競爭風險：CIF 與 Fine-Gray

這個方法在解決什麼問題

survival::mgus2 是 Mayo Clinic 的 MGUS（monoclonal gammopathy of undetermined significance，意義未明的單株免疫球蛋白病變）世代，1384 位病人，診斷時年齡中位數 72 歲，最長追蹤 35 年。

臨床問題是：MGUS 病人有多少比例最後會進展成漿細胞惡性腫瘤（plasma cell malignancy，多發性骨髓瘤那一類）？ 這個數字直接決定要不要追蹤、多久追一次。

資料的實際結局分布是：

看到問題了嗎：死於他因的人數（860）是進展人數（115）的 7.5 倍。 這個世代的病人診斷時年紀很大，多數人在 MGUS 有機會進展之前，就先死於心血管疾病、其他癌症或其他原因了。

一個已經死於心肌梗塞的病人，不可能再進展為骨髓瘤。這叫做競爭事件（competing event）：它不是把你的觀察打斷，而是永久排除了目標事件發生的可能。這就是設限與競爭事件的根本差別。

為什麼 1 − KM 會高估

Kaplan-Meier 的整套邏輯建立在非資訊性設限上：被設限的人「還在，只是我們看不到了」，他們未來仍然可能發生事件，而且發生率與那些留下來的人一樣。

把死亡當成設限來跑 KM，等於在假設：那些死掉的人，如果我們繼續追蹤，會以跟活人一樣的速率進展成骨髓瘤。 這是不可能的。

⚠️ 要注意這個假設是 1 − KM 這個估計量才需要的，不是「設限」這個動作本身需要的。cause-specific hazard 依定義就條件在「此刻還在風險集合裡的人」，它從來沒有要對死者的未來說什麼——所以把競爭事件設限來估 cause-specific hazard 是合法的。本頁後面談兩種 hazard 時會回到這件事。

結果是 1 − KM 系統性地把死者的「假想進展」算了進去。這不是小偏差：

在 20 年這個時點，正確答案是 10.0%，而 1 − KM 說 21.0%——超過兩倍。到了 30 年更誇張（3.17 倍），因為那時風險集合只剩 3 個人。

「死於他因」那一欄才是高估的來源。 同樣在 20 年，有 72.4% 的人已經死於其他原因——他們不可能再進展成惡性疾病了。1 − KM 把這些人當成「還在等著進展、只是我們沒看到」，於是把剩下那些真正還有風險的人的比例，錯誤地放大成整群人的比例。

把三個部分加起來就看得出這件事：進展 10.0% ＋ 死於他因 72.4% ＋ 仍然存活且未進展 17.6% ＝ 100%。每一個人都必然落在這三格之一，這是 CIF 的定義帶來的性質，而 1 − KM 沒有這個性質。

CIF：累積發生函數

累積發生函數（cumulative incidence function, CIF）回答的是一個直接的問題：到時間 $t$ 為止，實際上有多少比例的人已經因為原因 1 發生了事件？

拆開來看：$h_1(u)$ 是在 $u$ 這一刻因為原因 1 發生事件的瞬時速率，$S(u^-)$ 是能撐到 $u$ 這一刻、任何事件都還沒發生的機率。兩者相乘再沿時間累加。

關鍵在於 $S(u^-)$：它是整體的無事件存活機率，被所有原因共同壓低。所以競爭事件多的時候，$S(u^-)$ 掉得快，能夠進展的人越來越少，$F_1(t)$ 自然爬得慢。1 − KM 之所以錯，就是它在計算時把 $S(u^-)$ 換成了「只考慮原因 1」的存活機率——一個假裝沒有人會死於他因的世界。

組間比較：分組的 CIF 與 Gray’s test

前面兩節的圖都只有一條曲線。但論文裡的 Figure 1 幾乎從來不是一條曲線——它是分組的 CIF 加上一個 p 值。讀者到這裡的處境通常是：已經知道 1 − KM 不能拿來當累積發生率，卻不知道那兩條 CIF 曲線旁邊的 p 值該從哪裡來——log-rank 檢定的是另一件事。這一節補的就是那格空白。

把同一份 mgus2 依診斷時的年齡切成兩組：

兩個檢定問的不是同一件事

Gray’s test 是 CIF 的兩組檢定：它比較的是 subdistribution hazard，也就是讀者眼睛正在看的那兩條曲線。log-rank 比較的則是 cause-specific hazard——它把競爭事件當成設限移出風險集合，而那正是取得 cause-specific hazard 的標準做法（本頁下一節說的「跑一個普通的 Cox 模型，把競爭事件當成設限」用的是同一招）。

所以這不是一個有效、一個壞掉。本頁第二節反對的是另一件事：把設限之後算出來的 1 − KM 當成累積發生率。設限這個動作本身沒有問題，有問題的是拿設限之後的存活曲線去回答「最後有多少人出事」。

