数学いらずの医科統計学PART7 CHAPTER43で、サンプルサイズについて書いてある。
パッケージならpwr、デフォルトならpower.*.test関数群でできる。
最近はやりのびっぐでーたなら、メモリに乗ってPCで解析できる分だけとにかく持っているデータを使えばいいんじゃない?と思うが、医療系だと莫大なコストと倫理的観点でサンプルがあまり稼げないことが多いのでサンプルサイズを前もって考えておきましょうという話。
サンプルサイズを決めるのは
有意水準：帰無仮説検定の棄却水準。たいてい0.05。小さくするとサンプルサイズは大きくなる。
検出力：差があるというときに、差があると判断する確率。たいてい0.8。業界による。type II error とは
の関係。大きくなるとサンプルサイズは大きくなる。
効果量：見たい差の大きさ。小さい差を見ようとするとサンプルサイズは大きくなる。
ばらつき：データを眺めて決まる。ばらつきが大きいとサンプルサイズは大きくなる。
この中で最も決めにくいのは効果量である。有意水準と検出力は業界の水準を流用すればいいが、効果量は自分で設定しなければならない。こういう困ったときにCohenさんは、大中小の効果量を事前に決めてあげればよくね?ということで適当な変換式を導入している。しかしCohenは行動科学の専門家であり、その専門領域に限って経験的にこれくらいの効果量が役立つといっており、すべての領域に普遍的に使えるものではない。普通は科学的な考察を加えてそこそこな効果量を決定する。いやだからそれが難しいんだってばよ…
比率、平均値、相関、ANOVA、χ自乗などで適当な効果量を計算する。
library(pwr) mapply(function(x) mapply(function(y) cohen.ES(test=x, size=y)$effect.size, c("small", "medium", "large")), c("p", "t", "r", "anov", "chisq", "f2"))
p t r anov chisq f2 small 0.2 0.2 0.1 0.10 0.1 0.02 medium 0.5 0.5 0.3 0.25 0.3 0.15 large 0.8 0.8 0.5 0.40 0.5 0.35
本では大きい効果量で26、中程度で65、小さいときで400くらいのサンプルサイズが必要と言っているので、これをCohenさんで計算すると、たいていそんな感じになる。
library(pwr) es <- mapply(function(x) cohen.ES("t", x)$effect.size, c("small", "medium", "large")) mapply(function(x) pwr.t.test(d=x, power=0.8)$n, es)
small medium large 393.40570 63.76561 25.52457
真面目にするならpower.*.testを使うほうがいい気がする。ただ、この関数群で使いにくい感じなのが、2群からサンプリングするときでないと対応していないっぽい感じなので、1群サンプリングのときはCohenさんを使うか、シミュレーションするか、確率分布から真面目に計算し直すか。
今、問題設定として、確率60%で当たると言われているクジがあるが、挑戦した人の話をまとめると40%くらいしか当たりが出ない雰囲気が漂っている。自分がクジを引いてこれを確かめるとき、何回引いたらこのクジが40%であることを示せるか
という感じ…?
有意水準0.05、検出力0.8として、効果量はクジが当たると言われている0.6と、いやウソだろう…という0.4を考える。
シミュレーションで、0.4でサンプリングした時に帰無仮説0.6が棄却されるのが、検出力分だけあるかを、サンプルサイズを増やしながらやってみる。今回は
偶奇で振動してるんだが…面倒なので保留。
# クッソ時間かかるよ!! p1 <- 0.4 p2 <- 0.6 alpha <- 0.05 beta <- 0.2 # power = 1- beta n <- 1 res <- NULL hoge <- 0 niter <- 2000 # 1サンプルサイズでのシミュレーション回数 while(mean(hoge) <= 1 - beta){ # hoge は検定結果。差があることが前提なので、検出力分だけ有意な結果が得られるかどうか。 hoge <- replicate(10, mean(replicate(niter, binom.test(sum(sample(0:1, size=n, prob=c(1-p1, p1), replace=TRUE)), n, p2)$p.value) < alpha)) res <- rbind(res, hoge) n <- n + 1 } rownames(res) <- NULL boxplot(t(res), xlab="sample size", ylab="power") lines(seq(nrow(res)), rowMeans(res), lty=2, col=4, lwd=3)
pwrでもできる。片側検定だとけっこうずれてしまった。両側だといいかんじだった。
pwr.p.test(ES.h(p1, p2), n=NULL, sig.level=alpha, power=1-beta, alternative="less")
proportion power calculation for binomial distribution (arcsine transformation)
h = -0.4027158
n = 38.12156
sig.level = 0.05
power = 0.8
alternative = less
