r - 如何在 R 中运行多元回归时缩短 broom 中的代码？

Question

我正在按不同的组进行多重回归。我想让事情自动化一点。我最初尝试运行并保存回归模型 1、模型 2 和模型 3。然后我尝试将代码缩短如下:

 temp <- df %>%
      group_by(group) %>%
      do(model1, model2, model3, data = .))) %>%   
      gather(model_name, model, -group) %>%                        
      unnest() 

但这没有用。您可以在下面找到有效的长版本。有人可以指导我如何缩短它吗？

df <- tibble(
  a = rnorm(1000),
  b = rnorm(1000),
  c = rnorm(1000),
  d = rnorm(1000),
  group =sample.int(300,size=1000,replace=TRUE)-1)
)

df$group = as.factor(df$group)

temp1 <- df %>%
  group_by(group) %>%
  do(model2 = tidy(lm(a ~ b , data = .))) %>%   
  gather(model_name, model, -group) %>%                        
  unnest() 

temp2 <- df %>%
  group_by(group) %>%
  do(model2 = tidy(lm(a ~ c , data = .))) %>%   
  gather(model_name, model, -group) %>%                        
  unnest() 

temp3 <- df %>%
  group_by(group) %>%
  do(model3 = tidy(lm(a ~ d , data = .))) %>%   
  gather(model_name, model, -group) %>%                        
  unnest() 

Answer 1

这可能有效，使用 nest_by 和 purrr 中的 map。

尝试使用 nest_by(dplyr 版本 1.0.0)代替 group_by，并在嵌套数据的每一行上运行您的模型。使用 nest_by 将创建一个新的临时列表列 data。它类似于以前使用 nest 和 rowwise。该模型也需要在此处的列表中。

使用 map，您可以为每个自变量“b”、“c”和“d”建立模型。自变量也作为单独的列包含在内(也可以是特定模型的标签)。

library(tidyverse)
library(purrr)
library(broom)

df %>%
  nest_by(group) %>%
  mutate(model = list(map(c("b", "c", "d"), ~
                       cbind(independent = .x, 
                             tidy(lm(formula(paste0("a ~ ", .x)), data = data)))))) %>%
  summarise(bind_rows(model))

输出

   group independent term        estimate std.error statistic p.value
   <fct> <chr>       <chr>          <dbl>     <dbl>     <dbl>   <dbl>
 1 0     b           (Intercept)   0.0480   NaN       NaN     NaN    
 2 0     b           b             0.268    NaN       NaN     NaN    
 3 0     c           (Intercept)  -0.124    NaN       NaN     NaN    
 4 0     c           c            -0.447    NaN       NaN     NaN    
 5 0     d           (Intercept)  -0.107    NaN       NaN     NaN    
 6 0     d           d             0.377    NaN       NaN     NaN    
 7 1     b           (Intercept)   0.473      0.296     1.60    0.356
 8 1     b           b             0.383      0.261     1.47    0.380
 9 1     c           (Intercept)   0.547      0.544     1.01    0.498
10 1     c           c            -0.183      0.798    -0.229   0.857

编辑 (12/19/20):如果您想包含具有多个协变量和交互项的模型，您只需在字符串中提供公式即可。

例如，如果您想为每个 group 运行 3 个模型:

• 主效应“b”和“c”以及交互项“b*c”
• 主要影响“c”和“d”
• 主要影响“d”

您可以执行以下操作:

df %>%
  nest_by(group) %>%
  mutate(model = list(map(c("b + c + b*c", "c + d", "d"), ~
                       cbind(model = .x, 
                             tidy(lm(formula(paste0("a ~ ", .x)), data = data)))))) %>%
  summarise(bind_rows(model))

输出

   group model       term        estimate std.error statistic  p.value
   <fct> <chr>       <chr>          <dbl>     <dbl>     <dbl>    <dbl>
 1 0     b + c + b*c (Intercept)   0.718      0.281    2.56     0.0835
 2 0     b + c + b*c b             0.819      0.348    2.35     0.100 
 3 0     b + c + b*c c            -0.351      0.315   -1.11     0.346 
 4 0     b + c + b*c b:c           0.0444     0.461    0.0964   0.929 
 5 0     c + d       (Intercept)   0.614      0.409    1.50     0.208 
 6 0     c + d       c            -0.271      0.439   -0.618    0.570 
 7 0     c + d       d             0.182      0.487    0.374    0.727 
 8 0     d           (Intercept)   0.605      0.383    1.58     0.175 
 9 0     d           d             0.208      0.455    0.456    0.667 
10 1     b + c + b*c (Intercept)   0.590    NaN      NaN      NaN     
# … with 2,600 more rows

或者，如果您愿意，您可以像这样完整且单独地列出方程式:

my_models <- c(
  "a ~ b + c + b*c", 
  "a ~ c + d", 
  "a ~ d"
)

df %>%
  nest_by(group) %>%
  mutate(model = list(map(my_models, ~
                       cbind(model = .x, 
                             tidy(lm(formula(.x), data = data)))))) %>%
  summarise(bind_rows(model))

输出

   group model           term        estimate std.error statistic  p.value
   <fct> <chr>           <chr>          <dbl>     <dbl>     <dbl>    <dbl>
 1 0     a ~ b + c + b*c (Intercept)   0.718      0.281    2.56     0.0835
 2 0     a ~ b + c + b*c b             0.819      0.348    2.35     0.100 
 3 0     a ~ b + c + b*c c            -0.351      0.315   -1.11     0.346 
 4 0     a ~ b + c + b*c b:c           0.0444     0.461    0.0964   0.929 
 5 0     a ~ c + d       (Intercept)   0.614      0.409    1.50     0.208 
 6 0     a ~ c + d       c            -0.271      0.439   -0.618    0.570 
 7 0     a ~ c + d       d             0.182      0.487    0.374    0.727 
 8 0     a ~ d           (Intercept)   0.605      0.383    1.58     0.175 
 9 0     a ~ d           d             0.208      0.455    0.456    0.667 
10 1     a ~ b + c + b*c (Intercept)   0.590    NaN      NaN      NaN     
# … with 2,600 more rows

数据

set.seed(123)

df <- tibble(
  a = rnorm(1000),
  b = rnorm(1000),
  c = rnorm(1000),
  d = rnorm(1000),
  group =sample.int(300,size=1000,replace=TRUE)-1)
)

df$group = as.factor(df$group)