対応のある多群間のノンパラメトリックな比較のためのFriedmanの検定後の対比較について

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弾）============================= 

結局、タイがあるため正確なp値を計算できないための違いではなく（実際、wilcox.test()とwilcox.exact()の差はわずかなものです）、３つのソフトではFriedmanの検定は同じだけれどもペアワイズの比較（対比較）をするためのpost hoc検定の方法が違うため、違う結果になるのは当然だというのが適切だと思います。

その後の調べで、SPSSにおけるFriedmanのpost hoc検定はDunn-Bonferroni
であることがわかりました。
https://www.ibm.com/support/pages/source-post-hoc-tests-friedman-nonparametric-tests

Dunn-BonferroniはRでもできて、dunn.test()というパッケージで可能です。
https://cran.r-project.org/web/packages/dunn.test/index.html
https://cran.r-project.org/web/packages/dunn.test/dunn.test.pdf
日本語のブログで使い方を説明している記事もありました。
https://note.com/eiko_dokusho/n/naac00eaec2ba

ただ、Dunn-Bonferroniはクラスカル＝ウォリス検定のpost hoc検定として
使われることも多く、対応のある検定になっているのかどうかが良くわか
りません。

jamoviが採用しているのは、Durbin-Conoverです。先のメールにも書きましたが、内部的にはPMCMRplusパッケージのfrdAllPairsConoverTest()という関数を
使っています。jamoviで出てくるp値は検定の多重性の補正がされていないので、5群間のすべての2群ずつの対比較なら10回の比較がされているので、Bonferroni補正ならp値をすべて10倍することになります。ConoverはFriedman
のpost hocでしか使われない検定なので、適用の適切性からいえば最善でしょう。

https://stats.stackexchange.com/questions/603451/post-hoc-tests-after-friedman-test

では、ウィルコクソンの符号順位検定やDunnの検定をBonferroni補正するのは
Friedmanのpost hoc検定としては適切でないのでConoverにしろ、と書かれて
います。Conoverの考え方や式は

https://real-statistics.com/anova-repeated-measures/friedman-test/friedman-test-post-hoc-analysis/

の説明がわかりやすいと思いました。難しいですが。

EZRでは、https://waidai-csc.jp/updata/2018/11/EZRNyumon.pdf
のp.39に書かれているように、2群ずつウィルコクソンの符号順位検定を実行して、Bonferroni補正の場合は得られたp値をすべて比較回数倍（5群間の場合なら10倍）していますが、この方法を使っている場合もよくあります。ウィルコクソンの符号順位検定ならば対応のある検定であることは間違いないので、その点を重視するなら、Friedmanのpost hocとしては
Dunn-Bonferroniより良いと思います。
https://www.datanovia.com/en/lessons/friedman-test-in-r/
は、この方法を説明しています。考え方としてはとてもシンプルです。

コメント

  1. Minato Nakazawa2024年3月27日 23:17

    この問題をシミュレーションで調べている論文 [https://www.tandfonline.com/doi/pdf/10.1080/03610918.2014.931971] が見つかった。jamoviが採用しているConoverは検討されていなかったが、EZRが使っているWilcoxonの符号順位検定について多重性の調整をする方法は、SPSSが使っているDunn-Bonferroniより検出力が低いと書かれていた

    返信削除

コメントを投稿

このブログの人気の投稿

差がゼロという帰無仮説以外の仮説の下でのp値の計算

イメージ
『宇宙怪人しまりす　統計よりも重要なことを学ぶ』朝倉書店（以下しまりす本３）の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,
続きを読む

baseグラフィクスの方が楽

イメージ
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活用が便利という発表主旨には賛成で、スライドも素晴らしいと思うが、世の中で広まっているやり方に引きずられるのは欠点だと思う。  
続きを読む

統計教育サンプルデータをAIで自動生成するのは無理だった

小説や漫画の登場人物のスペックを読者が勝手に想像して入力することで、統計教育に使えるオープンデータを作ろうという思いつきについて、AIでできるかも？　と思いついたことから、試してみたが、失敗したという話である。詳細は以下。 Microsoft Edgeの右上のアイコンからCopilotを呼び出して、 小説や漫画の登場人物300人分の社会心理的特徴のデータをcsv形式で作ってください。 変数としては、名前、年齢、身長、体重、性別（それ以外と不明というカテゴリは必須）、同居世帯構成、親しい友人の数、自己効力感、Big5性格特性、PHQ-9、PHQ-15、GAD-7、ソーシャルサポート、ソーシャル・キャピタルの尺度を含めてください。 と質問を打ってみた。すると以下の回答が表示された。 もちろんです。小説や漫画の登場人物に関する社会心理的特徴のデータをCSV形式で作成します。以下にリストされた変数を含めて、300人分のデータを生成します。 **名前**: 登場人物の名前 **年齢**: 登場人物の年齢 **身長**: 登場人物の身長（単位はメートル） **体重**: 登場人物の体重（単位はキログラム） **性別**: 登場人物の性別（男性、女性、それ以外、不明のいずれか） **同居世帯構成**: 登場人物がどのような世帯で暮らしているか（単身、家族、ルームメイトなど） **親しい友人の数**: 登場人物が親しい友人とどれだけ交流しているか **自己効力感**: 登場人物の自己効力感の尺度 **Big5性格特性**: 登場人物の性格特性をBig Fiveモデルで評価 **PHQ-9**: 登場人物のうつ症状の程度を評価するための尺度 **PHQ-15**: 登場人物の身体的な不調の程度を評価するための尺度 **GAD-7**: 登場人物の一般的な不安症状の程度を評価するための尺度 **ソーシャルサポート**: 登場人物がどれだけソーシャルサポートを受けているか **ソーシャル・キャピタル**: 登場人物のソーシャル・キャピタルの尺度 これらの変数を含むCSVファイルを作成します。お待ちください。 ソース: Bing との会話 2024/2/3 (1) 物語はいかにして心を動かすのか： 物語説得研究の現状と態度 .... https://www.jstage.js
続きを読む