pwr.anova.test() と power.anova.test() 両方あって,面倒くさい。同じ条件で計算したとき,どちらも同じ答えにならないと困るので,調べてみた。
分散分析のサンプルサイズ計算―pwr.anova.test()を使う方法 でも取り上げられていたのだけど,結論がちょっと違ったので,まとめておく。
まず,両方のプログラムで p.body で定義される関数本体の違いを見る。
power.anova.test で群の個数を表す変数が groups なので,pwr.anova.test と揃えるために k に置換して両者を以下に示す。
power.anova.test
lambda <- (k - 1) * n * (between.var/within.var)
pf(qf(sig.level, k - 1, (n - 1) * k, lower.tail = FALSE),
k - 1, (n - 1) * k, lambda, lower.tail = FALSE)
pwr.anova.test
lambda <- k * n * f^2
pf(qf(sig.level, k - 1, (n - 1) * k, lower = FALSE),
k - 1, (n - 1) * k, lambda, lower = FALSE)
両者は lamda の計算式が異なる。しかし,power.anova.test で between.var と within.var を使おうが,pwr.anova.test で between.var と within.var から計算される f を使おうが,同じ結果になるためには between.var と within.var と f の関係式を調べればよい。
つまり,(k - 1) * n * (between.var/within.var) = k * n * f^2 ならば,
f = sqrt((between.var * (k-1)/k) / within.var)
さて,between.var とは何者かということもあるが,これは
各群の平均値を groupemean としたとき,between.var = var(groupmean) で,一元配置分散分析表の 群間の平均平方 Mean Sq を(各群で同じ)サンプルサイズ n で割ったものに等しい。(最後に書くけど,var() がくせ者)
前出の,分散分析のサンプルサイズ計算―pwr.anova.test()を使う方法 では,
群間の平均平方 Mean Sq 484.7, n = 4, 群内の平均平方 35.9 として,効果量 f=sqrt(484.7/4)/sqrt(35.9) とした計算結果が示されている。
> pwr.anova.test(k=3, f=sqrt(484.7/4)/sqrt(35.9), power=0.8)
Balanced one-way analysis of variance power calculation
k = 3
n = 2.289809 # 一番最後の結果に示した power.anova(... between.var=484.7/4 * (3/2), ...) の結果と同じ
f = 1.837212
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
これを power.anova.test() でやってみると,結果が異なる。
> power.anova.test(groups=3, between.var=484.7/4, within.var=35.9, power=0.8)
Balanced one-way analysis of variance power calculation
groups = 3
n = 2.712944 # 2.289809 と合わない
between.var = 121.175
within.var = 35.9
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
そこで,上に示したように効果量 f の計算に使う between.var を修正して分析する。
> pwr.anova.test(k=3, f=sqrt(484.7/4 * 2/3)/sqrt(35.9), power=0.8)
Balanced one-way analysis of variance power calculation
k = 3
n = 2.712935 # 誤差範囲で 2.712944 と同じになった
f = 1.500077
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
この比較は power.anova.test を基準にしたのだが,pwr.anova.test を基準にしても構わない。
> power.anova.test(groups=3, between.var=484.7/4 * (3/2), within.var=35.9, power=0.8)
Balanced one-way analysis of variance power calculation
groups = 3
n = 2.289779 # 最初に示した pwr.anova.test の結果と一致する
between.var = 181.7625
within.var = 35.9
sig.level = 0.05
power = 0.8
NOTE: n is number in each group
いずれを基準としてもよいのだけど,別々では困る。
between.var = var(groupmean) の自由度は R では k - 1 だから,これに基づいて計算する場合は (k-1)/k の修正が必要かなあ?と思う次第。つまり,power.anova.test() を基準とする。
※コメント投稿者のブログIDはブログ作成者のみに通知されます