首页 文章

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

提问于
浏览
8

我跑回了一个回归:

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

    指定 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 .

  • 27

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

相关问题