predict.lm（）是如何计算置信区间和预测区间的？

我跑了一个回归：

CopierDataRegression <- lm(V1~V2, data=CopierData1) 

我的任务是获得一个

• 给定V2=6和V2=6的平均响应的90％ 置信区间
• 当V2=6时预测间隔为 90％。

我使用了下面的代码：

 X6 <- data.frame(V2=6) predict(CopierDataRegression, X6, se.fit=TRUE, interval="confidence", level=0.90) predict(CopierDataRegression, X6, se.fit=TRUE, interval="prediction", level=0.90) 

我得到了(87.3, 91.9)和(74.5, 104.8) ，这似乎是正确的，因为PI应该更宽。

两者的输出也包括se.fit = 1.39这是相同的。 我不明白这个标准错误是什么。 PI与CI的标准误差不应该大于？ 如何在R中find这两个不同的标准错误？

---

数据：

 CopierData1 <- structure(list(V1 = c(20L, 60L, 46L, 41L, 12L, 137L, 68L, 89L, 4L, 32L, 144L, 156L, 93L, 36L, 72L, 100L, 105L, 131L, 127L, 57L, 66L, 101L, 109L, 74L, 134L, 112L, 18L, 73L, 111L, 96L, 123L, 90L, 20L, 28L, 3L, 57L, 86L, 132L, 112L, 27L, 131L, 34L, 27L, 61L, 77L), V2 = c(2L, 4L, 3L, 2L, 1L, 10L, 5L, 5L, 1L, 2L, 9L, 10L, 6L, 3L, 4L, 8L, 7L, 8L, 10L, 4L, 5L, 7L, 7L, 5L, 9L, 7L, 2L, 5L, 7L, 6L, 8L, 5L, 2L, 2L, 1L, 4L, 5L, 9L, 7L, 1L, 9L, 2L, 2L, 4L, 5L)), .Names = c("V1", "V2"), class = "data.frame", row.names = c(NA, -45L)) 

Answer 1

简短的回答

 ## no need to specify "interval"; even no need to specify "level" z <- predict(CopierDataRegression, X6, se.fit=TRUE) 

用于CI的标准错误是：

 se.CI <- z$se.fit # [1] 1.396411 

用于PI的标准错误是：

 se.PI <- sqrt(z$se.fit^2 + z$residual.scale^2) # [1] 9.022228 

要计算90％显着性水平下的置信度/预测间隔，请执行：

 alpha <- qt((1-0.9)/2, df = z$df) # [1] -1.681071 CI <- z$fit + c(alpha, -alpha) * se.CI # [1] 87.28387 91.97880 PI <- z$fit + c(alpha, -alpha) * se.PI # [1] 74.46433 104.79833 

---

更详细的答案与math的细节

当你拟合一个线性模型时，拟合模型用matrixforms表示：

y = X％*％beta.hat

你可以得到beta.hat

 beta.hat <- CopierDataRegression$coefficients # (Intercept) V2 # -0.5801567 15.0352480 

和它的协方差matrix

 V <- vcov(CopierDataRegression) # (Intercept) V2 # (Intercept) 7.862086 -1.1927966 # V2 -1.192797 0.2333733 

当我们用预测matrixXp进行预测时，我们有预测的平均值：

 Xp <- model.matrix(~ V2, X6) pred <- as.numeric(Xp %*% beta.hat) # [1] 89.63133 

和预测方差：

 se2 <- unname(rowSums((Xp %*% V) * Xp)) # [1] 1.949963 

对于90％的置信区间，我们做

 alpha <- qt((1-0.9)/2, df = CopierDataRegression$df.residual) # [1] -1.681071 CI <- pred + c(alpha, -alpha) * sqrt(se2) # [1] 87.28387 91.97880 

预测区间是一个更宽的区间，因为它进一步说明了噪声sigma2不确定性。 Pearson对sigma2估计是

 sigma2 <- sum(CopierDataRegression$residuals ^ 2) / CopierDataRegression$df.residual # [1] 79.45063 

因此，90％水平的预测间隔是：

 PI <- pred + c(alpha, -alpha) * sqrt(se2 + sigma2) # [1] 74.46433 104.79833 

---

最后， z <- predict(CopierDataRegression, X6, se.fit=TRUE)返回

• z$fit ： pred以上;
• z$se.fit sqrt(se2)以上;
• z$df ：上面的df.residuals ;
• z$residual.scale ： sqrt(sigma2)以上。

Answer 2

我不知道是否有一个快速的方法来提取标准错误的预测间隔，但你总是可以解决SE的间隔（即使它不是超级优雅的方法）：

 m <- lm(V1 ~ V2, data = d) newdat <- data.frame(V2=6) tcrit <- qt(0.95, m$df.residual) a <- predict(m, newdat, interval="confidence", level=0.90) cat("CI SE", (a[1, "upr"] - a[1, "fit"]) / tcrit, "\n") b <- predict(m, newdat, interval="prediction", level=0.90) cat("PI SE", (b[1, "upr"] - b[1, "fit"]) / tcrit, "\n") 

请注意，CI SE与se.fit值相同。