邊際估計與 G-computation

這一頁在還兩張支票

邏輯迴歸那一頁講到 OR 不可塌縮時，撂下了一句 「想要一個可以在族群層次解讀的效果量（例如要算 NNT），應該報邊際估計」， 常見誤用表裡也寫了「拿模型的 OR 去算 NNT 或族群層次的影響」是錯的， 需要的是「標準化 / G-computation」。

兩句話都沒有說怎麼做。這一頁做給你看。

用的是同一份 MASS::birthwt、同一個結果變項（新生兒體重是否低於 2500 g）， 共變項也刻意跟那一頁重疊，好讓兩頁的數字可以直接對照。

條件 OR 為什麼不能拿去算 NNT

邏輯迴歸的係數取指數之後是 OR，而且是條件（conditional）OR： 它的完整念法是「在年齡、產婦體重、族裔、高血壓病史、子宮易激性都固定不動的前提下， 抽菸讓低出生體重的勝算變成幾倍」。

問題就在「固定不動」這五個字。臨床上要回答的問題通常長這樣：

如果這個族群裡所有人都不抽菸，低出生體重會少多少？

這句話裡沒有任何東西被固定住。族群裡的人年齡不同、體重不同、族裔不同， 問的是整個族群的平均風險會怎麼變。條件 OR 沒有回答這個問題， 它回答的是「某一個特定共變項組合底下的那一小群人」會怎麼變。

再加上兩件事，條件 OR 就更不能直接拿來算 NNT：

• NNT 的定義是風險差的倒數，而 OR 不是風險差，也不能單獨換算成風險差—— 同一個 OR 在不同基線風險下對應完全不同的風險差（這件事在 邏輯迴歸那頁的 OR 對 RR 那張圖上已經看過）。
• OR 不可塌縮，所以就算你把它換算成某個「平均」的量，換出來的也不是族群層次的那個量。 這一點在談不可塌縮性那一節專門處理。

G-computation：三個步驟

標準化（standardisation）在迴歸模型上的做法叫 G-computation，只有三步：

1. 照常配適一個模型。 這裡就是普通的邏輯迴歸，沒有任何特別之處。
2. 把整份資料複製兩份。 第一份把每個人的暴露欄位都改成「有」，第二份都改成「沒有」， 其他欄位一個字都不動。所以每一位產婦在兩份資料裡都出現一次， 一次當作她抽菸、一次當作她沒抽菸。
3. 各自預測後取平均。 兩個平均值就是兩個邊際風險，相減是風險差，相除是相對風險。

第二步是整個方法的重點：不動其他欄位，等於用這份研究自己的共變項分布當作標準人口。 所以算出來的量是「這個族群如果全部暴露 / 全部不暴露」，不是任何假想族群。

library(MASS)
data(birthwt, package = "MASS")

bw <- birthwt
bw$race_f  <- factor(bw$race,  levels = 1:3, labels = c("White", "Black", "Other"))
bw$smoke_f <- factor(bw$smoke, levels = 0:1, labels = c("No", "Yes"))

FORM <- low ~ smoke_f + age + lwt + race_f + ht + ui
fit <- glm(FORM, data = bw, family = binomial)

# 第二、三步：複製兩份、各自預測、取平均
gcomp <- function(model, data) {
  d1 <- d0 <- data
  d1$smoke_f <- factor("Yes", levels = levels(data$smoke_f))
  d0$smoke_f <- factor("No",  levels = levels(data$smoke_f))
  r1 <- mean(predict(model, d1, type = "response"))
  r0 <- mean(predict(model, d0, type = "response"))
  c(risk1 = r1, risk0 = r0, RD = r1 - r0, RR = r1 / r0, NNT = 1 / (r1 - r0))
}
gcomp(fit, bw)

# 邊際 OR 要從兩個平均風險換算，不是拿模型係數
r <- gcomp(fit, bw)
(r["risk1"] / (1 - r["risk1"])) / (r["risk0"] / (1 - r["risk0"]))
exp(coef(fit))["smoke_fYes"]        # 條件 OR，數值不一樣

# 信賴區間：重抽個體，整個流程重跑一次
set.seed(20260823)
boot <- replicate(2000, {
  ix <- sample(nrow(bw), replace = TRUE)
  d  <- bw[ix, ]
  gcomp(glm(FORM, data = d, family = binomial), d)
})
apply(boot, 1, quantile, c(0.025, 0.975))

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

