Cox 比例風險模型

這個模型在解決什麼問題

Kaplan-Meier 與 log-rank 可以告訴你兩條曲線不一樣，但給不了兩件事：差多少，以及扣掉其他因素之後還差多少。

真實的臨床問題幾乎都是後者。女性肺癌病人存活比男性長——這是年齡的效果嗎？是體能狀態（performance status）的效果嗎？還是性別本身？而且分組變項如果是連續的（年齡、腫瘤大小、eGFR），log-rank 根本無從下手，你只能硬把它切成幾組，切點還是自己選的。

Cox 比例風險模型（Cox proportional hazards model）處理的就是這件事：同時放進多個變項（連續的、類別的都行），對每一個給出一個效果量，而且不需要假設存活時間服從任何特定分布。

核心概念：風險函數與 HR

先定義風險函數（hazard function）$h(t)$：一個到 $t$ 這一刻為止還沒發生事件的人，在這一瞬間發生事件的速率。

這個式子有兩個關鍵：分子裡的條件 $T \ge t$（只算還在風險集合裡的人），以及它是個速率不是機率，可以大於 1。

Cox 模型寫成：

拆開來看是兩塊乘在一起：

• $h_0(t)$ 是基線風險（baseline hazard），所有變項都等於 0 的那個人的風險隨時間怎麼跑。它可以是任何形狀——先高後低、有峰、有波動都行。
• $\exp(\beta_1 x_1 + \cdots)$ 是一個不隨時間變的倍數，把基線風險整條往上或往下拉。

於是兩個人的風險比值是：

$h_0(t)$ 上下抵消掉了。這就是 hazard ratio（HR，風險比），也解釋了兩件事：為什麼 $h_0(t)$ 不用指定（它在比值裡消失了，這叫半參數 semi-parametric——效果的部分是參數化的，時間的部分不是），以及為什麼這叫「比例風險」（proportional hazards）——這個比值被假設在整段時間裡都是同一個數。

動手跑一次

library(survival)
data(cancer, package = "survival")   # lung、colon、mgus2 都在這包裡

lung$sex_f <- factor(lung$sex, levels = c(1, 2),
                     labels = c("Male", "Female"))

fit <- coxph(Surv(time, status) ~ age + sex_f + ph.ecog + wt.loss,
             data = lung)
summary(fit)        # coef、exp(coef) = HR、95% CI、p、Concordance

# 年齡改成「每 10 歲」的對比：要在係數尺度上乘，不是把 HR 乘 10
b  <- coef(fit)["age"]; se <- sqrt(vcov(fit)["age", "age"])
exp(10 * c(b, b - 1.96 * se, b + 1.96 * se))

# 大一點的真實試驗：colon 輔助化療，etype 2 = 死亡
colon_d <- subset(colon, etype == 2)
colon_d$rx <- factor(colon_d$rx, levels = c("Obs", "Lev", "Lev+5FU"))
fit2 <- coxph(Surv(time, status) ~ rx + sex + age + nodes + extent + differ,
              data = colon_d)
summary(fit2)

# 校正後的存活曲線：三組都代入同一組共變項值
ref <- data.frame(rx = factor(levels(colon_d$rx), levels = levels(colon_d$rx)),
                  sex = 1, age = 60, nodes = 2, extent = 3, differ = 2)
plot(survfit(fit2, newdata = ref), lwd = 2,
     col = c("#8c6d4a", "#4d6a8c", "#c44d4d"),
     xlab = "Days", ylab = "Adjusted survival")

import statsmodels.api as sm
from lifelines import CoxPHFitter

# 命名陷阱：lung 在 Rdatasets 上掛在 "cancer" 這個 help page 底下
lung = sm.datasets.get_rdataset("cancer", "survival").data
d = lung[["time", "status", "age", "sex", "ph.ecog", "wt.loss"]].dropna().copy()
d["event"]  = (d["status"] == 2).astype(int)
d["female"] = (d["sex"] == 2).astype(int)
d = d.drop(columns=["status", "sex"])

cph = CoxPHFitter().fit(d, duration_col="time", event_col="event")
cph.print_summary()          # coef、exp(coef)、CI、p、concordance
cph.plot()                   # 就是 forest plot

# colon 與 mgus2 沒有掛在 Rdatasets 上（同樣是 help page 命名的問題），
# 要在 Python 用它們，先在 R 跑一次：
#   write.csv(subset(colon, etype == 2), "colon_death.csv", row.names = FALSE)

怎麼讀報表

summary(coxph(...)) 印出來的東西，對應到論文 Table 的每一欄。用上面那個 lung 模型逐格看——先說清楚這個模型是用多少人跑的：資料集有 228 人，其中 15 人因為共變項有缺失被完整個案分析排除，實際進入模型的是 213 人、151 個事件。

（掉了幾個人要講出來。coxph() 不會警告你，它只是安靜地少算幾列；而缺失往往與預後有關，所以這不是無害的四捨五入。世代研究那一章有整節在談這件事。）

逐欄怎麼講：

