1
$\begingroup$

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.

$\endgroup$
3
$\begingroup$

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)

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ It works perfectly! Thank you very much for the suggestion. $\endgroup$ – Flora Grappelli Apr 16 '19 at 12:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.