同一份資料、同一個目標事件、同一組分組，兩個檢定回答兩個問題：

在 α = 0.05 之下，兩者的結論並不一致：Gray’s test 偵測到兩組 CIF 的差異，cause-specific log-rank 則未偵測到差異。

CIF 與 1 − KM 會排出相反的順序

兩個檢定不一致有解釋——它們問的不是同一件事。真正該停下來的是另一件事：同一個問題被兩種算法排出相反的順序——哪一組比較高。

這個相反是有時點的，不是整條曲線的性質。 把表格依時點讀一遍就看得出來：

• 第 60 個月：兩種算法都把年長組排在上面。
• 第 180 個月：兩種算法又都把年輕組排在上面。
• 第 240 個月：正確的 CIF 說未滿 70 歲比較高（13.4% 對 7.8%），而 1 − KM 說70 歲以上比較高（27.8% 對 20.3%）。

所以不要把這件事記成「兩條曲線整段方向相反」。曲線在追蹤中途交叉過，相反的是第 240 個月那一格的排序——而那正好是很多論文選來報「20 年累積發生率」的那一格。

原因攤在同一份資料裡。同樣到第 240 個月，死於他因的 CIF 是 51.7%（未滿 70 歲）對 87.6%（70 歲以上），Gray’s test p < 0.001。年長組不是比較不會進展，是多數人先死於別的原因；1 − KM 把那些人記成「還在等著進展、只是我們沒看到」，於是把年長組的風險灌大到超過年輕組。

尾端那一段陡直的跳升，是個位數的人畫出來的

本頁前面已經說過，1 − KM 曲線「尾端的劇烈跳動是因為那時風險集合只剩個位數」。分組之後，那句話有了具體的數字：70 歲以上到第 240 個月只剩 4 個人還在風險集合裡（未滿 70 歲是 53 個）。

圖上紅色虛線在接近右緣處那一段近乎垂直的跳升，就是這幾個人裡有人進展造成的：那條線從第 180 個月的 13.3% 一路跳到第 240 個月的 27.8%，高估倍數也從 1.93 倍變成 3.55 倍。分母只剩個位數的時候，一個事件就能把曲線推上十幾個百分點。

這與本頁「怎麼讀報表」第 5 點、以及「常見誤用」最後幾列講的是同一件事：沒有風險人數表，你沒有辦法判斷尾端那一段是幾個人畫的。 只是在競爭風險的情境下，它還多一層後果——被推上去的那條線，正好是排序反轉的那一條。

堆疊 CIF：三個狀態填滿到 100%

本頁在「為什麼 1 − KM 會高估」那一節用文字論證過一件事：進展、死於他因、仍存活且未進展，三格加起來剛好是 100%。那是 CIF 的定義帶來的性質，也是 1 − KM 沒有的性質。當時只有文字，現在把它畫出來。

這不是示意圖。腳本在 538 個實際發生轉移的時點上檢查「三者相加是否等於 1」，最大偏差是 1.1e-15——那是浮點運算的捨入誤差量級，不是四捨五入湊出來的。（檢查跑在真實的轉移時點上，不是跑在匯出給這張圖的那 21 個時間網格點上；只查網格點會讓檢查變寬鬆，因為兩個轉移之間本來就不會出問題。）

而兩格擺在一起，上一節那個「排序反轉」就不需要再解釋了：

• 70 歲以上那一格，淺褐色到第 240 個月已經佔掉 87.6%，只剩 4.6% 的人還活著而且沒進展。深藍色不是被什麼力量壓下去的——是沒有剩下多少人可以進展。
• 未滿 70 歲那一格，綠色到第 240 個月還有 35.0%。藍色因此還有空間繼續往上長，最後長到 13.4%。

綠色那一塊剩多少，決定了藍色還能長多高。 這就是競爭風險的全部機制，畫成圖之後不需要公式。

也因為這樣，這張圖是分辨兩種算法最快的方法：能填滿到 1 的是 CIF，填不滿或滿出來的是 1 − KM。 本頁前面那個 Callout 說「各結局的 1 − KM 加起來可以超過 100%」，講的就是同一件事；「怎麼讀報表」第 2 點給的判準，操作起來就是這張圖。

兩種 hazard，兩個不同的問題

這是競爭風險最常被誤用的地方，也是這一頁真正的重點。有兩種模型，它們不是「一個對一個錯」，而是回答不同的問題。

Cause-specific hazard（因別風險）

分母的條件是 $T \ge t$：還沒發生任何事件的人。做法就是跑一個普通的 Cox 模型，把競爭事件當成設限。

