R语言线性回归

パーフェクト

212人浏览 · 2026-06-09 17:47:35

パーフェクト · 2026-06-09 17:47:35 发布

一元线性回归模型

一元线性回归模型为：

$Y=\beta_0+\beta_1X+\varepsilon$

其中： $\beta_0$ 是截距项； $\beta_1$ 是回归系数； $\varepsilon$ 是误差项。

对于建模，我们需要找到一条直线，使所有样本点到这条直线的残差平方和最小，也就是到这条直线的距离之和最小。也就是：

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

等价于：

$SSE=\sum_{i=1}^n(y_i-\beta_0-\beta_1x_i)^2$

通过对 $\beta_0$ 和 $\beta_1$ 分别求偏导，并令偏导等于0，得到最小二乘估计。

$\hat{\beta}_1= \frac{\sum_{i=1}^{n}(x_i-\bar{x})(y_i-\bar{y})} {\sum_{i=1}^{n}(x_i-\bar{x})^2}$

$\hat{\beta}_0=\bar{y}-\hat{\beta}_1\bar{x}$

# ---------------------------------------------------------
# 1. 构造模拟数据
# ---------------------------------------------------------

set.seed(123)

# 样本量
n <- 100

# 自变量 x
x <- 1:n

# 随机误差项 epsilon ~ N(0, 1)
epsilon <- rnorm(n, mean = 0, sd = 1)

# 因变量 y
# 真实模型：y = 10 + 3*x + epsilon
y <- 10 + 3 * x + epsilon


# ---------------------------------------------------------
# 2. 手动计算最小二乘估计
# ---------------------------------------------------------

# 样本均值
x_bar <- sum(x) / n
y_bar <- sum(y) / n

# 计算 Sxx 和 Sxy
Sxx <- sum((x - x_bar)^2)
Sxy <- sum((x - x_bar) * (y - y_bar))

# 回归系数 beta1
beta1_hat <- Sxy / Sxx

# 截距项 beta0
beta0_hat <- y_bar - beta1_hat * x_bar

# 拟合值
y_hat_manual <- beta0_hat + beta1_hat * x

# 残差
residuals_manual <- y - y_hat_manual

# 残差平方和 SSE
SSE <- sum(residuals_manual^2)

# 总平方和 SST
SST <- sum((y - y_bar)^2)

# 回归平方和 SSR
SSR <- sum((y_hat_manual - y_bar)^2)

# 决定系数 R^2
R_squared <- SSR / SST

# 残差方差估计
sigma2_hat <- SSE / (n - 2)

# 回归系数标准误
se_beta1 <- sqrt(sigma2_hat / Sxx)

# 截距标准误
se_beta0 <- sqrt(sigma2_hat * (1 / n + x_bar^2 / Sxx))

# t 统计量
t_beta1 <- beta1_hat / se_beta1
t_beta0 <- beta0_hat / se_beta0


# ---------------------------------------------------------
# 3. 使用 R 自带 lm() 进行回归
# ---------------------------------------------------------

model_lm <- lm(y ~ x)

lm_result <- summary(model_lm)

# lm() 的拟合值
y_hat_lm <- fitted(model_lm)


# ---------------------------------------------------------
# 4. 输出手动计算结果
# ---------------------------------------------------------

cat("\n========== 手动最小二乘估计结果 ==========\n")

cat("样本量 n =", n, "\n")
cat("x 的均值 =", x_bar, "\n")
cat("y 的均值 =", y_bar, "\n")

cat("\n手动估计的截距 beta0_hat =", beta0_hat, "\n")
cat("手动估计的回归系数 beta1_hat =", beta1_hat, "\n")

cat("\n手动回归方程为：\n")
cat("y_hat =", beta0_hat, "+", beta1_hat, "* x\n")

cat("\n残差平方和 SSE =", SSE, "\n")
cat("回归平方和 SSR =", SSR, "\n")
cat("总平方和 SST =", SST, "\n")
cat("决定系数 R^2 =", R_squared, "\n")

