Как я могу применить сгруппированные данные к сгруппированным моделям, используя метлу и dplyr?

Question

Я хотел бы сделать эквивалент подгонки модели gpm (галлонов на милю = 1/mpg) к wt в наборе данных mtcars. Это кажется легким:

data(mtcars)
library(dplyr)
library(tidyr)
library(broom)
library(ggplot2)
library(scales)

mtcars2 <-
    mtcars %>%
    mutate(gpm = 1 / mpg) %>%
    group_by(cyl, am)

lm1 <-
    mtcars2 %>%
    do(fit = lm(gpm ~ wt, data = .))

Это дает мне строку данных с шестью строками, как и ожидалось.

Этот график подтверждает, что существует шесть групп:

p1 <-
    qplot(wt, gpm, data = mtcars2) +
    facet_grid(cyl ~ am) +
    stat_smooth(method='lm',se=FALSE, fullrange = TRUE) +
    scale_x_continuous(limits = c(0,NA)) 

Я могу использовать augment () для получения подходящих результатов:

lm1 %>% augment(fit)

Это дает мне 32 строки, по одной на каждую строку в mtcars2, как и ожидалось.

Теперь задача: я хотел бы получить подходящие результаты, используя newdata, где я увеличил вес на cyl/4:

newdata <-
    mtcars2 %>%
    mutate(
        wt = wt + cyl/4)

Я ожидаю, что это создаст фрейм данных того же размера, что и lm1 %>% augment(fit): по одной строке для каждой строки в newdata, потому что метла будет сопоставлять модели и новые данные по групповым переменным cyl и am.

К несчастью,

pred1 <-
    lm1 %>%
    augment(
        fit,
        newdata = newdata)

дает мне фрейм данных с 192 строками (= 6 x 32), по-видимому, подгоняя каждую модель к каждому ряду новых данных.

Из других источников я понял, что фреймы данных group_by и rowwise не совместимы, поэтому lm1 не сгруппирован, и augment не может связать модели и новые данные. Есть ли другой шаблон дизайна, который позволяет мне сделать это? Было бы хорошо, если бы он был таким же простым и прозрачным, как и приведенная выше попытка, но более важно, чтобы он работал.

Вот мой sessionInfo():

> sessionInfo()
R version 3.3.1 (2016-06-21)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United States.1252 
[2] LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252
[4] LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] scales_0.4.0  ggplot2_2.1.0 broom_0.4.1   tidyr_0.6.0   dplyr_0.5.0  

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.7      magrittr_1.5     mnormt_1.5-4     munsell_0.4.3   
 [5] colorspace_1.2-6 lattice_0.20-34  R6_2.1.3         stringr_1.1.0   
 [9] plyr_1.8.4       tools_3.3.1      parallel_3.3.1   grid_3.3.1      
[13] nlme_3.1-128     gtable_0.2.0     psych_1.6.9      DBI_0.5-1       
[17] lazyeval_0.2.0   assertthat_0.1   tibble_1.2       reshape2_1.4.1  
[21] labeling_0.3     stringi_1.1.1    compiler_3.3.1   foreign_0.8-67  

РЕДАКТИРОВАТЬ:

@aosmith: я изучаю ваш второй вариант, и мне это нравится. Однако, когда я пробую это на моих реальных данных, у меня возникает проблема с командой mutate: она возвращает "Ошибка: augment не знает, как обращаться с данными списка классов".

Мой реальный код больше похож на:

newdata %>% 
dplyr::select(cyl, am, wt) %>% # wt holds new predictor values
group_by(cyl, am) %>%
nest() %>%
inner_join(regressions, .) %>% 
## looks like yours at this point
mutate(pred = list(augment(fit, newdata = data))) %>% # Error here
unnest(pred)

Когда я говорю, что это выглядит как ваш, я имею в виду, что у меня есть следующие столбцы (переименовано здесь для согласованности): ID (chr), attr1 (dbl), cyl (dbl), am (chr), fit (список) и данные (список). У вас есть цил, утра (DBL), подходят и данные. Я изменил свой я на DBL, но это не помогло.

Я думаю, что разница в том, что у меня есть 3 (ID... аналогично именам строк в mtcars) x 2 (цил) x 2 (am) единиц в этом образце (каждый образец имеет 12 измерений), в то время как пример mtcars имеет 3 (цил) х 2 (ам) ячеек х случайное количество типов автомобилей на ячейку. В моем анализе мне нужно увидеть значения идентификаторов, но новые данные в равной степени применяются ко всем единицам. Если это помогает, думайте об этом как о скорости встречного ветра, примененной к каждой машине в тесте. Означает ли это, что причина жалобы дополнения не связана с данными списка классов?

РЕДАКТИРОВАТЬ: Объединение идентификатора с новыми данными (используя full=TRUE) решило последнюю проблему. В настоящее время я использую ваше первое предложенное решение.

Answer 1

Я использовал map2 из пакета мурлыкать для такого рода ситуации. map2 проходит через элементы двух списков одновременно. Списки должны быть одинаковой длины и в одинаковом порядке.

Элементы списков используются в качестве аргументов для некоторой функции, которую вы хотите применить (augment, в твоем случае). Здесь ваши два списка будут список моделей и список наборов данных (один список для каждого cyl/am комбинация).

С помощью map2_df возвращает результаты в виде data.frame вместо списка.

library(purrr)

Я сделал список data.frames для прогнозирования с помощью split, Порядок факторов, на которые делятся, определил порядок списка, поэтому я убедился, что он был в том же порядке, что и список. lm1,

test_split = split(newdata, list(newdata$am, newdata$cyl)

map2_df(lm1$fit, test_split, ~augment(.x, newdata = .y))

Чтобы не беспокоиться о заказе так много, вы могли бы nest данные прогноза по группам, присоедините это к lm1и вернуть результаты augment как список для unnesting.

newdata %>%
    group_by(cyl, am) %>%
    nest() %>%
    inner_join(lm1, .) %>%
    mutate(pred = list(augment(fit, newdata = data))) %>%
    unnest(pred)