Package {glmmFEL}

Type:	Package
Title:	Generalized Linear Mixed Models via Fully Exponential Laplace in EM
Version:	1.0.6
Date:	2026-09-18
Description:	Fit generalized linear mixed models (GLMMs) with normal random effects using first-order Laplace, fully exponential Laplace (FEL) with mean-only corrections, and FEL with mean and variance-diagonal corrections in the E-step of an expectation-maximization (EM) algorithm. The current development version provides a matrix-based interface (y, X, Z) and supports binary logit and probit, and Poisson log-link models. An EM framework is used to update fixed effects, random effects, and a single variance component tau^2 for G = tau^2 I, with staged approximations (Laplace -> FEL mean-only -> FEL full) for efficiency and stability. A pseudo-likelihood engine glmmFEL_pl() implements the working-response / working-weights linearization approach of Wolfinger and O'Connell (1993) <doi:10.1080/00949659308811554>, and is adapted from the implementation used in the 'RealVAMS' package (Broatch, Green, and Karl (2018)) <doi:10.32614/RJ-2018-033>. The FEL implementation follows Karl, Yang, and Lohr (2014) <doi:10.1016/j.csda.2013.11.019> and related work (e.g., Tierney, Kass, and Kadane (1989) <doi:10.1080/01621459.1989.10478824>; Rizopoulos, Verbeke, and Lesaffre (2009) <doi:10.1111/j.1467-9868.2008.00704.x>; Steele (1996) <doi:10.2307/2532845>). Package code was drafted with assistance from generative AI tools.
License:	GPL-3
Encoding:	UTF-8
Imports:	Matrix, stats, methods
Suggests:	testthat (≥ 3.0.0), lme4, numDeriv
Config/testthat/edition:	3
NeedsCompilation:	no
Config/roxygen2/version:	8.1.0
Packaged:	2026-09-18 17:04:56 UTC; andre
Author:	Andrew T. Karl [cre, aut]
Maintainer:	Andrew T. Karl <akarl@asu.edu>
Repository:	CRAN
Date/Publication:	2026-09-18 17:20:02 UTC

glmmFEL: Generalized Linear Mixed Models via Fully Exponential Laplace in EM

Description

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 and variance-diagonal corrections), and includes pseudo-likelihood alternatives based on working-response / working-weights linearization.

The matrix interface supports multiple-membership designs with a single variance component. It does not implement all covariance structures in the cited papers. See glmmFEL() for the precise approximation and output limits.

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).

Approximations

glmmFEL() supports:

• "Laplace": first-order Laplace approximation,

• "FE_mean": fully exponential Laplace corrections to \widehat\eta only,

• "FE_full" (or "FE"): mean corrections plus the posterior variance-diagonal corrections needed for the scalar variance update,

• "RSPL" / "MSPL": restricted/marginal pseudo-likelihood (working response / working weights).

Output

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: approximate posterior covariance, with diagonal-only FEL covariance corrections in FE_full,

• vcov_beta: approximate covariance of beta when available,

• convergence: iteration counts and flags.

Author(s)

Maintainer: Andrew T. Karl akarl@asu.edu (ORCID)

Authors:

• Andrew T. Karl akarl@asu.edu (ORCID)

References

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. 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. 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. 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. 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. 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. 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. doi:10.1080/00949659308811554

See Also

glmmFEL_pl() for the pseudo-likelihood engines; glmmFEL-benchmarks for the simulation appraisal.

Extract model coefficients (fixed effects)

Description

Extract model coefficients (fixed effects)

Usage

## S3 method for class 'glmmFELMod'
coef(object, ...)

Arguments

Value

A named numeric vector of estimated fixed-effect regression coefficients (on the linear predictor scale). Names correspond to columns of the fixed-effects design matrix X when available.

Extract fitted values

Description

Extract fitted values

Usage

## S3 method for class 'glmmFELMod'
fitted(object, type = c("response", "link"), ...)

Arguments

Value

A numeric vector of fitted values. If type = "link", the fitted linear predictor values are returned. If type = "response", the fitted mean response values on the response scale are returned. The length equals the number of observations in the fitted object.

Fit GLMMs via Laplace and fully exponential Laplace (matrix interface)

Description

