進階已經雙重審閱，尚未人工抽查

Poisson 與負二項迴歸

計數資料為什麼不能用線性迴歸、Poisson 的等分散假設幾乎總是不成立、過度分散怎麼偵測、負二項與 quasi-Poisson 各修了什麼，以及 offset 這一步——次數與發生率的差別是醫學研究最常搞錯的地方。

這一類結果變項

第三大類的臨床結果是次數：一年發作幾次癲癇、一季跌倒幾次、住院期間感染幾次、一個月急診幾趟。它們的共同特徵是非負整數、沒有上限、而且通常右偏—— 大多數人是小數字，少數人非常大。

線性迴歸處理連續結果、邏輯迴歸處理二元結果，計數要的是第三種連結函數。

計數資料還有第二個特徵，而它比第一個更重要：次數幾乎永遠要配一段觀察時間才有意義。「發作 5 次」在追蹤 2 週與追蹤 2 年是完全不同的臨床狀況。第六節專門處理這件事，它是本頁最實用的一節。

為什麼不能直接用線性迴歸

MASS::epil 是 Thall 與 Vail 的抗癲癇藥 progabide 試驗： 59 位病人（28 位安慰劑、 31 位 progabide），隨機分配前有 8 週的基線發作次數，之後每 2 週記錄一次、共四次。

本頁的主要分析在病人層級（每人一列，結果是整段 8 週追蹤的總發作次數），這樣同一個人四次觀察之間的相關性就不會被偷偷忽略掉。

這個結果變項的樣子：中位數 16、最大值 302、平均 33.05、變異數 2,078。記住最後兩個數字，第四節會用到。

Poisson 迴歸

Poisson 迴歸把「平均次數」放進對數尺度：

\log \mathbb{E}[Y \mid x] = \beta_0 + \beta_1 x_1 + \cdots + \beta_k x_k

三個直接的後果：預測值取指數之後永遠是正的； $\exp(\beta)$ 是發生率比 （rate ratio，也叫 incidence rate ratio, IRR），與 OR、HR 一樣是個乘法效果—— 另外要注意，下面每張表的 RR 欄指的都是這個發生率比，不是邏輯迴歸那一頁同樣縮寫成 RR 的相對風險（risk ratio）；以及模型自帶了「變異數隨平均值增加」的性質。

變項	RR	95% CI	p
Progabide vs placebo	0.983	0.894–1.080	0.715
log(baseline count), per unit	3.405	3.195–3.629	< 0.001
log(age), per unit	1.800	1.451–2.233	< 0.001

兩個連續變項是取過對數才放進模型的，所以表上的「per unit」指的是 log 尺度上的一單位，不是原始尺度的一次或一歲：基線發作次數變成 e 倍（約 2.72 倍）時，追蹤期的發作率乘以 3.405；年齡變成 e 倍時乘以 1.800。把它讀成「基線每多發作一次，發作率變 3.405 倍」是錯的。

看起來很乾淨。但這張表的信賴區間全部是錯的，理由在下一節。

等分散：Poisson 唯一的、也幾乎總是不成立的假設

Poisson 分布只有一個參數，於是它的變異數等於平均值——這叫等分散（equidispersion）。這不是可選項，是分布本身的性質，而模型的標準誤完全建立在它上面。

真實的臨床計數資料幾乎不長這樣。原因很好懂：Poisson 描述的是「事件彼此獨立、每個人的發生率相同」，而病人不是同質的。有人天生發作頻繁、有人本來就少，這種未被模型解釋的個體差異會讓整體變異數超過平均值，這叫過度分散（overdispersion）。

本頁的資料：整體變異數 2,078 對平均 33.05，比值 62.9。當然這是沒有校正共變項的粗比值，正式的檢查看模型的Pearson 卡方除以自由度：

\hat{\phi} = \frac{1}{n - p}\sum_i \frac{(y_i - \hat{\mu}_i)^2}{\hat{\mu}_i}

這裡是 608.0 / 55 = 11.05 （p < 0.001）。等分散成立時這個值應該在 1 附近。

左圖是觀察到的發作次數分布長條圖，疊上 Poisson 與負二項兩條期望曲線；兩條曲線一路都很接近，最大的落差在 10 到 14 次那一格——實際 17 人，兩個模型分別只預期 9.2 與 8.7 人。右圖是分組平均對分組變異數的對數座標散布圖，五個點都遠高於 Poisson 的等分散直線，貼近負二項曲線。 — 左：兩個模型的期望次數分布其實差別不大——模型裡有 log(基線次數) 這個很強的預測項，Poisson 也預測得到發作數十次的病人，所以這張圖不是過度分散的證據。右：把病人依模型預測值分成五組，各組的實際變異數遠高於等分散直線（變異數 = 平均值）——過度分散的證據是這張圖與 φ̂ = 11.05，不是左圖。產圖腳本 `figures/scripts/B2-03-poisson.R`

