前の記事では、20人のCD4データから平均変化0.805と、対象者単位のブートストラップ標準誤差0.174896を求めました。ここから95%信頼区間を作る段になると、パーセンタイル、基本、正規近似、Studentized、BCaという五つの名前が並び、bootパッケージはそれらを一度に返してきます。
五つはどれも同じ10,000個の再抽出推定値から作られますが、その使い方が違うため、端点は同じ値になりません。名前を選ぶ前に、それぞれが何を分位点として扱い、どの尺度で区間を組み立てているのかを押さえておくと、出力のどこを見ればよいかが決まります。
この記事では、五つの区間を数式と実行値で対応させ、どれを主に報告するかの決め方、そして計算回数を増やしたときに端点がどこまで落ち着くのかまでを扱います。前提となる記号と実行環境は、経験分布から標本を作り直す仕組みと共通です。
パーセンタイル区間と基本区間
前の記事で得た10,000個の平均変化を小さい順に並べ、下側2.5%と上側97.5%の位置を使うと、パーセンタイル区間になります。作り方が直感的なので広く使われますが、ブートストラップ分布の中央95%を取れば、どの推定量でも母数を95%で覆えるという保証はありません。
再抽出推定値の$p$分位点を$q_p^*$、有意水準を$\alpha$とすると、パーセンタイル区間は次の形です。
$$\mathrm{CI}_{\mathrm{perc}} = \left[q_{\alpha/2}^*,\ q_{1-\alpha/2}^*\right] \tag{1}$$
式(1)は、ブートストラップ推定値の分布から、両側に$\alpha/2$ずつ残す境界を選んでいます。95%区間では$\alpha=0.05$なので、2.5%点と97.5%点を使います。
今回の区間は0.4750〜1.1510でした。これは20人の平均変化についての区間で、一人のCD4変化が95%の確率でこの範囲へ入るという意味ではありません。また、パラメータの事後分布から取ったベイズ信用区間とも区別します。
元の推定値を中心に反転する区間
ブートストラップ推定値そのものではなく、元の推定値からの誤差に注目すると、基本区間、basic intervalが得られます。
$$\mathrm{CI}_{\mathrm{basic}} = \left[ 2\widehat\theta-q_{1-\alpha/2}^*,\ 2\widehat\theta-q_{\alpha/2}^* \right] \tag{2}$$
式(2)は、ブートストラップ分位点を元の推定値$\widehat\theta$の周りで反転しています。推定誤差が大きい側へ同じ幅だけ区間を伸ばすのではなく、その誤差を元の推定値から差し引く向きで区間を作ります。
平均変化0.805の今回の例では、基本区間は0.4590〜1.1350でした。パーセンタイル区間と幅がほぼ同じでも、上下の位置が異なるため、論文では区間の種類まで書き分けます。
また、パーセンタイル区間は、単調変換した推定値の分位点を戻す操作と整合しやすい性質を持ちますが、基本区間や正規近似区間では、どの尺度で区間を作るかによって結果が変わります。オッズ比なら元の尺度と対数尺度、相関係数なら元の尺度とFisher変換後の尺度を区別する必要があります。
正規近似と標準化による区間
標準誤差を得られたら、元の推定値の左右へ約1.96倍の幅を取る方法が候補になります。ただし、ブートストラップ標準誤差を使うことと、区間で正規近似を避けることは別の選択です。
正規近似区間には、中心を元の推定値にする形と、ブートストラップのバイアス推定で補正する形があります。前の記事で示したbootパッケージの出力では、後者の形で計算されています。
$$\mathrm{CI}_{\mathrm{norm}} = \widehat\theta-\widehat{\operatorname{bias}}_{\mathrm{boot}} \ \pm\ z_{1-\alpha/2}\widehat{\operatorname{SE}}_{\mathrm{boot}} \tag{3}$$
式(3)は、バイアス補正した中心から、ブートストラップ標準誤差の一定倍だけ左右へ広げています。$z_p$は標準正規分布の$p$分位点で、95%区間では約1.96です。
今回の平均0.805、バイアス推定約−0.000692、標準誤差約0.174896を代入すると、0.4629〜1.1485となり、関数の出力と一致しました。ライブラリや引数によって分散の渡し方が異なることもあるため、「normal」と表示された区間がどの中心と分散を使っているかを確認します。
各標本の標準誤差で割る方法
ブートストラップ標本によって、推定値だけでなく標準誤差も変わります。この違いを使うのが、studentized bootstrap、あるいはbootstrap-tと呼ばれる方法です。
$$T_b^* = \frac{\widehat\theta_b^*-\widehat\theta} {\widehat{\operatorname{SE}}_b^*} \tag{4}$$
式(4)は、各再抽出標本での推定値のズレを、その標本自身の標準誤差で割っています。$\widehat{\operatorname{SE}}_b^*$は$b$番目の再抽出標本から得た標準誤差です。
この$T_b^*$の$p$分位点を$t_p^*$、元標本での標準誤差を$\widehat{\operatorname{SE}}$とすると、区間は次のようになります。
$$\mathrm{CI}_{\mathrm{stud}} = \left[ \widehat\theta-t_{1-\alpha/2}^*\widehat{\operatorname{SE}},\ \widehat\theta-t_{\alpha/2}^*\widehat{\operatorname{SE}} \right] \tag{5}$$
式(5)は、標準化した誤差の上下の分位点を使い、元標本の尺度へ戻して区間を作っています。上側の誤差分位点が区間の下限に入り、下側の誤差分位点が上限に入る点に注意します。
今回の平均変化では、各標本内の変化の不偏分散を20で割って平均の分散を求められるため、標準誤差を得るための内側のブートストラップは不要です。そして、元標本の標準誤差約0.179257を使って、区間0.4569〜1.2116を得ました。
「bootstrap-t」という名前でも、必ず既知の$t$分布表の分位点を使うわけではありません。今回のStudentized区間では再抽出から得た$T_b^*$の分位点を使い、対応のある通常の$t$区間では自由度19の$t$分布を使っています。
推定量の標準誤差を式やモデル出力から得られない場合には、各再抽出標本の内側で再びブートストラップする方法もあります。ただし、外側5,000回と内側200回なら最大で約100万回の内側の計算が必要になり、計算負荷と内側の標準誤差の安定性を別に確認しなければなりません。
また、再抽出標本の標準誤差が0へ近づくと、式(4)の値が極端になります。この不安定さは方法を複雑にしても消えず、標準誤差が安定して求まるかという問題として残ります。
BCa区間の補正
パーセンタイル区間は2.5%点と97.5%点をそのまま使いますが、BCa区間では、推定値の位置のズレと標準誤差の変わり方を反映して、使う分位点の位置を調整します。BCaはbias-corrected and acceleratedの略称です。
Efronの1987年の論文は、この調整と、成立条件の下での高次の近似精度を論じています。ただし、有限標本でどこまで正確になるかは、推定量の性質と標本サイズによって変わります。
位置のズレを表す補正値
元の推定値より小さいブートストラップ推定値の割合から、位置の補正値$\widehat z_0$を求めます。
$$\widehat z_0 = \Phi^{-1}\!\left\{ \frac1B\sum_{b=1}^{B} I(\widehat\theta_b^*<\widehat\theta) \right\} \tag{6}$$
式(6)は、元の推定値がブートストラップ分布のどの位置にあるかを、標準正規分布の尺度へ移しています。$\Phi$は標準正規分布の分布関数で、ちょうど半数が元の推定値より小さければ$\widehat z_0=0$になります。
これは、前の記事で求めたブートストラップ平均と元の推定値の差を、標準誤差で割ったものではありません。同じ「バイアス補正」という言葉でも、BCaで使う$\widehat z_0$と、ブートストラップ平均からのバイアス推定は計算が異なります。
今回の$\widehat z_0$は約0.0276でした。また、同じ値が多数並ぶ離散的な推定量では、同値をどう扱うかで結果が動く場合があるため、滑らかな分布を前提とする説明をそのまま移すのは注意が必要です。
加速度と対象者ごとの影響
加速度$\widehat a$は、推定値が変わった時に標準誤差がどう変わるかを調整する量で、通常の独立な対象者の設定では、対象者を一人ずつ除いたジャックナイフ推定値から近似できます。
一人除いた推定値を$\widehat\theta_{(-i)}$、その平均を$\overline\theta_{(-\cdot)}$とすると、よく使われる形は次のとおりです。
$$\widehat a= \frac{ \sum_{i=1}^{n} (\overline\theta_{(-\cdot)}-\widehat\theta_{(-i)})^3 }{ 6\left\{ \sum_{i=1}^{n} (\overline\theta_{(-\cdot)}-\widehat\theta_{(-i)})^2 \right\}^{3/2} } \tag{7}$$
式(7)は、一人ずつ除いた時の推定値の変化について、非対称な影響の出方を要約しています。ある一人を除くと推定値が大きく動く場合には、その対象者の影響が加速度の推定にも入ります。
今回の平均変化では$\widehat a$は約0.0137でした。ただし、加速度は、再抽出分布の歪度をそのまま代入する量ではなく、対象者単位の影響や推定量の性質に関わります。また、施設を独立な抽出単位とする場合、一人ずつ除くジャックナイフを使えるかは、どこを独立と考えたかで決まります。
使用する分位点の移動
両側の確率$p=\alpha/2$と$p=1-\alpha/2$について、BCaで使う確率$\widetilde p$を求めます。
$$\widetilde p= \Phi\!\left[ \widehat z_0+ \frac{\widehat z_0+z_p} {1-\widehat a(\widehat z_0+z_p)} \right] \tag{8}$$
式(8)は、位置の補正値と加速度を使い、元の$p$分位点を別の分位点へ移しています。BCa区間は、この調整後の下側と上側の確率に対応するブートストラップ分位点を使います。
今回の調整後の確率は約0.0318と0.9808で、通常の0.025と0.975より右側へ移りました。その結果、BCa区間は0.4925〜1.1715となっています。この右方向への移動は今回のデータでの結果で、BCaなら上方へ補正されるという一般的な性質ではありません。
また、補正後の確率が0や1へ近づきすぎると、少数の再抽出値へ区間端点が依存します。関数が極端な順序統計量について警告した場合には、反復回数だけでなく、推定量の離散性、影響の大きい対象者、加速度の不安定さも確認します。
信頼区間を選ぶ時の判断
今回の95%区間を並べると、正規近似、基本、パーセンタイル、BCa、Studentizedで上下限が異なります。また、対応のある$t$区間は0.4298〜1.1802でしたが、一つの標本で短い区間が得られたことから、方法の優劣は判定できません。
信頼区間の性能として確認したいのは、同じ研究を繰り返した時に真値を含む割合、つまり被覆率と、区間の長さや左右の偏りです。そして、真値が分からない今回の一つのデータだけでは、その被覆率を測れません。
平均と複雑な指標の違い
私なら、今回のような平均変化では、まず通常の標準誤差と$t$区間を計算し、ブートストラップの結果がどこで異なるかを確認します。元の差の分布、外れ値、一人ずつ除いた場合の動きを見れば、区間の違いを説明する根拠が増えるためです。
滑らかな推定量で、再抽出分布の非対称性を反映したい場合には、BCaやStudentizedが候補になります。ただし、BCaでは補正値と加速度、Studentizedでは各標本の標準誤差が安定して求まることを確認します。パーセンタイル区間を比較用に添える場合も、主に使う方法は結果を見る前に理由とともに決めます。
反対に、最もp値が小さくなる方法や、0を含まない方法を後から選ぶことは推奨しません。区間の作り方の選択が結果に依存すると、報告した95%という水準だけでは解析手順全体の性質を説明できなくなります。
平均の区間と予測区間
平均変化の95%信頼区間は、集団の平均をどこまで絞れるかに関する区間です。新しい一人の変化を予測する区間には、平均の推定誤差だけでなく、個人間の変動も必要になります。
したがって、10,000個の平均を並べて、新しい患者の予測区間として使うことはできません。予測を目的にするなら、モデルの推定の揺れに加え、アウトカムの条件付き分布からの変動も扱う方法が必要です。
また、Kaplan-Meier推定量のような生存曲線で各時点に付けた95%区間と、曲線全体を95%で覆う同時信頼帯は、別のものです。複数時点での最大のズレを統計量にするなど、曲線全体の不確実性に対応する構成を別に考えます。
反復回数と計算の安定性
20人を固定してブートストラップ回数を1,000回から10,000回へ増やすと、主に小さくなるのは乱数計算による誤差です。経験分布が母集団分布をどれくらい表せているかという問題は、回数だけでは改善しません。
再抽出推定値がほぼ正規分布に従い、反復が独立とみなせる場合、ブートストラップ標準誤差の相対的なモンテカルロ誤差は、おおむね次の形で捉えられます。
$$\frac{ \operatorname{MCSE} (\widehat{\operatorname{SE}}_{\mathrm{boot}}) }{ \widehat{\operatorname{SE}}_{\mathrm{boot}} } \approx \frac1{\sqrt{2(B-1)}} \tag{9}$$
式(9)は、標準誤差の推定そのものが乱数計算で揺れる相対的な大きさを表しています。再抽出分布が強く歪む場合や裾が重い場合には、この正規近似による目安をそのまま使えません。
この計算では、$B=200$で約5.0%、1,000で約2.2%、10,000で約0.7%です。ただし、信頼区間の端点や極端な分位点に、同じ相対誤差が当てはまるわけではありません。
95%区間の裾に残る個数
2.5%点を求める場合、200回の再抽出では、その点より下の値は期待値で約5個しかありません。1,000回なら約25個、5,000回なら約125個、10,000回なら約250個です。BCaで下側確率が0.002へ移るなら、10,000回でも期待値は約20個になります。
研究の最終的な95%区間を作る場合、計算負荷が許せば、5,000〜10,000回程度から始めて乱数を変えた時の端点の動きを確認する進め方を推奨します。これは本稿での作業上の提案で、5,000回を超えれば精度が保証されるという基準ではありません。
端点の揺れをより細かく考えるなら、分位点$p$での分布関数のモンテカルロ標準誤差は、独立反復の下でおおむね$\sqrt{p(1-p)/B}$です。$p=0.025$、$B=10,000$では約0.00156ですが、これは確率の尺度での揺れで、CD4変化の尺度での端点誤差ではありません。
値の尺度へ移す時には、分位点付近にどれくらい推定値が集まっているかも関わります。裾の密度が低いと、確率のわずかなズレが端点の大きな移動につながるため、計算回数だけでなく、端点付近の分布を確認します。
回数を増やした時の端点
同じ10,000回の結果について、最初の200回、500回、1,000回と使う範囲を広げ、2.5%点と97.5%点の推移を描きます。別々の標本や独立した六つの実験を比べているのではなく、同じ計算結果の先頭部分を順に増やしています。
checkpoints <- c(200, 500, 1000, 2000, 5000, 10000)
limits <- sapply(checkpoints, FUN = function(repetitions)
quantile(resampled$t[seq_len(repetitions), 1], probs = c(0.025, 0.975)))
png("bootstrap-repetitions.png", width = 1100, height = 688, res = 200, type = "cairo")
plot(checkpoints, y = limits[1, ], type = "b", log = "x", ylim = range(limits),
xlab = "Bootstrap replicates", ylab = "Percentile interval endpoints", pch = 19)
lines(checkpoints, y = limits[2, ], type = "b", pch = 17, col = "blue3")
legend("right", c("97.5% endpoint", "2.5% endpoint"), col = c("blue3", "black"), lty = 1, pch = c(17, 19), bty = "n")
dev.off()
この描画での上下限は、200回で約0.4766〜1.0976、10,000回で約0.4765〜1.1505でした。下限だけを見ると変化が小さく見えますが、途中では上下し、上限とは動きの向きも幅も違います。
なお、図1ではRのquantile関数の既定の補間を使っています。一方、前の記事で実行したboot.ci関数は分位点を正規分位点の尺度で補間するため、パーセンタイル区間0.4750〜1.1510とは端数が一致しません。乱数計算の違いだけでなく、分位点の定義や補間方法の違いも区別します。
この図が落ち着いて見えても、同じ乱数系列の中で安定していることを見ているだけです。最終解析では独立の乱数シードでも計算し、主張に影響する程度の差が残らないかを確認します。また、シード、反復数、ソフトウェアの版を記録して、どの計算から本文の数値を得たかを追えるようにします。
区間を選ぶ順序を先に決める
ここまでで五つの区間を並べましたが、今回の一つのデータからは、どれが最も被覆率の良い方法かは決まりません。決められるのは、推定量の性質と、補正値や標準誤差が安定して求まるかどうかという、計算前に確認できる条件です。
次の記事では、この区間の作り方を回帰分析、階層のあるデータ、時系列へ広げます。行を選ぶのか残差を選ぶのか、施設を選ぶのか施設内で選ぶのかで、同じ名前の方法が別の標本変動を再現することになります。
参考文献
Efron B. Better bootstrap confidence intervals. Journal of the American Statistical Association. 1987;82(397):171-185.
DiCiccio TJ, Efron B. Bootstrap confidence intervals (with discussion). Statistical Science. 1996;11(3):189-228.
Efron B. Bootstrap methods: Another look at the jackknife. The Annals of Statistics. 1979;7(1):1-26.
Davison AC, Hinkley DV. Bootstrap Methods and Their Application. Cambridge University Press; 1997.
