ねこすたっと

ねこの気持ちと統計について悩む筆者の備忘録的ページ。

サンプルサイズ計算:1群の連続アウトカム(正規分布)の平均を既知の値と比較する [R]

自分用のリファレンスとして、サンプルサイズ計算に関する記事をシリーズで書いています。 なるべく体裁を統一するために、以下のように決めています。

  • 群を示す添字について:
    • c, C = 対照群
    • t, T = 介入群*1
    • e, E = 要因群
    • 添字なし = 全体
  • 使用する記号:
    • p = 反応割合
    • μ = 平均値
    • σ = 標準偏差
    • Φ = 割り付け比 Nt/Nc(デフォルトは1)
    • α = αエラー(デフォルトは0.05)
    • β = βエラー(デフォルトは0.2)
    • z = 標準正規分布の累積分布関数

想定シナリオ

疾患Dに対する治療として、従来薬Cと新規薬Tのどちらが優れているだろうか?
従来薬Cについては既に多くの研究が行われており、治療後の疾患活動度スコアが60点となることが分かっている。 そこで今回の研究では、患者全員に新規薬Tを投与して、従来薬Cの疾患活動度スコア 60点と比較したい。 またこれまでの知見から、新規薬Tにおける疾患活動度スコアは50点となることが見込まれる。 両側有意水準 5%, 検出力 80%として、必要なサンプルサイズはどれくらいだろうか?

このシナリオでサンプルサイズを計算するのに必要な条件を整理してみると、以下のようになる。

  • 帰無仮説  H_0: \mu_T = \mu_C
  • 対立仮説  H_A: \mu_T \neq \mu_C
  • 有意水準:α = 0.05
  • 検出力:1-β = 0.80
  • 従来薬C群における治療後活動度スコアの平均値(既知の値): \mu_C = 60
  • 新規薬T群における治療後活動度スコアの平均値(想定値): \hat{\mu}_T = 50
  • 新規薬T群における活動度スコアの標準偏差(想定値): \hat{\sigma} = 10

新規薬Tの効果は、2群の平均値の差  \delta = \mu_T - \mu_S で推定され、効果サイズと呼ばれる。 効果サイズの大小は、スコアのバラツキ具合、つまり標準偏差  \sigma(standard deviation, SD)と比べて相対的に判断する必要があるので、

 
\begin{aligned}
\Delta = \frac{\delta}{\sigma}
\end{aligned}

で計算される「標準効果サイズ(別名:Cohen's d) 」が用いられる。
標準効果サイズは現実的には0.1-1.0の範囲で設定され、目安は以下のとおり。

  • Δ=0.2:効果が小さい
  • Δ=0.5:効果は中等度
  • Δ=0.8:効果が大きい

方法0:手計算

なるべく正確に計算しようとすると、対立仮説のもとでの効果サイズの分布が非心t分布(non-central t [NCT] distribution)*2に従うことをもとに計算しなければならないが、自由度(平たく言えばサンプルサイズ)が増えるほど正規分布に近づく。

非心t分布を正規分布に近似するの条件の下で、サンプルサイズを計算するには以下の式を使う。

 
\begin{aligned}
N = \frac{(z_{1-\alpha/2} + z_{1-\beta})^2}{\Delta ^2} + \frac{z_{1-\alpha/2}^2}{2}
\end{aligned}

上の式に、

  •  \Delta = \frac{(\mu_T - \mu_C)}{\sigma} = 10/10 = 1.0
  •  z_{1-\alpha/2} = 1.96
  •  z_{ 1 - \beta} = 0.8416

を代入すると、N = 9.76...となるので、サンプルサイズは10例必要なことが分かる。

方法1:pwr.t.test( )を使って計算する

まずはpwrパッケージの関数を使ってみる。この関数でサンプルサイズを計算するときは以下の引数を指定する。

  • d:Cohenのd(標準化効果サイズ)。ここでは(60-50)/10=1を指定。
  • sig.level:有意水準
  • power:検出力
  • type:"two.sample"(2群比較)、"one.sample"(1群の平均を既知の値と比較)、"paired"(対応のある比較)から選ぶ。ここでは"one.sample"を指定。
  • alternative:対立仮説の種類。ここでは"two.sided"(両側検定)を指定。
library(pwr)
pwr.t.test(d = (60-50)/10,
           type = "one.sample",alternative="two.sided",
           sig.level=0.05,
           power=0.8)
     One-sample t test power calculation 

              n = 9.93785
              d = 1
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

10例必要という結果になった。

方法2:TwoSampleMean.Equality( )を使って計算する

次にTrialSizeパッケージの関数を使ってみる。
この関数では対立仮説(alternative hypothesis)は関数名に含まれているので、引数として指定する必要はない。*3

指定する引数は以下のとおり。

  • alpha:αエラー
  • beta:βエラー
  • sigma:標準偏差
  • margin:両群の平均の差。標準化していない効果サイズ。ここでは60-50=10を指定。
library(TrialSize)
OneSampleMean.Equality(alpha = 0.05,
                       beta = 0.2,
                       sigma = 10,
                       margin = 10)
[1] 7.84888

必要数は8例という結果になった。

方法3:power.prop.test( )を使って計算する

このシナリオなら基本関数を使っても計算できる。
サンプルサイズ計算のときに指定する主な引数は以下のとおり。

  • delta:両群の平均の差
  • sd:標準偏差
  • sig.level:αエラー
  • power:1-βエラー
  • type:"two.sample"(2群比較)、"one.sample"(1群の平均を既知の値と比較)、"paired"(対応のある比較)から選ぶ。ここでは"one.sample"を指定。
  • alternative:対立仮説。"one.sided"(片側)か"two.sided"(両側)を指定。
power.t.test(delta = 10,
             sd = 10,
             sig.level = 0.05,
             power = 0.8,
             type = "one.sample",
             alternative = "two.sided")
     One-sample t test power calculation 

              n = 9.937864
          delta = 10
             sd = 10
      sig.level = 0.05
          power = 0.8
    alternative = two.sided

10例必要という結果になった。

おわりに

  • 家を一泊空けてから、ネコにご飯以外で呼ばれることが増えた気がします

参考資料

  • 奥村先生の効果サイズに関する記事です。

oku.edu.mie-u.ac.jp

  • サンプルサイズ計算の本なのにあまり眠くならない。手計算はこの本を参考にしました。

*1:interventionのIを使うと1と区別しにくいのでtrial

*2:t分布を一般化したもの。一般的には中心が0ではなく、左右対称でない。

*3:"equality"は"equivalent"(同等性)と混同しやすいので注意。この研究の仮説は「優越性(superiority)」なんだけど、このパッケージで"superiority"は非劣性試験と対応させた「十分な優越性(substantially superiority)」の意味で使われているので、これも注意。