予測モデルの評価指標

予測モデル評価指標それぞれの 意味づけと指標間の関係 東京大学 大学院医学系研究科 生物統計学分野 篠崎 智大 [email protected] 日本計量生物学会 計量生物セミナー 2018年12月7日 18:25–19:25 @キャンパスプラザ京都

イベント発症の予測 イベント予測モデル 対象者の状態（共変量 Xi）からイベント発症 Yi を予測したい イベント発症 二値アウトカム   Yi = Yi(t) = 0, 1 定まった追跡時点 t までの発症有無がデータとして得られている 生存時間アウトカム Ti   Yi(t) = I(Ti < t) （追跡期間内の）任意の t でのイベントが予測対象だが、打ち切りの可能性 アウトカムの型によらず t-年イベント確率 pr(Yi (t) = 1|Xi) を求めたい 2

Yi(t) の予測確率 Pi(t) （リスクスコアとも） 二値アウトカム Y = Y(t) ロジスティックモデル pr(Y = 1|X = x) = {1 + exp(–xTβ )}–1 Pi = Pi(t) = {1 + exp(–XiT β� )}–1 生存時間アウトカム T Y(t) = I(Ti ≤ t) Coxモデル h(t|X = x) = h0(t)exp(xTβ ) t Pi(t) = 1 – exp{–∫ h� 0(u)du exp(XiT β� )} 0 3

Framingham Risk Score http://cvdrisk.nhlbi.nih.gov/ 4

• http://cvdrisk.nhlbi.nih.gov/

JALS-ECC AMI risk score JALS 急性心筋梗塞リスクスコアシート 5年以内に心筋梗塞を発症する確率を推定します。（non-HDLCモデル） non-HDLコレステロールは自動的に計算されます。 ピンク色部分はプルダウンで選択して下さい。 黄色部分は値を入力して下さい。 項目 値 プルダウンで選んでください 性 年齢（歳） 総コレステロール（TC, mg/dL） HDLコレステロール（HDLC, mg/dL） 血圧 収縮期（上の血圧, mmHg） 拡張期（下の血圧, mmHg） プルダウンで選んでください 糖尿病 プルダウンで選んでください 喫煙 スコア 計 5年以内の心筋梗塞発症確率 http://jals.gr.jp/result/datas/2010_74_1346-56AMI03.xls 5

• http://jals.gr.jp/result/datas/2010_74_1346-56AMI03.xls

JJ Risk Engine http://www.med.niigata-u.ac.jp/emh/jjre.html 6

• http://www.med.niigata-u.ac.jp/emh/jjre.html

作業モデル（working model）としての予測モデル Pi を得るための予測モデル pr(Y = 1|X = x) = {1 + exp(–xTβ )}–1 h(t|X = x) = h0(t)exp(xTβ ) 予測モデルは「誤特定」されていても良い モデルに含まれる変数 X 選ばれた変数 X で条件付けた回帰に対するモデル形状   E(Y|X = x) = pr(Y = 1|X = x) h(t|X = x) 予測モデル評価は特定のモデルを仮定しない指標に基づくのが好ましい AIC（赤池情報量基準）によるモデル選択と同様の思想 ここではモデルそのものではなく、既に与えられている Pi の評価を考える 7

予測モデル（リスクスコア）の評価 アウトカム発症 Yi = Yi(t) 予測確率 Pi = Pi(t) 何らかの「近さ」を測る指標がほしい 連続アウトカムの場合：R2 ３つの解釈 8

