如何使用 broom 和 dplyr 将分组数据应用于分组模型?

Bil*_*llH 5 r dplyr broom

我想做相当于将 gpm(每英里加仑数 = 1/mpg)模型拟合到 mtcars 数据集中的 wt。这似乎很容易:

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 = .))
Run Code Online (Sandbox Code Playgroud)

正如预期的那样,这为我提供了一个 6 行的 rowwise 数据框。

此图确认有六组:

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)) 
Run Code Online (Sandbox Code Playgroud)

我可以使用 Augment() 来获得拟合的输出:

lm1 %>% augment(fit)
Run Code Online (Sandbox Code Playgroud)

正如预期的那样,这给了我 32 行,mtcars2 中的每一行。

现在的挑战是:我想使用 newdata 获得拟合输出,其中我已将 wt 增加了 cyl/4:

newdata <-
    mtcars2 %>%
    mutate(
        wt = wt + cyl/4)
Run Code Online (Sandbox Code Playgroud)

我希望这会产生一个与 lm1 %>% Augment(fit): 相同大小的数据框:newdata 中的每一行对应一行,因为 broom 将通过分组变量 cyl 和 am 匹配模型和新数据。

很遗憾,

pred1 <-
    lm1 %>%
    augment(
        fit,
        newdata = newdata)
Run Code Online (Sandbox Code Playgroud)

给了我一个 192 行(= 6 x 32)的数据框,显然每个模型都适合每一行新数据。

从别处阅读,我发现 group_by 和 rowwise 数据帧不兼容,因此 lm1 未分组,并且增加无法关联模型和新数据。是否有另一种设计模式可以让我这样做?如果它像上述尝试一样简单透明就好了,但更重要的是它可以工作。

这是我的 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  
Run Code Online (Sandbox Code Playgroud)

编辑:

@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)
Run Code Online (Sandbox Code Playgroud)

我说它看起来像你的,我的意思是我有以下列(为了一致性在这里重命名):ID (chr), attr1 (dbl), cyl (dbl), am (chr), fit (list), and data (列表)。你有 cyl、am (dbl)、fit 和 data。我将我的 am 改为 dbl,但这没有帮助。

我认为不同之处在于我在这个样本中有 3 个(ID ...类似于 mtcars 中的行名)x 2 (cyl) x 2 (am) 个单位(每个样本有 12 个测量值),而 mtcars 示例有 3 (cyl) x 2 (am) 个单元格 x 每个单元格的随机数量的汽车类型。在我的分析中,我需要查看 ID 值,但 newdata 同等适用于所有单位。如果有帮助,请将其视为测试中应用于每辆车的逆风速度。这是否暗示了增加无法处理类列表数据的原因?

编辑:将 ID 与新数据合并(使用 full=TRUE)解决了最后一个问题。我目前正在使用您提出的第一个解决方案。

aos*_*ith 5

我已经使用purrrmap2包来处理这种情况。 同时循环遍历两个列表的元素。这些列表必须具有相同的长度并且具有相同的顺序。map2

列表的元素用作您想要应用的某些函数的参数(augment在您的情况下)。这里您的两个列表将是模型列表和数据集列表(每个cyl/am组合一个列表)。

使用map2_df将结果作为 data.frame 而不是列表返回。

library(purrr)
Run Code Online (Sandbox Code Playgroud)

我制作了 data.frames 列表来使用 进行预测split。要拆分的因素的顺序决定了列表顺序,因此我确保它的顺序与 相同lm1

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

map2_df(lm1$fit, test_split, ~augment(.x, newdata = .y))
Run Code Online (Sandbox Code Playgroud)

为了避免过多担心顺序,您可以nest按组预测数据,将其连接到lm1,并将结果augment作为列表返回以进行取消嵌套。

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