AMA3602 应用线性模型教程：从简单回归、多元回归到模型诊断、变量选择与混合效应模型

这篇文章根据 AMA3602 Applied Linear Models 的讲义、tutorial、summary note 和 Boston Housing 项目资料整理。它不是考试提纲，而是一篇面向“真的会用回归建模”的教程：先把简单线性回归和多元线性回归讲清楚，再进入残差诊断、变量选择、多重共线性、ridge regression，最后接到随机效应模型和一个 Boston Housing 房价建模案例。

如果把机器学习理解成“用数据拟合一个可泛化的函数”，那么线性模型就是最值得先学扎实的一类模型。它简单，但不幼稚；它的假设透明，诊断方法成熟，而且很多复杂模型的思想都可以在这里找到原型。

原始资料与数据：

• Chapter 1: Simple Linear Regression
• Chapter 2: Multiple Linear Regression
• Chapter 3-I: Residual Analysis
• Chapter 3-II: Variable Selection
• Chapter 4: Linear Random-Effects Models
• Chapter 5: Multicollinearity and Ridge Regression
• Summary: Model Adequacy Checking
• Summary: Linear Mixed Model
• Boston Housing 数据

1. 线性模型到底在做什么

线性模型的核心问题是：

$$y = f(x) + \varepsilon$$

我们观测到的是 $y$ 和 $x$，真正想学的是 $f(x)$。在线性模型里，我们假设这个关系可以近似写成：

$$y = \beta_0 + \beta_1 x_1 + \cdots + \beta_p x_p + \varepsilon$$

这里的 $\beta$ 是模型参数，$\varepsilon$ 是误差项。一个好的回归模型至少要回答四件事：

1. 变量之间有没有稳定关系？
2. 每个变量的影响方向和大小是什么？
3. 这个关系是否显著，还是只是样本噪声？
4. 模型假设是否合理，预测能不能信？

所以线性回归不是简单地调用 lm() 得到一行公式。完整流程应该是：

看数据 -> 建模型 -> 看系数 -> 做显著性检验 -> 做诊断 -> 修模型 -> 解释结果

很多初学者会停在 summary(model)，但真正的统计建模能力在后半段：你要知道模型哪里可能错，为什么错，怎么改。

2. 简单线性回归：一条直线的完整含义

简单线性回归只有一个自变量：

$$y_i = \beta_0 + \beta_1 x_i + \varepsilon_i,\quad i=1,\ldots,n$$

其中：

• $\beta_0$ 是截距，表示 $x=0$ 时的理论平均响应。
• $\beta_1$ 是斜率，表示 $x$ 每增加 1 个单位，$y$ 的条件均值平均变化多少。
• $\varepsilon_i$ 是第 $i$ 个样本没有被直线解释掉的部分。

最小二乘法的目标是让残差平方和最小：

$$SS_{Res}=\sum_{i=1}^n (y_i-\hat y_i)^2$$

直观上，模型在所有可能的直线里，选择一条让点到直线的垂直误差平方和最小的直线。

对于简单线性回归，斜率估计量可以写成：

$$\hat\beta_1=\frac{S_{xy}}{S_{xx}} = \frac{\sum_i (x_i-\bar x)(y_i-\bar y)} {\sum_i (x_i-\bar x)^2}$$

截距则是：

$$\hat\beta_0=\bar y-\hat\beta_1\bar x$$

在 R 里可以直接写：

fit <- lm(y ~ x, data = df)
summary(fit)
anova(fit)

summary(fit) 主要看系数估计、标准误、t 值、p 值和 $R^2$；anova(fit) 主要看回归平方和、残差平方和、均方和 F 检验。

2.1 t 检验：斜率是不是显著

最常见的问题是：$x$ 对 $y$ 有没有线性影响？

这可以写成假设检验：

$$H_0:\beta_1=0,\quad H_1:\beta_1\ne 0$$

