R线性回归完全指南：lm 建模、系数解读与残差诊断

方差分析回答「组间有没有差别」，回归回答「一个变量变化时另一个变量怎么变」。两者在 R 里共用同一套底层机制，只是 aov() 帮你把系数重新整理成了方差分解表。

从一个问题开始

mtcars 里车重（wt，单位千磅）和油耗（mpg）的关系是最经典的例子。散点图当然该先看，这里直接从模型入手——summary() 一次给出的信息比一张图多：

fit <- lm(mpg ~ wt, data = mtcars)
summary(fit)

Call:
lm(formula = mpg ~ wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max
-4.5432 -2.3647 -0.1252  1.4096  6.8727

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  37.2851     1.8776  19.858  < 2e-16 ***
wt           -5.3445     0.5591  -9.559 1.29e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.046 on 30 degrees of freedom
Multiple R-squared:  0.7528,  Adjusted R-squared:  0.7446
F-statistic: 91.38 on 1 and 30 DF,  p-value: 1.294e-10

公式 mpg ~ wt 读作「用 wt 解释 mpg」。截距项默认包含，想强制过原点要写 mpg ~ wt - 1，但除非有明确的物理理由，否则不要这么做——去掉截距会让 R² 的计算方式改变，数值不再能和普通模型比较。

逐行读懂 summary 输出

这份输出有五块信息，很多人只看了中间的系数表。

Residuals 是残差（实际值减预测值）的五数概括。中位数 −0.1252 很接近 0，说明残差大体居中；第一四分位数 −2.36 与第三四分位数 1.41 量级接近，没有明显偏斜。如果中位数明显偏离 0，说明有系统性偏差没被模型捕捉。

Coefficients 表是最常看的部分，四列各有含义：

Estimate：点估计。截距 37.29 表示车重为 0 时预测油耗 37.29（现实中不存在 0 重量的车，所以截距在这里没有实际意义，只是让直线位置正确的数学常数）。wt 的系数 −5.34 是可解释的：车重每增加 1 千磅，油耗平均减少 5.34 英里/加仑。
Std. Error：估计量的标准误，衡量这个系数估得有多准。
t value：Estimate 除以 Std. Error，即 −5.3445 / 0.5591 = −9.56。
Pr(>|t|)：检验「该系数等于 0」的 p 值。wt 对应 1.29e-10，远小于 0.05。

星号那一行只是 p 值的可视化标记，别把 *** 当作「更重要」的证据——它和 p 值携带的信息完全一样。截距的 p 值通常没有解读价值。

想直接拿到置信区间：

confint(fit)

                2.5 %    97.5 %
(Intercept) 33.450500 41.119753
wt          -6.486308 -4.202635

wt 的 95% 置信区间是 [−6.49, −4.20]。这个区间的宽度比 p 值更有信息量：它告诉你效应的可能范围，而不是简单地二分成「有」或「没有」。

最后三行是模型整体的情况。Residual standard error: 3.046 是残差的标准差，含义是「用这个模型预测油耗，平均会差 3 个 mpg 左右」——直接对应业务上的预测精度。Multiple R-squared: 0.7528 表示模型解释了 mpg 变异的 75.3%。

自变量标准化：让系数可比

wt 的系数是 −5.34，hp 的系数是 −0.03，看起来前者重要得多。这个比较站不住：wt 的单位是千磅，hp 的单位是马力，两者的「一个单位」不在同一量级。把三个自变量各除以自己的标准差再看：

std <- mtcars
std[c("wt", "hp", "disp")] <- scale(std[c("wt", "hp", "disp")])

round(coef(lm(mpg ~ wt + hp + disp, data = mtcars)), 3)
round(coef(lm(mpg ~ wt + hp + disp, data = std)), 3)

(Intercept)          wt          hp        disp
     37.106      -3.801      -0.031      -0.001
(Intercept)          wt          hp        disp
     20.091      -3.719      -2.136      -0.116

标准化之后系数读作「该变量每增加 1 个标准差，油耗平均变化多少英里/加仑」。按这个口径，hp 的效应（−2.14）虽然仍小于 wt（−3.72），但已经从「几乎看不见」进入同一量级；disp 无论怎么换算都接近零。

scale() 对每一列减去均值再除以标准差，返回矩阵。因变量不标准化，保持原单位，预测值才仍然是可以解释的油耗。这一步只改变系数的呈现方式，模型本身没变：R²、F 统计量和 p 值完全相同，只有系数的单位和数值不同。需要比较多个自变量的相对重要性时做这一步，只关心某个变量每变动一个实际单位的影响时用原始尺度。

R² 与调整 R²

R² 的定义是「被模型解释的方差占总方差的比例」，取值 0 到 1。它有一个让人误用的特性：往模型里加任何变量，R² 都不会下降，哪怕新变量纯属噪声。

调整 R²（adjusted R²）对此做了惩罚，按自由度折算：

summary(lm(mpg ~ wt, data = mtcars))$adj.r.squared

[1] 0.7445939

当新变量带来的解释力提升小于惩罚时，调整 R² 反而会降低——这正是判断「该不该加这个变量」的一个信号。比较嵌套模型还有一个更正式的做法：

fit2 <- lm(mpg ~ wt + cyl, data = mtcars)
anova(fit, fit2)

Analysis of Variance Table

Model 1: mpg ~ wt
Model 2: mpg ~ wt + cyl
  Res.Df    RSS Df Sum of Sq     F   Pr(>F)
1     30 278.32
2     29 191.17  1     87.15 13.22 0.001064 **

F = 13.22、p = 0.001，加入气缸数确实显著改善了拟合。注意这一步要求两个模型是嵌套关系（小模型的变量是大模型的子集），并且用同一批数据拟合，否则 F 检验没有意义。

加变量之后，系数为什么会变

单变量模型里 disp 的系数是 −0.0412，加上 wt 之后变成 −0.0177，掉了一半多：

round(c(disp_only = coef(lm(mpg ~ disp, data = mtcars))["disp"],
        disp_and_wt = coef(lm(mpg ~ disp + wt, data = mtcars))["disp"]), 4)
round(cor(mtcars$disp, mtcars$wt), 3)

  disp_only.disp disp_and_wt.disp
         -0.0412          -0.0177
[1] 0.888

disp 和 wt 的相关系数是 0.888，两者携带大量重叠的信息。单独放进模型时，disp 把这份共享的解释力全算在自己头上；wt 一进来，重叠部分被分走，disp 只剩下独立贡献。**多元回归的系数含义是「控制其他变量不变时该变量的效应」，不是它单独出现时的效应。**把单变量系数和多元系数混着解释，结论会前后矛盾。

AIC 与 BIC：比较模型的另一把尺子

调整 R² 只衡量拟合优度，不惩罚模型复杂度。AIC 和 BIC 在拟合优度上加了参数个数的惩罚项，数值越小越好：

ms <- list(lm(mpg ~ wt, data = mtcars),
           lm(mpg ~ wt + cyl, data = mtcars),
           lm(mpg ~ wt + cyl + hp, data = mtcars))

data.frame(
  model  = c("wt", "wt + cyl", "wt + cyl + hp"),
  adj_r2 = round(sapply(ms, function(m) summary(m)$adj.r.squared), 4),
  AIC    = round(sapply(ms, AIC), 2),
  BIC    = round(sapply(ms, BIC), 2)
)

         model adj_r2    AIC    BIC
            wt 0.7446 166.03 170.43
      wt + cyl 0.8185 156.01 161.87
 wt + cyl + hp 0.8263 155.48 162.81

三个指标在第三个模型上出现分歧：调整 R²（0.8263）和 AIC（155.48）都认为继续加变量更好，BIC 却回升到 162.81，判定第二个模型最优。原因是 BIC 的惩罚项是 log(n)，比 AIC 的常数 2 重得多。两者不一致时看目的：在意预测精度用 AIC，想要精简可解释的模型用 BIC。

这些数值只在同一批数据、同一个因变量上可比。换数据集或对因变量做了变换，AIC 的绝对值不再有意义。

分类自变量：因子与基准组

mtcars$am 是 0/1 编码。直接放进模型，R 会把它当成数值，系数变成「从 0 变到 1 的斜率」；多数时候你希望它作为因子（factor）：

mt <- mtcars
mt$am <- factor(mt$am, levels = c(0, 1), labels = c("auto", "manual"))
summary(lm(mpg ~ am, data = mt))

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)   17.147      1.125  15.247 1.13e-15 ***
ammanual       7.245      1.764   4.106 0.000285 ***
---
Residual standard error: 4.902 on 30 degrees of freedom
Multiple R-squared:  0.3598,  Adjusted R-squared:  0.3385

