CochranのQ検定の話

以前投稿した、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="\t")


とすれば、Rの作業フォルダにcochran.txtというファイルができる。これは縦長形式なので、EZRやjamoviで分析するには横長形式に変換する必要がある。最も簡単なのは、EZRであればddというデータフレームにファイルを読み込んだ後で、

library(car)

ddw <- reshape(dd, direction="wide", idvar="block", timevar="fact")


として横長形式のddwに変換し、ddwをアクティブデータセットにしてから、統計解析＞名義変数の解析＞対応のある3群以上の比率の比較（Cochran Q検定）と進んで、時点を示す3つの変数を選んでOKをクリックするという方法である。X-squared = 10.889, df = 2, p-value = 0.00432と同じ結果が得られる。ノンパラメトリック検定＞対応のある3群以上の間の比較（Friedman検定）でも、同じp値が得られる。jamoviでも、横三本線のメニューの「開く」からブラウズしてデータを読み込んだ後、jReshapeモジュールのLong to wideから、Rows to columnsとしてresponse、indexing variableとしてfact、ID variableとしてblockを指定し、Reshapeというボタンをクリックすると、横長形式への変換が行われ、新しいjamoviのウィンドウとして自動的に開くので、その後で、分散分析＞反復測定分散分析（ノンパラメトリック）と進んで3つの変数を選ぶと、表示されるp値は0.004である（同じ0.00432だが、丸められて0.004までしか表示されない）。なお、RVAdeMemoireパッケージのcochran.qtest()関数ではfdr法で検定の多重性を調整した対比較結果が表示されるが、jamoviではダービン・コノバー法が可能である（逆にfdrやholmの結果を得るメニューはない）。記述統計にチェックを入れれば、平均値として1の割合を計算させることもできる。

 

コメント

コメントを投稿

このブログの人気の投稿

統計教育サンプルデータ作成プロジェクト

このところ保健学研究共通特講IV, VIIIにおける統計教育の改善方法についてつらつら考えていて（というか結果的には無駄になってしまった某業務のため、強制的に内容の再整理を考える必要があったのだが）、自由に使える適当なデータが欲しいと思うようになった。作図の説明のためには、量的なデータと質的なデータが混在している何百人規模のものが適していて、これまではMASSパッケージのsurveyという、アデレード大学の統計関係の学生のデータを使ってきた。自分のデータを公開するのは、匿名化したとしても倫理審査関係のクリアが難しいなあと思って、そこで思考が停まっていたが、ふと、架空の人物について想像したデータをバザール形式で作ってしまえば良いのではないかと思いついた。 質問紙はウェブのフォームで作り、完全匿名で（何の登録も必要ないゲストとして）ボランティアを募って（botに入力されないような仕組みは必要だと思うが、メール認証までしなくても、良くある CAPTCHA 、例えばPerlのcgiなら このモジュールを使う とか、あるいは Turstile で十分だろう）、自分が好きな小説や漫画の登場人物について作品中で明らかになっている情報や、性格的には読者が想像した情報で、登場人物の気になって埋めてもらうのだ。そうすれば、個人情報にも抵触しないし、倫理的な問題も生じないはずだ。データやプロジェクトのタイトルを 読者が勝手に想像した架空人物の社会心理的特徴 、とでもするか。 たぶん架空人物も成長するので、何歳時点のデータというのを想像した場合は何度でも入力できるようにしたい。例えば20代のシャーロック・ホームズと、ライヘンバッハの滝から生還後の40代のシャーロック・ホームズでは別データとなる。高村薫『李歐』の李歐や吉田一彰は、22歳の時と再会後では別データとなる。名前、年齢、身長、体重（を入れるためには、人間に限るとしなくてはデータにならないが、明らかに定量的なデータとしては身長と体重くらいしか想像で入力できるものが思いつかない）、性別（架空人物の場合は、新井素子『二分割幽霊奇譚』の主人公のようなインターセックスの登場人物は珍しくないので、それ以外と不明というカテゴリは必須だな）、同居世帯構成（ワンピースのサニー号乗組員のような場合もあるので、血縁関係なしが何人、というカテゴリが実在す...
続きを読む

差がゼロという帰無仮説以外の仮説の下での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, ...
続きを読む

対応のある多群間のノンパラメトリックな比較のための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弾）===...
続きを読む