統計相談に乗ったときのメモ

 院生から相談があったわけではないが、何回かのメモから採録する。    【第1807回】 国際保健医療学会西日本地方会など（2025年3月1-2日） いくつかの発表について三重大の谷村先生が指摘していたのは、 NEJMのStat Guideline にも書かれている prevalenceが高いときにはORは過大評価になるのでロジスティック回帰でなく修正ポアソン回帰を使ってCPRを出すべき （リンク先のプレゼン資料が大変わかりやすい）という話。参考に "Odds Ratio or Prevalence Ratio? An Overview of Reported Statistical Methods and Appropriateness of Interpretations in Cross-sectional Studies with Dichotomous Outcomes in Veterinary Medicine" (Frontiers in Veterinary Science, 2017)、 "Commentary: Calculating risk and prevalence ratios and differences in R: developing intuition with a hands-on tutorial and code" (Annals of Epidemiology, 2023)、 "Incorrect inference in prevalence trend analysis due to misuse of the odds ratio" (Annals of Epidemiology, 2016)、 "Overestimation of Relative Risk and Prevalence Ratio: Misuse of Logistic Modeling" (Diagnostics, 2022)、辺りは読んでおくべきか。ちなみに修正ポアソン回帰は、Rならsandwichパッケージかbootパッケージかrqlmパッケージを使うのが便利そう。   https://minato.sip21c.org/im3r/2025030...

以前投稿した、 https://stat-consul.blogspot.com/2024/03/friedman.html とも少し関連しているが、 対応のある3群以上の2値データの比較の話で相談に乗ったのでメモ。  jamoviを使っている院生からの相談に関連して。2時点間でMcNemarの検定をするような、0/1のカテゴリ変数の割合の変化がランダムではない一方向への変化なのかを、3時点以上で検定する方法としてはCochranのQ検定を使うのが普通だが、 jamoviでは特別なメニューが用意されていないので、Friedmanの検定をすればCochranのQ検定と同じp値が得られる 、という記述を見つけたので、確かめてみる。Rのパッケージとしては、RVAideMemoireパッケージに入っているcochran.qtest()関数が便利なので（DescToolsパッケージのCochranQTest()関数もほぼ同様である）、そのexample()の例で試してみると、 response <- c(0,1,1,0,0,1,0,1,1,1,1,1,0,0,1,1,0,1,0,1,1,0,0,1,0,1,1,0,0,1) fact <- gl(3,1,30,labels=LETTERS[1:3]) block <- gl(10,3,labels=letters[1:10]) library(RVAideMemoire) cochran.qtest(response ~ fact | block) とすれば、 Q = 10.8889, df = 2, p-value = 0.00432 というCochranのQ検定の結果、A、B、Cそれぞれにおける1の割合、AとB、BとC、CとAという対比較をfdr法で検定の多重性を調整した結果のp値が示される。反復測定デザインとしては、blockが個人、factが測定日に当たると考えれば良い。このデータをデータフレームddにしてタブ区切りテキスト形式でファイル出力するには、 dd <- data.frame(response=response, fact=fact, block=block) write.table(dd, "cochran.txt", row.names=FALSE, sep=...

