校正後的相對風險：log-binomial 與 modified Poisson

為什麼需要這一頁

邏輯迴歸那一頁的結論是：世代研究與 RCT 有分母， 可以直接算相對風險，這時候還報 OR 等於多要一個假設。那一頁也順手指出， 想要校正後的 RR 可以用 log-binomial 迴歸或 Poisson 配穩健標準誤。

聽起來像是換個 family 參數就好。實際上不是——log-binomial 很常跑不動， 而它跑不動的方式會讓第一次遇到的人以為是自己寫錯。

所以這一頁不按照教科書的順序寫（先講原理、再講方法、最後提一句「有時候不收斂」）， 而是按照你在自己的電腦上會遇到的順序寫：先成功一次、再失敗一次、再修好、 再換一條路。中間那次失敗是這一頁的重點之一，不是插曲。

這一頁的例子

同一份 MASS::birthwt：189 位產婦，59 個低出生體重， 盛行率 31.2%。抽菸組 30 / 74， 未抽菸組 29 / 115， 所以未暴露組的風險是 25.2%——這個數字整頁都會用到。

校正的變項是年齡、產婦體重、族裔，跟邏輯迴歸那一頁重疊， 好讓兩邊的估計可以直接對照。

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"))
ADJ <- low ~ smoke_f + age + lwt + race_f

# 第一步：單變項的 log-binomial，直接收斂
fit_uni <- glm(low ~ smoke_f, data = bw, family = binomial(link = "log"))
exp(coef(fit_uni))["smoke_fYes"]

# 第二步：加共變項，預設起始值會失敗
# try() 只是讓下面幾步接得下去；直接呼叫會中止，那正是這一步的重點。
try(glm(ADJ, data = bw, family = binomial(link = "log")))

# 第三步：給起始值。log(0.3) 是「基線風險約三成」，斜率都從 0 開始
fit_lb <- glm(ADJ, data = bw, family = binomial(link = "log"),
              start = c(log(0.3), rep(0, 5)))
exp(cbind(RR = coef(fit_lb), confint.default(fit_lb)))
max(fitted(fit_lb))                 # 收斂後有沒有貼在 1 附近

# 第四步：modified Poisson。三明治估計手寫版（HC0）
fit_p <- glm(ADJ, data = bw, family = poisson)
robust_vcov <- function(model) {
  X <- model.matrix(model)
  r <- model$y - fitted(model)
  bread <- solve(crossprod(X, X * model$weights))   # (X'WX)^-1
  bread %*% crossprod(X, X * r^2) %*% bread         # 夾住 X'diag(r^2)X
}
se <- sqrt(diag(robust_vcov(fit_p)))
cbind(RR = exp(coef(fit_p)),
      lcl = exp(coef(fit_p) - 1.96 * se),
      ucl = exp(coef(fit_p) + 1.96 * se))
summary(fit_p)$coefficients[, 2]     # 對照：Poisson 的裸標準誤

# 第五步：對照 logistic 的 OR
exp(coef(glm(ADJ, data = bw, family = binomial)))["smoke_fYes"]

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"})
ADJ = "low ~ C(smoke_f) + age + lwt + C(race_f)"

# log-binomial：statsmodels 同樣可能不收斂，用 start_params 給起始值
lb = smf.glm(ADJ, data=bw,
             family=sm.families.Binomial(sm.families.links.Log())
             ).fit(start_params=np.r_[np.log(0.3), np.zeros(5)])
print(np.exp(lb.params), np.exp(lb.conf_int()))

# modified Poisson：穩健標準誤一個參數就有
mp_ = smf.glm(ADJ, data=bw, family=sm.families.Poisson()).fit(cov_type="HC0")
print(np.exp(mp_.params), np.exp(mp_.conf_int()))

# 對照 logistic 的 OR
print(np.exp(smf.logit(ADJ, data=bw).fit().params))

第一步：單變項的 log-binomial，直接收斂

