Линейная регрессия с помощью `lm()`: интервал прогнозирования для агрегированных прогнозных значений
Я использую predict.lm(fit, newdata=newdata, interval="prediction") получить прогнозы и их интервалы прогнозирования (PI) для новых наблюдений. Теперь я хотел бы агрегировать (суммировать и означать) эти прогнозы и их ИП на основе дополнительной переменной (т. Е. Пространственной агрегации на уровне прогнозов почтового индекса для отдельных домохозяйств).
Я узнал из 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))
Я не смог найти хороший учебник или пакет, который описывает, как правильно агрегировать прогнозы и их PI в R при использовании predict.lm(), Там что-то есть? Буду очень признателен, если вы укажете мне правильное направление, как это сделать в R.
1 ответ
Ваш вопрос тесно связан с темой, на которую я ответил 2 года назад: линейная модель с `lm`: как получить прогнозную дисперсию суммы прогнозируемых значений. Он обеспечивает R-реализацию ответа Glen_b на Cross Validated. Спасибо за цитирование этой темы 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 нам нужно указать весовой вектор, длина которого равна количеству данных. Мы можем определить это, посоветовавшись с длиной fit в каждом oo[[i]]:
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
Это замечательно! Спасибо вам большое! Я забыл упомянуть одну вещь: в моем реальном приложении мне нужно суммировать ~300000 прогнозов, которые создали бы полную дисперсионно-ковариационную матрицу размером около ~700 ГБ. Есть ли у вас какие-либо идеи, если существует более эффективный в вычислительном отношении способ прямого получения суммы дисперсионно-ковариационной матрицы?
Использовать fast_agg_pred функция, представленная в ревизии оригинальной Q & A. Давайте начнем все сначала.
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,