Using purrr map for several linear regressions for each predictor in df

user_12746690 0
Viewed 623

I'm trying to run many individual linear regressions with one Y variable and many x variables. My data has 300+ x variables. I've been trying to do this with purrr and broom but cannot figure out how to get the output how I'd like.

Example:

iris <- iris %>% 
  select_if(is.numeric)

iris %>% 
  map(~lm(Sepal.Length ~ .x, data = iris)) %>% 
  map(summary) %>% 
  map_df(tidy)

This produces the following output:

# A tibble: 6 x 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 (Intercept)    0      3.79e-17   0.      1.00e+  0
2 .x             1      6.43e-18   1.56e17 0.       
3 (Intercept)    6.53   4.79e- 1   1.36e 1 6.47e- 28
4 .x            -0.223  1.55e- 1  -1.44e 0 1.52e-  1
5 (Intercept)    4.31   7.84e- 2   5.49e 1 2.43e-100
6 .x             0.409  1.89e- 2   2.16e 1 1.04e- 47

Which is close to what I'm looking for, but not quite! I want the variable names to in the 'term' column here and I don't want the intercept pasted for each model. The results I'm looking for a more like: 

# A tibble: 6 x 5
  term        estimate std.error statistic   p.value
  <chr>          <dbl>     <dbl>     <dbl>     <dbl>
1 Sepal.Width    0      3.79e-17   0.      1.00e+  0
2 Petal.Width    1      6.43e-18   1.56e17 0.       
3 Petal.Length   6.53   4.79e- 1   1.36e 1 6.47e- 28

Any help getting to that point would be greatly appreciated!! And of course extra appreciation for explanations of the process (i'm learning)

Cheers

broom purrr tidyverse regression r
1 Answers

lm() follows a special set of rules called non-standard evaluation, where it uses parts of the expression in its computation. Here's a simple example showing the difference:

a <- "purrr"

print(a)        # Standard evaluation - expression a is evaluated to its value
# [1] "purrr"

library(a)      # Non-standard evaluation - expression a is used as-is
# Error in library(a) : there is no package called ‘a’

Similarly, lm() uses parts of the expression Sepal.Length ~ .x, which is why you see .x in your output, instead of the value that map() places into .x. Here's a couple of options to solve this.

Option 1: Construct and evaluate expressions "by hand"

colnames(iris) %>%                                # Start with all column names
    setdiff( "Sepal.Length" ) %>%                 # ...that are not Sepal.Length
    rlang::syms() %>%                             # Convert them to symbols
    map( ~rlang::expr(lm(Sepal.Length ~ !!.x,
                         data=iris)) ) %>%        # Create expressions
    map( eval.parent ) %>%                        # Evaluate expressions
    map( broom::tidy ) %>%                        # Tidy up the output
    bind_rows() %>%                               # Combine into a single data frame
    filter( term != "(Intercept)" )               # Drop all (Intercept) entries

Here, !!.x will replace .x with the symbol stored inside of it. The output of this map() step will be a set of expressions that look exactly like what you want:

lm(Sepal.Length ~ Sepal.Width, data = iris)
lm(Sepal.Length ~ Petal.Length, data = iris)
lm(Sepal.Length ~ Petal.Width, data = iris)

Option 2: Annotate the rows yourself

iris %>%
    select( -Sepal.Length ) %>%                   
    map( ~lm(Sepal.Length ~ .x, data=iris) ) %>%    # As before
    map( broom::tidy ) %>%                          # Tidy up the output
    map( filter, term != "(Intercept)" ) %>%        # Remove (Intercept) entries
    map( select, -term ) %>%                        # Remove the default term column
    bind_rows( .id="term" )                         # Make your own from the list names
user_3001870
Related