在 R 中迭代多个回归模型和数据子集

Question

我正在尝试学习如何使用 R 中的 purrr 和 broom 包在数据集的子集上自动运行 3 个或更多回归模型。我正在使用 nest %>% mutate(map()) %>% unnest() 牢记在心。

当只有一个回归模型应用于多个数据子集时，我能够在线复制示例。但是，当我的函数中有多个回归模型时，我遇到了问题。

我想做什么

library(tidyverse)
library(broom)

estimate_model <- function(df) {
  model1 <- lm(mpg ~ wt, data = df)
  model2 <- lm(mpg ~ wt + gear, data = df)
  model3 <- lm(mpg ~ wt + gear + vs, data = df)
}

ols_1dep_3specs <- mtcars %>%
    nest(-cyl) %>%
    mutate(
       estimates = map(data, estimate_model), # want to run several models at once
       coef_wt = map(estimate,  ~pluck(coef(.), "wt")), # coefficient of wt only
       se_wt = map(estimate, ~pluck(tidy(.), "std.error")[[2]]), # se of wt only
       rsq = map(model, ~pluck(glance(.), "r.squared")),
       arsq = map(model, ~pluck(glance(.), "adj.r.squared"))
    ) %>%
    unnest(coef_wt, se_wt, rsq, arsq)

ols_1dep_3specs

不幸的是，这似乎只在函数estimate_model 仅包含一个回归模型时才有效。关于在有多个规范时如何编写代码的任何建议？对 nest() %>% mutate(map()) %>% nest() 框架之外的建议开放。

---

以下代码可以达到我希望达到的目的，但它涉及大量重复。

estimate_model1 <- function(df) {
  model1 <- lm(mpg ~ wt, data = df)
}
estimate_model2 <- function(df) {
  model2 <- lm(mpg ~ wt + gear, data = df)
}
estimate_model3 <- function(df) {
  model3 <- lm(mpg ~ wt + gear + vs, data = df)
}

ols_1dep_3specs <- mtcars %>%
  nest(-cyl) %>%
  mutate(model_1 = map(data, estimate_model1),
         model_2 = map(data, estimate_model2),
         model_3 = map(data, estimate_model3)) %>%
  mutate(coef_wt_1 = map_dbl(model_1, ~pluck(coef(.), "wt")),
         coef_wt_2 = map_dbl(model_2, ~pluck(coef(.), "wt")),
         coef_wt_3 = map_dbl(model_3, ~pluck(coef(.), "wt")),
         rsq_1 = map_dbl(model_1, ~pluck(glance(.), "r.squared")),
         rsq_2 = map_dbl(model_2, ~pluck(glance(.), "r.squared")),
         rsq_3 = map_dbl(model_3, ~pluck(glance(.), "r.squared"))) %>% 
  dplyr::select(starts_with("coef_wt"), starts_with("rsq")) 

Answer 1

在函数中，没有return调用，最好把所有模型放在一个list中

estimate_model <- function(df) {
        model1 <- lm(mpg ~ wt, data = df)
        model2 <- lm(mpg ~ wt + gear, data = df)
        model3 <- lm(mpg ~ wt + gear + vs, data = df)
        list(model1, model2, model3)
      }

然后通过循环遍历每个 list 元素来应用第一段代码

mtcars %>% 
  group_by(cyl) %>%
  nest() %>% 
  mutate(estimates = map(data, estimate_model),
         coef_wt = map(estimates,  ~map_dbl(.x, ~ pluck(coef(.x), "wt"))),
         se_wt = map(estimates, ~map_dbl(.x, ~pluck(tidy(.x), "std.error")[[2]])), 
         rsq = map(estimates, ~ map_dbl(.x, ~pluck(glance(.x), "r.squared"))),
          arsq = map(estimates, ~map_dbl(.x, ~ pluck(glance(.x), "adj.r.squared")))) %>%
  unnest(c(coef_wt, se_wt, rsq, arsq))
# A tibble: 9 x 7
# Groups:   cyl [3]
#    cyl            data estimates  coef_wt se_wt   rsq  arsq
#  <dbl> <list<df[,10]>> <list>       <dbl> <dbl> <dbl> <dbl>
#1     6        [7 × 10] <list [3]>   -2.78 1.33  0.465 0.357
#2     6        [7 × 10] <list [3]>   -3.92 1.41  0.660 0.489
#3     6        [7 × 10] <list [3]>   -6.19 4.49  0.690 0.379
#4     4       [11 × 10] <list [3]>   -5.65 1.85  0.509 0.454
#5     4       [11 × 10] <list [3]>   -5.38 2.08  0.517 0.396
#6     4       [11 × 10] <list [3]>   -5.13 2.16  0.555 0.364
#7     8       [14 × 10] <list [3]>   -2.19 0.739 0.423 0.375
#8     8       [14 × 10] <list [3]>   -2.43 0.798 0.459 0.361
#9     8       [14 × 10] <list [3]>   -2.43 0.798 0.459 0.361

Answer 2

使用purrr::lst 会自动根据其元素命名该列表，这有助于您以后跟踪模型。将函数应用到嵌套数据后，您可以拉出一列模型名称。

我选择在工作流程的早期将 plucking 替换为取消嵌套，并使用 2 个映射调用将值作为向量而不是列表获取。这只是一种偏好，但是当列嵌套的深度较低时，我会更轻松。

library(tidyverse)
library(broom)

estimate_model <- function(df) {
  model1 <- lm(mpg ~ wt, data = df)
  model2 <- lm(mpg ~ wt + gear, data = df)
  model3 <- lm(mpg ~ wt + gear + vs, data = df)
  lst(model1, model2, model3)
}

mtcars %>%
  group_by(cyl) %>%
  nest() %>%
  mutate(mods = map(data, estimate_model),
         mod_id = map(mods, names)) %>%
  unnest(c(mod_id, mods)) %>%
  mutate(coef_wt = map(mods, coef) %>% map_dbl("wt"),
         se_wt = map(mods, tidy) %>% map("std.error") %>% .[[2]],
         rsq = map(mods, glance) %>% map_dbl("r.squared"),
         arsq = map(mods, glance) %>% map_dbl("adj.r.squared"))
#> # A tibble: 9 x 8
#> # Groups:   cyl [3]
#>     cyl            data mods   mod_id coef_wt  se_wt   rsq  arsq
#>   <dbl> <list<df[,10]>> <list> <chr>    <dbl>  <dbl> <dbl> <dbl>
#> 1     6        [7 × 10] <lm>   model1   -2.78  6.36  0.465 0.357
#> 2     6        [7 × 10] <lm>   model2   -3.92  1.41  0.660 0.489
#> 3     6        [7 × 10] <lm>   model3   -6.19  0.727 0.690 0.379
#> 4     4       [11 × 10] <lm>   model1   -5.65 11.6   0.509 0.454
#> 5     4       [11 × 10] <lm>   model2   -5.38  2.08  0.517 0.396
#> 6     4       [11 × 10] <lm>   model3   -5.13  2.20  0.555 0.364
#> 7     8       [14 × 10] <lm>   model1   -2.19  4.91  0.423 0.375
#> 8     8       [14 × 10] <lm>   model2   -2.43  0.798 0.459 0.361
#> 9     8       [14 × 10] <lm>   model3   -2.43  0.835 0.459 0.361