迴歸模型診斷

診斷在檢查什麼

跑出一組係數不代表模型是對的。線性迴歸的每一個 p 值與信賴區間，都建立在四個假設上：

1. 設定正確——放進去的變項與結果的關係形狀對了（線性、有沒有交互作用）
2. 變異數固定（homoscedasticity）——殘差的散布不隨預測值改變
3. 殘差近似常態——影響的主要是小樣本時的區間，大樣本靠中央極限定理救得回來
4. 觀察彼此獨立——這一條看不出來，要靠設計（重複測量、群集抽樣、同一家醫院的病人都會違反）。 違反之後怎麼辦，見重複測量與群集資料

診斷就是逐條檢查前三條，外加兩個報表看不出來的問題：共線性與影響點。 本頁沿用 B2-01 的線性模型與 B2-02 的邏輯模型，你已經讀過它們的報表。

殘差圖要看什麼

殘差在低預測值那一半的標準差是 628， 高預測值那一半是 643——幾乎一樣，C 圖平坦是有數字支持的。

共線性與 VIF

共線性（collinearity）是指某個變項幾乎可以被其他變項的線性組合預測出來。 變異數膨脹因子（variance inflation factor, VIF）量的就是這件事：

其中 $R_j^2$ 是「把第 $j$ 個變項當結果、用其他所有變項去迴歸」得到的 R²。 $\mathrm{VIF} = 4$ 代表那個係數的標準誤被放大了 2 倍。

B2-01 那個模型的 VIF 全部很低（最大 1.34）。 為了示範，故意再放進一個「同一個體重換算成公斤」的欄位——這在真實資料合併之後很常見。 注意換算時帶了一點量測／四捨五入誤差（lwt * 0.4536 + rnorm(n, 0, 0.5)）， 這正是真實資料合併後的樣子；如果換算是完全等比例的，R 會判定它與 lwt 完全共線（aliased）、 係數直接給 NA、VIF 是無限大，跑不出下面這張表：

影響點

離群值（outlier）、槓桿（leverage）、影響力（influence）是三件不同的事：

• 離群值：$y$ 離預測值很遠 → 看標準化殘差
• 槓桿：$x$ 落在共變項空間的邊緣 → 看 hat value；平均值是 $p/n$ = 0.042
• 影響力：把它拿掉，係數會不會變 → Cook’s distance，它是前兩者的乘積效果

只有第三件事真的要緊。$x$ 很極端但落在迴歸線上的點，槓桿高、影響力低， 它其實在幫模型固定斜率。

本模型 Cook’s D 最大的五位：

這一頁的診斷還替 B2-01 收了一個尾。那頁的年齡平方項在統計上被支持（p 0.018）， 但把 3 位超過 35 歲的產婦拿掉重跑，證據就消失了（p 0.373）。 3 個人可以決定一個 189 人模型裡的一項結論—— 這就是為什麼「顯著」與「穩健」是兩件事，也是為什麼多項式不如樣條。

二元與計數結果的診斷不一樣

線性模型的那四張圖不能直接搬到邏輯迴歸上。 原因很簡單：結果只有 0 與 1 兩個值，原始殘差只能落在兩條曲線上， 畫出來永遠是兩條帶子，看不出任何東西（這個模型的原始殘差範圍是 -0.74 到 0.91， 但形狀完全由結果值決定）。

要看的東西換成三樣：

1. 校準（calibration）——預測機率與實際發生率對不對得上。把病人依預測機率分成十組：

Hosmer-Lemeshow 檢定統計量 4.68，df = 8， p 0.792——未偵測到系統性的校準偏差。

2. 區辨力（discrimination）——見ROC 與 AUC（B4 家族）。校準與區辨是兩件獨立的事： 一個模型可以排序很準但機率整體偏高，也可以機率很準但排不出順序。

3. 影響力仍然要看——glm 物件一樣有 cooks.distance() 與 hatvalues()，判讀方式相同。

至於Poisson 迴歸，最關鍵的診斷不是殘差圖而是分散度： Pearson 卡方除以自由度是否接近 1。那一頁的例子是 11.05，也就是說標準誤要放大約 3.32 倍。

逐步迴歸：醫學論文最常見的統計錯誤之一

