This tutorial provides an overview of the mathematical framework behind GLMMs, the likelihood function, numerical integration using AGQ, parameter estimation, and prediction. We demonstrate how to fit GLMMs using the {GLMMadaptive} package in R with examples of binary, Poisson, and negative binomial response data.
Generalized Linear Mixed Models (GLMMs) are a powerful statistical framework for analyzing complex data structures that involve hierarchical or clustered data with non-normal outcomes. GLMMs extend the concept of Generalized Linear Models (GLMs) by incorporating random effects to account for the correlation among observations within the same group. The Generalized Linear Mixed Models using Adaptive Gaussian Quadrature with {GLMMadaptive} package in R provides a flexible and efficient tool for fitting GLMMs using Adaptive Gaussian Quadrature (AGQ) for numerical integration using the adaptive Gauss-Hermite quadrature rule. This involves modeling fixed effects, random effects, and the response distribution using likelihood-based methods. The package supports various response distributions, complex random effects structures, and custom link functions. Multiple random effects terms can be included for the grouping factor (e.g., random intercepts, random linear slopes, random quadratic slopes), but currently only a single grouping factor is allowed.
The key features of the {GLMMadaptive} package include:
Supports various response distributions, such as:
tools for making predictions with confidence intervals.
Supports marginal predictions (accounting for random effects).
A GLMM can be expressed as:
\[ g(\mu_{ij}) = X{ij}\beta + Z_{ij}u_i \]
The conditional distribution of the response \(y_{ij}\), given the random effects \(u_i\), is part of the exponential family:
\[ f(y_{ij} \| u_i) = \exp\left(\frac{y_{ij}\theta_{ij} - b(\theta_{ij})}{\phi} + c(y_{ij}, \phi)\right), \] where:
The likelihood function involves integrating out the random effects \(u_i\):
\[ L(\beta, \theta) = \prod_{i=1}^m \int \left[ \prod_{j=1}^{n_i} f(y_{ij} | u_i) \right] f(u_i) , du_i \]
where: .
AGQ is a numerical integration technique used to approximate these integrals efficiently and accurately. It is particularly suitable for GLMMs because it:
Exploits the assumption that random effects are normally distributed.
the quadrature points to the mode of the integrand for better accuracy.
To approximate the likelihood, AGQ is used. The key steps are:
The integral is centered around the mode of the conditional random-effects density:
\[ f(u_i \| y_i) \propto f(y_i \| u_i) f(u_i), \] where:
\(y_i\) = \(y_{i1}, y_{i2}, \dots, y_{in_i})\): Observations for group \(i\).
The mode \(\hat{u}_i\) is found by maximizing \(\log f(u_i \| y_i)\).
The integral is approximated using Gaussian quadrature. Let \(\tilde{u}_i = \hat{u}_i + s z_k\), where \(z_k\) are quadrature points and \(s\) is a scaling factor based on the curvature (Hessian) at \(\hat{u}_i\). The integral becomes:
\[ \int f(y_i \| u_i) f(u_i) , du_i \approx \sum_{k=1}^Q w_k f(y_i \| \tilde{u}\_i) f(\tilde{u}_i), \]
where \(Q\) is the number of quadrature points, \(w_k\) are quadrature weights, and \(\tilde{u}_i\) are adjusted points.
The log-likelihood is computed as:
\[ \log L(\beta, \theta) = \sum_{i=1}^m \log \left( \sum{k=1}^Q w_k f(y_i \| \tilde{u}_i) f(\tilde{u}_i) \right). \]
Parameters \(\beta\) (fixed effects) and \(\theta\) (random effects and dispersion) are estimated by maximizing the log-likelihood:
\[ \hat{\beta}, \hat{\theta} = \arg \max_{\beta, \theta} \log L(\beta, \theta). \]
Optimization methods like Newton-Raphson or quasi-Newton are used.
After fitting the model, the posterior distribution of the random effects \(u_i\) can be estimated:
\[ f(u_i \| y_i) \propto f(y_i \| u_i) f(u_i). \] Predictions for random effects are often based on:
Predictions can be made for:
Marginal Means (without random effects): \(\hat{y}_{ij} = g^{-1}(X_{ij}\hat{\beta})\).
Subject-Specific Means (with random effects): \(\hat{y}_{ij} = g^{-1}(X_{ij}\hat{\beta} + Z_{ij}\hat{u}_i)\)
In this tutorial, we will demonstrate how to fit Generalized Linear Mixed Models (GLMMs) using Adaptive Gaussian Quadrature in R with the {GLMMadaptive} package. Basic features of the package include:
The package contains a single model-fitting function named mixed_model() with four required arguments, 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.
Methods for standard generics are provided, i.e., coef(), fixef(), ranef(), vcov(), logLik(), summary(), anova(), confint(), fitted(), residuals(), predict(), and simulate().
Negative binomial mixed models can be fitted using the negative.binomial() family object.
Zero-inflated Poisson and negative binomial models using the zi.poisson() and zi.negative.binomial() family objects.
Hurdle Poisson and negative binomial models using the hurdle.poisson() and hurdle.negative.binomial() family objects.
Zero-inflated binomial models using the zi.binomial() family objects.
Two-part/hurdle mixed models for semi-continuous normal data using the hurdle.lognormal() family object.
Mixed models for censored normal data using the censored.normal() family object.
Continuation ratio mixed models for ordinal data using functions cr_setup() and cr_marg_probs().
Beta and hurdle Beta mixed effects models using beta.fam() and hurdle.beta.fam() family objects.
Gamma mixed effects models using the Gamma() or Gamma.fam() family object.
Linear mixed effects models with right and left censored data using the censored.normal() family object.
Users may also specify their own log-density function for the repeated measurements response variable, and the internal algorithms will take care of the optimization.
Calculates the marginalized coefficients using the idea of Hedeker et al. (2017) using function marginal_coefs().
Predictions with confidence interval for constructing effects plots are provided by function effectPlotData().
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. nAGQ is an optional argument that specifies the number of quadrature points to use in the adaptive Gaussian quadrature approximation. The default is 11 when the number of random effects is one or two, and 7 otherwise. We will walk through the process of specifying the model, fitting the model, and interpreting the results.
Following R packages are required to run this notebook. If any of these packages are not installed, you can install them using the code below:
packages <- c('tidyverse',
'dlookr',
'sjPlot',
'jtools',
'lme4',
'glmmTMB',
'MASS',
'MuMIn',
'performance',
'GLMMadaptive',
'mlmRev'
)#| warning: false
#| error: false
# Install missing packages
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)
# Verify installation
cat("Installed packages:\n")
print(sapply(packages, requireNamespace, quietly = TRUE))# Load packages with suppressed messages
invisible(lapply(packages, function(pkg) {
suppressPackageStartupMessages(library(pkg, character.only = TRUE))
}))# Check loaded packages
cat("Successfully loaded packages:\n")Successfully loaded packages:print(search()[grepl("package:", search())]) [1] "package:mlmRev" "package:GLMMadaptive" "package:performance"
[4] "package:MuMIn" "package:MASS" "package:glmmTMB"
[7] "package:lme4" "package:Matrix" "package:jtools"
[10] "package:sjPlot" "package:dlookr" "package:lubridate"
[13] "package:forcats" "package:stringr" "package:dplyr"
[16] "package:purrr" "package:readr" "package:tidyr"
[19] "package:tibble" "package:ggplot2" "package:tidyverse"
[22] "package:stats" "package:graphics" "package:grDevices"
[25] "package:utils" "package:datasets" "package:methods"
[28] "package:base" We will use the Contraception dataset from the mlmRev package, which come from the 1988 Bangladesh Fertility Survey and available for download from here.
# load the data as tibble
mf.binomial<-as_tibble(Contraception, package = "mlmRev")
str(mf.binomial)tibble [1,934 × 6] (S3: tbl_df/tbl/data.frame)
$ woman : Factor w/ 1934 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
$ district: Factor w/ 60 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
$ use : Factor w/ 2 levels "N","Y": 1 1 1 1 1 1 1 1 1 1 ...
$ livch : Factor w/ 4 levels "0","1","2","3+": 4 1 3 4 1 1 4 4 2 4 ...
$ age : num [1:1934] 18.44 -5.56 1.44 8.44 -13.56 ...
$ urban : Factor w/ 2 levels "N","Y": 2 2 2 2 2 2 2 2 2 2 ...We fit a mixed effects logistic regression for use and age, urban, livch are fixed effect and district is random effect (random intercepts) random-effects part.
# covert the livch variable to factor
mf.binomial$livch<-as.factor(mf.binomial$livch)
# fit the mixed effects logistic regression model
model.logit_q11 <- mixed_model(fixed = use ~ age+urban+livch, random = ~1|district,
family = binomial(),
control = glmerControl(optimizer = "nAGQ = 11"),
data = mf.binomial)Warning in mixed_model(fixed = use ~ age + urban + livch, random = ~1 | :
unknown names in control: restart_edge, boundary.tol, calc.derivs,
use.last.params, checkControl, checkConv, optCtrl, tolPwrss, compDev,
nAGQ0initStepsummary(model.logit_q11)
Call:
mixed_model(fixed = use ~ age + urban + livch, random = ~1 |
district, data = mf.binomial, family = binomial(), control = glmerControl(optimizer = "nAGQ = 11"))
Data Descriptives:
Number of Observations: 1934
Number of Groups: 60
Model:
family: binomial
link: logit
Fit statistics:
log.Lik AIC BIC
-1206.674 2427.348 2442.009
Random effects covariance matrix:
StdDev
(Intercept) 0.4645476
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -1.6904 0.1478 -11.4406 < 1e-04
age -0.0266 0.0079 -3.3727 0.00074429
urbanY 0.7324 0.1195 6.1290 < 1e-04
livch1 1.1094 0.1580 7.0206 < 1e-04
livch2 1.3766 0.1748 7.8750 < 1e-04
livch3+ 1.3457 0.1796 7.4923 < 1e-04
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11
Optimization:
method: EM
converged: TRUE 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 (nAGQ). We fit then with 15, and 21 points:
model.logit_q15 <- stats::update(model.logit_q11, nAGQ = 15)
model.logit_q21 <- stats::update(model.logit_q11, nAGQ = 21)
models.logit_01 <- list("nAGQ=11" = model.logit_q11, "nAGQ=15" = model.logit_q15, "nAGQ=21" = model.logit_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.logit_01, extract) nAGQ=11 nAGQ=15 nAGQ=21
(Intercept) -1.690396e+00 -1.690390e+00 -1.690390e+00
age -2.660182e-02 -2.660177e-02 -2.660177e-02
urbanY 7.323517e-01 7.323533e-01 7.323533e-01
livch1 1.109383e+00 1.109382e+00 1.109382e+00
livch2 1.376598e+00 1.376596e+00 1.376596e+00
livch3+ 1.345699e+00 1.345696e+00 1.345696e+00
var_(Intercept) 2.158045e-01 2.157964e-01 2.157964e-01
logLik -1.206674e+03 -1.206674e+03 -1.206674e+03We observe a rather stable model with virtually no differences between the different choices of quadrature points.
In this section we use the Salamanders data set from the {glmmTMB} package. The data set containing counts of salamanders with site covariates and sampling covariates. Each of 23 sites was sampled 4 times. (Price et al. (2016); Price et al. 2015).
# load the data as tibble
mf.poisson<-as_tibble(Salamanders, package = "glmmTMB")
str(mf.poisson)tibble [644 × 9] (S3: tbl_df/tbl/data.frame)
$ site : Ord.factor w/ 23 levels "R-1"<"R-2"<"R-3"<..: 13 14 15 1 2 3 4 5 6 7 ...
$ mined : Factor w/ 2 levels "yes","no": 1 1 1 2 2 2 2 2 2 2 ...
$ cover : num [1:644] -1.442 0.298 0.398 -0.448 0.597 ...
$ sample: int [1:644] 1 1 1 1 1 1 1 1 1 1 ...
$ DOP : num [1:644] -0.596 -0.596 -1.191 0 0.596 ...
$ Wtemp : num [1:644] -1.2294 0.0848 1.0142 -3.0234 -0.1443 ...
$ DOY : num [1:644] -1.497 -1.497 -1.294 -2.712 -0.687 ...
$ spp : Factor w/ 7 levels "GP","PR","DM",..: 1 1 1 1 1 1 1 1 1 1 ...
$ count : int [1:644] 0 0 0 2 2 1 1 2 4 1 ...We fit a mixed effects Poisson regression for count and mined, spp, cover are fixed effect and site is random effect (random intercepts) random-effects part.
# fit the mixed effects Poisson regression model
model.poisson_01 <- mixed_model(fixed = count ~ mined + spp + cover, random = ~1|site,
family = poisson(),
data = mf.poisson)
summary(model.poisson_01)
Call:
mixed_model(fixed = count ~ mined + spp + cover, random = ~1 |
site, data = mf.poisson, family = poisson())
Data Descriptives:
Number of Observations: 644
Number of Groups: 23
Model:
family: poisson
link: log
Fit statistics:
log.Lik AIC BIC
-972.1232 1964.246 1975.601
Random effects covariance matrix:
StdDev
(Intercept) 0.5651011
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -1.6904 0.2580 -6.5513 < 1e-04
minedno 2.3960 0.3350 7.1531 < 1e-04
sppPR -1.3864 0.2152 -6.4432 < 1e-04
sppDM 0.2305 0.1289 1.7885 0.073702
sppEC-A -0.7701 0.1711 -4.5022 < 1e-04
sppEC-L 0.6211 0.1193 5.2060 < 1e-04
sppDES-L 0.6791 0.1181 5.7491 < 1e-04
sppDF 0.0800 0.1334 0.5997 0.548733
cover -0.1206 0.1649 -0.7313 0.464609
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE We extend model 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.,
# fit the mixed effects Poisson regression model
model.poisson_02 <- mixed_model(fixed = count ~ mined + spp + cover, random = ~ mined||site,
family = poisson(),
data = mf.poisson)
summary(model.poisson_02)
Call:
mixed_model(fixed = count ~ mined + spp + cover, random = ~mined ||
site, data = mf.poisson, family = poisson())
Data Descriptives:
Number of Observations: 644
Number of Groups: 23
Model:
family: poisson
link: log
Fit statistics:
log.Lik AIC BIC
-972.2275 1966.455 1978.946
Random effects covariance matrix:
StdDev
(Intercept) 0.5714
minedno 0.0874
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -1.6924 0.2604 -6.4998 < 1e-04
minedno 2.3964 0.3402 7.0449 < 1e-04
sppPR -1.3863 0.2152 -6.4429 < 1e-04
sppDM 0.2305 0.1289 1.7885 0.07370
sppEC-A -0.7701 0.1711 -4.5023 < 1e-04
sppEC-L 0.6212 0.1193 5.2065 < 1e-04
sppDES-L 0.6792 0.1181 5.7495 < 1e-04
sppDF 0.0800 0.1334 0.5998 0.54865
cover -0.1181 0.1679 -0.7031 0.48200
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE # compare the two models
anova(model.poisson_01, model.poisson_02)
AIC BIC log.Lik LRT df p.value
model.poisson_01 1964.25 1975.60 -972.12
model.poisson_02 1966.46 1978.95 -972.23 0.21 1 0.6478The results now suggest that indeed the covariance between the two random effects terms is not statistically different from zero.
{GLMMapative} 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.,
# fit the mixed effects Poisson regression model with penalized fixed effects
model.poisson_03 <- mixed_model(fixed = count ~ mined + spp + cover, random = ~ 1|site,
family = poisson(), penalized = TRUE,
data = mf.poisson)
summary(model.poisson_03)
Call:
mixed_model(fixed = count ~ mined + spp + cover, random = ~1 |
site, data = mf.poisson, family = poisson(), penalized = TRUE)
Data Descriptives:
Number of Observations: 644
Number of Groups: 23
Model:
family: poisson
link: log
Fit statistics:
log.Lik AIC BIC
-981.3343 1982.669 1994.024
Random effects covariance matrix:
StdDev
(Intercept) 0.5570291
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -1.4874 0.2432 -6.1158 < 1e-04
minedno 2.1016 0.3250 6.4670 < 1e-04
sppPR -1.3490 0.2065 -6.5319 < 1e-04
sppDM 0.2209 0.1260 1.7526 0.079678
sppEC-A -0.7641 0.1666 -4.5878 < 1e-04
sppEC-L 0.6102 0.1165 5.2371 < 1e-04
sppDES-L 0.6681 0.1153 5.7923 < 1e-04
sppDF 0.0713 0.1305 0.5467 0.584606
cover -0.0388 0.1637 -0.2372 0.812480
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE cbind('unpenalized' = fixef(model.poisson_01), 'penalized' = fixef(model.poisson_03)) unpenalized penalized
(Intercept) -1.69039072 -1.48738250
minedno 2.39603247 2.10164569
sppPR -1.38636329 -1.34901230
sppDM 0.23050847 0.22087622
sppEC-A -0.77010542 -0.76411050
sppEC-L 0.62111277 0.61024378
sppDES-L 0.67912236 0.66812308
sppDF 0.08001793 0.07134521
cover -0.12056806 -0.03882774In 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:
# fit the mixed effects Poisson regression model with penalized fixed effects
model.poisson_04 <- mixed_model(fixed = count ~ mined + spp + cover, random = ~ 1|site,
family = poisson(),
penalized = list(pen_mu = 0, pen_sigma = 1, pen_df = 200),
data = mf.poisson)
summary(model.poisson_04)
Call:
mixed_model(fixed = count ~ mined + spp + cover, random = ~1 |
site, data = mf.poisson, family = poisson(), penalized = list(pen_mu = 0,
pen_sigma = 1, pen_df = 200))
Data Descriptives:
Number of Observations: 644
Number of Groups: 23
Model:
family: poisson
link: log
Fit statistics:
log.Lik AIC BIC
-977.6296 1975.259 1986.614
Random effects covariance matrix:
StdDev
(Intercept) 0.5576869
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -1.4751 0.2310 -6.3845 < 1e-04
minedno 2.0835 0.3060 6.8093 < 1e-04
sppPR -1.3464 0.2053 -6.5571 < 1e-04
sppDM 0.2203 0.1258 1.7514 0.079885
sppEC-A -0.7636 0.1662 -4.5940 < 1e-04
sppEC-L 0.6096 0.1163 5.2428 < 1e-04
sppDES-L 0.6675 0.1151 5.7996 < 1e-04
sppDF 0.0708 0.1303 0.5438 0.586603
cover -0.0336 0.1609 -0.2091 0.834373
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE Negative binomial mixed models can be fitted using the negative.binomial() family object.
# fit the mixed effects Poisson regression model
model.nb <- mixed_model(fixed = count ~ mined + spp + cover, random = ~1|site,
family = negative.binomial(),
data = mf.poisson)
summary(model.nb)
Call:
mixed_model(fixed = count ~ mined + spp + cover, random = ~1 |
site, data = mf.poisson, family = negative.binomial())
Data Descriptives:
Number of Observations: 644
Number of Groups: 23
Model:
family: negative binomial
link: log
Fit statistics:
log.Lik AIC BIC
-825.9648 1673.93 1686.42
Random effects covariance matrix:
StdDev
(Intercept) 0.5356483
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -1.7395 0.2899 -6.0010 < 1e-04
minedno 2.3789 0.3408 6.9810 < 1e-04
sppPR -1.3237 0.2877 -4.6019 < 1e-04
sppDM 0.3694 0.2234 1.6532 0.09829453
sppEC-A -0.7107 0.2530 -2.8090 0.00496941
sppEC-L 0.5655 0.2192 2.5799 0.00988368
sppDES-L 0.7905 0.2166 3.6495 0.00026275
sppDF 0.3083 0.2329 1.3236 0.18562086
cover -0.1046 0.1677 -0.6236 0.53292062
log(dispersion) parameter:
Estimate Std.Err
-0.0592 0.135
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE Zero-inflated Poisson models can be fitted using the zi.poisson() family object.
# fit the mixed effects zero-inflated Poisson regression model
model.zi <- mixed_model(fixed = count ~ mined + spp + cover, random = ~1|site,
family = zi.poisson(),
zi_fixed = ~ mined + spp + cover,
data = mf.poisson)
summary(model.zi)
Call:
mixed_model(fixed = count ~ mined + spp + cover, random = ~1 |
site, data = mf.poisson, family = zi.poisson(), zi_fixed = ~mined +
spp + cover)
Data Descriptives:
Number of Observations: 644
Number of Groups: 23
Model:
family: zero-inflated poisson
link: log
Fit statistics:
log.Lik AIC BIC
-874.0457 1786.091 1807.666
Random effects covariance matrix:
StdDev
(Intercept) 0.1890082
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -0.2181 0.2413 -0.9038 0.366083
minedno 1.2575 0.2195 5.7297 < 1e-04
sppPR -0.5753 0.2901 -1.9829 0.047375
sppDM 0.2410 0.1418 1.6987 0.089378
sppEC-A -0.2023 0.2010 -1.0063 0.314253
sppEC-L 0.6389 0.1313 4.8675 < 1e-04
sppDES-L 0.6028 0.1284 4.6944 < 1e-04
sppDF 0.0502 0.1577 0.3181 0.750436
cover -0.1827 0.0789 -2.3172 0.020495
Zero-part coefficients:
Estimate Std.Err z-value p-value
(Intercept) 1.1298 0.3926 2.8779 0.0040032
minedno -2.0505 0.3174 -6.4605 < 1e-04
sppPR 1.5002 0.5053 2.9688 0.0029899
sppDM -0.3211 0.4144 -0.7748 0.4384351
sppEC-A 1.1196 0.4358 2.5689 0.0102030
sppEC-L -0.1397 0.3969 -0.3521 0.7247834
sppDES-L -0.4347 0.4042 -1.0754 0.2822063
sppDF -0.5028 0.4566 -1.1012 0.2707919
cover -0.0215 0.1437 -0.1499 0.8808374
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE Zero-inflated negative binomial models can be fitted using the zi.negative.binomial() family object.
# fit the mixed effects Poisson regression model
model.zi.nb <- mixed_model(fixed = count ~ mined + spp + cover, random = ~1|site,
family = zi.negative.binomial(),
zi_fixed = ~ mined + spp + cover,
data = mf.poisson)
summary(model.zi.nb)
Call:
mixed_model(fixed = count ~ mined + spp + cover, random = ~1 |
site, data = mf.poisson, family = zi.negative.binomial(),
zi_fixed = ~mined + spp + cover)
Data Descriptives:
Number of Observations: 644
Number of Groups: 23
Model:
family: zero-inflated negative binomial
link: log
Fit statistics:
log.Lik AIC BIC
-813.0351 1666.07 1688.78
Random effects covariance matrix:
StdDev
(Intercept) 0.3538294
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -0.6673 0.3824 -1.7451 0.0809676
minedno 1.5897 0.3768 4.2193 < 1e-04
sppPR -1.5691 0.3001 -5.2283 < 1e-04
sppDM 0.1465 0.2354 0.6224 0.5336926
sppEC-A -0.8771 0.2667 -3.2884 0.0010075
sppEC-L 0.3428 0.2355 1.4558 0.1454517
sppDES-L 0.5422 0.2287 2.3706 0.0177578
sppDF -0.0893 0.2439 -0.3659 0.7144036
cover -0.2303 0.1372 -1.6791 0.0931314
Zero-part coefficients:
Estimate Std.Err z-value p-value
(Intercept) 2.0512 1.0353 1.9814 0.0475472
minedno -6.2910 3.3492 -1.8783 0.0603345
sppPR -2.8438 1.6721 -1.7008 0.0889886
sppDM -2.5648 1.2060 -2.1267 0.0334453
sppEC-A -1.9382 1.4223 -1.3627 0.1729854
sppEC-L -2.6403 1.3086 -2.0176 0.0436297
sppDES-L -2.6414 1.1739 -2.2502 0.0244362
sppDF -3.7705 1.4058 -2.6822 0.0073135
cover -0.8901 0.5714 -1.5578 0.1192817
log(dispersion) parameter:
Estimate Std.Err
0.0973 0.1443
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE # compare the two models
anova(model.zi, model.zi.nb)
AIC BIC log.Lik LRT df p.value
model.zi 1786.09 1807.67 -874.05
model.zi.nb 1666.07 1688.78 -813.04 122.02 1 <0.0001Hurdle Poisson models can be fitted using the hurdle.poisson() family object.
# fit the mixed effects Poisson regression model
model.hardle <- mixed_model(fixed = count ~ mined + spp + cover, random = ~1|site,
family = hurdle.poisson(),
zi_fixed = ~ mined + spp + cover,
data = mf.poisson)
summary(model.hardle)
Call:
mixed_model(fixed = count ~ mined + spp + cover, random = ~1 |
site, data = mf.poisson, family = hurdle.poisson(), zi_fixed = ~mined +
spp + cover)
Data Descriptives:
Number of Observations: 644
Number of Groups: 23
Model:
family: hurdle poisson
link: log
Fit statistics:
log.Lik AIC BIC
-874.2537 1786.507 1808.082
Random effects covariance matrix:
StdDev
(Intercept) 0.1251081
Fixed effects:
Estimate Std.Err z-value p-value
(Intercept) -0.1345 0.1931 -0.6965 0.48613578
minedno 1.1990 0.1724 6.9550 < 1e-04
sppPR -0.5452 0.2781 -1.9609 0.04989212
sppDM 0.2235 0.1446 1.5456 0.12220016
sppEC-A -0.1989 0.2009 -0.9901 0.32211763
sppEC-L 0.6266 0.1320 4.7466 < 1e-04
sppDES-L 0.6059 0.1300 4.6594 < 1e-04
sppDF 0.0470 0.1534 0.3060 0.75959241
cover -0.2078 0.0621 -3.3481 0.00081362
Zero-part coefficients:
Estimate Std.Err z-value p-value
(Intercept) 1.8034 0.2936 6.1414 < 1e-04
minedno -2.4951 0.2551 -9.7827 < 1e-04
sppPR 1.6800 0.3966 4.2359 < 1e-04
sppDM -0.4272 0.3504 -1.2194 0.2227040
sppEC-A 1.1055 0.3686 2.9995 0.0027046
sppEC-L -0.4272 0.3504 -1.2194 0.2227040
sppDES-L -0.6720 0.3520 -1.9093 0.0562246
sppDF -0.4272 0.3504 -1.2194 0.2227040
cover 0.0771 0.1212 0.6364 0.5245160
Integration:
method: adaptive Gauss-Hermite quadrature rule
quadrature points: 11
Optimization:
method: hybrid EM and quasi-Newton
converged: TRUE This tutorial provided an overview of Generalized Linear Mixed Models (GLMMs) using Adaptive Gaussian Quadrature (AGQ) with the {GLMMadaptive} package in R. We discussed the mathematical framework behind GLMMs, the likelihood function, numerical integration using AGQ, parameter estimation, and prediction. We demonstrated how to fit GLMMs using the {GLMMadaptive} package with examples of binary, Poisson, and negative binomial response data. The package provides a flexible and efficient tool for fitting GLMMs with various response distributions, complex random effects structures, and custom link functions. The package supports multiple random effects terms, custom link functions, and robust estimation methods. The {GLMMadaptive} package is a valuable resource for researchers analyzing complex data structures with non-normal outcomes.