ここまでの記事では、推定値の標準誤差と信頼区間を求めるためにブートストラップを使ってきました。ただし、再抽出の使い道はそれだけではありません。予測モデルを作ったデータでそのモデルの性能を測ると、新しい対象者に使った場合より良い値が出ますが、その差そのものを推定する道具にもなります。
さらに、欠測を補う多重代入と組み合わせる場面、追跡時間とイベント指標の組を選び直す生存時間解析、傾向スコアの推定を含む複数段階の手順では、どこまでを繰り返すかという指定が結果を分けます。
この記事では、楽観性補正による内部検証、MICEとの組み合わせ、生存時間解析での再抽出、通常のブートストラップが成り立たない推定量、そして論文に何を書くかまでを扱います。近似が崩れる条件は、最大値とまれなイベントで扱いました。
予測モデルの内部検証
予測モデルを作ったデータで、そのモデルのAUCや予測誤差を測ると、新しい対象者に使った場合より良い性能が得られやすくなります。説明変数やハイパーパラメータを同じデータに合わせて選んだ分、そのデータ固有の特徴にも適合するためです。
この楽観性について、Steyerbergらは、心筋梗塞後の30日死亡を扱うロジスティック回帰モデルなどを使い、分割標本、交差検証、ブートストラップによる内部検証を比較しました。彼らの研究では、単純な標本分割の非効率性と、ブートストラップによる内部検証の利用価値が論じられています。
ここで推定したいのは、係数の標準誤差だけではなく、元データで測った性能がどれくらい楽観的かという量です。
楽観性を差し引く手順
元データを$D$、$b$番目の再抽出データを$D_b^*$、それぞれで作った予測モデルを$\widehat f$と$\widehat f_b^*$とします。また、$A(f;D)$を、値が大きいほど良い性能指標とすると、楽観性補正の基本形は次のように書けます。
$$O_b= A(\widehat f_b^*;D_b^*)-A(\widehat f_b^*;D), \qquad A_{\mathrm{corrected}} = A(\widehat f;D)-\frac1B\sum_{b=1}^{B}O_b \tag{1}$$
式(1)は、再抽出データで作ったモデルについて、学習に使った再抽出データと元データでの性能差を求め、その平均を元モデルの見かけの性能から引いています。$O_b$は$b$番目の楽観性の推定値です。
ここで、この計算では元データが再抽出データと重なっているので、独立した外部検証標本として扱う条件は満たしていません。同じモデルの二つのデータでの性能差を利用し、開発手順に伴う楽観性を近似しています。
なお、予測誤差や損失のように小さいほど良い指標では、差の向きに合わせて補正を定義します。AUC用の式を符号まで確認せずにコピーすると、誤差をさらに楽観的な方向へ補正してしまう場合があります。
変数選択と調整をやり直す範囲
元の解析で説明変数を選び、連続変数の変換を決め、正則化の強さを調整していたなら、その開発手順を再抽出の中でも繰り返します。元データで選ばれた変数を固定して係数だけ推定し直す方法では、変数選択による楽観性を十分に扱えません。
そのため、標準化、欠測処理、変数選択、内部の交差検証による調整なども、どのデータから学習した処理かを管理します。そして、各再抽出標本で学習した変換を、評価先のデータへ同じように適用します。評価対象のアウトカムを、実運用では利用できない予測入力として混ぜないことも必要です。
この問題は、Coxモデルの診断と変数選択で扱うモデル選択ともつながります。ただし、選択操作を全て繰り返せば、選択後の係数に対する通常のブートストラップ区間まで必ず正当化できる、という意味ではありません。
固定したモデルの評価との違い
既に固定したモデルを独立の外部検証データで評価し、そのAUCの標準誤差を求める場合には、検証対象者を選び直してAUCを計算する方法が候補になります。この場合、元の開発モデルを毎回学習し直す必要があるかは、固定済みモデルの性能を評価するのか、モデル開発の揺れまで含めるのかによって異なります。
また、二つの固定モデルを同じ対象者で比べる場合には、両モデルを同じ再抽出対象者で評価し、AUCなどの差を毎回計算します。モデルごとに別々の対象者を選ぶと、同じ検証対象者に基づく性能差の相関を失います。
そして、ブートストラップによる内部検証は、別施設、別時期、異なる診療環境での外部検証を代替しません。元の対象者に含まれない集団への移行で起きる変化は、その対象者を選び直すだけでは分かりません。
バギングとの違い
バギングは、再抽出標本ごとに予測モデルを作り、それらの予測を平均するなどして、一つの予測手順として利用する方法です。再抽出による学習を使う点は共通しても、標準誤差の推定や楽観性補正とは目的が異なります。
さらに、各モデルの予測値のばらつきをそのまま並べても、それだけで校正された予測区間にはなりません。モデルの推定の揺れと、新しいアウトカム自体の変動をどのように扱うかが別に必要です。
低出生体重モデルで楽観性を測る
MASSパッケージのbirthwtデータセットを使い、低出生体重の有無をアウトカムとするロジスティック回帰で楽観性を実測します。189人について、年齢、最終月経時体重、喫煙、早産歴、高血圧、子宮過敏、受診回数を説明変数に入れたモデルです。AUCは順位和から計算し、500回の再抽出ごとにモデルを当てはめ直しています。
set.seed(42)
births <- MASS::birthwt
auc <- function(outcome, score) {
ranks <- rank(score)
positives <- sum(outcome == 1)
(sum(ranks[outcome == 1]) - positives * (positives + 1) / 2) /
(positives * sum(outcome == 0))
}
develop <- function(data) glm(low ~ age + lwt + smoke + ptl + ht + ui + ftv,
family = binomial, data = data)
apparent <- auc(births$low, predict(develop(births), type = "link"))
optimism <- replicate(500, {
drawn <- births[sample.int(nrow(births), nrow(births), replace = TRUE), ]
refit <- develop(drawn)
auc(drawn$low, predict(refit, type = "link")) -
auc(births$low, predict(refit, newdata = births, type = "link"))
})
round(c(apparent = apparent, corrected = apparent - mean(optimism)), 4)#> apparent corrected
#> 0.7383 0.6934元のデータで測ったAUCは0.7383で、楽観性の平均は0.0449、差し引いた後は0.6934でした。189人に対して説明変数が7つという規模で、見かけの性能が0.04ほど高く出ていたことになります。
ここで注意したいのは、0.6934を外部検証で得た性能として報告できない点です。評価に使った189人は開発にも使った189人なので、別の施設や別の時期の患者でどうなるかは、この計算からは分かりません。
欠測データとMICEの組み合わせ
欠測を含むデータでは、誰が標本に入るかによる揺れに加えて、観測されなかった値についての不確実性があります。ブートストラップと多重代入は、この異なる不確実性を扱うため、組み合わせる順序と統合方法を明らかにする必要があります。
ここで、欠測を1回だけ補ったデータを完成済みの実測データとして扱い、その行をブートストラップしても、多重代入で表すべき補完の不確実性を十分に反映できません。また、30組の補完データを縦に積んで600人分の独立データとして抽出すれば、元の20人を30倍の人数として扱うことになってしまいます。
外側で抽出し、内側で代入する方法
一つの方式では、欠測を残した元データから対象者を抽出し、その再抽出データの中でMICEによる多重代入を行います。そして、補完データごとの推定値をその再抽出内で統合し、外側のブートストラップ間の分布を評価します。
ただし、外側の回数を$B$、各再抽出標本の補完データ数を$m$とすると、単純な構成では$B\times m$回の解析が必要になります。$B=5,000$、$m=20$なら100,000回で、MICEの連鎖を更新する反復回数はさらに別の設定です。
もっとも、全ての$B\times m$個の推定値を独立なブートストラップ反復として一列へ並べればよい、という意味ではありません。同じ外側標本を共有する補完結果の階層と、有限個の補完によるモンテカルロ誤差を扱う必要があります。
順序だけでは決まらない妥当性
BartlettとHughesは、補完モデルと解析モデルの不整合やモデルの誤指定の下で、ブートストラップと多重代入を組み合わせる方法の違いを調べています。特定の「抽出後に代入する」方式が有効な条件を扱っていますが、そこで示されているのは条件付きの成立であって、順序の選択が欠測モデルの誤りまで補うという結果ではありません。
もう一つ、多重代入の各データでブートストラップにより完全データ分散を推定し、Rubinの統合則へ渡す方法もありますが、その成立条件は補完と解析の関係に依存します。代入後に抽出する方法を全て同じと扱ったり、常に不適と説明したりすることも正確ではありません。
そして、欠測がMNARなのにMARのモデルだけを使っている場合や、推定したい効果そのものが識別できない場合、ブートストラップはその偏りを自動的に除きません。推定値が目的の量へ近づく条件と、分散を正しく評価できる条件は区別します。
欠測アウトカムを予測因子として使うこと、アウトカム自体を補うこと、補完値を最終解析へ含めることも、再抽出とは別に決めます。欠測機構については、欠測データと多重代入と合わせて、どの不確実性をどの段階で扱っているかを確認します。
生存時間解析での再抽出
生存時間解析では、一人分の情報は観測追跡時間、イベント指標、共変量の組です。独立した対象者を想定する場合には、この組をまとめて選び、各再抽出標本でKaplan–Meier推定やCox回帰などを計算し直します。
ここで、追跡時間だけを抽出し、イベント指標を別の人から選ぶと、誰がどこまでイベントなく追跡されたかという情報が崩れます。また、イベントが観測された人だけを選ぶ方法では、打ち切りを持つ研究を再現できません。
打ち切りを残した再計算
通常の対象者単位の再抽出では、打ち切り例も観測された時間とイベント指標の組のまま使います。打ち切り後の未知のイベント時刻を新しく確定する処理は不要で、打ち切りを扱う解析を各標本でやり直します。
一方、モデルに基づいて生存時間や打ち切り時間を生成する方式もあり、bootパッケージのcensboot関数は複数の方式を扱っています。ただし、イベント分布と打ち切り分布のどの関係を仮定するかを指定する必要があります。
そのうえ、情報的打ち切りに対する重みを推定しているなら、重みの推定も含めてどこまで繰り返すかを決めます。再抽出を加えても、無視できない打ち切りの仕組みは元のまま残ります。
複数行と推定できない生存中央値
同じ人がstart–stop形式で複数行を持つ場合には、行を独立に選ぶのではなく、人を選んでその人の全区間を持ち込みます。そして、同じ人のコピーを識別できるIDを付け、各再抽出標本でリスク集合を作り直します。
ここで注意が必要なのは、元標本で生存中央値が求まっていても、再抽出標本では生存曲線が0.5まで下がらず、中央値を推定できない場合があります。このような標本を黙って除いて区間を作ると、中央値を推定できた標本に条件を付けた分布へ変わってしまいます。
そのため、解析前に、中央値を支える追跡とイベントがあるかを確認し、研究の問いに合えば、追跡で支えられる時点の生存確率や、期間を事前に定めたRMSTなども考えます。ただし、推定しやすい指標へ事後的に変更した場合は、その経緯も区別して報告します。
論文に残す抽出と計算の情報
「ブートストラップを10,000回行いました」という一文だけでは、読者は何を抽出し、何を計算し直し、どの区間を報告したかを再現できません。また、対象者の標本変動を調べたのか、予測モデルの楽観性を補正したのかも分からないままです。
そこで、方法には、推定対象、抽出単位、復元の有無、層別やクラスター構造、再抽出ごとの計算手順、反復数、信頼水準、区間の種類、乱数シード、ソフトウェアの版を書きます。そして、推定に失敗した反復の件数と扱い、計算回数や乱数を変えて確認した内容を続けます。
今回の方法の記述例
今回のcd4データでの計算を文章にすると、次のようになります。
開始時から1年後のCD4数の平均変化を推定するため、開始時と1年後の測定値の組を保った対象者単位のノンパラメトリック・ブートストラップを行いました。20人から毎回20人分を復元抽出し、平均変化を10,000回計算しました。また、平均変化の標準誤差は再抽出推定値の標準偏差から求め、BCa法による両側95%信頼区間を計算しました。計算にはR4.4.2とbootパッケージ1.3-31を使用し、乱数シードを42としました。
なお、ここでBCaだけを書くのは、BCaを報告する方法の例です。このシリーズでは比較のために他の区間も計算しているので、実際に複数の結果を報告する論文では、その比較が事前の計画か追加解析かも明らかにします。
もう一点、反復回数を増やしながら端点の推移を見た図は、同じ乱数系列で使う範囲を変えた診断なので、独立した複数シードによる安定性確認を行ったと書き換えることはできません。診断の記述は、実際に行った範囲に合わせます。
今回の結果の記述例
結果は、平均変化と標準誤差、区間、元の人数を対応させて書けます。
20人におけるCD4数の平均変化は0.805で、対象者単位のブートストラップ標準誤差は0.175、BCa法による95%信頼区間は0.493〜1.172でした。値は元データの100個単位で記載しています。また、対応のある$t$区間は0.430〜1.180でしたが、これは比較のために求めた別方式の区間です。
そして、この結果から書けるのは、観測された平均変化と、置いた抽出・解析の仮定の下での推定の不確実性です。対照群のないデータから薬剤の因果効果を確定したり、別施設での再現性まで確認したりした結果としては報告しません。
自分のデータで判断が難しくなりやすいのは、回数の指定よりも、一人、施設、時点のどれを選び直すか、そして補完や変数選択までどこを推定し直すかという部分です。こうした解析手順と推定対象の対応については、Dr.データサイエンスでご相談いただけます。
再抽出の記録を残す
ブートストラップは、標準誤差を出す道具としても、開発手順の楽観性を測る道具としても使えます。どちらのつもりで回したのかは、出力された数字を見ても区別が付きません。
そのため、推定対象、抽出単位、繰り返した処理の範囲、反復数、区間の種類、乱数シードまでを方法へ書いておくことが、計算そのものと同じだけの意味を持ちます。
参考文献
Steyerberg EW, Harrell FE Jr, Borsboom GJJ, Eijkemans MJC, Vergouwe Y, Habbema JDF. Internal validation of predictive models: Efficiency of some procedures for logistic regression analysis. Journal of Clinical Epidemiology. 2001;54(8):774-781.
Bartlett JW, Hughes RA. Bootstrap inference for multiple imputation under uncongeniality and misspecification. Statistical Methods in Medical Research. 2020;29(12):3533-3546.
Abadie A, Imbens GW. On the failure of the bootstrap for matching estimators. Econometrica. 2008;76(6):1537-1557.
Efron B, Tibshirani RJ. An Introduction to the Bootstrap. Chapman and Hall; 1993.
Davison AC, Hinkley DV. Bootstrap Methods and Their Application. Cambridge University Press; 1997.
