predict.lm() 如何计算置信区间和预测区间？
How does predict.lm() compute confidence interval and prediction interval?
我运行一个回归：
CopierDataRegression <- lm(V1~V2, data=CopierData1)
我的任务是获得一个
给定 V2=6 和 的平均响应的 - 90% 置信区间
- 90% 预测区间 when
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))
我不知道是否有快速提取预测区间标准误差的方法，但您始终可以反解 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 的值相同。
当指定 interval 和 level 参数时，predict.lm 可以 return 置信区间 (CI) 或预测区间 (PI)。此答案显示了如何在不设置这些参数的情况下获取 CI 和 PI。有两种方式：
- 使用
predict.lm 的中间阶段结果；
- 一切从零开始。
了解如何使用这两种方法可以使您对预测过程有一个透彻的了解。
请注意，我们将只涵盖 predict.lm 的 type = "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
What is se.fit?
z$se.fit 是预测平均值 z$fit 的标准误差，用于为 z$fit 构造 CI。我们还需要具有自由度 z$df.
的 t 分布的分位数
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") 一致。
What is the standard error for PI?
PI 比 CI 宽，因为它考虑了残差：
variance_of_PI = variance_of_CI + variance_of_residual
请注意，这是逐点定义的。对于非加权线性回归（如您的示例），残差在任何地方都相等（称为 homoscedasticity），它是 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 更容易（即，在 predict.lm 中使用 type = "prediction" 时不设置 newdata），因为权重是已知的（您必须通过weight 使用 lm 时的参数）。对于样本外预测（即，您将 newdata 传递给 predict.lm），predict.lm 希望您告诉它应如何对残差方差进行加权。您需要在 predict.lm 中使用参数 pred.var 或 weights，否则您会收到来自 predict.lm 的警告，抱怨构建 PI 的信息不足。以下内容转自?predict.lm：
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 的构造不受回归类型的影响。
一切从零开始
基本上我们想知道如何在z.
中获取fit、se.fit、df和residual.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)，其中 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
"Do everything from scratch" 中的代码在本问答中被干净地组织成一个函数 lm_predict：.
我运行一个回归：
CopierDataRegression <- lm(V1~V2, data=CopierData1)
我的任务是获得一个
-
给定
- 90% 置信区间
- 90% 预测区间 when
V2=6.
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))
我不知道是否有快速提取预测区间标准误差的方法，但您始终可以反解 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 的值相同。
当指定 interval 和 level 参数时，predict.lm 可以 return 置信区间 (CI) 或预测区间 (PI)。此答案显示了如何在不设置这些参数的情况下获取 CI 和 PI。有两种方式：
- 使用
predict.lm的中间阶段结果； - 一切从零开始。
了解如何使用这两种方法可以使您对预测过程有一个透彻的了解。
请注意，我们将只涵盖 predict.lm 的 type = "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
What is
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") 一致。
What is the standard error for PI?
PI 比 CI 宽，因为它考虑了残差：
variance_of_PI = variance_of_CI + variance_of_residual
请注意，这是逐点定义的。对于非加权线性回归（如您的示例），残差在任何地方都相等（称为 homoscedasticity），它是 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 更容易（即，在 predict.lm 中使用 type = "prediction" 时不设置 newdata），因为权重是已知的（您必须通过weight 使用 lm 时的参数）。对于样本外预测（即，您将 newdata 传递给 predict.lm），predict.lm 希望您告诉它应如何对残差方差进行加权。您需要在 predict.lm 中使用参数 pred.var 或 weights，否则您会收到来自 predict.lm 的警告，抱怨构建 PI 的信息不足。以下内容转自?predict.lm：
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 的构造不受回归类型的影响。
一切从零开始
基本上我们想知道如何在z.
fit、se.fit、df和residual.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)，其中 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
"Do everything from scratch" 中的代码在本问答中被干净地组织成一个函数 lm_predict：