Use broom and tidyverse to run regressions on diff

Question

I'm looking for a Tidyverse / broom solution that can solve this puzzle:

Let's say I have different DVs and a specific set of IVS and I want to perform a regression that considers every DV and this specific set of IVs. I know I can use something like for i in or apply family, but I really want to run that using tidyverse.

The following code works as an example

ds <- data.frame(income = rnorm(100, mean=1000,sd=200),
                 happiness = rnorm(100, mean = 6, sd=1),
                 health = rnorm(100, mean=20, sd = 3),
                 sex = c(0,1),
                 faculty = c(0,1,2,3))

mod1 <- lm(income ~ sex + faculty, ds)
mod2 <- lm(happiness ~ sex + faculty, ds)
mod3 <- lm(health ~ sex + faculty, ds)
summary(mod1)
summary(mod2)
summary(mod3)

Income, happiness, and health are DVs. Sex and Faculty are IVs and they will be used for all regressions.

That was the closest I found

Let me know If I need to clarify my question. Thanks.

Answer 1

Another method is to gather the dependent variables and use a grouped data frame to fit the models with do. This is the method explained in the broom and dplyr vignette.

library(tidyverse)
library(broom)
ds <- data.frame(
  income = rnorm(100, mean = 1000, sd = 200),
  happiness = rnorm(100, mean = 6, sd = 1),
  health = rnorm(100, mean = 20, sd = 3),
  sex = c(0, 1),
  faculty = c(0, 1, 2, 3)
)
ds %>%
  gather(dv_name, dv_value, income:health) %>%
  group_by(dv_name) %>%
  do(tidy(lm(dv_value ~ sex + faculty, data = .)))
#> # A tibble: 9 x 6
#> # Groups:   dv_name [3]
#>   dv_name   term         estimate std.error statistic  p.value
#>   <chr>     <chr>           <dbl>     <dbl>     <dbl>    <dbl>
#> 1 happiness (Intercept)     6.25      0.191    32.7   3.14e-54
#> 2 happiness sex             0.163     0.246     0.663 5.09e- 1
#> 3 happiness faculty        -0.172     0.110    -1.56  1.23e- 1
#> 4 health    (Intercept)    20.1       0.524    38.4   1.95e-60
#> 5 health    sex             0.616     0.677     0.909 3.65e- 1
#> 6 health    faculty        -0.653     0.303    -2.16  3.36e- 2
#> 7 income    (Intercept)  1085.       32.8      33.0   1.43e-54
#> 8 income    sex           -12.9      42.4      -0.304 7.62e- 1
#> 9 income    faculty       -25.1      19.0      -1.32  1.89e- 1

Created on 2018-08-01 by the reprex package (v0.2.0).

Answer 2

As you have different dependent variables but the same independent, you can form a matrix of these and pass to lm.

mod = lm(cbind(income, happiness, health) ~ sex + faculty, ds)

And I think broom::tidy works

library(broom)
tidy(mod)

#    response        term      estimate  std.error  statistic      p.value
# 1    income (Intercept) 1019.35703873 31.0922529 32.7849205 2.779199e-54
# 2    income         sex  -54.40337314 40.1399258 -1.3553431 1.784559e-01
# 3    income     faculty   19.74808081 17.9511206  1.1001030 2.740100e-01
# 4 happiness (Intercept)    5.97334562  0.1675340 35.6545278 1.505026e-57
# 5 happiness         sex    0.05345555  0.2162855  0.2471528 8.053124e-01
# 6 happiness     faculty   -0.02525431  0.0967258 -0.2610918 7.945753e-01
# 7    health (Intercept)   19.76489553  0.5412676 36.5159396 1.741411e-58
# 8    health         sex    0.32399380  0.6987735  0.4636607 6.439296e-01
# 9    health     faculty    0.10808545  0.3125010  0.3458723 7.301877e-01

Answer 3

We can loop through the column names that are dependent variables, use paste to create the formula to be passed into lm and get the summary statistics with tidy (from broom)

library(tidyverse)
library(broom)
map(names(ds)[1:3], ~ 
            lm(formula(paste0(.x, "~", 
                paste(names(ds)[4:5], collapse=" + "))), data = ds) %>%
               tidy)

---

If we want it in a single data.frame with a column identifier for dependent variable,

map_df(set_names(names(ds)[1:3]), ~ 
      lm(formula(paste0(.x, "~", 
        paste(names(ds)[4:5], collapse=" + "))), data = ds) %>%
   tidy, .id = "Dep_Variable")