それが、power.t.test()です。
今回はこの関数を用いて、検出力と例数設計について検証してみたいと思います。
先日、身長平均身長175cm(±5cm）と他の集団との身長の平均値に関しての比較を扱いました。その時、有意差を検出するためのn数について直感的な検討を行いました。
今回は、パワーテストなので、より定量的に議論したいと思います。
以下のスクリプトは、1)サンプル数、標準偏差、α危険率、検出力をある値に想定したときに検出可能な平均値の差 2)想定される平均値の差、標準偏差、α危険率、検出力をある値に想定したときに要請されるn数 3)想定される平均値の差、n数、標準偏差、α危険率、検出力をある値に想定したときの検出力の値 について着目しています。
今回はベースに「平均身長175cm(±5cm)の集団と平均身長180cm(±5cm)の集団の平均値の比較を意識していて、サンプル数は１０、p値は0.05を想定している」状況があります。
$ vim command.in
#以下の条件を満たすのに必要なn数を求める。
n_null = power.t.test(delta=5,sd=5,sig.level=0.05,power=0.8,type=c("two.sample"),alternative="two.sided")
#以下の条件で有意差が検出可能な母平均の差を求める。
delta_null = power.t.test(n=10,sd=5,sig.level=0.05,power=0.8,type=c("two.sample"),alternative="two.sided")
#以下の条件での検出力を求める。
power_null = power.t.test(delta=5,n=10,sd=5,sig.level=0.05,type=c("two.sample"),alternative="two.sided")
$ R
> source("command.in")
#特定の項目を取り出す方法を調べる。
> names(delta_null)
[1] "n" "delta" "sd" "sig.level" "power"
[6] "alternative" "note" "method"
#以下、1)-3)のケースの計算結果を表示
> delta_null$delta
[1] 6.624728
> n_null$n
[1] 16.71477
> power_null$power
[1] 0.5619846
1)サンプル数、標準偏差、α危険率、検出力をある値に想定したときに検出可能な平均値の差
> delta_null$delta
[1] 6.624728
サンプル数は１０、p値は0.05、検出力0.8、標準偏差5を想定しているで検出可能な平均値の差は6.6cmであるということです。
> n_null$n
[1] 16.71477
平均値の差が5、p値は0.05、検出力0.8、標準偏差5の条件を満足させるために必要なn数は16.7である。
#実験計画を立てる上で特に重要である。
3)想定される平均値の差、n数、標準偏差、α危険率、検出力をある値に想定したときの検出力の値
> power_null$power
[1] 0.5619846
n数が10、平均値の差が5、p値は0.05、標準偏差5の条件の時の検出力は0.56である。
2種類の過誤を明確に意識し、さらには母集団の平均値の差を実験によりある程度想定しながら、サンプル数を決定したいと思いました。
・第一種過誤（α過誤、偽陽性）は、帰無仮説が実際には真であるのに棄却してしまう過誤
・第二種過誤（β過誤、偽陰性）は、対立仮説が実際には真であるのに帰無仮説を採用してしまう過誤
最後に実際に実験を行うことを意識して
「想定する母集団（標準偏差５）の差を1から0.5cmずつ10cmまで増加させていったときに、p値は0.05、検出力0.8の条件を満足させるために必要なn数を２次元プロットに描写してみます。
$ vim command_2.in
n <- c(1:100)
for(i in 1:100)
{
n[i] <- (power.t.test(delta=0.1*i,sd=5,sig.level=0.05,power=0.8,type=c("two.sample"),alternative="two.sided"))$n
}
delta <- c(1:100)*0.1
png("120801_delta_n_plot.png")
plot(delta,n)
dev.off()
$ fg
> source("command_2.in")
このシミュレーション結果は非常に深い意味があります。
平均値の差が小さい時(2ぐらいまで）は、非常に多くのサンプル数が必要となるのに対して、それよりも大きい平均値の差であるときには、サンプル数がほぼ一定（5程度）で済んでしまうというわけです。
deltaの値を0.11から10までの間に関してプロットを書かせるとよりわかりやすいです。
> plot(delta[11:100],n[11:100])
実験を計画するときはこのような事実を念頭に置いてサンプル数を設定したいものです。
