統計ER

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

t 検定のサンプル数計算

t 検定のサンプル数計算の方法

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


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

t 検定のサンプル数計算(1:1の場合)

R の power.t.test()を用いる。

必要な情報は、検出力をどうするかと、二群の平均値の差。

標準偏差で割った標準化された差が必要だ。

検出力は慣例で80%=0.8が多い。

90%=0.9もよくある。

標準化された差dが難しい。

先行研究は二群の差の標準偏差が計算されていないことがほとんどだ。

そんなときは、以下の目安を使う。

  • 小さな差を検出したければ d=0.1~0.2
  • 中くらいな差を検出したければ d=0.4~0.5
  • 大きな差を検出したければ d=0.8~0.9

検出力80%、中くらいな差を検出したいとして d=0.5 とする。

power.t.test()では、検出力をpower=で、標準化された差をdelta=で、指定する。

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

結果は一群64例必要と計算される。

二群合わせて128例が必要なサンプルサイズだ。

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

     Two-sample t test power calculation 

              n = 63.76576
          delta = 0.5
             sd = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

t 検定のサンプル数計算(1:1の場合)power.t.test() を使わない方法

power.t.test() を使わない簡易的な方法の R スクリプトを下記参考書籍を参照して自作した。

t.test.samplesize <- function(sig.level=0.05, power=0.8, 
delta, sd=1, alternative=c("two.sided","one.sided"))
{
 alternative <- match.arg(alternative)
 tside <- switch(alternative, one.sided=1, two.sided=2)

 Za <- qnorm(sig.level/tside, lower.tail=FALSE)
 Zb <- qnorm(power)
 n <- 2*((Za+Zb)/(delta/sd))^2

 NOTE <- "n is number in *each* group"
 METHOD <- "t-test sample size"
 structure(list(n=n,delta=delta, sd=sd, 
 sig.level=sig.level, power=power, 
 alternative=alternative, note=NOTE, method=METHOD),
 class="power.htest")
}

結果は一群63例と計算された。

power.t.test()を使った時は64例と計算された。

簡易的な方法でも大きな違いがない。

> t.test.samplesize(delta=0.5)

     t-test sample size 

              n = 62.79104
          delta = 0.5
             sd = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n is number in *each* group

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


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

t 検定のサンプル数計算(割り付け比対応)

R の標準では、1:2とか1:3に割り付けるサンプルサイズ計算のプログラムは用意されていない。

自前で作ったプログラムが以下の通り。

power.t.test() で1:1のサンプルサイズを計算した後、割り付けの比(ctrl.ratio)で調整している。

power.t.test.unequal <- function(n = NULL, delta = NULL, sd = 1, sig.level = 0.05, power = NULL, 
 type = "two.sample", alternative = "two.sided", strict = FALSE, ctrl.ratio = NULL) {
 res <- power.t.test(n = n, delta = delta, sd = sd, sig.level = sig.level, power = power,
 type = type, alternative = alternative, strict = strict)
 alternative <- match.arg(alternative)
 tside <- switch(alternative, one.sided=1, two.sided=2)
 nA <- 2 * res$n / (ctrl.ratio + 1)
 nB <- 2 * res$n * ctrl.ratio / (ctrl.ratio + 1)
 Za <- qnorm(sig.level/tside, lower.tail=FALSE)
 power.u <- 1 - pnorm(Za - delta/sd * (1 / sqrt(1 / nA + 1 / nB)))
 METHOD <- paste("Two-sample", "t test power calculation")
 structure(list(n.equal = res$n, n.unequal.1 = nA, n.unequal.2 = nB, delta = delta, sd = sd, 
 sig.level = sig.level, power.equal = res$power, power.unequal = power.u, alternative = alternative,
 method = METHOD), class = "power.htest")
}

powerとdeltaは先ほどと同じにして、1:1, 1:2, 1:3の場合を計算してみる。

まず1:1の場合。

1:1は前節と同じように一群64例と計算される。

> power.t.test.unequal(power=0.8, delta=0.5, ctrl.ratio=1)

     Two-sample t test power calculation 

        n.equal = 63.76576
    n.unequal.1 = 63.76576
    n.unequal.2 = 63.76576
          delta = 0.5
             sd = 1
      sig.level = 0.05
    power.equal = 0.8
  power.unequal = 0.8060089
    alternative = two.sided

1:2は43例と86例と計算される。

