[Paper Review] Parsimoniously Fitting Large Multivariate Random Effects in glmmTMB

Maeve McGillycuddy, Gordana Popović+2 more|arXiv (Cornell University)|Nov 7, 2024

Statistical Methods and Inference4 citations

TL;DR

This paper introduces a reduced-rank approach in the glmmTMB R package to fit large-dimensional multivariate random effects by expressing them as linear combinations of fewer latent variables (d < q), enabling efficient estimation of complex mixed models with high-dimensional correlation structures. The method allows for scalable modeling of multivariate responses in ecological and social science data, overcoming convergence and computational limits of unstructured covariance matrices.

ABSTRACT

Multivariate random effects with unstructured variance-covariance matrices of large dimensions, $q$, can be a major challenge to estimate. In this paper, we introduce a new implementation of a reduced-rank approach to fit large dimensional multivariate random effects by writing them as a linear combination of $d &lt; q$ latent variables. By adding reduced-rank functionality to the package glmmTMB, we enhance the mixed models available to include random effects of dimensions that were previously not possible. We apply the reduced-rank random effect to two examples, estimating a generalized latent variable model for multivariate abundance data and a random-slopes model.

Motivation & Objective

• To address the computational and estimation challenges of fitting unstructured multivariate random effects with large dimensions (q) in generalized linear mixed models.
• To extend the glmmTMB package with reduced-rank factor analytic structures to handle high-dimensional multivariate random effects efficiently.
• To enable flexible modeling of complex study designs—such as multi-level and repeated measures—while maintaining parsimony and computational feasibility.
• To provide a practical solution for fitting generalized latent variable models (GLVMs) with large numbers of correlated responses, such as species abundance or survey items.
• To allow estimation of high-dimensional random effects (e.g., hundreds or thousands of species) that were previously infeasible with standard mixed model packages.

Proposed method

• Represent large-dimensional multivariate random effects (q) as a linear combination of d < q latent variables, reducing the number of parameters to estimate from q(q+1)/2 to d(q + d)/2.
• Use a factor analytic model structure where the random effect vector is expressed as η = Lξ, with ξ ~ N(0, I_d) and L a q×d loading matrix.
• Implement the model in glmmTMB using automatic differentiation and Laplace approximation for fast estimation of the marginal likelihood.
• Integrate the reduced-rank structure into glmmTMB’s existing interface, allowing seamless use with standard mixed model syntax and support for various exponential family families.
• Apply optimization strategies such as multiple starting values and information criteria (e.g., AIC) to select the optimal rank d.
• Allow for singular variance-covariance matrices by design, avoiding convergence issues common in positive-definite estimation when q is large.

Figure 1: Boxplot of fish abundance ( $\log(y+1)$ scale) for each species at three zones, windfarm (WF, green), north (N, orange), and south (S, lilac), before (2003) and after (2010) construction of the offshore wind farm.

Experimental results

Research questions

• RQ1 Can a reduced-rank factor analytic structure be effectively integrated into glmmTMB to estimate high-dimensional multivariate random effects?
• RQ2 How does the reduced-rank approach improve estimation stability and computational feasibility compared to unstructured covariance matrices for large q?
• RQ3 To what extent do fixed effect estimates and inferences change as the rank d of the latent variable space varies?
• RQ4 Can this method be applied to complex study designs such as multi-level or repeated measures with high-dimensional responses?
• RQ5 What strategies are effective for selecting the optimal rank d in practice, especially when the true dimensionality is unknown?

Key findings

• The reduced-rank approach in glmmTMB enables estimation of multivariate random effects with dimensions in the hundreds or even thousands, which were previously infeasible with standard packages like lme4 or glmmTMB without rank reduction.
• For a nine-species fish abundance dataset, the unstructured model failed to converge, but the reduced-rank model with d=2 successfully estimated correlations across species.
• In the PIRLS education study with 15 outcome variables, the reduced-rank model provided stable estimates and robust fixed effect inference, with results insensitive to moderate changes in d.
• The method allows for singular variance-covariance matrices by design, circumventing convergence warnings common in positive-definite estimation when q is large.
• Information criteria and visual inspection (e.g., biplots) were effective tools for selecting the optimal rank d, with qualitative changes in interpretation occurring primarily at low d (e.g., d < 2).
• The approach maintains compatibility with glmmTMB’s existing interface, enabling use in multi-level, repeated measures, and complex sampling designs previously unsupported by dedicated GLVM packages.

Figure 2: Conditional estimates (95% confidence interval) of the Zone by Year interaction terms for species from the diagonal random effect. The contrast between a zone (North, or South), and the Wind Farm zone in 2010 is shown.