發作次數	實際人數	Poisson 預期	負二項預期
0-4	3	3.0	4.4
5-9	6	9.5	9.2
10-14	17	9.2	8.7
15-19	6	6.7	6.9
20-29	8	7.9	9.3
30-49	7	10.1	9.6
50-99	9	9.9	7.8
100+	3	2.7	3.0

兩條修法：quasi-Poisson 與負二項

三個模型的全部三個係數放在一起（同一份資料、同一條公式，只換變異數假設）：

係數	做法	變異數假設	RR	95% CI	p	標準誤
Progabide vs placebo	Poisson	μ	0.983	0.894–1.080	0.715	0.048
Progabide vs placebo	quasi-Poisson	φ·μ	0.983	0.718–1.345	0.913	0.160
Progabide vs placebo	負二項	μ + μ²/θ	0.769	0.574–1.031	0.079	0.149
log(baseline count), per unit	Poisson	μ	3.405	3.195–3.629	< 0.001	0.033
log(baseline count), per unit	quasi-Poisson	φ·μ	3.405	2.755–4.209	< 0.001	0.108
log(baseline count), per unit	負二項	μ + μ²/θ	2.823	2.312–3.447	< 0.001	0.102
log(age), per unit	Poisson	μ	1.800	1.451–2.233	< 0.001	0.110
log(age), per unit	quasi-Poisson	φ·μ	1.800	0.879–3.684	0.114	0.365
log(age), per unit	負二項	μ + μ²/θ	1.386	0.711–2.699	0.337	0.340

治療效果的結論沒有變，log(age) 的結論變了。 治療那三列在三個模型裡都未達統計顯著（p 0.715 / 0.913 / 0.079）。但 log(age) 在 Poisson 底下是 RR 1.800 （1.451–2.233，p < 0.001）—— 看起來是很強的預測因子；換成 quasi-Poisson，點估計一個字都沒動，信賴區間卻寬到跨過 1（0.879–3.684， p 0.114），負二項也一樣（p 0.337）。 這就是上一個 Callout 說的「假顯著」實際發生的樣子：一個在 Poisson 報表上 p < 0.001 的效應，在正確的變異數假設下退回未達統計顯著——資料沒有變，變的只有標準誤該有多大。 log(base) 則是三個模型都站得住（三列的 p 都 < 0.001），它是這份資料裡真正強的預測項。

quasi-Poisson 不換分布，只把所有標準誤乘上 $\sqrt{\hat{\phi}}$ 。點估計一個字都不動（表上每個係數的前兩列，RR 完全相同），變的只有不確定性。估計的分散參數是 11.05。

負二項（negative binomial）換了一個分布：它在 Poisson 上面再疊一層 gamma 分布來描述個體之間的差異，於是變異數變成 $\mu + \mu^2/\theta$ 。這裡估出來的 $\theta$ 是 3.69 （標準誤 0.80）； $\theta \to \infty$ 就回到 Poisson。在平均次數 33.05 這個位置，它容許的變異數是 329，而 Poisson 只容許 33.05。

正式檢定：概似比檢定統計量 377.8，p < 0.001。模型配適也差很多（AIC 848 對 472）。

兩者怎麼選：

點估計會不會變。 quasi-Poisson 不改點估計，負二項會（這裡治療效果從 0.983 移到 0.769），因為它換了加權方式—— 高計數的觀察不再被賦予那麼大的權重。這個差距本身就是過度分散嚴重程度的訊號。
負二項是一個真正的機率模型，所以有概似、有 AIC、可以做概似比檢定、可以模擬預測分布。 quasi-Poisson 不是完整的模型（沒有概似），AIC 算不出來。
quasi-Poisson 更穩健：它不假設過度分散的形狀，只假設變異數與平均值成正比。真實的變異數結構若不是負二項那個形狀，quasi 比較不會出錯。
實務上先跑負二項，若 $\theta$ 估得很不穩、或殘差顯示變異數結構不像 $\mu + \mu^2/\theta$ ，改用 quasi-Poisson 或穩健（sandwich）標準誤。

別做的事：看到過度分散就直接改用線性迴歸，或把次數切成有無二分再跑邏輯迴歸。 後者會丟掉「整段追蹤只發作一次」與「發作數十次」的差別，那通常正是臨床上最重要的資訊。殘差與影響點怎麼看見迴歸模型診斷。

