比例風險假設與 Schoenfeld residuals

這個假設在假設什麼

Cox 模型的形式是「基線風險乘上一個倍數」：

注意 $\beta$ 沒有帶 $t$。這代表那個倍數在整段追蹤期間都是同一個數：下面那個 lung 模型估出來的女性 HR 是 0.602，模型的意思是第 30 天、第 300 天、第 900 天，女性的死亡風險都是男性的 0.602 倍。這就是比例風險假設（proportional hazards assumption）。

它不是一個可以「附帶檢查一下」的技術細節，而是那個 HR 有沒有意義的前提。假設不成立時，coxph() 還是會吐出一個數字，但那個數字是把不同時期不同大小（甚至方向相反）的效果，用一個跟事件分布有關的權重平均起來的產物。它不對應到任何一個時點的真實風險比，而且它的大小會隨著追蹤期長短而變——同一個治療追蹤三年跟追蹤十年會得到不同的 HR，即使真實情況完全一樣。

臨床上假設會破功的典型場景：

怎麼檢查：Schoenfeld residuals

Schoenfeld 殘差（Schoenfeld residuals）的想法很直接：在每一個事件時點，模型會依各人的共變項算出「誰最可能是這次發生事件的人」的預測值；殘差就是實際發生事件那個人的共變項值，減掉模型在該時點的預測值。

關鍵在於：如果比例風險成立，這些殘差與時間應該沒有關係——早期跟晚期的殘差都該在 0 附近隨機散布。如果殘差隨時間有系統性的趨勢，代表 $\beta$ 其實在隨時間變化。

實務上用的是尺度化 Schoenfeld 殘差（scaled Schoenfeld residuals），它有一個很好的性質：把它對時間作圖，那條平滑曲線就是 $\hat{\beta}(t)$ 的估計——直接看得到係數隨時間跑去哪裡。cox.zph() 做的正式檢定，就是在檢定這條線的斜率是不是 0。

動手跑一次

library(survival)
data(cancer, package = "survival")
lung$sex_f <- factor(lung$sex, levels = c(1, 2),
                     labels = c("Male", "Female"))

fit <- coxph(Surv(time, status) ~ age + sex_f + ph.karno + wt.loss,
             data = lung)

zph <- cox.zph(fit)
zph                      # 每個變項一列 + GLOBAL
par(mfrow = c(2, 2)); plot(zph)   # 平滑線就是 beta(t)

# log(-log) 圖：類別變項的圖形檢查，平行 = 假設沒問題
plot(survfit(Surv(time, status) ~ sex_f, data = lung),
     fun = "cloglog", col = c("#4d6a8c", "#c44d4d"), lwd = 2)

# 補救一：把時間切開，讓係數在前後兩段各自估
sp <- survSplit(Surv(time, status) ~ ., data = lung, cut = 180,
                episode = "period")
coxph(Surv(tstart, time, status) ~ age + sex_f + wt.loss +
        ph.karno:strata(period), data = sp)

# 補救二：讓係數是時間的函數 beta(t) = b0 + b1*log(t)
coxph(Surv(time, status) ~ age + sex_f + wt.loss + ph.karno + tt(ph.karno),
      data = lung, tt = function(x, t, ...) x * log(t))

# 補救三：分層——每一層有自己的 baseline hazard
lung$kstrata <- cut(lung$ph.karno, c(-Inf, 70, 80, Inf))
coxph(Surv(time, status) ~ age + sex_f + wt.loss + strata(kstrata), data = lung)

# 補救四（換掉效果量）的程式碼在 B3-07，那一頁是 RMST 的完整說明

import statsmodels.api as sm
from lifelines import CoxPHFitter

lung = sm.datasets.get_rdataset("cancer", "survival").data
d = lung[["time", "status", "age", "sex", "ph.karno", "wt.loss"]].dropna().copy()
d["event"]  = (d["status"] == 2).astype(int)
d["female"] = (d["sex"] == 2).astype(int)
d = d.drop(columns=["status", "sex"]).rename(columns={"ph.karno": "ph_karno",
                                                      "wt.loss": "wt_loss"})

cph = CoxPHFitter().fit(d, duration_col="time", event_col="event")

# 尺度化 Schoenfeld 檢定 + 殘差圖，並印出補救建議
cph.check_assumptions(d, p_value_threshold=0.05, show_plots=True)

# 分層版本：strata 內各有自己的 baseline hazard
d["k_strata"] = (d["ph_karno"] <= 70).map({True: "low", False: "high"})
CoxPHFitter().fit(d, duration_col="time", event_col="event",
                  strata=["k_strata"]).print_summary()

cox.zph() 的表怎麼讀

拿 survival::lung 的四變項模型跑出來——資料集 228 人，完整個案分析後進入模型的是 213 人、151 個事件：

cox.zph() 對這個模型的結果：

這裡的虛無假設是「比例風險成立」，所以 p 小才是壞消息。 Karnofsky 分數的 p = 0.007，明確違反；整體檢定（GLOBAL）p = 0.038，也達到顯著。性別的 p = 0.080 落在邊界，年齡與體重減輕沒有問題。