cat("\n截距标准误 se_beta0 =", se_beta0, "\n")
cat("回归系数标准误 se_beta1 =", se_beta1, "\n")

cat("\n截距 t 统计量 =", t_beta0, "\n")
cat("回归系数 t 统计量 =", t_beta1, "\n")


# ---------------------------------------------------------
# 5. 输出 lm() 的结果
# ---------------------------------------------------------

cat("\n========== lm() 回归结果 ==========\n")
print(lm_result)

cat("\n========== lm() 提取的回归系数 ==========\n")
print(coef(model_lm))


# ---------------------------------------------------------
# 6. 比较手动结果和 lm() 结果
# ---------------------------------------------------------

cat("\n========== 手动结果与 lm() 结果比较 ==========\n")

cat("手动截距 beta0_hat =", beta0_hat, "\n")
cat("lm() 截距 =", coef(model_lm)[1], "\n")

cat("\n手动回归系数 beta1_hat =", beta1_hat, "\n")
cat("lm() 回归系数 =", coef(model_lm)[2], "\n")

cat("\n截距差异 =", beta0_hat - coef(model_lm)[1], "\n")
cat("回归系数差异 =", beta1_hat - coef(model_lm)[2], "\n")


# ---------------------------------------------------------
# 7. 绘制散点图和两条回归直线
# ---------------------------------------------------------

plot(
  x,
  y,
  main = "一元线性回归：手动OLS与lm()结果比较",
  xlab = "x",
  ylab = "y",
  pch = 19,
  col = "black"
)

# 手动OLS回归线：红色
abline(
  a = beta0_hat,
  b = beta1_hat,
  col = "red",
  lwd = 3
)

# lm()回归线：蓝色虚线
abline(
  model_lm,
  col = "blue",
  lwd = 2,
  lty = 2
)

# 添加图例
legend(
  "topleft",
  legend = c("样本点", "手动OLS回归线", "lm()回归线"),
  col = c("black", "red", "blue"),
  pch = c(19, NA, NA),
  lty = c(NA, 1, 2),
  lwd = c(NA, 3, 2),
  bty = "n"
)

手动编写和利用R包的效果几乎一致。

当维度增加，我们的一元线性回归就变成多元线性回归了。

多元线性回归模型

模型的一般形式为：

$Y=\beta_0+\beta_1X_1+\beta_2X_2+\cdots+\beta_pX_p+\varepsilon$

其中： $\beta_0$ 是截距项， $\beta_1,\beta_2,......,\beta_p$ 表示各个自变量的回归系数； $X_1,X_2,......,X_p$ 表示多个解释变量； $\varepsilon$ 表示随机误差项。

矩阵形式：

$\mathbf{Y}=\mathbf{X}\boldsymbol{\beta}+\boldsymbol{\varepsilon}$

其中：

$\mathbf{Y}= \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{pmatrix}$

$\mathbf{X}= \begin{pmatrix} 1 & x_{11} & x_{12} & \cdots & x_{1p}\\ 1 & x_{21} & x_{22} & \cdots & x_{2p}\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{pmatrix}$

其建模的核心思想与一元回归模型的建模核心思想是一样的。

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

矩阵形式：

$SSE=(\mathbf{Y}-\mathbf{X}\boldsymbol{\beta})^{T}(\mathbf{Y}-\mathbf{X}\boldsymbol{\beta})$

通过最小化 SSE，可以得到参数估计公式：

$\hat{\boldsymbol{\beta}}= (\mathbf{X}^{\mathrm{T}}\mathbf{X})^{-1} \mathbf{X}^{\mathrm{T}}\mathbf{Y}$

多元线性回归通常有几个重要假设：

第一，线性关系假设。因变量和自变量之间满足线性关系。

第二，误差项均值为 0：

$E({\varepsilon})=0$

第三，同方差性。误差项的方差保持不变：

$VAR(\varepsilon)=\sigma^2$

第四，误差项相互独立。

第五，不存在严重多重共线性。也就是说，自变量之间不能高度重复或高度相关。

