R逻辑回归：glm 建模、优势比解读与分类评估

因变量是「是/否」「患病/未患病」「手动挡/自动挡」这类二分类结果时，线性回归直接失效：它预测的是连续数值，会给出 −0.3 或 1.7 这种不可能的概率值，残差方差也随预测值变化，违反同方差假设。

逻辑回归的做法是先对线性组合做一次 logit 变换（logit transformation），把取值范围从整个实数轴压缩到 (0, 1)。

拟合模型

mtcars 里 am 是 0/1 二值变量，用 glm() 配合 family = binomial 拟合。注意这里 am 保持数值型即可，取值为 0 和 1 时 R 不会报错：

fit <- glm(am ~ wt, data = mtcars, family = binomial)
summary(fit)

Call:
glm(formula = am ~ wt, family = binomial, data = mtcars)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   12.040      4.510   2.670  0.00759 **
wt            -4.024      1.436  -2.801  0.00509 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 43.230  on 31  degrees of freedom
Residual deviance: 19.176  on 30  degrees of freedom
AIC: 23.176

Number of Fisher Scoring iterations: 6

输出结构和 lm() 很像，但有三处关键差别：

系数列的名字是 z value 而不是 t value。逻辑回归的参数估计靠迭代（Fisher Scoring iterations: 6 就是迭代次数），没有 lm 那样的精确 t 分布，用的是标准正态近似——所以精度依赖样本量，n 很小的时候这列 p 值偏乐观。

系数是 log odds（对数优势），不是概率的变化量。wt 的系数 −4.024 意思是「车重每增加 1 千磅，log odds 减少 4.024」。这个数字没法直接讲给人听。

没有 R²。因为逻辑回归用最大似然而非最小二乘估计，不存在「解释了多少方差」这个量。代替它的是 deviance：零模型（只含截距）的 Null deviance 是 43.230，加入 wt 后残差 deviance 降到 19.176。两者之比的减少量可以当作一个粗糙的拟合指标（McFadden 伪 R²）：

1 - fit$deviance / fit$null.deviance

[1] 0.5564145

这个 0.556 只能和同族模型比较，不能说成「解释了 55.6% 的变异」，写论文时不要按 R² 的口径描述。AIC: 23.176 用于模型比较，越小越好。

把系数转成优势比

exp() 是解读逻辑回归的唯一实用入口：

exp(coef(fit))

 (Intercept)           wt
1.694596e+05 1.788183e-02

优势比（odds ratio，OR）= 0.0179。含义是：车重每增加 1 千磅，手动挡的优势（odds）变为原来的 0.018 倍，也就是下降约 98%。系数为负对应 OR 小于 1，两者方向一致。

这里的「1 千磅」跨度太大，实际解释时通常换算到更合理的单位。每增加 500 磅：

exp(coef(fit)["wt"] * 0.5)

       wt
0.133723

OR ≈ 0.134，即优势下降到原来的 13.4%。报告时一定要写清楚「每增加多少单位」，只给一个 OR 数值是没有信息量的。

置信区间同样要取指数：

exp(confint(fit))

Waiting for profiling to be done...
                   2.5 %       97.5 %
(Intercept) 1.837902e+02 1.827703e+10
wt          4.533119e-04 1.598747e-01

wt 的 OR 区间是 [0.00045, 0.160]，整个区间远小于 1，和 p 值结论一致。注意 confint() 默认走的是 profile likelihood（似然剖面）方法，所以要打印一句 Waiting for profiling to be done...。想跳过这一步用快速的 Wald 区间，写 confint(fit, method = "Wald")，但小样本下 profile 区间更可靠。

截距的区间宽到 1.8e+10，这很正常——它对应的是车重为 0 时的优势，而这个场景在数据里根本不存在。

预测概率

默认情况下 predict() 返回的是 log odds，必须显式指定类型：

p <- predict(fit, type = "response")
round(head(p, 8), 3)

        Mazda RX4     Mazda RX4 Wag        Datsun 710    Hornet 4 Drive
            0.817             0.616             0.937             0.290
Hornet Sportabout           Valiant        Duster 360         Merc 240D
            0.142             0.132             0.089             0.311

type = "response" 会套一层 plogis() 把 log odds 转成概率。忘了写这个参数，得到的会是 −1.5、2.3 这样的数值，画出来的曲线也不在 0 到 1 之间——这是初学阶段最高频的错误之一。

想预测新车的概率，用 newdata：

predict(fit, newdata = data.frame(wt = c(2.0, 3.0, 4.0)), type = "response")

        1         2         3
0.9818796 0.4921156 0.0170315

2 千磅的车有 98% 的概率是手动挡，4 千磅的车只剩 1.7%。概率随车重单调下降，符合系数为负的预期。

混淆矩阵与准确率

概率要变成分类结论，需要选一个阈值。最简单的是 0.5：

pred <- ifelse(p > 0.5, 1, 0)
table(predicted = pred, actual = mtcars$am)
mean(pred == mtcars$am)

         actual
predicted  0  1
        0 18  2
        1  1 11
[1] 0.90625