R 用哑变量编码（dummy coding）：截距 17.147 是基准组（auto）的均值，ammanual 的系数 7.245 是手动挡相对自动挡的均值差。判断基准组看 levels() 的顺序，levels = c(0, 1) 让自动挡当基准。想改基准组就调换 levels 的顺序，或者用 relevel(mt$am, ref = "manual")——这不会改变模型的拟合优度，只会改变系数的呈现方式。

有 k 个水平的因子会产生 k−1 个系数，剩下的那个水平「藏」在截距里。模型里加第二个分类变量时，R 默认仍然是哑变量编码（contr.treatment），而不是很多人以为的方差分析式编码。

交互项：`x1:x2` 与 `x1*x2`

当两个变量的效应不独立时，需要交互项。* 会同时展开主效应和交互项，: 只给交互项。

summary(lm(Sepal.Length ~ Petal.Length * Species, data = iris))

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)
(Intercept)                      4.2132     0.4074  10.341  < 2e-16 ***
Petal.Length                     0.5423     0.2768   1.959  0.05200 .
Speciesversicolor               -1.8056     0.5984  -3.017  0.00302 **
Speciesvirginica                -3.1535     0.6341  -4.973 1.85e-06 ***
Petal.Length:Speciesversicolor   0.2860     0.2951   0.969  0.33405
Petal.Length:Speciesvirginica    0.4534     0.2901   1.563  0.12029
---
Residual standard error: 0.3365 on 144 degrees of freedom
Multiple R-squared:  0.8405,  Adjusted R-squared:  0.8349

