交互作用與次族群分析

交互作用在回答什麼問題

一篇試驗的主要結果給的是平均效果。看完那個數字，臨床上的下一個問題幾乎一定是： 「這個藥對誰比較有效？年紀大的人呢？本來就高風險的人呢？」

這個問題有一個正式的名字：交互作用（interaction）， 在流行病學的脈絡下也叫效果修飾（effect modification）。 它問的是——治療效果會不會隨著某個特徵而改變。

回答它的方式不是把資料切成幾塊、各自跑一次、然後比較誰的 p 值比較小。 那是本頁要拆掉的錯誤，而且它是臨床論文裡最常見的統計誤讀之一：

「女性顯著、男性不顯著」＝「女性和男性的治療效果不同」

這個等號不成立。左邊是兩個各自對照 0 的檢定，右邊是一個兩者互相對照的檢定。 兩個層別的樣本數不同、事件數不同，光是這一點就足以讓其中一邊跨過 0.05 而另一邊沒有， 即使底下的效果完全沒有隨層別改變。這一頁會用同一份資料把兩種讀法擺在一起， 讓落差看起來像它實際的樣子。

這一頁的例子

沿用隨機對照試驗那一章的 medicaldata::indo_rct： 直腸給予 indomethacin 預防 ERCP 術後胰臟炎（post-ERCP pancreatitis, PEP）的隨機試驗。 602 位受試者，共 79 個 PEP 事件。

整體效果（未校正的邏輯迴歸）：indomethacin 組 27／295（9.2%）， 安慰劑組 52／307（16.9%）， OR 0.49（0.30–0.81）， p 0.005。

這份資料被選來教這一頁，是因為它自帶一整排預先指定的次族群變項， 而其中好幾個的層別估計看起來差很多、交互作用檢定卻完全沒有支持。 那正是最難教、也最需要被看見的組合。

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

d <- indo_rct
d$y    <- as.integer(d$outcome == "1_yes")          # 1 = 發生 PEP
d$tx   <- as.integer(d$rx == "1_indomethacin")      # 1 = indomethacin
d$male <- as.integer(d$gender == "2_male")          # 0 = 女性（參考層）

# 一個含交互作用的模型，四個係數
fit <- glm(y ~ tx * male, data = d, family = binomial)
summary(fit)

# 層別 OR：不要把資料切開跑兩次，從同一個模型算回來
exp(coef(fit)["tx"])                                 # 女性層
exp(coef(fit)["tx"] + coef(fit)["tx:male"])          # 男性層

# 男性層的 CI 要用兩個係數的變異數「和」共變異數，不能只加兩個 SE
V  <- vcov(fit)
b  <- unname(coef(fit)["tx"] + coef(fit)["tx:male"])
se <- sqrt(V["tx", "tx"] + V["tx:male", "tx:male"] + 2 * V["tx", "tx:male"])
exp(c(OR = b, lcl = b - 1.96 * se, ucl = b + 1.96 * se))

# P for interaction：概似比檢定，不是 Wald
fit0 <- glm(y ~ tx + male, data = d, family = binomial)
anova(fit0, fit, test = "LRT")

# 對每個次族群變項重複同一件事
for (v in c("gender", "sod", "pep", "psphinc", "precut", "train", "bsphinc")) {
  f  <- droplevels(d[[v]])
  m0 <- glm(y ~ tx + f, data = d, family = binomial)
  m1 <- glm(y ~ tx * f, data = d, family = binomial)
  cat(v, anova(m0, m1, test = "LRT")[["Pr(>Chi)"]][2], "\n")
}

import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from scipy import stats

d = sm.datasets.get_rdataset("indo_rct", "medicaldata").data
d["y"]    = (d["outcome"] == "1_yes").astype(int)
d["tx"]   = (d["rx"] == "1_indomethacin").astype(int)
d["male"] = (d["gender"] == "2_male").astype(int)

fit = smf.logit("y ~ tx * male", data=d).fit()
print(fit.summary())

# 層別 OR
b_tx, b_int = fit.params["tx"], fit.params["tx:male"]
print(np.exp(b_tx), np.exp(b_tx + b_int))

