带有“lm()”的线性回归：聚合预测值的预测区间

Question

我正在使用 predict.lm(fit, newdata=newdata, interval="prediction") 来获取预测及其预测区间 (PI) 以进行新观察。现在，我想根据一个附加变量（即单个家庭的邮政编码级别的空间聚合）来汇总（求和和平均）这些预测及其 PI。

我了解到from StackExchange，您不能仅通过聚合预测区间的限制来聚合单个预测的预测区间。这篇文章非常有助于理解为什么不能这样做，但我很难将这一点翻译成实际的代码。答案是：

这是一个可重现的例子：

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 中执行此操作的正确方向，将不胜感激。

Answer 1

您的问题与我 2 年前回答的一个主题密切相关：linear model with `lm`: how to get prediction variance of sum of predicted values。它提供了Glen_b's answer on Cross Validated 的 R 实现。感谢您引用该交叉验证线程；我不知道；也许我可以在此处留下评论，链接 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)

我们按组Species 拆分pred，然后将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 

如果聚合操作是求和，我们做

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

---

这太棒了！太感谢了！我忘了提一件事：在我的实际应用中，我需要对约 300,000 个预测求和，这将创建一个大小约为 700GB 的完整方差-协方差矩阵。您是否知道是否有一种计算上更有效的方法可以直接得到方差-协方差矩阵的总和？

使用原问答修订版中提供的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 

如果聚合操作是求和，我们做

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 一起使用。