log-log plot：類別變項的圖形版本

另一個古典檢查是把 $\log(-\log \hat S(t))$ 對 $\log t$ 作圖。它的邏輯是：如果兩組符合比例風險（風險比為 $\theta$），那麼

也就是兩條線垂直相差一個常數——平行。曲線靠攏、發散或交叉，就是假設有問題。

違反了怎麼辦

先問一個問題：這個違反是發生在你關心的變項，還是只是某個共變項？

如果只是共變項（例如你要估治療效果，違反的是年齡），那麼把它分層掉就行，代價很小。如果違反的是主要暴露或治療變項，那你就不能再報單一 HR 了——因為那個數字沒有清楚的意義。以下四條路各有代價：

路一：分層 Cox（stratified Cox）

把違反的變項從「共變項」改成「層」（strata）：每一層有自己的基線風險 $h_0(t)$，不再要求層與層之間成比例。

把 Karnofsky 分成三層之後，cox.zph() 的 GLOBAL 從 p = 0.038 變成 p = 0.441——假設救回來了。性別的 HR 變成 0.572（0.405–0.809）。

代價：分層的那個變項不再有 HR。 你買回了假設，但失去了對它的估計。所以分層適合用在「這只是個需要控制的因素，我不在乎它的效果量」的變項上。

路二：把時間切開，各段各估一個 HR

用 survSplit() 在第 180 天把每個人的追蹤時間切成兩段，讓 Karnofsky 的係數在前後兩段各自估：

故事變得清楚了：Karnofsky 分數的保護效果集中在前 180 天（每高一分 HR 0.966，區間不跨 1），180 天之後未偵測到關聯（HR 0.998，區間跨過 1）。這比原本那個單一的 HR 0.988 有資訊得多——原本那個數字其實是這兩段的加權平均。

代價：切點是你選的。 先看資料再挑一個讓 p 值好看的切點，就是在做 data dredging。切點應該由臨床理由決定（例如「術後 30 天」「療程結束」），而且要在論文裡講清楚是怎麼決定的。

路三：讓係數變成時間的函數

用 tt() 直接假設 $\beta(t) = \beta_0 + \beta_1 \log t$，把 $\beta_1$ 估出來。這裡的 $\beta_1$ = 0.011（SE 0.006，p = 0.062）——正號代表係數隨時間往上跑，與 Schoenfeld 圖看到的方向一致。

代價：$\log t$ 這個形式也是你假設的。 換成 $t$、$\sqrt t$ 會得到不同的結果。而且輸出變成兩個係數，臨床讀者不容易解讀。

路四：換一個不需要比例風險假設的效果量

前面三條路都還在修 Cox 模型。第四條路是換掉效果量本身：限制平均存活時間（restricted mean survival time, RMST）是 KM 曲線在 0 到 $\tau$ 之間的面積——「在前 $\tau$ 這段期間，平均活了多久」。它不需要任何比例假設，單位是「天」，而且不必要求兩組的風險比在整段追蹤裡維持固定。

RMST 現在不只是 PH 違反時的補救，在腫瘤與心衰竭試驗裡它已經是常規的次要分析。它的完整說明——面積怎麼來的、$\tau$ 為什麼必須事先指定、換一個 $\tau$ 結論會變多少、以及怎麼把它講成「平均多活幾個月」——在 限制平均存活時間（RMST）。

怎麼讀報表

論文裡跟這一頁有關的訊息通常只有一句話，藏在 Statistical analysis 段落。要找的是三件事：

1. 有沒有提到檢查過。 完全沒提，而主要結果是 HR，就是一個實質的方法學缺口。寫法通常是「The proportional hazards assumption was assessed using scaled Schoenfeld residuals」。
2. 檢查結果怎麼處理。 「假設成立」很好；「假設違反，因此改採分層／時間分段／RMST」也很好；「假設違反但我們仍報告單一 HR」則要自己在心裡打折。
3. KM 圖有沒有交叉。 這是不必看 Methods 就能做的檢查。兩條曲線在追蹤中段明確交叉、而且交叉之後持續分開，是比例風險不成立的強烈訊號——真實的風險比若固定在某個不等於 1 的值，兩條真實的存活曲線不會交叉。但要記得你看到的是估計出來的曲線：尾端只剩下少數人的時候，兩條 KM 曲線光靠抽樣波動就可能互相穿過，這種交叉不構成違反的證據。所以要對照風險人數表看交叉發生在哪一段、幅度多大，再回到 Methods 找 Schoenfeld 殘差檢定或時間交互作用的結果來確認。確認之後，那個單一 HR 就是兩段相反效果的平均，log-rank 的檢定力也會被削弱。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B3-04-ph-assumption.R

延伸觀看

▶The Cox proportional hazards model explained

▶Survival Analysis Part 9 | Cox Proportional Hazards Model

▶Cox Proportional Hazard Models

▶【Lecture】L20 Survival Analysis (2)

回報這一頁的問題