Understanding Confidence Intervals via Hypothetical Replications in R

95%信頼区間の「95%」は、何を意味しているのでしょうか。仮想的反復をRで体験するcoffee-chat guide。娘が感じた信頼区間の感覚を再現してみてください。

Frequentist Experiments II − Understanding Confidence Intervals via Hypothetical Replications in R

Keywords: probability model, simulation, study design, survival & competing risks

お父さん、その95%信頼区間って本当に95%なの？

私「あのさ、もうひとつ聞きたいと思っていたことがあって。95%信頼区間（confidence interval）って、みんな当たり前のようにいうけど、ほんとに95%も当たってるの？」

お父さん「ん？そりゃ当たってるよ。仮定した確率モデルが正しければね」

私「ふーん。シミュレーションもできる？」

お父さん「ああ、そういうこと。じゃあ続けてハザード比の95%信頼区間の被覆確率（coverage）を出してみようか。Rでね」

私「coverage？つまり、本当に100回に95回当たるかどうか？」

お父さん「そう。信頼区間が真値を当てる頻度を、シミュレーションデータで数えるんだ。頻度論の約束が、どの程度守られているかを確かめる仮想実験だよ。前提条件を整理しよう。2群比較の生存曲線を比べたいんだよね」

私「ストーマあり群とストーマなし群ね」

お父さん「そうだね、ストーマあり群とストーマなし群のOSを比べよう。真のハザード比は1.5、つまりストーマあり群の方が寿命が1.5倍短い。ハザード比の95%信頼区間は、1000回シミュレーションしたら、1.5を約950回は含んでいるはず」

私「ん？私、調査は1回しかしないつもりだけど」

お父さん「あのね、聞いて。統計家が話している確率モデルは、実際にこれから行われる研究そのものじゃないと思った方がいい。たとえるなら”こういう集団から、こういうルールでデータを発生させて”と別世界に指示を出す手紙なんだ。現実と別世界はデータとRとシミュレーションでつながっている。その手紙に書かれたルールに従って別世界の誰かが研究を1000回繰り返す」

私「ふーん。そこで当たりがでた割合がcoverageね」

お父さん「そう。以前JCOG9502の話をしたとき、“p値が0.05だったら100回に5回は正しい”っていってたよ。自分の経験だと違って見えるものだけど同じこと。それでね、生存曲線の形は指数分布で決める。途中で追跡できなかった患者は打ち切り扱いになるけど、それも指数分布で発生させる。この設定だと、Cox回帰の仮定はすべて正しい。だから理論的にcoverageは95%になるはず」

私「じゃあ、その通りになるかどうか、Rで確かめればいいんだね」

coxph()を用いた95%信頼区間のシミュレーション

ハザード比を推定する方法はいくつかありますが、もっともポピュラーなのはsurvivalパッケージのcoxph()です。以前のデモで、generate_data(hr1, hr2)を用いて、ストーマの有無やOSのデータを生成しました。今回は同じデータにCox回帰を当てはめ、ストーマあり群とストーマなし群のハザード比を計算してみます。この関数では、死亡ハザード比の真値はhr2という引数で指定できます。