逐步迴歸（stepwise regression）是指讓演算法依 p 值或 AIC 自動增刪變項， 最後留下一組「顯著」的變項。R 的 step() / stepAIC()、SPSS 的 forward / backward 選項都是。

它在論文裡非常普遍，而且幾乎總是錯的。用模擬把它做給你看： 拿 B2-01 那個模型，額外塞進 10 個純粹的標準常態亂數 （與新生兒體重毫無關係），跑 stepAIC()，重複 300 次。

更狠的版本：把結果變項也換成亂數，只丟 15 個亂數預測項進去， 資料裡沒有任何東西可以找。結果是 300 次執行中：

• 平均留下 2.59 個變項
• 55% 的最終模型至少有一個變項 p < 0.05
• 平均 R² 是 0.047
• 64% 的最終模型整體 F 檢定也「顯著」

最後這一項要小心分母：這 300 次裡有 24 次 stepAIC 把變項全部刪光、只剩截距，那種模型沒有整體 F 檢定。 上面的 64% 是把這 24 次 算成「不顯著」的版本（分母 300）；若只算 276 次 真的有 F 檢定的執行，比例是 70%。 兩個數字都要看——只報後者等於把「逐步迴歸什麼都沒選中」這些最接近虛無的執行悄悄拿掉， 數字會往「逐步迴歸很會製造假象」的方向多推幾個百分點。順帶一提， 連純亂數資料都有 92% 的執行被 stepAIC 留下至少一個變項。

那變項該怎麼選

三種情況都有一個共同點：變項的去留不由同一份資料的 p 值決定。

動手跑一次

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"))

fit <- lm(bwt ~ age + lwt + smoke_f + race_f + ht + ui, data = bw)

par(mfrow = c(2, 2)); plot(fit)          # 四張標準診斷圖

# ⚠️ 標準四張圖的橫軸是預測值；單一變項的曲率會被稀釋掉
plot(bw$age, residuals(fit)); lines(lowess(bw$age, residuals(fit)))

# VIF，從定義算（不需要 car）
X <- model.matrix(fit)[, -1]
sapply(seq_len(ncol(X)), function(j) 1 / (1 - summary(lm(X[, j] ~ X[, -j]))$r.squared))

# 共線性示範：把體重換算成公斤再放一次
# ⚠️ 換算要帶一點誤差，否則兩欄完全等比例，R 會判為 aliased、係數給 NA、算不出 VIF
bw$lwt_kg <- bw$lwt * 0.4536 + rnorm(nrow(bw), 0, 0.5)
fit_col <- lm(bwt ~ age + lwt + lwt_kg + smoke_f + race_f + ht + ui, data = bw)
summary(fit_col)$coefficients[c("lwt", "lwt_kg"), ]

# 影響點：找出來，然後實際重跑看結論會不會變
cd <- cooks.distance(fit)
worst <- which.max(cd)
coef(fit)["smoke_fYes"]
coef(lm(formula(fit), data = bw[-worst, ]))["smoke_fYes"]

# 逐步迴歸會撿到什麼：塞 10 個亂數進去試一次
d <- bw[, c("bwt", "age", "lwt", "smoke_f", "race_f", "ht", "ui")]
for (j in 1:10) d[[paste0("z", j)]] <- rnorm(nrow(d))
summary(stepAIC(lm(bwt ~ ., data = d), direction = "both", trace = 0))

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

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"})

fit = smf.ols("bwt ~ age + lwt + C(smoke_f) + C(race_f) + ht + ui", data=bw).fit()

infl = fit.get_influence()
cooks = infl.cooks_distance[0]
lev   = infl.hat_matrix_diag

X = np.asarray(fit.model.exog)
vif = [variance_inflation_factor(X, i) for i in range(1, X.shape[1])]

# statsmodels 沒有內建 stepwise——這其實是件好事

常見誤用

重跑本頁的所有數字

/opt/homebrew/bin/Rscript figures/scripts/B2-04-model-diagnostics.R

延伸觀看

▶生物統計學一 93.【迴歸分析 (1)】Simple Linear Regression Model

▶Lec04 統計學(二) Ch11.1-11.5 簡單廻歸分析與相關分析

回報這一頁的問題