至於這個試驗的結論：兩個模型都未偵測到 progabide 與癲癇發作次數的統計顯著關聯（負二項 RR 0.769，0.574–1.031， p 0.079）。信賴區間的下界仍容得下相當程度的減少，所以正確的寫法是「本分析未達統計顯著，區間尚不足以排除中等程度的療效」，不是「progabide 無效」。

offset：次數與發生率是兩件事

這是本頁最實用的一節，也是醫學研究最常搞錯的地方。

同樣的次數配上不同的觀察時間，代表的臨床事實完全不同。 一年發作 5 次與一個月發作 5 次，數字一樣，病情差十二倍。所以模型要算的通常是發生率（rate）＝次數 ÷ 人時（person-time），而不是次數本身。

作法是在對數尺度上把觀察時間加成一個係數固定為 1 的項，這叫 offset：

\log \frac{\mathbb{E}[Y]}{T} = \beta_0 + \beta_1 x_1 + \cdots \quad\Longleftrightarrow\quad \log \mathbb{E}[Y] = \log T + \beta_0 + \beta_1 x_1 + \cdots

offset(log(T)) 與「把 log(T) 當成一個普通共變項」不一樣：後者會去估一個係數，等於讓資料決定「時間的影響有多大」，而時間與次數的關係在定義上就是一比一。

這份資料自己就示範了忽略 offset 的後果。基線期是 8 週一筆，追蹤期是 2 週一筆。把兩種紀錄疊在一起問「隨機分配之後，發作有沒有變少」：

兩張長條圖。左圖比較每筆紀錄的平均發作次數，基線期明顯高於追蹤期。右圖比較每人每週的發作率，兩者幾乎一樣高。 — 左：每筆紀錄的平均次數，基線 31.2 對追蹤 8.3。右：同一份資料換算成每人每週，基線 3.90 對追蹤 4.13。差別全部來自觀察窗長度。產圖腳本 `figures/scripts/B2-03-poisson.R`

模型	追蹤期 vs 基線期的比值	95% CI	解讀
沒有 offset（比次數）	0.265	0.248–0.282	「發作減少了 74%」
有 offset（比發生率）	1.059	0.993–1.128	未偵測到發作率的變化

第一列的「減少 74%」百分之百是觀察窗長度造成的假象，一次治療效果都沒有。基線每筆蓋 8 週、追蹤每筆蓋 2 週，比值 0.265 幾乎就是 2 除以 8。

而它的信賴區間非常窄、p 值極小——統計上「非常顯著」，臨床上完全是垃圾。 這正是為什麼「顯著」不能拿來判斷分析對不對。

真實研究裡這個錯誤長什麼樣：

各組追蹤時間不等（有人早退出、有人死亡），卻直接比事件次數
健保資料庫各人的納保月數不同，卻直接比就醫次數
各醫院收案期間不同，卻直接比通報件數
各科病床數不同，卻直接比感染件數（分母應該是床日數）

判斷法則很簡單：問這個結果「每什麼」。 每人年、每千床日、每百人月—— 只要答得出「每什麼」，那個東西就該進 offset。

動手跑一次

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

# 病人層級：一人一列，避免同一人四次觀察的相關性被忽略
pt <- aggregate(y ~ subject + trt + base + age, data = epil, FUN = sum)

# Poisson
m_p <- glm(y ~ trt + log(base) + log(age), data = pt, family = poisson)
summary(m_p)

# ⚠️ 一定要自己算分散度，summary() 不會給
sum(residuals(m_p, type = "pearson")^2) / m_p$df.residual

# quasi-Poisson：點估計不變，標準誤放大
summary(glm(y ~ trt + log(base) + log(age), data = pt, family = quasipoisson))

# 負二項
m_nb <- glm.nb(y ~ trt + log(base) + log(age), data = pt)
summary(m_nb)                       # theta 在最後
2 * (logLik(m_nb) - logLik(m_p))    # 概似比統計量，p 值要除以 2

# offset：把觀察時間放進去，比的就從次數變成發生率
long <- rbind(
  data.frame(subject = pt$subject, phase = "Baseline",  count = pt$base, weeks = 8),
  data.frame(subject = epil$subject, phase = "Follow-up", count = epil$y, weeks = 2))
glm(count ~ phase, data = long, family = poisson)                      # ❌ 比次數
glm(count ~ phase + offset(log(weeks)), data = long, family = poisson) # ✅ 比發生率

驗證環境：R 4.6.0 + MASS 7.3.65。glm.nb() 在 MASS 裡，不必另外裝套件。

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

epil = sm.datasets.get_rdataset("epil", "MASS").data
pt = (epil.groupby(["subject", "trt", "base", "age"], as_index=False)["y"].sum())

