predict.lm()如何计算置信区间和预测区间？
Mit*_*tty 11 regression r prediction linear-regression lm
我跑回了一个回归:
CopierDataRegression <- lm(V1~V2, data=CopierData1)
我的任务是获得一个
- 给出和 的平均响应的90%置信区间
V2=6 - 当预测间隔为 90%时
V2=6.
我使用了以下代码:
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中找到这两个不同的标准错误？
数据:
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))
李哲源*_*李哲源 36
指定interval和level参数时,predict.lm可以返回置信区间(CI)或预测区间(PI).此答案显示如何在不设置这些参数的情况下获取CI和PI.有两种方法:
- 使用中期结果
predict.lm; - 从头开始做一切.
了解如何使用这两种方式可以让您彻底了解预测过程.
请注意,我们仅涵盖type = "response"(默认)情况predict.lm.讨论type = "terms"超出了这个答案的范围.
建立
我在这里收集你的代码,以帮助其他读者复制,粘贴和运行.我还更改变量名称,以便它们具有更清晰的含义.另外,我扩展newdat为包含多行,以显示我们的计算是"矢量化"的.
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,用于构造CI z$fit.我们还需要具有一定自由度的分布的分位数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").
备注
如果你有权重线性回归,那么事情会更复杂,其中残差方差在任何地方都不相等,所以z$residual.scale ^ 2应该加权.这是比较容易构建拟合值PI(即,您没有设置newdata使用时type = "prediction"的predict.lm)(您必须通过提供它,因为权重被称为weight使用时的说法lm).对于样本外预测(即,您传递newdata给predict.lm),predict.lm期望您告诉它应如何加权残差方差.您需要使用参数pred.var或weightsin predict.lm,否则您会因predict.lm抱怨构建PI的信息不足而收到警告.以下是引用自?predict.lm:
Run Code Online (Sandbox Code Playgroud)The prediction intervals are for a single observation at each case in ‘newdata’ (or by default, the data used for the fit) with error variance(s) ‘pred.var’. This can be a multiple of ‘res.var’, the estimated value of sigma^2: the default is to assume that future observations have the same error variance as those used for fitting. If ‘weights’ is supplied, the inverse of this is used as a scale factor. For a weighted fit, if the prediction is for the original data frame, ‘weights’ defaults to the weights used for the model fit, with a warning since it might not be the intended result. If the fit was weighted and ‘newdata’ is given, the default is to assume constant prediction variance, with a warning.
请注意,CI的构造不受回归类型的影响.
从头开始做一切
基本上,我们想知道如何获取fit,se.fit,df和residual.scale在z.
预测平均值可以通过矩阵向量乘法来计算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),其中V在所述方差-协方差矩阵的b其可通过计算
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
备注
注意,在加权回归的情况下,sig2应计算为
sig2 <- c(crossprod(sqrt(lmObject$weights) * lmObject$residuals)) / dof
附录:模仿的自编函数 predict.lm
"从头开始做所有事情"中的代码已经完整地组织成lm_predict这个Q&A中的函数:线性模型lm:如何获得预测值之和的预测方差.
我不知道是否有一种快速方法来提取预测区间的标准误差，但您始终可以反求解 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。