R で多重比較に必要となるサンプルサイズを計算する方法

多重比較のサンプルサイズ計算を R で行う方法

ボンフェローニ型 多重比較のサンプルサイズ計算

Bonferroni（ボンフェローニ）型多重比較とは、比較する数で有意水準を割って、割った有意水準より小さい有意確率の場合、統計学的有意と考える方法。

三群あって、総当たりで三つの比較をする場合、有意水準たとえば0.05を3で割って、一つの比較に対して有意水準を0.0167とする方法。

詳しくは過去記事を参照。

より現実的で適切な方法として、HolmやHochbergの方法も紹介している。

多重比較の前に平均値の差の検定のサンプルサイズ計算

多重比較の前に二群の比較のサンプルサイズ計算を実行してみる。

統計ソフトRでは、power.t.test()で計算できる。検出力 powerと差（SDを指定しない場合標準化した差）delta を指定すると計算できる。

deltaについてはこちらも参照のこと。

検出力を0.8（80%）、標準化した差を0.8（大き目）とする。

power.t.test(power=0.8, delta=0.8)

計算結果は、一群25.5、（人数なので）切り上げると26例必要と計算された。

> power.t.test(power=0.8, delta=0.8)
Two-sample t test power calculation
n = 25.52463
delta = 0.8
sd = 1
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group

平均値の差を多重比較するのに必要なサンプルサイズの計算

平均値の差の検定を三群以上で行う場合のサンプルサイズ計算のスクリプトは以下の通り。

最低限、標準化した差 delta と比較ペア数 num.tests は決めないといけない。SDを1とし、有意水準0.05、検出力0.8、両側検定 two.sidedをデフォルトにしてある。

samplesize.bonferroni.mean <- function(sig.level=.05, power=.8, delta, sd=1, alternative=c("two.sided","one.sided"), num.tests)
{
alternative <- match.arg(alternative)
tside <- switch(alternative, one.sided=1, two.sided=2)
d <- delta/sd
Za <- qnorm(sig.level/num.tests/tside, lower.tail=FALSE)
Zb <- qnorm(power)
n <- 2*((Za+Zb)/d)^2
NOTE <- "n is number in *each* group"
METHOD <- "Sample size for Bonferroni-type multiple comparison of the mean"
structure(list(n = n, "Numbers of tests"=num.tests, "Difference of a pair" = delta,
SD = sd, sig.level = sig.level, "sig.level per test" = sig.level/num.tests,
power = power, alternative = alternative, note = NOTE,
method = METHOD), class = "power.htest")
}

二群比較で検算

二群比較でpower.t.test()の結果と一致するか確認してみる。delta=0.8、比較ペアは1つ num.tests=1だ。

# Two groups
samplesize.bonferroni.mean(delta=0.8, num.tests=1)

結果は、一群25例必要と計算された。power.t.test()は非心ｔ分布を使った本格的な計算だったはずで、上記のスクリプトのような近似の方法ではないので、若干のずれがある。とは言うものの実用には問題ない範囲のずれと思う。

> samplesize.bonferroni.mean(delta=0.8, num.tests=1)
Sample size for Bonferroni-type multiple comparison of the mean
n = 24.52775
Numbers of tests = 1
Difference of a pair = 0.8
SD = 1
sig.level = 0.05
sig.level per test = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group

三群比較（三通りの組み合わせ）

本題の三群比較のサンプルサイズ計算を実行してみる。delta=0.8は同じにして、num.tests=3とする。

deltaは三つの比較の平均的な差を意味している。それが0.8（大き目）と考えて計算するという意味になる。

# Three groups
samplesize.bonferroni.mean(delta=0.8, num.tests=3)

計算結果は、一群33例と計算された。二群比較よりも必要な人数が増える。かなり厳しい計算方法であると感じる。

> samplesize.bonferroni.mean(delta=0.8, num.tests=3)
Sample size for Bonferroni-type multiple comparison of the mean
n = 32.71598
Numbers of tests = 3
Difference of a pair = 0.8
SD = 1
sig.level = 0.05
sig.level per test = 0.01666667
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group

