Yatesの補正について

 院生の相談に乗っていて気付いたが、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の補正をすべきかどうかについては論戦が繰り広げられてきたところだし、独立性の検定よりも関連性を示した方が情報量が維持されるので良いのではないかとも思われる。

コメント

コメントを投稿

このブログの人気の投稿

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

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=...
続きを読む

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