32 辆车里判对 29 辆，准确率 90.6%。拆开看两个方向：19 辆自动挡里判对 18 辆（特异度 94.7%），13 辆手动挡里判对 11 辆（灵敏度 84.6%）。样本这么小，每个错误都值 7 个百分点，所以不要盯着单一的准确率数字做决定。

第三个指标是精确率（precision），它换了个方向问：判为手动挡的那些车，有多少真的是手动挡。从混淆矩阵的对应行算出来：

tab <- table(predicted = pred, actual = mtcars$am)
tab["1", "1"] / sum(tab["1", ])

[1] 0.9167

12 辆被判为手动挡，其中 11 辆判对，精确率 91.7%。精确率和灵敏度容易混：灵敏度问「真的手动挡有没有被找出来」，分母是实际的手动挡数量；精确率问「找出来的准不准」，分母是模型判为手动挡的数量。类别不平衡时这两个值会拉开很大，只看其中一个都会被误导。

阈值怎么选

0.5 不是定理，只是一个默认值。它假设两类同等重要、且模型给出的概率是校准过的。实际研究里往往不是这样：

for (th in c(0.3, 0.5, 0.7)) {
  pr <- ifelse(p > th, 1, 0)
  cat(sprintf("threshold %.1f  accuracy %.3f\n", th, mean(pr == mtcars$am)))
}

threshold 0.3  accuracy 0.875
threshold 0.5  accuracy 0.906
threshold 0.7  accuracy 0.875

提高阈值会让判为「手动挡」的样本变少，漏判增多；降低阈值则相反。选择依据应该来自研究目的：筛查疾病时宁可误诊也不能漏诊，取低阈值；确认性诊断则取高阈值。pROC 包把不同阈值下的灵敏度与特异度一次算出来：

library(pROC)

r <- roc(mtcars$am, p, quiet = TRUE)
coords(r, c(0.3, 0.5, 0.7), input = "threshold",
       ret = c("threshold", "specificity", "sensitivity"))

  threshold specificity sensitivity
1       0.3   0.8421053   0.9230769
2       0.5   0.9473684   0.8461538
3       0.7   0.9473684   0.7692308

阈值从 0.5 降到 0.3，灵敏度从 84.6% 升到 92.3%，特异度从 94.7% 掉到 84.2%——多找出一辆手动挡的代价是多误判两辆自动挡。阈值升到 0.7 只是继续牺牲灵敏度，特异度已经到顶，没有收益。这种表比准确率那张表有用，因为它把两个方向的错误分开摆出来了。

把整条 ROC 曲线（receiver operating characteristic curve）下的面积算出来，得到一个与阈值无关的指标：

auc(r)

Area under the curve: 0.9332

AUC 取值 0.5（随机猜测）到 1.0（完全分开）。0.9332 说明这个单变量模型已经有不错的区分能力。它的好处是不依赖阈值——衡量的是模型给手动挡车打出更高概率的排序能力，所以比较不同模型时比准确率稳健。

pROC 会自动判断哪一类算「正类」，依据是两类预测值的中位数。方向判断反了，AUC 就变成 1 − AUC：

r_wrong <- roc(mtcars$am, p, quiet = TRUE, direction = ">")
auc(r_wrong)

Area under the curve: 0.0668

0.9332 变成了 0.0668，正是 1 − 0.9332。拿到一个明显低于 0.5 的 AUC 时，先查 roc() 返回对象里的 $direction，而不是急着推翻模型。换数据、换正类的定义方式时方向都可能翻转。

coords() 还能直接给出精确率与召回率，和前面手算的混淆矩阵指标对上：

coords(r, 0.5, input = "threshold",
       ret = c("threshold", "specificity", "sensitivity", "precision", "recall"))

          threshold specificity sensitivity precision    recall
threshold       0.5   0.9473684   0.8461538 0.9166667 0.8461538

特异度 0.947、灵敏度 0.846、精确率 0.917，与 table() 手算的结果一致。召回率（recall）就是灵敏度的另一个名字，医学文献常用前者，机器学习文献常用后者。

要从曲线上挑一个阈值，常用的规则是 Youden 指数，也就是「灵敏度 + 特异度 − 1」最大的那个点：

coords(r, "best", best.method = "youden",
       ret = c("threshold", "specificity", "sensitivity"))

  threshold specificity sensitivity
1 0.3196104   0.8947368   0.9230769

结果是 0.32，和上面表格里 0.3 附近表现最好互相印证。Youden 指数对两类错误同等看待，如果研究场景里漏诊的代价远高于误诊，应该改用 best.method = "closest.topleft"，或者直接按目标灵敏度反查阈值。

另外，准确率在类别不平衡时几乎没有意义。mtcars 里自动挡占多数：

mean(mtcars$am == 0)

[1] 0.5938

一个把所有车都判为自动挡的模型准确率就是 59.4%，而它一辆手动挡也找不出来——灵敏度为 0。类别越不平衡这个数字越高，这时应该看灵敏度、特异度、精确率或者 AUC。

加变量与模型比较

