Brier score、NRI 與 IDI

三個分數，補的是 C-index 沒說的話

C-index 只看排序。 校準只看刻度。兩件事都看過之後，論文通常還會再放三個數字：

• Brier score —— 一個把鑑別力與校準揉在一起的整體分數
• NRI（net reclassification improvement，淨重分類改善）—— 新模型把多少人挪到「對的那一格」
• IDI（integrated discrimination improvement，整合鑑別改善）—— 事件組與非事件組的平均預測風險被拉開了多少

三個都算給你看，用的是與決策曲線那一頁完全相同的兩個模型與同一批病人： 在 survival::rotterdam（2982 人、1713 個事件）開發， 在 survival::gbsg（686 人、299 個事件）驗證，時間視野 5 年。

模型 A 的 C-index 高 0.020，但它沒有重新校準： 平均只預測 43.1%，而這群人實際的五年風險是 50.8%——系統性低估。整頁的張力都從這裡長出來。

Brier score：一個數字裡混了兩件事

Brier score 就是預測機率與實際結果之間的均方誤差：

$y_i$ 是 0 或 1，$\hat{p}_i$ 是模型給的機率。值愈小愈好，理論範圍是 0 到 1。

它可以拆成校準項與鑑別項（calibration-refinement 分解，DeGroot & Fienberg 1983；Murphy 1973 那個分解是三項，多一個只跟事件率有關的 uncertainty 項）。所以一個 Brier score 變好， 可能是排序變好，也可能只是刻度被修正——光看這個數字分不出來。而且它沒有尺度。0.228 是好還是壞？ 唯一能回答的方式是跟一個基準比。這一頁用的基準寫在 stats 檔裡： 把這群人自己的 Kaplan-Meier 五年風險 50.8% 發給每一個人。 連這個都贏不了的模型，等於沒有告訴臨床任何它從盛行率不知道的事。

有設限的資料直接算 Brier 是錯的

五年風險模型的麻煩是：有些人根本沒被追蹤到五年。這一批人裡， 686 位有 280 位在五年前就被設限， 佔 40.8%。 406 位五年狀態已知，其中 285 位發生事件、121 位無事件。

最直覺的做法是「只算五年狀態已知的人」。這叫完整病例計分，而它不是雜訊比較多的同一個估計—— 它估的是另一個量：追蹤得完整的人事件率是 70.2%， 而這個世代的五年風險是 50.8%。

正確的做法是 IPCW（inverse probability of censoring weighting，設限機率反向加權）： 先用反向 Kaplan-Meier 估設限分布 $G(t)$，再給每個還有資訊的人一個權重：

五年內出事的人權重是 $1/\hat{G}(T_i^-)$，五年後還在的人是 $1/\hat{G}(t)$， 提前設限的人權重 0——他們沒有被丟掉，是被其他同類的人替他們發言。

Scaled Brier：跟基準比，才有尺度

把 Brier score 除以基準的 Brier score 再拿 1 去減，就得到 scaled Brier（也叫 Brier skill score）：

0 代表「跟每個人都發基準值一樣好」，1 代表完美。它才是可以拿來說「這個模型加了多少值」的那個數字。

C-index 偏好模型 A， Brier score 偏好模型 B。 兩個 scaled Brier 的區間大幅重疊，但區間重疊本身不能拿來下結論——這兩個估計是在同一批病人身上算出來的， 誤差幾乎完全相關，所以要看的是差值自己的區間。把差值放進同一輪 bootstrap 重抽： A 減 B 是 −0.004，95% 區間 −0.045 到 +0.024，跨過零。 正確的說法因此是在這個世代上未偵測到兩個模型的整體準確度差異，不是「B 比較好」。

每個時點都算一次，順便看估計什麼時候不能信

Brier score 是「某一個時間視野」的分數，換一個視野答案就換一次。整條算出來長這樣：

Reclassification table：怎麼排、怎麼讀

NRI 與 IDI 都建立在重分類表（reclassification table）上。做法是先把預測風險切成幾組， 再看每個病人在新舊兩個模型下分別落在哪一組。

排法在這一頁固定成一種：舊模型（B）的分組沿著橫列走，一橫列就是舊模型的一組； 新模型（A）的分組沿著直行走。所以對角線上的人是兩個模型意見一致的， 對角線右上方是被新模型往上推的，左下方是被往下推的。