generate_data <- function(n=200, hr1, hr2) {
  stoma <- rbinom(n, size = 1, prob = 0.4)
  sex <- rbinom(n, size = 1, prob = 0.5)
  age <- rnorm(n, mean = 65 + 3 * stoma, sd = 8)
  hazard_relapse   <- ifelse(stoma == 1, hr1*0.10, 0.10) # 再発のハザード（大きいほど早く起こる）
  hazard_death     <- ifelse(stoma == 1, hr2*0.10, 0.10) # 死亡のハザード（大きいほど早く起こる）
  hazard_censoring <- 0.05                               # 打ち切りハザード（群に依存しない）
  
  t_relapse   <- rexp(n, rate = hazard_relapse)   # 再発までの潜在時間
  t_death     <- rexp(n, rate = hazard_death)     # 死亡までの潜在時間
  t_censoring <- rexp(n, rate = hazard_censoring) # 打ち切りまでの潜在時間
  
  ## --- 全生存期間（OS） ----------------------------------------
  time_os   <- pmin(t_death, t_censoring)
  status_os <- as.integer(t_death <= t_censoring)  # 1 = 死亡, 0 = 打ち切り
  
  ## --- 無再発生存期間（RFS） -----------------------------------
  time_rfs <- pmin(t_relapse, t_death, t_censoring)
  
  status_rfs <- integer(n)
  status_rfs[time_rfs == t_relapse & time_rfs < t_censoring] <- 1 # 再発
  status_rfs[time_rfs == t_death & time_rfs < t_censoring] <- 1   # 死亡
  
  ## --- 累積再発率（CIR） + 競合リスク --------------------------
  time_cir <- pmin(t_relapse, t_death, t_censoring)
  
  status_cir <- integer(n)
  status_cir[time_cir == t_relapse & time_cir < t_censoring] <- 1 # イベント1: 再発
  status_cir[time_cir == t_death & time_cir < t_censoring] <- 2   # イベント2: 競合リスクとしての死亡
  
  ## --- データフレームにまとめる --------------------------------
  dat <- data.frame(
    id         = 1:n,
    sex        = factor(sex, levels = c(0, 1), labels = c("WOMAN", "MAN")),
    age        = age,
    stoma      = factor(stoma, levels = c(0, 1),
                        labels = c("WITHOUT STOMA", "WITH STOMA")),
    time_os    = time_os,
    status_os  = status_os,
    time_rfs   = time_rfs,
    status_rfs = status_rfs,
    time_cir   = time_cir,
    status_cir = status_cir
  )
}
calculate_coverage <- function(model=c("coxph","finegray"), n, hr1, hr2, hr_true) {
  replications=1000
  covered <- numeric(replications)

  for (r in seq_len(replications)) {
    dat <- generate_data(n, hr1, hr2)
    if (identical(model, "coxph")) {
      fit <- coxph(Surv(time_os, status_os) ~ stoma, data = dat)
    } else if (identical(model, "finegray")) {
      dat$fstatus_cir <- factor(dat$status_cir,
                                levels = 0:2,
                                labels = c("censor",
                                           "relapse",
                                           "death"))
      fgdat <- finegray(Surv(time_cir, fstatus_cir) ~ ., data=dat)
      fit <- coxph(Surv(fgstart, fgstop, fgstatus) ~ stoma, weight=fgwt, cluster=id, data=fgdat)
    }
    confidence_interval <- confint(fit)
    covered[r] <- as.integer(confidence_interval[1] <= log(hr_true) && log(hr_true) <= confidence_interval[2])
  }
  coverage <- mean(covered)
  return(coverage)
}
# install.packages("survival") #インストールが必要なら実行
library(survival)
coverage_200 <- calculate_coverage(model="coxph", n=200, hr1=2, hr2=1.5, hr_true=1.5)
print(coverage_200)
[1] 0.938

私「だいたいcoverageが95%になるっていうのは確認した。でも結局、95%信頼区間ってどう理解すればいいの？」

お父さん「教科書的にいうとね。同じようなデータを何度も集めて、毎回95%信頼区間を作ったとき、“真の値を含む区間”の割合が0.95になるように設計された指標なんだ。これを”in the long run”の性能って言ったりする。でも、概念的に理解するよりも、シミュレーションを行った方がわかりやすいと思って」

私「いやわかりにくいだろ。1回の研究で”この区間が当たっている確率が95%“って意味じゃないんだね」

お父さん「そうだね」

私「じゃあ、人数が100人だったらどうなるの？人数が減ったら当たりにくくはならないの？」

お父さん「設計上は当たりにくくなったりはしないよ。95%信頼区間のcoverageはサンプルサイズによらず95%になるように作ってある。n=100、200、400、800のシミュレーション結果をみてみようか」

coverage_100 <- calculate_coverage(model="coxph", n=100, hr1=2, hr2=1.5, hr_true=1.5)
coverage_400 <- calculate_coverage(model="coxph", n=400, hr1=2, hr2=1.5, hr_true=1.5)
coverage_800 <- calculate_coverage(model="coxph", n=800, hr1=2, hr2=1.5, hr_true=1.5)
print(coverage_100)
[1] 0.958
print(coverage_200)
[1] 0.938
print(coverage_400)
[1] 0.953
print(coverage_800)
[1] 0.948

信頼区間アプローチによるサンプルサイズ設計

私「ほんとだ、だいたい95%。じゃあ調査人数は何人でもいいってことだ」

お父さん「ん？どういう意味？」

私「がんサバイバー何人に調査票を送るかってこと。だって信頼区間の性能は人数によらないんでしょ」

お父さん「あ、ごめん、誤解があったね。サンプルサイズが小さくなると、coverageは同じでも、信頼区間幅が広くなる。つまり、同じ95%でも、どれくらい曖昧な推定値で我慢しないといけないかが変わるんだ。今のところ、何人に調査するつもりだった？」

私「100人から調査票戻ってきたら上出来かなって思ってるけど」

お父さん「やっぱり必要サンプルサイズは計算してなかったか。がんを手術した後の復職率を推定したいってことでいいの？それなら100人でギリギリだと思う」

私「いや？前にお父さんにPECOを教えてもらったよね。あれから上司と相談して、ストーマを造設した患者の復職状況を造設しなかった患者と比べることにしたの」

お父さん「そうなんだ。そうすると計算が違ってくる。100人だと足りないはず」

私「へ？なんでそんなことわかるのよ。100人に協力してもらうのってたいへんなのよ」

お父さん「どれだけサンプルサイズが必要なのかについて、統計学的な根拠に基づいて見積もる方法があるんだよ。この表をみてみて。ある書籍から借りてきたものなんだけど（Machin, et al. 2022）」

