因果中介分析

這一頁接的是「那該怎麼做」

DAG 與該校正什麼那一頁已經把「不要這樣做」講完了：治療之後才測到的變項放進迴歸，會把你想估的那條路徑一起校正掉，報出來的既不是總效果、也不是乾淨的直接效果。

那該怎麼做？答案不是「別碰中介」。答案是：要問「這個治療有多少作用是透過 X 傳遞的」，就得先把「透過 X 傳遞的那一部分」定義成一個量，再去估那個量。 把中介變項丟進迴歸之所以不算中介分析，是因為它跳過了「定義」那一步——迴歸係數不會因為你心裡想的是直接效果就變成直接效果。

讀者最可能撞到這一頁的時機，是論文的 Methods 寫著 “mediation analysis was performed using the Baron and Kenny approach”，或是 “causal mediation analysis (VanderWeele)“。讀完這一頁之後，你應該要能問出四個問題：

1. 中介變項是治療之後才測的嗎——如果是治療前測的，它根本不是中介變項。
2. 中介與結果之間的干擾處理了嗎——即使暴露已經隨機分派，這一條也不會自動成立。
3. 有沒有考慮暴露與中介的交互作用——沒有交互作用項的分解，是一個更強的假設，不是一個更保守的做法。
4. proportion mediated 的分母穩不穩——總效果接近零時，這個比例可以是任何數字。

資料：一個 negative trial，而這正是重點

本頁的例子是 medicaldata::opt（Obstetrics and Periodontal Therapy 試驗：孕婦的牙周治療會不會降低早產）。原始資料 823 位孕婦。

• 暴露 Group：T（治療）對 C（對照），已隨機分派
• 中介變項 V3.PD.avg：治療後的平均牙周探測深度，以同一指標的基線值 BL.PD.avg 校正
• 結果 Preg.ended…37.wk：是否在滿 37 週之前生產

排除掉的人比你以為的多，而且主因是中介變項

排除的主因不是結果變項，是中介變項——只有 9 位缺結果，卻有 139 位缺治療後的探測深度。

這是中介分析特有的脆弱點：一般的療效分析只需要暴露與結果，中介分析需要暴露、中介、結果三個都有，於是完整個案分析掉的樣本必然更多。而中介變項是「治療之後才測」的東西，它的缺失往往與治療經驗本身有關（沒回診就量不到），也就更難假定是隨機缺失——見遺漏值。

分析集的兩組長這樣：

共 66 個早產事件，事件率 9.7%。未經任何校正的風險差是 −0.44 個百分點。

NDE 與 NIE：先把問題定義成一個量

用反事實記號。$Y(a, m)$ 是「把治療設成 $a$、把中介變項設成 $m$」時這個人的結果；$M(a)$ 是「把治療設成 $a$」時這個人的中介變項會是多少。

自然直接效果（natural direct effect, NDE）問的是：把治療打開，但讓中介變項停在沒有治療時的值，結果會差多少。

自然間接效果（natural indirect effect, NIE）問的是：治療固定打開，只把中介變項從「沒治療時的值」換成「有治療時的值」，結果會差多少。

兩者相加，中間那一項剛好抵銷：

這不是一個近似，也不是一個要去檢定的性質——在風險差這個尺度上，它是定義本身帶來的恆等式。本頁的分解實際算出來的加總殘差是 0，那只是浮點數的層級。

本頁三個反事實風險是這樣：

這份試驗的分解結果

主分析是帶交互作用項的 g-computation，尺度是風險差（單位：百分點）：

區間是 2000 次重抽的百分位區間，每一次重抽都把中介模型與結果模型一起重配。

所以這份試驗的正確寫法是：在這個分析集裡，未偵測到牙周治療經由治療後探測深度傳遞的間接效果，也未偵測到直接效果；總效果同樣未達統計顯著。

隨機化買到了什麼、沒買到什麼

要讓 NDE 與 NIE 從觀察資料裡被辨識出來，需要的假設是這六條：