• HR 的方向。 大於 1 是風險升高、小於 1 是降低。女性相對男性 HR = 0.55，讀作「在同年齡、同 ECOG、同體重減輕的條件下，女性任一時刻的死亡風險約為男性的 0.55 倍」。「在其他變項相同的條件下」這句話不是客套話，是 HR 的定義的一部分，拿掉就是另一個數字了。
• 連續變項的單位。 年齡的 HR = 1.013，看起來小到沒意義——因為那是每多一歲。換成每 10 歲是 HR 1.14（0.95–1.38）。換算要在係數尺度做（$\exp(10\beta)$），不能把 HR 乘以 10。論文沒寫單位的連續變項 HR 是無法解讀的。
• 信賴區間比 p 值有用。 年齡的區間是 0.995–1.033，跨過 1；體重減輕是 0.978–1.004，也跨過 1。
• Concordance（Harrell’s C）= 0.647。 隨機挑兩個人，模型把「先發生事件的那個」給較高風險的機率。0.5 是亂猜，1 是完美。0.647 意思是這個模型有訊號但區辨力普通——這在臨床預後模型裡很常見，也是為什麼「p 值很顯著」跟「模型能不能拿來預測個別病人」是兩件事。
• 整體模型檢定（likelihood ratio test）= 31.0，df = 4，p < 0.001。 問的是「全部係數同時為 0 嗎」，相當於迴歸裡的整體 F 檢定。

多變項校正實際上在做什麼

「校正年齡」聽起來像是把年齡的影響減掉，但模型裡發生的事更具體：在每一個事件時點，把該時點風險集合裡每個人的預測風險依各自的共變項加權，再問「實際發生事件的是不是預測風險較高的那個」。 係數就是讓這件事最一致的那組數字（技術上叫偏概似 partial likelihood）。

結果是：報出來的性別 HR，是「同年齡、同 ECOG、同體重減輕」的兩個人的比較。這帶來三個實務後果：

1. 校正的變項換了，HR 就換了。 兩篇論文報不同的 HR，可能只是模型裡的變項不同，不是資料衝突。看 HR 一定要連著看那個模型放了什麼。
2. 不能把中介因子（mediator）放進去。 如果治療是透過降低發炎指標起作用，把發炎指標校正進去，等於把治療的效果扣掉一部分。要決定放什麼，靠的是因果圖（DAG）不是 p 值。
3. 事件數決定你能放幾個變項。 傳統經驗法則是每個變項至少 10 個事件（events per variable, EPV），近年的模擬研究認為視情況可以放寬，但方向不變：151 個事件的模型放四個變項還好，放十五個就是在配雜訊。

用一個真的臨床試驗看校正後的樣子——survival::colon 是大腸癌術後輔助化療試驗。它的原始檔是雙端點的長格式，每位病人各有復發與死亡兩列，一共 1858 列；本頁先用 etype == 2 篩出死亡端點的 929 列，再扣掉共變項缺失，模型用的是 888 人、430 例死亡：

這張表把上一段那個警告變成了具體例子：Levamisole 單用的 HR 是 0.913，區間 0.731–1.141 跨過 1 —— 本試驗未偵測到單用 Levamisole 與死亡風險的關聯。而 Levamisole + 5-FU 的 HR 是 0.673（0.530–0.854），區間完全落在 1 的左邊。歷史上這個試驗正是 5-FU 併用方案成為大腸癌標準輔助治療的依據之一。

Ties：兩個人同一天發生事件

Cox 模型的推導假設事件時間是連續的，也就是不會有兩個人在完全同一刻發生事件。真實資料當然會——時間記到天、記到月、記到年，同時發生（ties）就一定出現。R 提供三種處理方式：

這一節換回完整的 colon 死亡端點資料（不做完整個案排除，因為 ties 的比較只牽涉時間與治療組，不需要那些有缺失的共變項），所以事件數是 452 而不是上面模型表的 430——同一頁上兩個不同的數字，是因為分析集不同，不是其中一個算錯。

時間以天為單位時，452 個死亡事件散在 409 個不同的日子，三種方法幾乎沒差：

但把同一份資料的時間粗化成「年」（452 個事件擠進 8 個時點），差距就出來了：

用之前要檢查的假設

Cox 模型有三個假設，其中最容易被忽略、後果最嚴重的是第一個：

1. 比例風險 —— HR 在整段追蹤期間是同一個數。違反的話，那個單一 HR 是把不同時期相反的效果平均掉的產物，本身沒有清楚的意義。檢查方式與補救辦法在比例風險假設與 Schoenfeld residuals。
2. 連續變項與 log hazard 呈線性 —— 年齡對死亡風險的作用未必是直線。檢查方式是 martingale residuals 對該變項作圖，或直接改用限制性立方樣條（restricted cubic spline）。
3. 非資訊性設限 —— 與所有存活方法共用的地基，見 Censoring 與存活資料的結構。

另外，有競爭風險時，cause-specific Cox 回答的問題與 Fine-Gray 不同，這一點在競爭風險：CIF 與 Fine-Gray 裡談。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B3-03-cox.R

延伸觀看

▶Hazard Ratios – Best explanation for beginners

▶Cox Regression [Cox Proportional Hazards]

▶COX REGRESSION and HAZARD RATIOS

▶存活分析（Survival Analysis）第二部分

回報這一頁的問題