这里面每个系数都是相对 setosa 而言的：versicolor 的截距比 setosa 低 1.8056，斜率比 setosa 陡 0.2860。把它换算成三组各自的回归方程：

setosa：截距 4.2132，斜率 0.5423
versicolor：截距 4.2132 − 1.8056 = 2.4075，斜率 0.5423 + 0.2860 = 0.8283
virginica：截距 4.2132 − 3.1535 = 1.0597，斜率 0.5423 + 0.4534 = 0.9957

这正是分组各做一次回归的结果，交互模型只是额外假定了三组的残差方差相同。主效应不能单独解读：Petal.Length 的 0.5423 只是 setosa 一组的斜率，p = 0.052 说的是「setosa 组的斜率是否为零」，不能说成「花瓣长度整体不显著」。报告交互模型时，务必给出各组的分组斜率或画一张分组散点加拟合线。

预测：`predict()` 与两种区间

predict(fit, newdata = data.frame(wt = c(2.0, 3.5)), interval = "confidence")
predict(fit, newdata = data.frame(wt = c(2.0, 3.5)), interval = "prediction")

       fit      lwr      upr
1 26.59618 24.82389 28.36848
2 18.57948 17.43342 19.72553
       fit      lwr      upr
1 26.59618 20.12811 33.06425
2 18.57948 12.25426 24.90469

两组输出点估计完全相同，区间宽度差很多。confidence 区间描述的是回归线本身的不确定性，prediction 区间描述的是单辆新车的不确定性，后者要额外加上个体随机误差，所以宽得多。论文里报「模型预测的均值」用前者，报「某辆车的油耗可能范围」必须用后者。

诊断图：`plot(fit)`

plot(fit) 会依次画出四张图，按提示回车切换。它们检查的是模型假设，不是「模型好不好看」。

Residuals vs Fitted：残差对拟合值。理想情况是围绕 y = 0 的随机散点。出现 U 形或倒 U 形说明线性假设不成立，需要加二次项 I(x^2)；出现喇叭口（残差随拟合值变大而扩散）说明方差不齐，可以考虑对 y 取对数。
Normal Q-Q：残差的分位数对正态分位数。点大致落在对角线上就行，两端轻微偏离可接受。严重弯曲说明残差非正态，此时置信区间和 p 值都不太可靠。
Scale-Location：标准化残差的平方根对拟合值，把「喇叭口」放大后更容易看出来。
Residuals vs Leverage：残差对杠杆值（leverage，衡量每个观测对回归线的影响力）。这里要看的是库克距离（Cook’s distance）较大的点。杠杆高、残差又大的点会明显拉偏回归线。

想一次看全四张，可以 par(mfrow = c(2, 2)) 后再 plot(fit)。

诊断图看的是「有没有严重违反假设」，不是「必须完美」；n = 32 这样的小样本，残差图歪一点很正常。判断异常点别只凭肉眼，用数字更稳：

sort(cooks.distance(fit), decreasing = TRUE)[1:3]

Chrysler Imperial    Toyota Corolla          Fiat 128
           0.5319            0.2599            0.1930

Chrysler Imperial 的库克距离 0.53 是最大的，一般把大于 1 的点视为需要检查的对象。注意：发现异常点不等于可以删掉它，先搞清楚它是不是数据录入错误，删点必须在方法部分写明。

异方差与稳健标准误

残差图出现喇叭口时，最小二乘的标准误不再可靠，因为它假定所有观测的误差方差相同。先做检验：

library(car)
library(lmtest)

ncvTest(lm(mpg ~ wt, data = mtcars))
coeftest(lm(mpg ~ wt, data = mtcars), vcov = hccm(lm(mpg ~ wt, data = mtcars)))

Non-constant Variance Score Test
Variance formula: ~ fitted.values
Chisquare = 0.03794177, Df = 1, p = 0.84556

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)
(Intercept) 37.28513    2.42678 15.3640 9.253e-16 ***
wt          -5.34447    0.73811 -7.2408 4.636e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