而且事件組與非事件組一定要分開兩張表。往上推對「後來真的出事的人」是對的， 對「後來沒出事的人」是錯的；混在一張表裡看不出誰被幫到、誰被害到。

這些表格是人數，所以只能包含五年狀態已知的 406 位（686 位裡的）； 提前設限的 280 位放不進整數格子裡。 這件事有後果，本頁最後一節會用 IPCW 加權版把它量出來。

NRI 的算法

把兩張表的「淨移動比例」相加，就是分組式 NRI：

第一項：事件組往上推是對的，所以「上減下」。第二項：非事件組往下推才是對的，所以「下減上」。 兩項各自的範圍是 −1 到 +1，加起來是 −2 到 +2。

NRI 為什麼在這裡是正的

先注意一件事：除了完全沒有人移動的那一組之外，每一組切點的事件組貢獻都是負的， 非事件組貢獻都是大的正數。這不是巧合，是模型 A 系統性低估的直接後果—— 它把幾乎每一個人都往下推一級。

以 15/30 這組切點為例。被往下推的人是 30 位事件者與 38 位非事件者。 前者是被害的（本來就會出事，卻被調降風險），後者是被幫到的。 NRI 的算法是把兩個比例直接相加：

• 事件組：30 除以 285，記成 −0.105
• 非事件組：38 除以 121，記成 +0.314
• 相加得到 +0.209

切點決定一切

同一對模型，同一批病人，換一組切點，NRI 從 +0.000 一路跑到 +0.209，而且最後一組是負的。

淨效益差的單位是「每千位病人，模型 A 比模型 B 多幾個淨真陽性」， 定義與決策曲線那一頁完全相同。

三件事值得逐一看：

第一組切點的 NRI 恰好是 +0.000，因為沒有任何人被移動。 回頭看那組的兩張表就知道： 285 位事件者與 121 位非事件者，兩個模型都把他們全部放進最高的那一組。 兩個切點都低於這群人的風險分布，於是這組切點什麼也沒問到——一組從別的疾病、別的族群抄過來 而沒有對照過本地資料的切點，長的就是這個樣子。

最後一組把 NRI 的正負號翻過來了，而它是看過答案之後才挑的。 stats 檔裡這一組帶著 chosen_to_demonstrate 旗標，上面那兩張表與圖上的星號都是從那個旗標渲染出來的。 把它留在頁面上而不掩飾，是因為這一頁正在批評的就是這種挑法——只報一組切點的論文， 讀者沒有辦法知道作者試過幾組。

同一組切點下，NRI 與淨效益不必同號。 在 50/60 這組的低切點處， 模型 A 的淨效益比模型 B 少 23.55 每千人， 而在高切點處反而多 3.04。

IDI：不用切點，但也不是免費的

IDI 把切點整個拿掉，改比兩組人的平均預測風險：

也可以讀成「鑑別斜率（discrimination slope，事件組平均預測減非事件組平均預測）的增加量」。

IDI 是 0.0400，bootstrap 95% CI 0.0139 到 0.0563。這個區間不跨過零， 所以可以直接說：模型 A 把事件組與非事件組拉得比模型 B 開。

完整病例的表格帶著偏誤：IPCW 加權版對照

上面所有的重分類表都是整數人數，所以放不進提前設限的 280 位。把同樣的 NRI 與 IDI 用 IPCW 權重重算一次， 就看得到那件事付出的代價：

IDI 也一樣：完整病例是 0.0400， IPCW 加權是 0.0371。

差距不大，但方向一致而且不是四捨五入：除了沒有任何人被移動的那一組（兩邊都是零）以外，完整病例版都比加權版更樂觀一點點。加權之後的有效樣本數是 事件組 348.701、非事件組 337.107，比整數表格的 285 與 121 都大——被丟掉的人回來了，而他們回來之後把答案往回拉。

在這份資料上這個偏誤小到不改變結論。它會不會小，取決於設限有多重、以及設限跟預測值有沒有關係， 兩者都是要在論文裡交代的事，不是可以預設的。

動手跑一次

library(survival)
data(cancer, package = "survival")
# 與 B5-05 相同的世代調和與模型配適。攤開寫，這一段才自己跑得動。


H <- 5 * 365.25

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)

