Quick reprex for @ jess_carilli

same_model_many_responses.R

#some libraries
library(dplyr)
library(tidyr)
library(broom)
library(purrr)

#let's use the npk data
head(npk)

#but make up a few extra columns
set.seed(2019)
npk <- npk %>%
  mutate(yield2 = rnorm(n()),
         some_variable = rnorm(n()))

head(npk)

#now, let's automate some modeling

npk_aov <- npk %>%
  #reshape so all responses are in the same column
  gather(response_type, value, yield:some_variable) %>%
  
  #group by response
  group_by(response_type) %>%
  
  #nest the data
  nest() %>%
  
  #fit models for each response
  mutate(mod = map(data, ~lm(value ~ block + N * P + K, data = .))) %>%
           
  #extract the anova info
  mutate(anova_tables = map(mod, car::Anova)) %>%
  
  #make it pretty
  mutate(nice_anova = map(anova_tables, tidy)) %>%
  
  #unnest!
  unnest(nice_anova)

#what did we get?
npk_aov


#let's just look at the N effect
npk_aov %>% filter(term == "N")

#what had a p value < 0.05
npk_aov %>% filter(p.value <= 0.05)