glmmFEL() fits a generalized linear mixed model (GLMM) with multivariate normal random effects using EM-type algorithms and likelihood approximations:

• first-order Laplace (approx = "Laplace"),

• fully exponential corrections to the random-effects mean (approx = "FE_mean"),

• fully exponential corrections to both mean and variance diagonals (approx = "FE_full" / "FE"),

• pseudo-likelihood / PL linearization (approx = "RSPL" or "MSPL") via glmmFEL_pl().

The interface is matrix-based: provide the response y, fixed-effects design matrix X, and random-effects design matrix Z.

Random effects are assumed \eta \sim N(0, G) with a single variance component

G = \tau^2 I_q,

allowing arbitrary (including multi-membership) Z while keeping the variance update simple and stable.

Usage

glmmFEL(
  y,
  X,
  Z,
  family = stats::binomial(link = "probit"),
  approx = c("FE", "Laplace", "FE_mean", "FE_full", "RSPL", "MSPL"),
  max_iter = 200,
  tol = 1e-06,
  control = list()
)

Arguments

Details

Laplace denotes an EM algorithm using posterior modes in its expected score and Gaussian posterior second moments. It is not direct maximization of the first-order Laplace marginal likelihood (as in lme4::glmer). FE_mean adds mean and expected-score corrections. FE_full also corrects the posterior variance diagonals needed to update the single variance component; off-diagonal entries of var_eta remain at their Laplace values. Thus "full" distinguishes the two implemented EM stages, not a correction of every entry of an arbitrary posterior covariance matrix.

Returned moments are recomputed at the final parameter estimates. A failed inner solve or exhausted iteration budget produces a warning and a false convergence flag. Always inspect convergence before interpreting estimates. Approximate posterior variances can become invalid in sparse/extreme settings; such fits are flagged rather than reported as converged. A strictly separating fixed-effect direction found at initialization is rejected for binary data, since no finite estimate exists. This sufficient check is not a complete detector of all forms of quasi-separation.

vcov_beta is a joint-Hessian approximation conditional on the fitted variance component; it is not a fully exponential observed-information estimate and does not incorporate variance-component estimation uncertainty. For EM fits, logLik evaluates the first-order Laplace marginal likelihood at the returned estimates, which need not maximize it. For PL fits it is unavailable (NA); working_logLik stores the working Gaussian objective. These objectives should not be used for cross-method AIC comparisons. Fitted response values plug the estimated random effects into the inverse link; they are not posterior predictive averages over random effects.

See glmmFEL-benchmarks for a reproducible simulation appraisal and its limits.

Value

A fitted model object of class glmmFELMod, with estimates beta, tau2, random-effect predictions eta, approximate covariance matrices, and convergence. EM convergence diagnostics include em_converged, reason, mode_converged, mode_gradient, beta_converged, and moments_valid. PL fits report both pql_converged and em_converged.

References

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. 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. 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. 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. 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. 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. 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. doi:10.1080/00949659308811554

Examples

## Example 1: Simulated probit random-intercept GLMM (matrix interface)
set.seed(1)
n_id <- 30
m_per_id <- 6
n <- n_id * m_per_id
id <- factor(rep(seq_len(n_id), each = m_per_id))
x  <- rnorm(n)
X  <- model.matrix(~ x)
Z  <- Matrix::sparseMatrix(i = seq_len(n),
                           j = as.integer(id),
                           x = 1,
                           dims = c(n, n_id))
beta_true <- c(0.2, 0.7)
tau2_true <- 0.5
eta_true  <- rnorm(n_id, sd = sqrt(tau2_true))
lp <- as.vector(X %*% beta_true + Z %*% eta_true)
y  <- rbinom(n, 1, pnorm(lp))

fit <- glmmFEL(y, X, Z, family = "binomial_probit", approx = "Laplace")
fit$beta
fit$tau2

## Example 2: Get X, y, Z from an lme4 formula (without a glmmFEL formula wrapper)

if (requireNamespace("lme4", quietly = TRUE)) {
  dat <- data.frame(y = y, x = x, id = id)
  lf  <- lme4::lFormula(y ~ x + (1 | id), data = dat)
  X_lme4 <- lf$X
  Z_lme4 <- Matrix::t(lf$reTrms$Zt)
  y_lme4 <- lf$fr$y

  fit2 <- glmmFEL(y_lme4, X_lme4, Z_lme4, family = "binomial_probit", approx = "Laplace")
}