正規線形モデルでの R2 1. Mittlebock & Shemper, 1996 � i）で説明された Yi の平方和の割合（explained variation） 予測値 Pi（= Y � i)2 /∑i(Yi – Y) �2 1 – ∑i(Yi – Y 予測誤差の指標 2. � i） の相関係数の2乗 結果変数 Yi と予測値 Pi（= Y �Y)2/√{var(Y) � � � � Y)} cov(Y, var( 判別、較正の指標 3. 無情報モデルとの尤度比の（データあたりの）幾何平均 1 – {L(0)/L(β� )}2/n 作業モデルの正しい特定を要するので扱わない 9

今日のトピック（1）：指標どうしのつながり 予測誤差（prediction error） Brierスコア 判別（discrimination） Se(v), Sp(v), ROC曲線, C 較正（calibration） Calibration-in-the-large, calibration slope 較正プロット, Hosmer-Lemeshow統計量 再分類（reclassification） 再分類表, NRI(r1,...,rk), cNRI（NRI>0）, IDI 意思決定分析（decision analysis） NB(v) 意思決定曲線 10

今日のトピック（2）：生存時間アウトカムへの拡張 打ち切りの考慮 Tobs = min(T, U)  U：打ち切りまでの時間 Y(t) = I(T < t) が分からない対象者がいる 生存時間特有の指標 Cτ （overall C）  T そのものを判別 11

データ (Yi, Pi), i = 1,..., n Yi ： 時点 t でのイベント有無（1: あり、0：なし） Pi ： 時点 t でのイベント予測確率（P ∈ [0,1]） � イベント割合 Y � = ∑iYi /n Y � 予測値の平均 P � = ∑iPi /n P Yi 1 0 0 1 Pi 12

予測誤差 prediction error SAS 対象者 i での二乗誤差 squared error (Yi – Pi)2 予測誤差：Yi と Pi との「距離」を測る指標の総称  例）絶対誤差、カルバックライブラー（エントロピー）、... （期待）Brierスコア BS = E{(Y – P)2} = E{Y(1 – P)2 + (1 – Y)P2} データから推定 2 ∑ (Y – P ) i � = i i BS n 13

Brierスコアはそのままでは解釈しにくい 数式上の上下限 Pi = Yi （for all i）と完全に予測できれば下限 0（best） Pi = 1 – Yi （for all i）のとき上限 1（worst）  各対象者に予測確率 Pi = 1 または Pi = 0 を割り当てるような予測モデル より現実的な「無情報」モデル 各対象者におなじ予測確率 Pi = P を与えるモデル 予測にまったく寄与しない 2 ∑ (Y – P) � = i i � – Y) � � 2 + Y(1 BS = (P – Y) n P と Y との系統的なずれ （calibration-in-the-large） Y の分散 14

SAS Scaled Brier score 無情報モデルのBrierスコア 2 ∑ (Y – P) � = i i � – Y) � � 2 +Y(1 = (P – Y) BS n � = 0.25 � によらず） BS Pi = 0.5 を与えたとき、（Y � = Y(1 � を与えると、BS � – Y) � = P(1 – P) （ Y � = P = 0.5 のとき、最大値 0.25） Pi = Y スケール化（⇒ 1：best, 0: worst） �s = 1 – BS � BS � – P) � P(1 � 予測モデル Pi に対する「無情報」モデル： Pi = P � – Y）が無い状況での「最悪値」 � 無情報モデルに系統的なずれ（P Steyerberg, et al. 2010 15

スケール化したBrierスコア vs. 平方和にもとづく R2 平方和の分解にもとづく R2 （= Efron の疑似 R2） � 2 = ∑i(Yi – Y � i)2 + ∑i(Y � i –Y) �2 ∑i(Yi – Y) � i = Pi かつ P �=Y � おなじデータから予測 Pi を得ると Y Brierスコア � 2 ∑ (Y –P ) BS REfron2 = 1 – i i �i 2 = 1 – � � ∑i(Yi – Y) Y(1 – Y) �s = 1 – BS � BS � – P) � P(1 平方和の分解を モデル評価の文脈に置き換えた指標と見なせる 16

Brierの指摘 多値アウトカムに対するBrierスコアの原型 観測 Ei1,..., Eir 予報 fi1,..., fir ∑i∑j=1,...,r(Eij – fij)2 n Brierスコアを高めるための予測のあり方 おなじ予測値を与える状況では、予測値は「較正」されるべき  予測降水確率がおなじ P の日を集めると、確率 P で雨天となるように 確率 1 か 0 を割り当てる予測より、不確実な場合は確率 Pi で個々を予測  0 < Pi < 1 確率的な予測においては、「判別」が大事   なるべく Pi の分布が極端になるように（できるだけ 0 か 1 に近い方が良い） Eir と fir との「相関」が高くなる Brier, 1950 17

SAS 判別 discrimination イベント発症者（Y = 1）と非発症者（Y = 0）を正しく弁別する特性 Harrell et al. 1996 P をカットポイント v で上下に分ける 感度 sensitivity Se(v) = pr(P > v|Y = 1) 特異度 specificity Sp(v) = pr(P ≤ v|Y = 0) データから推定 ∑ I(P > v)Yi � Se(v) = i i ∑iYi ∑ I(P ≤ v)(1 – Yi) � Sp(v) = i i ∑i(1 – Yi) Yi = 0 Yi = 1 0 v 1 Pi Yi 1 0 0 v 1 Pi 18

SAS Se(v) と Sp(v) の統合：ROC曲線 Se(v) = pr(P > v|Y = 1) 1 – Sp(v) = pr(P > v|Y = 0) P に対する「生存関数」 単調減少 Se(v) 1 0 ROC（receiver operating characteristics）曲線 横軸： 1 – Sp(v) 縦軸： Se(v) 1 0 P 1 – Sp(v) 0 1 v Se(v) 1 曲線下面積 AUC（area under the curve） ROC曲線の要約 0 0 1 1 – Sp(v) 19

ROC曲線AUC “concordance” index に一致 C = pr(Pi > Pj |Yi = 1, Yj = 0) Y = 1からランダムサンプル Y = 0からランダムサンプル 0 1 Pi Hanley & McNeil, 1982 20

SAS C “concordance” index C = pr(Pi > Pj |Yi = 1, Yj = 0) データから推定 Yi = 1, Yj = 0 から成る全てのペア (i, j) で数え上げ（n1×n0組） n1 = ∑iYi n0 = ∑i(1 – Yi) ∑i,j I(Pi > Pj)I(Yi = 1, Yj = 0) �= C ∑i,j I(Yi = 1, Yj = 0) Mann-Whitney統計量に等しい  順序統計量 21

R2 再び：Y と P との相関としての解釈 Pearsonの積率相関 r P が Y の線形変換なら「平方和にもとづく R2 」に一致 スケール化したBrierスコアに類似 Spearmanの順位相関 ρ P を順位変換した積率相関 Kendallの順位相関 τ 全てのペアに対して「Y と P の大小一致割合 pc」 – 「不一致の割合1 – pc」  τ = 2pc – 1 ⇒ 順序一致の割合 pc = (τ + 1)/2 「Y が異なるペア」に限定したものが、Mann-Whitney統計量  Yi = 1, Yj = 0 から成るペア (i, j) 22

P から Y への「回帰」としての判別 相関というよりも... Y = 1 vs. Y = 0 での P の比較 C = pr(Pi > Pj |Yi = 1, Yj = 0) C – (1 – C) = 2C –1 を Somer の D という Kendallの順位相関τ の「回帰」バージョンとみることができる Mittlebock & Schemper, 1996 E(P|Y = 1) – E(P|Y = 0) Yates の判別係数（discrimination slope） P を Y に回帰 = r×√{var(P)/var(Y)} 23

較正 calibration 判別の評価指標 Se(v), Sp(v), ROC曲線, C いずれも P のスケールに対して不変、順序のみを評価している 裏を返せば、P の値自体と Y との「近さ」を直接評価できていない 較正 予測リスクと実際のリスクの近さを測る 1. 2. 予測リスクと実際のリスクとの系統的なずれ： E(P) – pr(Y = 1) 予測リスク Pi と実際のリスク pr(Y = 1|Pi) のずれ 24

系統的な較正：calibration-in-the-large 予測リスク Pi の平均値 イベント Y = 1 の割合 SAS の差：pr(Y = 1) – E(P) 比  pr(Y = 1)/E(P) オッズ比  pr(Y = 1){1 – E(P)}/{pr(Y = 0)E(P)} calibration-in-the-largeの検定 作業モデル log{pr(Y = 1|P)/pr(Y = 0|P)} = α0 + 1∙log P/(1 – P) α0 = 0（オッズ比 = 1）の検定が calibration-in-the-large = 0 の検定  ただし、exp(� α0) はオッズ比の推定値とはならない（オッズ比のnon-collapsibilityより） Miller et al. 1993 25

個々の P の予測：calibration slope SAS 個々の予測が個々のリスクに近いか Pi vs. pr(Y = 1|Pi) Pi = pr(Y = 1|Pi) の検定 作業モデル α1 = 1 の検定  log{pr(Y = 1|P)/pr(Y = 0|P)} = α1 log P/(1 – P) α1：refinement または calibration slope という Cox, 1958 Miller et al., 1993 注意 calibration-in-the-large と calibration slope は一意な指標ではない 26

calibration-in-the-large + calibration slope SAS 2つの較正指標をまとめて検定 log{pr(Y = 1|P)/pr(Y = 0|P)} = α0 + α1 log P/(1 – P)  (α0, α1) = (0, 1) の自由度2の検定 Cox, 1958 検定に関する注意 最尤法で Pi = {1 + exp(–XiT β� )}–1 を得たとき 1. 同一データでは必ず (� α0, α �1) = (0, 1) 2. 内部検証法では検定のサイズが保たれる  (a) 外部検証か、(b) 最尤法を用いずに作られたリスクスコア評価を前提 ロジスティックモデルの正しい特定は仮定しない  ただし、logit pr(Y = 1|P) と logit P に線形関係がないと検出力低下 27

SAS Hosmer-Lemeshow統計量 予測値 P に応じて集団を K 分割（通常は K = 10） 各集団 k のサイズ nk �k 各集団 k でのイベント数 Ok = nkY �k 各集団 k での予測値平均 P K K k=1 k=1 2 2 � � � P (O –n ) ( Y – P ) χ2HL = � �k k �k = � � k � k Pk(1 – Pk)/nk nkPk(1 – Pk) 自由度 K – 2 のカイ二乗検定が Pi = pr(Y = 1|Pi) の検定  � k = ∑kOk の制約より ∑knkP 大きい p-値の方がよく較正された予測 28

較正プロット calibration plot Hosmer-Lemeshow統計量を図示 �k 横軸に予測値平均 P �k 縦軸にイベント割合 Y  0.12 割合の信頼区間をつけることも k = 1,..., K をプロット Y と P との「一致度」 �k = P � k の直線 Y 「相関」を表す判別とは補完的に 0.08 �k Y 0.04 0.00 0.00 0.04 0.08 �k P 0.12 29

異なる予測どうしの比較 Pnew vs. Pold 1. Pnew と Pold それぞれで評価指標を計算、比較 例： Cnew – Cold とその検定・信頼区間 「予測力の平均的指標」を比較 2. Pnew と Pold による予測の変化の指標 「個々の予測力比較」の平均的指標 30

再分類 reclassification 予測値 Pi によって個人 i をリスク分類 カットオフ値 r1,..., rJ   例：J = 1  P < r1：経過観察群、 r1 ≤ P：要治療群 例：J = 2  P < r1：低リスク、 r1 ≤ P < r2：中リスク、 r2 ≤ P：高リスク Pi ⇒ Ri  R ∈ {1,..., J + 1} 異なる予測によるリスク分類 Pnew,i, Pold,i ⇒ Rnew,i, Rold,i Rnew と Rold でクロス集計 31

再分類表 reclassification table Pold Pnew 0%--<5% Rnew = 1 5%--<10% Rnew = 2 10%-Rnew = 3 0%--<5% Rold = 1 542 84 1 5%--<10% Rold = 2 32 320 51 10%-- Rold = 3 2 14 211 対角線上はリスク分類が変わらなかった対象者 右上は新スコアで予測リスク上昇 イベントあり（Y = 1）で増えていて欲しい 左下は新スコアで予測リスク下降 イベントなし（Y = 0）で増えていて欲しい 32

SAS イベントごとの再分類表 Y = 1（n = 219） Y = 0（n = 1038） Rnew = 1 Rnew = 2 Rnew = 3 Rold = 1 510 20 0 40 Rold = 2 23 292 11 42 Rold = 3 2 11 169 Rnew = 1 Rnew = 2 Rnew = 3 Rold = 1 32 64 1 Rold = 2 9 28 Rold = 3 0 3 好ましい再分類 pr(Rnew > Rold|Y = 1) 好ましくない再分類 pr(Rnew < Rold|Y = 1) 好ましくない再分類 pr(Rnew > Rold|Y = 0) 好ましい再分類 pr(Rnew < Rold|Y = 0) 33

SAS NRI（net reclassification improvement） 「ネット」再分類指標 予測 Pnew vs. Pold カットオフ値 r1,..., rJ NRI(r1,..., rJ) = pr(Rnew > Rold|Y = 1) – pr(Rnew < Rold|Y = 1) + pr(Rnew < Rold|Y = 0) – pr(Rnew > Rold|Y = 0) 前頁の例  NRI(5%, 10%) = 105/219 – 12/219 + 36/1038 – 31/1038 = 43.0% Pencina et al., 2008 34

SAS cNRI “continuous” NRI category-free NRI, NRI>0とも 「すべての」カットオフ値で再分類を（= P そのものの大小） cNRI = pr(Pnew > Pold|Y = 1) – pr(Pnew < Pold|Y = 1) + pr(Pnew < Pold|Y = 0) – pr(Pnew > Pold|Y = 0) 書き直すと Vi = I(Pnew,i > Pold,i) – I(Pnew,i < Pold,i) cNRI = E(V|Y = 1) – E(V|Y = 0) 35

IDI（integrated discrimination improvement） SAS cNRI = E(V|Y = 1) – E(V|Y = 0) V は3値しかとらない Pnew と Pold の乖離度合いを反映していない 重み付けcNRI wi = |Pnew,i – Pold,i| IDI = E(wV|Y = 1) – E(wV|Y = 0) = {E(Pnew|Y = 1) – E(Pnew|Y = 0)} – {E(Pold|Y = 1) – E(Pold|Y = 0)} 36

再分類から再び判別へ IDI = {E(Pnew|Y = 1) – E(Pnew|Y = 0)} – {E(Pold|Y = 1) – E(Pold|Y = 0)} 新旧モデルのYatesの判別係数（discrimination slope; p. 23）の差 “integrated” discrimination   E(P|Y = 1) = ∫1 Se(v)dv 0 1 E(P|Y = 0) = ∫ {1 – Sp(v)}dv Se(v) と 1 – Sp(v) はP に対する「生存関数」 0 再分類指標とも判別指標とも見ることができる Pencina et al., 2008 37

カットオフ値による「意思決定」 実用上は P の値を見たあとで医療行為をするかどうかを判断 1つのカットオフ値 v による二元的な意思決定   判別 再分類 意思決定に伴う「コスト」 P > v の場合  イベントあり（Y = 1）を治療する便益 + イベントなし（Y = 0）を治療するコスト P < v の場合  イベントあり（Y = 1）を治療しないコスト + イベントなし（Y = 0）を治療しない便益 38

カットオフ値とは 治療するかどうかの判断が均衡する予測値 P = v の場合   イベントあり（Y = 1）を治療する効用 O1T + イベントなし（Y = 0）を治療する効用 O0T イベントあり（Y = 1）を治療しない効用 O1U + イベントなし（Y = 0）を治療しない効用 O0U その人のイベント確率を Pi = v と予測しているので vO1T + (1 – v)O0T = vO1U + (1 – v)O0U (O1T – O1U)/(O0T – O0U) = v/(1 – v) Y=1を 治療した便益 Y=0を 治療したコスト 39

ネットベネフィット カットオフ値 v での判断（P > v）に伴う、治療を行う便益 vs. コストの差 便益： Y = 1 で治療できた集団  pr(Y = 1, P > v) コスト： Y = 0 で治療してしまった集団  pr(Y = 0, P > v) この比が v/(1 – v) なので v NB(v) = pr(Y = 1, P > v) – pr(Y = 0, P > v) 1–v Vickers & Elkin, 2006 Baker et al., 2014 40

意思決定曲線 decision curve NB(v) を v に対してプロット NB(v) = pr(Y = 1, P > v) – v pr(Y = 0, P > v) 1–v 0.2 2つの予測モデル Pold と Pnew 0.1 NB(v) 0.0 Vickers & Elkin, 2006 Baker et al., 2014 -0.1 0.0 0.1 0.2 v 0.3 0.4 41

SAS 意思決定曲線 decision curve NB(v) を v に対してプロット NB(v) = pr(Y = 1, P > v) – v pr(Y = 0, P > v) 1–v 「誰も治療しない方針」（“treat none”） 0.2 Pnone,i = 0 v NBnone(v) = 0 – 0=0 1–v 「全員治療する方針」（“treat all”） 0.1 0.0 Pall,i = 1 v NBall(v) = pr(Y = 1) – pr(Y = 0) 1–v -0.1 0.0 0.1 0.2 0.3 0.4 42

（再）今日のトピック（1）：指標どうしのつながり 予測誤差（prediction error） Brierスコア 判別（discrimination） Se(v), Sp(v), ROC曲線, C 較正（calibration） Calibration-in-the-large, calibration slope 較正プロット, Hosmer-Lemeshow統計量 再分類（reclassification） 再分類表, NRI(r1,...,rk), cNRI（NRI>0）, IDI 意思決定分析（decision analysis） NB(v) 意思決定曲線 43

（再）今日のトピック（2）：生存時間アウトカムへの拡張 打ち切りの考慮 Tobs = min(T, U)  U：打ち切りまでの時間 Y(t) = I(T ≤ t) が分からない対象者がいる 生存時間特有の指標 Cτ （overall C）  T そのものを判別 44

データ (Ti*, Di, Pi(t)), i = 1,..., n; t > 0 Ti* = Tobs,i = min(Ti, Ui)   Ti ：イベントまでの時間 Ui ：打ち切りまでの時間 Di = I(Ti* = Ti) Pi(t)：Yi(t) = I(Ti ≤ t) の予測確率 t  例：Pi(t) = 1 – exp{–∫ h� 0(u)du exp(XiT β� )} 0 Yi(t) が観察されない人がいる 45

Yi(t) が観察されないのは タイプ1： Ti* ≤ t で Di = 1 の対象者 Yi(t) = 1 タイプ2： Ti* > t の対象者 Yi(t) = 0 タイプ3： Ti* ≤ t で Di = 0 の対象者 Yi(t) = ?? 46

時点 t での予測誤差 期待Brier score（p. 13） BS(t) = E[Y(t){1 – P(t)}2 + {1 – Y(t)}P(t)2] タイプ1 タイプ2 タイプ1 & 2 は Y(t) – P(t) の計算に使える タイプ3 は使えない ⇒ 「各時点での打ち切り確率」の推定で利用 IPCW（inverse probability censoring weighting） タイプ1：イベント発症時点 Ti* まで打ち切られない確率 タイプ2：イベントなし確認時点 t まで打ち切られない確率 の逆数 Graf et al., 1999 47

打ち切られない確率 SAS 打ち切り分布 G(s) = pr(U > s) � Kaplan-Meier推定：G(s) (Ti*, Di) を利用 Di = 1 を「打ち切り」として扱う 48

SAS IPCW BS(t) 推定量 タイプ1の重み タイプ1 タイプ2の重み タイプ2 � i∗)–1I(Ti∗ ≤ t)Di {1 – Pi(t)}2 +G(t) � –1I(Ti∗ > t)Pi(t)2] ∑i[G(T � = BS(t) � i∗)–1I(Ti∗ ≤ t)Di + G(t) � –1I(Ti∗ > t)} ∑i {G(T スケール化 P(t) = ∑i Pi(t)/n � BS s(t) = 1 – � BS(t) P(t){1 – P(t)} Graf et al., 1999 49

時点 t での判別 Se(v, t) = pr{P(t) > v | Y(t) = 1} {1 – S(t|P(t) > v)} pr{P(t) > v} = 1 – S(t) Sp(v, t) = pr{P(t) ≤ v | Y(t) = 0} ROC(t) S(t) = pr(T > t) S(t|∙) = pr(T > t|∙) S(t|P(t) ≤ v) pr{P(t) ≤ v} = S(t) 1 – Sp(v, t) vs. Se(v, t) pr(∙) は割合から推定、 S(t) と S(t|∙) は K-M推定 Heagerty et al., 2000 50

SAS

いくつかのノンパラメトリック法
1.

条件付き Kaplan-Meier 法（前頁）
簡便だが、Se(v, t) と Sp(v, t) が v に関して単調にならない

2.

最近傍法 （nearest neighbor method）
単調性を保つよう S(t|∙) を推定
Heagerty et al., 2000; French et al., 2016

3.

再帰法 （recursive method）
スムージング手法を利用
Chambless & Diao, 2006

4.

IPCW法
タイプ1から pr{P(t) > v | Y(t) = 1}、タイプ2から pr{P(t) ≤ v | Y(t) = 0}
Uno et al., 2007
51

IPCW Se(v, t) & Sp(v, t) 推定量 タイプ1 は Y(t) = 1 � i∗)–1Di I{Ti∗ ≤ t, Pi(t) >v} ∑iG(T � = Se(v,t) � i∗)–1Di I(Ti∗ ≤ t) ∑iG(T タイプ1の重み タイプ1 タイプ2 は Y(t) = 0 � –1I{Ti∗ > t, Pi(t) ≤ v} ∑iG(t) � = Sp(v,t) � –1I(Ti∗ > t) ∑iG(t) タイプ2の重み タイプ2 （キャンセル） 52

C(t) 時点 t でのイベント発症者と非発症者を弁別 C(t) = pr{Pi(t) > Pj(t) |Yi(t) = 1, Yj(t) = 0} Se(v, t) と Sp(v, t) のどの推定方法からでも ROC(t) を描ける AUC(t) として C(t) を推定 （補）cumulative/dynamic ROC(t) という t までの全イベント（cumulative case）から感度 Se(v, t) を定義 t ごとに変化するリスクセット（dynamic control）から特異度 Sp(v, t) を定義 Heagerty & Zheng, 2005 53

生存時間アウトカムの判別 discrimination
時点 t でのイベント発症ありなしを弁別する特性
Y(t) = 1 vs. Y(t) = 0

C(t) = pr{Pi(t) > Pj(t) |Yi(t) = 1, Yj(t) = 0}
イベント発症までの時間の長短を弁別する特性
T そのもの

τ ：制限時間

Cτ = pr{Pi(t) > Pj(t) |Ti < Tj, Ti < τ}

Overall C という
Harrell et al., 1996; Pencina & D’Agostino, 2004
54

2つの C C(t) = pr{Pi > Pj|Yi(t) = 1, Yj(t) = 0} Yi(t) = 1, Yj(t) = 0 から成るペア (i, j) から選ぶ  n1×n0 ペア Mann-Whitney の U-統計量 Cτ = pr(Pi > Pj|Ti < Tj, Ti < τ) （τ = ∞ にとると）すべてのペア  n(n – 1)/2 ペア Kendall の τ-統計量 55

SAS

Cτ の推定量
Harrellの推定量

comparable のうち T と P の大小がそろっている（concordant）ペアの数

∑i,j Di I(Ti∗ <Tj∗) I{Pi(t) > Pj(t)}
�τ =
C
∑i,j Di I(Ti∗ <Tj∗)
比較可能な（comparable）ペアの数

Harrell et al., 1996

comparableである必要十分条件
ペア内で観察時間 T * の短い方のイベントが観察されていること（D = 1）

打ち切りの存在下では U の分布 G(∙) を含む指標に収束
Uno et al., 2011; Pencina et al., 2012
56

SAS

Cτ のIPCW推定量
「ペア」を観察確率（打ち切られない確率）の逆数で重み付け

ペア内で T * の短い方のイベントが観察されていれば（= comparable）利用可
ペア (i, j) の短い方の Ti* までどちらも打ち切られない確率 = G(Ti*)2
Kaplan-Meier法で推定
concordant

� i∗)–2Di I(Ti∗ <Tj∗) I{Pi(t) > Pj(t)}
∑i,jG(T
�τ,IPCW =
C
� i∗)–2Di I(Ti∗ <Tj∗)
∑i,jG(T
comparable

Uno et al., 2011

57

Cτ の回帰モデルを使った推定量 incident/dynamic ROC(t)（p. 53）の AUC(t) を重み付け積分 incident：t でのイベント（incident case）のみで感度 Se(v, t) を定義 dynamic： t ごとに変化するリスクセット（dynamic control）で特異度 Sp(v, t) を定義 Heagerty & Zheng, 2005 Coxモデルの線形予測子 XiT β� の変換 非常に簡便 Gönen & Heller, 2006 いずれもCoxモデルの正しい特定が必要 「作業モデル」を評価することができない 58

較正：calibration-in-the-large, calibration slope 時点 t での較正は難しい pr{Y(t) = 1} vs. E{P(t)} pr{Y(t) = 1|P(t)} vs. P(t) 代わりに、観察イベントインディケータ Di = Yi(Ti*) vs. Pi(Ti*) Poisson回帰（作業モデル） log pr{Di = 1| Ti*, Pi(Ti*)} = α0 + log Pi(Ti*) log pr{Di = 1| Ti*, Pi(Ti*)} = α1 log Pi(Ti*) α0 = 0 の検定（calibration-in-the-large）、α1 = 1 の検定（calibration slope）  log Pi(Ti*)：マルチンゲール残差 Di – exp{–∫ Ti ∗ 0 h� 0(u)du exp(XiT β� )} から得られる Crowson et al., 2016 （注：α1 = 1 の検定はモデル構築と同一サンプルでも検定のサイズを保てない） 59

Hosmer-Lemeshow統計量（再） イベント有無が分かっている場合 P で分割したグループ k でのKaplan-Meier推定値 k = 1,…, K （通常10） K �k – P � k)2 ( Y χ2HL = � � � k)/nk Pk(1 – P k=1 60

SAS D’Agostino & Nam統計量 イベント割合をKaplan-Meier推定 � k(t) � k(t) = 1 – S K K � k(t) – P � k(t)}2 { K χ2DN(t) = � � � k(t)}/nk Pk(t){1 – P 自由度 K – 1 のカイ二乗検定 打ち切りが多いとき k=1 D’Agostino & Nam, 2005 K � k(t) – P � k(t)}2 { K χ2DNG(t) = � � k(t)} � K var{ k=1 Greenwood推定量 Demler et al., 2015 61

較正プロット D’Agostino & Nam統計量を図示 � k(t) 横軸に予測値平均 P � k(t) 縦軸に K-M推定値 K  � k(t)} で信頼区間をつけることも � K Greenwood公式 var{ k = 1,..., K をプロット 0.12 0.08 � k(t) K 0.04 0.00 0.00 0.04 0.08 0.12 � k(t) P 62

再分類表の作り方 Pnew(t) vs. Pold(t) カットオフ値 r1,..., rJ Pold(t) Pnew(t) <r1 Rnew(t) = 1 r1--<r2 Rnew(t) = 2 r2-Rnew(t) = 3 <r1 Rold(t) = 1 n11 n12 n13 r1--<r2 Rold(t) = 2 n21 n22 n23 r2-- Rold(t) = 3 n31 n32 n33 � ij(t) 各セル (Rold(t), Rnew(t)) = (i, j) ごとにKaplan-Meier推定 K  � ij(t) � ij(t) = 1 – S K 63

イベントごとの再分類表と NRI(t; r1,…, rJ) Y(t) ごとの「期待」再分類表 Y(t) = 0 Y(t) = 1 Rnew = 1 Rold = 1 Rold = 2 Rold = 3 Rnew = 2 Rnew = 1 Rnew = 3 � 11(t) n12K � 12(t) n13K � 13(t) n11K � 21(t) n22K � 22(t) n23K � 23(t) n21K � 31(t) n32K � 32(t) n33K � 33(t) n31K Rold = 1 Rold = 2 Rold = 3 Rnew = 2 Rnew = 3 � 11(t) n12S � 12(t) n13S � 13(t) n11S � 21(t) n22S � 22(t) n23S � 23(t) n21S � 31(t) n32S � 32(t) n33S � 33(t) n31S NRI(t; r1,…, rJ) = pr{Rnew(t) > Rold(t)|Y(t) = 1} – pr{Rnew(t) < Rold(t)|Y(t) = 1} + pr{Rnew(t) < Rold(t)|Y(t) = 0} – pr{Rnew(t) > Rold(t)|Y(t) = 0} Steyerberg & Pencina, 2010 64

SAS

NRI(t)
再分類表を作らなくても
up(t):
“Rnew(t) > Rold(t)” or “Pnew(t) > Pold(t)”
down(t): “Rnew(t) < Rold(t)” or “Pnew(t) < Pold(t)”
NRI(t; r1,…, rJ)

cNRI(t)

[1 – S{t|up(t)}] pr{up(t)} – [1 – S{t|down(t)}] pr{down(t)}
NRI(t) =
1 – S(t)
S{t|down(t)} pr{down t } – S{t|up(t)} pr{up t }
+
S(t)
Pencina et al., 2011
65

SAS IDI(t) Pnew(t) と Pold(t) との判別係数の差 IDI(t) = [E{Pnew(t)|Y(t) = 1} – E{Pnew(t)|Y(t) = 0}] – [E{Pold(t)|Y(t) = 1} – E{Pold(t)|Y(t) = 0}] 注 var{Pnew(t)} var{Pold(t)} − = S(t){1 – S(t)} S(t){1 – S(t)} Chambless et al., 2011 ただし、E{Pnew(t)} = 1 – S(t)、 S{t|P(t)} = 1 – P(t) などを仮定（同一標本でのみ成立） 外部検証ではモデルを仮定しない手法としてIPCW法が適用可 66

NB(v, t) SAS v NB(v, t) = pr{Y(t) = 1, P(t) > v} – pr{Y(t) = 0, P(t) > v} 1–v = Se(v, t) {1 – S(t)} = {1 – Sp(v, t)} S(t) S(t) Kaplan-Meier推定 Se(v, t) と Sp(v, t) 4つのノンパラメトリック推定法（p. 51） 条件付きKaplan-Meier推定 = Vickers et al.（2008）の方法 67

後半の内容 生存時間（打ち切り）データ IPCWによる重み付け Kaplan-Meier法による生存関数推定 生存時間特有の指標も Overall C 68

今日お話しできなかったトピック 各指標の臨床的解釈および注意点 特に再分類指標 モデルの過適合（overfit）とオプティミズム 内部検証 外部検証 モデルのアップデート、再較正 あしたのお話 井上先生 再分類の指標 NRI/IDI 宇野先生 “Moving Beyond Association” 69

まとめ イベント予測モデル 二値アウトカム 生存時間アウトカム イベント予測の評価 イベント予測の指標には R2 のような「統一指標」がない R2 のいくつかの解釈への類比と、それを超えて指標が発展してきたとも見える  予測誤差、判別、較正、再分類 評価はなるべくモデルの正しい特定を前提としない指標と手法を 謝辞 横田勲先生（北海道大学）と大庭幸治先生（東京大学）には一部内容の整理とともに文献情報を、 上妻佳代子さん（東京大学）には付録の統計パッケージ情報を提供いただきました 70

文献 Baker, S.G., Schuit, E., Steyerberg, E.W., Pencina, M.J., Vickers, A., Moons, K.G., Mol, B.W. & Lindeman, K.S. (2014). How to interpret a small increase in AUC with an additional risk prediction marker: decision analysis comes through. Statistics in Medicine 33, 3946–3959. Brier, G.W. (1950). Verification of forecasts expressed in terms of probability. Monthly Weather Review 78, 1–3. Chambless, L.E. & Diao, G. (2006). Estimation of time-dependent area under the ROC curve for long-term risk prediction. Statistics in Medicine 25, 3474–3486. Chambless, L.E., Cummiskey, C.P. & Cui, G. (2011). Several methods to assess improvement in risk prediction models: extension to survival analysis. Statistics in Medicine 30, 22–38. Cox, D.R. (1958). Two further applications of a model for binary regression. Biometrika 45, 562–565. Crowson, C.S., Atkinson, E.J. & Therneau, T.M. (2016). Assessing calibration of prognostic risk scores. Statistical Methods in Medical Research 25, 1692–1706. D’Agostino, R.B. & Nam, B.H. (2005). Evaluation of the performance of survival analysis models: discrimination and calibration measures. In: Balakrishan, N. and Rao, C.R. (editors). Handbook of Statistics, 23. Philadelphia, PA: Elsevier, pp. 1–25. Demler, O.V., Paynter, N.P. & Cook, N.R. (2015). Tests of calibration and goodness-of-fit in the survival setting. Statistics in Medicine 34, 1659–1680. French, B., Saha-Chaudhuri, P., Ky, B., Cappola, T.P. & Heagerty, P.J. (2016). Development and evaluation of multi-marker risk scores for clinical prognosis. Statistical Methods in Medical Research 25, 255–271. 71

文献 Gönen, M. & Heller, G. (2005). Concordance probability and discriminatory power in proportional hazards regression. Biometrika 92, 965–970. Graf, E., Schmoor, C., Sauerbrei, W. & Schumacher, M. (1999). Assessment and comparison of prognostic classification schemes for survival data. Statistics in Mediciine;18:2529–2545. Hanley, J.A. & McNeil, B.J. (1982). The meaning and use of the area under a receiver operating characteristic (ROC) curve. Radiology 143, 29–36. Harrell, F.E., Lee, K.L. & Mark, D.B. (1996). Tutorial in Biostatistics. Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors. Statistics in Medicine 15, 361–387. Heagerty, P.J., Lumley, T. & Pepe, M.S. (2000). Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics 56, 337–344. Heagerty, P.J. & Zheng, Y. (2005). Survival model predictive accuracy and ROC curves. Biometrics 61, 92–105. Miller, M.E., Langefeld, C.D., Tierney, W.M., Hui, S.L. & McDonald, C.J. (1993). Validation of probabilistic predictions. Medical Decision Making 13, 49–58. Mittlböck, M. & Schemper, M. (1996). Explained variation for logistic regression. Statistics in Medicine 15, 1987–1997. Pencina, M.J. & D’Agostino, R.B. (2004). Overall C as a measure of discrimination in survival analysis: model specific population value and confidence interval estimation. Statistics in Medicine 23, 2109–2123. Pencina, M.J., D’Agostino, R.B. Sr, D’Agostino, R.B. Jr & Vasan, R.S. (2008). Evaluating the added predictive ability of a new marker: from area under the ROC curve to reclassification and beyond. Statistics in Medicine 27, 157–172. 72

文献 Pencina, M.J., D'Agostino, R.B. Sr & Steyerberg, E.W. (2011). Extensions of net reclassification improvement calculations to measure usefulness of new biomarkers. Statistics in Medicine 30, 11–21. Pencina, M.J., D'Agostino, R.B. Sr & Song, L. (2012). Quantifying discrimination of Framingham risk functions with different survival C statistics. Statistics in Medicine 31, 1543–1553. Steyerberg, E.W. & Pencina, M.J. (2010). Reclassification calculations for persons with incomplete follow-up. Annals of Internal Medicine 152, 195–196. Steyerberg, E.W., Vickers, A.J., Cook, N.R., Gerds, T., Gönen, M., Obuchowski, N., Pencina, M.J. & Kattan, M.W. (2010). Assessing the performance of prediction models: a framework for traditional and novel measures. Epidemiology 21, 128– 138. Uno, H., Cai, T., Tian, L. & Wei, L.J. (2007). Evaluating prediction rules for t-year survivors with censored regression models. Journal of American Statistical Association 102, 527–537. Uno, H., Cai, T., Pencina, M.J., D'Agostino, R.B. & Wei, L.J. (2011). On the C-statistics for evaluating overall adequacy of risk prediction procedures with censored survival data. Statistics in Medicine 30, 1105–1117. Vickers, A.J. & Elkin, E.B. (2006). Decision curve analysis: a novel method for evaluating prediction models. Medical Decision Making 26, 565–574. Vickers, A.J., Cronin, A.M., Elkin, E.B. & Gönen, M. (2008). Extensions to decision curve analysis, a novel method for evaluating diagnostic tests, prediction models and molecular markers. BMC Medical Informatics and Decision Making 8, 53. 73

付録 SASコード例 SASマクロ、Rパッケージ一覧 ウェブサイト上にアップします http://www.epistat.m.u-tokyo.ac.jp/ 74

• http://www.epistat.m.u-tokyo.ac.jp/