Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Quick reprex for @ jess_carilli
#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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.