• 一致性：觀察到的結果就是「實際接受到的那組治療與那個中介值」下的反事實結果
• 暴露與結果之間沒有干擾（本例由隨機分派保證）
• 暴露與中介變項之間沒有干擾（本例由隨機分派保證）
• 中介變項與結果之間沒有干擾（隨機分派保證不了；這一頁講的就是它）
• 沒有被暴露影響、又同時干擾中介與結果的變項
• 中介模型與結果模型都設定正確

隨機分派把其中兩條變成設計層面的事實：暴露與結果之間沒有干擾、暴露與中介之間沒有干擾。 這是真的買到了，而且買得很乾淨——治療是抽籤決定的，不會有任何東西同時決定「誰被治療」與「誰的探測深度比較深」。

但第四條完全沒有被碰到。 中介變項不是隨機分派的：誰的治療後探測深度比較深，是自然發生的。任何同時影響牙周狀態與早產的東西——吸菸、社經地位、全身性的感染負荷——都會讓「探測深度較深的人早產較多」這個關聯裡混進非因果的成分，而那正好是 $b$ 這條路。抽籤決定的是治療，不是牙周狀態。

還有第五條值得單獨講：被暴露影響、又同時干擾中介與結果的變項。 這種變項校正也不是、不校正也不是——不校正，$b$ 那條路上留著干擾；校正了，它是暴露的後代，會把 NDE 的一部分擋掉。這種情形下自然效果不再可辨識，要改用 interventional effects 那一類的估計量，而那已經超出這一頁。

第六條「模型設定正確」在中介分析裡比在一般迴歸裡更貴：兩個模型都要對，而且中介模型的錯設會經由積分傳遞到兩個效果上。中介變項本身的測量誤差同樣會把 NIE 往零壓——見測量誤差。

Baron-Kenny 為什麼不夠

Baron 與 Kenny 的逐步法是這樣：配 $M \sim A$ 拿到 $a$、配 $Y \sim A + M$ 拿到 $b$ 與 $c'$，然後宣稱間接效果是 $a \times b$、直接效果是 $c'$、總效果是 $c' + ab$。它在線性模型、而且沒有交互作用的世界裡是對的。臨床論文的結果通常兩個條件都不滿足。

不夠的第一個地方：它處理不了暴露與中介的交互作用

逐步法的 $b$ 是一個數字，意思是「中介變項每高一單位，結果的變化」——不論是治療組還是對照組。可是治療本身可能改變中介變項的作用強度，那就是一個 $A \times M$ 的交互作用項。

本頁資料裡這個交互作用項是 +0.837（95% CI −0.442 到 +2.098，log odds 尺度），區間跨過零，未達統計顯著。但「未達統計顯著」不是「可以當它是零」的授權——尤其這裡的區間寬到兩端符號相反。把交互作用項留著與拿掉，分解結果會這樣變：

單位都是百分點。總效果幾乎不動（+0.10），但直接與間接各自被移動了兩個百分點以上，而且方向相反。這正是交互作用項的作用：它不改變「總共」，它改變「怎麼分」。

所以「模型裡不放交互作用項」不是一個比較保守的選擇，而是一個額外的假設：你在斷言治療不改變中介變項的作用強度。這個假設要嘛有理由，要嘛就把兩組結果都報出來。

不夠的第二個地方：非線性模型下係數不可相乘

$c = c' + ab$ 的推導用到「係數可以相加相乘」，那是線性模型的性質。這裡的結果是二元的，結果模型是 logistic 迴歸，而 log odds ratio 不可塌縮（non-collapsible）：$Y \sim A$ 的係數與 $Y \sim A + M$ 的係數就算完全沒有干擾也不會相等，因為它們回答的是不同的問題（邊際 vs 條件）。相關的完整說明在邊際效果與校正後的相對風險。

四種算法放在一起看：

前兩列與後兩列不能直接比大小，尺度不同。要看「捷徑錯了多少」，就得在同一個尺度上比——也就是後兩列。兩者都是 log OR，都在算同一件事，差別只在 VanderWeele 的閉合式把交互作用項算了進去：Baron-Kenny 的 NIE 是 −0.147，VanderWeele 的是 −0.352，後者是前者的 2.39 倍。在同一份資料、同一個尺度上，捷徑把間接效果縮成了不到一半。

