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

自分が指導しているのではない院生から表題の計算をしてくれるツールがないか尋ねられた。 Rのfisher.test()では2x2に限らずカテゴリ３つ以上の変数が含まれていても独立性の確率を計算してくれるが、パッケージを探してもサンプルサイズと効果量と有意水準を与えて検出力を計算してくれる関数は見つからなかった。 2x2の場合ならEZRでもG*Powerでもjamoviのモジュールでもできるのだが、2x3以上についてはシミュレーションプログラムを書くしかなさそうだった。 Googleで"power.fisher.test multiple categories"と検索してみたら、AIによる回答が以下のコードを提案してきたが、これは正しいのだろうか？　効果量を与えるのではなく、各セルの期待確率を与える必要があるようだが。とりあえず後で考えるためメモだけしておく。 # Define a function to simulate power for an r x c table sim_power_fisher rc <- function (matrix_p, n_total, alpha = 0.05 , nsim = 1000 ) {    # matrix_p is a matrix of expected cell probabilities    p_vals <- replicate(nsim, {     # Generate multinomial counts     sim_data <- rmultinom( 1 , size = n_total, prob = matrix_p)   dim(sim_data) <- dim(matrix_p)    # Run Fisher's exact test (uses hybrid or Monte Carlo for > 2x2)     res <- fisher.test(sim_data, simulate.p.value = TRUE , B = 2000 )   res$p.value  }) ...

カテゴリ変数について割合の信頼区間を計算するには、通常、二項分布を仮定してClopper and Pearsonの方法を使う（ 『Rによる保健医療データ解析演習』 の第9章、9.3に書いた通り）。しかし、カテゴリが3つ以上ある場合、各カテゴリの割合は独立ではない（もちろん2つの場合も片方のカテゴリの割合が決まれば、残りも自動的に決まってしまうわけだが、1つのカテゴリが生起する確率、と考えれば二項分布を普通に適用して差し支えない）ので、別のアプローチをした方が良いと思われた。 自分の本にも書かなかったけれども、カテゴリ３つ以上の場合の、各カテゴリの割合の信頼区間の推定についても、既に多くの研究があり、Rでも実装しているパッケージがある。 DescToolsパッケージのMutinomCI()関数 である（もっとも、 DescToolsパッケージのvignettesの１つ の13.1の解説にあるように、カテゴリ３つ以上でも二項分布を当てはめて信頼区間を求める、BinomCI()という関数も用意されているので、必ずしも各カテゴリについて二項分布で信頼区間を推定して悪いわけではないと思われる→正しくは、個々のカテゴリの割合についての95％信頼区間を求めるのか、すべてのカテゴリの割合について95％同時信頼区間を求めるのか、という違いで、どちらも間違いではない。黒木さん、ご指摘ありがとうございました）。 MultinomCI()関数にはさまざまな信頼区間推定方法が実装されていて、古くからあり他のパッケージでも実装されているものもあるらしいGoodman (1965)の方法も、method="goodman"というオプションを指定すれば使えるが、デフォルトがSison and Glaz (1995)の方法（明示的にオプション指定するにはmethod="sisonglaz"だが、省略できる）なので、特別な理由がなければそれで計算するのが良いと思う。DescToolsパッケージを引用するには、Signorell A (2025). _DescTools: Tools for Descriptive Statistics_.  doi:10.32614/CRAN.package.DescTools , R package version 0...

 院生の相談に乗っていて気付いたが、chisq.test()やprop.test()は、2行2列でないと、correct=TRUEをつけてもYatesの連続性の修正をしてくれない。強引にやらせる方法としては、例えば Hitchcock, David B. (2009). Yates and Contingency Tables: 75 Years Later. Retrieved 4/8/2015 from: University of South Carolina に書かれている式を使うと、   chisqY.test <- function(.X) {  .EX <- outer(rowSums(.X), colSums(.X))/sum(.X)  # same as chisq.test(.X)$expected  X2Y <- sum((abs(.X-.EX)-0.5)^2/.EX)  DF <-  (NROW(.X)-1)*(NCOL(.X)-1)  # same as chisq.test(.X)$parameter list(chisq.Yates=X2Y, df= DF,  p.value=1-pchisq(X2Y, DF)) } と関数定義してから、chisqY.test()を使えば、行数か列数、あるいはその両方が3以上でも、強引にYatesの補正をしてくれる。 なお、当初は内部で期待度数と自由度を得るためにchisq.test()を呼び出すとき、期待度数が5より小さいセルがセル数の20％以上あると警告メッセージが表示されたので、chisq.test()を呼ばない形に書き換えた。 また、chisq.test()やprop.test()のYatesの補正では、0.5ではなく、min(0.5, abs(.X - .EX)を補正項としているので、2行2列の場合でも結果は必ずしも一致しない点にも注意。 けれども、前掲Hitchcock (2009)にあるように、そういう場合は素直にfisher.test()を使えば良いので、無理にYatesの補正をするのは筋が良くない。そもそも、Yatesの補正をすべきかどうかについては論戦が繰り広げられ...

 院生から相談があったわけではないが、何回かのメモから採録する。    【第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(...