1

I want to do the same thing as described here but for glmer models: How do I store lm object in a data frame in R

example model

fit1 <- glmer(y ~ scale(x) + (1 | ID), data=dat, family = binomial, control=glmerControl(optimizer="bobyqa", optCtrl=list(maxfun=100000)))

the code described in the above post works well and also gives model summary with

list_models

But does not show specific model elements when requested as below

list_models$fit1$coefficents

#Error in list_models$vasb3$coefficients : $ operator not defined for this S4 class

Question 1(unsolved): I would like to make it work so I can request any elements (e.g. coefficients) from list_models

Question 2 (solved): how should I store it so next time I do not need to re-run models and get the model fit elements from a file?

for question 2 solution:

#save objects to R
saveRDS(list_models, file="list_models.RData")
#read object
readRDS("list_models.RData")

Thanking you!

nnet
  • 39
  • 4

4 Answers4

0

Using $ to access the coefficients doesn't work in the first place; you need to use coef() to access them. Saving and loading the data works as usual; once you load the file, it is in the original variable name that it was saved from.

gm1 <- glmer(cbind(incidence, size - incidence) ~ period + (1 | herd),
              data = cbpp, family = binomial)
list.models = list(gm1)
coef(list.models[[1]])

setwd("pathname")
save(list=c("list.models"), file="model.RData")
load("model.RData")
coef(list.models[[1]])

Output

$herd
   (Intercept)   period2   period3   period4
1   -0.8087134 -0.991925 -1.128216 -1.579745
2   -1.6974362 -0.991925 -1.128216 -1.579745
3   -0.9924771 -0.991925 -1.128216 -1.579745
4   -1.3593198 -0.991925 -1.128216 -1.579745
...
Vons
  • 3,277
  • 2
  • 16
  • 19
  • this works! however, how do we request coef by name of model than specifying element no. [1]? e.g., coef(list.models[[gm1]]). This is important as there will be several models in the list and specifying them by name would be easier than element number e.g. [1] – nnet Jun 05 '23 at 14:40
  • As @rps1227 suggested, list_models$fit1@beta works same as fixef(list_models[[1]]) but without column headings. Mainly I was able to pull only beta, what other elements could we request this way? e.g., list_models$vasb3@residuals did not work. how to retrieve confidence interval than manually exponentiating coefficients? Thanks to both of you! – nnet Jun 05 '23 at 14:49
  • The goal is to use this stored list/rdata to retrieve everything about models so we don't have to go back and rerun individual models. – nnet Jun 05 '23 at 14:56
0

See ?coef

get_coef <- function(list_of_models, variable) {
  lapply(list_of_models, function(x) { coef(x)[["id"]][[variable]] }) 
}

use

get_coef(list_of_models, "x")

result

[[1]]
[1] -0.7613817 -1.2056582 -1.5324857

[[2]]
[1]  1.0493066 -0.3348213  1.9110284
jyr
  • 690
  • 6
  • 20
0

merMod objects (the results of glmer() and lmer()) are S4 objects so to access individual elements you have to use @ and not $ which is for S3 objects such as lists. list_models$fit1@beta will give you the fixed effects. But using coef() or fixef() as suggested by the other answers is better and safer.

rps1227
  • 472
  • 2
  • 5
0

We could use tidy from broom.mixed package:

Here is an example:

library(lme4)
library(dplyr)
library(broom.mixed)

# Fit the glmer model
fit1 <- glmer(y ~ scale(x) + (1 | ID), data = dat, family = binomial, control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 100000)))

tidy(fit1)

  effect   group term            estimate std.error statistic  p.value
 <chr>    <chr> <chr>              <dbl>     <dbl>     <dbl>    <dbl>
1 fixed    NA    (Intercept)        0.632     0.212      2.97  0.00295
2 fixed    NA    scale(x)          -0.295     0.215     -1.37  0.170  
3 ran_pars ID    sd__(Intercept)    0        NA         NA    NA     

# To get access we could use pull: 

tidy(fit1) %>% 
  pull(std.error) # estimate or p.value etc. 
[1] 0.2124558 0.2151282        NA
TarJae
  • 72,363
  • 6
  • 19
  • 66