如果 $\beta_1=0$，说明直线没有斜率，$x$ 的变化不能解释 $y$ 的平均变化。R 输出里的 t value 和 Pr(>|t|) 就是在做这个判断。

但要小心：显著不等于因果。线性回归能告诉你变量之间有统计关系，不自动说明一个变量导致另一个变量变化。

2.2 置信区间和预测区间

回归里有两种区间经常混淆：

• 置信区间：估计平均响应 $E(y|x_0)$ 的不确定性。
• 预测区间：预测一个新观测值 $y_0$ 的不确定性。

预测区间通常更宽，因为它不仅包含平均函数的不确定性，还包含新样本自身的随机误差。

predict(fit, newdata = data.frame(x = 10), interval = "confidence")
predict(fit, newdata = data.frame(x = 10), interval = "prediction")

如果你是在解释平均趋势，用 confidence interval；如果你是在预测一个新的样本，用 prediction interval。

3. 多元线性回归：从一条线到一个平面

多元线性回归把一个自变量扩展到多个自变量：

$$y_i = \beta_0+\beta_1x_{i1}+\beta_2x_{i2}+\cdots+\beta_px_{ip}+\varepsilon_i$$

矩阵形式更清楚：

$$\mathbf y = X\beta+\varepsilon$$

其中 $X$ 是 design matrix，第一列通常是全 1 的截距列。最小二乘估计是：

$$\hat\beta=(X^TX)^{-1}X^T\mathbf y$$

这条公式非常重要。它说明多元回归本质上是一个线性代数问题：只要 $X^TX$ 可逆，模型参数就可以被估计。

在 R 中：

fit <- lm(y ~ x1 + x2 + x3, data = df)
summary(fit)
anova(fit)

3.1 多元回归系数怎么解释

在简单回归里，$\beta_1$ 表示 $x$ 对 $y$ 的平均影响。多元回归里，$\beta_j$ 的解释要加上一句非常关键的话：

在其他变量保持不变的情况下，$x_j$ 增加 1 个单位，$y$ 的条件均值平均变化 $\beta_j$。

这就是 controlled effect。比如房价模型里，如果 RM 的系数为正，意思不是“房间数多的房子一定更贵”，而是在其他变量如 LSTAT、NOX、PTRATIO 控制住之后，RM 增加通常对应更高的 MEDV。

3.2 为什么系数可能符号反常

多元回归里常见一个现象：单独看某个变量时它和响应变量正相关，但放进多元模型后系数变成负的，或者反过来。

原因通常是：

• 自变量之间存在相关性。
• 某个变量在单变量分析中吸收了其他变量的影响。
• 模型遗漏了重要变量。
• 数据里存在高杠杆点或异常点。

所以多元回归不能只看散点图，也不能只看系数。你需要同时看相关性、VIF、残差图、变量选择结果和领域解释。

4. 回归模型的五个基本假设

线性模型常见假设包括：

1. 线性关系：$y$ 和 regressors 之间的关系至少近似线性。
2. 零均值误差：$E(\varepsilon_i)=0$。
3. 同方差：$Var(\varepsilon_i)=\sigma^2$。
4. 误差不相关：$Cov(\varepsilon_i,\varepsilon_j)=0$。
5. 正态性：$\varepsilon_i\sim N(0,\sigma^2)$，主要用于小样本推断。

这些假设不是装饰品。如果它们严重不成立，系数、标准误、p 值、置信区间和预测区间都会受到影响。

5. 残差诊断：模型错在哪里

残差定义为：

$$e_i=y_i-\hat y_i$$

它可以理解成“模型没解释掉的部分”。残差诊断的目标不是追求每个残差都很小，而是看残差有没有系统模式。

常用图形包括：

诊断图	看什么	可能问题
residuals vs fitted	是否随机水平散布	非线性、异方差
normal Q-Q plot	点是否接近直线	非正态、重尾、偏态
scale-location plot	残差波动是否稳定	异方差
residuals vs leverage	是否有高影响点	influential points

R 里最基本的诊断方式：

