No.21941　分散分析のpower analysisについて　　【Kei】　2016/03/23(Wed) 20:41

分散分析で有意差が出なかった研究結果について，post hocでパワーアナリシスをするように求められており質問させていただきます。

標準パッケージに含まれているpower.anova.test（群間分散，群内分散を使用）と
pwrパッケージに含まれているpwr.anova.test（効果量を使用）の結果が大きく異なり
使い方が間違えているのではないかと思い質問させていただきます。

例えば
group <- c("a","a","a","a","a","b","b","b","b","b","c","c","c","c","c")
data <- rnorm(15)
testdata <- data.frame(group,data)
anova(aov(data~group,data=testdata))
として結果で表示される

Df Sum Sq Mean Sq F value Pr(>F)
group 2 3.0252 1.5126 1.1924 0.337
Residuals 12 15.2231 1.2686

にあるgroupのMean Sqを群間分散，ResidualsのMean Sqを群内分散として
用いてよいのでしょうか。この結果から

a <- anova(aov(data~group,data=testdata))
power.anova.test(group = 3, n = 5, between.var = a[1,3], within.var = a[2,3], sig.level = 0.05, power = NULL)

として計算されるPowerと，以下の方法の方法で計算される結果が大きく異なってしまいます

# 効果量es = sqrt(群間の平方和/群内の平方和)

es <- sqrt(a[1,2]/a[2,2])
pwr.anova.test(k = 3, n = 5, f = es, sig.level = 0.05, power = NULL)

ご多忙のところ恐縮ですがご教授いただけると助かります。
よろしくお願いいたします。
No.21942　Re: 分散分析のpower analysisについて　　【青木繁伸】　2016/03/24(Thu) 13:47

power.anova.test の中で，F 分布の非心度 lambda の計算をするところがあります

lambda <- (k - 1) * n * (between.var/within.var)

pwr.anova.test では，

f <- sqrt(between.var / within.var)
lambda <- k * n * f^2

となっております

(k-1) を掛けるか k を掛けるかの違いになっているわけです

どちらが正しいのか，参考書を全て処分したので，私にはわかりません
No.21943　Re: 分散分析のpower analysisについて　　【Kei】　2016/03/24(Thu) 18:27

青木先生

早々のお返事ありがとうございます。先生のご返信を見て気づいたのですが，

within.var(群内分散)と
Analysis of Variance Tableにある
residualsのMean Sq(群内平均平方)，Sum Sq(群内平方和)
を私は混同してしまっているかもしれません。

power.anova.testを使うとき私は
between.varに群間のMean Sqを，within.varに群内のMean Sqを
代入しました。

また，効果量を求めるときに
青木先生はf <- sqrt(between.var / within.var)とされていますが
私はf <- sqrt(群間のSum Sq/群内のSum Sq)としました。

このどちらかが間違っていると思うのですが，ご教授いただけませんでしょうか。
No.21944　Re: 分散分析のpower analysisについて　　【青木繁伸】　2016/03/24(Thu) 19:15

> また，効果量を求めるときに
> 青木先生はf <- sqrt(between.var / within.var)とされていますが
> 私はf <- sqrt(群間のSum Sq/群内のSum Sq)としました。

f <- sqrt(between.var / within.var) は，プログラムにそう書かれていて，pwr.anova.test で実際に行われている演算だからそう書いたまで。

ところで，ここでいうところの between.var，within.var は，データから計算された分散分析表での MSb, MSw とは本来別のものであるということがそもそもの勘違い。

k 群あって，各群のサンプルサイズが同一の n であるとき，各群の平均値を μi，全体の平均値を μ としたとき，between.var は Σ(μi-μ)^2 / k。すなわち，
f = sqrt{ (Σ(μi-μ)^2 / k) / within.var } 。
within.var は見積もられるものですから，いろいろな値をとりうるけど，MSw ではない。

example(power.anova.test) の最初の例で，
power.anova.test(groups=3, n=4, var(c(59, 66, 42)), 12^2, sig.level=0.05)
が挙げられている。12^2 が within.var。これは天下り的に来ているわけではないが MSw とは無関係。
between.var は var(c(59, 66, 42)) とされているが，上の f の定義でいうと var(c(59, 66, 42))*(k-1)/k
lambda = k * n * var(c(59, 66, 42))*(k-1)/k / 12^2
= (k-1) * n * var(c(59, 66, 42)) / 12^2
となる。
一方，pwr.anova.test は，引数に f を与える。
繰り返すが，f = sqrt(var(c(59, 66, 42))*(k-1)/k / 12^2) で，内部では lambda の計算に f^2 が使われる。
lambda = k * n * f^2
= k * n * var(c(59, 66, 42))*(k-1)/k / 12^2
= (k-1)*n*var(c(59, 66, 42)) / 12^2

