回归模型的诊断，主要是在模型建立之后检查：主要检查这个回归模型是否可靠？参数估计是否稳定？显著性检验和预测结果是否可信？它不是重新估计参数，而是检查模型是否满足线性回归的基本假设。

为什么要做回归诊断？

通过最小二乘法得到

但是，得到以后，并不代表模型一定可靠。因为线性回归的 t 检验、F 检验、置信区间、预测结果都依赖一些基本假设。（线性关系、误差均值为0、残差近似正态......）

如果这些假设不满足，可能出现：

模型显著但解释错误；
回归系数不稳定；
p 值不可信；
预测误差很大；
少数异常点严重影响模型结果。

所以回归诊断的核心目的就是检查模型是否存在这些问题。

回归诊断主要检查什么？

回归诊断一般检查五个方面：

线性关系→残差正态性→同方差性→异常点和强影响点→多重共线性

线性关系诊断

线性回归的基本假设是：

也就是说，因变量 Y 与自变量之间应该近似满足线性关系。

如果真实关系是非线性的，例如：

但我们却强行使用线性模型：

就会导致模型拟合不充分。

常用诊断图是：

残差—拟合值图 Residuals vs Fitted

如果模型合理，残差应该围绕 0 随机分布，没有明显规律。

如果残差图出现弧形、U 形、波浪形，就说明可能存在非线性关系。

残差定义为：

残差正态性诊断

线性回归中，通常假设误差项服从正态分布：

这个假设主要影响 t 检验、F 检验和置信区间的可靠性。

如果样本量较大，正态性要求可以适当放宽；但如果样本量较小，残差严重偏离正态分布，就会影响显著性检验。

常用诊断图是：

Normal Q-Q 图

如果残差近似正态分布，点应该大致落在一条直线上。

如果两端明显偏离直线，说明残差可能存在偏态或厚尾问题。

常用残差形式是标准化残差：

其中，是杠杆值， 是残差标准差。

同方差性诊断

线性回归要求误差项具有相同方差：

这叫同方差性。

如果误差方差随着拟合值或某些自变量变化而变化，就叫异方差。

例如，收入越高，消费波动越大；企业规模越大，利润波动越大。这些情况都容易出现异方差。

常用诊断图是：

Scale-Location 图

也可以看残差—拟合值图。

如果残差的波动范围大致稳定，说明同方差性较好。

如果残差呈现“漏斗形”，例如左边窄、右边宽，就说明可能存在异方差。

异方差会导致：

回归系数估计仍可能无偏；
但标准误可能错误；
进一步导致 t 检验和 p 值不可靠。

异常点诊断

异常点是指某些样本的残差特别大，也就是实际值和预测值相差很远。

判断异常点常用标准化残差：

一般经验上：

说明该点可能异常；

说明该点高度异常，需要重点检查。

异常点可能来自：

数据录入错误；
极端样本；
模型遗漏重要变量；
模型形式不正确。

但是要注意：异常点不一定必须删除，需要结合实际背景判断。

多重共线性诊断

多元回归中，自变量之间如果高度相关，就会产生多重共线性。

例如：

这两个变量可能高度相关，因为学习时间越长，刷题数量往往也越多。

多重共线性会导致：

回归系数不稳定；
标准误变大；
变量本来有影响，但 t 检验不显著；
系数方向可能反常。

常用诊断指标是 VIF，方差膨胀因子：

其中，是用第 个自变量作为因变量，用其他自变量进行回归得到的决定系数。

一般来说：

说明可能存在较强共线性；

说明共线性问题比较严重。

高杠杆点诊断

杠杆值用于判断某个样本在自变量空间中是否特殊。

Hat 矩阵为：

第 个样本的杠杆值为：

它是 Hat 矩阵对角线上的元素。

如果某个样本的自变量取值非常特殊，比如远离大多数样本，它就可能有较高杠杆值。

常用经验阈值是：

或者

其中，是自变量个数， 是样本量。

高杠杆点本身不一定是坏点。如果它的残差很小，说明它虽然特殊，但仍然符合模型规律；如果它同时残差很大，就可能严重影响回归结果。

如果诊断发现问题怎么办？

如果发现非线性关系，可以加入二次项、交互项，或者进行变量变换。

如果发现异方差，可以考虑对取对数，或者使用加权最小二乘法。

如果发现残差严重非正态，可以检查异常值，或者考虑稳健回归。

如果发现强影响点，要回到原始数据检查是否存在录入错误。如果不是错误，就需要结合实际背景决定是否保留。

