带 lm() 的线性回归：聚合预测值的预测区间
Linear regression with `lm()`: prediction interval for aggregated predicted values
我正在使用 predict.lm(fit, newdata=newdata, interval="prediction") 来获取新观察的预测及其预测区间 (PI)。现在我想根据一个附加变量（即单个家庭预测的邮政编码级别的空间聚合）聚合（求和和平均）这些预测及其 PI。
我了解到 from StackExchange，您不能仅通过聚合预测区间的限制来聚合单个预测的预测区间。 post 非常有助于理解为什么不能这样做，但我很难将这一点转换为实际代码。答案是：
这是一个可重现的例子：
library(dplyr)
set.seed(123)
data(iris)
#Split dataset in training and prediction set
smp_size <- floor(0.75 * nrow(iris))
train_ind <- sample(seq_len(nrow(iris)), size = smp_size)
train <- iris[train_ind, ]
pred <- iris[-train_ind, ]
#Fit regression model
fit1 <- lm(Petal.Width ~ Petal.Length, data=train)
#Fit multiple linear regression model
fit2 <- lm(Petal.Width ~ Petal.Length + Sepal.Width + Sepal.Length, data=train)
#Predict Pedal.Width for new data incl prediction intervals for each prediction
predictions1<-predict(fit1, newdata=pred, interval="prediction")
predictions2<-predict(fit2, newdata=pred, interval="prediction")
# Aggregate data by summing predictions for species
#NOT correct for prediction intervals
predictions_agg1<-data.frame(predictions1,Species=pred$Species) %>%
group_by(Species) %>%
summarise_all(funs(sum,mean))
predictions_agg2<-data.frame(predictions2,Species=pred$Species) %>%
group_by(Species) %>%
summarise_all(funs(sum,mean))
我找不到描述使用 predict.lm() 时如何在 R 中正确聚合预测及其 PI 的好的教程或包。外面有东西吗？如果您能为我指出正确的方向，说明如何在 R 中执行此操作，我将不胜感激。
您的问题与我 2 年前回答的一个话题密切相关：. It provides an R implementation of Glen_b's answer on Cross Validated。感谢您引用该交叉验证线程；我不知道；也许我可以在那里发表评论，链接 Stack Overflow 线程。
我已经完善了我的原始答案，将逐行代码干净地包装成易于使用的函数 lm_predict 和 agg_pred。然后将解决您的问题简化为按组应用这些功能。
考虑问题中的 iris 示例，并使用第二个模型 fit2 进行演示。
set.seed(123)
data(iris)
#Split dataset in training and prediction set
smp_size <- floor(0.75 * nrow(iris))
train_ind <- sample(seq_len(nrow(iris)), size = smp_size)
train <- iris[train_ind, ]
pred <- iris[-train_ind, ]
#Fit multiple linear regression model
fit2 <- lm(Petal.Width ~ Petal.Length + Sepal.Width + Sepal.Length, data=train)
我们将 pred 按组 Species 拆分，然后在所有子数据帧上应用 lm_predict（使用 diag = FALSE）。
oo <- lapply(split(pred, pred$Species), lm_predict, lmObject = fit2, diag = FALSE)
要使用agg_pred，我们需要指定一个权重向量，其长度等于数据的数量。我们可以通过查询每个 oo[[i]]:
中 fit 的长度来确定这一点
n <- lengths(lapply(oo, "[[", 1))
#setosa versicolor virginica
# 11 13 14
如果聚合操作是sum，我们做
w <- lapply(n, rep.int, x = 1)
#List of 3
# $ setosa : num [1:11] 1 1 1 1 1 1 1 1 1 1 ...
# $ versicolor: num [1:13] 1 1 1 1 1 1 1 1 1 1 ...
# $ virginica : num [1:14] 1 1 1 1 1 1 1 1 1 1 ...
SUM <- Map(agg_pred, w, oo)
SUM[[1]] ## result for the first group, for example
#$mean
#[1] 2.499728
#
#$var
#[1] 0.1271554
#
#$CI
# lower upper
#1.792908 3.206549
#
#$PI
# lower upper
#0.999764 3.999693
sapply(SUM, "[[", "CI") ## some nice presentation for CI, for example
# setosa versicolor virginica
#lower 1.792908 16.41526 26.55839
#upper 3.206549 17.63953 28.10812
如果聚合操作是平均的，我们将 w 重新缩放 n 并调用 agg_pred.
w <- mapply("/", w, n)
#List of 3
# $ setosa : num [1:11] 0.0909 0.0909 0.0909 0.0909 0.0909 ...
# $ versicolor: num [1:13] 0.0769 0.0769 0.0769 0.0769 0.0769 ...
# $ virginica : num [1:14] 0.0714 0.0714 0.0714 0.0714 0.0714 ...
AVE <- Map(agg_pred, w, oo)
AVE[[2]] ## result for the second group, for example
#$mean
#[1] 1.3098
#
#$var
#[1] 0.0005643196
#
#$CI
# lower upper
#1.262712 1.356887
#
#$PI
# lower upper
#1.189562 1.430037
sapply(AVE, "[[", "PI") ## some nice presentation for CI, for example
# setosa versicolor virginica
#lower 0.09088764 1.189562 1.832255
#upper 0.36360845 1.430037 2.072496
This is great! Thank you so much! There is one thing I forgot to mention: in my actual application I need to sum ~300,000 predictions which would create a full variance-covariance matrix which is about ~700GB in size. Do you have any idea if there is a computationally more efficient way to directly get to the sum of the variance-covariance matrix?
使用原版问答改版中提供的fast_agg_pred功能，从头再来
set.seed(123)
data(iris)
#Split dataset in training and prediction set
smp_size <- floor(0.75 * nrow(iris))
train_ind <- sample(seq_len(nrow(iris)), size = smp_size)
train <- iris[train_ind, ]
pred <- iris[-train_ind, ]
#Fit multiple linear regression model
fit2 <- lm(Petal.Width ~ Petal.Length + Sepal.Width + Sepal.Length, data=train)
## list of new data
newdatlist <- split(pred, pred$Species)
n <- sapply(newdatlist, nrow)
#setosa versicolor virginica
# 11 13 14
如果聚合操作是sum，我们做
w <- lapply(n, rep.int, x = 1)
SUM <- mapply(fast_agg_pred, w, newdatlist,
MoreArgs = list(lmObject = fit2, alpha = 0.95),
SIMPLIFY = FALSE)
如果聚合操作是平均的，我们做
w <- mapply("/", w, n)
AVE <- mapply(fast_agg_pred, w, newdatlist,
MoreArgs = list(lmObject = fit2, alpha = 0.95),
SIMPLIFY = FALSE)
请注意，在这种情况下我们不能使用 Map，因为我们需要为 fast_agg_pred 提供更多参数。在这种情况下使用 mapply，MoreArgs 和 SIMPLIFY。
我正在使用 predict.lm(fit, newdata=newdata, interval="prediction") 来获取新观察的预测及其预测区间 (PI)。现在我想根据一个附加变量（即单个家庭预测的邮政编码级别的空间聚合）聚合（求和和平均）这些预测及其 PI。
我了解到 from StackExchange，您不能仅通过聚合预测区间的限制来聚合单个预测的预测区间。 post 非常有助于理解为什么不能这样做，但我很难将这一点转换为实际代码。答案是：
这是一个可重现的例子：
library(dplyr)
set.seed(123)
data(iris)
#Split dataset in training and prediction set
smp_size <- floor(0.75 * nrow(iris))
train_ind <- sample(seq_len(nrow(iris)), size = smp_size)
train <- iris[train_ind, ]
pred <- iris[-train_ind, ]
#Fit regression model
fit1 <- lm(Petal.Width ~ Petal.Length, data=train)
#Fit multiple linear regression model
fit2 <- lm(Petal.Width ~ Petal.Length + Sepal.Width + Sepal.Length, data=train)
#Predict Pedal.Width for new data incl prediction intervals for each prediction
predictions1<-predict(fit1, newdata=pred, interval="prediction")
predictions2<-predict(fit2, newdata=pred, interval="prediction")
# Aggregate data by summing predictions for species
#NOT correct for prediction intervals
predictions_agg1<-data.frame(predictions1,Species=pred$Species) %>%
group_by(Species) %>%
summarise_all(funs(sum,mean))
predictions_agg2<-data.frame(predictions2,Species=pred$Species) %>%
group_by(Species) %>%
summarise_all(funs(sum,mean))
我找不到描述使用 predict.lm() 时如何在 R 中正确聚合预测及其 PI 的好的教程或包。外面有东西吗？如果您能为我指出正确的方向，说明如何在 R 中执行此操作，我将不胜感激。
您的问题与我 2 年前回答的一个话题密切相关：
我已经完善了我的原始答案，将逐行代码干净地包装成易于使用的函数 lm_predict 和 agg_pred。然后将解决您的问题简化为按组应用这些功能。
考虑问题中的 iris 示例，并使用第二个模型 fit2 进行演示。
set.seed(123)
data(iris)
#Split dataset in training and prediction set
smp_size <- floor(0.75 * nrow(iris))
train_ind <- sample(seq_len(nrow(iris)), size = smp_size)
train <- iris[train_ind, ]
pred <- iris[-train_ind, ]
#Fit multiple linear regression model
fit2 <- lm(Petal.Width ~ Petal.Length + Sepal.Width + Sepal.Length, data=train)
我们将 pred 按组 Species 拆分，然后在所有子数据帧上应用 lm_predict（使用 diag = FALSE）。
oo <- lapply(split(pred, pred$Species), lm_predict, lmObject = fit2, diag = FALSE)
要使用agg_pred，我们需要指定一个权重向量，其长度等于数据的数量。我们可以通过查询每个 oo[[i]]:
fit 的长度来确定这一点
n <- lengths(lapply(oo, "[[", 1))
#setosa versicolor virginica
# 11 13 14
如果聚合操作是sum，我们做
w <- lapply(n, rep.int, x = 1)
#List of 3
# $ setosa : num [1:11] 1 1 1 1 1 1 1 1 1 1 ...
# $ versicolor: num [1:13] 1 1 1 1 1 1 1 1 1 1 ...
# $ virginica : num [1:14] 1 1 1 1 1 1 1 1 1 1 ...
SUM <- Map(agg_pred, w, oo)
SUM[[1]] ## result for the first group, for example
#$mean
#[1] 2.499728
#
#$var
#[1] 0.1271554
#
#$CI
# lower upper
#1.792908 3.206549
#
#$PI
# lower upper
#0.999764 3.999693
sapply(SUM, "[[", "CI") ## some nice presentation for CI, for example
# setosa versicolor virginica
#lower 1.792908 16.41526 26.55839
#upper 3.206549 17.63953 28.10812
如果聚合操作是平均的，我们将 w 重新缩放 n 并调用 agg_pred.
w <- mapply("/", w, n)
#List of 3
# $ setosa : num [1:11] 0.0909 0.0909 0.0909 0.0909 0.0909 ...
# $ versicolor: num [1:13] 0.0769 0.0769 0.0769 0.0769 0.0769 ...
# $ virginica : num [1:14] 0.0714 0.0714 0.0714 0.0714 0.0714 ...
AVE <- Map(agg_pred, w, oo)
AVE[[2]] ## result for the second group, for example
#$mean
#[1] 1.3098
#
#$var
#[1] 0.0005643196
#
#$CI
# lower upper
#1.262712 1.356887
#
#$PI
# lower upper
#1.189562 1.430037
sapply(AVE, "[[", "PI") ## some nice presentation for CI, for example
# setosa versicolor virginica
#lower 0.09088764 1.189562 1.832255
#upper 0.36360845 1.430037 2.072496
This is great! Thank you so much! There is one thing I forgot to mention: in my actual application I need to sum ~300,000 predictions which would create a full variance-covariance matrix which is about ~700GB in size. Do you have any idea if there is a computationally more efficient way to directly get to the sum of the variance-covariance matrix?
使用原版问答改版中提供的fast_agg_pred功能，从头再来
set.seed(123)
data(iris)
#Split dataset in training and prediction set
smp_size <- floor(0.75 * nrow(iris))
train_ind <- sample(seq_len(nrow(iris)), size = smp_size)
train <- iris[train_ind, ]
pred <- iris[-train_ind, ]
#Fit multiple linear regression model
fit2 <- lm(Petal.Width ~ Petal.Length + Sepal.Width + Sepal.Length, data=train)
## list of new data
newdatlist <- split(pred, pred$Species)
n <- sapply(newdatlist, nrow)
#setosa versicolor virginica
# 11 13 14
如果聚合操作是sum，我们做
w <- lapply(n, rep.int, x = 1)
SUM <- mapply(fast_agg_pred, w, newdatlist,
MoreArgs = list(lmObject = fit2, alpha = 0.95),
SIMPLIFY = FALSE)
如果聚合操作是平均的，我们做
w <- mapply("/", w, n)
AVE <- mapply(fast_agg_pred, w, newdatlist,
MoreArgs = list(lmObject = fit2, alpha = 0.95),
SIMPLIFY = FALSE)
请注意，在这种情况下我们不能使用 Map，因为我们需要为 fast_agg_pred 提供更多参数。在这种情况下使用 mapply，MoreArgs 和 SIMPLIFY。