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:
| Id | rand | visit | x1 | x2 | ... |
|---|---|---|---|---|---|
| 1 | 0 | 0 | 178 | 5,2 | |
| 2 | 0 | 0 | 165 | NA | |
| 3 | 2 | 0 | 142 | 1,3 | |
| 4 | 1 | 0 | 198 | 2,7 | |
| 1 | 0 | 1 | 191 | 9,5 | |
| 2 | 0 | 1 | 183 | 3,9 |
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."
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