登録した患者のデータから予測モデルを作り、同じ患者で予測値と実測値を比べると良い成績が得られても、別の患者では誤差が大きくなり、開発時の成績を再現できない場合があります。また、クロスバリデーションを使っていても、変数選択や欠測値の補完を分割前に済ませると、評価用の情報が学習へ入ります。
そのため、予測性能の評価では分割数だけでなく、評価対象者の情報をどの段階でモデルへ渡したかを確認します。ここでは、見かけの性能、分割ごとの予測、評価対象がモデルか学習手順かという区別、内側で選び外側で評価する構成を、300人のランダムに生成したデータセットで確かめます。
学習に使った患者での成績
血圧や検査値の予測を想定すると、回帰モデルの係数は、学習データの予測誤差を小さくする方向へ推定されます。そして、変数を追加したり、曲線の形を柔軟にしたりすると、対象者に共通する傾向だけでなく、その標本に偶然含まれた変動まで拾うことがあります。
ただし、学習データで計算した誤差が小さくても、新しい患者を同じ精度で予測できるとは限りません。見かけの性能は、同じ情報でモデルを合わせて評価するため楽観的になりやすく、この補正は予測モデルの検証でも扱っています。
ただし、学習時に直接最小化していない指標や偶然の影響まで含めると、単一の標本で見かけの誤差が常に最小になるわけではありません。Stone（1974）が論じたように、データを分けた選択と評価は、機械学習だけでなく回帰モデルにも共通する予測評価の考え方です。
回帰係数の推論との区別
ある検査値の回帰係数と、その信頼区間を知りたい研究と、将来の患者の検査値を予測したい研究では、解析の目的が異なります。前者では推定対象、交絡、モデルの仮定などが問題になり、後者では、新しい対象者に対する予測のズレが中心になります。
したがって、クロスバリデーションで良い成績が得られたという理由だけで、選ばれた変数を原因と解釈したり、選択後の通常のp値をそのまま因果的な根拠にしたりすることは推奨しません。Coxモデルの診断と変数選択でも、変数の選び方と、選んだ後の解釈を分けて扱っています。
分割ごとに予測を作る仕組み
ここで、クロスバリデーションは交差検証とも呼ばれ、英語ではcross-validation、以下ではCVと表記します。K-fold CVでは、データを重複しないK個の集合に分け、そのうちの一つを評価用に残し、残りでモデルを学習します。
また、5-fold CVなら、評価用に残す集合を順に変えながら、学習と予測を5回繰り返します。各対象者は一度だけ評価側へ入り、それ以外の反復では学習側に入りますが、自分の予測値を作るモデルの学習には参加しません。
ここで、対象者$i$のアウトカムを$y_i$、予測因子のベクトルを$\mathbf{x}_i$とし、第$k$評価集合の添字を$I_k$とします。また、$I_k$を除いたデータで学習した予測関数を$\widehat f^{(-k)}$、予測の誤りを数値化する損失関数を$L$、対象者数を$n$とすると、CVの平均損失は次のように書けます。
$$\widehat R_{\mathrm{CV}} =\frac{1}{n}\sum_{k=1}^{K}\sum_{i\in I_k} L\!\left(y_i,\widehat f^{(-k)}(\mathbf{x}_i)\right) \tag{1}$$
そのため、式(1)は、各対象者について、その人を学習から外したモデルの予測誤差を求め、全対象者で平均しています。同じ対象者の学習内予測と学習外予測を混ぜて平均する式ではありません。
分割の人数が異なる場合
第$k$集合の人数を$n_k$、その集合の平均損失を$\widehat R_k$とすれば、式(1)は次の加重平均と同じです。
$$\widehat R_{\mathrm{CV}} =\sum_{k=1}^{K}\frac{n_k}{n}\widehat R_k \tag{2}$$
この加重平均では、式(2)は、対象者を等しく重み付けするために、人数の多い評価集合ほど大きな重みを与えています。各集合の人数が同じ場合に限り、集合ごとの平均損失を単純平均した値と一致します。
施設単位の分割では、評価集合の人数が大きく異なることがあります。この場合、患者を等しく扱う平均と、施設を等しく扱う平均は異なる内容を表すため、どちらを報告するかも研究目的から決めます。
そして、二乗誤差を損失とする場合、式(1)から平均二乗誤差、MSEが得られ、その平方根が二乗平均平方根誤差、RMSEです。
$$\mathrm{RMSE}_{\mathrm{CV}} =\sqrt{\frac{1}{n}\sum_{k=1}^{K}\sum_{i\in I_k} \left\{y_i-\widehat f^{(-k)}(\mathbf{x}_i)\right\}^{2}} \tag{3}$$
ここで、式(3)は、学習外の予測誤差を二乗して平均し、その後に平方根を取っています。RMSEの単位はアウトカムと同じになり、誤差の大きさを元の測定尺度で読めます。
評価されるモデルと学習手順
CVを終えるとK個の学習済みモデルが残りますが、実際に使うのは全データで再度学習した別のモデルであることが多いため、CVの成績を最終モデルの性能そのものとみなすのは不適切です。
そこで、新しい対象者のデータが、学習データと同じ母集団から独立に得られる場合を考えます。学習データを$\mathcal D$、それを使ってモデルを作る手順を$\mathcal A$、新しい対象者を$(\mathbf X_{\mathrm{new}},Y_{\mathrm{new}})$とすると、固定した最終モデルの予測リスクは次の量です。
$$R(\mathcal D) =\mathbb E\!\left[ L\!\left(Y_{\mathrm{new}}, \mathcal A(\mathcal D)(\mathbf X_{\mathrm{new}})\right) \mid\mathcal D \right] \tag{4}$$
この式(4)は、手元の学習データと学習済みモデルを固定し、新しい対象者だけを入れ替えた場合の平均損失です。学習データ自体を集め直し、モデルも作り直した場合の平均性能とは区別しています。
通常のCVは、少し小さい学習集合で同じ手順を繰り返した成績を集めており、その手順に変数選択や調整値の選択が含まれるなら、それらも評価対象です。Batesらの理論解析も、手元の一つのモデルの誤差と、学習データを変えて作り直す手順の平均誤差を区別しています。
学習人数の違い
300人を均等に5分割すると、各モデルの学習人数は240人です。したがって、CVで直接予測に使われるのは240人で学習したモデルで、300人で学習した最終モデルではありません。
そのため、学習人数が増えるほど平均的な性能が改善する状況では、CVが最終モデルを悲観的に評価する場合がありますが、有限標本の変動や正則化の選択も重なるため、誤差の大小関係は条件によって変わります。
最小の誤差を選ぶ際の偏り
たとえば、LASSOの正則化係数を決める場合、複数の候補についてCVの誤差を計算し、その値が最も小さい候補を選ぶ方法があります。正則化の強さを$\lambda$、事前に定めた候補集合を$\Lambda$とすると、この選択は次の形になります。
$$\widehat\lambda =\operatorname*{arg\,min}_{\lambda\in\Lambda} \widehat R_{\mathrm{CV}}(\lambda) \tag{5}$$
式(5)は、候補ごとのCV誤差を比べ、その標本で最も小さかった値に対応する正則化係数を選んでいます。
ここで、各候補のCV誤差には偶然の変動が含まれるため、最小値を選ぶと、本当に性能が良い候補に加えて偶然小さく見えた候補も選ばれやすくなります。この関係は、固定した候補集合について次の不等式で表せます。
$$\mathbb E\!\left[\min_{\lambda\in\Lambda} \widehat R_{\mathrm{CV}}(\lambda)\right] \leq \min_{\lambda\in\Lambda} \mathbb E\!\left[\widehat R_{\mathrm{CV}}(\lambda)\right] \tag{6}$$
そのため、式(6)は、各データで最小値を選んでから平均すると、各候補の平均を求めてから最小値を選ぶ場合より小さくなりやすいことを表しています。候補同士の誤差が独立という仮定は、この不等式には不要です。
ただし、式(6)だけで、実際の研究の楽観的な偏りが何点になるかは分かりません。候補の数に加えて、候補同士の類似性、学習の不安定さ、評価標本の大きさが関わります。
CawleyとTalbot（2010）およびVarmaとSimon（2006）は、モデルを選ぶ評価基準にも過剰適合が起こることを示しているため、モデル選択に使った最小CV誤差は選択のための情報として読みます。
内側で選択し外側で評価する構成
300人のうち60人を外側の評価用に残し、残る240人だけで正則化の強さを選ぶ場面を考えます。この240人の中でさらに5分割し、内側では192人で学習、48人で評価する操作を繰り返します。
そして、内側の成績から$\lambda$を選んだ後に、240人全員を使って、その$\lambda$のモデルを学習し直します。最後に、ここまで一度もモデル選択へ使っていない外側の60人を予測します。
この外側の役割を入れ替え、300人全員について学習外の予測値を得る構成が、ネストしたCV、nested CVです。外側の第$k$学習集合だけから選ばれた設定を$\widehat\lambda_{-k}$とすると、平均損失は次の形です。
$$\widehat R_{\mathrm{nested}} =\frac{1}{n}\sum_{k=1}^{K}\sum_{i\in I_k} L\!\left(y_i, \widehat f_{\widehat\lambda_{-k}}^{(-k)}(\mathbf{x}_i)\right) \tag{7}$$
ここでも、式(7)は、対象者を外して係数を推定するだけでなく、その対象者を外したデータの中で正則化係数も選んでいます。内側で選んだ値は、外側の分割ごとに異なって構いません。
評価の外に置く処理
正則化係数だけを内側へ入れても、全データで選んだ変数を固定しているなら、変数選択の情報は外側へ漏れています。また、外側の結果を読んでモデルの種類や候補範囲を繰り返し変更すれば、その外側の成績も開発へ使ったことになります。
したがって、候補変数、前処理、モデルの種類、探索範囲、採用基準をどこまで外側の学習集合だけで再現したかが評価の意味を決めるため、これらを一続きの学習手順として記載します。
一方、モデルの構造も前処理も事前に固定され、データからチューニングする設定がない場合には、その固定手順を通常のCVで評価できます。ネストしたCVが必要になる理由は、分割を二重にすること自体ではなく、データに基づく選択と最終評価を分けることにあります。
乱数で生成した連続アウトカム
CVの操作を確認するために、独立した300人と、60個の候補予測因子を持つランダムに生成したデータセットを使います。各予測因子は互いに独立な標準正規分布から生成し、アウトカムとの関連を持つのは最初の3個だけとします。
ただし、この設定は多数の検査指標から連続的な測定値を予測する場面を表した生成モデルであり、実際の患者の検査値や特定の疾患、測定単位を再現したものではありません。
$$Y_i=10+2X_{i1}-1.5X_{i2}+X_{i3}+\varepsilon_i, \qquad \varepsilon_i\sim N(0,3^2) \tag{8}$$
式(8)は、最初の3個の予測因子から決まる条件付き平均に、標準偏差3の誤差を加えてアウトカムを生成しています。誤差は予測因子と独立で、対象者同士も独立とします。
そのうえで、候補変数を全て使ったLASSOなどの正則化を学習し、$\alpha=1$に固定して、$\lambda$だけを選びます。探索範囲は0.01から10までの対数等間隔の60候補とし、この範囲をCVの結果から変更する操作は行いません。
また、実行にはnestedcvパッケージを使い、内部の正則化回帰には依存パッケージのglmnetを用います。実行環境はR4.4.2、nestedcv0.8.0、glmnet4.1-8です。
set.seed(42)
library(nestedcv)
predictors <- matrix(rnorm(300 * 60), nrow = 300, ncol = 60)
colnames(predictors) <- paste0("marker", seq_len(60))
outcome <- 10 + 2 * predictors[, 1] - 1.5 * predictors[, 2] + predictors[, 3] + rnorm(300, sd = 3)
fold_id <- sample(rep(seq_len(5), length.out = 300))
outer_folds <- split(seq_len(300), f = fold_id)
nested <- nestcv.glmnet(outcome, x = predictors, family = "gaussian", alphaSet = 1, outer_folds = outer_folds, n_inner_folds = 5, min_1se = 0, cv.cores = 1,
lambda = exp(seq(log(10), to = log(0.01), length.out = 60)), type.measure = "mse")outer_foldsには、外側の評価集合に入る対象者の行番号を渡しています。n_inner_folds = 5は、その外側学習集合の内部を5分割する指定で、min_1se = 0は内側CVの平均二乗誤差が最小になる$\lambda$を採用する指定です。
また、この例には欠測値がなく、事前の単変量スクリーニングもありません。glmnetによる標準化は各学習時に行われますが、欠測処理や特徴量選択まで、このコードが自動的に安全な位置へ移してくれるわけではありません。
外側の評価集合の入れ替え
次に、最初に作成した300人の分割を図にすると、どの反復でも、60人が評価用として残っていることが分かります。図では読みやすくするために、対象者を割り当てられたfoldの順へ並べ替えていますが、実際の割り当ては乱数で決めています。
png("cv-fold-allocation.png", width = 1100, height = 800, res = 200, type = "cairo")
plot(rep(seq_len(300), times = 5), y = rep(seq_len(5), each = 300), pch = 15, cex = 0.5, col = ifelse(rep(sort(fold_id), times = 5) == rep(seq_len(5), each = 300), yes = "#0072B2", no = "#D8DEE4"), xlab = "Participants ordered by assigned fold", ylab = "Outer iteration", ylim = c(0.5, 6), yaxp = c(1, 5, 4))
legend("top", legend = c("Outer test (60)", "Outer training (240)"), col = c("#0072B2", "#D8DEE4"), pch = 15, bty = "n", cex = 0.8)
dev.off()
灰色の部分では、さらに内側の5-fold CVが動いています。図1は外側だけを表しているため、内側の48人の評価集合を、外側の60人の評価集合と取り違えないようにします。
また、CVの反復間では、同じ人が別のモデルの学習に使われます。この情報の共有があるため、外側で得た5個の誤差を、独立した研究を5回行った結果と同じようには扱えません。
評価の設計は分割数より先に決まる
ここまでで決めたのは、外側の60人へ何を渡さないかという一点です。分割数を5にしたことよりも、正則化係数の選択を内側へ入れたことの方が、得られる誤差の意味を大きく変えます。
次の記事では、この構成を動かし、内側の誤差曲線でどの正則化係数が選ばれたか、外側の予測で測った誤差が見かけの誤差とどれだけ違うかを、数値と図で確かめます。
参考文献
Stone M. Cross-validatory choice and assessment of statistical predictions. Journal of the Royal Statistical Society: Series B (Methodological). 1974;36(2):111-133.
Varma S, Simon R. Bias in error estimation when using cross-validation for model selection. BMC Bioinformatics. 2006;7(1):91.
Cawley GC, Talbot NLC. On over-fitting in model selection and subsequent selection bias in performance evaluation. Journal of Machine Learning Research. 2010;11:2079-2107.
Lewis MJ, Spiliopoulou A, Goldmann K, Pitzalis C, McKeigue P, Barnes MR. nestedcv: an R package for fast implementation of nested cross-validation with embedded feature selection designed for transcriptomics and high-dimensional data. Bioinformatics Advances. 2023;3(1):vbad048.
Friedman J, Hastie T, Tibshirani R. Regularization paths for generalized linear models via coordinate descent. Journal of Statistical Software. 2010;33(1):1-22.