四群比較（六通りの組み合わせ）

四群はどうか？四群になると4つの中から2つを選んで組み合わせるので、コンビネーション4の2は、(4×3)/(2×1)=6通り。

# Four groups
samplesize.bonferroni.mean(delta=0.8, num.tests=6)

結果は、一群38例必要と計算された。

> samplesize.bonferroni.mean(delta=0.8, num.tests=6)
Sample size for Bonferroni-type multiple comparison of the mean
n = 37.84236
Numbers of tests = 6
Difference of a pair = 0.8
SD = 1
sig.level = 0.05
sig.level per test = 0.008333333
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group

平均値の差を多重比較するのに必要なサンプルサイズを計算できるエクセルファイル

エクセルファイルでBonferroni型 平均値の多重比較のためのサンプルサイズ計算ができるようにした。

よければどうぞ。

ボンフェローニ型 平均値の多重比較 サンプルサイズ計算【エクセルでサンプルサイズ】 | TKER SHOP

動画による使い方解説

よければこちらもどうぞ。

割合の差の検定のサンプルサイズ計算

割合の差の検定はどうか？

まず二群のサンプルサイズ計算から。

群1を0.8、群2を0.5として、検出力を0.8とする。

power.prop.test(p1=0.8, p2=0.5, power=0.8)

結果は、一群39例必要と計算された。

> power.prop.test(p1=0.8, p2=0.5, power=0.8)
Two-sample comparison of proportions power calculation
n = 38.48004
p1 = 0.8
p2 = 0.5
sig.level = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group

割合の差を多重比較するのに必要なサンプルサイズの計算

割合の差の検定を三群以上で行うときのサンプルサイズ計算スクリプトは以下の通り。

群Aと群Bの割合と比較ペア数が必須項目だ。

有意水準0.05、検出力0.8、両側検定はデフォルト値としている。

samplesize.bonferroni.prop <- non.inferior.sample.size <- function(pA, pB, power=.8, sig.level=.05, alternative=c("two.sided","one.sided"), num.tests)
{
alternative <- match.arg(alternative)
tside <- switch(alternative, one.sided = 1, two.sided = 2)
p.bar <- (pA+pB)/2
R <- sqrt(2*(p.bar)*(1-p.bar))
S <- sqrt(pA*(1-pA)+pB*(1-pB))
Za <- qnorm(sig.level/num.tests/tside, lower.tail=FALSE)
Zb <- qnorm(power)
n <- ((Za*R+Zb*S)/abs(pA-pB))^2
n.alt <- ((Za*S+Zb*S)/abs(pA-pB))^2
NOTE <- "n is number in *each* group"
METHOD <- "Sample size for Bonferroni-type multiple comparison of the proportion"
structure(list(n = n, "n calc. by H1" = n.alt, "Numbers of tests"=num.tests,
"Difference of a pair" = abs(pA-pB),
sig.level = sig.level, "sig.level per test" = sig.level/num.tests,
power = power, alternative = alternative, note = NOTE,
method = METHOD), class = "power.htest")
}

二群の割合比較で検算

二群比較で、power.prop.test()の結果と一致するか確認してみる。

群Aが0.8、群Bが0.5、比較ペア数は1だ。

# Two groups
samplesize.bonferroni.prop(pA=0.8, pB=0.5, num.tests=1)

結果は、38.48つまり一群39例必要と計算され、power.prop.test()の結果と一致した。

ちなみにn calc. by H1は、平均割合の分散を使わずに、群ごとの割合の分散を使って、対立仮説寄りにした計算。

ゆえに症例数も緩め（少な目）。

> samplesize.bonferroni.prop(pA=0.8, pB=0.5, num.tests=1)
Sample size for Bonferroni-type multiple comparison of the proportion
n = 38.48004
n calc. by H1 = 35.75601
Numbers of tests = 1
Difference of a pair = 0.3
sig.level = 0.05
sig.level per test = 0.05
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group

三群の割合比較（三通りの組み合わせ）

三群比較は、num.tests=3とする。

# Three groups
samplesize.bonferroni.prop(pA=0.8, pB=0.5, num.tests=3)

