§4.6 廣義線性模型(Generalized Linear Models)

★★★☆☆ 廣義線性模型 GLM 卜瓦松迴歸 計數資料 連結函數 課本 pp. 167–172

📖 前言:當 Y 不是數字也不是類別

到目前為止,我們學了兩種 Y 的型態:連續數值(第三章的線性迴歸) 和類別(第四章前半的分類方法)。但真實世界還有一種尷尬的 Y:計數(counts)—— 它不是連續的實數,也不是「是/否」的類別,而是 0, 1, 2, 3, … 的非負整數。

想像你在管理 華盛頓特區的公共自行車系統(Bikeshare)。 你想預測「某個小時會有幾個人租車」。這個數字:

線性迴歸在這三個點上都會出問題。本節將介紹 卜瓦松迴歸(Poisson Regression), 並將其統一進 廣義線性模型(GLM) 的框架中。

📚 理論基礎:本節內容來自 James, Witten, Hastie, Tibshirani (2023) An Introduction to Statistical Learning with Python §4.6。GLM 框架最早由 Nelder & Wedderburn (1972) 在 Journal of the Royal Statistical Society 提出,是現代統計建模的基石之一。

4.6.1 用線性迴歸硬幹 Bikeshare 資料——然後翻車

我們拿 Bikeshare 資料直接跑最小平方法線性迴歸。結果看起來還行:

係數估計值標準誤t 值p 值
截距73.605.1314.340.00
workingday1.271.780.710.48
temp157.2110.2615.320.00
weathersit[陰/霧]-12.891.96-6.560.00
weathersit[小雨/雪]-66.492.97-22.430.00

表 4.10 摘要:線性迴歸預測 Bikeshare 每小時租車人數。

乍看合理:溫度高 → 租車多、天氣差 → 租車少。但仔細看有三個致命問題:

問題一:9.6% 的預測值是負的 🔴

線性迴歸模型不「知道」租車人數不能為負。在大雪紛飛的凌晨兩點, 模型可能預測「-15 個人租車」——這在現實中毫無意義。

問題二:變異數不是常數(異質變異性)🔴

線性模型假設誤差項的變異數 σ² 是常數,但自行車資料明顯違反: 凌晨雨天平均 5 人、標準差 3.7;晴天上午尖峰平均 244 人、標準差 132。 平均值越大,變異數也越大。這違反了同質變異性假設。

問題三:Y 是整數,模型卻給連續值 🔴

線性模型 Y = β₀ + Σ Xⱼβⱼ + ε,其中 ε 是連續誤差 → Y 必然是連續的。 但租車人數只能是整數(0, 1, 2, …),連續預測在概念上就不對。

取 log 可以救嗎?

對 Y 取 log 再跑迴歸可以避免負預測,也部分緩解異質變異性(見課本 Figure 4.14 右圖)。 但解釋變得很尷尬:「Xⱼ 增加一單位,log(Y) 的平均增加 βⱼ」——誰聽得懂? 而且 Y=0 時無法取 log。所以 取 log 只是權宜之計,不是正解。

🏪 應用場景:零售店每日來客數預測

便利商店想預測「明天會有幾個人走進店裡」。Y 是計數(0, 1, 2, …), 天氣好的週末人多(平均高、變異大),颱風天人少(平均低、變異小)。 用線性迴歸會預測出「負的來客數」,明顯不合理。Poisson 迴歸是正確選擇。

4.6.2 卜瓦松迴歸——為計數而生的模型

卜瓦松分佈(Poisson Distribution)

如果隨機變數 Y 取值為非負整數 {0, 1, 2, …},且服從卜瓦松分佈,則:

\[ \Pr(Y = k) = \frac{e^{-\lambda} \lambda^k}{k!}, \quad k = 0, 1, 2, \ldots \]

這裡 λ > 0 是 Y 的期望值,而且 λ = E(Y) = Var(Y)。 這就是卜瓦松分佈的關鍵特性:平均值等於變異數。 平均值大 → 變異數也大,完美匹配 Bikeshare 的「尖峰波動大、離峰波動小」現象。

📐 生活化理解:想像一個櫃台,平均每小時有 λ=5 個人來。 有時多來幾個(8 人)、有時少來幾個(2 人),這個「來客數」的波動幅度跟 λ 有關。 λ=5 時波動有限,λ=100 時波動自然就大——卜瓦松用 λ = Var(Y) 捕捉了這個自然的規律。

卜瓦松迴歸模型

當然,λ 不該是常數——它應該隨著天氣、時間、月份而變化。 因此我們讓 log(λ) 是預測變數的線性組合

\[ \log\big(\lambda(X_1, \ldots, X_p)\big) = \beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p \]

等價於:

\[ \lambda(X_1, \ldots, X_p) = e^{\beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p} \]

