How do I run a mixed linear regression model on several outcomes variables and get presentable results?

I finally gave up and admitted I need help. I have this data set with 3 different groups, measured at 2 time points and 49 outcome variables. I would like to do a mixed linear regression analysis on each outcome variable for within group change between time points. As shown in table below:

Naturally, I’d rather not do all the 147 analysis manually (even though at this stage it would have saved me a lot of time).

So after scouring forums for answers this is what I’ve tried so far:

library(lme4)
library(lmerTest)
library(tidyverse)

df <- data.frame(
  id = rep(1:66, each = 2),
  visit = 0:1,
  rand = rep(0:2, each = 2),
  x1 = sample(4000:9000, 132),
  x2 = sample(1200:3400, 132),
  x3 = sample(220:400, 132)
)

df_rand0 <- df %>%
  filter(rand == "0")
df_rand1 <- df %>%
  filter(rand == "1")
df_rand2 <- df %>%
  filter(rand == "2")

depVarList <- colnames(df_rand0[4:6])
allModels <- lapply(depVarList, function(x){
  lmer(formula = paste0("`", x, "` ~ visit + (1| id)"),
       data = df_rand0, na.action = na.omit)
})

Which does generate a list of results but I’m missing p-values and with 49 variables it generates a large list. I would like to get a better overview as well as get the p-values from the tests. I tried loading the tidymodels package and running tidy() but it returns "Error: No tidy method recognized for this list."

>Solution :

broom.mixed provides tidiers for mixed models. You can also use purrr::map_dfr() instead of lapply() to get all coefficients in one dataframe.

library(lme4)
library(lmerTest)
library(tidyverse)
library(broom.mixed)
set.seed(13)

allModels <- map_dfr(
  set_names(depVarList), 
  \(x) {
    tidy(lmer(
       formula = paste0(x, " ~ visit + (1| id)"),
       data = df_rand0, 
       na.action = na.omit
   ))
  },
  .id = "dv"
)

allModels

# A tibble: 12 × 9
   dv    effect   group    term          estim…¹ std.e…² stati…³    df   p.value
   <chr> <chr>    <chr>    <chr>           <dbl>   <dbl>   <dbl> <dbl>     <dbl>
 1 x1    fixed    <NA>     (Intercept)   6372.     286.   22.3    32.9  2.00e-21
 2 x1    fixed    <NA>     visit          229.     278.    0.821  21.0  4.21e- 1
 3 x1    ran_pars id       sd__(Interce…  973.      NA    NA      NA   NA       
 4 x1    ran_pars Residual sd__Observat…  923.      NA    NA      NA   NA       
 5 x2    fixed    <NA>     (Intercept)   2278.     123.   18.5    42.0  8.84e-22
 6 x2    fixed    <NA>     visit          -19.2    174.   -0.110  42.0  9.13e- 1
 7 x2    ran_pars id       sd__(Interce…    0       NA    NA      NA   NA       
 8 x2    ran_pars Residual sd__Observat…  578.      NA    NA      NA   NA       
 9 x3    fixed    <NA>     (Intercept)    314.      12.1  26.0    42.0  1.63e-27
10 x3    fixed    <NA>     visit           -8.82    17.1  -0.516  42.0  6.09e- 1
11 x3    ran_pars id       sd__(Interce…    0       NA    NA      NA   NA       
12 x3    ran_pars Residual sd__Observat…   56.7     NA    NA      NA   NA       
# … with abbreviated variable names ¹​estimate, ²​std.error, ³​statistic