尖度(kurtosis)と歪度(skewness)について質問があり、R/EZR/jamoviでの実装を調べてみた。 R本体には関数として実装されていないが、e1071パッケージに入っていて、helpを見たら、 Joanes and Gill (1988) に基づいて、それぞれ3通りの計算方法が実装されているとわかった。 古典的定義（歪度を3次のモーメントを2次のモーメントの平方根の3乗で割った値、尖度を4次のモーメントを2次のモーメントの二乗で割って３を引いた値とする定義。ちなみにr次のモーメントとは平均からの偏差のr乗和をサンプルサイズで割った値である）に基づくのが type=1 で、SASやSPSSと同じ実装が type=2 （e1071パッケージのヘルプでは表示されないが、Joanes and Gill, 1988によると、EXCELもこの実装を採用しているとのこと）、MINITABやBMDPと同じ実装が type=3 で、どちらも デフォルトは type=3 とのこと。 正規分布の下で不偏推定量であることがわかっているのは、歪度の方は３つともだが、尖度は type=2 だけとのことなので（なぜ type=2 をデフォルトにしなかったのか謎だが）、R本体で library(e1071) から skewness() や kurtosis() を使う場合は、特別な理由がない限り、 type=2 を指定すべきと思う。 EZRでは統計解析＞連続変数の解析＞正規性の検定でアウトプットウィンドウの中に表示されるが、 type=2 の結果と同じようだ。 jamoviは探索(Exploration)の記述統計(Descriptives)で変数を選んで計算させた後で統計量(Statistics)の分布(Distribution)の下にある歪度(skewness)と尖度(kurtosis)にチェックを入れると type=2 と同じ結果が表示される。 医学統計ではSASやSPSSと同じ結果になるのが無難な選択であろうから、EZRの方針は合理的だと思う。元々心理統計のために開発されたらしいjamoviも、最近は医学論文でも使われているから、妥当な選択であろう。 なお、EZRやjamoviでは欠損値があっても自動的に除去して計算してくれるが、e1071パッケージの skewness(...

Friedmanの検定後のpost hocな多重比較の結果がEZRとSPSSとjamoviで食い違うのは何故かという問い合わせメールをもらったので調べてみた。以下メールの返事を貼り付けて若干改変。 （返信第1弾）============================= jamoviは、 https://bookdown.org/sbtseiji/jamovi_complete_guide/sec-ANOVA-friedman.html に書かれている通り、Friedmanの後の対比較ではbonferroni補正なしのDurbin-Conover法を使っています。  このDurbin-Conover法というのは検定の多重性の調整をしていません。内部的にはPMCMRplusというパッケージに含まれている、frdAllPairsConoverTest()という関数を使っているようで、 この関数自体にはオプションでp.adj="bon"を指定できる のですが、jamoviではそこをいじれない実装になっています。 もっとも、Bonferroniの補正という手法は、対比較で出てきたp値に比較回数を掛けるだけなので、5群間の2群ずつの総当たり比較なら10回の比較がされているので、得られたp値を10倍すれば良いわけです。 （データを使った検証部分は省略） EZRが何故合わないかですが、EZRは対比較にWilcoxonの 符号順位検定を使っています（とメッセージにでてきます）。EZRでは正規近似でp値を出しているので、同順位の値が含まれていると、 wilcox.test.default(pairwise.response1, pairwise.response2, alternative = "two.sided", で警告がありました: タイがあるため、正確な p 値を計算することができません という警告が出ます。 ここで表示されているp値は、wilcox.test(..., paired=TRUE)で得られるp値の10倍になっています。同順位があっても正確なp値が欲しい場合は、exactRankTestsパッケージのwilcox.exact()関数を使えばOKです（もちろん個々のp値は10倍する必要があります）。 （返信第2弾）===...

『宇宙怪人しまりす　統計よりも重要なことを学ぶ』朝倉書店（以下しまりす本３）のpp.12「回復割合の差は20%」のp値を計算しようとしたが合わない。何かぼくのミスだと思うが原因がわからない。 x <- c(58, 62); n <- c(80, 100); X <- cbind(x, n-x) colnames(X) <- c("回復", "未回復") rownames(X) <- c("ヨクナール", "経過観察") print(X) s <- c(colSums(X)[1]/sum(X), colSums(X)[2]/sum(X)) # 帰無仮説：回復割合が処置と独立の場合の期待度数E print(E <- rbind(ヨクナール=rowSums(X)[1]*s, 経過観察=rowSums(X)[2]*s)) E/rowSums(E) # 確認：回復割合にヨクナールと経過観察で差はない 1-pchisq(sum((X-E)^2/E), 1) # 帰無仮説の下でのp値 Z <- c(0.1, -0.1) # ヨクナール群の回復割合が経過観察群より20%良い場合の期待度数EE print(EE <- rbind(ヨクナール=rowSums(X)[1]*(s+Z), 経過観察=rowSums(X)[2]*(s-Z))) EE/rowSums(EE) # 確認：回復割合がヨクナールの方が20％良い 1-pchisq(sum((X-EE)^2/EE), 1) # ヨクナールの回復割合が20%高い仮説の下でのp値 というRコードを実行すると、 差がないという帰無仮説（しまりす本３では「ゼロ仮説」）の下でのp値は本と同じく13.8%になるのだけれども、最後に表示される「回復割合の差が20%」という仮説の下でのp値は16.4%になってしまって、しまりす本３に書かれている17.2%とは合わない。しまりす本３ではpp.13-14に17.2%の算出方法が書かれていて、上のコードで使ったカイ二乗適合度検定とは違う方法（下記）なのだけれども、どちらでも一致するはずなんだがなあ。自分の頭の回転が鈍くて嫌になる。 rdp <- function(a, ...

gridグラフィクスが流行している理由が、やっぱり良くわからない。 例えば、 bob3bob3さんの『データ可視化学入門』をPythonからRに翻訳した話 で、日本の人口変化をChatGPT4ベースのCopilotを使ってPythonからRに翻訳させるとggplot2を使ったコードが出てくるようだが、この程度ならデータ定義部分だけ元コードからコピペして X <- data.frame(year = 1990:2022, pop = c(1.23611, 1.24101, 1.24567, 1.24938, 1.25265, 1.2557, 1.25859, 1.26157, 1.26472, 1.26667, 1.26926, 1.27316, 1.27486, 1.27694, 1.27787, 1.27768, 1.27901, 1.28033, 1.28084, 1.28032, 1.28057, 1.27834, 1.27593, 1.27414, 1.27237, 1.27095, 1.27042, 1.26919, 1.26749, 1.26555, 1.26146, 1.25502, 1.24947))   とデータフレーム X を定義しさえすれば、コードとしては、 plot(pop ~ year, data=X, xlab="西暦", ylab="日本の総人口[億人]", type="b", col="blue", pch=16, las=1); grid() で済む。パイプを使うよりも、データはいったんオブジェクトとして定義（付値）しておく方がRの思想的にもわかりやすいし、すべて標準のbaseグラフィクスで書けば、間違いなくどのR環境でも動くわけだから、Copilotすらggplot2を使うコードを出してくるのが謎でならない。念の為、bob3bob3さんのコード移植にはAI活用が便利という発表主旨には賛成で、スライドも素晴らしいと思うが、世の中で広まっているやり方に引きずられるのは欠点だと思う。