邊際結構模型與時間相依 IPTW

為什麼這個方法值得有自己的一頁

你在論文裡撞到這一頁的時機大概是三種：Methods 寫著 patients who switched treatment、 寫著 time-updated covariates、或者寫著 per-protocol analysis 而且不只是把不服藥的人刪掉。 自己手上那份資料撞到的時機只有一種：病人會中途開始用藥、停藥、換藥，而不是在第一天就被分好組、 從此不動。

這種資料有一個標準迴歸處理不了的形狀。決定病人第二次要不要治療的那個指標—— CD4 數值、血壓、腎功能、疾病活動度——本身就是第一次治療造成的。它同時是後續治療的干擾因子， 又是先前治療的後果。標準迴歸對這種變項只有兩個選項：放進去，或不放。 這一頁要說明的是為什麼兩個都錯，以及第三條路長什麼樣。

站上有三個地方把讀者送到這裡：時間相依共變項那頁講到 「時間相依的中介兼混淆不能用標準 Cox 處理」、IPTW 那頁講到 「暴露會隨時間改變時只有時間相依的加權走得通」、目標試驗模擬 那頁講到「per-protocol 要處理隨時間改變的暴露」。這一頁是那三條線的終點， 所以它自己不再把你轉介回去——底下有完整的權重形式、可以照跑的程式碼、 以及一份真值已知的模擬拿來對答案。

治療與干擾互相餵養：feedback 是什麼

把時間軸縮到最短——兩次門診——就足以把問題整個攤開：

• visit 0：病人接受或不接受治療，記作 $A_0$。在這份模擬裡它是隨機分派的。
• visit 1：先量到一個臨床指標 $L_1$，再根據它決定第二次治療 $A_1$。
• 結局 $Y$：例如一年內的某個事件。

$L_1$ 站在一個尷尬的位置。它是 $A_0$ 的後果（治療改變了指標）， 又是 $A_1$ 的干擾因子（指標決定了第二次要不要治療，而指標本身也預測結局）。 這就是 treatment-confounder feedback：治療改變共變項，共變項再回頭決定下一次治療。

圖上還有一個沒被觀察到的 $U$，它同時造成 $L_1$ 與 $Y$。臨床上這種東西到處都是： 沒被記錄下來的疾病嚴重度、體能狀態、社經條件、家屬支持——它們影響指標，也影響結局， 而資料庫裡沒有欄位。

模擬用的生成過程把上面每一條箭頭都寫成一行式子：

實際生出來的世代長這樣：20,000 個人，事件率 13.4%。 $A_0$ 沒治療的人裡有 40.4% 的 $L_1$ 是陽性，有治療的則是 68.5%——這是 $A_0 \to L_1$ 那支箭頭在資料裡的樣子。 而 $L_1$ 陰性的人有 40.3% 接受第二次治療，陽性的則是 70.6%——這是 $L_1 \to A_1$。兩件事同時成立，feedback 就成立了。

「放進去」與「不放」，兩個都錯

不放進去也沒有比較好：$L_1$ 是 $A_1$ 真正的干擾因子（它同時預測 $A_1$ 與 $Y$， 後者經由 $U$），不處理就是 confounding by indication 原封不動留在估計裡。

三條路跑出來的合計 log odds ratio（$A_0$ 與 $A_1$ 兩期都治療對兩期都不治療）是這樣：

方向是可以事先推出來的，而資料也照著走：把 $L_1$ 放進迴歸偏過頭 （-1.427，比真值離虛無值更遠），完全不校正偏不足（-0.975，比真值更靠近零）， 加權落在中間而且最靠近真值（-1.173，差 0.074）。 三條路各偏一個方向，這正是這份模擬要示範的東西。

加權是第三條路：MSM 在估什麼

先把估計目標寫清楚。用潛在結果的記號，$Y^{a_0, a_1}$ 是「假如這個人在 visit 0 被指定 $a_0$、在 visit 1 被指定 $a_1$，他會不會發生事件」。四種靜態策略各對應一個全世代的風險：