# 男性層的 CI：用 contrast 向量取變異數，statsmodels 會自動帶入共變異數
c = np.zeros(len(fit.params))
c[list(fit.params.index).index("tx")] = 1
c[list(fit.params.index).index("tx:male")] = 1
se = np.sqrt(c @ fit.cov_params().values @ c)
print(np.exp([b_tx + b_int - 1.96 * se, b_tx + b_int + 1.96 * se]))

# P for interaction：LRT 要自己做
fit0 = smf.logit("y ~ tx + male", data=d).fit()
lr = 2 * (fit.llf - fit0.llf)
print(lr, stats.chi2.sf(lr, df=1))

四個係數各代表什麼

glm(y ~ tx * male) 展開之後是這樣一個模型：

跑出來的四個係數：

一列一列看：

• 截距。把 tx 與 male 都設成 0，剩下的就是截距。 取指數之後 0.2108 是安慰劑組女性的 PEP 勝算—— 可以直接用手驗：那一格有 43 個事件、 204 個非事件，相除就是這個數字。

• tx 的係數，是女性層的治療效果，不是「整體的治療效果」。 取指數 0.4540 就是女性層的 OR。 這是全篇最常被誤讀的一格。

• male 的係數，是安慰劑組裡男性相對女性的效果， 不是「男性整體的風險」。取指數 0.8372。

• tx:male 才是交互作用：男性層的治療 log OR 減去女性層的治療 log OR， 也就是兩個效果之間的差，而不是任何一個效果本身。

exp(交互作用項) 是「比值的比值」

把上面那個模型的定義展開，兩個層別的治療 OR 分別是：

兩式相除，$\beta_1$ 消掉：

所以交互作用項取指數之後，是一個勝算比的勝算比（ratio of odds ratios）。 產圖腳本把等號兩邊都算出來對過：

這個恆等式是代數上的必然，不是這份資料的巧合——腳本裡用 stopifnot() 卡住，不成立就不會產出。

知道它是「比值的比值」之後，判讀方式就跟著決定了。 它的虛無值是 1（不是 0），要在 log 尺度上做檢定與算信賴區間， 而且它跟一般的 OR 一樣需要看區間： 這裡是 1.48 （0.45–4.91）， 這個區間要翻回 OR 才讀得懂它有多寬，因為它是比值不是效果：把兩端各乘上女性層的 OR 0.45，就得到與資料相容的男性層 OR 範圍。

下界 0.447 對應男性層 OR 0.20—— 比女性層的 0.45 更低，也就是男性得到的保護反而更強。 上界 4.906 對應男性層 OR 2.23—— 已經越過 1，也就是在男性層可能反而偏向有害。

同一個區間同時容得下「男性得益更大」與「男性受害」這兩個相反的方向， 本研究並未偵測到治療效果隨性別而異。

從模型算回層別估計

有了四個係數，各層的 OR 直接由定義得到：參考層是 exp(第二個係數)， 另一層是 exp(第二個係數 + 交互作用項)。

CI 就沒那麼直接了。男性層的 log OR 是兩個係數相加，它的變異數是

漏掉那個共變異數項，區間就會算錯，而錯的方向由它的正負決定。 在這個模型裡它是負的，所以漏掉它會把區間算得太寬。 上面 R 程式碼裡的 vcov(fit) 就是在取它。

次族群 forest plot 怎麼讀

這份試驗有 7 個二元的預先指定次族群變項。 把每一層的治療 OR 與 P for interaction 一起畫出來，就是論文裡那張圖：

同一批數字列成表：

讀這張圖的順序：

1. 先看整體效果那一列，它是所有層別的參考點。這裡是 OR 0.49。
2. 看每個變項的 P for interaction，判斷有沒有證據支持效果隨這個特徵改變。 本頁 8 個次族群變項的 P for interaction 落在 0.16 到 0.95 之間； 低於 0.05 的有 0 個。
3. 最後才看各層的點估計與區間，而且是看它們有沒有一致地落在整體效果附近， 不是看誰的區間有沒有壓到 1。
4. 順便看每一層的 n 與事件數。Precut 括約肌切開術「是」那一層只有 32 人、5 個事件， 區間 0.02–2.36 涵蓋了從「幾乎完全防住」到「風險增為兩倍以上」的範圍。這種層別的點估計不含什麼資訊， 而它在圖上跟其他列長得一樣顯眼。

顯著性的差異，不是差異的顯著性