這裡用 log 連結 保證 λ 永遠為正(指數函數的輸出恆正), 而且用最大概似法(MLE)估計係數,跟邏輯迴歸的估計方式相同。

程式碼示範:Poisson 迴歸實作

# Google Colab / 本機相容資料讀取
try:
    from google.colab import drive
    drive.mount('/content/drive')
    DATA_PATH = '/content/drive/MyDrive/ISLP_data/'
except ImportError:
    DATA_PATH = '/tmp/'

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

# 載入 Bikeshare 資料(來自 ISLP 套件)
from ISLP import load_data
Bike = load_data('Bikeshare')

# 查看資料結構
print("資料維度:", Bike.shape)
print(Bike.head(3))
print("\n租車人數分布:")
print(Bike['bikers'].describe())

# ---- Poisson 迴歸 ----
# 使用 statsmodels 的 GLM 搭配 Poisson family
model_pois = smf.glm(
    'bikers ~ mnth + hr + workingday + temp + weathersit',
    data=Bike,
    family=sm.families.Poisson()
).fit()

print("\n=== Poisson 迴歸結果 ===")
print(model_pois.summary().tables[1])

# ---- 比較:線性迴歸(OLS)----
model_lm = smf.ols(
    'bikers ~ mnth + hr + workingday + temp + weathersit',
    data=Bike
).fit()

# 計算線性迴歸的負預測比例
neg_pred = (model_lm.fittedvalues < 0).mean() * 100
print(f"\n線性迴歸負預測比例: {neg_pred:.1f}%")
print(f"Poisson 迴歸負預測比例: 0.0% (模型本身禁止負值)")

plt.show()

Bikeshare 資料上的 Poisson 迴歸結果

係數估計值標準誤z 值p 值
截距4.120.01683.960.00
workingday0.010.007.50.00
temp0.790.0168.430.00
weathersit[陰/霧]-0.080.00-34.530.00
weathersit[小雨/雪]-0.580.00-141.910.00

表 4.11:Poisson 迴歸預測 Bikeshare 租車人數。

係數解讀:乘法思維 🔑

Poisson 迴歸的係數要用 乘法 來解讀,而不是加法。 因為 λ = e^(β₀ + β₁X₁ + …),所以 Xⱼ 增加 1 單位,λ 會 乘以 e^(βⱼ)

📊 應用場景:客服中心來電量預測

電信公司的客服中心每小時接到 Y 通電話。Y 是計數: 週一早上 9 點來電多(λ 大、波動大)、週日凌晨 3 點來電少(λ 小、波動小)。 用 Poisson 迴歸建模 λ = f(星期幾, 時段, 節假日),就能準確安排客服人力。 線性迴歸預測的「負來電數」在這裡同樣荒謬。

Poisson 迴歸 vs. 線性迴歸——三個關鍵差異

面向線性迴歸 (OLS)Poisson 迴歸
係數解讀 加法:Xⱼ +1 → E(Y) +βⱼ 乘法:Xⱼ +1 → E(Y) × e^(βⱼ)
變異數假設 Var(Y) = σ²(常數) Var(Y) = E(Y)(隨平均值變化)
預測範圍 (-∞, +∞),可能有負值 [0, +∞),永遠非負

✅ Poisson 迴歸優點

  • 預測值永遠非負,不會出現「負的租車人數」
  • 自動捕捉均值-變異數關係(mean ↑ → variance ↑)
  • 係數有乘法解讀,在許多領域(流行病學、保險)更自然
  • MLE 估計,具備良好的大樣本性質

⚠️ Poisson 迴歸缺點

  • 假設 Var(Y) = E(Y),實務上常有過度離散(overdispersion)
  • 過度離散會使 z 值膨脹,p 值過度樂觀
  • 需要較大的 λ 才能用常態近似做推論
  • 只能處理計數,不能處理連續或有界資料

4.6.3 GLM:一個框架統治三種迴歸

現在我們學了三種迴歸——線性、邏輯斯、卜瓦松。它們看似不同,但其實共享一個 統一結構。這就是 廣義線性模型(Generalized Linear Model, GLM) 的核心洞見。

GLM 的三個組件

每個 GLM 由三個部分組成:

  1. 隨機成分(Random Component):Y 服從某個指數族分佈
    • 線性迴歸:Y ~ Normal(μ, σ²)
    • 邏輯斯迴歸:Y ~ Bernoulli(p)
    • 卜瓦松迴歸:Y ~ Poisson(λ)
  2. 系統成分(Systematic Component):線性預測子 η = β₀ + β₁X₁ + … + βₚXₚ
  3. 連結函數(Link Function):把 E(Y) 和線性預測子 η 連起來
