r - 从 predict.lm() 匹配 R 的置信区间

Question

我在匹配 R 的线性模型置信度和预测区间以进行点估计时遇到问题：

set.seed(3)
x <- rnorm(392, 20, 5)
y <- 2*x + 3 + rnorm(392, sd=3)

lm.fit <- lm(y~x)

(est <- as.vector(matrix(c(1,20),nrow=1,ncol=2) %*% lm.fit$coefficients))
(ci <- est + c(-1,1) * qt(.975, df=lm.fit$df.residual)*sqrt(var(lm.fit$residuals)/(lm.fit$df.residual)))
(pred <- predict(lm.fit,data.frame(x=(c(20))), interval="confidence", se.fit = TRUE))

sqrt(var(lm.fit$residuals)/(lm.fit$df.residual)) 
pred$se.fit

如果您运行上面的代码，您会看到点估计值是相同的，但我得到的标准误差略有不同，因此置信区间略有偏差。R如何匹配标准误差？

我在 predict.lm() 的文档（如下）中看到 R 在 predict.lm 函数中使用 pred.var，但 pred.var 使用 res.var，并且没有说明 res.var 的来源。res.var 也不匹配残差的方差：

var(lm.fit$residuals)

（http://127.0.0.1:25345/library/stats/html/predict.lm.html）

谢谢！

标签： rlinear-regression

Answer 1

您混淆了预测区间标准误差（对于新观察）与置信区间标准误差（对于平均响应）。

具体来说，对于预测，SE 在标准误差的表达式中多了一个 1。我在下面展示了完整的计算。

请看下面的代码：

set.seed(3)
x <- rnorm(392, 20, 5)
y <- 2*x + 3 + rnorm(392, sd=3)

lm.fit <- lm(y~x)

如果这是一个新的观察，你需要一个 1 的比例因子（见第一项

pred <- predict(lm.fit,data.frame(x=(c(20))), interval="prediction", se.fit = TRUE)
MSE <- sum(lm.fit$residuals^2)/(length(x)-2)
Scaling_factor = (1 + 1/length(x) + ((20 - mean(x))^2) / sum( (x- mean(x))^2     ))
est - qt(.975, length(x)-2)*sqrt(MSE*Scaling_factor)
est + qt(.975, length(x)-2)*sqrt(MSE*Scaling_factor)

如果您要预测平均响应（这是线性回归通常所做的），则不需要该术语。

pred <- predict(lm.fit,data.frame(x=(c(20))), interval="confidence", se.fit = TRUE)
MSE <- sum(lm.fit$residuals^2)/(length(x)-2)
Scaling_factor = (1/length(x) + ((20 - mean(x))^2) / sum( (x- mean(x))^2     ))
est - qt(.975, length(x)-2)*sqrt(MSE*Scaling_factor)
est + qt(.975, length(x)-2)*sqrt(MSE*Scaling_factor)

注意它们现在是如何完美匹配的。对于数学分析，看看这篇很棒的文章：

平均响应和单个响应的线性回归

具体来说，对于平均响应：

但是，对于个人观察：

那个额外的“1”就是让你与众不同的地方。基本思想是，当您预测单个响应时，存在额外的随机性，因为您对 MSE 的估计计算的是平均响应值周围的变异性，而不是单个响应值。

直觉上，这很有意义。由于您的观察次数 -> 无穷大，1 使标准误差不会变为 0（因为任何人总是存在一些可变性）。

希望有帮助！