从predict.lm（）匹配R的置信区间

詹姆斯·万豪

我在匹配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）

谢谢！

用户名

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

具体来说，对于预测，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)

注意它们现在如何完美匹配。要进行数学分析，请看一下以下出色的文章：

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

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