\[ \eta = g\big(E(Y|X)\big) = \beta_0 + \beta_1 X_1 + \cdots + \beta_p X_p \]

三種迴歸的連結函數對照

模型Y 的分佈連結函數 g(μ)連結名稱反連結 g⁻¹(η)
線性迴歸 Normal(μ, σ²) g(μ) = μ 恆等連結 (identity) μ = η
邏輯斯迴歸 Bernoulli(p) g(μ) = log(μ/(1-μ)) 邏輯連結 (logit) μ = e^η / (1+e^η)
卜瓦松迴歸 Poisson(λ) g(μ) = log(μ) 對數連結 (log) μ = e^η

程式碼示範:GLM 統一介面

import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
import numpy as np
import statsmodels.api as sm
import statsmodels.formula.api as smf
from ISLP import load_data
Bike = load_data('Bikeshare')

# GLM 統一介面:只需換 family 參數
families = {
    'Gaussian (線性迴歸)': sm.families.Gaussian(),
    'Poisson (卜瓦松迴歸)': sm.families.Poisson(),
}

print("=== 同一個公式,不同 family,不同結果 ===\n")
for name, fam in families.items():
    model = smf.glm(
        'bikers ~ temp + workingday',
        data=Bike,
        family=fam
    ).fit()
    print(f"[{name}]")
    print(f"  截距: {model.params.iloc[0]:.3f}")
    print(f"  temp 係數: {model.params.iloc[1]:.3f}")
    print(f"  AIC: {model.aic:.1f}\n")

# 視覺化:比較三種連結函數
eta = np.linspace(-3, 3, 200)
fig, axes = plt.subplots(1, 3, figsize=(14, 4))

# Identity link
axes[0].plot(eta, eta, 'b-', lw=2)
axes[0].axhline(y=0, color='gray', ls='--', alpha=0.5)
axes[0].set_title('Identity: μ = η')
axes[0].set_xlabel('η')
axes[0].set_ylabel('μ')

# Logit link
mu_logit = 1 / (1 + np.exp(-eta))
axes[1].plot(eta, mu_logit, 'r-', lw=2)
axes[1].axhline(y=0.5, color='gray', ls='--', alpha=0.5)
axes[1].set_title('Logit: μ = e^η / (1+e^η)')
axes[1].set_xlabel('η')

# Log link
mu_log = np.exp(eta)
axes[2].plot(eta, mu_log, 'g-', lw=2)
axes[2].axhline(y=0, color='gray', ls='--', alpha=0.5)
axes[2].set_title('Log: μ = e^η')
axes[2].set_xlabel('η')

for ax in axes:
    ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('/tmp/glm_links.png', dpi=100, bbox_inches='tight')
print("✅ 已儲存連結函數比較圖")
plt.show()
📚 延伸閱讀:GLM 框架遠不止這三種。其他常用家族包括 Gamma(正連續資料,如保險理賠金額)、 Inverse Gaussian(偏態嚴重的正資料)、Negative Binomial(處理過度離散的計數資料)。 所有這些都可以在同一個 GLM 框架內用 MLE 估計。

🔬 GLM 的深層洞見:對 Hermes 系統設計的啟發

GLM 給我們的最大教訓不是數學,而是設計哲學

統一的介面 + 可替換的引擎 = 最大彈性
GLM 的精髓在於:η = Xβ(線性預測子)不變,只換 g(·)(連結函數)和分佈家族, 就能從線性迴歸切換到邏輯斯迴歸再切到卜瓦松迴歸。

這對 Hermes 子代理(sub-agent)系統的啟發:

這正是 subagent-driven-development skill 設計中的核心原則: 定義統一的任務 schema → 子代理自由選擇實作方式 → 結果經過審查層(critic)驗證。

📊 模型比較總覽

特性 線性迴歸 邏輯斯迴歸 卜瓦松迴歸
Y 型態 連續數值 二元類別 (0/1) 非負整數(計數)
Y 分佈 Normal Bernoulli Poisson
連結函數 Identity: g(μ)=μ Logit: g(μ)=log(μ/(1-μ)) Log: g(μ)=log(μ)
預測範圍 (-∞, +∞) [0, 1] [0, +∞)
估計方法 最小平方法 / MLE MLE MLE
變異數結構 常數 σ² μ(1-μ) μ (= λ)
適用場景 房價、溫度、銷量 違約與否、點擊與否 來客數、事故數、電話量
💡 「GLM 不只是三種迴歸的集合——它是一種分層設計的勝利: 把『線性預測子』和『Y 的分佈假設』拆開, 換一個連結函數就是一個全新的模型,但核心估計邏輯完全一樣。」

ISLP §4.6 Generalized Linear Models · 課本 pp. 167–172 · 下一節:§4.7 Lab: Logistic Regression, LDA, QDA, and KNN