結局の所，以上を理解した上で使えば，power.anova.test も pwr.anova.test も同じ結果になる。
> power.anova.test(groups=3, n=4, var(c(59, 66, 42)), 12^2, sig.level=0.05)

     Balanced one-way analysis of variance power calculation 

         groups = 3
              n = 4
    between.var = 152.3333
     within.var = 144
      sig.level = 0.05
          power = 0.5846256

NOTE: n is number in each group

> 
> f <- sqrt(var(c(59, 66, 42))*(k-1)/k / 12^2)
> f
[1] 0.8397898
> pwr.anova.test(k=3, n=4, f=f, sig.level=0.05)

     Balanced one-way analysis of variance power calculation 

              k = 3
              n = 4
              f = 0.8397898
      sig.level = 0.05
          power = 0.5846256

NOTE: n is number in each group
No.21945　Re: 分散分析のpower analysisについて　　【青木繁伸】　2016/03/24(Thu) 19:37

以下のように，MSw = n * between.var なので，within.var = MSe / n と指定してやると，データに基づく power の予測になるでしょう
この場合は比例するわけだから，MSb，MSw をそのまま使っても結果は同じになる
> k <- groups <- 3
> n <- 10
> x <- c(42, 44, 58, 48, 61, 62, 47, 52, 46, 50,
+        43, 40, 56, 49, 54, 39, 42, 47, 42, 54,
+        37, 62, 59, 37, 58, 27, 56, 50, 26, 50)
> g <- rep(factor(1:k), each=n)
> 
> summary(aov(x ~ g))
            Df Sum Sq Mean Sq F value Pr(>F)
g            2  141.9   70.93   0.786  0.466
Residuals   27 2436.0   90.22               
> 
> #######
> 
> ( between.var <- var(by(x, g, mean)) )
[1] 7.093333
> ( within.var <- 16)
[1] 16
> power.anova.test(groups=groups, n=n, between.var=between.var, within.var=within.var)

     Balanced one-way analysis of variance power calculation 

         groups = 3
              n = 10
    between.var = 7.093333　　★★ MSb の 1/n
     within.var = 16
      sig.level = 0.05
          power = 0.713057

NOTE: n is number in each group

> 
> ( f <- sqrt(var(by(x, g, mean))*(k-1)/k / within.var) )
[1] 0.5436502
> pwr.anova.test(k=k, n=n, f=f)

     Balanced one-way analysis of variance power calculation 

              k = 3
              n = 10
              f = 0.5436502
      sig.level = 0.05
          power = 0.713057

NOTE: n is number in each group

> 
> #######
> 
> ( within.var <- 9.022)
[1] 9.022
> power.anova.test(groups=groups, n=n, between.var=between.var, within.var=within.var)

     Balanced one-way analysis of variance power calculation 

         groups = 3
              n = 10
    between.var = 7.093333 　　★★ MSb の 1/n
     within.var = 9.022        ★★ MSw の 1/n を指定
      sig.level = 0.05
          power = 0.9285187

NOTE: n is number in each group

> 
> ( f <- sqrt(var(by(x, g, mean))*(k-1)/k / within.var) )
[1] 0.7239826
> pwr.anova.test(k=k, n=n, f=f)

     Balanced one-way analysis of variance power calculation 

              k = 3
              n = 10
              f = 0.7239826
      sig.level = 0.05
          power = 0.9285187

NOTE: n is number in each group
No.21947　Re: 分散分析のpower analysisについて　　【Kei】　2016/03/25(Fri) 21:57

ありがとうございます。

f = sqrt{ (Σ(μi-μ)^2 / k) / within.var } と定義したとき，
power.anova.testとpwr.anova.testが求めるlambdaが同一であり，結果も同じになることは理解できました。

また，
データ個数をn，グループ数をkとすると
群間平方和が
Σ((群内平均 - 全体平均)^2 * n)
それを群間の自由度　k – 1で割ったものが群間平均平方なので
群間平均平方　= Σ( (群内平均 - 全体平均)^2 * n) / (k - 1 )

