knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) library("GLMMadaptive")

Package **GLMMadaptive** provides a suit of functions for fitting and post-processing
mixed effects models for grouped/clustered outcomes which have a distribution other than a
normal distribution. In particular, let $y_i$ denote a vector of grouped/clustered outcome for
the $i$-th sample unit ($i = 1, \ldots, n$). The conditional distribution of $y_i$ given a
vector of random effects $b_i$ is assumed to be a member of the extended exponential
family, with linear predictor given by $$g{E(y_i \mid b_i)} = X_i \beta + Z_i b_i,$$
where $g(\cdot)$ denotes a monotonic link function, $X_i$ a design matrix for the fixed
effects coefficients $\beta$, and $Z_i$ a design matrix for the random effects
coefficients $b_i$. Typically, matrix $Z_i$ is assumed to be a subset of $X_i$. The random
effects are assumed to follow a normal distribution with mean 0 and variance-covariance
matrix $D$. In addition, the distribution $[y_i \mid b_i]$ may potentially have
extra dispersion/shape parameters $\phi$.

The package focuses in settings in which the distribution $[y_i \mid b_i]$ is not normal
and/or the link function $g(\cdot)$ is not the identity. In these settings, the estimation
of the model is complicated by the fact that the marginal log-likelihood function of the
observed $y_i$ cannot be derived analytically. In particular, the log-likelihood function
has the form:
$$\begin{eqnarray*}
\ell(\theta) & = & \sum_{i = 1}^n \log p(y_i; \theta)\
& = & \sum_{i = 1}^n \log \int p(y_i \mid b_i; \theta) \, p(b_i; \theta) \, db_i,
\end{eqnarray*}$$
where $\theta$ denotes the full parameter vector including the fixed effects, the extra
potential dispersion/shape parameters $\phi$ and the unique element of the covariance
matrix $D$, and $p(\cdot)$ denotes a probability density or probability mass function. The
integral in the definition of $\ell(\theta)$ does not have a closed-form solution, and
numerical approximations are required to obtain the maximum likelihood estimates.

In the literature several approaches have been proposed to approximate such integrals, and
a nice overview is given in
Pinheiro and Chao (2006). A typical approach to
approximate these integrals is the Laplace approximation. However, the general consensus
has been that in the standard but difficult cases of binary/dichotomous data and count data
with small counts and few repeated measurements, the accuracy of this approximation is
rather low. Due to this fact, the general consensus is that the gold standard numerical
approximation method is the adaptive Gaussian quadrature rule (*note: we focus here on
maximum likelihood estimation; under the Bayesian paradigm, approaches, such as, MCMC and
Hamiltonian Monte Carlo also provide accurate evaluation of the integrals*). This is more
computationally intensive but also more accurate. This package provides an efficient
implementation of the adaptive Gaussian quadrature rule, allowing for multiple correlated
random effects (e.g., random intercepts, linear and quadratic random slopes) but currently
a single grouping factor (i.e., no nested or crossed random effects designs).

A hybrid optimization procedure is implemented starting with an EM algorithm, treating the random effects as 'missing data', followed by a direct optimization procedure with a quasi-Newton algorithm. Fine control of this procedure is allowed with a series of control arguments.

We illustrate the use of the package in the standard case of a mixed effects logistic regression. That is, the distribution of $[y_i \mid b_i]$ is binomial, and the distribution of $[b_i]$ multivariate normal.

We start by simulating some data for a binary longitudinal outcome:

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 fit a mixed effects logistic regression for `y`

