時間相依共變項

兩件常被混為一談的事

中文都叫「時間相依」，但它們是完全不同的兩個問題：

(start, stop] 長表：一個人可以有很多列

一般的存活資料一人一列：一個時間、一個事件指標、若干個固定共變項。暴露會變動的時候這個形狀就不夠了，因為「這個人的 $x$ 是多少」在不同時刻有不同答案。

解法是把每個人切成若干段，每一段內部 $x$ 是常數。每一段一列，欄位是 (tstart, tstop] 這個左開右閉的區間、該段結束時有沒有發生事件、以及該段適用的共變項值。

關鍵在於，這個結構沒有把一個人變成很多人。Cox 模型的估計完全發生在「每一個事件時點的風險集合」上（見B3-01）：在時間 $t$ 這一刻，模型只問「此刻還在追蹤的這些人，各自的 $x$ 現在是多少」。長表的每一列只是回答這個問題的查詢表；一個人在時間 $t$ 只會有一列涵蓋那一刻，所以他只會被算進去一次。

tmerge()：把寬表變成長表

手動切區間很容易出錯（重疊、漏段、事件記在錯的那一列），survival 套件的 tmerge() 就是專門做這件事的。它的心智模型是：先建一個只有追蹤區間與結果的骨架，再一層一層把事件與時間相依變項疊上去。

三個主要動詞：

用史丹佛換心資料（103 位等待名單上的病人，其中 69 位接受了移植）跑一次，會得到 170 列：34 位從未移植的人各一列，67 位移植的人各兩列，剩下 2 位是在第 0 天就移植的——tdc() 在時間零切不出區間，所以他們也只有一列，而且那一列的 transplant 已經是 1。

34 + 2 = 36 位只有一列的人，加上 67 位有兩列的人，正好是 103 位病人；而 67 + 2 = 69 位接受移植。建好長表就把這兩條加法算一次——本頁下面的「常見誤用」表說長表建錯不會報錯、只會靜靜給出錯的答案，這兩位第 0 天移植的病人就是那句話的活教材。

重複測量：基線值還是當下值

pbcseq 是原發性膽汁性膽管炎的追蹤資料：312 位病人、定期回診驗血，中位回診 5 次、最多 16 次，最長追蹤 12.5 年，期間 125 人死亡。轉成長表後有 1807 列。

膽紅素是 PBC 最重要的預後指標，而它會隨疾病進展上升。用基線那一次的值，跟用「當下最近一次測到的值」，跑出來差很多：

基線膽紅素的 HR 比較小，理由不神祕：基線值是當下值的一個逐漸失效的代理。 追蹤越久，基線那一次抽血離「此刻的病情」越遠，於是關聯被稀釋（regression dilution）。用當下值就把那層稀釋拿掉了。

重複事件：Andersen-Gill 與穩健變異數

有些結果會反覆發生：感染、氣喘發作、心衰竭再住院、癲癇。一人一列的做法只能留下第一次，其餘全部丟掉。

cgd 是慢性肉芽腫病（chronic granulomatous disease）的隨機試驗，比較干擾素 gamma 與安慰劑對嚴重感染的預防效果。128 位孩子，追蹤期間共發生 76 次感染，其中 17 位發生不只一次，最多的一位發生了 7 次。

只取第一次事件，等於丟掉 32 次事件，佔全部的 42.1%。

Andersen-Gill（AG）模型的做法非常樸素：把每個人的追蹤切成「上一次事件之後到下一次事件」的區間，全部丟進同一個 Cox 模型，就當作是不同的觀察。資料結構與上一節的長表完全一樣。

Landmark 分析：更笨但更容易講清楚的替代方案

時間相依模型不是唯一的解。landmark 分析（landmark analysis）的做法是：挑一個時點 $L$，只保留活過 $L$ 的人，依他們在 $L$ 那一刻的狀態分組，時間從 $L$ 重新起算。

它為什麼避開了 immortal time bias：分組所需的資訊在 $L$ 那一刻全部已知，$L$ 之前的時間對兩組一視同仁地被排除。

用史丹佛換心資料，在第 30 天設 landmark：

三個 landmark 的結果都未達統計顯著（信賴區間都跨過 1），與 B3-01 那個時間相依模型的結論一致——一致的是「都未達統計顯著」這一點，不是點估計的方向（B3-01 的 HR 大於 1，這裡三個裡有兩個小於 1）。但點估計從第 14 天的 1.441 移到第 60 天的 0.808，而且保留的人數從 88 掉到 64。

解讀上的限制：時間相依的 HR 不能拿來預測

這是本頁最重要、也最容易被跳過的一節。

一個時間相依共變項的 HR 是同時期的關聯（contemporaneous association）：它說的是「在時間 $t$ 這一刻，$x(t)$ 高的人在該刻的風險較高」。它不是一個預測陳述。

理由是：要用這個模型算「某人未來五年的存活機率」，你必須知道他未來五年每一刻的 $x(t)$——而那正是你沒有的東西。膽紅素會怎麼變，本身就是未知的；如果你知道，你多半也已經知道結局了。

動手跑一次

library(survival)

# ── 1. 長表長什麼樣 ──
data(cgd, package = "survival")
head(cgd[cgd$id %in% c(1, 2), c("id", "tstart", "tstop", "status", "treat", "enum")])

