[Paper Review] Generalised Linear Mixed Model Specification, Analysis, Fitting, and Optimal Design in R with the glmmr Packages

Sam Watson|arXiv (Cornell University)|Mar 22, 2023

Optimal Experimental Design MethodsDecision Sciences3 citations

TL;DR

This paper introduces the R packages glmmrBase and glmmrOptim, offering a flexible, object-oriented framework for specifying, fitting, and analyzing generalised linear mixed models (GLMMs) using Markov Chain Monte Carlo Maximum Likelihood (MCML), Laplace approximation, and Hamiltonian Monte Carlo. The key contribution is a unified, extensible system that supports complex covariance structures, non-linear mean functions, and c-optimal experimental design, with benchmarked performance and accurate parameter estimates comparable to established packages like lme4 and glmmTMB.

ABSTRACT

We describe the \proglang{R} package \pkg{glmmrBase} and an extension \pkg{glmmrOptim}. \pkg{glmmrBase} provides a flexible approach to specifying, fitting, and analysing generalised linear mixed models. We use an object-orientated class system within \proglang{R} to provide methods for a wide range of covariance and mean functions, including specification of non-linear functions of data and parameters, relevant to multiple applications including cluster randomised trials, cohort studies, spatial and spatio-temporal modelling, and split-plot designs. The class generates relevant matrices and statistics and a wide range of methods including full likelihood estimation of generalised linear mixed models using stochastic Maximum Likelihood, Laplace approximation, power calculation, and access to relevant calculations. The class also includes Hamiltonian Monte Carlo simulation of random effects, sparse matrix methods, and other functionality to support efficient estimation. The \pkg{glmmrOptim} package implements a set of algorithms to identify c-optimal experimental designs where observations are correlated and can be specified using the generalised linear mixed model classes. Several examples and comparisons to existing packages are provided to illustrate use of the packages.

Motivation & Objective

• To address the fragmentation in existing GLMM software, where users often need multiple packages with inconsistent interfaces for a single analysis.
• To develop a unified, extensible R framework that supports complex GLMM specifications, including non-linear functions of data and parameters.
• To enable efficient estimation via MCML, Laplace approximation, and Hamiltonian Monte Carlo, with support for sparse matrices and high-performance C++ backends.
• To implement c-optimal experimental design for correlated observations using the GLMM framework.
• To provide a consistent, object-oriented interface using R6 classes that allows shared covariance and mean structures across multiple models.

Proposed method

• The glmmrBase package implements a Model class using R6's object-oriented system, encapsulating model components such as fixed and random effects, link functions, and distributional families.
• The framework supports user-defined non-linear functions for both the linear predictor and covariance structures, enabling flexible modeling of clustered, cohort, spatial, and spatio-temporal data.
• Model fitting is performed using Laplace approximation for initial parameter estimates, followed by Markov Chain Monte Carlo Maximum Likelihood (MCML) with Hamiltonian Monte Carlo for posterior sampling.
• The system leverages the Matrix package for sparse matrix operations and integrates C++ code for performance-critical computations.
• The glmmrOptim package extends the framework to compute c-optimal designs for experiments with correlated observations, using the GLMM structure to identify optimal allocation of observations.
• The implementation includes methods for full likelihood estimation, standard error computation (robust and bias-corrected), and power calculations.

Figure 1: Illustrative example of formula parsing and calculation including autodifferention with respect to parameter \code b_1

Experimental results

Research questions

• RQ1 How can a unified, extensible framework be designed in R to support complex GLMM specifications, including non-linear functions of data and parameters, across diverse applications?
• RQ2 What is the performance and accuracy of MCML and Laplace approximation methods when implemented in a modular, object-oriented R/C++ system for GLMMs?
• RQ3 Can the same GLMM framework be used to derive c-optimal experimental designs for correlated observations, and how does it compare to existing design methods?
• RQ4 How do parameter estimates and standard errors from the glmmrBase package compare to those from established packages like lme4 and glmmTMB in simulated settings?
• RQ5 To what extent does the object-oriented design of glmmrBase enable code reuse, modularity, and extensibility in GLMM-based statistical workflows?

Key findings

• Parameter estimates from glmmrMCML using MCML and Laplace approximation methods were within 0.01 of estimates from lme4 and glmmTMB, demonstrating high accuracy.
• The glmmrMCML implementation achieved parameter estimates consistent with reference methods, with standard errors differing by less than 0.02 across all tested parameters.
• The Laplace approximation with Newton-Raphson and BOBYQA optimisation methods provided stable starting values for MCML, reducing convergence issues.
• The glmmrOptim package successfully computed c-optimal designs for correlated observations, enabling efficient experimental planning in complex settings.
• Model fitting times varied significantly: MCNR took 7.4 seconds, while MCEM with BOBYQA reduced time to 0.36 seconds, showing the impact of optimisation choice on efficiency.
• The object-oriented design allowed shared covariance and mean structures across multiple models, improving code maintainability and reducing redundancy.