Search code examples
rloopsfor-loopgtsummary

code to automate multiple glmer model and comparison for different dependent variables


I want to:

  1. to automate running the following 3 models again for different outcomes from yone, ytwo, ythree, zone, ztwo etc. (in model 1st variable should change from yone0 to ytwo0 and so on..)

2)lastly I want to export model summary and comparison to word

#I tried reading the following but do not know how to get summary and modify as per my requirement: Unexpected error when looping the glmer function with the effects package

dat3_long <- structure(list(
  ID = c("1", "1", "1", "1", "2", "2", "2", "2", "3", "3", "3", "3", "4", "4", "4", "4"),
  visit = c("0", "1", "3", "4", "0", "1", "3", "4", "0", "1", "3", "4", "0", "1", "3","4"),
  x1 = c("5.6", "1.5", "0.5", NA, "6", NA, NA, NA, "3.4","2.4", "2.5", "1", NA, 0, NA, "3.3"),
  x2 = c("0", "0", "0", "0", "1", "1", "1", "1", "2","2", "2", "2", "1", "1", "1", "1"),
  yone = c("0", NA, NA, NA, "1", "0", "0", "0", "0", "0", "0", "1", 1, NA, NA, "0"),
  yone0 = c("0", 1, NA, NA, "1", "0", "0", "0", "0", "0", "0", "1", NA, NA, 0, "0"),
  ytwo = c("0", NA, 1, NA, "1", "0", "0", "0", "0", "0", "0", "1", NA, 0, NA, "0"),
  ytwo0 = c("0", NA, NA, 1, "1", "0", "0", "0", "0", "0", "0", "1", NA, NA, NA, "0"),
  ythree = c("0", NA, 1, NA, "1", "0", "0", "0", "0", "0", "0", "1", NA, 0, NA, "0"),
  ythree0 = c("0", 1, NA, NA, "1", "0", "0", "0", "0", "0", "0", "1", NA, NA, NA, "0"),
  zone = c("0", NA, NA, 1, "1", "0", "0", "0", "0", "0", "0", "1", NA, 1, NA, "0"),
  zone0 = c("0", NA, NA, NA, "1", "0", "0", "0", "0", "0", "0", "1", NA, NA, 1, "0"),
  ztwo = c("0", NA, 1, NA, "1", "0", "0", "0", "0", "0", "0", "1", NA, NA, NA, "0"),
  ztwo0 = c("0", NA, NA, 1, "1", "0", "0", "0", "0", "0", "0", "1", NA, 1, NA, "0")),
row.names = 2:17, class = "data.frame")

#
#character to factor
names <- c(4:14)
dat3_long[,names] <- lapply(dat3_long[,names] , factor)
str(dat3_long)

# I want 1) to automate following 3 models again for different outcomes from yone, ytwo, ythree, zone, ztwo etc. (in model 1st variab changed from yone0 to ytwo 0 and so on)
#models
#1
yone_1 <- glmer(yone ~ yone0 + visit + x1+ x2 + (1 | ID), data=dat3_long, family = binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)))
summary(yone_1)
#2 + x2*visit
yone_2 <- glmer(yone ~ yone0 + visit + x1 + x2 + x2*visit +(1 | ID), data=dat3_long, family = binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)))
summary(yone_2)
#3 
yone_3 <- glmer(yone ~ yone0 + visit + x1 + x2*visit + (1 | ID), data=dat3_long, family = binomial)#, 
summary(yone_3)

#
anova(yone_1,yone_2)
#
anova(yone_3,yone_1)

#2)lastly I want to export model summary and comparison to word



#Thanking you!

Solution

  • update() and reformulate() are your friends.

    Let's say you've fitted the model once:

    base_model <- glmer(yone ~ yone0 + visit + x1+ x2 + (1 | ID), 
                    data=dat3_long, family = binomial,  ...)
    

    Set up names for the responses and predictors you want to try:

    responses <- grep("(one|two|three)$", names(dat3_long), value = TRUE)
    preds <- c("x2", "x2+x2:visit", "x2:visit")
    
    results <- list()
    for (r in responses) {
       for (p in preds) {
          form <- reformulate(response = r,
                c(paste0(r, "0"), "x1", p, "(1 | ID)"))
          results[[paste(p, r, sep = ":")]] <- update(base_model, formula = form)
       }
    }
    

    This will give you a list of fitted models; I'll let someone else explain how to export them to Word. (If I were doing it I would probably write a Quarto or RMarkdown file that looped/printed the various summaries, and compile it to a .docx file, but I guess the officer package is good for this stuff too?

    PS you should be aware that x2*visit expands to x2 + x2:visit, so x2 + x2*visit and x2*visit are equivalent as far as R is concerned. I've adjusted this in my answer.


    • if you want summaries including p-values for the fixed effects you should load library(lmerTest) instead of/after you load lme4.
    • lapply(results, summary) or library(broom.mixed); lapply(results, tidy, effects = "fixed") will get you summaries with p-values
    • to get anova() results, you can use anova() within the inner for loop (or run a separate for loop) to compare the corresponding fitted models
    • the warning method you're getting about running gtsummary indicates that you need to update your version of dplyr; install.packages("dplyr") should probably work.