「統計学関連なんでもあり」の過去ログ--- 039

No.01897　パワーアナリシス　　【tr】　2006/12/14(Thu) 00:19

パワーアナリシスをしようとしています。
例えば，青木先生がお書きになった Rstat.pdf の p.56 からの
「6.6 独立k 標本の平均値の差の検定（一元配置分散分析）」
の例で行くと，oneway.test ではなく，

> summary(aov(df$結果 ~ df$群))
            Df Sum Sq Mean Sq F value Pr(>F)
df$群        2   84.1    42.1  0.4469  0.641
Residuals   97 9129.7    94.1

として，結果の df$群と Residuals の Mean Sq を使い，

> power.anova.test(groups=3, n=33, between.var=42.1/33, within.var=94.1, sig.level=0.05)

     Balanced one-way analysis of variance power calculation 

         groups = 3
              n = 33
    between.var = 1.275758
     within.var = 94.1
      sig.level = 0.05
          power = 0.1209878

 NOTE: n is number in each group 

のようにして power = 0.1209878 と求めればよいのでしょうか。

No.01902　Re: パワーアナリシス　　【青木繁伸】　2006/12/14(Thu) 09:50

それでたぶんよろしいと思います。

No.01909　Re: パワーアナリシス　　【tr】　2006/12/14(Thu) 21:03

ありがとうございます。

最初に R で power.anova.test を使ったところ，
G*Powerで求めた場合と大幅に違う値が出ていたので悩んでいたところでした。

http://www.okada.jp.org/RWiki/index.php?cmd=read&page=%C4%B6%CC%F5%A1%A7%A5%D1%A5%EF%A1%BC%A5%A2%A5%CA%A5%EA%A5%B7%A5%B9%A1%CA%B0%EC%B8%B5%C7%DB%C3%D6%CA%AC%BB%B6%CA%AC%C0%CF%A1%CB
の情報では between.varには群間分散そのものを与えてはダメで，
群間分散/n を与えなければならないとのことで，このように計算してみました。

＃そのページにも書いてありますが，一体どういうわけなんでしょうかね。

重ねての質問になってしまいますが，p. 53
「6.4 独立2 標本の平均値の差の検定（t 検定）」
の場合では，
> d <- abs(mean(df$結果[df$群=="g1"])-mean(df$結果[df$群=="g2"]))
> s <- sd(df$結果)
> power.t.test(n=50, delta=d, sd=s, sig.level=0.05, type="two.sample", alternative="two.sided")

     Two-sample t test power calculation 

              n = 50
          delta = 3.208625
             sd = 9.64723
      sig.level = 0.05
          power = 0.3770338
    alternative = two.sided

 NOTE: n is number in *each* group 
で，power = 0.3770338 というような感じでしょうか。
一標本の場合や対応のある場合は type が "one.sample" なり "paired" になるだけ，
と思っております。

No.01910　Re: パワーアナリシス　　【青木繁伸】　2006/12/14(Thu) 21:33

s を求めるときに，t検定のときと同じように両群の標準偏差をプールして求めるやり方もありますね。
いずれにしろ，パワーアナリシスに用いる数値は，若干の誤差を持っていても差し支えないというか，この数値がこれくらいのときにパワーはどれくらいになるのかなみたいに使われるので。

対応のある場合は最終的には一標本の場合に帰着できる（つまり対応のあるデータの差を取ってそれを使うとすれば一標本だから）わけですが，パラメータを変えるといろいろな場合に対応しているというのはありがたいですね。

No.01916　Re: パワーアナリシス　　【マスオ】　2006/12/15(Fri) 00:39

> ＃そのページにも書いてありますが，一体どういうわけなんでしょうかね。

between.var は群平均の分散です．
提示例は群内の繰り返し数 nが不揃いなのでややこしいのですが，nが各群同じであれば，
群平均の分散 = 平均平方和（MS）/ n
になるから．
> set.seed(123)
> ng <- 3
> n <- 30
> g <- gl(ng, n)
> x <- rnorm(n*ng)*10+80
> (m <- tapply(x, g, mean)) #群平均
       1        2        3 
79.52896 81.78338 80.24420 
> var(m) #群平均の分散
[1] 1.327176
> 
> (s <- anova(lm(x ~ g)))
Analysis of Variance Table