第六，误差项服从正态分布。

第六条假设主要用于后续的 t 检验、F 检验和置信区间推断。

多元回归中常见两个检验。

1. 单个回归系数的 t 检验

检验某个自变量是否对因变量有显著影响。

假设检验为：

$H_0:\beta_j=0\quad VS \quad H_1:\beta_j\neq0$

t 统计量为：

$t=\frac{\hat{\beta}_j}{SE(\hat{\beta}_j)}$

如果 p 值小于 0.05，说明该变量对因变量具有显著影响。

2. 模型整体的 F 检验

检验整个回归模型是否显著。

$H_0:\beta_1=\beta_2=......=\beta_p=0\quad VS \quad H_1:$ 至少有一个 $\beta_j\neq0$

如果 F 检验显著，说明模型整体有解释能力。

模型评价指标

常用的是决定系数 $R^2$ ：

$R^2=1-\frac{SSE}{SST}$

它表示自变量对因变量变异的解释程度。

但是在多元回归中，加入更多变量通常会让 $R^2$ 增大，所以更常用调整后的 $R^2$ ：

$R_a^2=1-\frac{SSE/(n-p-1)}{SST/(n-1)}$

调整后的 $R^2$ 会考虑自变量个数，更适合比较不同复杂度的模型。

总结就是：多元线性回归模型是在多个自变量共同作用下解释或预测因变量的方法，其核心是利用最小二乘法估计参数，并通过 t 检验、F 检验和 $R^2$ 等指标判断变量影响和模型效果。

案列

$Y=5+2X_1-1.5X_2+3X_3+0X_4+1.2X_5-2X_6+0X_7+1.5X_9-X_{10}+\varepsilon$

# ---------------------------------------------------------
# 1. 构造模拟数据
# ---------------------------------------------------------

set.seed(123)

# 样本量
n <- 200

# 自变量个数
p <- 10

# 生成 n × p 的自变量矩阵
X_raw <- matrix(rnorm(n * p, mean = 0, sd = 1), nrow = n, ncol = p)

# 给自变量命名
colnames(X_raw) <- paste0("X", 1:p)

# 真实参数
beta0_true <- 5

beta_true <- c(
  2.0,    # X1
  -1.5,   # X2
  3.0,    # X3
  0.0,    # X4
  1.2,    # X5
  -2.0,   # X6
  0.0,    # X7
  0.8,    # X8
  1.5,    # X9
  -1.0    # X10
)

# 随机误差项
epsilon <- rnorm(n, mean = 0, sd = 1)

# 生成因变量 Y
Y <- beta0_true + X_raw %*% beta_true + epsilon

# 转成向量
Y <- as.vector(Y)

# 整理成数据框，方便后面与 lm() 比较
data <- data.frame(Y = Y, X_raw)

cat("\n========== 模拟数据前 6 行 ==========\n")
print(head(data))


# ---------------------------------------------------------
# 2. 手动编写多元线性回归函数
# ---------------------------------------------------------

