校準

為什麼校準比鑑別力更值得擔心

鑑別力（discrimination）問的是排序：模型有沒有把出事的人排在沒出事的人前面。 校準（calibration）問的是數值：模型說 30% 的那群人，是不是真的有三成出事。

只要臨床決策長成「排序然後挑前面幾個」（誰先做電腦斷層、器官怎麼分配），鑑別力就夠。 但只要決策綁在一個絕對機率門檻上——「十年心血管風險超過 20% 就給 statin」、 「復發風險低於某個數就不加輔助治療」——校準不良的模型會系統性地把病人分錯邊， 而 C-index 完全看不出來。

還有一個現實層面的理由：校準比鑑別力更容易壞，也更容易修。 族群的基準風險會因為轉診模式、篩檢政策、治療進展而改變， 而這些改變幾乎不會改變變項之間的排序關係。所以外部驗證最常見的畫面就是本頁的畫面—— 在此驗證世代，C-index 為 0.642；但平均預測五年風險比觀察風險低 7.21 個百分點（O/E 1.165）， 校準斜率為 0.691（95% CI 0.539–0.842）。

四個層次，愈後面愈嚴格

在真實資料上逐層量一次

模型在 survival::rotterdam（2982 人）開發， 在 survival::gbsg（686 人、299 個事件）驗證， 時間視野 5 年。鑑別力是 C-index 0.642。

第一層：平均校準

平均預測風險 43.6%， 實際觀察（Kaplan-Meier）50.8%。 O/E 比值 1.165，差距 7.21 個百分點。

O/E 大於 1 代表這群人出的事比模型預期的多，也就是模型低估了風險。

第二層：弱校準

校準斜率 0.691（95% CI 0.539–0.842）。

存活模型沒有一個字面上的「截距」，它的對應物是基準累積風險。 把斜率固定在 1、只讓基準風險自由，新的基準風險是原來的 1.211 倍，取對數是 0.192——這個數字扮演的角色就是校準截距， 正值代表低估。

第三層：中度校準

校準斜率小於 1 代表什麼

斜率是把驗證世代的線性預測值當成唯一的預測變項重配一次模型，得到的那個係數。

• 斜率 ＝ 1：模型給出的風險差異幅度剛好。
• 斜率 < 1：風險差異太極端——高風險的人被推得太高、低風險的人被壓得太低。
• 斜率 > 1：反過來，模型給出的差異太保守、太集中在平均值附近。

小於 1 是外部驗證裡最常見的形態，它有兩個來源，而且校準斜率本身分不出是哪一個：

1. 開發時的過度配適——係數本來就太大（見收縮那一頁）
2. 兩個族群的風險結構不同——某個預測變項在新族群裡的效應本來就比較小

不要用 Hosmer-Lemeshow 檢定

它是最常出現在論文裡的「校準檢定」，而方法學文獻建議不要再用它。三個理由：

1. 它靠人為分組——組數與切法會改變結果。
2. 它給的 p 值不帶方向也不帶幅度——不顯著不代表校準好，只代表沒偵測到偏離； 顯著也不告訴你是高估還是低估、差多少。
3. 它的檢定力很低——樣本小的時候幾乎一定不顯著，於是「p 值不顯著」被寫成「校準良好」， 而那正好是樣本不足時最不該下的結論。

那 calibration belt 呢

它在義大利與重症醫學文獻裡很常見（GiViTI calibration belt，R 套件 givitiR）， 形式上是把觀察結果對預測風險配一條多項式，在整段預測風險上畫出那條曲線的 同時信賴區域（simultaneous confidence region），再從那個區域導出一個 p 值—— 區域在任何一段離開對角線就顯著。

它與這一頁第三層的彈性校準曲線問的是同一個問題，但它多給的不只是一個 p 值， 而是曲線本身的不確定性——而這一頁的圖完全沒有畫這一層。 上面那條 spline 曲線是點估計，旁邊沒有任何帶子， 所以讀者無法分辨圖上那個凹陷是真的校準偏離，還是幾百個事件配出來的 spline 在晃。 belt 的加值就落在這個缺口上。p 值只是那個區域的副產品，而上面三點勸退的是那個副產品。

這一頁沒有畫 belt，理由不是「belt 沒有加值」。 一部分是範圍： 這一頁講的是四個層次與重新校準，而「彈性曲線的不確定性」這一層，本頁根本沒有畫—— 上面那個校準斜率是有給信賴區間的，但曲線本身只是一條點估計。 另一部分是對不上——belt 在原始文獻與 givitiR 裡都是為二元結果與其預測機率定義的， 而這一頁驗證的是帶設限追蹤的 Cox 模型，這份驗證世代裡有相當一部分人根本沒被追滿五年。 把它壓成「五年有沒有事件」的二元變項，等於不是把這些人丟掉、就是把他們算成沒事件。 所以看到論文報 belt 不必反對，但要看的是那個區域在哪一段離開對角線、離開多少，不是旁邊那個 p。

重新校準：三個層次

校準壞掉的預設反應不是丟掉模型，是更新它。從輕到重：

在整個驗證世代上各做一次：

不論做哪一層，更新之後都需要再驗證一次，而且不能在拿來更新的那批人身上驗。

動手跑一次

library(survival); library(splines)
data(cancer, package = "survival")

rot <- rotterdam
rot$rfs_time  <- pmin(rot$rtime, rot$dtime)
rot$rfs_event <- as.integer(rot$recur == 1 | rot$death == 1)
rot$size_mm   <- c("<=20" = 15, "20-50" = 35, ">50" = 60)[as.character(rot$size)]
rot$grade3    <- as.integer(rot$grade >= 3)