模型只換一個地方——連結函數從 logit 換成 log：

於是 $\exp(\beta)$ 直接就是相對風險，不再需要任何換算。 只放抽菸一個變項時，R 完全不抗議：

• RR 1.608 （95% CI 1.06–2.44， p 0.026）， 6 次疊代收斂。

第二步：加共變項，預設就失敗

把年齡、產婦體重、族裔加進去，同樣的寫法，R 直接報錯：

兩句是同一個錯誤。 R 會依照系統語系翻譯自己的訊息，所以你在網路上搜尋時 要用英文那句才找得到東西；如果你的環境是中文，看到的會是下面那句。 （順帶一提，這也是為什麼把錯誤訊息貼到 issue 或問人時最好附上英文版。）

錯誤的內容是「找不到有效的起始係數」。原因是 log 連結沒有把預測值限制在 1 以內： logit 的反函數怎麼樣都落在 0 與 1 之間，而 $\exp(\beta_0 + \beta_1 x_1 + \cdots)$ 可以大於 1，那就不是機率了。R 用來猜起始值的預設規則在這裡猜出一組會讓某些人的預測機率超過 1 的係數， 疊代第一步就無處可去。

第三步：給起始值就收斂了

起始值只要落在合法範圍內就行，不必接近答案。一組很好用的通用起始值是： 截距設成「大約的基線風險取對數」，所有斜率設成 0—— 意思是「先假設沒有任何關聯，基線風險大概這麼多」。

這樣設之後，7 次疊代就收斂：

• 校正後 RR 1.786 （95% CI 1.18–2.71， p 0.006）

第四步：modified Poisson，不必調任何參數

另一條路更省事：明知故犯地用錯的概似。把 0/1 結果丟給 Poisson 迴歸配 log 連結， 點估計仍然是一致的相對風險估計，然後用穩健（三明治）標準誤去修變異數。 這個組合就是文獻上的 modified Poisson（Zou 2004）。

跑起來完全不用調參：

• 校正後 RR 1.915 （95% CI 1.26–2.92， p 0.002）， 5 次疊代，沒有給任何起始值。

穩健標準誤在這裡修的是什麼

三明治估計的形狀是「兩片麵包夾一塊肉」：

麵包是模型假設下的資訊矩陣，肉是實際觀察到的殘差平方。 所以它的作用是：把「變異數等於模型說的那樣」這個假設，換成「變異數等於資料實際表現的那樣」。

第五步：對照 logistic 的 OR

同一份資料、同一組共變項，改用 logistic：

• 校正後 OR 2.870 （95% CI 1.36–6.05， p 0.006）

在盛行率 31.2% 的這份資料上， OR 比 log-binomial 的 RR 大了約 61% （比 modified Poisson 的 RR 大約 50%）。

「罕見疾病假設」在哪裡開始失效

右邊那條線是嚴格的直線，不是碰巧看起來像。把換算式整理一下就看得出來： 給定 OR 與未暴露組風險 $p_0$，

也就是說，「OR 比 RR 大多少」與未暴露組風險成正比，比例常數就是 $\mathrm{OR} - 1$。 所以沒有任何一個基線風險是「安全門檻」——落差是從 0 開始一路線性長上去的， 只是低基線風險那一段長得慢。$p_0 \to 0$ 時落差趨近 0，這才是罕見疾病假設真正的內容。

把校正後的 OR 固定住，看它在不同基線風險下對應的 RR：

該用哪一個

常見誤用

延伸

• 邏輯迴歸與勝算比——OR 從哪來、為什麼結果常見時它一定誇大
• 邊際估計與 G-computation——不換模型也能得到風險差與 NNT
• Poisson 與負二項迴歸——Poisson 迴歸原本要解的問題（計數結果）
• 群集與重複測量——穩健標準誤的另一個用途，方向相反

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B2-09-adjusted-rr.R

回報這一頁的問題