同じ対象者の開始時と1年後の検査値を比べると、平均変化は計算できても、その値が対象者の入れ替わりでどれくらい動くかは、平均だけでは分かりません。また、中央値、回帰モデルから得る予測値、複雑な効果指標になると、標準誤差を求める式がすぐには見つからない場合があります。
そこで、観測した対象者から繰り返し標本を作り直し、毎回同じ解析を行う方法がブートストラップです。ただし、行をランダムに選ぶだけでは、研究で知りたい不確実性を再現できないこともあります。同じ人の測定値を別々に選ぶのか、人ごとに選ぶのかで、作り直している研究が変わってしまうためです。
そして、反復回数を10,000回へ増やすことと、対象者を10,000人へ増やすことは異なります。計算による揺れと、元のデータに残る情報不足を区別すると、信頼区間の読み方や、ブートストラップが役立たない場面まで見通しやすくなります。
推定値と標本分布の違い
20人の平均変化が0.805だったとしても、同じ集団から別の20人を集めれば、通常は別の平均になります。このように、同じ手順で研究を繰り返した場合の推定値の分布が、標本分布です。
対象者ごとの測定値の分布と、20人から求める平均の標本分布は、別のばらつきを表しています。前者の広がりは人同士の違い、後者の広がりは平均という推定値の揺れに関わります。また、標準偏差と標準誤差を区別するのも、この違いを表現するためです。
平均なら標本標準偏差を人数の平方根で割る方法が知られていますが、二つの回帰モデルを組み合わせた効果指標や、複雑な予測手順では、同じ程度に簡単な式が使えるとは限りません。そのような場合でも、対象者が入れ替わったデータへ解析を繰り返せるなら、推定値の揺れを計算で追う余地があります。
ただし、ブートストラップが成立するかどうかは、解析式の複雑さとは別の条件で決まります。どの集団から、どの単位が独立に得られたと考えるか、そして推定したい量が再抽出によって安定して近似できるかが前提になります。
観測分布から標本を作る考え方
Efronが1979年に提案したブートストラップでは、未知の母集団分布を、観測データが作る経験分布へ置き換えて考えます。観測された値へ同じ確率を割り当て、その分布から標本を作り直すことが、通常のノンパラメトリック・ブートストラップの出発点です。
独立で同じ分布に従う観測値を$X_1,\ldots,X_n$とすると、一変量の経験分布関数は次のように書けます。
$$\widehat F_n(x)=\frac{1}{n}\sum_{i=1}^{n}I(X_i\le x) \tag{1}$$
式(1)は、観測された$n$個の値のうち、$x$以下に入るものの割合を求めています。$I$は条件が成り立つ場合に1、成り立たない場合に0となる指示関数で、各観測者へ$1/n$ずつの重みを付けています。
実際の解析では、経験分布関数を画面に描いてから抽出する必要はありません。対象者番号を1から$n$までの間で復元抽出すると、経験分布からの抽出と同じ処理になります。また、一人分の観測が複数の変数を持つ場合には、その変数の組へ重みを割り当てます。
ここで、母集団で知りたい量を$\theta$、元標本からの推定値を$\widehat\theta$、再抽出標本からの推定値を$\widehat\theta^*$とします。通常の平方根$n$で尺度をそろえられる推定量では、近似したい関係は次の形です。
$$\mathcal L_F\!\left\{\sqrt n(\widehat\theta-\theta)\right\} \approx \mathcal L_*\!\left\{\sqrt n(\widehat\theta^*-\widehat\theta) \mid X_1,\ldots,X_n\right\} \tag{2}$$
式(2)は、母集団からの抽出で生じる推定誤差を、元の観測データを固定した再抽出で生じる推定値のズレによって近似することを表しています。$\mathcal L_F$は母集団からの標本分布、$\mathcal L_*$は観測データを条件としたブートストラップ分布です。
左側は未知の$\theta$を中心に考えますが、右側では計算できる$\widehat\theta$を中心にします。そのため、ブートストラップの推定値が元の推定値の周辺に集まることは、真値を発見したことを意味しません。
BickelとFreedmanは、多くの設定での漸近的な成立条件とともに、近似が成立しない例も論じています。特に、最大値のように別の尺度で誤差を考える量へ、式(2)をそのまま当てはめることはできません。
復元抽出と対象者の重複
20人から20人分を復元抽出すると、同じ人が2回以上現れたり、一度も現れない人がいたりします。一つの標本内で対象者番号が重複しても、入力を間違えたわけではなく、観測者ごとの重みを変えながら推定値の揺れを作っています。
反対に、20人を非復元で20人分選ぶと、全員が1回ずつ現れ、順番しか変わりません。順番に依存しない平均や通常の回帰係数は同じ値になるので、その操作では標本変動を再現できません。
$n$人から$n$回復元抽出する場合、ある一人が一度も選ばれない確率と、少なくとも一度選ばれる人数の期待値は、次のようになります。
$$\Pr_*(i\text{ が未抽出})=\left(1-\frac1n\right)^n, \qquad \mathbb E_*(N_{\mathrm{unique}}) =n\left\{1-\left(1-\frac1n\right)^n\right\} \tag{3}$$
式(3)は、各回でその人以外を選ぶ確率を$n$回分掛け、そこから一度以上選ばれる割合を求めています。$N_{\mathrm{unique}}$は再抽出標本に含まれる異なる対象者の人数です。
$n=20$では異なる対象者の期待人数は約12.83人で、大きな$n$では元の人数の約63.2%へ近づきます。約63.2%という値は「各標本が元データの63.2%の情報しか持たない」と一律に評価するための割合ではなく、異なる対象者が含まれる割合の計算です。
また、一つの再抽出標本の行数は20のままで、12.83という数字が指しているのは異なる対象者の期待人数です。重複を削除すると、意図した復元抽出とは異なる標本になるため、解析へ渡す前に対象者番号の重複を取り除く操作は推奨しません。
同じ人の測定値を一緒に選ぶ理由
開始時と1年後の測定値を比較する研究では、一人分の情報は二つの測定値の組です。開始時だけを20回選び、1年後を別に20回選ぶと、元の人同士の対応が崩れ、対応のある変化を調べる研究とは異なる標本を作ってしまいます。
同じ人の開始時を$X_i$、1年後を$Y_i$、変化を$D_i=Y_i-X_i$とすると、変化の分散には個人内の共分散が入ります。
$$\operatorname{Var}(D_i) = \operatorname{Var}(Y_i)+\operatorname{Var}(X_i) -2\operatorname{Cov}(X_i,Y_i) \tag{4}$$
式(4)は、変化のばらつきが各時点のばらつきだけでなく、同じ人の値同士がどれくらい一緒に動くかにも依存することを表しています。$\operatorname{Cov}$は開始時と1年後の共分散です。
後で使うデータでは、両時点の相関係数は約0.723でした。平均変化の通常の標準誤差は約0.179ですが、二つの時点を独立として分散を足すと約0.316になります。この例では、対応を無視すると標準誤差が大きくなる方向です。
依存を無視した時に標準誤差がどちらへ動くかは、共分散の符号と、それが分散にどう入るかで決まります。個人内の変化と、多施設から集めた人全体の平均では共分散の入る形が異なるため、何を計算する時の依存かを確認する必要があります。
今回は一人分の行を選び、その行にある開始時と1年後を一緒に使います。同じ人が2回選ばれれば、その人の測定値の組も2回使うことになり、元の対応関係を保った復元抽出になります。
20人のデータによる実行例
bootパッケージのcd4データセットを使用して、平均変化のブートストラップを確かめます。HIV陽性の20人について、試験開始時と薬剤使用1年後のCD4数を記録したデータで、値は100個を単位として表されています。
DiCiccioとEfronの1996年の論文でも、このデータが信頼区間の説明に使われています。ただし、ここで計算するのは平均変化で、原著の相関係数の区間をそのまま再現する例ではありません。
また、このデータ表には無治療の対照群がないため、開始時からの平均変化を、そのまま薬剤の因果効果と解釈することはできません。ブートストラップは比較群を追加する方法ではなく、ここでは観測された人の変化を要約する推定量の不確実性を扱います。
以下のコードはR4.4.2、bootパッケージ1.3-31で実行しました。後に掲載する描画コードも含め、上から順に同じRセッションで実行すると、図は作業フォルダへPNG形式で保存されます。関数に渡す処理は、選ばれた人の平均変化と、その平均変化の分散推定値を返すものです。復元抽出そのものはbootパッケージが行います。
set.seed(42)
library(boot)
study <- cd4
change_statistic <- function(data, indices) {
change <- data$oneyear[indices] - data$baseline[indices]
c(mean(change), var(change) / length(change))
}
resampled <- boot(study, statistic = change_statistic, R = 10000)
resampled
intervals <- boot.ci(resampled, type = "all", index = c(1, 2))
intervals
t.test(study$oneyear, y = study$baseline, paired = TRUE)$conf.int#> 元標本の平均変化 0.805000
#> ブートストラップbias -0.000692
#> ブートストラップSE 0.174896
#>
#> 95%信頼区間
#> Normal 0.4629 ～ 1.1485
#> Basic 0.4590 ～ 1.1350
#> Studentized 0.4569 ～ 1.2116
#> Percentile 0.4750 ～ 1.1510
#> BCa 0.4925 ～ 1.1715
#>
#> 対応のあるt区間
#> 0.4298 ～ 1.1802開始時の平均は3.288、1年後は4.093で、平均変化は0.805でした。100個を1単位とした尺度なので、0.805は80.5個分の変化に相当しますが、その医学的な意味や治療選択について結論付ける計算ではありません。
出力のbiasは、ブートストラップ平均から元の平均を引いた値で、臨床研究の交絡や選択バイアスを評価した数字ではありません。また、Studentizedの計算に必要な分散を2番目の戻り値として渡しているため、元の全出力にはその分散についての行も現れます。上には平均変化に関する部分を抜粋しています。
五つの区間はどれも0をまたぎませんが、そこから区間の作り方まで同じだと読むのは飛躍です。そして、最も短い区間を選ぶことが、最も信頼できる方法を選ぶことにもなりません。各区間の構成は次の記事で、数式と実行値を対応させます。
観測された対応関係
散布図では、一つの点が一人の開始時と1年後の組です。破線は両時点が同じ値になる位置なので、その上にある点は1年後の値が高い対象者を表します。
png("bootstrap-paired-data.png", width = 1100, height = 688, res = 200, type = "cairo")
plot(study$baseline, y = study$oneyear, xlim = c(1, 7), ylim = c(1, 7),
xlab = "Baseline CD4 (hundreds)", ylab = "One-year CD4 (hundreds)", pch = 19)
abline(a = 0, b = 1, lty = 2, col = "grey40")
dev.off()
図1の右上の対象者が再抽出標本で複数回選ばれる場合もあれば、一度も選ばれない場合もあります。その重みの変化によって平均変化が動きますが、動かしているのは観測済みの20人へ配る重みだけです。
標準誤差とバイアスの計算
復元抽出を$B$回繰り返して得た推定値を$\widehat\theta_1^*,\ldots,\widehat\theta_B^*$とします。ここでの$B$は計算回数で、元の対象者数$n$とは別の記号です。
$$\overline\theta^*=\frac1B\sum_{b=1}^{B}\widehat\theta_b^*, \qquad \widehat{\operatorname{SE}}_{\mathrm{boot}} = \sqrt{\frac1{B-1} \sum_{b=1}^{B}(\widehat\theta_b^*-\overline\theta^*)^2} \tag{5}$$
式(5)は、再抽出のたびに得た推定値の平均と標準偏差を求めています。この推定値の標準偏差を、元の推定量の標準誤差の近似として使います。
ブートストラップ標準誤差を、もう一度$\sqrt B$で割ってはいけません。そうすると、推定量の標本変動ではなく、$B$回の平均を計算したことによるモンテカルロ誤差へ近い量を計算することになります。
今回の平均変化の標準誤差は約0.174896です。元の20人の変化そのものの標準偏差は約0.801659なので、両者が表すばらつきの対象は異なります。
元の推定値を残す判断
ブートストラップによるバイアス推定は、次の差として計算します。
$$\widehat{\operatorname{bias}}_{\mathrm{boot}} = \overline\theta^*-\widehat\theta \tag{6}$$
式(6)は、再抽出による推定値の平均が、元標本の推定値からどちらへどれくらいズレたかを求めています。母集団の真値を直接使っていないため、ここで分かるのは再抽出による推定値の偏りだけで、実際の研究で生じる交絡や選択の偏りは別の話になります。
今回は約−0.000692でした。そして、標本平均の通常の復元抽出では、観測データを固定した時のブートストラップ平均の期待値は元の標本平均と一致するため、この程度の差には有限回の乱数計算による揺れが含まれます。
バイアスを引いた$2\widehat\theta-\overline\theta^*$を点推定値とする方法もありますが、補正量自体も推定しているので、常に平均二乗誤差が小さくなるとは限りません。今回のように標本平均を扱う例で、わずかな計算上の差を理由に点推定値を置き換える必要はなく、元の平均0.805を報告し、標準誤差や区間を追加します。
式による標準誤差との照合
標本平均については、観測データを固定した通常の$n$個復元抽出での分散を、式から直接求められます。元標本の不偏分散を$s^2$とすると、その分散は$(n-1)s^2/n^2$です。
$n=20$の今回のデータでは、式からのブートストラップ標準誤差は約0.174718で、10,000回から得た約0.174896と近い値でした。また、通常の平均の標準誤差$s/\sqrt n$は約0.179257で、有限標本では両者が完全には一致しません。
ここでのズレは、経験分布の分散が分母$n$で定義されることと、不偏分散が分母$n-1$を使うことに関わります。ブートストラップ標準誤差が少し小さかった理由は、この分母の取り方にあります。対象者が増えたわけでも、欠測が埋まったわけでもありません。
再抽出による推定値の分布
平均変化を10,000回計算すると、一つの平均値だけでは見えなかった広がりを描けます。図2では、その10,000個の平均変化に密度曲線を当てています。
png("bootstrap-mean-distribution.png", width = 1100, height = 688, res = 200, type = "cairo")
plot(density(resampled$t[, 1]), main = "", xlab = "Bootstrap mean change", lwd = 2)
abline(v = resampled$t0[1], lty = 2, col = "red3")
dev.off()
このデータでは、平均変化の分布は単峰に見えます。ただし、密度曲線には平滑化が入っているため、細かな離散性や少数の極端な再抽出結果は、曲線だけでは見落とす場合があります。
また、元データの分布が非正規でも、平均の標本分布は比較的対称に見えることがあります。反対に、比や極端な分位点などの推定量では、観測値が整った分布でも、推定値の分布が強く歪む場合があります。確認したいのは、解析する推定量の分布です。
抽出単位を決めてから回数を決める
ここまでで動かした設定は、何を一単位として選び直すかと、それを何回繰り返すかの二つでした。平均変化0.805に対する標準誤差0.174896は、一人分の測定値の組をそのまま選ぶという指定から出ています。
開始時と1年後を別々に選べば、回数を100,000回へ増やしても、対応のある研究の標準誤差には近づきません。回数が小さくするのは乱数計算による揺れであって、何を再現しているかというズレは回数では埋まらないためです。
欠測を含むデータで、どの段階の不確実性をどこで扱っているかという問題は、MICEによる多重代入とも共通します。また、開発したモデルの性能を同じデータで測るときの揺れについては、予測モデルの検証で扱っています。
次の記事では、ここで得た10,000個の平均変化から信頼区間を作ります。パーセンタイル、基本、正規近似、Studentized、BCaの五つで、同じデータの端点がどこまで動くかを比べます。
参考文献
Efron B. Bootstrap methods: Another look at the jackknife. The Annals of Statistics. 1979;7(1):1-26.
Bickel PJ, Freedman DA. Some asymptotic theory for the bootstrap. The Annals of Statistics. 1981;9(6):1196-1217.
DiCiccio TJ, Efron B. Bootstrap confidence intervals (with discussion). Statistical Science. 1996;11(3):189-228.
