glmmFEL fits generalized linear mixed models (GLMMs) with normal random effects using a matrix-based interface where users supply \((y, X, Z)\) directly. Model fitting is performed with an EM algorithm whose E-step can be approximated using first-order Laplace or fully exponential Laplace (mean-only or mean + covariance corrections), and includes pseudo-likelihood alternatives based on working-response / working-weights linearization.
This development branch is intentionally matrix-first and supports a flexible (including multiple-membership) random-effects design matrix \(Z\); see the package NOTE for details on current covariance support and planned extensions.
Supported families in this branch:
family = stats::binomial(link = "probit") (binary probit),
family = stats::binomial(link = "logit") (binary logit),
family = stats::poisson(link = "log") (Poisson log-link).
This matrix-only development branch is a streamlined rewrite based on the binary/Poisson engines in mvglmmRank (which were based on Karl, Yang, and Lohr, 2014). The fully exponential Laplace EM strategy is described in Karl, Yang, and Lohr (2014) and related fully exponential Laplace work such as Tierney, Kass, and Kadane (1989) and Rizopoulos, Verbeke, and Lesaffre (2009); the FE-within-EM lineage is attributed to Steele (1996). Karl, Yang, and Lohr (2014) referenced the source code of package JM when writing their own code.
At present, the package supports a single random-effect vector (i.e., one variance component). However, the theory in the references above—and the joint Poisson-binary model already implemented in mvglmmRank—extends naturally to a block-diagonal random-effects covariance matrix \(G\). This would allow multiple independent random effects and/or random intercept–slope specifications, with intercepts and slopes either correlated or uncorrelated (depending on the block structure). Enabling this functionality primarily involves allowing the user to specify a \(G\) structure in the main fitting function. The maintainer has working prototype code for selected \(G\) structures, but it has not yet been tested sufficiently for release; future support may be added. If you would like a particular \(G\) structure supported, please email the package maintainer with a request.
Importantly, the current implementation already supports a multiple-membership (or otherwise arbitrary) random-effects design matrix \(Z\), which is more general than what is available in some other implementations. In particular, the multi-membership setting allows more than one random effect from the same variance component to be active on the same observation, possibly with different weight entries per random effect; see Karl, Yang, and Lohr (2014) for details.
The RSPL/MSPL pseudo-likelihood code paths are adapted from the
RealVAMS implementation described in Broatch, Green, and Karl (2018), which
follows Wolfinger and O'Connell (1993) closely (working response + working
weights).
OpenAI's GPT models (such as GPT-5 Pro) were used to assist with coding and roxygen documentation; all content was reviewed and finalized by the author.
glmmFEL() supports:
"Laplace": first-order Laplace approximation,
"FE_mean": fully exponential Laplace corrections to \(\widehat\eta\) only,
"FE_full" (or "FE"): fully exponential Laplace corrections to both
\(\widehat\eta\) and \(\widehat{\mathrm{Var}}(\eta\mid y)\),
"RSPL" / "MSPL": restricted/marginal pseudo-likelihood
(working response / working weights).
glmmFEL() returns an object of class "glmmFELMod" containing:
beta: fixed-effect estimates,
eta: empirical Bayes predictions of random effects,
tau2: the scalar variance component,
G: \(q\times q\) covariance matrix (diagonal in this branch),
var_eta: prediction-error covariance for eta (approx.),
vcov_beta: approximate covariance of beta when available,
convergence: iteration counts and flags.
Maintainer: Andrew T. Karl akarl@asu.edu (ORCID)
Broatch, J., Green, J. G., & Karl, A. T. (2018). RealVAMS: An R Package for Fitting a Multivariate Value-added Model (VAM). The R Journal, 10(1), 22--30. tools:::Rd_expr_doi("10.32614/RJ-2018-033")
Karl, A. T., Yang, Y., & Lohr, S. L. (2014). Computation of maximum likelihood estimates for multiresponse generalized linear mixed models with non-nested, correlated random effects. Computational Statistics & Data Analysis, 73, 146--162. tools:::Rd_expr_doi("10.1016/j.csda.2013.11.019")
Rizopoulos, D., Verbeke, G., & Lesaffre, E. (2009). Fully exponential Laplace approximations in joint models for longitudinal and survival data. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 71(3), 637--654. tools:::Rd_expr_doi("10.1111/j.1467-9868.2008.00704.x")
Rizopoulos, D. (2010). JM: An R package for the joint modeling of longitudinal and time-to-event data. Journal of Statistical Software, 35(9), 1--33. tools:::Rd_expr_doi("10.18637/jss.v035.i09")
Steele, B. M. (1996). A modified EM algorithm for estimation in generalized mixed models. Biometrics, 52(4), 1295--1310. tools:::Rd_expr_doi("10.2307/2532845")
Tierney, L., Kass, R. E., & Kadane, J. B. (1989). Fully exponential Laplace approximations to expectations and variances of nonpositive functions. Journal of the American Statistical Association, 84(407), 710--716. tools:::Rd_expr_doi("10.1080/01621459.1989.10478824")
Wolfinger, R., & O'Connell, M. (1993). Generalized linear mixed models: a pseudo-likelihood approach. Journal of Statistical Computation and Simulation, 48(3--4), 233--243. tools:::Rd_expr_doi("10.1080/00949659308811554")
glmmFEL_pl() for
the pseudo-likelihood engines.