Simulation appraisal of glmmFEL approximations

Description

A fixed-design simulation comparison of estimation accuracy, convergence, and computing time. The results describe the settings studied, not a general ranking of GLMM estimators.

Details

These static results were generated for version 1.0.6 on 2026-09-18 from a prespecified design: 2,400 datasets and 22,400 applicable fits. The two glmmML methods have no probit implementation; their 1,600 non-applicable rows are recorded separately. No simulation runs during installation or package checks.

Design

Each family/setting combination has 200 independent datasets, paired across methods. The covariate is standard normal; Gaussian random intercepts have variance 0.5. The true slope is 0.5, and the intercept is -0.5 for binary responses and log(2) for Poisson. Binary responses are Bernoulli observations.

Methods and diagnostics

Modal EM, FEL mean/full, and MSPL/RSPL are the package's five algorithms. The package's Laplace uses modes in EM; it is not direct marginal Laplace optimization. lme4 uses glmer with nAGQ=1 or 9. glmmML uses Laplace or 20-point Gaussian quadrature. Quadrature is an implementation reference, not simulated truth. MASS::glmmPQL uses its ML working-model convention with residual scale fixed at one.

glmmFEL EM has a 400-iteration budget and 1e-6 tolerance for all three stages. Its PL fits have 100 outer and 100 inner iterations. MASS PQL has 100 outer iterations; its reported convergence is a proxy based on returning before the cap, because the API exposes no outer convergence flag. Boundary fits are reported separately from convergence. An iteration cap counts as failure, even if the numerical estimates might already be practically close.

Interpretation

FEL was not an overall winner in these settings. It often reduced variance-component bias relative to modal EM or pseudo-likelihood, but the accompanying sampling variance, iteration cost, and failures matter. Bias reduction did not imply a smaller RMSE. With well-informed clusters, FEL estimates were generally close to quadrature estimates, including binary settings.

Sparse binary data were unfavorable. For logit responses, FEL full converged in 80/200 small-cluster datasets and 53/200 many-small-cluster datasets at the 400-iteration budget. Its variance squared error exceeded AGQ9 on their common converged datasets. One sparse probit fit produced a variance estimate around 1.7e155 before being flagged for invalid posterior moments. Failed fits must not be used for inference, even when an estimate is returned.

Rich Poisson data exposed slow staged EM: only 72/200 FEL-full fits converged within 400 iterations. The post-hoc sensitivity increased only the budget to 2,000 on the first 40 seeds of every cell. In rich Poisson data, FEL-full convergence rose from 12/40 to 40/40. Its variance RMSE was 0.10449, versus 0.10450 for both lme4 AGQ9 and lme4 Laplace, and 0.10434 for MASS PQL. Slope RMSE was approximately 0.02545 for all four. The FEL-full minus AGQ9 variance MSE difference was -2.09e-6 (MCSE 1.11e-6): a very small practical difference in a 40-replicate sensitivity, not evidence of broad superiority.

Thus the study supports accurate FEL approximation in the tested Poisson settings, but does not establish an advantage over these competing solvers, even with rich data. Increasing the budget also improved binary convergence, yet sparse-logit failures and excess variance error remained. The independent random-intercept structure particularly favors available quadrature solvers; the small crossed-design diagnostic does not settle accuracy for the broader nonnested problems motivating the literature.

Convergence across all settings

Counts below are out of 200 applicable datasets per cell. NA means that the method does not implement the link, not optimizer failure.

Bernoulli-logit

Bernoulli-probit

Poisson-log

Variance estimation with rich clusters

For the 30-by-20 setting, the following are bias, empirical sampling variance, and RMSE of the estimated random-intercept variance (truth 0.5). Each method uses its own converged subset; use the counts above and the paired contrasts below. Complete results for every setting and all three parameters, including Monte Carlo SEs, are in the accompanying CSV files.

Bernoulli-logit

Bernoulli-probit

Poisson-log

Paired variance comparisons

