統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

分散分析 ANOVAのサンプルサイズ計算はどうやる?

Toukei Consul Banner

KH Coder Consul Banner

一元配置分散分析、ANOVA(アノーバ) One-way Analysis of Variance のサンプルサイズ計算はどのようにすればいいか?

統計ソフトRでやってみた。

ANOVAのサンプルサイズ計算(誤差分散がわかっている場合)

群の数をk、
群間で最大の差(母平均の差の最大値)をd、
誤差の母分散の平方根をsigma0

とする。

検出力が80%か90%、
有意水準が5%か1%で、

合計4通りに対応できるスクリプトが以下の通り。

anova.sample.size <- function(k, d, sigma0, alpha=0.05, power=.8){
DELTA <- d^2/(2*sigma0^2)
phiA <- k-1
if (power==.8 & alpha==.05) {lambda <- 4.860+3.584*sqrt(phiA)}
else if (power==.9 & alpha==.05) {lambda <- 7.049+4.244*sqrt(phiA)}
else if (power==.8 & alpha==.01) {lambda <- 7.736+4.551*sqrt(phiA)}
else if (power==.9 & alpha==.01) {lambda <- 10.439+5.213*sqrt(phiA)}
n <- lambda/DELTA
METHOD <- "ANOVA Sample Size Calculation"
structure(list(n=n, delta=d, sd=sigma0, sig.level=alpha, power=power, method=METHOD),class="power.htest")
}

群の数が4、最大の群間差が2、誤差分散の平方根を1、 検出力80%、有意水準5%とすると、 一群6例と計算される。

> anova.sample.size(k=4, d=2, sigma0=1)

     ANOVA Sample Size Calculation 

              n = 5.533835
          delta = 2
             sd = 1
      sig.level = 0.05
          power = 0.8

検出力を90%にすると、 一群8例と計算される。

> anova.sample.size(k=4, d=2, sigma0=1, power=.9)

     ANOVA Sample Size Calculation 

              n = 7.199912
          delta = 2
             sd = 1
      sig.level = 0.05
          power = 0.9

ANOVAのサンプルサイズ計算(誤差分散が未知の場合)

通常、誤差分散は未知なので、 こちらを使うべき。

群の数をk、
群間で最大の差(母平均の差の最大値)をd、
誤差の母分散の平方根をsigma0

とするのは前節と同じ。

誤差分散が既知とみなしてサンプルサイズ計算をして、 そのサンプルサイズnで誤差分散未知の場合の検出力を計算し、 検出力が希望通りになるようにサンプルサイズを微調整して、 最終的なサンプルサイズを求めている。

anova.power.michi <- function(k, n, d, sigma0, alpha=.05){
DELTA <- d^2/(2*sigma0^2)
phiA <- k-1
phiE <- k*(n-1)
lambda <- n*DELTA
cA <- (phiA+2*lambda)/(phiA+lambda)
phiA.star <- (phiA+lambda)^2/(phiA+2*lambda)
w <- qf(p=alpha, df1=phiA, df2=phiE, lower.tail=FALSE)
u <- (sqrt(w/phiE)*sqrt(2*phiE-1)-sqrt(cA/phiA)*sqrt(2*phiA.star-1))/(sqrt(cA/phiA+w/phiE))
power <- 1-pnorm(u)
METHOD <- "ANOVA Power Calculation Variance unknown"
structure(list(n=n, delta=d, sd=sigma0, sig.level=alpha, power=power, method=METHOD), class="power.htest")
}

前節の例だと、 群の数が4で、差が2、誤差分散平方根が1だとすると、 一群6例と計算されたので、その通り計算すると、 検出力は80%に到達せず76%になる。

> anova.power.michi(k=4, n=6, d=2, sigma0=1)

     ANOVA Power Calculation Variance unknown 

              n = 6
          delta = 2
             sd = 1
      sig.level = 0.05
          power = 0.7579354

一群に1例ずつ足して、 一群7例とすると検出力はどうなるか?

84%になり、検出力80%を超える。

ゆえに一群7例が適切といえる。

> anova.power.michi(k=4, n=7, d=2, sigma0=1)

     ANOVA Power Calculation Variance unknown 

              n = 7
          delta = 2
             sd = 1
      sig.level = 0.05
          power = 0.8388784

もう一つ前節の例で、 群の数が4、差が2、誤差分散平方根が1、検出力90%のときは、 一群8例と計算される。

> anova.sample.size(k=4, d=2, sigma0=1, power=.9)

     ANOVA Sample Size Calculation 

              n = 7.199912
          delta = 2
             sd = 1
      sig.level = 0.05
          power = 0.9

