Convex estimation of Gaussian graphical regression models with covariates

1. Problem Statement

The paper addresses the estimation of Gaussian Graphical Models (GGMs) in the presence of covariates. In standard GGMs, the goal is to recover the conditional independence structure among random variables $X \in \mathbb{R}^p$ (represented by the sparsity pattern of the precision matrix $\Omega = \Sigma^{-1}$).

In many applications, such as co-expression quantitative trait locus (eQTL) studies, the network structure itself may vary depending on external covariates $U \in \mathbb{R}^q$ (e.g., genetic variants). The authors consider a general model where both the mean vector $\mu(u)$ and the precision matrix $\Omega(u)$ depend on covariates $u$:
$$X | U = u \sim \mathcal{N}(\mu(u), \Omega(u)^{-1})$$

The Core Challenge:
Existing methods for covariate-adjusted GGMs (e.g., Zhang and Li [2023]) typically model the dependence of $\Omega(u)$ on $U$ via a linear expansion (e.g., $\Omega(u) = B_0 + \sum B_h u_h$). However, estimating the parameters jointly leads to a non-convex optimization problem due to the interaction between regression coefficients (for the mean) and precision matrix parameters.

• Current approaches rely on two-stage procedures (first estimating the mean, then centering the data to estimate the precision).
• These two-stage methods suffer from error propagation: inaccuracies in the first stage (mean estimation) degrade the performance of the second stage (graph estimation), leading to weaker theoretical guarantees and sensitivity to tuning parameters.

2. Methodology: The "cspine" Framework

The authors propose a new method called cspine (Covariate-adjusted Sparse Precision with Natural Estimation). The core innovation is a natural parametrization of the multivariate Gaussian likelihood that renders the joint estimation problem convex.

A. Natural Parametrization

Instead of estimating $\mu(u)$ and $\Omega(u)$ directly, the authors reparameterize the model using the natural parameters of the exponential family:

1. Define $\tilde{\theta} = \Omega \mu$ and $\tilde{\Theta} = -\Omega$.
2. The authors introduce a specific scaling to decouple the mean and precision effects:
   $$\theta(u) = \text{diag}(\Omega)^{-1} \Omega \mu(u) = \Gamma u$$
   $$\Theta(u) = -\text{diag}(\Omega)^{-1} \Omega(u) = B_0 + \sum_{h=1}^q B_h u_h$$
   Here, $\Gamma$ captures covariate effects on the mean, and $B_h$ capture covariate effects on the graph structure.

B. Convex Nodewise Regression

Using this parametrization, the conditional distribution of $X_j$ given $X_{-j}$ and $U$ can be written as a linear regression model:
$$X_j = u^\top \gamma_j + \sum_{k \neq j} \beta_{jk0} X_k + \sum_{k \neq j} \sum_{h=1}^q \beta_{jkh} u_h X_k + \epsilon_j$$
Crucially, the least squares criterion for this equation is jointly convex in the parameters $\gamma_j$ (mean effects) and $\beta_{jkh}$ (graph effects). This eliminates the need for a two-stage procedure.

C. Estimation Procedure

The method estimates parameters via nodewise regression with structured sparsity penalties:

1. Objective Function: For each node $j$, minimize:
   $$\frac{1}{2n} \| x_j - U\gamma_j - W_{-j}\beta_j \|_2^2 + \lambda \|\gamma_j\|_1 + \lambda \|\beta_j\|_1 + \lambda_g \|\beta_{j, -0}\|_{1,2}$$
   • $\ell_1$ penalties induce element-wise sparsity (for baseline connections and mean effects).
   • Group Lasso penalty ($\|\cdot\|_{1,2}$) induces group-wise sparsity, allowing the model to determine if a specific covariate $u_h$ affects any edge connected to node $j$.
2. Symmetrization: Since nodewise regressions are performed independently, the resulting precision matrix estimate is not symmetric. The authors apply a post-hoc "and-rule" symmetrization (an edge exists only if both directions are detected).
3. Recovery: The final estimates for $\mu(u)$ and $\Omega(u)$ are recovered via the inverse transformation of the natural parameters.

3. Key Contributions

1. Joint Convex Formulation: The primary theoretical contribution is the derivation of a parametrization that makes the joint estimation of mean and precision matrix parameters a convex optimization problem. This avoids the non-convexity inherent in previous formulations.
2. Improved Theoretical Guarantees:
   • The authors establish $\ell_2$ estimation error bounds and support recovery guarantees.
   • Scaling Conditions: The proposed method allows for a wider range of scaling regimes between the sample size $n$, dimension $p$, and covariate dimension $q$. Specifically, it permits $s \log(pq) = o(\sqrt{n})$ (where $s$ is sparsity), which is strictly weaker than the conditions required by the two-stage method of Zhang and Li [2023] (which required $s \log(pq) = o(n^{1/6})$).
   • Error Propagation: By estimating jointly, the method avoids the error accumulation seen in two-stage approaches, leading to tighter probability bounds.
3. Structured Sparsity: The use of sparse-group lasso penalties allows for flexible modeling where covariates can affect the global network structure or specific local edges.

4. Results

A. Simulation Studies

The authors compared cspine against RegGMM (Zhang and Li [2023]) under two scenarios:

1. Natural Model (Correctly Specified for cspine): cspine significantly outperformed RegGMM in True Positive Rate (TPR) for both population and covariate-dependent edges, and achieved lower estimation error ($\beta_{err}$).
2. Original Model (Correctly Specified for RegGMM): Even when the data was generated under the model assumed by RegGMM, cspine remained competitive, demonstrating robustness to model misspecification.
3. Sparsity Sensitivity: In simulations where the mean structure ($\Gamma$) became dense, RegGMM's performance degraded significantly due to first-stage estimation errors. cspine maintained stable performance, validating the advantage of joint convex optimization.

B. Real Data Application: Glioblastoma Multiforme (GBM)

The method was applied to re-analyze an eQTL dataset involving 73 genes and 118 SNPs from 401 GBM patients.

• Network Structure: cspine identified a sparser, more interpretable population network (476 edges) compared to RegGMM (750 edges).
• Covariate Effects: cspine identified 28 SNPs with significant effects on gene co-expression networks, whereas RegGMM identified only 12.
• Biological Insights:
  • rs10505980 (near KRAS): Both methods detected effects on the PI3K/MTOR pathway (PIK3CA-SHC2). cspine additionally identified an effect on PIK3CB-CAMK1G, suggesting a role in regulating oncogenic pathways.
  • rs1267622 (near BRAF): Both methods detected effects on the PDGF receptor pathway. cspine provided additional insights into GADD45A/B co-expression.

5. Significance and Conclusion

This paper provides a significant advancement in high-dimensional statistical learning for graphical models with covariates. By transforming a non-convex problem into a convex one via natural parametrization, the authors:

• Eliminate the need for error-prone two-stage estimation.
• Provide stronger theoretical guarantees regarding consistency and support recovery in high-dimensional settings.
• Offer a practical, scalable tool (cspine) that outperforms existing state-of-the-art methods in both simulations and real-world genomic applications.

The framework is particularly valuable for genomics and other fields where understanding how external factors (like genetics or environment) dynamically reshape interaction networks is critical. The method's ability to handle heterogeneous network effects with rigorous convex optimization makes it a robust choice for complex, high-dimensional data analysis.