# 兩個候選模型，與本頁產圖腳本同一份定義。
BIG   <- c("age", "meno", "size_mm", "grade3", "nodes", "log_pgr", "log_er", "hormon")
SMALL <- c("age", "size_mm", "nodes", "grade3")
rot$log_pgr <- log1p(rot$pgr); rot$log_er <- log1p(rot$er)
ext$log_pgr <- log1p(ext$pgr); ext$log_er <- log1p(ext$er)
keep <- c("rfs_time", "rfs_event", BIG)
dev  <- rot[complete.cases(rot[, keep]), keep]
ext  <- ext[complete.cases(ext[, keep]), keep]

# Breslow 的累積基線風險，直接寫出來——basehaz(centered = FALSE) 與
# predict(type = "lp") 不匹配（後者是中心化的），配錯會讓每個預測風險都少乘一個常數。
breslow_h0 <- function(time, event, lp, h = H) {
  o <- order(time); time <- time[o]; event <- event[o]; r <- exp(lp[o])
  at_risk <- rev(cumsum(rev(r)))
  sum(1 / at_risk[event == 1 & time <= h])
}
form_of <- function(v) as.formula(paste("Surv(rfs_time, rfs_event) ~",
                                        paste(v, collapse = " + ")))
fit_big   <- coxph(form_of(BIG),   data = dev)
fit_small <- coxph(form_of(SMALL), data = dev)
h0_big   <- breslow_h0(dev$rfs_time, dev$rfs_event,
                       as.numeric(as.matrix(dev[, BIG]) %*% coef(fit_big)))
lp_big   <- as.numeric(as.matrix(ext[, BIG])   %*% coef(fit_big))
lp_small <- as.numeric(as.matrix(ext[, SMALL]) %*% coef(fit_small))

# A：八個變項，原封不動搬過來。B：四個變項，再對這個世代重新校準。
pred_A <- 1 - exp(-h0_big * exp(lp_big))
sl_B   <- unname(coef(coxph(Surv(rfs_time, rfs_event) ~ lp_small, data = ext))[1])
h0_B   <- breslow_h0(ext$rfs_time, ext$rfs_event, sl_B * lp_small)
pred_B <- 1 - exp(-h0_B * exp(sl_B * lp_small))


# 1. 設限分布：把 event 指標反過來做 Kaplan-Meier
Gfit <- survfit(Surv(rfs_time, 1 - rfs_event) ~ 1, data = ext)
G_at    <- function(u) summary(Gfit, times = u, extend = TRUE)$surv
G_minus <- function(u) sapply(u, function(x)          # 左極限，給 case 用
  summary(Gfit, times = max(0, x - 1e-8), extend = TRUE)$surv)

# 2. IPCW 權重。提前設限的人權重 0：不是被丟掉，是被同類的人代表
ipcw <- function(d, t) {
  is_case <- d$rfs_time <= t & d$rfs_event == 1
  is_ctrl <- d$rfs_time >  t
  w <- numeric(nrow(d))
  w[is_case] <- 1 / G_minus(d$rfs_time[is_case])
  w[is_ctrl] <- 1 / G_at(t)
  list(w = w, case = is_case, ctrl = is_ctrl)
}

# 3. IPCW Brier，以及它的基準
brier_ipcw <- function(pred, t = H) {
  z <- ipcw(ext, t)
  contrib <- numeric(nrow(ext))
  contrib[z$case] <- (1 - pred[z$case])^2
  contrib[z$ctrl] <- pred[z$ctrl]^2
  mean(z$w * contrib)
}
base <- 1 - summary(survfit(Surv(rfs_time, rfs_event) ~ 1, ext),
                    times = H)$surv
1 - brier_ipcw(pred_A) / brier_ipcw(rep(base, nrow(ext)))   # scaled Brier

# 4. 完整病例版，拿來對照。差別只有 known 這一行
known <- ext$rfs_time > H | (ext$rfs_time <= H & ext$rfs_event == 1)
y     <- as.integer(ext$rfs_time <= H & ext$rfs_event == 1)
mean((y[known] - pred_A[known])^2)      # 基準變成「追蹤得到的人的事件率」，不是世代的風險

