ブートストラップ法とは ― 信頼区間・バイアス補正・臨床試験での活用をRで実装する ―
この記事でわかること
- ブートストラップの原理：母集団からの標本抽出を手元の標本からの復元抽出で置き換えるだけで、分布の仮定なしに標準誤差が出せる理由
- 4種類の信頼区間の使い分け：正規近似法・基本法・パーセンタイル法・BCa法の定義式と長所短所、実務での既定の選び方
- bootパッケージの実装：
boot()とboot.ci()で標準誤差・バイアス・4種類の信頼区間を得る手順と層別ブートストラップ - パーセンタイル法とBCa法で結論が変わる実例：同じデータ・同じ再標本回数なのに「差がある」かどうかの判定が割れる理由
- 被覆確率のシミュレーションで見る限界：小標本では名目95%に届かないこと、そもそも破綻する統計量があること
はじめに
臨床試験の統計解析計画書（SAP）を書いていると、「この推定量の標準誤差はどう計算するのか」で手が止まる場面に必ず出会います。2群の平均値の差なら教科書どおりの公式があります。ところが実務で本当に報告したい量は、中央値、RMST（制限付き平均生存時間）、Win Ratio、ROC曲線下面積（AUC）、発生率の比といった、分散の閉じた公式が存在しない、あるいは存在しても近似の当てにならない統計量ばかりです。
t検定やWald型の信頼区間が寄りかかっているのは、「推定量がおおむね正規分布に従う」という前提でした。標本が十分に大きく分布がそれほど歪んでいなければうまく働きますが、炎症マーカーのように右へ大きく裾を引くデータや数十例規模の初期相試験では、この正規近似はしばしば裏切ります。区間が名目どおりの確率で真値を捕まえてくれないのは、推測統計としてかなり深刻な事態です。
そこで効いてくるのがブートストラップ法です。分布の形についての仮定を置かず、手元のデータだけを使って標準誤差と信頼区間を数値的に作ってしまう手法で、Efron が1979年に提案しました。名前の由来は「自分の靴ひも（bootstraps）を引っ張って自分自身を持ち上げる」という、ありえない離れ業を指す英語の慣用句です。たった1組の標本から、試験を何度も繰り返したかのような情報を引き出す、その少し反則めいた感触がそのまま名前になっています。
この記事は、製薬企業やCROで解析を担当されている方、統計検定準1級で再標本法を整理したい方に向けて書いています。前半では原理と4種類の信頼区間の作り方を、後半ではRの boot パッケージによる実装と限界を扱います。信頼区間そのものの考え方は 区間推定入門：数式と図解で理解する信頼区間の世界 を、正規近似が成り立つ背景は 中心極限定理とは を先にご覧いただくと、この先がぐっと読みやすくなります。
ブートストラップ法とは
考え方の核は「母集団を標本で置き換える」の一点
統計的推測の基本構図を思い出してください。背後に母集団分布 \( F \) があり、そこから \( n \) 個の観測値からなる標本が得られ、その標本から統計量 \( \hat{\theta} \) を計算します。同じ試験を何度も繰り返して標本を何組も手に入れられるなら、\( \hat{\theta} \) がどれくらいばらつくか（標本分布）を直接観察できます。しかし現実に手元にあるのは、たった1組の標本だけです。
ブートストラップは、ここで一点だけ大胆な置き換えを行います。\( F \) が分からないのなら、手元の標本そのものを \( F \) の代わりに使ってしまえばよい、というものです。観測値 \( x_1,\dots,x_n \) それぞれに確率 \( 1/n \) を与えた分布を経験分布 \( \hat{F}_n \) と呼びます。\( F \) から標本を取り直すことはできませんが、\( \hat{F}_n \) から取り直すことなら計算機でいくらでもできます。それが「元のデータから \( n \) 個を復元抽出する」という操作です。つまり、\( F \) からの標本抽出で生じるばらつきを、\( \hat{F}_n \) からの再標本抽出で生じるばらつきで近似する、という対応関係です。
プラグイン原理
知りたいのは母集団の量 \( \theta=T(F) \)、手元にあるのはその経験分布版 \( \hat{\theta}=T(\hat{F}_n) \) です。推測の本体は誤差 \( \hat{\theta}-\theta \) がどう散らばるかですが、これは \( F \) に依存するので分かりません。ブートストラップはここを次のように置き換えます。
\[ \mathcal{L}(\hat{\theta}-\theta\mid F)\ \approx\ \mathcal{L}(\hat{\theta}^{*}-\hat{\theta}\mid\hat{F}_n) \]
ここで \( \mathcal{L}(\cdot) \) は「〜の分布」を、\( \hat{\theta}^{*} \) は再標本から計算した統計量を表します。左辺は観測できませんが、右辺は「自分で作れる母集団」のもとでの誤差の分布なので、再標本を何度でも生成して数値的に評価できます。この差し替えをプラグイン原理と呼びます。
アルゴリズムは4ステップ
手順そのものは拍子抜けするほど単純です。ステップ1として、観測された \( n \) 個のデータから復元抽出で \( n \) 個を選び直し、再標本（ブートストラップ標本）を作ります。同じ観測値が2回3回と選ばれることもあれば、1回も選ばれない観測値も出ます。ステップ2として、その再標本に、元データとまったく同じ手順で統計量を計算します。ステップ3として、これを \( B \) 回（通常1,000〜10,000回）繰り返し、\( \hat{\theta}^{*1},\dots,\hat{\theta}^{*B} \) を得ます。ステップ4として、この \( B \) 個の散らばり方をそのまま推定量のばらつきの推定に使います。標準偏差が標準誤差の推定値、分位点が信頼区間です。
標準誤差とバイアスの定義式
標準誤差の推定値は、\( B \) 個の再標本統計量の標本標準偏差そのものです。
\[ \widehat{\mathrm{se}}_B=\sqrt{\frac{1}{B-1}\sum_{b=1}^{B}(\hat{\theta}^{*b}-\bar{\theta}^{*})^{2}},\qquad \bar{\theta}^{*}=\frac{1}{B}\sum_{b=1}^{B}\hat{\theta}^{*b} \]
ここで \( \hat{\theta}^{*b} \) は \( b \) 番目の再標本から計算した統計量、\( \bar{\theta}^{*} \) はその平均です。分散の公式を導出することなく、「何度も計算してばらつきを見る」だけで標準誤差が得られる点が、実用上の値打ちです。
バイアス（偏り）の推定値は、再標本統計量の平均が観測された推定値からどちらへずれているかで測ります。
\[ \widehat{\mathrm{bias}}_B=\bar{\theta}^{*}-\hat{\theta} \]
この値が正なら、再標本での推定値は観測値より系統的に大きめに出ているということです。平均値のような不偏推定量ではほぼ0ですが、中央値・分散・相関係数・比の推定量では無視できない大きさになり、後で見るBCa法はまさにこのずれを補正するために設計されています。
なぜこの置き換えでうまくいくのか
根拠は素朴です。大数の法則により、経験分布 \( \hat{F}_n \) は \( n \) が増えるにつれて真の分布 \( F \) に近づきます。そして統計量 \( T \) が分布の「滑らかな汎関数」、つまりデータをほんの少し動かしたとき出力もほんの少ししか動かない性質を持っていれば、\( \hat{F}_n \) のまわりでの \( T \) の振る舞いは \( F \) のまわりでの振る舞いをよく真似ます。その結果、ブートストラップ分布は真の標本分布に収束します（一致性）。裏を返せばこの「滑らかであること」が崩れる統計量では話が成り立たなくなるわけで、ここが後半で扱う破綻条件の伏線になります。
中心極限定理との関係
中心極限定理は「\( n \) を大きくすれば標本平均の分布は正規分布に近づく」という漸近論で、近似の良し悪しは \( n \) 次第です。対してブートストラップは、その漸近近似が当てにならない有限標本の局面で、標本分布そのものを数値的に組み立ててしまうアプローチです。両者は対立せず、正規近似が使えるところでは古典的な区間を、怪しいところではブートストラップを、という補完関係にあります。
ブートストラップは、無からデータを生み出す魔法ではありません。増えているのは情報量ではなく計算回数だけで、再標本はあくまで手元の \( n \) 個の観測値の組み合わせにすぎません。その標本が母集団を代表していなければ、再標本も同じ偏りを忠実に引き継ぎます。選択バイアスや測定バイアスは、何万回再標本を取ろうと一切救われません。ブートストラップが対処するのは「標本抽出に伴う偶然のばらつき」だけだ、とお考えください。
標準誤差という量そのものの意味は 【やさしく解説】標準偏差（SD）と標準誤差（SE）の違い が、正規性の仮定に依拠した古典的な区間の導出は t検定の数理的導出とRによる実装例 がよい対比材料になります。
ブートストラップ信頼区間の4つの作り方
\( B \) 個の再標本統計量が手元に並んだとして、そこから95%信頼区間をどう切り出すかには複数の流儀があります。Rの boot.ci() が返す代表的な4つを順に見ていきましょう。以下、\( \alpha=0.05 \)、\( \hat{\theta}^{*}_{(q)} \) は再標本統計量の \( q \) 分位点とします。
正規近似法（normal）
\[ \hat{\theta}-\widehat{\mathrm{bias}}\pm z_{1-\alpha/2}\,\widehat{\mathrm{se}} \]
形は古典的なWald区間そのままで、標準誤差だけをブートストラップで置き換えたものです。ただし単なるWald区間と違い、推定値からバイアス推定値を引いたうえで区間の中心に据えます。ブートストラップ分布がほぼ左右対称ならこれで十分ですが、歪みがあると当てが外れます。
基本法（basic / pivotal）
\[ \left(2\hat{\theta}-\hat{\theta}^{*}_{(1-\alpha/2)},\ 2\hat{\theta}-\hat{\theta}^{*}_{(\alpha/2)}\right) \]
直感は「ブートストラップ分布を観測値 \( \hat{\theta} \) のまわりで左右に反転させる」ことです。誤差 \( \hat{\theta}^{*}-\hat{\theta} \) の分布を \( \hat{\theta}-\theta \) の分布とみなし、\( \theta \) について解き直すとこの形になります。反転するので、再標本分布の上側分位点が区間の下限に効きます。歪みが強いと、割合や生存時間のように負を取りえない量の区間が負まではみ出すことがあります。
パーセンタイル法（percentile）
再標本分布の2.5%点と97.5%点をそのまま区間の両端に採用する、最も直感的で実装も簡単な方法です。長所は変換不変性にあります。対数尺度で区間を作ってから戻しても元の尺度で直接作った区間と一致するため、比やオッズと相性がよいのです。短所は、推定量にバイアスがあっても素通ししてしまう点で、ずれた分布の分位点をそのまま使う以上、区間全体がずれます。
BCa法（bias-corrected and accelerated）
パーセンタイル法の弱点を2つの補正量で埋めた方法です。ひとつはバイアス補正量 \( z_0 \) で、再標本統計量のうち観測値 \( \hat{\theta} \) を下回るものの割合を標準正規分位点に変換して求めます。もうひとつは加速度定数 \( a \) で、観測値を1つずつ抜いて推定値を計算し直すジャックナイフから分布の歪みを推定します。この2つで、実際に読み取る分位点をパーセンタイル法から非対称にずらします。理論的には二次の精度（second-order accurate）を持ち、被覆誤差が \( O(1/n) \) の速さで小さくなります。パーセンタイル法と基本法は一次の精度（\( O(1/\sqrt{n}) \)）にとどまるため、精度ではBCaが一段上です。
| 手法 | 考え方 | 長所 | 短所 | 使いどころ |
|---|---|---|---|---|
| 正規近似法 （normal） | 標準誤差とバイアスだけを使い、Wald型の対称な区間を作る | 解釈が容易で、標準誤差の報告と整合する | 正規性を仮定するため、歪んだ分布で崩れる | 分布がほぼ対称で、標準誤差を主役に報告したいとき |
| 基本法 （basic） | ブートストラップ分布を観測値のまわりで反転させる | 枢軸量の考え方に忠実で、正規性を仮定しない | 変換不変性がなく、定義域を飛び出すことがある | パーセンタイル法との比較・感度確認の材料として |
| パーセンタイル法 （percentile） | 再標本分布の2.5%点と97.5%点をそのまま使う | 実装が最も簡単で、変換不変性がある | バイアスを補正せず素通しするため区間がずれる | 1回の推定が重くBCaまで回せないとき、探索的解析 |
| BCa法 （bca） | バイアス補正量 \( z_0 \) と加速度定数 \( a \) で分位点をずらす | 二次の精度を持ち、バイアスと歪みの両方を補正する | 計算量が大きく、再標本回数が少ないと不安定 | 既定の選択。中央値・比・歪んだ指標の区間推定 |
実務でどれを使うかという結論を先に述べておくと、BCa法を既定にするのが最も無難です。バイアスと歪みという、ブートストラップを使いたくなる場面でまさに問題になる2つの要因を同時に補正してくれるからです。そのうえで、1回の推定にモデル当てはめが必要で計算コストが厳しい場合はパーセンタイル法に退く、という判断になります。ただしBCaは裾の分位点を非対称にずらして読む以上、再標本回数 \( R \) が小さいと補正後の分位点が数個の再標本値に振り回されます。\( R \) は最低でも2,000、可能であれば10,000は取ってください。
4つの区間は「同じデータから作った別の答え」ではなく、精度の階層が付いた選択肢です。迷ったらBCa、計算コストが許さなければパーセンタイル、標準誤差と対にして報告したいときだけ正規近似、と決めておくと迷いません。そして最も大切なのは、どの手法を使うかを解析計画書に事前に書いておくことです。4つ全部を計算してから都合のよいものを選ぶのは、検定の多重性と同じ構造の問題を招きます。
なお、信頼区間が0をまたぐかどうかで有意性を語る場面は多いのですが、その対応関係は自明ではありません。区間と検定の関係は 信頼区間とp値の関係を図解で理解する と p値とは何か？ を併せてご覧ください。後半では、この「0をまたぐかどうか」が手法の選択で入れ替わる実例を扱います。
Rでブートストラップを実装する
以下のコードはすべて R 4.5.1（boot パッケージ 1.3.31）で実行し、冒頭に set.seed(4869) を置いて乱数を固定していますので、同じコードをそのまま実行すれば本文に貼った出力と同じ結果が再現できます。なお boot は追加インストールが不要で、R本体に同梱されている推奨パッケージです。解析環境へのパッケージ追加に社内承認が必要な現場でも今日から使える点は、実務上かなり効いてきます。
ステップ1：歪んだ臨床データを用意する
想定するのは対照群・実薬群がそれぞれ50例の試験で、評価項目は炎症マーカー（CRP, mg/L を想定）です。CRPのような検査値は多くの患者が健常域に集まる一方で一部が桁違いの高値を示すため、右に大きく裾を引きます。そこで対数正規分布からデータを発生させます。真の中央値は対照群が \( \exp(1.6)=4.95 \)、実薬群が \( \exp(0.9)=2.46 \) です。
set.seed(4869)
n <- 50
ctrl <- rlnorm(n, meanlog = 1.6, sdlog = 0.8)
trt <- rlnorm(n, meanlog = 0.9, sdlog = 0.8)
dat <- data.frame(
group = rep(c("Control", "Treatment"), each = n),
crp = c(ctrl, trt)
)
round(do.call(rbind, tapply(dat$crp, dat$group, function(x)
c(n = length(x), mean = mean(x), median = median(x), sd = sd(x)))), 2)
> n mean median sd
> Control 50 6.46 5.0 6.91
> Treatment 50 3.32 2.7 2.49
両群とも平均が中央値より明確に大きくなっています。対照群は平均 6.46 に対して中央値 5.0、実薬群は平均 3.32 に対して中央値 2.7 で、少数の高値例が平均を引き上げている典型的な形です。しかも対照群は標準偏差 6.91 が平均 6.46 を上回っており、「平均±SD」で区間を作ると下限が負になります。CRPは定義上ゼロ未満になりませんから、正規分布を前提にした区間が現実と噛み合っていないことがこの時点でわかります。
それでも機械的に平均±SDとt検定を当てはめることは可能で、実際よく行われています。しかしそこで比較しているのは「典型的な患者の値」ではなく、外れ値に引きずられた平均です。臨床検査値・バイオマーカー・医療費・在院日数は軒並みこの形になりますので、まず分布を見てから要約統計量を選びたいところです。描画の具体は データ可視化の中級：箱ひげ図・バイオリンプロット・散布図行列 をご覧ください。
ステップ2：ブートストラップを自分で書く
boot を使う前に、ブートストラップが sample() と replicate() だけで原理どおり書けることを確認します。一度自分で書いておくと、パッケージが裏で何をしているのかが見通せます。やることは単純で、対照群の50例から「50例を復元抽出して中央値を計算する」操作を2000回繰り返すだけです。
set.seed(4869)
B <- 2000
boot_med <- replicate(B, median(sample(ctrl, size = length(ctrl), replace = TRUE)))
round(c(observed = median(ctrl),
boot_mean = mean(boot_med),
bias = mean(boot_med) - median(ctrl),
se = sd(boot_med)), 3)
quantile(boot_med, c(0.025, 0.975))
> observed boot_mean bias se
> 4.998 4.814 -0.184 0.854
> 2.5% 97.5%
> 2.914 6.039
観測された中央値は 4.998、2000回の再標本から得た中央値の平均は 4.814、その差であるバイアスは -0.184 です。再標本での中央値は観測値よりわずかに小さい側へ系統的にずれています。そして再標本分布のばらつき、すなわち中央値の標準誤差は 0.854 と推定され、パーセンタイル法による95%信頼区間は (2.914, 6.039) となりました。
ここで立ち止まっていただきたいのは、いま求めた「中央値の標準誤差」には、平均の標準誤差 \( s/\sqrt{n} \) のような閉じた式が（漸近的な近似式を除けば）存在しないという事実です。公式集を探しても出てきません。それが10行足らずのコードで、しかも分布の形に何も仮定を置かずに出てしまう。これがブートストラップの威力です。標準偏差と標準誤差の区別があいまいだとこの有り難みが伝わりにくいので、不安な方は 【やさしく解説】標準偏差（SD）と標準誤差（SE）の違い で整理しておくとよいと思います。
ステップ3：bootパッケージで4種類の信頼区間を作る
原理がわかったところで boot に任せます。ここが最大のつまずきポイントです。boot() に渡す統計量の関数は、必ず function(data, indices) という2引数の形で書かなければなりません。第1引数が元データ、第2引数が「今回の再標本でどの行を選んだか」を表す添字ベクトルで、関数の中では元データをそのまま使わず data[indices] と添字で抜き出した部分を使う、というのがお作法です。
function(data) median(data) と1引数で書くと unused argument というエラーで止まりますので、これはすぐ気づけます。むしろ厄介なのは function(data, indices) median(data) のように2引数で受けておきながら indices を使い忘れる場合です。この書き方はエラーにならず、毎回まったく同じ元データの中央値を返してしまうため、標準誤差が 0 と表示されて初めて誤りに気づくという事故になります。逆に関数の中に自分で sample() を書くのも誤りで、復元抽出は boot() 側が済ませた結果が indices として渡ってきます。「抽出はbootの仕事、計算は自分の仕事」と覚えると間違えません。
library(boot)
med_fun <- function(data, indices) median(data[indices])
set.seed(4869)
boot_out <- boot(data = ctrl, statistic = med_fun, R = 2000)
boot_out
boot.ci(boot_out, type = c("norm", "basic", "perc", "bca"))
> ORDINARY NONPARAMETRIC BOOTSTRAP
>
>
> Call:
> boot(data = ctrl, statistic = med_fun, R = 2000)
>
>
> Bootstrap Statistics :
> original bias std. error
> t1* 4.998115 -0.1873139 0.8458154
>
> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
> Based on 2000 bootstrap replicates
>
> CALL :
> boot.ci(boot.out = boot_out, type = c("norm", "basic", "perc",
> "bca"))
>
> Intervals :
> Level Normal Basic
> 95% ( 3.528, 6.843 ) ( 3.958, 6.942 )
>
> Level Percentile BCa
> 95% ( 3.054, 6.038 ) ( 2.914, 5.908 )
> Calculations and Intervals on Original Scale
4つの手法を並べると次のようになります。区間幅は上限から下限を引いた値です。
| 手法 | 95%信頼区間 | 区間幅 | 特徴 |
|---|---|---|---|
| Normal（正規近似法） | (3.528, 6.843) | 3.315 | 再標本分布の正規性を仮定。4手法で最も幅が広く、歪んだ統計量には勧めにくい |
| Basic（基本法） | (3.958, 6.942) | 2.984 | パーセンタイル区間を観測値の周りに折り返す。幅は同じで位置が上側にずれる |
| Percentile（パーセンタイル法） | (3.054, 6.038) | 2.984 | 再標本分布の2.5%点と97.5%点をそのまま使う。直感的だがバイアスは補正しない |
| BCa（バイアス修正加速法） | (2.914, 5.908) | 2.994 | バイアスと歪みの両方を補正して端点をずらす。実務の既定として使いやすい |
Bootstrap Statistics の行では、観測された中央値（original）が 4.998115、バイアスの推定値が -0.1873139、標準誤差が 0.8458154 です。ステップ2で自分で書いた結果（バイアス -0.184、標準誤差 0.854）とほぼ一致しており、パッケージが特別なことをしているわけではないと確認できます。区間で目を引くのは、正規近似法だけが上振れしていることです。下限が 3.528 とほかの手法より明らかに高く、上限も 6.843 と最も高い。区間幅も 3.315 で最大です。理由ははっきりしていて、バイアスの推定値が負であるために区間の中心が \( \hat{\theta}-\widehat{\mathrm{bias}} \) すなわち 5.185 と観測値 4.998 より上へ置かれ、そこから左右対称に幅を取っているからです。歪んだ再標本分布に対称な区間を当てはめる無理が、そのまま位置のずれとして出た格好です。
一方、パーセンタイル法 (3.054, 6.038) とBCa法 (2.914, 5.908) は観測値 4.998 に対して左右非対称で、どちらも下限までの距離のほうが長くなっています。歪んだデータの信頼区間が観測値を中心に対称にならないのはむしろ自然で、この非対称性を素直に表現できることがブートストラップの美点です。BCa法は再標本分布の偏り（バイアス補正量 \( z_0 \)）と歪み（加速度定数 \( a \)）を踏まえ、読み取る分位点をパーセンタイル法よりさらに下側へずらしています。
ステップ4：反復回数 R をいくつにするか
実務でよく聞かれるのが「Rはいくつにすればよいのか」です。実際に動かして確かめましょう。
set.seed(4869)
for (R in c(200, 1000, 2000, 10000)) {
ci <- boot.ci(boot(ctrl, med_fun, R = R), type = "perc")$percent[4:5]
cat(sprintf("R = %5d 95%% percentile CI = (%.3f, %.3f)\n", R, ci[1], ci[2]))
}
> R = 200 95% percentile CI = (2.917, 6.039)
> R = 1000 95% percentile CI = (2.914, 6.038)
> R = 2000 95% percentile CI = (2.914, 6.038)
> R = 10000 95% percentile CI = (3.054, 6.056)
R=200 の結果が (2.917, 6.039)、50倍の R=10000 でも (3.054, 6.056) です。下限は 2.917 から 3.054 へ、上限は 6.039 から 6.056 へ動いただけで、臨床的な解釈が変わる違いはどこにもありません。区間の端点を得るという意味では、必要な反復回数は思われているほど多くないのです。ただしこのあと述べるとおり、裾の分位点を非対称にずらして読むBCa法では話が別になります。
一方、Rを増やしても下限が滑らかに収束せず、2.917 → 2.914 → 2.914 → 3.054 と段階的に飛んでいる点も見逃せません。中央値は元の観測値そのもの（あるいは隣り合う2値の平均）しか取り得ないため、再標本分布が飛び飛びの値になります。したがって区間の端点も順序統計量の刻みで飛び、Rを増やしても中間の値には落ち着きません。平均や回帰係数では起きない、中央値・分位点に特有の癖です。
目安としては、標準誤差だけなら R=200〜1,000 で十分、信頼区間を作るなら R≧2,000、BCa法を使うなら R≧2,000（できれば10,000）と考えておけば大きく外しません。BCaは分布の両端に近い分位点を使うため、端の推定を安定させる分だけ多めの反復が要ります。
set.seed() を書かずに実行すると、同じデータ・同じコードでも実行のたびに区間の端点が変わります。「昨日の結果と今日の結果が微妙に違う」という事故は、ブートストラップで最も起きやすいトラブルです。SAP（統計解析計画書）や解析仕様書の段階で乱数シードと反復回数 R を数値として事前に規定し、誰がいつ実行しても同じ数字が出る状態にしておいてください。提出資料では再現性の担保そのものが要求事項です。ステップ5：2群の中央値の差に層別ブートストラップを使う
実際の臨床試験で知りたいのは片方の群の中央値ではなく、2群の差、すなわち治療効果の大きさです。中央値の差の信頼区間は解析的に求めるのが難しい代表例ですが、ブートストラップならそのまま計算できます。
ただし必ず押さえるべき引数が1つあります。strata = factor(dat$group) です。これを付けると、群ごとに例数を固定したまま群内で復元抽出が行われます。付け忘れると100例全体から無差別に復元抽出することになり、再標本ごとに「対照群53例・実薬群47例」のように群サイズが変動します。極端な場合には片方の群がほとんど選ばれず、中央値の差が異常な値になったり計算が破綻したりします。そもそも「各群50例」は偶然決まった数ではなく事前に固定された設計ですから、その構造を再標本でも保つのが自然です。
diff_med <- function(data, indices) {
d <- data[indices, ]
median(d$crp[d$group == "Treatment"]) - median(d$crp[d$group == "Control"])
}
set.seed(4869)
boot_diff <- boot(data = dat, statistic = diff_med, R = 2000,
strata = factor(dat$group))
boot_diff
boot.ci(boot_diff, type = c("perc", "bca"))
> STRATIFIED BOOTSTRAP
>
>
> Call:
> boot(data = dat, statistic = diff_med, R = 2000, strata = factor(dat$group))
>
>
> Bootstrap Statistics :
> original bias std. error
> t1* -2.299143 0.3293268 0.9571412
>
> BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
> Based on 2000 bootstrap replicates
>
> CALL :
> boot.ci(boot.out = boot_diff, type = c("perc", "bca"))
>
> Intervals :
> Level Percentile BCa
> 95% (-3.586, 0.021 ) (-3.772, -0.408 )
> Calculations and Intervals on Original Scale
出力の先頭が
STRATIFIED BOOTSTRAP に変わっており、層別が効いていることが確認できます。観測された中央値の差は -2.299143 で、実薬群のほうがCRPの中央値が低いという結果です。問題はここからです。パーセンタイル法の95%信頼区間は (-3.586, 0.021) で、上限がわずかに0を超えています。つまり「差があるとは言えない」。ところがBCa法は (-3.772, -0.408) で0を含みません。つまり「差がある」。まったく同じデータ、まったく同じ2000回の再標本から出発しているのに、区間の作り方を変えただけで結論が正反対になりました。上限 0.021 は境目をほんのわずかに跨いだにすぎず、実務者としては最も居心地の悪い場面です。
なぜ割れたのか。鍵は再標本分布が観測値のどちら側へ偏っているかです。
bias が 0.3293268 と正であることが示すとおり、再標本での中央値の差は観測値 -2.299143 より系統的に大きい側（0に近い側）へずれています。パーセンタイル法はこの偏りを補正せず2.5%点と97.5%点をそのまま使うため、区間が0側へ引きずられます。BCa法は偏り（バイアス補正量 \( z_0 \)）と分布の歪み（加速度定数 \( a \)）の両方で読み取る分位点をずらすため、区間が本来の位置へ押し戻され、上限が -0.408 まで下がりました。標準誤差 0.9571412 に対してバイアス 0.3293268 は無視できる大きさではなく、補正の有無が結論を左右したわけです。なお、この差が臨床的にどの程度の意味を持つかは統計的有意性とは別に評価すべき論点です。差の大きさを解釈する視点は 効果量（Effect Size）を理解すると統計が一気に実務的になる が参考になります。
信頼区間の種類は必ず解析前に決めておいてください。いま見たとおり、パーセンタイル法とBCa法で結論が割れることは現実に起こります。結果を見てから「BCaのほうが理論的に良いから」と選び直すのは、実質的に有意になる方法を探しているのと同じで、多重性の観点からも規制上の観点からも許されません。SAPには「中央値の差の95%信頼区間はブートストラップBCa法により算出する。反復回数 R=2000、乱数シード=…、層別変数=投与群」というレベルまで書き切っておく。そうすればどちらの結果が出ても堂々と報告できます。
ステップ6：並べ替え検定・ウィルコクソン検定と突き合わせる
結論が割れたとき、どう判断すればよいのか。別の道具と突き合わせます。その前に整理しておきたいのは、ブートストラップと並べ替え検定（permutation test）は目的が違うということです。ブートストラップは観測データを母集団の代理とみなして再標本を作り、「差はどれくらいか」という推定と区間推定に答えます。並べ替え検定は帰無仮説が正しいと仮定して群ラベルをシャッフルし、「差がないと言えるか」という検定に答えます。前者は推定の道具、後者は検定の道具です。
obs <- median(dat$crp[dat$group == "Treatment"]) -
median(dat$crp[dat$group == "Control"])
set.seed(4869)
perm <- replicate(2000, {
g <- sample(dat$group)
median(dat$crp[g == "Treatment"]) - median(dat$crp[g == "Control"])
})
p_perm <- (sum(abs(perm) >= abs(obs)) + 1) / (2000 + 1)
round(c(observed_diff = obs, p_permutation = p_perm), 4)
wilcox.test(crp ~ group, data = dat)
> observed_diff p_permutation
> -2.2991 0.0005
>
> Wilcoxon rank sum test with continuity correction
>
> data: crp by group
> W = 1716, p-value = 0.001332
> alternative hypothesis: true location shift is not equal to 0
p値の計算で分子と分母にそれぞれ 1 を足しているのは、観測値以上に極端な並べ替えが1回も現れなかったときに p=0 と報告してしまうのを避けるための標準的な補正です。観測データそのものを並べ替えの1つとして数える、という解釈に対応します。この補正により、2000回の並べ替えで得られるp値には \( 1/2001 \) という下限が入ります。
並べ替え検定のp値は 0.0005 で、これは上記の下限 \( 1/2001 \) に張り付いた値です。つまり2000回の並べ替えの中に、観測された差 -2.2991 と同じかそれ以上に極端なものが1回も現れなかったということです。ウィルコクソン検定も W = 1716、p-value = 0.001332 で明確に有意です。2つの検定はいずれも「差がある」を支持しており、ステップ5でBCa法が出した結論と整合します。パーセンタイル法だけが0をまたいだ格好で、BCa法を既定として採用すべき実務的な根拠がここで得られます。
ただし3つの手法が答えている問いは微妙に違うので、p値や区間が完全に一致する必要はありません。並べ替え検定の帰無仮説は「2群の分布が同一である」という強い主張で、位置だけでなくばらつきや形の違いにも反応します。ウィルコクソン検定も厳密には同じ帰無仮説を置いていますが、2群の分布が形とばらつきを共有するという位置シフトモデルを仮定すれば、中央値のずれを順位で検出する検定として解釈できます。そしてBCa法が答えているのは「中央値の差はどの範囲にありそうか」という推定の問いです。数字が食い違うのは当然で、大事なのは複数の道具が同じ方向を指しているかどうかです。順位に基づく検定の中身は ノンパラメトリック検定の代表格：ウィルコクソン検定の数理と実装 で詳しく解説しています。
実務的には、主要な結果はSAPで事前に規定した方法（ここではBCa法）で報告し、並べ替え検定やウィルコクソン検定は感度分析として並べる整理が素直です。方法を後から選び直すのではなく、複数の方法が同じ結論に至ることを示す。それが結論の頑健性を主張する正しい作法です。
ブートストラップは万能ではない
「分布形を仮定しなくても信頼区間が作れる」を「どんな状況でも正しい区間が出る」と読み替えてしまうと、実務では手痛い誤りにつながります。分布の仮定から自由であることと、有限の標本で正しく振る舞うことは別の話だからです。ここでは真値がわかる人工データで、作った区間が実際に真値を捉えるかを数えて確かめます。
被覆確率を数えてみる
信頼区間の性能を測るもっとも素直な指標が被覆確率（coverage probability）です。95%と名乗る以上、実験を繰り返したときに区間が真値を含む割合は0.95に近くなければなりません。対数正規分布（\( \mu=1, \sigma=1 \)）から n=25 の標本を1,000回発生させ、真の平均 \( \exp(1+1^2/2)=4.4817 \) を t区間・パーセンタイル法・BCa法の95%信頼区間が何回捉えたかを数えます。あわせて平均の区間幅も記録します。なお1,000回のシミュレーションそれぞれでブートストラップを回すため、計算時間の都合から再標本回数は R=999 としています。被覆確率の水準がこの反復回数で大きく変わることはありません。
set.seed(4869)
nsim <- 1000
n <- 25
true_mean <- exp(1 + 1^2 / 2)
mean_fun <- function(d, i) mean(d[i])
hit <- matrix(FALSE, nsim, 3,
dimnames = list(NULL, c("t", "percentile", "BCa")))
width <- matrix(NA_real_, nsim, 3,
dimnames = list(NULL, c("t", "percentile", "BCa")))
for (s in 1:nsim) {
x <- rlnorm(n, meanlog = 1, sdlog = 1)
ci_t <- t.test(x)$conf.int
bo <- boot(x, mean_fun, R = 999)
ci_b <- boot.ci(bo, type = c("perc", "bca"))
hit[s, "t"] <- ci_t[1] <= true_mean && true_mean <= ci_t[2]
hit[s, "percentile"] <- ci_b$percent[4] <= true_mean && true_mean <= ci_b$percent[5]
hit[s, "BCa"] <- ci_b$bca[4] <= true_mean && true_mean <= ci_b$bca[5]
width[s, "t"] <- ci_t[2] - ci_t[1]
width[s, "percentile"] <- ci_b$percent[5] - ci_b$percent[4]
width[s, "BCa"] <- ci_b$bca[5] - ci_b$bca[4]
}
round(colMeans(hit), 3)
round(colMeans(width), 3)
出力は次のとおりです。1行目が被覆確率、2行目が平均の区間幅です。
> t percentile BCa
> 0.890 0.886 0.897
> t percentile BCa
> 4.218 3.865 4.483
被覆確率は t区間 0.890、パーセンタイル法 0.886、BCa 0.897 で、3手法とも名目の95%を下回りました。10回に1回以上は真値を外している計算です。平均区間幅は t 4.218、パーセンタイル 3.865、BCa 4.483 で、被覆の順序と幅の順序がぴたりと対応しています。パーセンタイル法は区間がもっとも狭いぶん被覆がもっとも低く（3.865／0.886）、BCaはもっとも広いぶん被覆がもっとも高い（4.483／0.897）という取引関係です。狭くて当たる区間という都合のよいものは存在せず、「分布の仮定が要らない」は「小標本でも正しい」を意味しません。
標本サイズを増やすとどうなるか
では例数を増やせば解決するのでしょうか。まったく同じ設定で n だけを100に変えて測り直します。
set.seed(4869)
nsim <- 1000
n <- 100
hit2 <- matrix(FALSE, nsim, 3,
dimnames = list(NULL, c("t", "percentile", "BCa")))
for (s in 1:nsim) {
x <- rlnorm(n, meanlog = 1, sdlog = 1)
ci_t <- t.test(x)$conf.int
bo <- boot(x, mean_fun, R = 999)
ci_b <- boot.ci(bo, type = c("perc", "bca"))
hit2[s, "t"] <- ci_t[1] <= true_mean && true_mean <= ci_t[2]
hit2[s, "percentile"] <- ci_b$percent[4] <= true_mean && true_mean <= ci_b$percent[5]
hit2[s, "BCa"] <- ci_b$bca[4] <= true_mean && true_mean <= ci_b$bca[5]
}
round(colMeans(hit2), 3)
> t percentile BCa
> 0.920 0.929 0.922
n=100 では t区間 0.920、パーセンタイル法 0.929、BCa 0.922 と、n=25 のとき（0.890／0.886／0.897）から確かに改善しました。しかしそれでも3手法とも95%には届いていません。ブートストラップは標本サイズが大きくなるにつれて正しくなる漸近的な手法であり、対数正規分布のように歪みの強い分布では n=100 でも名目水準に届かないことがあるわけです。近似の精度が例数だけでなく分布の歪みに左右される点は、中心極限定理とはで扱った収束の速さの問題です。
この結果を「ブートストラップは使えない」と読まないでください。同じ条件で t区間も0.890であり、どの手法も等しく苦しいのが実態です。歪んだ分布から少数例を取ったときに名目どおりの区間を出せる万能の手法は存在しません。重要なのは「どこかに万能の手法があるはずだ」という期待を捨て、想定する分布と症例数で事前にシミュレーションを回し、使う手法の被覆と区間幅を確認しておくことです。解析計画の段階で済ませられる作業であり、結果を見てから慌てて選び直すよりはるかに健全です。
原理的に破綻する統計量がある
被覆が名目より少し低いのはまだ程度問題です。しかしブートストラップには、例数をいくら増やしても原理的に機能しない統計量があります。代表例が標本最大値です。一様分布から n=50 の標本を取り、その最大値をブートストラップしてみます。
set.seed(4869)
x <- runif(50, min = 0, max = 10)
boot_max <- replicate(2000, max(sample(x, replace = TRUE)))
round(c(observed_max = max(x),
n_unique = length(unique(boot_max)),
prob_equal_max = mean(boot_max == max(x))), 3)
quantile(boot_max, c(0.025, 0.5, 0.975))
> observed_max n_unique prob_equal_max
> 9.943 7.000 0.625
> 2.5% 50% 97.5%
> 9.612 9.943 9.943
再標本は観測データからの復元抽出ですから、再標本の最大値は観測された最大値 9.943 を原理的に超えられません。2000回の再標本で得られた最大値は7通りの値しか取らず、そのうち62.5%は観測最大値ちょうどに一致しています。理論的にも、n個から復元抽出したときに観測最大値が選ばれる確率は \( 1-(1-1/n)^n \approx 1-e^{-1} \approx 0.632 \) となり、n=50 での実測 0.625 とよく一致します。結果として97.5%点も中央値もどちらも 9.943 に潰れ、区間としての体をなしていません。
一般化すると、ブートストラップがうまく働くのは、対象の統計量が経験分布の滑らかな汎関数である場合です。平均や中央値、中央付近の分位点はこれを満たします。一方、最大値・最小値のような端点に依存する統計量、真のパラメータが定義域の境界にある場合（分散成分がゼロ近傍など）、極値統計量では一致性が成り立ちません。m-out-of-n ブートストラップやサブサンプリングという修正法はありますが、標準的な boot() をそのまま回して済む話ではありません。
| 統計量・状況 | ブートストラップの可否 | 理由・注意 |
|---|---|---|
| 平均・中央値・中央付近の分位点 | ◎ 使える | 経験分布の滑らかな汎関数にあたる。ただし小標本では被覆が名目を下回る（本記事のn=25で0.886〜0.897） |
| 相関係数・回帰係数 | ◎ 使える | 例数が十分なら問題ない。回帰では個体ごと抽出（ケース法）と残差抽出のどちらを使うかを明示する |
| RMST・Win Ratio などの派生指標 | ◎ 主戦場 | 分散の閉じた式が無い、あるいは近似が怪しい推定量。ブートストラップの価値がもっとも出る場面 |
| 標本最大値・最小値 | × 使えない | 分布の端に依存し一致性を持たない。m-out-of-n ブートストラップやサブサンプリング、極値理論を検討する |
| 分散成分がゼロ近傍（境界パラメータ） | × 原則不可 | 真値が定義域の境界にあると再標本分布が歪み、通常のブートストラップは破綻する |
| 反復測定・クラスター化データ | △ 抽出単位を変えれば可 | 個票の単純復元抽出は独立性の仮定を壊す。被験者・施設単位で抽出するクラスターブートストラップを使う |
最後の行は実務でもっとも起こりやすい落とし穴です。単純な復元抽出は、観測値が互いに独立であることを前提にしています。反復測定デザインや施設ごとに患者が束になったクラスター化データで個票をそのまま抽出すると、相関構造が壊れて標準誤差を過小評価し、区間が不当に狭くなります。対処は抽出の単位をデータの独立単位に合わせることで、被験者内相関があるなら被験者を丸ごと1つの塊として抽出するクラスターブートストラップを用います。相関構造をモデル側に組み込む選択肢もあり、その代表格はGEE（一般化推定方程式）とはで解説しています。
実務でのポイント
最後に、「では自分の試験のどこで使うのか」に答える形で整理します。
どんな場面でブートストラップが効くか
判断基準はシンプルで、分散の閉じた式が無いか、あっても近似が怪しい推定量かどうかです。平均の差のように公式が確立していて正規近似も効く場面なら、わざわざ持ち出す必要はありません。
第一に、中央値や分位点です。歪んだ検査値では臨床的に意味を持ちますが、標準誤差に扱いやすい閉じた式がありません。第二に、RMST（制限付き平均生存時間）解析ガイドで扱ったRMSTの差や比、Win Ratio法による複合エンドポイント解析のWin Ratioのように、複数の情報を組み合わせて作る指標です。デルタ法による近似の妥当性が症例数やイベント数に強く依存します。第三に、ROC曲線から得られるAUCのような曲線由来の要約量。第四に、逆確率重み付け（IPTW）によるATE推定のように、傾向スコアを推定してから重み付き推定量を作る二段構えの手続きです。第一段階の推定誤差まで含めた分散を書き下すのは容易ではありませんが、傾向スコアの推定ごと再標本に含めてしまえばその不確実性が自然に区間へ反映されます。比や差の対数変換のような非線形変換も、素直に扱える対象です。
層別・クラスターを再標本にも反映させる
再標本の取り方はデータの構造を写していなければなりません。意識すべきはランダム化の構造と測定の構造です。層別因子でブロックを分けて割り付けたのであれば、再標本でも各層の人数が保たれるよう strata= を指定します。層別割付の考え方はランダム化割付の方法とR実装にまとめています。試験デザインで層別したなら解析でも層別するという原則が、再標本の取り方として現れるわけです。
SAPへの書き方
自由度が高いぶん、事前規定の重要性は通常の手法より高くなります。統計解析計画書（SAP）には、信頼区間の種類（たとえばBCa）、反復回数（R=2000など）、乱数シードの値、層別・クラスターの抽出単位、使用ソフトとパッケージのバージョン（R 4.5.1、boot 1.3.31 のように）まで書き切ってください。
ステップ5で見たとおり結論が割れることは現実に起こりますから、結果を見てから区間の種類を選び直さないという規律は徹底してください。あわせて、監査や当局照会に備えてシードとスクリプトだけで同じ数字が再現できる状態を保つことも必須です。
規制当局から見たときの位置づけ
主要評価項目の主解析に据えるより、①分散の式が存在しない推定量に区間推定を与える用途、②正規近似の妥当性が疑わしいときに主解析の頑健性を確かめる感度分析、この二つで使うほうが受け入れられやすい傾向があります。いずれも事前に規定されていることが何より重要で、事後に持ち出した区間は解釈の重みが下がります。
並べ替え検定との使い分け
ステップ6で整理したとおり、帰無仮説のもとでの確率を評価したいなら並べ替え検定、推定量のばらつきを評価して区間を作りたいならブートストラップが基本の切り分けです。臨床試験の報告ではp値だけでなく効果の大きさと信頼区間が求められますから、実務的には両方を併記するのが自然です。
計算コストの現実
中央値や平均なら R=2000 は一瞬で終わります。しかし各反復の中でモデルを1本推定する入れ子構造、たとえば傾向スコアの推定から重み付け、効果推定までを丸ごと再標本する場合は計算時間が跳ね上がり、やり直しやレビュー対応が現実的でなくなります。並列化（parallel= の指定）を前提にした設計と反復回数の見積もりを、計画段階で済ませておくことをおすすめします。
| 状況 | 推奨されるアプローチ | 判断のポイント |
|---|---|---|
| 正規近似が妥当な平均の比較（例数が十分） | t検定・ANCOVAなど標準的手法 | 閉じた式があるならそれを使う。ブートストラップを足しても情報は増えない |
| 歪んだ分布の中央値・分位点 | ブートストラップ（BCa）＋ノンパラ検定 | 区間と検定の結論が整合するかを必ず突き合わせる。小標本では被覆低下を前提に解釈する |
| RMST・Win Ratio・AUCなどの派生指標 | ブートストラップ（層別・BCa） | 分散の閉じた式が無い典型例。反復回数・シード・抽出単位をSAPに明記する |
| 反復測定・クラスターデータ | クラスターブートストラップ／GEE・混合モデル | 抽出単位は被験者・施設。個票をそのまま抽出すると標準誤差を過小評価する |
| 標本最大値など端点統計量 | 通常のブートストラップは使わない | 一致性が成り立たない。極値理論、m-out-of-n、サブサンプリングを検討する |
・使いどころは「分散の閉じた式が無い／近似が怪しい」推定量。中央値・分位点、RMST、Win Ratio、AUC、傾向スコアを経由した推定量が典型
・再標本の取り方はデザインを写す。層別割付なら
strata=、反復測定・クラスターなら被験者・施設単位で抽出する・SAPには区間の種類・反復回数・乱数シード・抽出単位・ソフトとパッケージのバージョンまで事前に規定する
・結果を見てから区間の種類を選び直さない。シードとスクリプトだけで同じ数字が再現できる状態を保つ
・主解析に据えるより、区間推定の手段または感度分析としての位置づけが受け入れられやすい
参考書籍
ブートストラップは、数式を眺めるよりも「実際に手を動かして分布を作ってみる」ほうが理解の進みが早い手法です。ここでは、リサンプリングそのものを正面から解説した本、ブートストラップを計算統計という広い文脈に置き直してくれる本、そして臨床試験の評価指標の側からブートストラップの必要性を実感できる本を3冊選びました。いずれも本記事のRコードと行き来しながら読むと、BCa法や被覆確率の議論がぐっと具体的になります。
まとめ
本記事では、ブートストラップ法の原理から4種類の信頼区間の作り方、boot パッケージによるR実装、そして手法としての限界までを紹介しました。中心にあるのは「母集団からの標本抽出を、標本からの復元抽出で置き換える」という一つの発想です。
記事の中で示した数値を振り返っておきます。2群の中央値の差では、同じデータ・同じ2,000回の再標本から、パーセンタイル法が (-3.586, 0.021) と0をまたぐ一方、BCa法は (-3.772, -0.408) と0を含まず、区間の作り方だけで結論が割れました。被覆確率のシミュレーションでは、n=25 のとき t区間・パーセンタイル法・BCa がそれぞれ 0.890 / 0.886 / 0.897 と名目95%を下回り、n=100 に増やしても 0.920 / 0.929 / 0.922 にとどまりました。さらに標本最大値では、2,000回の再標本のうち62.5%が観測最大値 9.943 ちょうどに張り付き、区間が意味を失いました。
ここから引き出せる位置づけは明確です。ブートストラップは分布の仮定から自由になる道具ではなく、正規近似という一つの近似を、再標本という別の近似で置き換える道具です。置き換えた先にも、標本サイズや分布の歪み、統計量の滑らかさ、観測値の独立性という前提が残っています。だからこそ事前規定と性能確認が要るのであり、その準備さえできていれば、分散の式が書けない推定量にも堂々と信頼区間を付けられる強力な武器になります。
関連するトピックとしては、生存時間解析の派生指標をRで動かす手順をまとめたRMSTをRで実装する方法、複合エンドポイントの新しい要約方法を扱ったWin Ratio法による複合エンドポイント解析、信頼区間そのものの考え方を数式と図で整理した区間推定入門：数式と図解で理解する信頼区間の世界もあわせてご覧ください。
まずは手元のデータで boot() を一度回し、パーセンタイル法とBCaの区間を並べて眺めてみることをおすすめします。二つの区間がどれくらいずれるかをご自身の目で確かめるところから、ブートストラップを実務で使う感覚をつかんでいただければと思います。