如果发现多重共线性，可以删除高度相关变量，或者使用岭回归、LASSO 回归。

总结

回归诊断的核心不是看模型有没有显著，而是判断：

模型结果是否可信，参数是否稳定，残差是否满足基本假设，是否存在异常样本严重影响模型。

可以概括为：

建模→估计参数→显著性检验→回归诊断→模型修正

所以，回归诊断是从“会建模”走向“会判断模型质量”的关键步骤。

案例演示

# =========================================================
# 多元线性回归模型诊断
# 包括：
# 1. 线性关系诊断
# 2. 残差正态性检验
# 3. 同方差性检验
# 4. 异常点诊断
# 5. 多重共线性诊断
# 6. 高杠杆点诊断
# =========================================================


# =========================================================
# 1. 构造十维多元线性回归模拟数据
# =========================================================

set.seed(123)

n <- 200      # 样本量
p <- 10       # 自变量个数

# 生成 10 个自变量
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 <- beta0_true + X_raw %*% beta_true + epsilon
Y <- as.vector(Y)


# =========================================================
# 2. 手动建立多元线性回归模型
# =========================================================

# 加入截距项
X <- cbind(1, X_raw)
colnames(X) <- c("Intercept", colnames(X_raw))

Y_matrix <- matrix(Y, ncol = 1)

# OLS估计：beta_hat = (X'X)^(-1)X'Y
XtX <- t(X) %*% X
XtY <- t(X) %*% Y_matrix

beta_hat <- solve(XtX) %*% XtY

# 拟合值
Y_hat <- X %*% beta_hat

# 残差
e <- Y_matrix - Y_hat

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

# 总平方和
SST <- sum((Y_matrix - mean(Y_matrix))^2)

# 回归平方和
SSR <- SST - SSE

# 残差自由度
df_residual <- n - p - 1

# 残差方差估计
sigma2_hat <- SSE / df_residual

# 残差标准差
sigma_hat <- sqrt(sigma2_hat)

# R方
R_squared <- 1 - SSE / SST

# 调整后的R方
adj_R_squared <- 1 - (SSE / df_residual) / (SST / (n - 1))


# =========================================================
# 3. 回归诊断基础指标
# =========================================================

# Hat矩阵 H = X(X'X)^(-1)X'
H <- X %*% solve(t(X) %*% X) %*% t(X)

# 杠杆值 h_ii
h <- diag(H)

# 标准化残差
standardized_residuals <- as.vector(e) / (sigma_hat * sqrt(1 - h))

# Cook距离，这里作为辅助诊断指标
cooks_distance <- (as.vector(e)^2 / ((p + 1) * sigma2_hat)) *
  (h / (1 - h)^2)

# 常用阈值
outlier_threshold_2 <- 2
outlier_threshold_3 <- 3
leverage_threshold <- 2 * (p + 1) / n
cook_threshold <- 4 / n


# =========================================================
# 4. 手动计算卡方分布函数
# 用于 Jarque-Bera 和 Breusch-Pagan 检验的 p 值
# 不调用 pchisq()
# =========================================================

chisq_pdf <- function(x, df) {
  
  value <- rep(0, length(x))
  idx <- which(x > 0)
  
  if (length(idx) > 0) {
    xp <- x[idx]
    
    log_value <- (df / 2 - 1) * log(xp) -
      xp / 2 -
      (df / 2) * log(2) -
      lgamma(df / 2)
    
    value[idx] <- exp(log_value)
  }
  
  return(value)
}


simpson_integral_chisq <- function(f, a, b, df, grid_n = 8000) {
  
  if (a == b) {
    return(0)
  }
  
  if (grid_n %% 2 == 1) {
    grid_n <- grid_n + 1
  }
  
  h_grid <- (b - a) / grid_n
  x_grid <- seq(a, b, length.out = grid_n + 1)
  y_grid <- f(x_grid, df)
  
  result <- h_grid / 3 * (
    y_grid[1] +
      y_grid[grid_n + 1] +
      4 * sum(y_grid[seq(2, grid_n, by = 2)]) +
      2 * sum(y_grid[seq(3, grid_n - 1, by = 2)])
  )
  
  return(result)
}


chisq_cdf <- function(q, df) {
  
  if (q <= 0) {
    return(0)
  }
  
  prob <- simpson_integral_chisq(chisq_pdf, 0, q, df)
  
  if (prob < 0) {
    prob <- 0
  }
  
  if (prob > 1) {
    prob <- 1
  }
  
  return(prob)
}