邊際結構模型就是直接對這四個風險配一個模型，而不是對「觀察到的 Y 條件在共變項上」配模型：

「邊際」的意思是左邊沒有條件在任何共變項上——它講的是整個世代的風險，不是某一種病人的風險。 「結構」的意思是左邊是潛在結果，不是觀察到的結果。$\beta_1 + \beta_2$ 就是上面表裡 「兩期都治療」對「兩期都不治療」的 log odds ratio，也就是這一頁一路在追的那個數字。

那要怎麼在只有觀察資料的情況下估這個模型？把每個人依他實際接受的治療機率的倒數加權， 造出一個假想世代（pseudo-population）：在那個世代裡，每個時點的治療分派都跟當時的共變項無關， 就像有人在每次門診擲了一次硬幣。

穩定化權重長什麼樣

未穩定化的權重是每個時點治療機率倒數的連乘。它會爆炸：只要某個人在某個時點接受的是 「以他的狀況幾乎不會有人給」的治療，分母就趨近零，那一個人就可以主宰整個分析。

穩定化權重在分子放回一個只用過去治療、不用共變項的機率：

上面那條橫槓表示「到該時點為止的全部歷史」。分母負責把共變項造成的干擾拿掉， 分子負責把「純粹由過去治療決定的那一部分」放回去，於是權重會集中在 1 附近而不是散開。 分子裡不可以放共變項，放了就等於把要拿掉的東西又加回來。

這份模擬的兩個模型是：分子 A1 ~ A0，分母 A1 ~ A0 + L1。

權重的樣子：平均 1.000、標準差 0.284、範圍 0.653 到 1.517，有效樣本數（ESS）18,505， 是名目 n 的 92.5%。穩定化是它這麼乖的原因。 同一組權重的未穩定化版本——分子換成 1、不是邊際治療機率——落在 2.76 到 7.43 之間，平均 4.00， 最大值是穩定化版本的 4.90 倍、標準差是 5.20 倍。其中 2.00 倍只是記帳問題 （$A_0$ 在這裡是隨機分派，visit 0 的分子是個常數），但有效樣本數還是掉到 17,597，是名目 n 的 88.0%， 對比穩定化版本的 92.5%。

這張表值得多看兩眼，因為它把「加權到底在做什麼」變成看得見的東西： 權重大於 1 的四群，全部都是 $L_1$ 與 $A_1$ 不一致的人——指標陽性卻沒治療、 或指標陰性卻治療了，共 6,881 人。他們在真實世界裡是少數， 但在假想世代裡被放大到跟其他人一樣多，好讓「治療」與「指標」脫鉤。 權重小於 1 的四群則相反，是照著指標走的多數人，他們被縮小。

換句話說，加權是靠那些沒有照著臨床邏輯走的病人把干擾撬開的。 這也解釋了為什麼正性（positivity）是這個方法的命門：如果指標陽性的人沒有一個不治療， 那一格就是空的，權重會趨近無限大或者根本估不出來，而你在平均權重上看不見這件事—— 要看 ESS、看最大值、看每一個型態的人數，也就是上面這張表。

動手跑一次：長表、兩階段

實作分兩個階段，中間隔著一份權重。這是 MSM 跟一般迴歸最不一樣的地方： 你會配兩次模型，而第一次配的那個模型的係數，你一眼都不會看—— 它只是拿來算權重的中間產物。

# 產生寬表 d。U 是未測量的干擾——少了它，把 L1 條件化只會過度校正，
# naive 模型看起來就會是對的，整頁的論點會被自己的數字推翻。
set.seed(20260823)
n  <- 20000
U  <- rnorm(n)                                   # 未測量的干擾，非有不可
A0 <- rbinom(n, 1, 0.5)
L1 <- rbinom(n, 1, plogis(-0.5 + 1.5 * A0 + 1.2 * U))
A1 <- rbinom(n, 1, plogis(-0.5 + 1.2 * L1 + 0.3 * A0))
Y  <- rbinom(n, 1, plogis(-1.5 - 0.7 * A0 - 0.7 * A1 + 1.0 * U))
d  <- data.frame(A0, L1, A1, Y)

