R中手动计算多元线性回归预测区间与predict.lm()结果不符

嘿，我之前手动算预测区间的时候也踩过类似的坑！结合你的代码和无截距多元线性模型的场景，大概率是这几个细节没处理对，咱们一步步捋清楚：

先明确核心公式（无截距模型适用） #

预测区间是针对单个新观测的，公式是：
$$\hat{y} \pm t_{\alpha/2, n-p} \times \sqrt{\sigma^2 \left(1 + X_0 (X^T X)^{-1} X_0^T\right)}$$
其中：

• $\hat{y}$：预测值（$X_0 \hat{\beta}$）
• $t_{\alpha/2, n-p}$：t分布分位数，自由度是样本量$n$减去参数个数$p$（你的模型无截距，$p=6$）
• $\sigma^2$：模型的残差均方（MSE）
• $X_0$：新观测的自变量行向量
• $X$：模型的设计矩阵

最容易踩的几个坑 #

1. 混淆「置信区间」和「预测区间」 #

很多人手动计算时会漏掉公式里的+1——这是区分两者的关键：

• 置信区间是对总体均值的估计，公式里没有+1
• 预测区间是对单个新观测值的估计，需要加上新观测的随机误差项，所以要+1
  如果你的手动计算用了置信区间的公式，结果肯定和predict.lm(interval="prediction")差很多。

2. 无截距模型的自由度计算错误 #

你的模型没有截距，参数个数$p=6$（crim、nox、age、dis、tax、black），所以自由度是：
$$df = n - p = 506 - 6 = 500$$
如果错误地按有截距模型算（$df=506-7=499$），t分位数会有细微差异，最终区间也会偏差。

3. 残差均方（MSE）计算错误 #

MSE的正确计算是残差平方和除以自由度，而不是除以样本量：

# 正确的MSE
mse <- sum(residuals(model)^2) / (nrow(df) - length(coef(model)))
# 错误的MSE（很多人会犯）
wrong_mse <- mean(residuals(model)^2)

用错误的MSE会导致预测方差计算偏小，区间变窄。

4. 矩阵运算的转置/维度错误 #

$X_0$是行向量，计算$X_0 (X^T X)^{-1} X_0^T$时要注意转置，否则矩阵乘法会报错或者得到错误结果。R里用as.matrix()处理单行数据时默认是行矩阵，直接运算即可。

完整手动计算代码（和R结果对比） #

咱们用你的场景写完整代码，验证结果：

library(MASS)
data("Boston")
df <- Boston
model <- lm(medv ~ crim + nox + age + dis + tax + black + 0, data = df)
n_obs <- 3 # 第3个观测

# 1. 用R的predict得到基准结果
pred_r <- predict(model, newdata = df[n_obs,], interval = "prediction")
cat("R的预测结果：\n")
print(pred_r)

# 2. 手动计算
# 提取新观测的自变量矩阵（1×6）
X0 <- as.matrix(df[n_obs, c("crim", "nox", "age", "dis", "tax", "black")])
# 提取模型系数
beta_hat <- coef(model)
# 计算预测值
y_hat <- X0 %*% beta_hat
# 计算正确的MSE
mse <- sum(residuals(model)^2) / (nrow(df) - length(beta_hat))
# 计算(X^T X)的逆矩阵
XTX_inv <- solve(t(model.matrix(model)) %*% model.matrix(model))
# 计算预测方差
var_pred <- mse * (1 + X0 %*% XTX_inv %*% t(X0))
# 计算t分位数（95%置信水平，自由度500）
t_crit <- qt(0.975, df = nrow(df) - length(beta_hat))
# 计算上下限
lower <- y_hat - t_crit * sqrt(var_pred)
upper <- y_hat + t_crit * sqrt(var_pred)

# 3. 对比结果
cat("\n手动计算 vs R结果：\n")
print(data.frame(
  R_pred = pred_r[,"fit"],
  Manual_pred = as.numeric(y_hat),
  R_lower = pred_r[,"lwr"],
  Manual_lower = as.numeric(lower),
  R_upper = pred_r[,"upr"],
  Manual_upper = as.numeric(upper)
))

运行后你会发现，手动计算的结果和R的输出完全一致——如果之前有差异，对照上面的坑点检查就行。

内容的提问来源于stack exchange，提问作者Legacy