Multilevel Regression Modeling of Covariance Matrix Outcomes

Technical Summary: Multilevel Regression Modeling of Covariance Matrix Outcomes

Problem Statement
Covariance matrix outcomes are central to neuroimaging research, particularly in studying brain functional connectivity (FC) via fMRI. A primary challenge in analyzing these outcomes is the hierarchical structure of modern datasets, such as the Human Connectome Project (HCP) Lifespan Study, where subjects are nested within clusters (e.g., age cohorts). Existing covariance regression methods, such as Covariate-Assisted Principal (CAP) regression, typically operate in a single-level framework. These methods either vectorize the covariance matrix (ignoring positive-definiteness constraints and increasing dimensionality) or assume a common projection direction across all subjects. The latter assumption is often too restrictive for lifespan studies, where the underlying neural subnetworks governing connectivity may vary structurally across different developmental stages (childhood, adulthood, aging). There is a need for a framework that can model covariance matrices as outcomes while accommodating between-cluster heterogeneity in the covariance-associated subspaces and borrowing strength across clusters.

Methodology: Multilevel Covariate-Assisted Principal Regression (MCAP)
The authors propose the Multilevel Covariate-Assisted Principal Regression (MCAP) framework to address these challenges. The core of the method involves modeling the covariance matrix $\Sigma_{ij}$ for unit $j$ in cluster $i$ through a linear projection $\gamma_i$ on the unit sphere.

• Model Structure: The model assumes that for each cluster $i$, there exists a unit-norm projection vector $\gamma_i$ such that the log of the projected variance follows a generalized linear mixed effects model:
  $$\log(\gamma_i^\top \Sigma_{ij} \gamma_i) = \beta_0 + x_{1i(j)}^\top \beta_1 + x_{2ij}^\top (\beta_2 + \vartheta_i) + \varepsilon_i$$
  Here, $x_{1i(j)}$ represents fixed-effect covariates (cluster or unit level), and $x_{2ij}$ represents random-effect covariates. The term $\vartheta_i$ captures cluster-specific deviations in slopes.
• Hierarchical Modeling of Projections: Unlike single-level CAP, MCAP allows the projection direction $\gamma_i$ to vary across clusters. These cluster-specific directions are modeled as random variables on the unit sphere following a von Mises-Fisher (vMF) distribution with a mean direction $\gamma$ and concentration parameter $\kappa$. This formulation enables principled information borrowing across clusters while permitting heterogeneity.
• Estimation: Parameters are estimated by maximizing a hierarchical likelihood function. The likelihood comprises four components: the conditional likelihood of unit-level observations, the likelihood of random intercepts, the likelihood of random slopes, and the likelihood of the random rotations (vMF distribution).
  • Algorithm: A block coordinate-descent algorithm is employed. Updates for regression coefficients and variance components use the Newton-Raphson method. Updates for the projection directions $\gamma_i$ involve solving a generalized eigenvalue problem with a data-adaptive normalization constraint ($\tilde{\gamma}_i^\top H_i \tilde{\gamma}_i = 1$), where $H_i$ is the cluster's average sample covariance matrix. This geometric constraint aligns the feasible set with the data cloud.
  • Component Identification: Higher-order components are identified sequentially via a deflation strategy. To ensure computational tractability when sample covariance matrices become rank-deficient after deflation, a rank-completion step replaces zero singular values with estimated random intercepts.
• Inference: A two-stage bootstrap procedure is proposed for inference on regression coefficients. In Stage 0, the full model is fitted to obtain estimates of $\gamma_i$. In Stage 1, clusters are resampled with replacement, and the regression sub-problem is solved with $\gamma_i$ fixed at their Stage 0 estimates. This approach respects the hierarchical data structure and avoids the computational burden of re-estimating directions in every bootstrap replicate.
• Asymptotic Theory: The paper establishes the asymptotic normality of regression parameter estimators when projection directions are known. When directions are unknown, the consistency of all estimators (including $\gamma_i$, $\gamma$, and regression coefficients) is proven under specific identifiability conditions, assuming the number of clusters ($m$) and observations per cluster ($n, T$) grow.

Key Results

• Simulation Studies: MCAP was compared against a single-level CAP approach (SCAP) that averages cluster-specific estimates. In simulations with dimensions $p=5$ and $p=20$, MCAP demonstrated substantially lower bias and mean squared error (MSE) in estimating regression coefficients, particularly as sample sizes increased. While SCAP performed adequately in low dimensions, it failed to accurately estimate coefficients in higher dimensions ($p=20$). MCAP effectively leveraged the hierarchical structure to improve estimation stability.
• HCP Lifespan Study Application: Applied to 1,500 subjects from the HCP Lifespan Study (ages 5–90), MCAP identified a single dominant spectral brain network capturing age and sex effects.
  • Lifespan Trajectory: The projected functional connectivity (FC) score exhibited a non-linear trajectory: declining from childhood into the late 20s, rising through midlife to a peak in the early 70s, and declining again in late adulthood (70–90).
  • Sex Differences: Males consistently showed higher projected FC values than females across the age range.
  • Neural Reorganization: The trajectory suggested a "mirroring" effect, where the decline in late adulthood partially echoed the developmental decline observed in early life. Furthermore, the spatial organization of the subnetwork became more homogeneous across neighboring age bins in later life (post-50), suggesting a convergence of neural reorganization patterns in late adulthood.
  • Regional Contributions: The identified subnetwork was driven by regions in the rostral middle frontal gyrus (rMFG), middle temporal gyrus (MTG), inferior parietal cortex (IPC), and isthmus cingulate gyrus (ICG), spanning frontoparietal, default mode, ventral-attention, and lateral temporal systems.

Significance and Claims
The authors claim that MCAP provides a flexible and statistically rigorous framework for analyzing hierarchically nested covariance matrix outcomes. By relaxing the common-projection assumption, MCAP allows for the modeling of cluster-specific covariance subspaces while still estimating a population-level direction. This is critical for lifespan studies where neurobiological mechanisms differ across developmental stages.

The paper asserts that MCAP outperforms single-level competitors in estimating regression coefficients, particularly in moderate-to-high dimensional settings. The application to the HCP study demonstrates the method's ability to uncover complex, non-linear lifespan profiles of brain connectivity and identify distributed cross-network patterns that would be difficult to characterize using independent cluster analyses or fixed projection directions. The authors note that the framework is specifically designed to handle the positive-definiteness constraint of covariance matrices without resorting to vectorization, thereby preserving interpretability and statistical power.

Limitations Acknowledged
The authors modestly acknowledge several limitations:

1. The model assumes Gaussian observations; robust or quasi-likelihood versions would be needed for non-Gaussian data.
2. The vMF distribution is unimodal; mixture models might be required if projection directions form multiple distinct clusters.
3. Asymptotic normality for the projection directions ($\gamma_i$ and $\gamma$) is not established, only consistency.
4. Higher-order components are identified sequentially via deflation, which may be less efficient than joint multi-component procedures.