検出力は76%に落ちる。

> power.t.test.unequal(power=0.8, delta=0.5, ctrl.ratio=2)

     Two-sample t test power calculation 

        n.equal = 63.76576
    n.unequal.1 = 42.51051
    n.unequal.2 = 85.02102
          delta = 0.5
             sd = 1
      sig.level = 0.05
    power.equal = 0.8
  power.unequal = 0.7586038
    alternative = two.sided

1:3は32例と96例と計算される。

検出力はさらに落ちて69%になる。

> power.t.test.unequal(power=0.8, delta=0.5, ctrl.ratio=3)

     Two-sample t test power calculation 

        n.equal = 63.76576
    n.unequal.1 = 31.88288
    n.unequal.2 = 95.64865
          delta = 0.5
             sd = 1
      sig.level = 0.05
    power.equal = 0.8
  power.unequal = 0.6861757
    alternative = two.sided

ここからは1:1が最も検出力が高いことがわかる。

1:1の割り付けがおすすめなのは、検出力が高いからだ。

1:1割り付け出ない場合、最終的な検出力を80%に保つためには最初の検出力を上げておかねばならない。

最初の検出力を85%にしておいて、1:2割り付けのサンプルサイズを計算すると、どうなるか?

> power.t.test.unequal(power=0.85, delta=0.5, ctrl.ratio=2)

     Two-sample t test power calculation 

        n.equal = 72.80053
    n.unequal.1 = 48.53369
    n.unequal.2 = 97.06738
          delta = 0.5
             sd = 1
      sig.level = 0.05
    power.equal = 0.85
  power.unequal = 0.8116907
    alternative = two.sided

結果は、49例と98例で、最終の検出力が81%になる。

これならば、必要症例数は64例から49例に15例減り、検出力は80%を保っているため、有用かつ適切なデザイン変更だ。

ただし、合計の必要例数例数は49+98=147となり、1:1のときの合計例数128例より19例増加する。

1:3の場合は、1:1の検出力を90%にすると、1:3の検出力が80%を超える。

> power.t.test.unequal(power=0.90, delta=0.5, ctrl.ratio=3)

     Two-sample t test power calculation 

        n.equal = 85.03129
    n.unequal.1 = 42.51564
    n.unequal.2 = 127.5469
          delta = 0.5
             sd = 1
      sig.level = 0.05
    power.equal = 0.9
  power.unequal = 0.8060558
    alternative = two.sided

必要例数は43例と128例で、合計171例必要になる。1:1の時より43例増加する。

症例を集めるのが大変な場合は致し方ないが、できれば1:1にしたほうが、合計の必要なサンプルサイズは小さくなる。

t 検定のサンプル数計算(割り付け比対応)power.t.test() を使わない方法

1:1の場合と同様に、power.t.test() を使わない簡易的な方法のスクリプトを作成して、power.t.test() を使った場合と比較した。

power.t.test.unequal.approx <- function(
 delta = NULL, sd = 1, sig.level = 0.05, power = 0.8, 
 alternative = c("two.sided","one.sided"), 
 ctrl.ratio = NULL) {

 alternative <- match.arg(alternative)
 tside <- switch(alternative, one.sided=1, two.sided=2)

 Za0 <- qnorm(sig.level/tside, lower.tail=FALSE)
 Zb0 <- qnorm(power)
 n <- 2*((Za0+Zb0)/(delta/sd))^2

 nA <- 2 * n / (ctrl.ratio + 1)
 nB <- 2 * n * ctrl.ratio / (ctrl.ratio + 1)
 Za <- qnorm(sig.level/tside, lower.tail=FALSE)
 power.u <- 1 - pnorm(Za - delta/sd * (1 / sqrt(1 / nA + 1 / nB)))

 METHOD <- paste("Two-sample", "t test power calculation")
 structure(list(n.equal = n, n.unequal.1 = nA, 
 n.unequal.2 = nB, delta = delta, sd = sd, 
 sig.level = sig.level, power.equal = power, 
 power.unequal = power.u,  alternative = alternative,
 method = METHOD), class = "power.htest")
}

まず1:1の場合。

power.t.test()の結果と比べ、各群1例少ないのは上述した通り。

