R で分散分析に必要なサンプル数を計算する方法

2018年7月27日2024年10月13日

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

＞＞もう統計で悩むのは終わりにしませんか？

↑期間・数量限定で無料プレゼント中！

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

群の数を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万人以上の医療従事者が購読中

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

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

あわせて読みたい

R で一元配置分散分析を行う方法 Rで、一元配置分散分析を step by step で計算してみた。 lm() と Anova() を使えばあっという間だが、具体的な一つ一つの計算を自分で組み立ててみるとどうか？教科書…

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

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

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

群の数は３。

群間の最大の差は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) は以下のように約６になる。

> sqrt(322.75/9)
[1] 5.988415

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

> 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

母分散未知とすると、一群２の場合は検出力が50%を下回る。

一群３とすると検出力80%を上回る。

実際のデータと同じ一群４とすると検出力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

つまり、今回と同じ傾向のデータの場合、一群３で統計学的有意を検出することができる。

実際のデータは一群４例だったので、問題なしと言える。

事実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