現在把倒數第二欄拿出來看。在 7 個二元次族群變項裡， 有 6 個 出現「一層的 p 低於 0.05、另一層不低於」的組合。 如果拿層別 p 值當判準，這份試驗可以生出 6 種不同的 「這個藥只對某某人有效」的故事。

而同一批資料的 P for interaction 全部不顯著。下面這張圖把兩種讀法擺在一起：

以性別為例，把左邊面板的兩列拆開來看：

點估計 0.45 與 0.67 看起來確實有距離。 但是：

• 女性層的整個信賴區間，被男性層的信賴區間完全包住。 重疊區間是 0.26–0.80， 佔較窄那條區間的 100%。 也就是說，凡是與女性層資料相容的 OR，沒有一個與男性層資料不相容。
• 兩層的樣本數相差 3.78 倍 （476 對 126），事件數相差 3.94 倍 （63 對 16）。 男性層的區間比女性層寬得多，光靠這一點就足以讓它跨過 1。
• 把「差異」本身估出來：兩個 OR 的比值是 1.48 （0.45–4.91）， P for interaction 0.52。

你其實做了幾次檢定

還有一件事讓次族群的層別 p 值更不可靠：它們不只一個。

本頁一共測了 8 個次族群變項。 假設治療效果完全不隨任何一個特徵改變，每個檢定各有 0.05 的機率誤報， 至少出現一個 p 小於 0.05 的機率是

代入之後是 33.7%—— 超過三分之一。而這還只算了交互作用檢定； 如果拿的是層別 p 值，檢定次數要再乘以層數，機率更高。

實務上真實論文的次族群數量往往不只這些，多重比較的代價也就更大。 完整的機制、校正方法與「預先指定」為什麼是唯一有效的防線， 見多重比較與次族群分析。

層別越多，先崩掉的是層別估計，不是檢定

這份試驗還有一個 4 層的次族群變項：收案中心。 它沒有被畫進上面那張 forest plot，理由值得單獨看一次：

最後一個中心只收了 3 位受試者、 0 個事件，2×2 表有整整一列是零，OR 根本不存在。 倒數第二個中心有 22 位、2 個事件， OR 算得出來，但區間 0.07–22.40 的上界是下界的 336 倍，同樣不含資訊。

而交互作用檢定仍然算得出來，也不被那一層拖累：它是在 3 個自由度上做的整體檢定， P for interaction 0.89。

但要知道那個自由度裡有一格是空的。4_Case 只有 3 位受試者、0 個事件， 與它有關的交互作用參數在概似函數上不可辨識——資料裡沒有任何東西可以決定它的值。 R 是照名目上的參數個數報 3 個自由度， 不是照「資料真的支撐得起幾個參數」報。實際有效的自由度比這個數字小。

在這個例子上結論不受影響：P for interaction 0.89 不論放在 3 還是少一個自由度下都遠離顯著。 但若某個檢定剛好落在門檻附近，這件事就會變成結論的關鍵—— 看到多層次族群的整體檢定，先去數每一層的事件數，再決定那個 df 該不該相信。

LRT 還是 Wald

交互作用檢定有兩條路：

• Wald 檢定：直接看 summary() 裡交互作用項那一列的 p 值。
• 概似比檢定（likelihood ratio test, LRT）：比較「有交互作用」與「沒有交互作用」兩個模型的概似， 也就是 anova(fit0, fit, test = "LRT")。

在性別這個例子上，兩者幾乎沒有分別： Wald 0.5205、 LRT 0.5222。 本頁全程用 LRT，理由有三個：

1. 多層變項只有 LRT 做得到。 收案中心有 4 層， 交互作用是 3 個係數， 要問的是「這 3 個係數是不是同時為零」。 Wald 只給你單一係數的 p 值，逐個看等於又做了一次多重比較。
2. Wald 對參數化方式敏感，LRT 不敏感。 換一個參考層、把變項重新編碼， Wald 的 p 值會變，LRT 不會。
3. 樣本數小或格子稀疏時 Wald 會失真，而且失真的方向是低估顯著性 （這就是邏輯迴歸那頁完全分離的極端版本： 標準誤爆炸、Wald 的 p 值趨近 1）。次族群分析正好是最容易踩到這件事的場合。

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B2-06-interaction.R

回報這一頁的問題