一群52例と計算された。

緩めの結果だと48例だ。

> samplesize.bonferroni.prop(pA=0.8, pB=0.5, num.tests=3)
Sample size for Bonferroni-type multiple comparison of the proportion
n = 51.53939
n calc. by H1 = 47.69263
Numbers of tests = 3
Difference of a pair = 0.3
sig.level = 0.05
sig.level per test = 0.01666667
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group

四群の割合比較（六通りの組み合わせ）

四群比較はnum.tests=6とする。

# Four groups
samplesize.bonferroni.prop(pA=0.8, pB=0.5, num.tests=6)

結果は、一群60例（もしくは56例）と計算された。

> samplesize.bonferroni.prop(pA=0.8, pB=0.5, num.tests=6)
Sample size for Bonferroni-type multiple comparison of the proportion
n = 59.72726
n calc. by H1 = 55.16575
Numbers of tests = 6
Difference of a pair = 0.3
sig.level = 0.05
sig.level per test = 0.008333333
power = 0.8
alternative = two.sided
NOTE: n is number in *each* group

R の TrialSize パッケージを使った多重比較に必要なサンプルサイズ計算

R には TrialSize というパッケージがあり、臨床試験のサンプルサイズ計算が行える。

最初の一回だけインストールする。install.packages(“”)でインストールだ。

install.packages("TrialSize")

こちらも参照のこと。

library()で呼び出すと使えるようになる。

library(TrialSize)

三群以上の平均値の多重比較 OneWayANOVA.pairwise

三群以上の平均値の多重比較のためのサンプルサイズ計算は、OneWayANOVA.paierwise()を使う。

alphaが有意水準、betaが1-検出力、tauが比較ペア数、sigmaはSD、marginは差（SDが1なら標準化した差）。

三群比較（tau=3）、四群比較（tau=6）をそれぞれ計算してみる。

OneWayANOVA.pairwise(alpha=0.05, beta=0.2, tau=3, sigma=1, margin=0.8)
OneWayANOVA.pairwise(alpha=0.05, beta=0.2, tau=6, sigma=1, margin=0.8)

三群比較は一群33例と計算され、四群比較は、一群38例と計算された。

上記の結果と一致した。

# 三群比較（三通りの組み合わせ）
> OneWayANOVA.pairwise(alpha=0.05, beta=0.2, tau=3, sigma=1, margin=0.8)
[1] 32.71598
# 四群比較（四通りの組み合わせ）
> OneWayANOVA.pairwise(alpha=0.05, beta=0.2, tau=6, sigma=1, margin=0.8)
[1] 37.84236

三群以上の割合の多重比較 OneWayANOVA.PairwiseComparison

三群以上の割合の多重比較のためのサンプルサイズ計算は、OneWayANOVA.PairwiseComparison()を使う。

p1が「ある群」の割合、p2が「別のある群」の割合、deltaが群間差。

deltaは、p1とp2の差だが、指定が必須。

OneWayANOVA.PairwiseComparison(alpha=0.05, beta=0.2, tau=3, p1=0.8, p2=0.5, delta=0.3)
OneWayANOVA.PairwiseComparison(alpha=0.05, beta=0.2, tau=6, p1=0.8, p2=0.5, delta=0.3)

三群比較は48例、四群比較は56例と計算された。

上記の緩めの計算結果と一致する。

# 三群比較（三通りの組み合わせ）
> OneWayANOVA.PairwiseComparison(alpha=0.05, beta=0.2, tau=3, p1=0.8, p2=0.5, delta=0.3)
[1] 47.69263
# 四群比較（六通りの組み合わせ）
> OneWayANOVA.PairwiseComparison(alpha=0.05, beta=0.2, tau=6, p1=0.8, p2=0.5, delta=0.3)
[1] 55.16575

まとめ

Bonferroni型 p値 多重比較調整のためのサンプルサイズ計算のスクリプトを書いてみた。

Bonferroniのp値調整は厳しめな多重調整法なので、今回のサンプルサイズ計算は保守的でより確実な結果が得られる方法と言える。

R の TrialSize パッケージも利用可能である。

参考になれば。