fit2 <- glm(am ~ wt + hp, data = mtcars, family = binomial)
summary(fit2)

Call:
glm(formula = am ~ wt + hp, family = binomial, data = mtcars)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) 18.86630    7.44356   2.535  0.01126 *
wt          -8.08348    3.06868  -2.634  0.00843 **
hp           0.03626    0.01773   2.044  0.04091 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 43.230  on 31  degrees of freedom
Residual deviance: 10.059  on 29  degrees of freedom
AIC: 16.059

Number of Fisher Scoring iterations: 8

AIC 从 23.176 降到 16.059，加入马力确实改善了模型。但两个自变量之间的相关系数是 0.659，存在中等程度的共线性，wt 的系数从 −4.024 变成 −8.083 就是这个原因——系数含义变了（现在是「控制马力不变时车重的效应」），不能直接和前一个模型比较大小。

AIC 差多少才算改善没有硬标准，嵌套模型之间可以用似然比检验给出一个正式的 p 值。R 里就是 anova()：

anova(fit, fit2, test = "Chisq")

Analysis of Deviance Table

Model 1: am ~ wt
Model 2: am ~ wt + hp
  Resid. Df Resid. Dev Df Deviance Pr(>Chi)
1        30     19.176
2        29     10.059  1    9.117 0.002532 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Deviance 一列就是两个模型的残差偏差之差 19.176 − 10.059 = 9.117，自由度差 1（少了一个参数）。在原假设「hp 的系数为 0」下，这个统计量服从自由度为 1 的卡方分布，Pr(>Chi) 是 0.002532。这就是似然比检验（likelihood ratio test），与 lm 里 anova() 给出的 F 检验对应，只是从 F 分布换成了卡方分布。

用似然比检验有两个前提：两个模型必须嵌套（小模型的自变量是大模型的子集），且用同一批数据拟合。拿两个自变量集合互不包含的模型做这个检验是没有意义的。非嵌套模型之间只能比 AIC 一类的信息准则。

候选模型多于两个时，逐个跑 summary() 再手工抄数字容易出错。broom 包的 glance() 把一个模型压成一行：

library(broom)

glance(fit)[, c("deviance", "df.residual", "AIC", "BIC")]
glance(fit2)[, c("deviance", "df.residual", "AIC", "BIC")]

# A tibble: 1 × 4
  deviance df.residual   AIC   BIC
     <dbl>       <int> <dbl> <dbl>
1     19.2          30  23.2  26.1
# A tibble: 1 × 4
  deviance df.residual   AIC   BIC
     <dbl>       <int> <dbl> <dbl>
1     10.1          29  16.1  20.5

AIC 与 BIC 都从 23.2、26.1 降到 16.1、20.5，两个准则方向一致。把几个模型的 glance() 结果用 do.call(rbind, ...) 拼起来，就是一张可以贴进报告的模型比较表。

两个必须知道的坑

完全分离（complete separation）。如果某个自变量能把两类完美分开，系数的最大似然估计不存在，会朝着正负无穷发散。构造一组完全分开的数据就能看到：

sep <- data.frame(x = 1:8, y = c(0, 0, 0, 0, 1, 1, 1, 1))
m <- glm(y ~ x, data = sep, family = binomial)
round(summary(m)$coefficients, 2)

            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -206.12   365130.9       0        1
x              45.80    80643.9       0        1

同时会打印一句警告：

Warning message:
glm.fit: fitted probabilities numerically 0 or 1 occurred

8 个样本、一个再普通不过的自变量，系数却涨到 45.80，标准误 80643.9——比系数本身大三个数量级，p 值算出来是 1。这些数字没有统计意义。警告信息是唯一的线索，而它不会中断运行，模型照常返回结果。看到系数绝对值异常大、标准误离谱，先怀疑分离问题。处理办法是减少自变量、增加样本，或改用惩罚似然方法（logistf 包的 logistf()，或 brglm2 包），后者在分离情况下仍能给出有限估计。

样本量。逻辑回归的经验法则是「每个自变量至少对应 10 到 20 个较少发生的那一类事件」。这里手动挡只有 13 辆，只够放一个自变量。硬塞第三个进去：

fit3 <- glm(am ~ wt + hp + qsec, data = mtcars, family = binomial)

Warning messages:
1: glm.fit: algorithm did not converge
2: glm.fit: fitted probabilities numerically 0 or 1 occurred

13 个事件配 3 个自变量，平均每个变量只有 4.3 个事件，求解直接不收敛。上面那个双变量模型虽然收敛了（13 ÷ 2 = 6.5 个事件），也已经低于 10 的门槛，wt 的系数 −8.083 和单变量时的 −4.024 差了一倍，那个结论别当真写进论文。样本不够时优先减少自变量，而不是硬拟合。

下一步

因变量是计数或者多分类时，glm() 的 family 还能换成 poisson、quasipoisson 等，语法完全一致。分类型自变量的分组汇总可以在 /python/pandas/groupby-agg 里看到 Python 的对照写法；如果你的数据是重复测量，别用普通 glm，见 /r/modeling/mixed-models。