隨時間變動的區辨力

一個 C 統計量看不出它自己漏掉了什麼

預測模型那一章與它底下的 內部驗證、 外部驗證、 校準三頁，整條線都建立在同一個 Cox 模型上， 而區辨力那一欄從頭到尾只有一個數字：Harrell’s C。

那個數字是這樣算出來的：把每一對「可以比較先後」的病人拿出來，看模型有沒有把先發生事件的 那位排在前面，然後對所有這樣的配對取平均。問題就在「所有」兩個字—— 它把兩年後的區辨力、五年後的區辨力、以及一對都拖到第六年才發病的病人，全部平均成一個數。

但臨床決策不發生在「所有可比配對的平均」上。 你要跟一位剛開完刀的乳癌病人講的是 「五年之內復發或死亡的機率」，或者是「兩年之內」——那是一個特定時點的問題。 審稿人問「五年的 AUC 是多少」，問的也是同一件事。

這一頁把那一個數字攤成一條曲線。

AUC(t) 在問的是一個時點上的問題

先把時間軸上的某一個 $t$ 釘住，例如五年。釘住之後，每個人在這個時點只有兩種身分：

• case：在 $t$ 之前已經發生事件（本頁是復發或死亡，取先發生者）
• control：到了 $t$ 仍然沒有發生事件

有了 case 與 control，就有了一張標準的二分類問題，於是也就有了 ROC 曲線與 AUC 的全套機器。這個版本叫 cumulative/dynamic AUC：case 是「累積到 $t$ 為止」發生事件的人（cumulative）， control 是「$t$ 這個時刻仍然無事件」的人（dynamic）。寫成機率就是

其中 $\hat r$ 是模型給的風險分數（Cox 模型就是線性預測值）。這與 那一頁的機率解釋是同一句話， 差別只在「誰是病人、誰不是」這件事現在隨 $t$ 改變：同一個人在兩年時是 control， 在五年時可能已經變成 case。

這一頁沿用外部驗證那一頁的同一個模型

不重新認識一份資料，是刻意的：這一頁的整個論證是「一條曲線 vs 一個數字」， 而那個數字必須是讀者在外部驗證那一頁 已經見過的同一個 Harrell’s C，比較才誠實。

開發集是 survival::rotterdam（2982 位病人、 1713 個 recurrence-free survival 事件）， 驗證集是 survival::gbsg（686 位、 299 個事件，最長追蹤 7.28 年）。 Cox 模型的四個共變項是：

• 年齡（手術當時，連續變項）
• 腫瘤大小（毫米）；rotterdam 記的是分級區間，取各級中點以配合 gbsg 的毫米
• 陽性淋巴結數（連續變項）
• 腫瘤分級（第三級以上為 1，否則為 0）；兩個世代分級的細緻度不同，所以併成二分變項

把區辨力攤開來看

為什麼整條曲線都在 Harrell’s C 上面

讀者看這張圖的第一個反應，通常是「那兩個數字哪一個算錯了」。 網格上 11 個時點裡，有 11 個時點的 AUC(t) 點估計高於同一個模型的 Harrell’s C（0.642）；連曲線最低的那一點（4.5 年） 都還高出 0.028，落差最大處是 0.069。

沒有人算錯，因為兩者估的根本不是同一個量。

• AUC(t) 只問一個問題：一位在 $t$ 之前發病的病人，分數有沒有比一位在 $t$ 時仍無事件的 病人高。它從來不比較兩位都在 $t$ 之前發病的人誰先誰後——在 $t$ 這個時點上，他們是同一類。
• Harrell’s C 會比。它對所有可比配對取平均，其中包含「一位在第五年發病、一位在第六年發病」 這種配對。要把兩位預後都很差、只差幾個月的病人排對，接近擲銅板——而這類配對在長追蹤的 世代裡數量很多，它們把平均值往 0.5 拉。

還有第二個來源：Harrell’s C 的計算會受到設限分布的影響（哪些配對「可比」取決於誰何時被設限）， 而本頁的 AUC(t) 已經用 IPCW 把設限的影響校正掉了。

不過方向一致不等於每一點都分得開：11 個時點裡， 有 8 個時點的 95% 信賴區間仍然把 Harrell’s C 包在裡面。 所以這一段話的正確版本是「AUC(t) 系統性地落在 Harrell’s C 上方」， 而不是「每一個時點都與 Harrell’s C 有統計上可分辨的差距」——後者本資料撐不起來。

平坦本身就是這一頁的結論

現在誠實面對這條曲線的形狀。

在這個模型、這個外部世代上，11 個時點的 AUC(t) 落在 0.669（4.5 年）到 0.711（6.0 年）之間，全距 0.042， 標準差 0.014。

但只引全距會高估這條曲線的起伏。 最低點與最高點的 95% 信賴區間是重疊的（0.620–0.719 對上 0.630–0.792）—— 也就是說整條曲線的上下擺動落在它自己的不確定度之內，本資料未偵測到 AUC(t) 隨時間變動。 這句話不等於「AUC(t) 不隨時間變動」：11 個時點的區間都不窄， 中等幅度的變化它一樣排除不掉。

臨床論文實際會引用的三個時點差距更小：

這三個時點彼此相差 0.020。

右端的信賴區間為什麼變寬

整條網格攤開來看，會看到全距是怎麼來的：

最後一欄與倒數第二欄解釋了圖上那個往右變寬的喇叭口：時間愈往後， 還留在風險集合裡的人愈少（從 602 位掉到 36 位）， control 那一邊可用的人就愈少，而 IPCW 的權重同時愈來愈大—— 兩件事都會讓估計變不穩。整條曲線最寬的區間出現在 6.0 年，寬 0.162， 比前面那個全距 0.042 還大了好幾倍。