ここで
between.var = var (群内平均) = Σ( (群内平均 - 全体平均)^2) )/ (k - 1 )
なので
群間平均平方 = n * between.var　つまり　between.var = 群間平均平方 / n
となるところまではなんとか理解できました。

ここで質問なのですが，
within.var に群内平均平方 / n を指定するのはなぜなのですか。
within.var は見積もられるものなので，仮に群内平均平方 / nを指定するということでしょうか。

あと，もう1つ質問なのですが，
いくつかのサイトで異なるeffect sizeの定義を見つけました。
http://www.unt.edu/rss/class/Jon/ISSS_SC/Module009/isss_m91_onewayanova/node8.html
https://en.wikipedia.org/wiki/Effect_size

Eta^2 = SSb /SSt
f^2 = Eta^2 / (1 - Eta^2) もしくは R^2 / (1- R^2)　など

例えば以下の式で
k <- groups <- 3
n <- 10
x <- c(42, 44, 58, 48, 61, 62, 47, 52, 46, 50,43, 40, 56, 49, 54, 39, 42, 47, 42, 54,37, 62, 59, 37, 58, 27, 56, 50, 26, 50)
g <- rep(factor(1:k), each=n)
summary(aov(x ~ g))
の結果から

SSb <- 141.9
SSt <- sum((x - mean(x))^2)
Eta.squared <- SSb / SSt
f <- sqrt(Eta.squared / (1 - Eta.squared))

ここから求まったfは
within.var <- 9.022　
f <- sqrt(var(by(x, g, mean))*(k-1)/k / within.var)
で求めたfとかなり異なります。
f <- sqrt(Eta.squared / (1 - Eta.squared))
で求めたfはpwr.anova.testに使用するべきではないということでしょうか。

お忙しいところ恐縮ですがまたご指導下さいますようお願い申し上げます。
No.21949　Re: 分散分析のpower analysisについて　　【青木繁伸】　2016/03/26(Sat) 05:04

> within.var に群内平均平方 / n を指定するのはなぜなのですか。
> within.var は見積もられるものなので，仮に群内平均平方 / nを指定するということでしょうか。

「within.var = MSe / n と指定してやると，データに基づく power の予測になるでしょう」と書いておきましたが。

> f <- sqrt(Eta.squared / (1 - Eta.squared))
> で求めたfはpwr.anova.testに使用するべきではないということでしょうか。

要するに，oneway ANOVA の power を計算するには，非心 F 分布が必要で，その非心度（power.anova.test, pwr.anova.test では lambda）をどのように設定するかということです。

そ れを計算するための統計量は between.var, within.var を使う（power.anova.test），および，それを使って計算される f を使う（pwr.anova.test）以外にも，sqrt(Eta.squared / (1 - Eta.squared)) でも，F 値でも何でも計算の元にしてよいのだけれど，非心度として計算される lambda は同じでなければならない。逆に言えば，同じ lambda が得られように計算するなら，どのようなものを与えてもよいということ。pwr.anova.test はそのような要望を満たすために作られた。そして，power.anova.test と同じ答えが出る。

sqrt(Eta.squared / (1 - Eta.squared)) を引数に与えて power.anova.test と同じ答えが出るように lambda を計算して答えを出すプログラム（lambda を計算する部分だけが異なるプログラム）を書いて，それを pwr2.anova.test なりの名前で公表すればよいだけの話（自分で使うだけでも良いのだけど）。
No.21950　Re: 分散分析のpower analysisについて　　【Kei】　2016/03/27(Sun) 08:48

青木先生

ご指導ありがとうございます。
データに基づく power の予測にはwithin.var = MSe / n と指定するとよい点は理解しました。
今回の研究に関してはこれでパワーアナリシスを行うようにします。

また，
f <- sqrt(Eta.squared / (1 - Eta.squared))
をそのままpwr.anova.testに代入してしまうとlambdaの値が正しく求められないので適切でないということは理解できました。
勉強のためにf <- sqrt(Eta.squared / (1 - Eta.squared))から正しいlambdaを求めてパワーアナリシスを行うようなプログラムも書いてみるようにします。

● 「統計学関連なんでもあり」の過去ログ--- 048 の目次へジャンプ
● 「統計学関連なんでもあり」の目次へジャンプ
● 直前のページへ戻る