Proportion mediated：問題出在分母

proportion mediated 的定義很簡單：

本頁算出來的點估計是 9.89，百分位區間 -24.67 到 23.33。一個叫「比例」的量，點估計遠在 1 以上、區間還橫跨負值——這不是算錯，是這個量在這種情況下本來就沒有意義。

原因在分母。這份試驗的總效果是 −0.35 個百分點，本身未達統計顯著，數值上非常接近零。分母接近零時，比值對分子與分母的任何一點擾動都會爆炸，而且分母只要換個號，比值就整個翻過去。

這件事在本頁是被量化過的，不是形容詞。2000 次重抽裡：

將近一半的重抽讓分母變號，超過八成的重抽算出來的「比例」不在 0 到 1 之間。所以本頁那個區間不可以被讀成「中介比例介於某某之間」——它連一個比例都不是。

怎麼算

library(medicaldata)
data(opt, package = "medicaldata")

# 欄位名是三個點的 Preg.ended...37.wk。打錯不會報錯：opt$Preg.ends 回 NULL，
# 而 sum(is.na(NULL)) 是 0，讀起來像「沒有遺漏值」。所以這一行是承重的。
stopifnot("Preg.ended...37.wk" %in% names(opt))

# 這個 factor 的 level 補過空白："   "、"No "、"Yes"。沒有 trimws() 的
# `== "No"` 會靜默配不到任何一列，而且不會有任何警告。
yraw <- trimws(as.character(opt[["Preg.ended...37.wk"]]))

d <- data.frame(
  a  = as.integer(opt$Group == "T"),                            # 暴露，已隨機分派
  y  = ifelse(yraw == "Yes", 1L, ifelse(yraw == "No", 0L, NA)), # 結果
  m  = opt$V3.PD.avg,                                           # 中介：治療後探測深度
  bl = opt$BL.PD.avg                                            # 同一指標的基線值
)
d <- d[complete.cases(d), ]

# 兩個模型。結果模型帶 a:m 交互作用項——拿掉它是一個額外的假設，不是保守。
mfit <- lm(m ~ a + bl, data = d)
ofit <- glm(y ~ a * m + bl, data = d, family = binomial())

# g-computation：把中介變項的分布積掉，直接得到反事實風險。
# astar 決定「中介變項用哪一組治療下的分布」，a 決定「結果模型用哪一組治療」。
# 兩者可以不一樣，而那個不一樣正是跨世界的量 E{Y(1, M(0))}。
set.seed(20260823)
R <- 500
ey <- function(a, astar) {
  mu <- predict(mfit, newdata = transform(d, a = astar))
  md <- rnorm(R * nrow(d), rep(mu, R), sigma(mfit))
  mean(predict(ofit, type = "response",
               newdata = data.frame(a = a, m = md, bl = rep(d$bl, R))))
}

y11 <- ey(1, 1); y10 <- ey(1, 0); y00 <- ey(0, 0)
c(NDE = y10 - y00, NIE = y11 - y10, TE = y11 - y00)   # 前兩者相加恆等於第三者

# 信賴區間沒有閉合式，用 bootstrap。每一次重抽都要把「兩個模型」一起重配，
# 只重配結果模型會低估變異——中介模型的不確定性也會傳到 NIE 上。
# 產生本頁數字的腳本用的是 30 點 Gauss-Hermite 求積分而不是上面的隨機抽樣，
# 這樣 bootstrap 重抽的就只有資料，不會再疊一層蒙地卡羅雜訊。

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

opt = sm.datasets.get_rdataset("opt", "medicaldata").data

# 結果的 level 帶著空白填充（"   "、"No "、"Yes"），一定要先 strip，
# 否則 map 出來整欄都是 NaN，而且不會有任何警告。
y = opt["Preg.ended...37.wk"].str.strip()

d = pd.DataFrame({
    "a": (opt["Group"] == "T").astype(int),
    "y": y.map({"Yes": 1.0, "No": 0.0}),
    "m": opt["V3.PD.avg"],
    "bl": opt["BL.PD.avg"],
}).dropna()