# =========================================================
# 5. 残差正态性检验：Jarque-Bera 检验
# H0：残差服从正态分布
# H1：残差不服从正态分布
# =========================================================

resid_vec <- as.vector(e)

resid_centered <- resid_vec - mean(resid_vec)

m2 <- mean(resid_centered^2)
m3 <- mean(resid_centered^3)
m4 <- mean(resid_centered^4)

skewness <- m3 / (m2^(3 / 2))
kurtosis <- m4 / (m2^2)

JB_stat <- n / 6 * (skewness^2 + ((kurtosis - 3)^2) / 4)

JB_p_value <- 1 - chisq_cdf(JB_stat, df = 2)


# =========================================================
# 6. 同方差性检验：Breusch-Pagan 检验
# H0：同方差
# H1：异方差
# 思路：将残差平方 e^2 对所有自变量回归，取辅助回归 R^2
# BP = n * R_aux^2
# =========================================================

manual_R2 <- function(X_input, Y_input) {
  
  X_input <- as.matrix(X_input)
  Y_input <- matrix(Y_input, ncol = 1)
  
  n_input <- nrow(X_input)
  
  X_design <- cbind(1, X_input)
  
  beta_aux <- solve(t(X_design) %*% X_design) %*% t(X_design) %*% Y_input
  
  Y_hat_aux <- X_design %*% beta_aux
  
  SSE_aux <- sum((Y_input - Y_hat_aux)^2)
  SST_aux <- sum((Y_input - mean(Y_input))^2)
  
  R2_aux <- 1 - SSE_aux / SST_aux
  
  return(R2_aux)
}

e_squared <- resid_vec^2

BP_R2 <- manual_R2(X_raw, e_squared)

BP_stat <- n * BP_R2

BP_p_value <- 1 - chisq_cdf(BP_stat, df = p)


# =========================================================
# 7. 多重共线性诊断：手动计算 VIF
# VIF_j = 1 / (1 - R_j^2)
# =========================================================

compute_vif <- function(X_raw) {
  
  X_raw <- as.matrix(X_raw)
  
  p <- ncol(X_raw)
  
  vif_values <- numeric(p)
  
  for (j in 1:p) {
    
    # 第 j 个变量作为因变量
    Y_j <- X_raw[, j]
    
    # 其他变量作为自变量
    X_other <- X_raw[, -j, drop = FALSE]
    
    R2_j <- manual_R2(X_other, Y_j)
    
    vif_values[j] <- 1 / (1 - R2_j)
  }
  
  names(vif_values) <- colnames(X_raw)
  
  return(vif_values)
}

vif_values <- compute_vif(X_raw)


# =========================================================
# 8. 异常点、高杠杆点、强影响点识别
# =========================================================

# 标准化残差绝对值大于 2 的点
outlier_index_2 <- which(abs(standardized_residuals) > 2)

# 标准化残差绝对值大于 3 的点
outlier_index_3 <- which(abs(standardized_residuals) > 3)

# 高杠杆点
high_leverage_index <- which(h > leverage_threshold)

# Cook距离较大的点，辅助观察
high_cook_index <- which(cooks_distance > cook_threshold)


# =========================================================
# 9. 输出模型基本结果
# =========================================================

cat("\n========== 多元线性回归模型基本结果 ==========\n")
cat("样本量 n =", n, "\n")
cat("自变量个数 p =", p, "\n")
cat("残差自由度 =", df_residual, "\n")
cat("残差平方和 SSE =", SSE, "\n")
cat("回归平方和 SSR =", SSR, "\n")
cat("总平方和 SST =", SST, "\n")
cat("残差标准误 sigma_hat =", sigma_hat, "\n")
cat("R^2 =", R_squared, "\n")
cat("调整后的 R^2 =", adj_R_squared, "\n")

cat("\n========== 回归系数估计 ==========\n")
coef_table <- data.frame(
  variable = colnames(X),
  estimate = as.vector(beta_hat)
)
print(coef_table)


# =========================================================
# 10. 输出六类诊断结果
# =========================================================

cat("\n========== 1. 线性关系诊断 ==========\n")
cat("主要查看 Residuals vs Fitted 图。\n")
cat("如果残差围绕 0 随机分布，没有明显曲线或漏斗形，说明线性关系基本合理。\n")

cat("\n========== 2. 残差正态性检验 ==========\n")
cat("Jarque-Bera 统计量 =", JB_stat, "\n")
cat("Jarque-Bera p 值 =", JB_p_value, "\n")
if (JB_p_value > 0.05) {
  cat("结论：不能拒绝残差正态性假设，残差近似服从正态分布。\n")
} else {
  cat("结论：拒绝残差正态性假设，残差可能不服从正态分布。\n")
}