fit <- lm(MEDV ~ RM + LSTAT + NOX + PTRATIO, data = housing)
par(mfrow = c(2, 2))
plot(fit)

5.1 标准化残差与学生化残差

原始残差 $e_i$ 的尺度受响应变量单位影响，不方便比较，所以需要 scaled residuals。

标准化残差常写作：

$$d_i=\frac{e_i}{\sqrt{MS_{Res}}}$$

学生化残差考虑了杠杆值 $h_{ii}$：

$$r_i=\frac{e_i}{\sqrt{MS_{Res}(1-h_{ii})}}$$

如果某个点的学生化残差绝对值特别大，比如 $|r_i|>3$，它可能是 outlier。

R 中可以这样看：

rstandard(fit)
rstudent(fit)
hatvalues(fit)
cooks.distance(fit)

5.2 残差图怎么读

如果 residuals vs fitted 图是一条随机水平带，通常说明模型没有明显结构性问题。

如果图像像漏斗，残差随着 fitted value 增大而变大，说明可能有异方差。常见处理包括：

• 对响应变量做 log transform。
• 使用 weighted least squares。
• 换一个更合理的模型形式。

如果图像出现曲线，说明线性关系不够，可能需要：

• 加入二次项或交互项。
• 对变量做 transformation。
• 使用 GAM 等更灵活的模型。

6. 变量选择：不是变量越多越好

多元回归最容易犯的错是“把所有变量都扔进去”。变量太多会带来几个问题：

• 模型解释变差。
• 标准误变大。
• 多重共线性变严重。
• 训练集拟合看起来更好，但泛化更差。

变量选择的目标不是找到唯一正确的模型，而是在解释力、稳定性和简洁性之间取得平衡。

常见准则包括：

准则	思想
adjusted $R^2$	惩罚无意义增加变量
AIC	在拟合优度和模型复杂度之间折中
BIC	比 AIC 更强地惩罚复杂模型
Mallows’ $C_p$	比较候选模型的偏差和方差
cross-validation	直接看预测泛化表现

R 中可以用逐步回归：

full_model <- lm(MEDV ~ ., data = housing)
null_model <- lm(MEDV ~ 1, data = housing)

step_model <- step(
  null_model,
  scope = list(lower = null_model, upper = full_model),
  direction = "both"
)

summary(step_model)

但逐步回归不是魔法。它依赖样本，容易不稳定。更好的做法是结合：

领域理解 + 相关性分析 + VIF + 信息准则 + 交叉验证 + 残差诊断

7. 多重共线性：为什么变量相关会伤害模型

多重共线性指 regressors 之间存在近似线性依赖。直观说，就是有些变量在传递相似的信息。

在线性代数上，如果 $X$ 的列高度相关，那么 $X^TX$ 接近奇异，$(X^TX)^{-1}$ 会变得不稳定。结果是：

• 系数估计对数据扰动很敏感。
• 标准误变大。
• 单个变量的 t 检验不显著，但整体模型可能显著。
• 系数符号可能反常。

常用检查方式是 VIF：

$$VIF_j=\frac{1}{1-R_j^2}$$

其中 $R_j^2$ 是用其他自变量回归 $x_j$ 得到的 $R^2$。如果 $x_j$ 能被其他变量很好解释，那么 $R_j^2$ 接近 1，VIF 就会很大。

R 中：

library(car)
fit <- lm(MEDV ~ LSTAT + B + PTRATIO + TAX + RAD + DIS + AGE +
            RM + NOX + CHAS + INDUS + ZN + CRIM,
          data = housing)

vif(fit)

经验上，VIF 超过 5 或 10 就值得警惕，但不要机械套阈值。真正要问的是：这些变量是否在解释同一个现象？是否需要删变量、合成变量、标准化变量，或者改用 ridge regression？

8. Ridge Regression：用一点偏差换稳定性

普通最小二乘最小化：

$$\sum_{i=1}^n (y_i-\hat y_i)^2$$

Ridge regression 在此基础上加入 L2 penalty：

