- パッケージとデータセットの準備
- 事前準備
- bootstraps( )を使ってブートストラップサンプルを作成する
- 関数を定義する
- optimismを計算する
- optimism補正性能指標を計算する
- おわりに
- 参考資料
予測性能には、識別能と較正能があって、内的妥当性検証の1つとしてoptimism補正がある。
以前、ROC曲線下面積（AUROC）のoptimism補正についてはまとめてみた。
また別の記事で、予測の較正能についてもまとめてみた。
今回は、予測性能指標をまとめてoptimism補正してしまおうとおもう。これにはrmsパッケージのvalidate( )を使えば簡単そうだが、サンプルサイズが大きい症例（>5000）でやってみたらエラーになってしまった。 そこで、purrrパッケージmap( )の練習を兼ねて自力でやってみる。
パッケージとデータセットの準備
pROCパッケージのaSAHデータセットを使う。関数として必要なのはtidymodelsパッケージとrmsパッケージ。
library(tidymodels) library(rms) library(pROC) data(aSAH)
変数outcomeをアウトカムに使う。カテゴリーが"Good"と"Poor"なので、それぞれ1と0に、変数名をyに変更した。
sample_o <- aSAH %>% mutate(y = if_else(outcome == "Good", 1, 0))
事前準備
この後のコードは自作なので、必要になるオブジェクトを事前に指定しておく。
sample_o：元データ。既に指定済み。model_o：元データで開発モデル。ここでは予測因子としてageを使う。後で使うのはformulaの部分だけ。B：ブートストラップの回数。本当は200回以上必要だが、ここでは内容確認がしやすいようにB=10にした。
model_o <- glm(y ~ age, data = sample_o, family = binomial(link="logit")) B <- 10 set.seed(1234)
bootstraps( )を使ってブートストラップサンプルを作成する
tidymodelsパッケージのbootstraps( )を使ってB=10個のブートストラップサンプルと末尾に元データがくっついたデータリストを作成する。
d_boot_app <- bootstraps(sample_o, times = B, apparent = TRUE)
内容はこんな感じ。
> d_boot_app # Bootstrap sampling with apparent sample # A tibble: 11 × 2 splits id <list> <chr> 1 <split [113/43]> Bootstrap01 2 <split [113/41]> Bootstrap02 3 <split [113/38]> Bootstrap03 4 <split [113/49]> Bootstrap04 5 <split [113/40]> Bootstrap05 6 <split [113/45]> Bootstrap06 7 <split [113/43]> Bootstrap07 8 <split [113/49]> Bootstrap08 9 <split [113/44]> Bootstrap09 10 <split [113/43]> Bootstrap10 11 <split [113/113]> Apparent
関数を定義する
上記の各データセットに対して適用する関数を定義する。
fun <- function(split){ sample_b <- analysis(split) model_b <- glm(formula(model_o), data = sample_b, family=binomial(link="logit")) fv_boot <- predict.glm(model_b, type="response") fv_test <- predict.glm(model_b, newdata=sample_o, type="response") pfm_boot <- val.prob(p = fv_boot, y = sample_b$y, pl=FALSE) pfm_test <- val.prob(p = fv_test, y = sample_o$y, pl=FALSE) result <- list(model_b = model_b, pfm_boot = pfm_boot, pfm_test = pfm_test) return(result) }
中身は以下のとおり。
sample_b：各データセット（B個のブートストラップデータ + 元データ）。この関数の唯一の引数になる。model_b：引数にとったデータセットをもとに作成されたモデル。前述のmodel_oと同じ式を指定する。fv_boot：model_bにおける予測値（y=1の予測確率）fv_test：元データに対してmodel_bを使って計算した予測値pfm_boot：各データセットを使って開発したモデル（model_b）の予測性能指標（= bootstrap performance*1）pfm_test：元データに対して計算されたmodel_bの予測性能指標（= test performance*2）
model_b, pfm_boot, pfm_testの3つを結果として返す。
ちなみに、pfm_bootおよびpfm_testの最後（=B+1番目）はともに「元データを使って開発したモデル」の「元データに対して計算された予測性能指標」なので一致する（=apparent performance）。
optimismを計算する
上記で作成した関数をデータセットのリストに対して適用する。
d_boot_app <- d_boot_app %>% mutate(result = map(splits, fun), model_b = map(result, ~.$model_b), pfm_boot = map(result, ~.$pfm_boot), pfm_test = map(result, ~.$pfm_test), opt = map2(pfm_boot, pfm_test, ~.x-.y))
前述のとおり、resultは3つの結果のリストなので、3つの結果を取り出して元の名前で保持した（チルダ式（~.）については別の機会に勉強します）。
pfm_***には色々な指標がまとめて含まれている。
optでは対応する指標についてoptimismの計算に使う差を計算している（ここらへんの詳細は、較正プロット（calibration plot）で較正能を評価する（rmsパッケージvar.prob関数）を参照のこと）。
出来上がるものは次のとおり。
> d_boot_app # Bootstrap sampling with apparent sample # A tibble: 11 × 7 splits id result model_b pfm_boot pfm_test opt <list> <chr> <list> <list> <list> <list> <list> 1 <split [113/43]> Bootstrap01 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 2 <split [113/41]> Bootstrap02 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 3 <split [113/38]> Bootstrap03 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 4 <split [113/49]> Bootstrap04 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 5 <split [113/40]> Bootstrap05 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 6 <split [113/45]> Bootstrap06 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 7 <split [113/43]> Bootstrap07 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 8 <split [113/49]> Bootstrap08 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 9 <split [113/44]> Bootstrap09 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 10 <split [113/43]> Bootstrap10 <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl> 11 <split [113/113]> Apparent <named list [3]> <glm> <dbl [18]> <dbl [18]> <dbl>
optには18種類の予測指標に対して計算されたbootstrap performanceとtest performanceの差が含まれる。
元データ（idが"Apparent"のもの）はpfm_bootとpfm_testが同じもなのので、optは（計算可能であれば）0になっている。
optimism補正性能指標を計算する
まず、opt_resultとして単に計算対象となる部分を取り出した。
次に、データセットを超えて各指標のoptimismの平均を取り、optimism（opt_m）とした
（do.call( )について末尾リンク参照）。
opt_result <- d_boot_app %>% filter(id != "Apparent") %>% select(opt) opt_m <- do.call(rbind, unlist(opt_result, recursive = FALSE)) %>% as.data.frame() %>% summarize_all(mean)
計算結果は以下のとおり。
> opt_m Dxy C (ROC) R2 D D:Chi-sq D:p U U:Chi-sq 1 0.0293955 0.01469775 0.06099864 0.04731428 5.346513 NA -0.02343513 -2.648169 U:p Q Brier Intercept Slope Emax E90 1 0.5989901 0.0707494 -0.01738962 0.2270676 -0.147162 0.01542379 -0.006755463 Eavg S:z S:p 1 -0.003548597 -1.000607 0.6156088
Apparent performanceは以下のコードで取り出した。これがapparent performanceになるのは前述のとおり。
pfm_app <- d_boot_app %>% filter(id=="Apparent") %>% select(pfm_boot) %>% unlist()
> pfm_app pfm_boot.Dxy pfm_boot.C (ROC) pfm_boot.R2 pfm_boot.D 2.300136e-01 6.150068e-01 6.050485e-02 3.633655e-02 pfm_boot.D:Chi-sq pfm_boot.D:p pfm_boot.U pfm_boot.U:Chi-sq 5.106030e+00 NA -1.769912e-02 2.842171e-14 pfm_boot.U:p pfm_boot.Q pfm_boot.Brier pfm_boot.Intercept 1.000000e+00 5.403566e-02 2.212515e-01 1.320229e-08 pfm_boot.Slope pfm_boot.Emax pfm_boot.E90 pfm_boot.Eavg 9.999999e-01 9.125128e-02 4.154254e-02 2.322127e-02 pfm_boot.S:z pfm_boot.S:p 2.643858e-02 9.789075e-01
結果の体裁は悪いけど、apparent AUROCは0.615くらいみたい。
最後にapparent performanceからoptimismを引いて、optimism-corrected performance（pfm_oc）を得る。
pfm_oc <- pfm_app - opt_m
> pfm_oc Dxy C (ROC) R2 D D:Chi-sq D:p U U:Chi-sq 1 0.200618 0.600309 -0.000493784 -0.01097773 -0.2404836 NA 0.005736012 2.648169 U:p Q Brier Intercept Slope Emax E90 Eavg 1 0.4010099 -0.01671374 0.2386411 -0.2270676 1.147162 0.07582749 0.048298 0.02676987 S:z S:p 1 1.027046 0.3632988
optimism-corrected AUROCは0.600と少しだけ下がった。
おわりに
- 次は、ブートストラップサンプルごとに計算したROC曲線や較正プロットを全て重ねて描くコードを考えてみる予定。
- 「猫なで声」の由来を実感している今日この頃です。
参考資料
- do.call( )という見慣れない関数を調べるのに参考にさせてもらいました。
- このリンクでの質問にあるように、本当はdo.call( )を使わずにbind_rows( )で行きたかったけど、何故かうまくいかず...
- rmsパッケージの使い方がひととおり書かれているので参考にしてます。