它回答的是：在還沒出事的人裡面，這個變項怎麼影響「發生目標事件」的速率？ 這是關於生物機轉的問題——這個因子是不是真的在推動疾病進展。

Subdistribution hazard（次分布風險，Fine-Gray）

分母多了一塊：已經發生競爭事件的人仍留在風險集合裡（技術上是用逐漸下降的權重留著）。這聽起來很怪——一個已經死掉的人怎麼還在風險集合裡？但它的目的很明確：讓模型的係數直接對應到 CIF 的高低。Fine-Gray 模型的 sHR 大於 1，就代表那個變項會讓 CIF 曲線爬得更高。

它回答的是：這個變項怎麼影響「最後實際有多少人發生目標事件」？ 這是關於負擔與預後的問題——要不要追蹤、要追多久、資源怎麼配。

在 MGUS 資料上，兩者會給相反的答案

年齡的兩個估計指向相反方向，而且兩邊都是對的。

• Cause-specific HR = 1.016（1.000–1.033，p = 0.051）：在還沒出事的人裡，年紀大的人進展速率略高——不過區間貼著 1，本研究未達統計顯著。
• Fine-Gray sHR = 0.983（0.973–0.994，p = 0.003）：年紀越大，最後實際進展的比例越低，而且區間完全在 1 的左邊。

矛盾嗎？一點也不。看第三個模型就懂了——同一份資料、同樣的變項，但結局改成「死於他因」：

年齡對死於他因的 HR 是 1.067（p < 0.001）——每多一歲，死於他因的速率上升約 6.7%。 老年病人不是比較不會進展，而是還來不及進展就先走了。生物層次上年齡對進展沒有保護作用（甚至略往反方向），但在「最後有多少人真的進展」這個層次上，年齡確實把 CIF 壓低了。

相對地，M-spike 的兩個估計幾乎一樣（2.421 vs 2.431），因為它對死於他因幾乎沒有作用（HR 0.942，p = 0.353，未偵測到關聯）。兩種 HR 的差距，來自那個變項對競爭事件的影響有多大。

動手跑一次

library(survival)
data(cancer, package = "survival")   # mgus2 在這包資料裡

# 1. 把兩個時間軸併成一個，結局做成 factor
#    第一個 level 必須是「設限」，survival 是這樣讀 factor endpoint 的
mg <- mgus2
mg$etime <- with(mg, ifelse(pstat == 0, futime, ptime))
mg$event <- with(mg, ifelse(pstat == 0, 2 * death, 1))
mg$event <- factor(mg$event, 0:2, c("censored", "progression", "death"))
table(mg$event)

# 2. CIF：Surv() 吃 factor 就自動走多狀態（Aalen-Johansen）
cif <- survfit(Surv(etime, event) ~ 1, data = mg)
summary(cif, times = c(60, 120, 240, 360))   # pstate 欄就是各原因的 CIF
plot(cif, col = c("#4d6a8c", "#8c6d4a"), lwd = 2, noplot = "(s0)",
     xlab = "Months", ylab = "Cumulative incidence")

# 3. 錯誤做法的對照：把死亡當設限之後，再把 1 - KM 當成累積發生率
naive <- survfit(Surv(etime, pstat) ~ 1, data = mg)
1 - summary(naive, times = 240)$surv          # 明顯高於上面的 CIF

# 4. Cause-specific Cox：競爭事件當設限，就是普通的 coxph
coxph(Surv(etime, event == "progression") ~ age + sex + mspike, data = mg)
coxph(Surv(etime, event == "death")       ~ age + sex + mspike, data = mg)

# 5. Fine-Gray：finegray() 產生加權長表，再跑加權 Cox
fgd <- finegray(Surv(etime, event) ~ ., data = mg, etype = "progression")
coxph(Surv(fgstart, fgstop, fgstatus) ~ age + sex + mspike,
      weight = fgwt, data = fgd)

# ── 組間比較 ────────────────────────────────────────────────────────────────
mg$agegrp <- factor(ifelse(mg$age >= 70, "Age >= 70", "Age < 70"),
                    levels = c("Age < 70", "Age >= 70"))
mg$ecode <- as.integer(mg$event) - 1L        # 0 設限, 1 進展, 2 死於他因

# 6. Gray's test：CIF 的兩組檢定。它住在 cmprsk 不在 survival
#    （survival 只內建 Fine-Gray 的 finegray()，沒有 Gray's test）。
#    這裡用 cmprsk::cuminc() 的 namespace 呼叫，不用 library(cmprsk)：
#    附加套件會把它排在 survival 前面，一旦 mask 掉同名函式，上面那三個
#    已經發表在頁面上的模型可能悄悄換成別的數字。:: 呼叫 mask 不到東西。
gray <- cmprsk::cuminc(ftime = mg$etime, fstatus = mg$ecode,
                       group = as.character(mg$agegrp), cencode = 0)
