Mixed Models with Adaptive Gaussian Quadrature

[This article was first published on iProgn: Interactive Prediction Tools based on Joint Models, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Overview

In this post, I would like to introduce my new R package GLMMadaptive for fitting mixed-effects models for non-Gaussian grouped/clustered outcomes using marginal maximum likelihood.

Admittedly, there is a number of packages available for fitting similar models, e.g., lme4, glmmsr, glmmTMB, glmmEP, and glmmML among others; more information on other available packages can also be found in GLMM-FAQ. GLMMadaptive differs from these packages in approximating the integrals over the random effects in the definition of the marginal log-likelihood function using an adaptive Gaussian quadrature rule, while allowing for multiple correlated random effects.

An Example: Mixed Effects Logistic Regression

We illustrate the use of the package in the simple case of a mixed effects logistic regression. We start by simulating some data:

set.seed(1234)
n <- 100 # number of subjects
K <- 8 # number of measurements per subject
t_max <- 15 # maximum follow-up time
# we constuct a data frame with the design:
# everyone has a baseline measurment, and then measurements at random follow-up times
DF <- data.frame(id = rep(seq_len(n), each = K),
time = c(replicate(n, c(0, sort(runif(K - 1, 0, t_max))))),
sex = rep(gl(2, n/2, labels = c("male", "female")), each = K))
# design matrices for the fixed and random effects
X <- model.matrix(~ sex * time, data = DF)
Z <- model.matrix(~ time, data = DF)
betas <- c(-2.13, -0.25, 0.24, -0.05) # fixed effects coefficients
D11 <- 0.48 # variance of random intercepts
D22 <- 0.1 # variance of random slopes
# we simulate random effects
b <- cbind(rnorm(n, sd = sqrt(D11)), rnorm(n, sd = sqrt(D22)))
# linear predictor
eta_y <- drop(X %*% betas + rowSums(Z * b[DF$id, ]))
# we simulate binary longitudinal data
DF$y <- rbinom(n * K, 1, plogis(eta_y))


We continue by fitting the mixed effects logistic regression for the longitudinal outcome y assuming random intercepts for the random-effects part. The primary model-fitting function in the package is the mixed_model(), and has four required arguments, namely,
  • fixed: a formula for the fixed effects,  
  • random: a formula for the random effects,  
  • family: a family object specifying the type of response variable, and  
  • data: a data frame containing the variables in the previously mentioned formulas.
 Assuming that the package has been loaded using library(“GLMMadaptive”), the call to fit the random intercepts logistic regression is:

fm1 <- mixed_model(fixed = y ~ sex * time, random = ~ 1 | id, data = DF,
family = binomial())
summary(fm1)
#>
#> Call:
#> mixed_model(fixed = y ~ sex * time, random = ~1 | id, data = DF,
#> family = binomial())
#>
#> Data Descriptives:
#> Number of Observations: 800
#> Number of Groups: 100
#>
#> Model:
#> family: binomial
#> link: logit
#>
#> Fit statistics:
#> log.Lik AIC BIC
#> -374.5197 759.0395 772.0653
#>
#> Random effects covariance matrix:
#> StdDev
#> (Intercept) 2.064578
#>
#> Fixed effects:
#> Value Std.Err z-value p-value
#> (Intercept) -2.6162 0.4447 -5.8837 < 1e-04
#> sexfemale -1.3234 0.6435 -2.0567 0.039714
#> time 0.2850 0.0395 7.2077 < 1e-04
#> sexfemale:time 0.0401 0.0545 0.7356 0.461982
#>
#> Integration:
#> method: adaptive Gauss-Hermite quadrature rule
#> quadrature points: 11
#>
#> Optimization:
#> method: hybrid EM and quasi-Newton
#> converged: TRUE


By default, 11 quadrature points are used, but this can be adjusted using the nAGQ control argument. We extend model fm1 by also including a random slopes term; however, we assume that the covariance between the random intercepts and random slopes is zero. This is achieved by using the `||` symbol in the specification of the random argument. We fit the new model and compare it with fm1 using the anova() function that performs a likelihood ratio test:

fm2 <- mixed_model(fixed = y ~ sex * time, random = ~ time || id, data = DF,
family = binomial())
anova(fm1, fm2)
#>
#> AIC BIC log.Lik LRT df p.value
#> fm1 759.04 772.07 -374.52
#> fm2 730.94 730.94 -359.47 30.1 1 <0.0001


 We further extend the model by estimating the covariance between the random intercepts and random slopes, and we use 15 quadrature points for the numerical integration:

fm3 <- mixed_model(fixed = y ~ sex * time, random = ~ time | id, data = DF,
family = binomial(), nAGQ = 15)
summary(fm3)
#>
#> Call:
#> mixed_model(fixed = y ~ sex * time, random = ~time | id, data = DF,
#> family = binomial(), nAGQ = 15)
#>
#> Data Descriptives:
#> Number of Observations: 800
#> Number of Groups: 100
#>
#> Model:
#> family: binomial
#> link: logit
#>
#> Fit statistics:
#> log.Lik AIC BIC
#> -358.8617 731.7235 749.9597
#>
#> Random effects covariance matrix:
#> StdDev Corr
#> (Intercept) 0.7729
#> time 0.2539 0.5181
#>
#> Fixed effects:
#> Value Std.Err z-value p-value
#> (Intercept) -2.0841 0.3218 -6.4761 < 1e-04
#> sexfemale -1.1759 0.4763 -2.4690 0.013549
#> time 0.2640 0.0569 4.6424 < 1e-04
#> sexfemale:time 0.0030 0.0790 0.0380 0.969696
#>
#> Integration:
#> method: adaptive Gauss-Hermite quadrature rule
#> quadrature points: 15
#>
#> Optimization:
#> method: hybrid EM and quasi-Newton
#> converged: TRUE


Capabilities and Further Reading

The package offers a wide range of methods for standard generic functions in R applicable to regression models objects and mixed models objects in particular (e.g., fixef(), ranef(), etc.); more information can be found in the Methods for MixMod Objects vignette. In addition, some highlights of its capabilities:
  • It allows for user-defined family objects implementing not standardly available outcome distributions; more information can be found in the Custom Models vignette.
  • It can calculate fixed effects coefficients with a marginal / population averaged interpretation using the function marginal_coefs().
  • Function effectPlotData() calculates predictions and confidence intervals for constructing effect plots.
The development version of the package is available on the GitHub repository.

To leave a comment for the author, please follow the link and comment on their blog: iProgn: Interactive Prediction Tools based on Joint Models.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)