mfit = smf.ols("m ~ a + bl", data=d).fit()
ofit = smf.glm("y ~ a * m + bl", data=d, family=sm.families.Binomial()).fit()
sigma = np.sqrt(mfit.scale)

rng = np.random.default_rng(20260823)
R = 500

def ey(a, astar):
    mu = mfit.predict(d.assign(a=astar))
    md = rng.normal(np.tile(mu, R), sigma)
    nd = pd.DataFrame({"a": a, "m": md, "bl": np.tile(d.bl, R)})
    return ofit.predict(nd).mean()

y11, y10, y00 = ey(1, 1), ey(1, 0), ey(0, 0)
print({"NDE": y10 - y00, "NIE": y11 - y10, "TE": y11 - y00})

敏感度分析：換一個時點的中介變項

同一個試驗還量了另一個時點的探測深度 V5.PD.avg。用它重跑同一套分析：

單位是百分點。NDE 與 NIE 的方向與主分析一致、量級較小（間接效果約為主分析的一半），而且同樣三個都未達統計顯著（間接效果的區間是 −6.96 到 +2.27）。

但總效果的點估計換了正負號：主分析是 −0.35 個百分點，這一組是 +0.86 個百分點。這不是兩個互相矛盾的結論，而是上一節那件事的另一個面貌——總效果貼著零的時候，連它的正負號都是不穩的，換一個時點的中介變項就足以讓它翻過去。這也順帶說明了為什麼這份資料上的 proportion mediated 不能報：分母的符號本身都還沒定下來。

怎麼讀報表

論文裡的中介分析通常給你一張路徑圖加一張表，表上有 NDE、NIE、總效果與 proportion mediated。回到開頭那四個問題，每一個都對應報表上一個具體的位置：

1. 中介變項的測量時點。 Methods 應該寫得出來「治療後第幾週測的」。如果中介變項是基線測的，它不可能是中介變項；如果測量時點與結果的發生時間重疊，那可能是反向因果。
2. 中介與結果之間的干擾。 找有沒有一句「adjusted for … in the mediator-outcome relationship」，以及有沒有未測量干擾的敏感度分析。只寫「本研究為隨機試驗」是不夠的，理由見上面那一節。
3. 有沒有交互作用項。 現代的做法（VanderWeele 那一系）預設放；只報 $a \times b$ 的論文等於宣稱交互作用是零。兩種都給的論文，看兩組數字差多少。
4. proportion mediated 的分母。 總效果的區間如果跨過虛無值，那個百分比就不該被引用。看到「本研究顯示 XX% 的效果由 Y 中介」時，先回頭看總效果那一列。

還有一個容易漏掉的：尺度。風險差、相對風險、log OR 三種尺度上的 NDE 與 NIE 是三組不同的數字，proportion mediated 也是。表格若沒有寫尺度，那張表沒辦法被複製，也沒辦法跟別篇比。

常見誤用

這一頁與其他頁的關係

• 「不要這樣做」在那一頁——見 干擾、DAG 與該校正什麼。那一頁講的是把中介變項丟進迴歸會發生什麼事，這一頁講的是那該怎麼做。
• 交互作用項本身——見 交互作用。中介分析裡的 $A \times M$ 與那一頁的交互作用是同一件事，只是它在這裡改變的是「效果怎麼分」。
• 為什麼係數不可相乘——見 logistic 迴歸與邊際效果。不可塌縮性是這一頁不用 Baron-Kenny 當主分析的技術理由。
• g-computation 的一般形式——見 邊際效果與目標試驗模擬。本頁只是把它用在一個帶跨世界量的估計上。
• 缺失值——見遺漏值。中介分析要三個變項同時有值，完整個案分析掉的樣本必然比一般分析多。
• 中介變項的測量誤差——見測量誤差。它會把 NIE 往零壓，於是「未偵測到中介」有一部分可能來自量得不準。
• 這份資料的臨床脈絡——見設計章 隨機對照試驗。

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B6-09-mediation.R

回報這一頁的問題