ncvTest() 的 p = 0.846，没有异方差的证据。即便这样，稳健标准误也比普通标准误大：wt 从 0.559 变成 0.738。系数点估计一动不动（仍是 −5.3445），异方差只影响标准误，不影响系数本身。

hccm() 来自 car 包，返回异方差一致的协方差矩阵，lmtest::coeftest() 用它重算标准误。报告时要写明用了哪一种标准误，否则读者会按普通标准误理解置信区间。

残差自相关：Durbin-Watson 检验

按时间采集的数据，相邻观测的残差可能相关，这违反独立性假设，会让标准误偏小、p 值偏乐观。airquality 的 Month 和 Day 给出了真实的时间顺序：

aq <- na.omit(airquality)

set.seed(1)
car::durbinWatsonTest(lm(Ozone ~ Solar.R + Wind + Temp, data = aq))

 lag Autocorrelation D-W Statistic p-value
   1      0.03150895      1.935476    0.64
 Alternative hypothesis: rho != 0

D-W 统计量取值 0 到 4，2 表示无自相关，小于 2 是正自相关，大于 2 是负自相关。1.935 接近 2，p = 0.64，这里没有问题。

这个检验的前提是数据的行顺序有实际含义。同样的检验用在 mtcars 上会得到 D-W = 1.252、p = 0.026，看着「显著」——但 mtcars 的 32 行按车型名称排列，行与行之间没有时间或空间关系，这个显著只是排列顺序碰巧造成的。数据没有天然顺序时不要做这个检验。该检验的 p 值由自助法（bootstrap）算出，用 set.seed() 固定随机数才能得到可复现的结果。

多重共线性

自变量之间高度相关时，系数估计会变得不稳——标准误膨胀，甚至符号反转。先看相关矩阵：

round(cor(mtcars[, c("mpg", "disp", "wt", "hp")]), 3)

        mpg   disp     wt     hp
mpg   1.000 -0.848 -0.868 -0.776
disp -0.848  1.000  0.888  0.791
wt   -0.868  0.888  1.000  0.659
hp   -0.776  0.791  0.659  1.000

disp 和 wt 的相关系数 0.888，同时放进模型就要小心。更直接的指标是方差膨胀因子（VIF，variance inflation factor），需要 car 包：

library(car)
vif(lm(mpg ~ disp + wt + hp, data = mtcars))

    disp       wt       hp
7.324517 4.844618 2.736633

经验上按下面的档位判断：

VIF	判断
< 5	基本不用担心
5 – 10	系数不稳，报结果时要说明
> 10	严重共线，必须处理

disp 的 7.32 落在中间一档，已经偏高。处理办法有三种：只留物理意义更明确的那个变量（这里 wt 通常比 disp 好解释）、把两者合并成一个指标、或者用主成分回归。不要因为共线性就随便删变量——如果两个变量都是研究问题需要的，删掉会引入遗漏变量偏误。

交叉验证：模型的样本外表现

调整 R² 和 AIC 都是样本内指标。模型对新数据的预测能力要用交叉验证估计。把数据分成 5 份，每次留一份做测试：

set.seed(42)
folds <- sample(rep(1:5, length.out = nrow(mtcars)))

cv_r2 <- sapply(1:5, function(k) {
  train <- mtcars[folds != k, ]
  test  <- mtcars[folds == k, ]
  pred  <- predict(lm(mpg ~ wt + hp, data = train), newdata = test)
  1 - sum((test$mpg - pred)^2) / sum((test$mpg - mean(train$mpg))^2)
})

round(cv_r2, 4)
round(c(mean_cv = mean(cv_r2),
        in_sample = summary(lm(mpg ~ wt + hp, data = mtcars))$r.squared), 4)

[1] 0.8304 0.7582 0.8315 0.7336 0.7609
mean_cv in_sample
 0.7829    0.8268

样本外 R² 平均 0.783，比样本内的 0.827 低 0.044，这个差额就是过拟合的量。折与折之间在 0.734 到 0.832 之间波动，32 个样本下结果对数据划分的依赖很强。论文里报交叉验证结果要一并给出折间标准差，只报平均值会掩盖这种不确定性。

set.seed() 必须写在划分之前：每次跑出来的折不同，结果就不可复现。样本量小时单次五折波动大，更稳妥的做法是重复若干轮再取平均。

下一步

线性回归的因变量必须是连续变量。因变量是「是/否」这类二分类结果时，要换成逻辑回归，见 /r/modeling/logistic-regression；三组以上的均值比较见 /r/modeling/anova；数据有重复测量或嵌套结构时，普通 lm 的独立性假设不成立，见 /r/modeling/mixed-models。如果你也用 Python，同一套回归逻辑在 statsmodels 里的写法可以对照 /python/modeling/linear-regression/。