An R package for GLMM estimation with two random effects?

Question

I am looking for an R package to make an estimate on a general linear model with two random effects. I was used to lme4 (functionglmer), but from what I understood, adaptive Gaussian quadrature is only available for one random effects model.

Do you have another suggestion for me?

Thanks !

PS: I looked on the page https://rpubs.com/kaz_yos/glmm1 which contains many suggestions, but this page is almost 3 years old and I was wondering if something new and not based on Bayesian (time calculation too important) is available.

Answer 1

You could try mixed_model() from the GLMMadaptive package, which uses adaptive Gaussian quadrature, and which fits models that lme4::glmer() doesn't.

Taking the data from the website you linked to:

library(tidyr)
library(dplyr)
library(GLMMadaptive)
library(lme4)

epi <- haven::read_sas("http://www.hsph.harvard.edu/fitzmaur/ala2e/epilepsy.sas7bdat")
names(epi) <- tolower(names(epi))

epiL <- gather_(epi, key_col = "time", value_col = "y", gather_cols = paste0("y", 0:4)) %>%
  mutate(time = as.numeric(gsub("y", "", .$time))) %>%
  arrange(id, time)

epiL$T[epiL$time == 0] <- 8
#> Warning: Unknown or uninitialised column: 'T'.
epiL$T[epiL$time != 0] <- 2
epiL$logT <- log(epiL$T)

m1 <- glmer(
  formula = y ~ time * trt + (1 + time | id),
  data = epiL,
  family = poisson(link = "log"),
  nAGQ = 7
)
#> Error in updateGlmerDevfun(devfun, glmod$reTrms, nAGQ = nAGQ): nAGQ > 1 is only available for models with a single, scalar random-effects term

m2 <- mixed_model(
  fixed = y ~ time * trt,
  random = ~ 1 + time | id,
  data = epiL,
  family = poisson(link = "log")
)

summary(m2)
#> 
#> Call:
#> mixed_model(fixed = y ~ time * trt, random = ~1 + time | id, 
#>     data = epiL, family = poisson(link = "log"))
#> 
#> Data Descriptives:
#> Number of Observations: 295
#> Number of Groups: 59 
#> 
#> Model:
#>  family: poisson
#>  link: log 
#> 
#> Fit statistics:
#>    log.Lik      AIC      BIC
#>  -1142.416 2298.831 2313.374
#> 
#> Random effects covariance matrix:
#>              StdDev    Corr
#> (Intercept)  0.7217        
#> time         0.2074  0.0983
#> 
#> Fixed effects:
#>             Estimate Std.Err z-value  p-value
#> (Intercept)   2.9140  0.1428 20.4115  < 1e-04
#> time         -0.4301  0.0460 -9.3513  < 1e-04
#> trt           0.0590  0.1964  0.3005 0.763801
#> time:trt     -0.1377  0.0643 -2.1417 0.032216
#> 
#> Integration:
#> method: adaptive Gauss-Hermite quadrature rule
#> quadrature points: 11
#> 
#> Optimization:
#> method: hybrid EM and quasi-Newton
#> converged: TRUE

^{Created on 2019-04-15 by the reprex package (v0.2.1)}