Optimismとは
予測モデルの性能を高めようとすればするほど、予測モデルを開発するために使用するデータに当てはまりすぎてしまって、新規のデータに対しては役に立たなくなってしまうことがある。これを過剰学習（overfitting）、あるいは過学習（overtraining）という。
予測モデル開発に使用したデータセットを使った予測性能指標は、この過剰学習の分だけ楽観的な数値になってしまうので、その"optimism"*1を差っ引いて補正した指標（optimism-corrected performance, OCP）を算出しなさい、と言われている。内的妥当性検証の1つ。
代表的な予測性能指標であるROC曲線下面積（AUROC）については、以前1回まとめたことがあります。
Optimism補正予測性能の求め方
元データから再サンプリングを繰り返すbootstrap法を用いる方法が一般的。「どのデータセットを使って開発したモデルか」と「どのデータセットを使って算出した性能指標か」に注意しておかないと混乱しがち。
手順は次のとおり。
- 元データ
を使ってモデルを開発して（
）、
）を計算する。
- 元データ
を作成する。
- ブートストラップサンプル
）、
）を計算する。
- 元データ
）を計算する。
- このブートスラップサンプルにおけるoptimism
=
- 上記を十分な回数（200回以上）繰り返し、
- Optimism補正予測性能（OPC） =
ブートストラップについても以前まとめたことがあるので参考までに。
tidymodelsパッケージを使ってoptimism-corrected AUROCを計算する
Rで実装する場合、predbootパッケージが簡単でしたが、何かをアップデートをしてから使えなくなってしまった。開発者のGithubには使用されている関数の定義が非常に読みやすく書かれているので、これを使えば実装できそう。
しかし、tidymodelsパッケージとmap( )の練習のため、自分で書いてみる。
まずは必要なパッケージの読み込み。使用するデータセットは下の記事で使ったものと同じく、pROCパッケージのaSAHデータを使う。
library(tidymodels) library(pROC) data(aSAH) sample_o <- aSAH
後で定義する自作関数の都合上、使用する元データの名前をsample_o、予測したい変数をoutcomeという名前にしておく（aSAHデータでは元々outcomeという名前になっている）。
bootstraps( )を使ってブートストラップサンプルを作成する
繰り返し処理にはforループ文やapply系関数を用いるが、purrrパッケージは繰り返し処理を簡単に行うことができるtidyverse系パッケージ。tidymodelsパッケージはpurrrの要領で、モデル妥当性検証*2のときにデータセット分割およびその解析が便利に行えるもの。
ここではbootstraps( )という関数を使って、B=5個のブートストラップサンプルを作成してみる（実際はB=200以上にする）。
ちなみにパイプ演算子%>%を使って書くことも可。
set.seed(1234) B=5 bt_rsample <- bootstraps(sample_o, times=B) bt_rsample_app <- bootstraps(sample_o, times=B, apparent=TRUE)
引数apparent=TRUEは、末尾に元データを追加するかどうかを指定するためのもの。
ここでは元データを追加したものの中身を確認。
> bt_rsample_app # Bootstrap sampling with apparent sample # A tibble: 6 × 2 splits id <list> <chr> 1 <split [113/45]> Bootstrap1 2 <split [113/43]> Bootstrap2 3 <split [113/49]> Bootstrap3 4 <split [113/44]> Bootstrap4 5 <split [113/43]> Bootstrap5 6 <split [113/113]> Apparent
bootstraps( )で作成されるオブジェクトはrsetオブジェクトと呼ばれ、splitsとidの2つの要素からなる。idを見ると、ここではBootstrap1からBootstrap5まで5個のブートストラップサンプルが抽出されていて、末尾に元データがApparentという名前で保持されているのが分かる。
splitsには、それぞれのデータセットがリストとして保持されている。リストの各要素には、「解析用」と「評価用」のデータセットが含まれていて、例えば[113/45]の2つの数字はそれぞれのデータセットのデータサイズを表している。
ブートストラップ法を行うときは、「解析用」データセットをブートストラップサンプルとして用いる（ブートストラップサンプルのサイズは元データと同じなので間違えないと思う）。
「評価用」はブートストラップサンプルとして抽出されなかった観測値。ブートストラップでは「評価用」は使う機会はない（と思う）。
ブートストラップサンプルを使って妥当性を「評価」するのに、なぜ「解析用」と名付けられているの？と思ったが、bootstraps( )が含まれているrsampleパッケージには他にも
- initial_split( )：データを開発用と検証用に分ける
- vfold_cv( )：k分割交差検証
などのデータ分割を目的とした関数が含まれている。
交差検証（cross-validation）では、例えばデータの90%を使って開発し、残りの10%で検証する。なので、前者を「解析用」、後者を「評価用」と呼ぶことは理解できる。
それぞれのデータセットを取り出したいときは、
- 解析用を取り出したいときはanalysis( )、評価用を取り出したいときはassessment( )を使う。
- 上記の関数にわたす中身は、rsetオブジェクト（bt_rsample）のsplitsリスト（$splits）の何番目かを指定したもの（リストなので二重四角カッコで指定する）。
例えば、1つ目のブートストラップサンプルの解析用データセットを抽出したいときは、下のコードを参照。
> analysis(bt_rsample$splits[[1]]) gos6 outcome gender age wfns s100b ndka 125 5 Good Male 64 1 0.10 46.83 140 5 Good Female 39 1 0.50 6.79 ~~~~~（以下省略）~~~~~
ブートストラップサンプルは元データと同じサイズなので、analysis( )で「解析用」を取り出す。
optimismを計算する自作関数
次に、データセットとモデル式を入力すればoptimismを出力してくれる関数を定義する。 まず、(1) 元データを使ってロジスティック回帰を行い、次に(2) それぞれの症例に対する予測確率を計算してデータセットに追加する。最後にpROCパッケージのroc( )を使ってAUROCを算出する。
ここでは予測因子としてageしか使っていないので、ロジスティック回帰をせずにそのままroc( )に入れたらいいのだが、複数の予測因子を使った予測モデルへ拡張することを考えて、このようにしておく（fvはfitted valueのつもり）。
model_o <- glm(outcome~age, data = sample_o, family = binomial(link="logit")) sample_o$fv_app <- predict(model_o, type="response") c_app <- roc(outcome ~ fv_app, data=sample_o)$auc
ちなみに、apparent performanceは次のとおり。
> c_app Area under the curve: 0.615
次に、先程bootstraps( )で作成したデータセットリストを引数にして、optimism（)を計算する関数を次のように定義する。
fun_opt <- function(split){ d = analysis(split) model_b <- glm(formula(model_o), data=d, family=binomial(link="logit")) fv_boot <- predict(model_b, type="response") fv_test <- predict.glm(model_b, newdata=sample_o, type="response") c_boot <- roc(d$outcome, fv_boot)$auc c_test <- roc(sample_o$outcome, fv_test)$auc delta <- c_boot - c_test return(delta) }
元データと当てはめたモデルは、それぞれsample_o, model_oという名前で保持していることを前提にしている。
これらも引数にして関数に含めてしまえばいいんだろうけど、それならpredboot( )など、既に公開されているものを使ってもらう方がいいだろう。
map( )を使ってデータセットのリストに一括して適用する
purrrパッケージのmap( )はリストに対して繰り返し処理をスムーズに適用できて、（使いこなせれば）とっても便利な関数。
ブートストラップサンプルのsplitsを引数にして、自作関数fun_optを使ってdeltaを算出して、最後に平均を取ってoptimism（opt）を求めている。
Apparent performanceからoptimismを引いたものが、optimism-corrected performance（c_oc）になる。
opt <- bt_rsample %>% mutate(delta = map(splits, fun_opt)) %>% select(delta) %>% unlist() %>% mean() c_oc <- c_app - opt
Apparent AUROCが0.615だったので、optimismを補正すると少し低い値になった。
> c_oc [1] 0.5987854
おわりに
- （私の中で）打ち切りも心配されましたが、なんとか3rdシーズンにこぎつけました。
- 暑さのせいか、フローリングに「落ち猫」していることが増えてきました。
参考資料
- 野間先生のpredbootパッケージのGithub
- rsampleパッケージの関数がわかりやすくまとめられています
- 同じ方のpurrrパッケージ講座。初学者の私でも日本語として読めるように書かれているのがすごい。
- rsampleパッケージのvignette