bw = sm.datasets.get_rdataset("birthwt", "MASS").data
bw["race_f"]  = bw["race"].map({1: "White", 2: "Black", 3: "Other"})
bw["smoke_f"] = bw["smoke"].map({0: "No", 1: "Yes"})

FORM = "low ~ C(smoke_f) + age + lwt + C(race_f) + ht + ui"
fit = smf.glm(FORM, data=bw, family=sm.families.Binomial()).fit()

def gcomp(model, d):
    r1 = model.predict(d.assign(smoke_f="Yes")).mean()
    r0 = model.predict(d.assign(smoke_f="No")).mean()
    return dict(risk1=r1, risk0=r0, RD=r1 - r0, RR=r1 / r0, NNT=1 / (r1 - r0))

print(gcomp(fit, bw))

rng = np.random.default_rng(20260823)
boot = []
for _ in range(2000):
    d = bw.iloc[rng.integers(0, len(bw), len(bw))]
    f = smf.glm(FORM, data=d, family=sm.families.Binomial()).fit()
    g = gcomp(f, d)
    boot.append([g["RD"], g["RR"]])
print(np.percentile(np.array(boot), [2.5, 97.5], axis=0))

算出來是什麼

NNT 念成「這個族群裡每有這麼多位產婦戒菸，就少一個低出生體重的新生兒」。 這句話裡沒有「其他變項固定」，因為 G-computation 本來就沒有固定任何東西。

Non-collapsibility：兩個 OR 不相等，而且不是誰算錯

同一個模型、同一批人、同一個對比，卻有兩個 OR：

• 條件 OR 2.794
• 邊際 OR 2.419

差距是 0.375，比值 1.155。 條件 OR 離 1 比較遠——永遠是這個方向。

為什麼會這樣？因為勝算比是不可塌縮的效果量：各次群勝算比的平均，不等於平均起來那一群的勝算比。理由不是「非線性」——log 連結一樣是非線性，它的發生率比卻塌縮得好好的。

把一群風險高低不同的人各自算出勝算、再平均，跟先把他們的風險平均起來、再換成勝算， 本來就不是同一件事。模型給每個人的預測風險散得很開—— 抽菸情境下從 6.9% 到 90.7%， 不抽菸情境下從 2.6% 到 77.7%—— 散得越開，兩個 OR 差越多。

信賴區間為什麼要 bootstrap

邊際風險差是所有共變項與所有係數的一個非線性函數，不是任何一個係數。 模型報表裡沒有它的標準誤，summary() 那一欄也不是它的。

標準做法是非參數 bootstrap：重抽個體、重新配適模型、把整個 G-computation 流程再跑一次， 重複 2000 次，取百分位數。重點在「整個流程再跑一次」—— 只重抽而不重新配適模型，等於假裝係數是已知的。

那如果偷懶呢？下面三種做法在同一份資料上跑出來的風險差區間：

NNT 的信賴區間是歪的，不要寫成正負

NNT 是風險差的倒數，而倒數會把「接近 0」放大成「接近無限大」。 本頁的風險差下限是 4.2 個百分點，離 0 不遠， 於是 NNT 的上限被推到 21.7：

• 點估計 5.2
• 95% CI 2.9 到 21.7

下界離點估計只有 2.3， 上界卻離了 16.5。

預測機率圖：它是標準化曲線，不是 partial plot

右格容易被誤認成一般的 partial dependence plot（「其他變項固定在平均值，只讓體重動」）。 它不是。 那條曲線的每一個點是這樣算出來的：把所有人的體重都改成同一個值， 其他欄位維持各自原樣，預測後取平均。也就是把 G-computation 再套用一次到體重上。

差別不是吹毛求疵：

• 固定在平均值的做法會產生一個現實中可能不存在的「平均病人」—— 平均年齡、平均體重、族裔卻是某個小數點。共變項越多、越是類別變項，這個虛構人物越荒謬。
• 標準化沒有虛構任何人，它保留了整個族群的異質性，只換掉被指定的那一個變項。

兩條曲線之間的垂直距離，就是「在那個體重水準下的邊際風險差」。 體重越大，兩條線的絕對距離越小——從左端約 22.0 個百分點 縮到右端約 9.6 個百分點。

什麼時候該報哪一個

常見誤用

延伸

• 邏輯迴歸與勝算比——條件 OR 從哪來、不可塌縮的第一次登場
• 校正後的相對風險——不想繞過 OR 的另一條路：直接估 RR
• IPTW 與傾向分數加權——另一種得到邊際估計的做法， 靠加權而不是靠複製資料

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B2-08-marginal.R

回報這一頁的問題