という質問を同級生から受けた。
氏いわく、リスクの有無と疾患の有無で分割表をフィッシャー検定したら、R ではp=1 でエクセルではp=0.7 だった。たぶんp=0.7 っぽいからエクセルの解析を信じたいのだけど、なぜR でp=1 になったのか、とのこと。
そしてその分割表データとRスクリプトが添付されてた。
b2 <- matrix(c(14, 21, 10, 17), 2) fisher.test(b2)
Fisher's Exact Test for Count Data data: b2 p-value = 1 alternative hypothesis: true odds ratio is not equal to 1 95 percent confidence interval: 0.3588741 3.6373238 sample estimates: odds ratio 1.131049
R でやったんならR が正しいだろ…と思ったが、それはそれでR 盲信者だということに反省したので、原点に立ち返って定義からやってみる。
確率の基本は、受験数学では全通りの数え上げである。フィッシャーの正確確率検定は
という分割表において
で与えられる。
上の分割表は、自由度1、すなわち、 という整数がひとつ動くと、周辺度数を満たすためには他すべてが連動して動く。
すべてのセルは0 以上だから、考えられる表のパターンはいかの通り。
idx <- matrix(c(1, -1, -1, 1), 2) i <- -14:10 bs <- mapply(function(z) b2 + idx*z, i, SIMPLIFY=FALSE)
各々の分割表の確率を計算すれば、「考えられるすべての分割表」は数え上げられているので、確率の和は1になる。
f <- function(mat) choose(sum(mat[1,]), mat[1,1])*choose(sum(mat[2,]), mat[2,1])/choose(sum(mat), sum(mat[,1])) p <- sapply(bs, f) # Fisher exact probability plot(p, type="s", ylab="probability", main="Fisher exact probability", lwd=5) abline(v=15, lty=3, lwd=3)
b2 の分割表は15番目で、これはもっとも「起こりうりそう」な分割表なので、これより生起確率の小さい分割表の確率(b2 自身を含む)を足すと、確率1 になる。
というわけで、たとえば片側検定をしようとして、fisher.test のalternative にless もしくはgreater を指定して、15番目より小さい確率もしくは大きい確率を足してみると、一致していることがわかる。
fisher.test(b2, alternative="less")$p.value; sum(p[1:15]) fisher.test(b2, alternative="greater")$p.value; sum(p[15:length(p)])
というわけで、R でやったのがやはりよかった。
逆に、エクセルでなぜp=0.7 になったのかが再現性がまったくないのでわからない。エクセルでフィッシャー検定をしようとするとエクセル統計のアドインが必要っぽいので、氏が持っていたとは到底思えないので謎。
