How to write a for loop which creates a model and has a function which references that same model

user_9829297 0
Viewed 386

I am trying to run a post hoc analysis on an unbalanced two way anova using the anova_test funciton in the rstatix package. I need to run this post hoc test iteratively, as I have ~26 response (y) variables. My first step is to create models of all my y variables with relation to group and treatment. I have successfully managed to do this, creating a single list with 26 models:

models <- map(data[,y1:y26], ~(lm(.x ~data$group*data$treatment)))

Now comes the part I'm stuck on. Referring to these models iteratively. I would like to run the following code for every y variable I have:

group_by(group) %>%
anova_test(y ~ treatment, error = models(y), type = 3)

where my y changes every time and as it does, the "model" (referred to in the error = term) is updated accordingly. I'm struggling with this bit since first set of models I make is used to inform the second set of models.

However, if I run just one y variable through this whole bit of code at one time, I get the appropriate results. 

model <- lm(y ~ group*treatment, data = data)

data %>%
group_by(group) %>%
anova_test(y ~ treatment, error = model, type = 3)

I have tried creating a for loop as well as using the map function in the purrr package but I have been unsuccessful. I am new to for loops and purrr so I am sure it's a simple fix I just can't see it.

Basically I want a way to run 

data %>%
group_by(group) %>%
anova_test(y ~ treatment, error = model, type = 3)

iteratively for different y variables (y1, y2, ..., y26) while also referring to the approprite model (model$y1, model$y2, ..., model$26).

Thanks for your help!

purrr anova r for-loop
2 Answers

Well you didn't give any data so let's use toothgrowth. You seem to like the model format, so let's build a list of models. You could do this in an automated fashion but to make it clear lets do it by hand. The call purrr::map with the anova_test function. You'll get a list back. Since you're in charge of naming the list elements go to town.

Updated answer May 18th. Now using map2 since you want two different models passed build a list for each...

library(rstatix)
library(purrr)

ToothGrowth$len2 <- ToothGrowth$len^2 # for variety

models <- list(model1 = lm(len ~ supp*dose, ToothGrowth), 
               model2 = lm(len ~ dose*supp, ToothGrowth),
               model3 = lm(len2 ~ dose*supp, ToothGrowth),
               model4 = lm(len2 ~ supp*dose, ToothGrowth))

models2 <- list(model1 = lm(len ~ supp, ToothGrowth), 
               model2 = lm(len ~ dose, ToothGrowth),
               model3 = lm(len2 ~ dose, ToothGrowth),
               model4 = lm(len2 ~ supp, ToothGrowth))


# one model
purrr::map(models, ~ anova_test(.x, type = 3))


# now with model for error term
purrr::map2(models, models2, ~ anova_test(.x, error = .y, type = 3))
#> Coefficient covariances computed by hccm()
#> Coefficient covariances computed by hccm()
#> Coefficient covariances computed by hccm()
#> Coefficient covariances computed by hccm()
#> $model1
#> ANOVA Table (type III tests)
#> 
#>      Effect DFn DFd      F        p p<.05   ges
#> 1      supp   1  58  4.058 0.049000     * 0.065
#> 2      dose   1  58 12.717 0.000734     * 0.180
#> 3 supp:dose   1  58  1.588 0.213000       0.027
#> 
#> $model2
#> ANOVA Table (type III tests)
#> 
#>      Effect DFn DFd      F        p p<.05   ges
#> 1      dose   1  58 33.626 2.92e-07     * 0.367
#> 2      supp   1  58 10.729 2.00e-03     * 0.156
#> 3 dose:supp   1  58  4.200 4.50e-02     * 0.068
#> 
#> $model3
#> ANOVA Table (type III tests)
#> 
#>      Effect DFn DFd      F        p p<.05   ges
#> 1      dose   1  58 36.028 1.35e-07     * 0.383
#> 2      supp   1  58  7.128 1.00e-02     * 0.109
#> 3 dose:supp   1  58  2.709 1.05e-01       0.045
#> 
#> $model4
#> ANOVA Table (type III tests)
#> 
#>      Effect DFn DFd      F        p p<.05   ges
#> 1      supp   1  58  2.684 0.107000       0.044
#> 2      dose   1  58 13.566 0.000508     * 0.190
#> 3 supp:dose   1  58  1.020 0.317000       0.017
user_91900340
user_98292970
Related