CODEWITHSUNDEEP
rloopslme4mixed-models
Sergio SV
Sergio SV

Reputation: 15

How to loop over columns to evaluate different fixed effects in consecutive lme4 mixed models and extract the coefficients and P values?

I am new to R and am trying to loop a mixed model across 90 columns in a dataset. My dataset looks like the following one but has 90 predictors instead of 7 that I need to evaluate as fixed effects in consecutive models. I then need to store the model output (coefficients and P values) to finally construct a figure summarizing the size effects of each predictor. I know the discussion of P value estimates from lme4 mixed models.

For example:

set.seed(101)

mydata <- tibble(id = rep(1:32, times=25), 
   time = sample(1:800), 
   experiment = rep(1:4, times=200),
   Y = sample(1:800), 
   predictor_1 = runif(800), 
   predictor_2 = rnorm(800), 
   predictor_3 = sample(1:800), 
   predictor_4 = sample(1:800),
   predictor_5 = seq(1:800),
   predictor_6 = sample(1:800),
   predictor_7 = runif(800)) %>% arrange (id, time)

The model to iterate across the N predictors is:

library(lme4)
library(lmerTest)  # To obtain new values

mixed.model <- lmer(Y ~ predictor_1 + time + (1|id) + (1|experiment), data = mydata)

summary(mixed.model)

My coding skills are far from being able to set a loop to repeat the model across the N predictors in my dataset and store the coefficients and P values in a dataframe.

I have been able to iterate across all the predictors fitting linear models instead of mixed models using lapply. But I have failed to apply this strategy with mixed models.

varlist <- names(mydata)[5:11]

lm_models <- lapply(varlist, function(x) {
  lm(substitute(Y ~ i, list(i = as.name(x))), data = mydata)
})

Upvotes: 1

Views: 403

Answers (1)

jay.sf
jay.sf

Reputation: 72683

One option is to update the formula of a restricted model (w/o predictor) in an lapply loop over the predictors. Then summaryze the resulting list and subset the coefficient matrix using a Vectorized function.

library(lmerTest)
mixed.model <- lmer(Y ~ time + (1|id) + (1|experiment), data = mydata)

preds <- grep('pred', names(mydata), value=TRUE)
fits <- lapply(preds, \(x) update(mixed.model, paste('. ~ . + ', x)))

extract_coef_p <- Vectorize(\(x) x |> summary() |> coef() |> {\(.) .[3, c(1, 5)]}())

res <- `rownames<-`(t(extract_coef_p(fits)), preds)
res
#                 Estimate  Pr(>|t|)
# predictor_1 -7.177579138 0.8002737
# predictor_2 -5.010342111 0.5377551
# predictor_3 -0.013030513 0.7126500
# predictor_4 -0.041702039 0.2383835
# predictor_5 -0.001437124 0.9676346
# predictor_6  0.005259293 0.8818644
# predictor_7 31.304496255 0.2511275

Upvotes: 1

Related Questions