FEL full minus lme4 AGQ9 squared-error differences, on common converged datasets only. Negative differences favor FEL full. MCSE is the Monte Carlo standard error of that difference. Small differences relative to MCSE do not support a stable ranking; these are exploratory, multiple comparisons.

Post-hoc iteration-budget sensitivity

This smaller extension uses the same first 40 seeds in every cell, chosen by replicate number rather than outcome. It increases only the three EM methods' budgets to 2,000 and reuses unchanged main-run comparator fits. Below are FEL-full convergence counts (out of 40) and paired variance MSE differences versus AGQ9. It is exploratory and does not replace the 200-replicate primary study. Other methods and parameters are in the separate sensitivity CSV files.

Limits and reproducibility

This study varies cluster information and balance at one effect size and one random-effect variance. It does not establish performance for rare outcomes, larger variances, misspecified random effects, grouped binomial data, or arbitrary nonnested designs. A separate 24-dataset crossed/signed-membership diagnostic had 100 converged fits among 120 attempts; all 20 failures were iteration caps. Four seeds per family/design do not support accuracy claims.

Main accuracy summaries condition on convergence. Pairing successful fits does not remove selection bias from failures. The files also report all finite returned estimates separately, including nonconverged values. Wald coverage is exploratory: the package's Hessian standard errors condition on the estimated variance and are not full FEL observed-information standard errors. Random-effect prediction loss compares the means or modes each method returns.

Timings include contention from eight local workers; they are not isolated microbenchmarks. The static files under system.file("benchmarks", package="glmmFEL") include full summaries, paired contrasts, failure counts, package versions, the scenario grid, source/seed provenance, and column definitions. The maintainer retains raw per-fit results and the frozen simulation scripts.

See Also

glmmFEL()

Examples

results <- read.csv(system.file("benchmarks", "summary.csv", package = "glmmFEL"))
subset(results, family == "poisson_log" & setting == "rich" & parameter == "tau2",
       select = c(method, converged, bias, variance, rmse, bias_mcse))

Pseudo-likelihood engine for RSPL/MSPL (Wolfinger-style, simplified R = I)

Description

This follows the vp_cp / RealVAMS structure:

• Outer PQL loop updates working response (z_work) and weights (w_num).

• Inner EM loop (with z_work, w_num fixed) updates (beta, eta, tau2).

Usage

glmmFEL_pl(
  y,
  X,
  Z,
  family = "binomial_probit",
  approx = c("RSPL", "MSPL"),
  max_iter = 200L,
  tol = 1e-06,
  control = list()
)

Details

RSPL vs MSPL: #' - MSPL uses \mathrm{Var}(\eta \mid \beta, y) = (Z^\top W Z + G^{-1})^{-1} (called H.inv in vp_cp).

• RSPL uses the eta-eta block of the inverse of the full augmented system (called C.mat in vp_cp) for the variance-component moment update.

Coerce a fixed-effects design matrix to a numeric base matrix

Description

Coerce a fixed-effects design matrix to a numeric base matrix

Usage

glmmfe_as_X(X)

Arguments

Value

A base numeric matrix.

Coerce a random-effects design matrix to a sparse dgCMatrix

Description

Coerce a random-effects design matrix to a sparse dgCMatrix

Usage

glmmfe_as_Z(Z, n = NULL)

Arguments

Value

A sparse dgCMatrix with numeric storage.

Internal Gaussian inner fit for PL / weighted LMM with G = tau2 * I

Description

Model (conditional on working quantities):

z_{\mathrm{work}} = X\beta + Z\eta + e,\qquad e \sim N(0, W^{-1})

\eta \sim N(0, \tau^2 I_q)

Usage

glmmfe_lmm_inner_fit(
  z_work,
  w_num,
  X,
  Z,
  tau2,
  approx = c("RSPL", "MSPL"),
  vc_eps = 1e-12,
  ridge_init = 1e-08
)

Details

where W = \mathrm{diag}(w_{\mathrm{num}}).

Returns beta, eta and covariance blocks needed by RSPL/MSPL updates. Does NOT form n \times n matrices.

Construct a glmmFEL fitted-model object

Description

Internal helper to standardize creation of the fitted-model object returned by glmmFEL() and related engines.

