predict.lm()如何计算置信区间和预测区间?

我跑回了一个回归:

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

我的任务是获得一个

  • 90% confidence intervalV2=6V2=6 的平均响应
    V2=6
  • 90% prediction interval .

我使用了以下代码:

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 ,这是相同的 . I don't understand what this standard error is. Shouldn't the standard error be larger for the PI vs. the CI? How do I find these two different standard errors in R?
enter image description here


数据:

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))

回答(2)

2 years ago

指定 intervallevel 参数时, predict.lm 可以返回置信区间(CI)或预测区间(PI) . 此答案显示如何在不设置这些参数的情况下获取CI和PI . 有两种方法:

  • 使用 predict.lm 的中期结果;

  • 从头开始做一切 .

了解如何使用这两种方式可以让您全面了解预测过程 .

请注意,我们只会涵盖 predict.lmtype = "response" (默认)大小写 . 讨论 type = "terms" 超出了这个答案的范围 .


设置

我在这里收集你的代码,以帮助其他读者复制,粘贴和运行 . 我还更改变量名称,以便它们具有更清晰的含义 . 另外,我扩展 newdat 以包含多行,以显示我们的计算是"vectorized" .

dat <- 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))

lmObject <- lm(V1 ~ V2, data = dat)

newdat <- data.frame(V2 = c(6, 7))

以下是 predict.lm 的输出,稍后将与我们的手动计算进行比较 .

predict(lmObject, newdat, se.fit = TRUE, interval = "confidence", level = 0.90)
#$fit
#        fit       lwr      upr
#1  89.63133  87.28387  91.9788
#2 104.66658 101.95686 107.3763
#
#$se.fit
#       1        2 
#1.396411 1.611900 
#
#$df
#[1] 43
#
#$residual.scale
#[1] 8.913508

predict(lmObject, newdat, se.fit = TRUE, interval = "prediction", level = 0.90)
#$fit
#        fit      lwr      upr
#1  89.63133 74.46433 104.7983
#2 104.66658 89.43930 119.8939
#
#$se.fit
#       1        2 
#1.396411 1.611900 
#
#$df
#[1] 43
#
#$residual.scale
#[1] 8.913508

使用来自predict.lm的中间阶段结果

## use `se.fit = TRUE`
z <- predict(lmObject, newdat, se.fit = TRUE)
#$fit
#        1         2 
# 89.63133 104.66658 
#
#$se.fit
#       1        2 
#1.396411 1.611900 
#
#$df
#[1] 43
#
#$residual.scale
#[1] 8.913508

se.fit是什么?

z$se.fit 是预测平均值 z$fit 的标准误差,用于构造 z$fit 的CI . 我们还需要具有自由度的分布的分位数 z$df .

alpha <- 0.90  ## 90%
Qt <- c(-1, 1) * qt((1 - alpha) / 2, z$df, lower.tail = FALSE)
#[1] -1.681071  1.681071

## 90% confidence interval
CI <- z$fit + outer(z$se.fit, Qt)
colnames(CI) <- c("lwr", "upr")
CI
#        lwr      upr
#1  87.28387  91.9788
#2 101.95686 107.3763

我们看到这与 predict.lm(, interval = "confidence") 一致 .

PI的标准错误是什么?

PI比CI更宽,因为它考虑了剩余方差:

variance_of_PI = variance_of_CI + variance_of_residual

请注意,这是逐点定义的 . 对于非加权线性回归(如在您的示例中),残差方差在任何地方都相等(称为同方差性),并且它是 z$residual.scale ^ 2 . 因此,PI的标准误差是

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

PI构建为

PI <- z$fit + outer(se.PI, Qt)
colnames(PI) <- c("lwr", "upr")
PI
#       lwr      upr
#1 74.46433 104.7983
#2 89.43930 119.8939

我们看到这与 predict.lm(, interval = "prediction") 一致 .

remark

如果你有一个权重线性回归,那么事情会更复杂,其中残差方差在任何地方都不相等,因此应该加权.756734_ . 为拟合值构造PI更容易(也就是说,在 predict.lm 中使用 type = "prediction" 时不设置 newdata ),因为权重是已知的(使用 lm 时必须通过 weight 参数提供) . 对于样本外预测(即,您将 newdata 传递给 predict.lm ), predict.lm 希望您告诉它应如何对残差方差进行加权 . 您需要在 predict.lm 中使用参数 pred.varweights ,否则您会收到来自 predict.lm 的警告,抱怨构建PI的信息不足 . 以下引用自 ?predict.lm

预测间隔适用于每种情况下的单次观察
在'newdata'(或默认情况下,用于拟合的数据)中有错误
方差'pred.var' . 这可以是'res.var'的倍数
sigma的估计值^ 2:默认是假设未来
观察结果具有与用于的观察结果相同的误差方差
配件 . 如果提供'权重',则将其反向用作
比例因子 . 对于加权拟合,如果预测是针对
原始数据框,'权重'默认为用于的权重
模型适合,带有警告,因为它可能不是预期的
结果 . 如果拟合是加权的并且给出了'newdata',那么
默认是假设持续预测方差,并带有警告 .

请注意,CI的构造不受回归类型的影响 .


从头开始做一切

基本上我们想知道如何在 z 中获得 fitse.fitdfresidual.scale .

预测平均值可以通过矩阵向量乘法 Xp %*% b 来计算,其中 Xp 是线性预测矩阵, b 是回归系数向量 .

Xp <- model.matrix(delete.response(terms(lmObject)), newdat)
b <- coef(lmObject)
yh <- c(Xp %*% b)  ## c() reshape the single-column matrix to a vector
#[1]  89.63133 104.66658

我们看到这与 z$fit 一致 . yh 的方差 - 协方差是 Xp %*% V %*% t(Xp) ,其中 Vb 的方差 - 协方差矩阵,可以通过

V <- vcov(lmObject)  ## use `vcov` function in R
#             (Intercept)         V2
# (Intercept)    7.862086 -1.1927966
# V2            -1.192797  0.2333733

不需要 yh 的完全方差 - 协方差矩阵来计算逐点CI或PI . 我们只需要它的主对角线 . 因此,我们可以通过提高效率来实现 diag(Xp %*% V %*% t(Xp))

var.fit <- rowSums((Xp %*% V) * Xp)  ## point-wise variance for predicted mean
#       1        2 
#1.949963 2.598222 

sqrt(var.fit)  ## this agrees with `z$se.fit`
#       1        2 
#1.396411 1.611900

在拟合模型中可以轻松获得剩余自由度:

dof <- df.residual(lmObject)
#[1] 43

最后,要计算残差方差,请使用Pearson估算器:

sig2 <- c(crossprod(lmObject$residuals)) / dof
# [1] 79.45063

sqrt(sig2)  ## this agrees with `z$residual.scale`
#[1] 8.913508

remark

请注意,在加权回归的情况下, sig2 应计算为

sig2 <- c(crossprod(sqrt(lmObject$weights) * lmObject$residuals)) / dof

附录:一个模仿predict.lm的自编函数

"Do everything from scratch"中的代码已完整地组织到此问答中的函数 lm_predictlinear model with lm: how to get prediction variance of sum of predicted values .

2 years ago

我不知道是否有一种快速的方法来提取预测间隔的标准误差,但是你可以随时回溯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 中的值相同 .