, assuming random
intercepts for the random-effects part. The basic model-fitting function in
**GLMMadaptive** is called `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. Hence, the call to fit the
random intercepts logistic regression is:

fm1 <- mixed_model(fixed = y ~ sex * time, random = ~ 1 | id, data = DF, family = binomial())

The summary method gives a detailed output of the model:

```
summary(fm1)
```

We continue by checking the impact of the chosen number of quadrature points to the parameters estimates and the log-likelihood value at convergence. First, we refit the model with an increasing number of quadrature points. The default when the number of random effects is smaller or equal to two is 11 points. We fit then with 15, and 21 points:

fm1_q11 <- fm1 fm1_q15 <- update(fm1_q11, nAGQ = 15) fm1_q21 <- update(fm1_q11, nAGQ = 21) models <- list("nAGQ=11" = fm1_q11, "nAGQ=15" = fm1_q15, "nAGQ=21" = fm1_q21)

We now extract from the model the estimated parameter for the fixed effects (using
function `fixef()`

), for the random effects, and the log-likelihood (using function
`logLik()`

):

extract <- function (obj) { c(fixef(obj), "var_(Intercept)" = obj$D[1, 1], "logLik" = logLik(obj)) } sapply(models, extract)

We observe a rather stable model with virtually no differences between the different choices of quadrature points.

We first compare the model with a simple logistic regression that does not include any random effects:

km <- glm(y ~ sex * time, data = DF, family = binomial()) anova(fm1, km)

We obtain a highly significant p-value suggesting that there are correlations in the data
that cannot be ignored. *Note:* the `anova()`

method that performs the likelihood ratio
test calculates the p-value using the standard $\chi^2$ distribution, here with one
degree of freedom. However, because the null hypothesis for testing variance parameters is
on the boundary of the corresponding parameter space, it would be more appropriate to use
a mixture of $\chi^2$ distributions.

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, i.e.,

fm2 <- mixed_model(fixed = y ~ sex * time, random = ~ time || id, data = DF, family = binomial())

The likelihood ratio test between the two models is computed with function `anova()`

. When
two `"MixMod"`

objects are provided, the function assumes that the first object represents
the model under the null hypothesis, and the second object the model under the
alternative, i.e.,

```
anova(fm1, fm2)
```

The results suggest that we need the random slopes term. We continue by testing whether the covariance between the random effects terms is zero. The model under the alternative hypothesis is (results not shown):

fm3 <- mixed_model(fixed = y ~ sex * time, random = ~ time | id, data = DF, family = binomial())

And again the likelihood ratio test is performed by a call to `anova()`

(results not shown):

```
anova(fm2, fm3)
```

The results now suggest that indeed the covariance between the two random effects terms is not statistically different from zero.

We continue our illustration with a Poisson mixed effects model. We start again by simulating some data for a Poisson longitudinal outcome:

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) betas <- c(2.13, -0.25, 0.24, -0.05) # fixed effects coefficients D11 <- 0.48 # variance of random intercepts # we simulate random effects b <- rnorm(n, sd = sqrt(D11)) # linear predictor eta_y <- drop(X %*% betas + b[DF$id]) # we simulate Poisson longitudinal data DF$y <- rpois(n * K, exp(eta_y))

We fit the mixed effects Poisson regression for `y`

assuming random intercepts for the
random-effects part.

gm1 <- mixed_model(fixed = y ~ sex * time, random = ~ 1 | id, data = DF, family = poisson())

As previously, the summary method gives a detailed output of the model:

summary(gm1) #> #> Call: #> mixed_model(fixed = y ~ sex * time, random = ~1 | id, data = DF, #> family = poisson()) #> #> Data Descriptives: #> Number of Observations: 800 #> Number of Groups: 100 #> #> Model: #> family: poisson #> link: log #> #> Fit statistics: #> log.Lik AIC BIC #> -2782.331 5574.663 5587.689 #> #> Random effects covariance matrix: #> StdDev #> (Intercept) 0.9785518 #> #> Fixed effects: #> Estimate Std.Err z-value p-value #> (Intercept) 2.9461 0.3091 9.5316 < 1e-04 #> sexfemale -0.8651 0.2711 -3.1913 0.0014162 #> time 0.2404 0.0016 151.3470 < 1e-04 #> sexfemale:time -0.0511 0.0028 -18.3994 < 1e-04 #> #> Integration: #> method: adaptive Gauss-Hermite quadrature rule #> quadrature points: 11 #> #> Optimization: #> method: hybrid EM and quasi-Newton #> converged: TRUE

In settings with few subjects, repeated measurements of separation effects, the package
also allows to place a penalty in the fixed effects regression coefficients $\beta$. The
penalty/prior is in the form of a Student's t distribution with mean 0, scale parameter 1,
and 3 degrees of freedom, and it is placed in all $\beta$ coefficients except from the
intercept. The penalized model can be fitted by setting argument `penalized`

to `TRUE`

,
i.e.,

gm2 <- mixed_model(fixed = y ~ sex * time, random = ~ 1 | id, data = DF, family = poisson(), penalized = TRUE)

cbind('unpenalized' = fixef(gm1), 'penalized' = fixef(gm2)) #> unpenalized penalized #> (Intercept) 2.94614349 2.94238986 #> sexfemale -0.86510298 -0.82073289 #> time 0.24039622 0.24043929 #> sexfemale:time -0.05106264 -0.05118503

In this example we observe small differences between the penalized and unpenalized models. The users have the option to alter the specification of the Student's t penalty by directly specifying the mean, scale and degrees of freedom arguments of the distribution. For example, a ridge penalty could be placed by setting the degrees of freedom to a high value. The call in this case should be:

gm3 <- mixed_model(fixed = y ~ sex * time, random = ~ 1 | id, data = DF, family = poisson(), penalized = list(pen_mu = 0, pen_sigma = 1, pen_df = 200))

**Any scripts or data that you put into this service are public.**

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.