pwr.anova.test() と power.anova.test()

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() を基準とする。