This branch stores only a single variance component tau2 with G = \tau^2 I_q. More structured covariance parameterizations and the formula wrapper are intentionally removed to reduce complexity.

Usage

glmmfe_new_fit(
  y,
  X,
  Z,
  beta,
  eta,
  tau2,
  G = NULL,
  vcov_beta = NULL,
  vcov_eta = NULL,
  cov_beta_eta = NULL,
  var_eta = NULL,
  family,
  approx,
  control,
  convergence,
  logLik = NA_real_,
  call = NULL,
  reml = NULL
)

Arguments

Value

A list of class c("glmmFELMod", "glmmFEL").

vp_cp-style PL objective (includes constants)

Description

vp_cp-style PL objective (includes constants)

Usage

glmmfe_pl_objective(
  z_work,
  w_num,
  X,
  Z,
  beta,
  eta,
  tau2,
  inner,
  approx,
  vc_eps = 1e-12
)

Resolve approximation labels

Description

Internal helper that normalizes approximation input to one of:

• "Laplace"

• "FE_mean"

• "FE_full" (alias "FE")

• "RSPL"

• "MSPL"

Usage

glmmfe_resolve_approx(approx)

Arguments

Value

Canonical approximation label.

Resolve supported family specifications

Description

Internal helper that normalizes user family input into one of the canonical labels used by the internal FE/Laplace code:

• "binomial_probit"

• "binomial_logit"

• "poisson_log"

Users may pass either the canonical character labels above, or common stats::family() objects:

• stats::binomial(link = "probit")

• stats::binomial(link = "logit")

• stats::poisson(link = "log")

Usage

glmmfe_resolve_family(family)

Arguments

Value

A canonical family label.

Fast trace of a matrix product

Description

Internal helper to compute \mathrm{tr}(A B) without forming A %*% B:

\mathrm{tr}(A B) = \sum (A \circ B^\top).

This identity is used repeatedly in the fully exponential (FE) trace corrections (Karl, Yang, and Lohr, 2014, Appendix B), where a naive diag(A %*% B) would allocate the full product.

Usage

glmmfe_trAB(A, B)

Arguments

Value

A single numeric scalar equal to tr(A %*% B).

Extract log-likelihood (approximate)

Description

Extract log-likelihood (approximate)

Usage

## S3 method for class 'glmmFELMod'
logLik(object, ...)

Arguments

Value

An object of class "logLik" giving the approximate log-likelihood stored in object$logLik (NA for PL fits), with attributes "df" (effective number of parameters, taken as length(beta) + 1 for the variance component) and "nobs" (number of observations).

Predict from a fitted glmmFEL model

Description

Predict from a fitted glmmFEL model

Usage

## S3 method for class 'glmmFELMod'
predict(object, newdata = NULL, type = c("response", "link"), ...)

Arguments

Value

A numeric vector of predictions. If type = "link", predictions are returned on the linear predictor scale. If type = "response", predictions are returned on the response scale. The length equals the number of observations used to fit the model. Supplying newdata triggers an error in this matrix-only branch.

Print a glmmFEL model object

Description

Print a glmmFEL model object

Usage

## S3 method for class 'glmmFELMod'
print(x, ...)

Arguments

Value

Returns x invisibly (a glmmFELMod object), called for its side effect of printing model information to the console.

Print a summary.glmmFELMod object

Description

Print a summary.glmmFELMod object

Usage

## S3 method for class 'summary.glmmFELMod'
print(x, ...)

Arguments

Value

Returns x invisibly (a summary.glmmFELMod object), called for its side effect of printing the summary to the console.

Summary for a glmmFEL model object

Description

Summary for a glmmFEL model object

Usage

## S3 method for class 'glmmFELMod'
summary(object, ...)

Arguments

Value

An object of class summary.glmmFELMod (a list) containing summary information for the fitted model, including the call, family, approximation label, approximate log-likelihood, estimated variance component tau2, a coefficient table for fixed effects (and standard errors when available), number of observations, and convergence information.

Extract the covariance matrix of the fixed effects

Description

Extract the covariance matrix of the fixed effects

Usage

## S3 method for class 'glmmFELMod'
vcov(object, ...)

Arguments

Value

A variance-covariance matrix for the estimated fixed-effect regression coefficients. Row and column names correspond to coefficient names when available.