$$\sum_{i=1}^n (y_i-\hat y_i)^2+\lambda\sum_{j=1}^p\beta_j^2$$

它的效果是把系数往 0 收缩，但通常不会直接变成 0。为什么这有用？因为在多重共线性严重时，OLS 系数可能非常不稳定。Ridge 牺牲一点无偏性，换来更低的方差和更稳定的预测。

R 中可以用 MASS::lm.ridge：

library(MASS)

ridge_model <- lm.ridge(
  MEDV ~ LSTAT + B + PTRATIO + TAX + RAD + DIS + AGE +
    RM + NOX + CHAS + INDUS + ZN + CRIM,
  data = housing,
  lambda = seq(0, 1, by = 0.01)
)

matplot(ridge_model$lambda, t(coef(ridge_model)), type = "l",
        xlab = "lambda", ylab = "coefficients")

也可以用 glmnet 做更标准的交叉验证：

library(glmnet)

x <- model.matrix(MEDV ~ . - 1, data = housing)
y <- housing$MEDV

cv_ridge <- cv.glmnet(x, y, alpha = 0)
plot(cv_ridge)
coef(cv_ridge, s = "lambda.min")

其中 alpha = 0 表示 ridge，alpha = 1 表示 lasso。Ridge 更适合处理共线性，Lasso 更适合自动做稀疏变量选择。

9. 随机效应模型：当数据天然分组

普通线性回归默认所有样本相互独立。但很多数据不是这样。比如：

• 学生嵌套在学校里。
• 病人嵌套在医院里。
• 房屋嵌套在区域里。
• 多次测量嵌套在同一个个体里。

如果忽略这种分组结构，模型会低估不确定性，也可能把 group-level variation 错当成普通误差。

线性混合模型可以写成：

$$y_i=X_i\beta+Z_ib_i+\varepsilon_i$$

其中：

• $X_i\beta$ 是 fixed effects，表示所有组共享的平均规律。
• $Z_ib_i$ 是 random effects，表示第 $i$ 个组自己的偏移。
• $\varepsilon_i$ 是组内误差。

最简单的 random intercept model 是：

$$y_{ij}=\beta_0+\beta_1x_{ij}+b_{0i}+\varepsilon_{ij}$$

这里 $b_{0i}$ 表示第 $i$ 个组的截距偏移。它回答的问题是：不同组是否有不同 baseline？

R 中用 lme4：

library(lme4)

mix1 <- lmer(MEDV ~ RM + LSTAT + NOX + PTRATIO + (1 | CHAS),
             data = housing)

summary(mix1)

如果怀疑某个变量的斜率也随组变化，可以写 random slope：

mix2 <- lmer(MEDV ~ RM + LSTAT + NOX + PTRATIO + (1 + DIS | CHAS),
             data = housing)

这里 (1 + DIS | CHAS) 表示：不同 CHAS 组不仅可以有不同截距，也可以有不同的 DIS 斜率。

9.1 ICC：组间差异占多少

Intra-class correlation, ICC, 衡量响应变量有多少比例的变异来自组间差异。对于 random intercept model：

$$\rho=\frac{\psi^2}{\psi^2+\sigma^2}$$

其中 $\psi^2$ 是随机截距方差，$\sigma^2$ 是残差方差。

如果 ICC 很高，说明同一组内样本更相似，分组结构很重要。如果 ICC 很低，说明 random effect 可能没有太大必要。

10. Boston Housing 案例：从数据到模型

Boston Housing 数据有 506 条记录，响应变量是：

• MEDV: median value of owner-occupied homes，以千美元为单位。

常见解释变量包括：

变量	含义
RM	平均房间数
LSTAT	低收入人口比例
TAX	房产税率
CHAS	是否靠近 Charles River
CRIM	城镇犯罪率
ZN	大面积住宅用地比例
NOX	一氧化氮浓度
PTRATIO	学生教师比例
DIS	到就业中心的加权距离
RAD	到高速公路的便利程度
AGE	1940 年前建成的自住房比例
INDUS	非零售商业用地比例
B	原数据中的人口统计变量

