Convex estimation of Gaussian graphical regression models with covariates

Gaussian graphical models (GGMs) are widely used to recover the conditional independence structure among random variables. Recent work has sought to incorporate auxiliary covariates to improve estimation, particularly in applications such as co-expression quantitative trait locus (eQTL) studies, where both gene expression levels and their conditional dependence structure may be influenced by genetic variants. Existing approaches to covariate-adjusted GGMs either restrict covariate effects to the mean structure or lead to nonconvex formulations when jointly estimating the mean and precision matrix. In this paper, we propose a convex framework that simultaneously estimates the covariate-adjusted mean and precision matrix via a natural parametrization of the multivariate Gaussian likelihood. The resulting formulation enables joint convex optimization and yields improved theoretical guarantees under high-dimensional scaling, where the sparsity and dimension of covariates grow with the sample size. We support our theoretical findings with numerical simulations and demonstrate the practical utility of the proposed method through a reanalysis of an eQTL study of glioblastoma multiforme (GBM), an aggressive form of brain cancer.

💡 Research Summary

This paper addresses the challenging problem of estimating Gaussian graphical models (GGMs) in which both the mean vector and the precision matrix depend on external covariates. In many modern applications, such as expression quantitative trait locus (eQTL) studies, genetic variants can influence not only gene expression levels (the mean) but also the conditional dependence structure among genes (the precision). Existing methodologies either restrict covariate effects to the mean, or, when allowing covariate‑dependent precision, lead to a non‑convex joint likelihood that forces practitioners to adopt alternating updates or two‑stage procedures. These approaches suffer from bias propagation, weakened theoretical guarantees, and computational instability.

The authors propose a novel convex framework, named Covariate‑adjusted Sparse Precision with Natural Estimation (cspine), that simultaneously estimates covariate‑adjusted means and precision matrices. The key insight is to re‑parameterize the Gaussian distribution in terms of its natural parameters: θ = diag(Ω)⁻¹Ωµ and Θ = ‑diag(Ω)⁻¹Ω. Under this transformation, the negative log‑likelihood becomes jointly convex in (θ, Θ). The covariate dependence is then modeled linearly as θ(u) = Γu and Θ(u) = B₀ + ∑_{h=1}^q B_h u_h, where Γ captures mean effects and each B_h encodes how covariate u_h modulates the network.

By substituting these expressions into the node‑wise regression representation of a Gaussian graphical model, the authors derive for each node j the regression equation
X_j = uᵀγ_j + ∑{k≠j}β{jk0}X_k + ∑{h=1}^q u_h∑{k≠j}β_{jkh}X_k + ε_j,
where γ_j corresponds to the mean coefficients, β_{jk0} to the baseline precision, and β_{jkh} to covariate‑specific adjustments. Crucially, this formulation eliminates the product terms that caused non‑convexity in earlier work, rendering the least‑squares loss jointly convex in all unknowns. Consequently, a single convex optimization problem can be solved for each node, without any alternating or staged procedures.

To induce interpretable sparsity, the authors employ a sparse‑group lasso penalty on the β coefficients: an ℓ₁ penalty promotes element‑wise sparsity, while a group ℓ₂,₁ penalty (applied across covariates for each node) encourages entire covariate effects to be either retained or discarded. An ℓ₁ penalty is also placed on γ_j to allow selective mean adjustment. This structured regularization aligns with biological expectations in eQTL settings, where a SNP may affect only a subset of gene‑gene interactions.

Theoretical analysis under high‑dimensional scaling (p, q ≫ n) demonstrates that the proposed estimator attains ℓ₁‑error bounds comparable to the oracle case where Γ is known, improving upon the rates of existing non‑convex or two‑stage methods. Because Θ and Ω share the same sparsity pattern, recovering Θ directly yields the desired graph structure. The authors show that the required sample size scales with √log p rather than log p, reflecting the advantage of the convex formulation.

Algorithmically, each nodewise problem is a standard convex least‑squares with mixed ℓ₁/ℓ₂,₁ penalties, solvable by proximal gradient methods (e.g., FISTA) or ADMM. The problems are embarrassingly parallel across nodes, leading to fast computation and guaranteed convergence.

Extensive simulations confirm that cspine outperforms competing methods in terms of edge‑selection accuracy, precision‑matrix estimation error, and recovery of covariate‑specific edge effects. Finally, the method is applied to a glioblastoma multiforme (GBM) eQTL dataset. The analysis uncovers SNPs that modulate specific gene‑gene conditional correlations, producing a biologically plausible network that aligns with prior findings while revealing novel regulatory relationships.

In summary, by leveraging a natural‑parameter re‑parameterization, the authors transform a fundamentally non‑convex covariate‑adjusted GGM estimation problem into a tractable convex one, delivering stronger theoretical guarantees, superior empirical performance, and practical utility for high‑dimensional biological data. Future directions include extensions to non‑linear covariate effects, dynamic networks, and Bayesian formulations.