ext <- gbsg
ext$rfs_time  <- ext$rfstime; ext$rfs_event <- ext$status
ext$size_mm   <- ext$size;    ext$grade3    <- as.integer(ext$grade >= 3)

v   <- c("age", "size_mm", "nodes", "grade3")
H   <- 5 * 365.25
fit <- coxph(Surv(rfs_time, rfs_event) ~ age + size_mm + nodes + grade3, data = rot)

# Breslow 基準累積風險，自己寫。理由見上面的紅字提醒
breslow_h0 <- function(time, event, lp, h) {
  o <- order(time); time <- time[o]; event <- event[o]; r <- exp(lp[o])
  sum(1 / rev(cumsum(rev(r)))[event == 1 & time <= h])
}
lp_dev <- as.matrix(rot[, v]) %*% coef(fit)
h0     <- breslow_h0(rot$rfs_time, rot$rfs_event, lp_dev, H)

lp   <- as.numeric(as.matrix(ext[, v]) %*% coef(fit))
pred <- 1 - exp(-h0 * exp(lp))

# 第一層：平均校準
km  <- survfit(Surv(rfs_time, rfs_event) ~ 1, data = ext)
obs <- 1 - summary(km, times = H, extend = TRUE)$surv
c(mean_predicted = mean(pred), observed = obs, OE = obs / mean(pred))

# 第二層：斜率
cal <- coxph(Surv(rfs_time, rfs_event) ~ lp, data = ext)
coef(cal); confint(cal)

# 第三層：彈性校準曲線
flex <- coxph(Surv(rfs_time, rfs_event) ~ ns(lp, df = 3), data = ext)
h0f  <- breslow_h0(ext$rfs_time, ext$rfs_event,
                   as.numeric(predict(ns(lp, df = 3), newx = lp) %*% coef(flex)), H)
g    <- seq(quantile(lp, .02), quantile(lp, .98), length.out = 100)
obs_curve  <- 1 - exp(-h0f * exp(as.numeric(predict(ns(lp, df = 3), newx = g) %*% coef(flex))))
pred_curve <- 1 - exp(-h0 * exp(g))
plot(pred_curve, obs_curve, type = "l", xlim = c(0, 1), ylim = c(0, 1)); abline(0, 1, lty = 2)

# 重新校準第一層：只換基準風險
h0_new <- breslow_h0(ext$rfs_time, ext$rfs_event, lp, H)
mean(1 - exp(-h0_new * exp(lp)))     # 現在應該貼近 obs

import numpy as np, pandas as pd
from lifelines import CoxPHFitter, KaplanMeierFitter

RD = "https://vincentarelbundock.github.io/Rdatasets/csv/"
rot = pd.read_csv(RD + "survival/rotterdam.csv")
ext = pd.read_csv(RD + "survival/gbsg.csv")

# 與上面 R 端相同的衍生欄位。rotterdam 原始欄位沒有這些，少了這一段
# 下面每一行都會 KeyError。
rot["rfs_time"]  = rot[["rtime", "dtime"]].min(axis=1)
rot["rfs_event"] = ((rot["recur"] == 1) | (rot["death"] == 1)).astype(int)
rot["size_mm"]   = rot["size"].map({"<=20": 15, "20-50": 35, ">50": 60})
rot["grade3"]    = (rot["grade"] >= 3).astype(int)
ext["rfs_time"]  = ext["rfstime"]
ext["rfs_event"] = ext["status"]
ext["size_mm"]   = ext["size"]
ext["grade3"]    = (ext["grade"] >= 3).astype(int)
v = ["age", "size_mm", "nodes", "grade3"]; H = 5 * 365.25

fit = CoxPHFitter().fit(rot[v + ["rfs_time", "rfs_event"]], "rfs_time", "rfs_event")
lp  = ext[v].to_numpy() @ fit.params_[v].to_numpy()

def breslow_h0(time, event, lp, h):
    o = np.argsort(time); t, e, r = time[o], event[o], np.exp(lp[o])
    at_risk = np.cumsum(r[::-1])[::-1]
    return np.sum(1.0 / at_risk[(e == 1) & (t <= h)])

lp_dev = rot[v].to_numpy() @ fit.params_[v].to_numpy()
h0 = breslow_h0(rot["rfs_time"].to_numpy(), rot["rfs_event"].to_numpy(), lp_dev, H)
pred = 1 - np.exp(-h0 * np.exp(lp))

km  = KaplanMeierFitter().fit(ext["rfs_time"], ext["rfs_event"])
obs = 1 - float(km.predict(H))
print(pred.mean(), obs, obs / pred.mean())

tmp = ext[["rfs_time", "rfs_event"]].copy(); tmp["lp"] = lp
print(CoxPHFitter().fit(tmp, "rfs_time", "rfs_event").params_["lp"])   # 校準斜率

讀論文時要問的四件事

1. 有沒有校準圖？ 只報 C-index 或 AUC 的預測模型論文，等於只交了一半。
2. 報的是哪一層？ 只報 O/E 是第一層；只報斜率是第二層；有彈性曲線才到第三層。
3. 驗證世代有多少事件？ 少於一百時，校準的結論不論正負都要打折看。
4. 有沒有用 Hosmer-Lemeshow 的 p 值宣稱「校準良好」？ 那是不成立的推論。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B5-06-calibration.R

延伸觀看

▶VALIDATING PREDICTION MODELS – what is discrimination and calibration?

▶SEER 數據之臨床預測模型 課時09 利用校準圖評價模型

回報這一頁的問題