manual_multiple_lm <- function(X_raw, Y) {
  
  # 样本量
  n <- nrow(X_raw)
  
  # 自变量个数
  p <- ncol(X_raw)
  
  # 加入截距项，第一列全为 1
  X <- cbind(1, X_raw)
  
  colnames(X) <- c("Intercept", colnames(X_raw))
  
  # 转成矩阵形式
  Y <- matrix(Y, ncol = 1)
  
  # -------------------------------------------------------
  # OLS 核心公式：
  # beta_hat = (X'X)^(-1) X'Y
  # -------------------------------------------------------
  
  XtX <- t(X) %*% X
  XtY <- t(X) %*% Y
  
  beta_hat <- solve(XtX) %*% XtY
  
  # 拟合值
  Y_hat <- X %*% beta_hat
  
  # 残差
  residuals <- Y - Y_hat
  
  # Y 的均值
  Y_bar <- mean(Y)
  
  # 残差平方和 SSE
  SSE <- sum(residuals^2)
  
  # 总平方和 SST
  SST <- sum((Y - Y_bar)^2)
  
  # 回归平方和 SSR
  SSR <- sum((Y_hat - Y_bar)^2)
  
  # 残差方差估计
  sigma2_hat <- SSE / (n - p - 1)
  
  # 决定系数 R^2
  R_squared <- 1 - SSE / SST
  
  # 调整后的 R^2
  adj_R_squared <- 1 - (SSE / (n - p - 1)) / (SST / (n - 1))
  
  # 系数协方差矩阵
  cov_beta <- sigma2_hat * solve(XtX)
  
  # 回归系数标准误
  se_beta <- sqrt(diag(cov_beta))
  
  # t 统计量
  t_value <- as.vector(beta_hat) / se_beta
  
  # p 值，这里使用 R 自带 pt()，不是外部包
  p_value <- 2 * (1 - pt(abs(t_value), df = n - p - 1))
  
  # 输出结果表
  coef_table <- data.frame(
    variable = colnames(X),
    estimate = as.vector(beta_hat),
    std_error = se_beta,
    t_value = t_value,
    p_value = p_value
  )
  
  result <- list(
    n = n,
    p = p,
    X = X,
    Y = Y,
    beta_hat = beta_hat,
    fitted_values = Y_hat,
    residuals = residuals,
    SSE = SSE,
    SSR = SSR,
    SST = SST,
    sigma2_hat = sigma2_hat,
    R_squared = R_squared,
    adj_R_squared = adj_R_squared,
    coef_table = coef_table
  )
  
  return(result)
}


# ---------------------------------------------------------
# 3. 运行手动多元线性回归
# ---------------------------------------------------------

manual_result <- manual_multiple_lm(X_raw, Y)


# ---------------------------------------------------------
# 4. 输出手动回归结果
# ---------------------------------------------------------

cat("\n========== 手动多元线性回归结果 ==========\n")

cat("样本量 n =", manual_result$n, "\n")
cat("自变量个数 p =", manual_result$p, "\n")

cat("\n残差平方和 SSE =", manual_result$SSE, "\n")
cat("回归平方和 SSR =", manual_result$SSR, "\n")
cat("总平方和 SST =", manual_result$SST, "\n")

cat("\n残差方差估计 sigma^2 =", manual_result$sigma2_hat, "\n")
cat("R^2 =", manual_result$R_squared, "\n")
cat("调整后的 R^2 =", manual_result$adj_R_squared, "\n")

cat("\n========== 手动估计的回归系数表 ==========\n")
print(manual_result$coef_table)


# ---------------------------------------------------------
# 5. 使用 R 自带 lm() 进行比较
# ---------------------------------------------------------

lm_model <- lm(Y ~ ., data = data)

cat("\n========== lm() 回归结果 ==========\n")
print(summary(lm_model))


# ---------------------------------------------------------
# 6. 比较手动结果和 lm() 的回归系数
# ---------------------------------------------------------

manual_coef <- as.vector(manual_result$beta_hat)
lm_coef <- coef(lm_model)

comparison <- data.frame(
  variable = names(lm_coef),
  manual_estimate = manual_coef,
  lm_estimate = as.vector(lm_coef),
  difference = manual_coef - as.vector(lm_coef)
)

cat("\n========== 手动 OLS 与 lm() 系数比较 ==========\n")
print(comparison)


# ---------------------------------------------------------
# 7. 画真实值与拟合值对比图
# ---------------------------------------------------------

plot(
  Y,
  manual_result$fitted_values,
  main = "多元线性回归：真实值与拟合值对比",
  xlab = "真实 Y",
  ylab = "拟合 Y",
  pch = 19,
  col = "blue"
)

abline(
  a = 0,
  b = 1,
  col = "red",
  lwd = 2
)

legend(
  "topleft",
  legend = c("拟合点", "理想拟合线 y=x"),
  col = c("blue", "red"),
  pch = c(19, NA),
  lty = c(NA, 1),
  lwd = c(NA, 2),
  bty = "n"
)