Response: x
          Df Sum Sq Mean Sq F value Pr(>F)
g          2   79.6    39.8  0.4943 0.6117
Residuals 87 7007.6    80.5               
> s["g", "Mean Sq"]/n
[1] 1.327176

No.01917　Re: パワーアナリシス　　【マスオ】　2006/12/15(Fri) 00:44
> s <- sd(df$結果)
さすがにこれはまずいのでは？
試しに平均の差をいろいろ変えたサンプルを作ってみて，プールした分散と比べてみてください．

No.01923　Re: パワーアナリシス　　【tr】　2006/12/15(Fri) 21:40

between.var は「群間分散÷標本数」ではなくて「群平均の分散」であり，
たまたま各群の標本数が一致する場合には「群間分散÷標本数」と同じになるだけ，
ということですね。

そうすると最初の計算の例でいくと，

between.var には

> var(tapply(df$結果, df$群, mean))
[1] 1.254219

として 1.254912 を使い，
within.var には

> summary(aov(df$結果 ~ df$群))
Df Sum Sq Mean Sq F value Pr(>F)
df$群 2 84.1 42.1 0.4469 0.641
Residuals 97 9129.7 94.1

として，結果の Residuals の Mean Sq である 94.1 を使い，

> power.anova.test(groups=3, n=33, between.var=1.254219, within.var=94.1, sig.level=0.05)

Balanced one-way analysis of variance power calculation

groups = 3
n = 33
between.var = 1.254219
within.var = 94.1
sig.level = 0.05
power = 0.1197049

NOTE: n is number in each group

といったところでしょうか。
もちろん，標本数が 33 ではない群がありますので，
実際には多少違うと思っております。

No.01924　Re: パワーアナリシス　　【tr】　2006/12/15(Fri) 21:41

ご指摘の通り，s を求めるために全体の標準偏差をつかってはマズイですね。
例の場合は，同一の母集団から 2 群を取り出していたので，
全体の標準偏差でもいいのかもしれませんが，
普通は同一の母集団から取り出したと分かっている
ものの差を検定するなんてことはないですよね。

普通に考えても，
等分散で平均が違う 2 群をゴッチャにして分散なり標準偏差なりを求めると
元の分散なり標準偏差なりより大きく出てしまいますね。

t 検定のように両群の標準偏差をプールするやり方，
ということで以下のようにしてみました。
> d <- abs(mean(df$結果[df$群=="g1"])-mean(df$結果[df$群=="g2"]))
> vg1=var(df$結果[df$群=="g1"])
> vg2=var(df$結果[df$群=="g2"])
> ng1=length(df$結果[df$群=="g1"])
> ng2=length(df$結果[df$群=="g2"])
> v=((ng1-1)*vg1+(ng2-1)*vg2)/(ng1+ng2-2)
> s=sqrt(v)
> power.t.test(n=50, delta=d, sd=s, sig.level=0.05, type="two.sample", alternative="two.sided")

     Two-sample t test power calculation 

              n = 50
          delta = 3.208625
             sd = 9.560432
      sig.level = 0.05
          power = 0.3827252
    alternative = two.sided

 NOTE: n is number in *each* group

No.01925　Re: パワーアナリシス　　【青木繁伸】　2006/12/15(Fri) 22:18

パワーアナリシスの使い方としては，検定や推定みたいに，ある特定の状況での結果を得るというのではなく，注目しているパラメータを少しずつ変化させたときに結果がどのように変化するかをシミュレーションしてみるというような使い方もあるのではないかと思います。
G*powerみたいに，グラフで表して，妥協点を探すみたいな使い方もできる（というか，実験デザインをたてるときの参考にする）

No.01978　Re: パワーアナリシス　　【tr】　2006/12/20(Wed) 22:16

遅くなりましたが，ありがとうございます。
標本数をどれぐらいにするのか，とか
どれぐらいの差を検出することを考えるのか，とか
グラフで表すとわかりやすいですね。

不勉強だからか，
この手の話が書いてある文献を見つけることができなくて，
こちらで質問させていただきました。
大変役に立つご回答をいただくことができました。

青木先生，マスオ様，ありがとうございました。