# 5. NRI：分組，數上下移動，兩個比例相加
band <- function(p, cuts) cut(p, breaks = c(-Inf, cuts, Inf), labels = FALSE)
nri <- function(cuts) {
  a <- band(pred_A, cuts); b <- band(pred_B, cuts)
  ev <- known & y == 1; ne <- known & y == 0
  (sum(a[ev] > b[ev]) - sum(a[ev] < b[ev])) / sum(ev) +   # 事件組：上減下
  (sum(a[ne] < b[ne]) - sum(a[ne] > b[ne])) / sum(ne)     # 非事件組：下減上
}
sapply(list(c(.10, .20), c(.15, .30), c(.30, .60), c(.50, .60)), nri)
# 換一組切點答案就換一次 —— 只報一組的論文，讀者無從知道試過幾組

# 6. IDI：不用切點，但也不看絕對水準
ev <- known & y == 1; ne <- known & y == 0
(mean(pred_A[ev]) - mean(pred_B[ev])) - (mean(pred_A[ne]) - mean(pred_B[ne]))

import numpy as np
from lifelines import KaplanMeierFitter
# 與 B5-05 相同的世代調和與模型配適。攤開寫，這一段才自己跑得動。

import pandas as pd
from lifelines import CoxPHFitter

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

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"]
beta = (CoxPHFitter()
        .fit(rot[v + ["rfs_time", "rfs_event"]], "rfs_time", "rfs_event")
        .params_[v].to_numpy())

H = 5 * 365.25
G = KaplanMeierFitter().fit(ext["rfs_time"], 1 - ext["rfs_event"])   # 設限分布

def ipcw(t=H):
    case = (ext["rfs_time"].values <= t) & (ext["rfs_event"].values == 1)
    ctrl = ext["rfs_time"].values > t
    w = np.zeros(len(ext))                       # 提前設限的人留 0
    w[case] = 1 / G.predict(ext["rfs_time"].values[case] - 1e-8).values
    w[ctrl] = 1 / float(G.predict(t))
    return w, case, ctrl

def brier_ipcw(pred, t=H):
    w, case, ctrl = ipcw(t)
    c = np.zeros(len(ext))
    c[case] = (1 - pred[case]) ** 2
    c[ctrl] = pred[ctrl] ** 2
    return float(np.mean(w * c))

def nri(cuts, pred_new, pred_old, y, sel):
    a, b = np.digitize(pred_new, cuts), np.digitize(pred_old, cuts)
    ev, ne = sel & (y == 1), sel & (y == 0)
    return ((a[ev] > b[ev]).sum() - (a[ev] < b[ev]).sum()) / ev.sum() + \
           ((a[ne] < b[ne]).sum() - (a[ne] > b[ne]).sum()) / ne.sum()

讀論文時要問的五件事

1. Brier score 旁邊有沒有基準？ 沒有基準的 Brier score 沒有尺度。 最低限度要有 scaled Brier，或是「每個人都給盛行率」那個參考模型的分數。
2. 有設限的資料，作者怎麼處理？ 沒提 IPCW 而直接報 Brier score， 多半就是完整病例計分，而那可能連正負號都不對。
3. NRI 用的切點哪裡來的？ 事前指定、有臨床理由、而且只報一組嗎？ 報了好幾組而只討論最好看的那一組，跟只報一組但試過好幾組，在紙面上長得一樣。
4. NRI 的兩個組成分有沒有分開報？ 只報總和會藏掉「事件組其實變差」這種事， 而那正是本頁四組切點裡的三組。
5. 有沒有一起報決策曲線或淨效益？ NRI 是一個沒有單位的加總， 淨效益有單位而且綁在一個明講出來的臨床交換條件上。

常見誤用

這一頁的結論

同一對模型、同一批病人、同一個五年視野：

• C-index 偏好模型 A
• Brier score 與 scaled Brier 偏好模型 B
• NRI 在四組切點下從 +0.000 到 +0.209，事後挑的那組是 −0.045
• 淨效益差在整條門檻網格上變號 10 次，而且區間普遍跨過零

沒有一個分數是錯的，它們只是在回答不同的問題，而其中大部分不是臨床要問的那一個。 臨床要問的是「照這個模型做決定，比現在的做法好嗎」， 那個問題由決策曲線回答， 而它的前提是模型的刻度先是對的——那是校準那一頁的事。 Brier、NRI 與 IDI 適合放在那兩件事旁邊，當作補充的描述，不適合拿來取代它們。

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B5-10-brier-nri-idi.R

回報這一頁的問題