# ── 階段一：權重模型（長表，一人一個 visit 一列）──────────────────────
# d 是寬表：一人一列，欄位 A0 / L1 / A1 / Y。先攤成長表，
# 因為權重是「每個時點各算一次再連乘」，長表讓迴圈只寫一次。
n <- nrow(d)
long <- rbind(
  data.frame(id = 1:n, visit = 0, A = d$A0, L = 0,     lagA = 0),
  data.frame(id = 1:n, visit = 1, A = d$A1, L = d$L1,  lagA = d$A0)
)

K  <- 1                    # 最後一個 visit 的編號
sw <- rep(1, n)            # 每個人一個累積權重，從 1 開始連乘

for (k in 0:K) {
  dk <- long[long$visit == k, ]
  dk <- dk[order(dk$id), ]

  # 分子：只放「過去的治療」。分母：再加上「當下的共變項」。
  # 分子絕對不能放共變項——放了就等於把要拿掉的干擾又加回來。
  num <- glm(A ~ lagA,     data = dk, family = binomial())
  den <- glm(A ~ lagA + L, data = dk, family = binomial())

  pn <- predict(num, type = "response")   # P(A = 1 | 過去治療)
  pd <- predict(den, type = "response")   # P(A = 1 | 過去治療 + 共變項)

  # 每個人拿的是「他實際接受的那個治療」的機率，不是 P(A = 1)
  sw <- sw * ifelse(dk$A == 1, pn, 1 - pn) / ifelse(dk$A == 1, pd, 1 - pd)
}

# ⚠️ 在這份模擬裡 visit 0 的 lagA 與 L 都是常數，兩個模型都退化成截距模型，
# 那一期的權重因此恰好是 1（R 會把常數項的係數報成 NA，那是預期的）。
# 真實資料的 visit 0 有基線共變項 L0，分母是 A ~ L0，權重不會是 1。

summary(sw)                       # 平均要接近 1；離 1 很遠就是模型設錯了
sum(sw)^2 / sum(sw^2)             # 有效樣本數 ESS

# ── 階段二：加權的結果模型，就是 MSM ─────────────────────────────────
# 左邊是結局，右邊只有治療——沒有 L1。共變項只出現在階段一的分母裡。
# ⚠️ suppressWarnings 是刻意的，不是在蓋掉配壞的模型：加權後
# w_i * y_i 不再是整數，glm 因此丟 non-integer #successes in a binomial glm!
# 點估計是正確的加權 score 解，只有 glm 內部的二項式帳本不高興。
msm <- suppressWarnings(
  glm(Y ~ A0 + A1, data = d, family = binomial(), weights = sw)
)
coef(msm)["A0"] + coef(msm)["A1"]   # 兩期都治療 vs 兩期都不治療

# ⚠️ 不要讀 summary(msm) 的標準誤。理由見下一節。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
# 產生寬表 d。U 是未測量的干擾——少了它，把 L1 條件化只會過度校正，
# naive 模型看起來就會是對的，整頁的論點會被自己的數字推翻。

rng = np.random.default_rng(20260823)
n   = 20000
# numpy 的亂數流與 R 不同，同一個 seed 不會給出同一批人，
# 所以這一側的數字會接近但不等於上面那張表。
expit = lambda x: 1 / (1 + np.exp(-x))
U  = rng.standard_normal(n)
A0 = rng.binomial(1, 0.5, n)
L1 = rng.binomial(1, expit(-0.5 + 1.5 * A0 + 1.2 * U))
A1 = rng.binomial(1, expit(-0.5 + 1.2 * L1 + 0.3 * A0))
Y  = rng.binomial(1, expit(-1.5 - 0.7 * A0 - 0.7 * A1 + 1.0 * U))
d  = pd.DataFrame({"A0": A0, "L1": L1, "A1": A1, "Y": Y})

