統計ER

R, EZR, SPSS, KH Coder を使ったデータ分析方法を紹介するブログ。ニッチな内容が多め

分散分析のサンプルサイズ、サンプル数、n数の計算

分散分析のサンプルサイズ、サンプル数、n 数の計算を R でやってみた。

>>もう統計で悩むのを終わりにしませんか?


↑1万人以上の医療従事者が購読中

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

群の数を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

実際には、母分散はわかっていることはないため、過去の文献などから仮定して、以下の「誤差分散がわからない場合」の手順を用いるほうがよい。

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

通常、誤差分散は未知なので、こちらを使うほうがよい。

群の数を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

>>もう統計で悩むのを終わりにしませんか?


↑1万人以上の医療従事者が購読中

分散分析のサンプルサイズを別の例で計算してみる

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

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

分散分析のサンプルサイズ計算エクセルファイル

エクセルで計算できるようにした。

よければどうぞ。

分散分析 ANOVA サンプルサイズ計算【エクセルでサンプルサイズ】 | TKER SHOP

分散分析のサンプルサイズ計算エクセルファイル使い方解説動画

youtu.be

まとめ

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

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

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

また、エクセルでも計算できるようにした。

よければご利用を。

参考書籍

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