cat("\n========== 3. 同方差性检验 ==========\n")
cat("Breusch-Pagan 统计量 =", BP_stat, "\n")
cat("Breusch-Pagan p 值 =", BP_p_value, "\n")
if (BP_p_value > 0.05) {
  cat("结论：不能拒绝同方差假设，暂未发现明显异方差问题。\n")
} else {
  cat("结论：拒绝同方差假设，模型可能存在异方差问题。\n")
}

cat("\n========== 4. 异常点诊断 ==========\n")
cat("标准化残差阈值 |r_i| > 2 的样本点：\n")
print(outlier_index_2)

cat("标准化残差阈值 |r_i| > 3 的样本点：\n")
print(outlier_index_3)

if (length(outlier_index_3) == 0) {
  cat("结论：未发现高度异常点。\n")
} else {
  cat("结论：存在高度异常点，需要结合原始数据进一步检查。\n")
}

cat("\n========== 5. 多重共线性诊断 ==========\n")
print(vif_values)

if (max(vif_values) < 5) {
  cat("结论：所有 VIF 均小于 5，未发现明显多重共线性问题。\n")
} else if (max(vif_values) < 10) {
  cat("结论：存在一定多重共线性风险，需要关注 VIF 较大的变量。\n")
} else {
  cat("结论：存在严重多重共线性问题，需要考虑删除变量或使用岭回归、LASSO。\n")
}

cat("\n========== 6. 高杠杆点诊断 ==========\n")
cat("高杠杆点阈值 2(p+1)/n =", leverage_threshold, "\n")
cat("最大杠杆值 =", max(h), "\n")
cat("高杠杆点样本编号：\n")
print(high_leverage_index)

if (length(high_leverage_index) == 0) {
  cat("结论：未发现明显高杠杆点。\n")
} else {
  cat("结论：存在高杠杆点，这些样本在自变量空间中较特殊。\n")
}

cat("\n========== Cook 距离辅助诊断 ==========\n")
cat("Cook距离阈值 4/n =", cook_threshold, "\n")
cat("最大 Cook 距离 =", max(cooks_distance), "\n")
cat("Cook距离较大的样本编号：\n")
print(high_cook_index)


# =========================================================
# 11. 绘制六类回归诊断图
# =========================================================

par(mfrow = c(2, 3))

# 图1：线性关系诊断，残差 vs 拟合值
plot(
  as.vector(Y_hat),
  resid_vec,
  main = "1. Residuals vs Fitted",
  xlab = "Fitted values",
  ylab = "Residuals",
  pch = 19,
  col = "blue"
)
abline(h = 0, col = "red", lwd = 2)

# 图2：残差正态性诊断，QQ图
qqnorm(
  standardized_residuals,
  main = "2. Normal Q-Q",
  pch = 19,
  col = "blue"
)
qqline(
  standardized_residuals,
  col = "red",
  lwd = 2
)

# 图3：同方差性诊断，Scale-Location图
plot(
  as.vector(Y_hat),
  sqrt(abs(standardized_residuals)),
  main = "3. Scale-Location",
  xlab = "Fitted values",
  ylab = "Sqrt(|standardized residuals|)",
  pch = 19,
  col = "blue"
)
abline(
  h = mean(sqrt(abs(standardized_residuals))),
  col = "red",
  lwd = 2
)

# 图4：异常点诊断，标准化残差图
plot(
  standardized_residuals,
  main = "4. Standardized Residuals",
  xlab = "Observation index",
  ylab = "Standardized residuals",
  pch = 19,
  col = "blue"
)
abline(h = 0, col = "black", lwd = 1)
abline(h = c(-2, 2), col = "red", lwd = 2, lty = 2)
abline(h = c(-3, 3), col = "darkred", lwd = 2, lty = 3)

# 图5：多重共线性诊断，VIF图
barplot(
  vif_values,
  main = "5. VIF",
  ylab = "VIF value",
  col = "blue",
  ylim = c(0, max(5, max(vif_values) + 1))
)
abline(h = 5, col = "red", lwd = 2, lty = 2)
abline(h = 10, col = "darkred", lwd = 2, lty = 3)

# 图6：高杠杆点诊断
plot(
  h,
  main = "6. Leverage",
  xlab = "Observation index",
  ylab = "Leverage value",
  pch = 19,
  col = "blue"
)
abline(h = leverage_threshold, col = "red", lwd = 2, lty = 2)

par(mfrow = c(1, 1))