> power.t.test.unequal.approx(power=0.8, delta=0.5, ctrl.ratio=1)

     Two-sample t test power calculation 

        n.equal = 62.79104
    n.unequal.1 = 62.79104
    n.unequal.2 = 62.79104
          delta = 0.5
             sd = 1
      sig.level = 0.05
    power.equal = 0.8
  power.unequal = 0.8
    alternative = two.sided

1:2の場合は、1:1の場合の検出力を85%にすると、検出力が80%を超えるのは同じ。

> power.t.test.unequal.approx(power=0.85, delta=0.5, ctrl.ratio=2)

     Two-sample t test power calculation 

        n.equal = 71.82718
    n.unequal.1 = 47.88479
    n.unequal.2 = 95.76957
          delta = 0.5
             sd = 1
      sig.level = 0.05
    power.equal = 0.85
  power.unequal = 0.8064989
    alternative = two.sided

小さい例数の群は49から48例に1例減少、大きい例数の群は98例から96例に2例減少した。

大きな違いはないと言える。

1:3の場合は、1:1の場合の検出力を90%にすると、検出力が80%を超えるのはこちらも同様。

> power.t.test.unequal.approx(power=0.90, delta=0.5, ctrl.ratio=3)

     Two-sample t test power calculation 

        n.equal = 84.05938
    n.unequal.1 = 42.02969
    n.unequal.2 = 126.0891
          delta = 0.5
             sd = 1
      sig.level = 0.05
    power.equal = 0.9
  power.unequal = 0.8015779
    alternative = two.sided

小さい例数の群は、43例から変化なし。

大きい例数の群は、128例から127例に1例減少。

こちらも大きな違いはないと言える。

よって、簡易的な(近似的な)方法は、power.t.test() を使った方法と、実用上は差がないと言えるだろう。

t 検定のサンプル数計算(割り付け比対応)power.t.test() を使わないもう少し簡単な方法

EZR の「2群の平均値の比較のためのサンプルサイズの計算」を参考にして、もう少し簡単に割り付け比ありのサンプルサイズ計算ができるように改造した。

Caseの例数 n1 を計算したときに切り上げて、そこに割り付け比をかけるように変更した。

t.test.samplesize.unequal.approx2 <- function(sig.level=0.05, power=0.8, 
delta, sd=1, alternative=c("two.sided","one.sided"),ctrl.ratio=NULL)
{
 alternative <- match.arg(alternative)
 tside <- switch(alternative, one.sided=1, two.sided=2)

 Za <- qnorm(sig.level/tside, lower.tail=FALSE)
 Zb <- qnorm(power)
 n1 <- ceiling((1+1/ctrl.ratio)*((Za+Zb)/(delta/sd))^2)
 n2 <- n1*ctrl.ratio

 NOTE <- "n1 or n2 is number in *each* group"
 METHOD <- "t-test sample size"
 structure(list(n1=n1, n2=n2, "control ratio"=ctrl.ratio, delta=delta, sd=sd, 
 sig.level=sig.level, power=power, 
 alternative=alternative, note=NOTE, method=METHOD),
 class="power.htest")
}

sig.level=0.05, power=0.8, sd=1のデフォルトのままとして、delta=0.5, ctrl.ratio=1,2,3のとき、結果は以下の通り。

> t.test.samplesize.unequal.approx2(delta=0.5,ctrl.ratio=1)

     t-test sample size 

             n1 = 63
             n2 = 63
  control ratio = 1
          delta = 0.5
             sd = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n1 or n2 is number in *each* group

> t.test.samplesize.unequal.approx2(delta=0.5,ctrl.ratio=2)

     t-test sample size 

             n1 = 48
             n2 = 96
  control ratio = 2
          delta = 0.5
             sd = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n1 or n2 is number in *each* group

> t.test.samplesize.unequal.approx2(delta=0.5,ctrl.ratio=3)

     t-test sample size 

             n1 = 42
             n2 = 126
  control ratio = 3
          delta = 0.5
             sd = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

NOTE: n1 or n2 is number in *each* group

t 検定のサンプル数計算をエクセルで

エクセルファイル販売中。

よければどうぞ。

t検定 平均値の差の検定(割り付け比にも対応)サンプルサイズ計算【改良版】【エクセルでサンプルサイズ】 | TKER SHOP

使い方解説動画。

こちらもよければ、ぜひ。

youtu.be

まとめ

t 検定のサンプル数計算を R で行う方法を解説した。

1:1 以外の割り付け比にも対応。

エクセルでも計算可能。

参考になれば。

参考書籍