m_p = smf.glm("y ~ trt + np.log(base) + np.log(age)", data=pt,
              family=sm.families.Poisson()).fit()
print(m_p.summary())
print(m_p.pearson_chi2 / m_p.df_resid)        # 分散度

m_nb = smf.glm("y ~ trt + np.log(base) + np.log(age)", data=pt,
               family=sm.families.NegativeBinomial(alpha=1 / 3.69)).fit()

# offset 用 exposure=（statsmodels 會自己取 log）或 offset=np.log(weeks)
# smf.glm("count ~ phase", data=long, family=sm.families.Poisson(),
#         exposure=long["weeks"]).fit()

statsmodels 的負二項要先估 alpha（= 1/theta）再代回去；R 的 glm.nb() 是同時估的，所以兩邊的 theta 可能有小差異。

常見誤用

誤用	為什麼錯
用線性迴歸跑次數	預測值可能為負、變異數不是常數
取 log 後跑 OLS	有 0 就沒定義；加常數會讓係數依你加多少而變
跑完 Poisson 不檢查分散度	過度分散會讓標準誤被低估、假顯著
各單位觀察時間不同卻直接比次數	要用 offset 比發生率，否則量到的是觀察窗長度
把 `log(T)` 當普通共變項而不是 offset	時間與次數的關係在定義上是一比一，不該去估它
把次數二分成「有沒有」再跑邏輯迴歸	丟掉「只發作一次」與「發作數十次」的差別
用 Poisson 的 p 值下結論而分散度很大	那個 p 值是在錯誤的變異數假設下算的
把 quasi-Poisson 的 AIC 拿來比模型	quasi 沒有概似，AIC 沒有定義
零特別多就直接套零膨脹模型	先確認零是「結構性的」（有些人根本不可能發生）還是只是過度分散；後者負二項就夠
重複測量的資料當成獨立列	組間比較的標準誤會太小；要用GEE 或混合效應模型
把 IRR 當成 OR 或 HR 解讀	分母不同：IRR 比人時、OR 比勝算、HR 比瞬時風險
未達顯著寫成「沒有效果」	應寫「未偵測到差異」並附上區間

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B2-03-poisson.R

讀讀看這張圖

答案取自產生本頁圖表的同一份統計輸出，不是另外打上去的。

要判斷這份癲癇發作次數資料需不需要處理過度離散，該看哪一個數字？

看答案與解析

正確答案: 11.1——Poisson 模型配適之後殘餘的離散比

該看的是 11.1，也就是模型配適之後剩下多少離散（Pearson 卡方除以自由度）。原始的 62.9 只是第一個訊號：共變項本來就會解釋掉一部分變異，所以配適前的比值一定偏大，用它來判斷會過度反應。33.1 是平均發作次數，它跟離散程度無關——過度離散問的是變異數與平均的關係，不是平均本身大不大。這裡解釋完仍然遠大於一，所以標準誤不能照 Poisson 的算法報。

如果改用一般線性迴歸去配適發作次數，輸出裡哪一個數字最能說明這個模型形式不對？

看答案與解析

正確答案: 最小的預測值是 -5.7——模型預測出了不可能存在的負次數

-5.7 是所有預測值裡最小的那個，而發作次數不可能是負的——線性迴歸沒有任何機制阻止它預測負數，所以這不是精度不夠，是模型的形式錯了。而 -7.2 只是文中那個示例個案，一個個案偏掉確實不能說明什麼，但這裡的問題是整個配適值的下緣都在零以下。0.0 是實際觀察到的最小值：資料合法不代表模型合法，看到預測區間下界是負數就該回頭想結果變數的分布。

這個 Poisson 模型有三項，其中只有一項回答了這個試驗的問題。Progabide 相對於安慰劑的發生率比是哪一個？

看答案與解析

正確答案: 0.983——治療分組那一項，信賴區間跨過一，未達統計顯著

0.983 是治療那一項，它的信賴區間跨過一，未達統計顯著——這代表這份資料沒有偵測到 Progabide 與安慰劑的差別，不代表兩者相同。3.405 是 log(baseline count)、1.800 是 log(age)：這兩個大得多的數字很容易被讀成「效果很大」，但它們是共變項，回答的不是試驗的問題。報表上唯一回答試驗問題的那一項，正好夾在那兩個顯眼的數字中間。

用到這個方法的章節

延伸觀看

ENVery Normal· 12 min把線性、邏輯、Poisson 收進同一個框架。看完會發現這三章其實是同一件事換連結函數。

繁中MeDA School（洪弘）· 41 min繁中、研究所級的廣義線性模型講義。本站盤點的中文影片裡沒有專講計數迴歸的一集，這支是最接近的。

素材來源與授權

本頁為原創內容