表1. 割合の推定において指定した 95%信頼区間幅のため必要なサンプルサイズ

割合 \(\pi\) 0.05 0.10 0.15 0.20
0.01 108 43 25 17
0.02 152 52 29 19
0.03 201 61 33 21
0.04 252 72 37 23
0.05 304 83 41 25
0.06 356 95 46 28
0.07 407 107 50 30
0.08 458 118 55 32
0.09 508 130 60 35
0.10 554 139 62 35
0.11 604 153 69 40
0.12 651 164 74 42
0.13 696 175 78 44
0.14 741 186 83 47
0.15 784 196 87 49

表1. 割合の推定において指定した 95%信頼区間幅のため必要なサンプルサイズ

割合 \(\pi\) 0.05 0.10 0.15 0.20
0.16 826 206 92 51
0.17 867 216 96 54
0.18 907 226 100 56
0.19 945 236 104 58
0.20 982 245 108 60
0.25 1150 286 126 70
0.30 1288 320 141 78
0.35 1395 347 152 84
0.40 1472 366 161 89
0.45 1518 377 166 92
0.50 1533 381 167 93

お父さん「たとえばなんだけどね。復職率を調べたいとするでしょ。そうすると復職率をどれくらい精確に推定できたかを表す指標が、95%信頼区間だよね。仮に復職率が50%だったとしてみよう。そうすると、50%±10%程度の推定精度がほしいよね。ここでいう±の後の数字が、95%信頼区間の幅になるんだ」

私「ああ、”50%（95%CI 40～60%）”みたいな感じで論文に書くことになるわけだ。40〜60%の場合、95%信頼区間幅は20%ってことね」

お父さん「そういうこと。表では復職率を\(\pi\)という変数で表しているんだけど、”\(\pi=0.5\)”と”信頼区間幅0.2”に対応する数字をみて。93って書いてあるでしょ。これが必要サンプルサイズ」

私「どうだろう。復職率ってもっと高くない？手術したときの年齢にもよるけど、たとえば患者さんが50歳前後だったら、80%はいってほしいなあ」

お父さん「そうなんだ。そうすると復職率80%は非復職率20%と同じことだから、”\(\pi=0.2\)”と”信頼区間幅0.2”のところをみればいい。サンプルサイズは60人」

私「100人で足りてるじゃん」

お父さん「いやいやいや、これは復職率の信頼区間からサンプルサイズを計算したからこうなったんだ。“信頼区間アプローチ”といって、主に探索的研究で使うものなんだよ。別の計算の仕方もある。ストーマ造設ありとなしの復職率を比べるような仮説検証を目指した研究なら、“仮説検定アプローチ”といって、仮説検定の考え方にそってサンプルサイズを計算しないといけない」

私「そっか、だいたいわかった。私、午後から出かけるから。また今度教えてよ、じゃあね」

サンプルサイズ設計の2つのアプローチ

研究計画を立てるときには、サンプルサイズを見積もる必要があります。サンプルサイズの計算は大まかにいえば2種類あって、信頼区間アプローチ（たとえば表1）と仮説検定アプローチ（次回）を用いることができます。

信頼区間アプローチは、調査や探索的研究で用いられます。このアプローチによって割合を推定するために必要なサンプルサイズを求めるとき、次の数値を設定する必要があります。

  • 割合の真値（復職率など）
  • 信頼区間幅

もうひとつの仮説検定アプローチは、臨床試験や仮説検証型の研究に適しています。このアプローチでは、以下の3つの要素が最低でも必要です。

  • 有意水準（通常5%）
  • 検出力（通常80～90%）
  • 検出したい効果サイズ（2群間の生存確率の差やハザード比など）

サンプルサイズ設計は、βエラーを制御するという意味で、p値と対になる統計手法といえます。

このエピソードに関係するクイズです

臨床検査の正常範囲の意味として正しいのはどれでしょうか。

  1. ある集団からのランダムサンプルにおいて検査値の平均±1.96×標準偏差
  2. 健常人集団における検査値の95%が含まれる範囲
  3. 健常人集団における検査値の平均±1.96×標準誤差
  4. ある集団からのランダムサンプルにおいて検査値の95%が含まれる範囲
  • 正解は2です

標準偏差（SD）と標準誤差（SE）の区別が重要です。前者は「検査値自体が個人間でどのくらいばらつくか」を、後者は「平均の推定がどのくらい精確か」を表しています。そして、「平均±1.96×標準偏差」は、正規分布における「95%が含まれる範囲」に相当しています。

なお、似て非なる範囲に「平均±1.96×標準誤差」があります。これは「仮想的反復の下で真の平均を95%の割合で包む区間」、すなわち頻度論統計学でいう95%信頼区間に対応します。

文献

次のエピソードとRスクリプト