而 6.0 年剛好就是全域最大值所在。 換句話說，那個全距有相當一部分是這一點的雜訊撐出來的，不是模型在後期真的變厲害了。 讀 AUC(t) 曲線的第一個習慣，就是先看那個時點還剩幾個人，再看曲線的高低。

固定一個時點的 ROC 長什麼樣子

AUC(t) 是一個面積，面積底下有一條曲線。把 $t$ 釘在 5.0 年畫出來：

這張圖與一般的 ROC 曲線看起來一模一樣， 但底下的資料結構差了三件事，而三件都會影響你怎麼讀它：

動手跑一次

library(survival)
library(timeROC)
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")
fit <- coxph(Surv(rfs_time, rfs_event) ~ age + size_mm + nodes + grade3, data = rot)

# 風險分數：自己算 X %*% beta，不用 predict() 的置中版本（與 B5-05 一致）
lp <- as.numeric(as.matrix(ext[, v]) %*% coef(fit))

# 一個數字
concordance(Surv(ext$rfs_time, ext$rfs_event) ~ lp, reverse = TRUE)

# 一條曲線。weighting = "marginal" 就是 Kaplan-Meier 的設限權重（IPCW）
grid <- seq(1, 6, by = 0.5) * 365.25
roc  <- timeROC(T = ext$rfs_time, delta = ext$rfs_event, marker = lp,
                cause = 1, times = grid, weighting = "marginal", iid = TRUE)

auc <- as.numeric(roc$AUC)
se  <- as.numeric(roc$inference$vect_sd_1)   # iid = TRUE 才有
cbind(years = grid / 365.25, auc = auc,
      lcl = auc - 1.96 * se, ucl = auc + 1.96 * se)

# 每個時點還剩幾個人——讀曲線之前先看這個
sapply(grid, function(t) sum(ext$rfs_time >= t))

# 某個時點的 ROC 座標，可以自己畫
k <- which(roc$times == 5 * 365.25)
plot(roc$FP[, k], roc$TP[, k], type = "l", xlim = c(0, 1), ylim = c(0, 1))
abline(0, 1, lty = 2)

import numpy as np
from sksurv.util import Surv
from sksurv.metrics import cumulative_dynamic_auc, concordance_index_ipcw
import pandas as pd
from lifelines import CoxPHFitter

# 與 B5-05 相同的世代調和與模型配適。攤開寫，這一段才自己跑得動。
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())

# scikit-survival 吃的是結構化陣列：(event 為 bool, time)
y_train = Surv.from_arrays(event=rot["rfs_event"].astype(bool), time=rot["rfs_time"])
y_test  = Surv.from_arrays(event=ext["rfs_event"].astype(bool), time=ext["rfs_time"])

v  = ["age", "size_mm", "nodes", "grade3"]
lp = ext[v].to_numpy() @ beta          # beta 來自任一個配好的 Cox 模型；分數愈大風險愈高

times = np.arange(1, 6.5, 0.5) * 365.25    # 時點必須落在驗證世代的追蹤範圍內
auc, mean_auc = cumulative_dynamic_auc(y_train, y_test, lp, times)

# Uno's C：同樣走 IPCW 的單一摘要數字
print(concordance_index_ipcw(y_train, y_test, lp))

這一頁與其他頁的關係

• Harrell’s C 本身怎麼算、為什麼表面表現一定偏高 —— 見 內部驗證。本頁只是把同一個量攤到時間軸上， 樂觀偏誤那一整套問題並沒有因此消失：在開發集自己身上算 AUC(t)，一樣是偏高的。
• 本頁的模型、世代與那個 Harrell’s C 從哪裡來 —— 見 外部驗證。那一頁還示範了另一件事： 光看鑑別力分不出三種驗證設計的優劣，差別要看校準才顯現。
• 區辨力好不等於機率算得準 —— 見校準。 AUC(t) 與 Harrell’s C 一樣只看排序：把每個人的預測風險同乘一個常數， 這一頁的曲線一位小數都不會動，而每個人的預測機率都錯了。 一條漂亮的 AUC(t) 曲線不能替校準背書。
• 一般的 ROC 與 AUC —— 見 ROC 曲線與 AUC。 兩者的機率解釋是同一句話，差別在 case 與 control 的身分會不會隨時間改變。
• 設限到底是什麼、獨立設限假設在假設什麼 —— 見 設限與追蹤時間。IPCW 是把那條假設用進估計裡， 不是把它拿掉。
• 整章的臨床脈絡 —— 見預測模型研究。

讀論文裡的 AUC(t) 時要問的四件事

1. 時點是哪些，為什麼是這些？ 事後挑一個 AUC 最漂亮的時點來報， 與挑切點是同一種問題（見切點選擇）。時點應該由臨床決定。
2. 那個時點還剩幾個人在風險集合裡？ 論文很少寫，但這決定了信賴區間的寬度。 沒有附風險人數就把最右端那一點當結論，是這類圖最常見的誤讀。
3. 設限怎麼處理的？ 有沒有寫 IPCW（或等價的加權）。只寫「以 ROC 分析評估」的， 多半是把設限的人當 control 或直接刪掉。
4. 有沒有一起報校準？ 只報 AUC(t) 的外部驗證，仍然只報了 三件套的三分之一。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B5-08-time-dependent-auc.R

回報這一頁的問題