誤差分散が未知のときのスクリプトで、 一群8例と9例のときの検出力を計算すると、 8例の時は検出力90%をわずかに下回り、9例の時は検出力90%を上回る。

よって、検出力90%で検定を行うためには、 一群9例必要と考える。

> anova.power.michi(k=4, n=c(8,9), d=2, sigma0=1)

     ANOVA Power Calculation Variance unknown 

              n = 8, 9
          delta = 2
             sd = 1
      sig.level = 0.05
          power = 0.8955750, 0.9338735

エクセルでサンプルサイズ計算(2020年4月18日追記)

エクセルで計算できるようにした。よければどうぞ。

https://happyhappygk.base.ec/items/28166707

まとめ(2020年4月18日追記)

一元配置分散分析のサンプルサイズ計算スクリプトを統計ソフトRで書いてみた。

母分散が既知とみなして、予備的にサンプルサイズを計算しておき、母分散未知の計算式で、必要な検出力を確保するように、nを微調整する方法。

微調整は一群当たり1例追加するかどうかという程度。

また、エクセルでも計算できるようにした。よければご利用を。

参考書籍

永田靖著 サンプルサイズの決め方 朝倉書店 pp.129-164
9.3 (誤差分散が既知の場合の)サンプルサイズ設計方法 pp.136-139
10 1元配置分散分析-誤差分散が未知の場合
10.2 検出力の計算方法 pp.150-155
10.3 サンプルサイズの設計方法 pp.156-158

分散分析の実際例に近い値を使ってサンプルサイズ計算をしてみる(2020年5月23日追記)

過去記事で用いたデータを使ってサンプルサイズ計算を行ってみる。

toukeier.hatenablog.com

計算に必要なものは、群の数、群間の最大の差、誤差分散の平方根の3つ。

データは以下の通りとてもシンプル。

A1 <- c(64,72,68,77)
A2 <- c(82,78,77,85)
A3 <- c(55,64,66,49)

群の数は3。群間の最大の差はA2とA3の間で22。

> mean(A1)-mean(A2)
[1] -10.25
> mean(A2)-mean(A3)
[1] 22
> mean(A3)-mean(A1)
[1] -11.75

誤差分散の平方根は、分散分析結果から計算できる「誤差平均平方の平方根」である。

分散分析の結果は以下の通り。

> Anova(lm1)
Anova Table (Type II tests)

Response: X.all
          Sum Sq Df F value   Pr(>F)   
A.all     969.50  2  13.517 0.001945 **
Residuals 322.75  9                    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

誤差平均平方の平方根 sqrt(Sum Sq of Residuals / Df) は以下のように約6になる。

> sqrt(322.75/9) 
[1] 5.988415

群の数が3、最大の群間差が22、誤差分散の平方根は6として、母分散既知のサンプルサイズは以下のように一群2となる。

> anova.sample.size(k=3, d=22, sigma0=6)

     ANOVA Sample Size Calculation 

              n = 1.476973
          delta = 22
             sd = 6
      sig.level = 0.05
          power = 0.8

母分散未知とすると、一群2の場合は検出力が50%を下回る。一群3とすると検出力80%を上回る。実際のデータと同じ一群4とすると検出力90%を上回るという結果になった。

> anova.power.michi(k=3, n=2, d=22, sigma0=6)

     ANOVA Power Calculation Variance unknown 

              n = 2
          delta = 22
             sd = 6
      sig.level = 0.05
          power = 0.4644863

> anova.power.michi(k=3, n=3, d=22, sigma0=6)

     ANOVA Power Calculation Variance unknown 

              n = 3
          delta = 22
             sd = 6
      sig.level = 0.05
          power = 0.8729045

> anova.power.michi(k=3, n=4, d=22, sigma0=6)

     ANOVA Power Calculation Variance unknown 

              n = 4
          delta = 22
             sd = 6
      sig.level = 0.05
          power = 0.9794157

つまり、今回と同じ傾向のデータの場合、一群3で統計学的有意を検出することができる。実際のデータは一群4例だったので、問題なしと言える。事実p値は0.001945で、有意水準5%とすれば統計学的有意であった。

以下再掲。

> Anova(lm1)
Anova Table (Type II tests)

Response: X.all
          Sum Sq Df F value   Pr(>F)   
A.all     969.50  2  13.517 0.001945 **
Residuals 322.75  9                    
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

エクセルでサンプルサイズ使い方解説動画(2020年5月25日追記)

youtu.be