gray$Tests                                   # 每個 cause 一列：stat / pv / df

# cmprsk 自己算的 pv 是 1 - pchisq(stat, df)，統計量夠大時會 underflow 成
# 剛好的 0。頁面上印出 "p = 0" 是讀者看不出來的錯，所以改算上尾。
pchisq(gray$Tests[, "stat"], gray$Tests[, "df"], lower.tail = FALSE)

# 7. Cause-specific 的組間檢定：把競爭死亡當設限，跑普通 log-rank。設限是
#    取得 cause-specific hazard 的標準做法，所以這是有效的檢定——它只是問了
#    機轉層次的問題，與第 6 步的 Gray's test 不是同一件事。
survdiff(Surv(etime, event == "progression") ~ agegrp, data = mg)

# 8. 分組 CIF 與堆疊圖的資料來源是同一個物件。factor endpoint 的 survfit()
#    回傳 pstate：一個狀態一欄，而且每一列加起來恰好是 1——堆疊圖畫的就是
#    這個性質。
cif_g <- survfit(Surv(etime, event) ~ agegrp, data = mg)
sg <- summary(cif_g, times = c(60, 120, 180, 240))
sg$pstate                                    # (s0) / progression / death
rowSums(sg$pstate)                           # 每一列都是 1

import pandas as pd
from lifelines import AalenJohansenFitter, KaplanMeierFitter

# statsmodels 的 get_rdataset 查不到它（Rdatasets 把它登記在別的 bundle 底下），
# 所以直接讀 CSV。
RD = "https://vincentarelbundock.github.io/Rdatasets/csv/"
mg = pd.read_csv(RD + "survival/mgus2.csv")
mg["etime"] = mg["futime"].where(mg["pstat"] == 0, mg["ptime"])
mg["event"] = (2 * mg["death"]).where(mg["pstat"] == 0, 1)   # 0 censor,1 prog,2 death

ajf = AalenJohansenFitter()
ajf.fit(mg["etime"], mg["event"], event_of_interest=1)
print(ajf.cumulative_density_.loc[:240].tail(1))    # CIF at 240 months

# 錯誤做法的對照
kmf = KaplanMeierFitter().fit(mg["etime"], (mg["event"] == 1).astype(int))
print(1 - kmf.predict(240))

# Fine-Gray 在 Python 沒有成熟實作；需要 sHR 時請回到 R 的 survival::finegray()

# 分組 CIF：對每一組各配一次 AalenJohansenFitter
mg["agegrp"] = (mg["age"] >= 70).map({True: "Age >= 70", False: "Age < 70"})
for g, sub in mg.groupby("agegrp"):
    ajf = AalenJohansenFitter().fit(sub["etime"], sub["event"], event_of_interest=1)
    print(g, float(ajf.cumulative_density_.loc[:240].iloc[-1]))

# 組間檢定的部分同樣沒有：lifelines 有 statistics.logrank_test，但那是
# cause-specific（要把競爭事件當設限），不是 Gray's test。Python 生態裡
# 沒有與 cmprsk::cuminc() 同等成熟的 CIF 組間檢定，要 p 值請回到 R。

怎麼讀報表

看到一篇有競爭風險的論文（老年腫瘤、移植、器材登錄、失智世代都算），依序檢查：

1. 有沒有承認競爭事件存在。 主要結局是「疾病進展」「翻修手術」「特定死因」而族群年紀大，報表上卻只有 1 − KM 曲線與一個 log-rank p 值，那份分析就只做了機轉那一半：累積發生率被高估，而「最後有多少人出事」根本沒被檢定過。
2. 累積發生率是 CIF 還是 1 − KM。 圖的 y 軸寫 “cumulative incidence” 不保證是 CIF——有些論文把 1 − KM 也這樣標。判斷方式：如果各結局的曲線在最後加起來明顯超過 1，那就是 1 − KM。
3. HR 是哪一種，有沒有講清楚。 報表寫 “HR” 通常是 cause-specific，寫 “sHR” 或 “subdistribution hazard ratio” 才是 Fine-Gray。沒寫的話只能從 Methods 裡的函式名（coxph vs crr／finegray）推。
4. 解讀的語言與模型對不對得上。 用 Fine-Gray 的 sHR 講「這個因子促進疾病進展」是錯配——sHR 講的是最終比例，不是機轉。
5. 風險人數表。 與 KM 一樣重要，理由也一樣：尾端的曲線是幾個人畫的。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B3-05-competing-risk.R

延伸觀看

▶COMPETING RISK EXPLAINED

▶Competing Risk Analysis

▶Easily Perform Competing Risks Survival Analysis with SAS Studio

回報這一頁的問題