統計ER

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

Bonferroni 型 多重比較のサンプルサイズ計算 - 統計ソフトRのスクリプト

ブログランキングに参加しています。
まずはぽちぽちっとお願いします。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 数学へ

三群以上の平均値や割合の差を比較したいデザインの場合、サンプルサイズはどのようにすればよいか?

Bonferroni (ボンフェローニ)型の比較は、保守的過ぎて有意水準が厳しすぎるとは言われるが、逆にBonferroniの条件を満たすサンプルサイズなら、HolmやHochbergの検定もクリアできるはず。

Bonferroni型多重比較のサンプルサイズ計算を統計ソフトRで行ってみた。

Bonferroni型多重比較とは?

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

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

詳しくは過去記事を参照。より現実的で適切な方法として、HolmやHochbergの方法も紹介している。

toukeier.hatenablog.com

toukeier.hatenablog.com

平均値の差の検定のサンプルサイズ計算

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

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

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

toukeier.hatenablog.com

検出力を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

平均値の差の検定 Bonferroni型多重比較

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

最低限、標準化した差 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()は非心t分布を使った本格的な計算だったはずで、上記のスクリプトのような近似の方法ではないので、若干のずれがある。とは言うものの実用には問題ない範囲のずれと思う。

> 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は、(4x3)/(2x1)=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型 平均値の多重比較 サンプルサイズ計算 エクセルファイル(2020年5月21日追記)

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

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

動画による使い方解説(20205月27日追記)

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

youtu.be

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

割合の差の検定はどうか?まず二群のサンプルサイズ計算から。群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

割合の差の検定 Bonferroni型多重比較

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

群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")

こちらも参照のこと。

toukeier.hatenablog.com

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値調整は実際には使用しないが、今回のサンプルサイズ計算は保守的で厳しめな方法と言える。