# ── 2. tmerge()：從寬表建長表 ──
# 骨架：一人一列的追蹤區間 + 結果
base <- subset(pbc, id <= 312, select = c(id:sex, stage))
d <- tmerge(base, base, id = id, death = event(time, status == 2))
# 疊上時間相依的檢驗值；tdc() 預設「該時刻之後」才生效
d <- tmerge(d, pbcseq, id = id,
            bili = tdc(day, bili), albumin = tdc(day, albumin))
head(d[d$id == 2, c("id", "tstart", "tstop", "death", "bili", "albumin")])

# ── 3. 基線值 vs 當下值 ──
b <- merge(base, pbcseq[!duplicated(pbcseq$id), c("id", "bili", "albumin")], by = "id")
coxph(Surv(time, status == 2) ~ log(bili) + log(albumin) + age, data = b)
coxph(Surv(tstart, tstop, death) ~ log(bili) + log(albumin) + age, data = d)

# ── 4. 重複事件：Andersen-Gill ──
coxph(Surv(tstart, tstop, status) ~ treat, data = cgd)                  # 標準誤是錯的
coxph(Surv(tstart, tstop, status) ~ treat + cluster(id), data = cgd)    # 穩健變異數
# PWP：依事件次序分層，第 k 次事件有自己的基線風險
coxph(Surv(tstart, tstop, status) ~ treat + strata(enum) + cluster(id), data = cgd)

# ── 5. landmark 分析 ──
data(heart, package = "survival")
L <- 30
lm <- subset(jasa, futime >= L)
lm$g <- factor(ifelse(!is.na(lm$wait.time) & lm$wait.time <= L, "tx", "no"))
coxph(Surv(futime - L, fustat) ~ g, data = lm)
plot(survfit(Surv(futime - L, fustat) ~ g, data = lm))

# ── 6. 檢查資料整理有沒有出錯（一定要做）──
# 區間不可重疊、不可有 tstop <= tstart、每個人的區間要接得起來
with(cgd, tapply(seq_along(id), id, function(i)
  all(tstart[i][-1] == tstop[i][-length(i)]))) |> table()

import numpy as np
import pandas as pd
from lifelines import CoxTimeVaryingFitter

# 長表的來源：pbcseq 的逐次檢驗值，接上 pbc 的追蹤終點。
RD  = "https://vincentarelbundock.github.io/Rdatasets/csv/"
pbc = pd.read_csv(RD + "survival/pbc.csv")
seq = pd.read_csv(RD + "survival/pbcseq.csv")
base = (pbc.loc[pbc["id"] <= 312, ["id", "time", "status", "age"]]
           .rename(columns={"time": "futime"}))
measurements = (seq[["id", "day", "bili"]].merge(base, on="id")
                   .assign(log_bili=lambda d: np.log(d["bili"])))

# 長表欄位：id / start / stop / event / 共變項；一個人多列
# lifelines 沒有 tmerge()，區間要自己切（pandas 的 merge_asof 很好用）
long = (measurements
        .sort_values(["id", "day"])
        .assign(start=lambda d: d["day"],
                stop=lambda d: d.groupby("id")["day"].shift(-1)))
long["stop"] = long["stop"].fillna(long["futime"])
long["event"] = ((long["stop"] == long["futime"]) & (long["status"] == 2)).astype(int)
long = long[long["stop"] > long["start"]]          # 空區間會讓模型報錯

ctv = CoxTimeVaryingFitter()
ctv.fit(long[["id", "start", "stop", "event", "log_bili", "age"]],
        id_col="id", event_col="event", start_col="start", stop_col="stop")
ctv.print_summary()

# 重複事件：R 那邊是 coxph(..., cluster = id)。lifelines 這一側做不到——
# CoxTimeVaryingFitter 的 robust=True 目前直接丟 NotImplementedError，
# 而不做 robust 會低估標準誤。這一段要回 R，或自己算三明治估計。

怎麼讀報表

1. 暴露是什麼時候決定的。 只要分組依據要等到基線之後才知道（接受手術、有無反應、完成療程、開始某個藥），就去 Methods 找他們用了時間相依變項、landmark、還是什麼都沒做。
2. 時間相依的共變項是什麼時候測的、怎麼帶進模型。 「最近一次的值」「該區間的平均」「lag 一週的值」是三種不同的定義，結果會不同。都沒寫，就是不可複製。
3. 重複事件有沒有用穩健變異數。 找 “robust”、“sandwich”、“cluster”、“covs(aggregate)” 這些字。沒有，信賴區間就要自己在心裡放寬約兩成。
4. 重複事件用的是哪一種模型。 AG、PWP、WLW、frailty 的假設不同，只寫「Cox model for recurrent events」不夠。
5. landmark 的 $L$ 是怎麼決定的、丟掉了多少人。 有沒有寫在計畫書裡；被排除的那些人（死於 $L$ 之前）的特徵有沒有交代。
6. 有沒有拿時間相依模型去做預測。 出現「本模型的 C 統計量為 X，優於基線模型」而模型含時間相依共變項，那個比較不公平。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B3-06-time-varying.R

延伸觀看

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

▶【Hands-on】L12 R：Survival Analysis

▶The Statistics of Life and Death | Survival Analysis

回報這一頁的問題