# d 是寬表 DataFrame，欄位 A0 / L1 / A1 / Y
n = len(d)
long = pd.concat([
    pd.DataFrame({"id": range(n), "visit": 0, "A": d.A0, "L": 0,    "lagA": 0}),
    pd.DataFrame({"id": range(n), "visit": 1, "A": d.A1, "L": d.L1, "lagA": d.A0}),
])

sw = np.ones(n)
for k in (0, 1):
    dk = long[long.visit == k].sort_values("id")
    # visit 0 的 lagA 與 L 是常數，patsy 會直接把它們吃掉；
    # 真實資料的 visit 0 要換成基線共變項。
    rhs_num = "lagA" if dk.lagA.nunique() > 1 else "1"
    rhs_den = " + ".join(c for c in ("lagA", "L") if dk[c].nunique() > 1) or "1"
    pn = smf.glm(f"A ~ {rhs_num}", dk, family=sm.families.Binomial()).fit().fittedvalues
    pd_ = smf.glm(f"A ~ {rhs_den}", dk, family=sm.families.Binomial()).fit().fittedvalues
    a = dk.A.to_numpy()
    sw = sw * np.where(a == 1, pn, 1 - pn) / np.where(a == 1, pd_, 1 - pd_)

# var_weights，不是 freq_weights：權重不是「這一列重複幾次」
msm = smf.glm("Y ~ A0 + A1", d, family=sm.families.Binomial(),
              var_weights=sw).fit()
print(msm.params["A0"] + msm.params["A1"])
# 這裡印出來的 bse 同樣不能用，要自己 bootstrap

真值不是那個結構係數：非塌縮性

生成 $Y$ 的那條式子裡，$A_0$ 與 $A_1$ 的係數各是 -0.70， 合計 -1.40。但這一頁一路在對的真值是 -1.247，不是 -1.40。 兩個數字不一樣，而且不是誰算錯。

差別在於它們是兩個不同的參數：

• -0.70 是條件參數——「在 $U$ 固定的前提下，治療讓 log odds 變多少」。
• -1.247 是邊際參數——「整個世代都治療 vs 整個世代都不治療，全世代風險的 log odds 差多少」。

odds ratio 是不可塌縮（non-collapsible）的：即使沒有任何干擾， 把一個有預測力的變項（這裡是 $U$）從模型裡拿掉，係數也會往零縮。 這是 odds ratio 這個尺度本身的性質，不是偏誤——風險差與相對風險就沒有這個問題。

MSM 不是魔法，它也是一個模型

上面那張表裡，MSM 離真值還差 0.074。這個殘餘落差有一部分不是雜訊， 是模型設錯——而且錯在這一頁自己寫的那個工作模型上。

那個模型假設 $A_0$ 的效果與 $A_1$ 的取值無關（右邊只有 $a_0$ 與 $a_1$ 兩個主效果，沒有交互作用項）。 但真值算出來的兩個對比其實不相等：$A_1$ 固定在不治療時，$A_0$ 的對比是 -0.612；$A_1$ 固定在治療時是 -0.635， 差了 0.022。工作模型把這兩個不同的量當成同一個參數在估， 所以它報出來的東西是兩者的一種折衷。

這件事值得寫出來，因為 MSM 最常被誤解的地方就是這裡：加權處理的是干擾， 它不會替你把結果模型設對。 你仍然要決定要不要放交互作用項、要不要放時間的函數、 劑量要不要當連續變項。權重模型設錯會讓權重失真，結果模型設錯會讓估計目標本身跑掉， 兩個是各自獨立的失敗管道。

變異數為什麼不能直接讀

兩條可用的路：

• Bootstrap（本頁採用）。 重抽個體，而且每一次重抽都要重配權重模型—— 只重跑結果模型、把權重當成算好的常數帶過去，就等於把上面第 2 點的錯誤原封不動搬進 bootstrap 裡。 本頁跑 500 次，表上所有的區間都來自它。
• 穩健（sandwich）標準誤。 便宜、多數軟體按一個選項就有，但它把權重當成已知， 對穩定化權重通常偏保守。本頁也一併算了：MSM 的 $A_0$ 是 0.045， bootstrap 是 0.045；$A_1$ 是 0.044 對 0.046。兩者在這份資料上很接近，但這是資料給的，不是通則。

