Как применить сгруппированные данные к сгруппированным моделям с помощью метлы и 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)) 

можно использовать увеличивающие(), чтобы получить встроенные выходы:

lm1 %>% augment(fit) 

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

Теперь задача: Я хотел бы получить встроенные выходы с помощью NewData, где я увеличивается вес на цил/4:

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

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

К сожалению,

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

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

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

Вот мой 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 

EDIT:

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

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

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 (ДВМ), cyl (dbl), am (chr), fit (список) и данные (список). У вас есть цил, am (dbl), подгонка и данные. Я изменил свой ам на dbl, но это не помогло.

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

EDIT: Объединение идентификатора с помощью newdata (с использованием full = TRUE) решило последнюю проблему.В настоящее время я использую ваше первое предложенное решение.

Answer 1

Я использовал map2 из упаковки purrr для такого рода ситуации. 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)