Matrix-Response Generalized Linear Mixed Model with Applications to Longitudinal Brain Images

Longitudinal brain imaging data facilitate the monitoring of structural and functional alterations in individual brains across time, offering essential understanding of dynamic neurobiological mechanisms. Such data improve sensitivity for detecting early biomarkers of disease progression and enhance the evaluation of intervention effects. While recent matrix-response regression models can relate static brain networks to external predictors, there remain few statistical methods for longitudinal brain networks, especially those derived from high-dimensional imaging data. We introduce a matrix-response generalized linear mixed model that accommodates longitudinal brain networks and identifies edges whose connectivity is influenced by external predictors. An efficient Monte Carlo Expectation-Maximization algorithm is developed for parameter estimation. Extensive simulations demonstrate effective identification of covariate-related network components and accurate parameter estimation. We further demonstrate the usage of the proposed method through applications to diffusion tensor imaging (DTI) and functional MRI (fMRI) datasets.

💡 Research Summary

The paper introduces a novel statistical framework called the Matrix‑Response Generalized Linear Mixed Model (MR‑GLMM) designed specifically for longitudinal brain imaging data where the response at each time point is a square adjacency matrix representing a brain network. Traditional generalized linear mixed models (GLMMs) handle scalar or vector outcomes and cannot directly accommodate matrix‑valued responses, while existing matrix‑response regression methods are limited to cross‑sectional data and ignore subject‑specific random effects. MR‑GLMM fills this gap by modeling the conditional mean of each adjacency matrix through a canonical link function applied element‑wise, yielding the linear predictor

 η_it = Θ + θ_i + B ×₃ x_it,

where Θ ∈ ℝ^{n×n} is a population‑level fixed‑effect matrix assumed to be low‑rank, θ_i ∈ ℝ^{n×n} is a subject‑specific random‑intercept matrix with independent Gaussian entries, B ∈ ℝ^{n×n×p} is a three‑dimensional tensor encoding the effect of p covariates on each edge, and x_it is the covariate vector at time t for subject i. The low‑rank assumption reflects the empirical observation that brain connectivity can be summarized by a few latent network modes, while the sparsity assumption on B reflects the belief that only a subset of edges are truly associated with the covariates.

To enforce low‑rankness, the authors adopt the Burer‑Monteiro factorization Θ = UVᵀ with U, V ∈ ℝ^{n×r}, avoiding costly singular‑value decompositions at each iteration. Uniqueness of the factorization is encouraged by adding a penalty term ‖UᵀU − VᵀV‖_F². Sparsity of B is imposed via an ℓ₀‑constraint limiting the number of non‑zero entries to a pre‑specified proportion of the total edges. Although ℓ₀ optimization is NP‑hard, the authors embed it within a Monte Carlo Expectation‑Maximization (MCEM) algorithm that yields a practical approximation.

The MCEM algorithm proceeds as follows. In the E‑step, the posterior distribution of the high‑dimensional random effects θ_i given the observed data and current parameter estimates is intractable. The authors therefore generate M independent draws from this posterior using a Metropolis‑within‑Gibbs sampler, which efficiently explores the conditional Gaussian structure of θ_i while accounting for the non‑Gaussian likelihood contributed by the exponential‑family response. These draws are used to approximate the expected complete‑data log‑likelihood (the Q‑function). In the M‑step, the approximated Q‑function is minimized with respect to (U, V, B). The minimization is performed by alternating updates: (i) a gradient‑based update for U and V that includes the low‑rank penalty, and (ii) a hard‑thresholding step on B that retains only the s·n² largest absolute entries, thereby enforcing the ℓ₀ constraint. The algorithm iteratively increases the Monte Carlo sample size and reduces the step size to ensure stability and convergence.

Extensive simulations evaluate the method under varying network sizes (n = 50–100), numbers of time points (T = 3–5), covariate dimensions (p = 10–20), rank (r = 3–5), and sparsity levels (2–10% non‑zero entries). Performance metrics include variable‑selection F1‑score, Frobenius norm error for Θ, and Frobenius norm error for B. Across all scenarios, MR‑GLMM outperforms element‑wise penalized GLMMs that treat each edge independently, achieving 15–30 % higher F1‑scores and 20–40 % lower estimation errors, especially when random‑effect variance is large or the true B is highly sparse. The method’s robustness to the number of Monte Carlo samples is demonstrated by negligible performance change when M is increased from 500 to 2000.

The authors apply MR‑GLMM to two publicly available neuroimaging datasets. In the ADNI diffusion tensor imaging (DTI) study, 200 older adults are scanned up to three times over three years, yielding 68 × 68 structural connectivity matrices. Covariates include age, sex, APOE4 carrier status, and diagnostic group (normal, MCI, AD). The estimated low‑rank Θ captures dominant fronto‑temporal long‑range connections, while the sparse B identifies specific fronto‑hippocampal edges whose strength is modulated by APOE4 and disease status, corroborating known Alzheimer’s disease biomarkers. In the Human Connectome Project (HCP) functional MRI (fMRI) study, 300 young adults undergo resting‑state scans before and after a cognitive task, producing 100 × 100 functional connectivity matrices. The model reveals task‑induced strengthening of fronto‑striatal connections, aligning with prior neurophysiological findings about task‑related network reconfiguration.

Overall, the paper makes three major contributions: (1) it extends GLMM methodology to matrix‑valued longitudinal outcomes, thereby accommodating subject‑specific random effects and time‑varying covariates; (2) it integrates low‑rank and sparsity structures within a computationally efficient MCEM framework that leverages modern MCMC techniques; and (3) it demonstrates superior statistical performance and neuroscientific interpretability on both simulated and real brain‑imaging data. The authors suggest future extensions such as nonlinear link functions, multimodal integration (e.g., simultaneous structural and functional networks), and scalable implementations using parallel computing or GPU acceleration to handle thousands of subjects.