IPCW：另一個加權，而且可以跟 IPTW 相乘

上面處理的是「誰接受了治療」。真實的追蹤資料還有第二個選擇性的問題：誰還留在追蹤裡。 病人失聯、退出、換醫院，而退出的人往往跟留下來的人不一樣——如果退出的機率跟共變項和治療有關， 那麼「還在追蹤裡」本身就是一個被條件化的變項，偏誤的結構跟干擾一模一樣。

處理方式也一模一樣：設限加權（inverse probability of censoring weighting, IPCW）。 形式跟 IPTW 完全平行，只是把「接受治療的機率」換成「在該時點仍未被設限的機率」：

實作上就是多配一個 logistic 模型（結果是「這一期有沒有被設限」）， 把每個還在追蹤裡的人依他留下來的機率的倒數加權——留在追蹤裡的人要替那些跟他相似、 但已經退出的人發言。

本頁的模擬沒有設限，所以只示範了 IPTW 那一半。真實資料上少做這一半， 等於默默假設了「退出是完全隨機的」——那是設限與截切那一頁 講的非資訊性設限假設，而在會換藥、會停藥的族群裡，它通常是最站不住的一條。

怎麼讀報表

論文裡的 MSM 通常只留下三、四句話加一張表。要檢查的是五件事：

1. 權重模型的變項清單有沒有列出來，分子與分母分開列。 只寫「用 IPTW 校正」等於沒寫。 分子裡混進共變項是常見的實作錯誤，而它不會在任何診斷圖上現形。
2. 每個時點都算了權重嗎，還是只算了基線？ 只有基線權重的分析不是 MSM， 它是基線 IPTW，處理不了 feedback。
3. 權重的分布有沒有報。 至少要有平均值（穩定化權重應該接近 1）、最大值、 以及有效樣本數。只報平均值看不出極端權重；本頁的 ESS 是名目 n 的 92.5%， 實務上掉到五成以下就要回頭檢查重疊。
4. 有沒有截尾（truncation），截在哪裡，有沒有做敏感度分析。 截尾是用一點偏誤換變異數， 是可接受的決定，但必須事先講明並附上不截尾的結果。
5. 標準誤是怎麼來的。 要看到 bootstrap 或 robust／sandwich 字樣。 什麼都沒寫的加權迴歸，八成報的是那個不能用的標準誤。

還有一件常被跳過的：結果模型的右邊應該只有治療（與必要的基線變項），不應該有時間相依共變項。 看到 MSM 的結果模型裡放著 $L_k$，那不是 MSM，那是這一頁第三節在警告的那個分析。

常見誤用

這一頁與其他頁的關係

• 權重這個想法的起點——見 propensity score 加權（IPTW）。 那一頁講的是只有一個時點的版本；權重診斷那一套（平衡表、權重分布、ESS、截尾） 在這裡完全適用，只是要對每一個時點各做一次。
• 箭頭要怎麼畫、要條件在哪些變項上——見 DAG 與干擾的結構。 本頁的碰撞子偏誤就是那一頁的 collider 規則在時間軸上的一個實例。
• Cox 模型裡的時間相依共變項——見 時間相依共變項與長表。 那一頁教的是長表怎麼組、tmerge() 怎麼用，那份資料整理工作是本頁階段一的前置； 兩頁的分歧點在於共變項受不受暴露影響——不受影響就用那一頁的標準 Cox， 受影響就要走這一頁的加權。
• per-protocol 分析的落點——見 目標試驗模擬。 那一頁定義了估計目標（ITT 還是 per-protocol）與時間零點，本頁提供 per-protocol 的估計手段。
• 設限的假設——見設限與截切。IPCW 那一節換掉的是 非資訊性設限這個假設的處理方式，不是取消它。
• 臨床落點——見資料庫研究那一章， 健保申報資料裡的用藥紀錄天生就是時間相依暴露。

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B6-08-msm.R

回報這一頁的問題