10.1 读入数据并做初步检查

housing <- read.csv("boston-housing.csv")

dim(housing)
head(housing)
summary(housing)

第一步不要急着建模。先看：

• 有没有缺失值。
• 响应变量分布是否偏态。
• 自变量是否量纲差异很大。
• 是否有明显异常点。

colSums(is.na(housing))
hist(housing$MEDV, breaks = 30, main = "MEDV distribution")

10.2 相关性与 VIF

先看相关性矩阵：

library(ggplot2)
library(reshape2)

cor_matrix <- cor(housing)
cor_melt <- melt(cor_matrix)

ggplot(cor_melt, aes(Var1, Var2, fill = value)) +
  geom_tile() +
  scale_fill_gradient2(low = "blue", mid = "white", high = "red", midpoint = 0) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Correlation Heatmap", x = "", y = "")

再看 VIF：

library(car)

full_model <- lm(MEDV ~ ., data = housing)
vif(full_model)

在项目资料中，TAX 和 RAD 的 VIF 较高，说明它们和其他变量存在明显共线性。这个现象很合理：交通便利程度、税率、土地用途、污染程度这些变量往往不是独立变化的。

10.3 单变量探索

在建多元模型之前，可以先看几个变量和 MEDV 的关系：

plot(housing$RM, housing$MEDV,
     xlab = "Average number of rooms",
     ylab = "Median home value",
     main = "MEDV vs RM")

plot(housing$LSTAT, housing$MEDV,
     xlab = "Lower status population percentage",
     ylab = "Median home value",
     main = "MEDV vs LSTAT")

boxplot(MEDV ~ CHAS, data = housing,
        xlab = "Charles River dummy",
        ylab = "Median home value",
        main = "MEDV by CHAS")

通常会看到：

• RM 和 MEDV 大致正相关。
• LSTAT 和 MEDV 明显负相关，而且可能非线性。
• CHAS=1 的区域房价分布可能更高，但样本量较少，需要谨慎解释。
• DIS、NOX、PTRATIO 等变量的影响需要放在多元模型里看。

10.4 一个基础多元线性模型

可以先从可解释性较强的变量开始：

fit_lm <- lm(MEDV ~ RM + LSTAT + NOX + PTRATIO + DIS + CHAS,
             data = housing)

summary(fit_lm)
par(mfrow = c(2, 2))
plot(fit_lm)

解释时重点看：

• RM 是否显著为正。
• LSTAT 是否显著为负。
• PTRATIO 是否为负。
• NOX 是否在控制其他变量后仍显著。
• 残差图是否出现非线性或异方差。

如果 Q-Q plot 两端偏离明显，说明尾部样本拟合不好。如果 residuals vs fitted 有曲线，说明模型可能需要非线性项，例如：

fit_poly <- lm(MEDV ~ RM + LSTAT + I(LSTAT^2) + NOX + PTRATIO + DIS + CHAS,
               data = housing)

anova(fit_lm, fit_poly)

I(LSTAT^2) 表示加入 LSTAT 的二次项。

10.5 Ridge 辅助变量选择

如果变量很多而且共线性明显，可以用 ridge trace 看系数稳定过程：

library(MASS)

ridge_fit <- lm.ridge(MEDV ~ ., data = housing,
                      lambda = seq(0, 1, by = 0.01))

matplot(ridge_fit$lambda, t(coef(ridge_fit)), type = "l",
        xlab = "lambda",
        ylab = "ridge coefficients",
        main = "Ridge trace")

项目资料里通过 ridge 思路筛出了 CHAS、NOX、RM、DIS、PTRATIO、LSTAT 这几个变量。这个结果可以理解成：这些变量在房价建模中保留了较强解释力，同时比全模型更简洁。

但需要注意：用系数大小做变量筛选前最好标准化变量，否则不同量纲会影响比较。

