如何使用 broom 和 dplyr 将分组数据应用于分组模型?
How can I apply grouped data to grouped models using broom and dplyr?
我想做相当于将 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 = .))
正如预期的那样,我得到了一个包含 6 行的按行数据框。
此图确认有六个组:
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 中的每一行。
现在是挑战:我想使用新数据获得拟合输出,其中我将 wt 增加 cyl/4:
newdata <-
mtcars2 %>%
mutate(
wt = wt + cyl/4)
我希望这会产生一个与 lm1 %>% augment(fit): newdata 中的每一行对应一行大小的数据框,因为 broom 将通过分组变量 cyl 和 newdata 匹配模型和 newdata上午
不幸的是,
pred1 <-
lm1 %>%
augment(
fit,
newdata = newdata)
给我一个包含 192 行 (= 6 x 32) 的数据框,显然每个模型都适合新数据的每一行。
从其他地方阅读,我了解到 group_by 和行向数据帧不兼容,因此 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 命令中遇到了一个问题:它 returns "Error: augment doesn't know how to deal with data of class list"。
我的真实代码更像是:
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 (list ) 和数据(列表)。您有 cyl、am (dbl)、fit 和数据。我将我的 am 更改为 dbl,但这没有帮助。
我认为不同之处在于我在这个样本中有 3 (ID ... 类似于 mtcars 中的行名) x 2 (cyl) x 2 (am) 个单位(每个样本有 12 个测量值),而mtcars 示例有 3 (cyl) x 2 (am) 个单元 x 每个单元随机数量的汽车类型。在我的分析中,我需要查看 ID 值,但 newdata 同样适用于所有单元。如果有帮助,请将其视为应用于测试中每辆汽车的逆风速度。这是否表明 augment 抱怨它无法处理 class 列表的数据?
编辑:将 ID 与新数据合并(使用 full=TRUE)解决了最后一个问题。我目前正在使用您提出的第一个解决方案。
对于这种情况,我使用了包 purrr 中的 map2
。 map2
同时循环遍历两个列表的元素。列表的长度和顺序必须相同。
列表的元素用作您要应用的某些函数的参数(augment
,在您的例子中)。这里你的两个列表将是一个模型列表和一个数据集列表(每个 cyl
/am
组合一个列表)。
使用 map2_df
return 将结果作为 data.frame 而不是列表。
library(purrr)
我制作了 data.frame 的列表以使用 split
进行预测。要拆分的因素的顺序决定了列表顺序,因此我确保它与 lm1
的顺序相同。
test_split = split(newdata, list(newdata$am, newdata$cyl)
map2_df(lm1$fit, test_split, ~augment(.x, newdata = .y))
为了避免太担心顺序,您可以按组 nest
预测数据,将其加入 lm1
,return augment
的结果作为取消嵌套的列表。
newdata %>%
group_by(cyl, am) %>%
nest() %>%
inner_join(lm1, .) %>%
mutate(pred = list(augment(fit, newdata = data))) %>%
unnest(pred)
我想做相当于将 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 = .))
正如预期的那样,我得到了一个包含 6 行的按行数据框。
此图确认有六个组:
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 中的每一行。
现在是挑战:我想使用新数据获得拟合输出,其中我将 wt 增加 cyl/4:
newdata <-
mtcars2 %>%
mutate(
wt = wt + cyl/4)
我希望这会产生一个与 lm1 %>% augment(fit): newdata 中的每一行对应一行大小的数据框,因为 broom 将通过分组变量 cyl 和 newdata 匹配模型和 newdata上午
不幸的是,
pred1 <-
lm1 %>%
augment(
fit,
newdata = newdata)
给我一个包含 192 行 (= 6 x 32) 的数据框,显然每个模型都适合新数据的每一行。
从其他地方阅读,我了解到 group_by 和行向数据帧不兼容,因此 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 命令中遇到了一个问题:它 returns "Error: augment doesn't know how to deal with data of class list"。
我的真实代码更像是:
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 (list ) 和数据(列表)。您有 cyl、am (dbl)、fit 和数据。我将我的 am 更改为 dbl,但这没有帮助。
我认为不同之处在于我在这个样本中有 3 (ID ... 类似于 mtcars 中的行名) x 2 (cyl) x 2 (am) 个单位(每个样本有 12 个测量值),而mtcars 示例有 3 (cyl) x 2 (am) 个单元 x 每个单元随机数量的汽车类型。在我的分析中,我需要查看 ID 值,但 newdata 同样适用于所有单元。如果有帮助,请将其视为应用于测试中每辆汽车的逆风速度。这是否表明 augment 抱怨它无法处理 class 列表的数据?
编辑:将 ID 与新数据合并(使用 full=TRUE)解决了最后一个问题。我目前正在使用您提出的第一个解决方案。
对于这种情况,我使用了包 purrr 中的 map2
。 map2
同时循环遍历两个列表的元素。列表的长度和顺序必须相同。
列表的元素用作您要应用的某些函数的参数(augment
,在您的例子中)。这里你的两个列表将是一个模型列表和一个数据集列表(每个 cyl
/am
组合一个列表)。
使用 map2_df
return 将结果作为 data.frame 而不是列表。
library(purrr)
我制作了 data.frame 的列表以使用 split
进行预测。要拆分的因素的顺序决定了列表顺序,因此我确保它与 lm1
的顺序相同。
test_split = split(newdata, list(newdata$am, newdata$cyl)
map2_df(lm1$fit, test_split, ~augment(.x, newdata = .y))
为了避免太担心顺序,您可以按组 nest
预测数据,将其加入 lm1
,return augment
的结果作为取消嵌套的列表。
newdata %>%
group_by(cyl, am) %>%
nest() %>%
inner_join(lm1, .) %>%
mutate(pred = list(augment(fit, newdata = data))) %>%
unnest(pred)