更稳妥的版本：

library(glmnet)

x <- model.matrix(MEDV ~ . - 1, data = housing)
y <- housing$MEDV

cv_lasso <- cv.glmnet(x, y, alpha = 1)
coef(cv_lasso, s = "lambda.min")

Lasso 可以直接把一些系数压到 0，更适合做变量选择；Ridge 更适合缓解共线性。

10.6 混合效应模型是否合适

项目资料尝试把 CHAS 作为分组变量：

library(lme4)

mix1 <- lmer(MEDV ~ RM + LSTAT + NOX + PTRATIO + (1 | CHAS),
             data = housing)

mix2 <- lmer(MEDV ~ RM + LSTAT + NOX + PTRATIO + (1 + DIS | CHAS),
             data = housing)

anova(mix1, mix2)

这里的想法是：靠近 Charles River 和不靠近 Charles River 的区域，可能不仅平均房价不同，而且某些变量的影响斜率也不同。

但是这里有一个很重要的统计建模提醒：CHAS 只有两个组，作为 random effect 的依据并不强。Random effect 更适合“组数较多，每组有多个样本”的结构，例如很多学校、很多医院、很多地区。如果只有 0/1 两组，通常把 CHAS 当作 fixed effect 更自然。

所以这部分更适合当作 mixed model 语法和思想练习，而不是说 Boston Housing 一定需要混合效应模型。

11. 一个更可靠的建模流程

如果我要把这个项目改成一份更成熟的数据科学作品，我会按下面流程做：

1. 数据理解：解释变量含义，检查缺失、分布、异常点。
2. EDA：画 MEDV 分布、相关性热力图、关键变量散点图。
3. 基础模型：先做一个可解释的 OLS baseline。
4. 诊断：检查残差、杠杆点、Cook’s distance、异方差和非线性。
5. 改进：加入 transformation、二次项或交互项。
6. 共线性处理：检查 VIF，用 ridge/lasso 做稳健比较。
7. 验证：划分 train/test 或做 cross-validation。
8. 解释：用系数、标准化系数、误差指标和图形共同说明结论。

一个简化但比较完整的 R 模板：

library(tidyverse)
library(car)
library(glmnet)

housing <- read.csv("boston-housing.csv")

set.seed(42)
idx <- sample(seq_len(nrow(housing)), size = 0.8 * nrow(housing))
train <- housing[idx, ]
test <- housing[-idx, ]

baseline <- lm(MEDV ~ RM + LSTAT + NOX + PTRATIO + DIS + CHAS,
               data = train)

summary(baseline)
vif(baseline)

pred <- predict(baseline, newdata = test)
rmse <- sqrt(mean((test$MEDV - pred)^2))
mae <- mean(abs(test$MEDV - pred))

c(RMSE = rmse, MAE = mae)

如果进一步加正则化：

x_train <- model.matrix(MEDV ~ . - 1, data = train)
y_train <- train$MEDV
x_test <- model.matrix(MEDV ~ . - 1, data = test)

cv_ridge <- cv.glmnet(x_train, y_train, alpha = 0)
ridge_pred <- predict(cv_ridge, newx = x_test, s = "lambda.min")

sqrt(mean((test$MEDV - ridge_pred)^2))

这样你就不只是“跑了一个回归”，而是在比较模型的解释性和预测能力。

12. 学完 AMA3602 应该带走什么

这门课最值得带走的不是某个公式，而是一种建模习惯：

模型不是结果，模型是一个需要被检查、解释和迭代的假设。

简单线性回归教你理解关系，多元线性回归教你控制变量，残差诊断教你怀疑模型，变量选择教你控制复杂度，ridge regression 教你处理不稳定，mixed model 教你面对分组数据。

如果以后做机器学习、推荐系统、金融建模、A/B test 或实验分析，这些思想都会反复出现。深度学习模型可以更复杂，但“先建立 baseline、检查误差、理解变量、验证泛化”的习惯，依然是从线性模型这里长出来的。