missoNet: Convex Framework for mQTL Analysis

• missoNet is a convex, three-stage estimation framework that jointly models cis SNP effects and conditional CpG dependencies for regional mQTL analysis.
• It replaces unavailable empirical covariances with unbiased surrogate statistics under an MCAR missingness model, preserving convex optimization.
• The method employs sequential convex programs using Lasso, GLASSO, and FISTA, with strong theoretical guarantees and favorable empirical performance.

Searching arXiv for the specified paper and closely related context papers so the article can be grounded in current literature. missoNet is a convex, three-stage estimation framework for multivariate regression with missing response data that was introduced for sequencing-based regional DNA methylation QTL analysis. It jointly estimates a sparse regression coefficient matrix for cis SNP effects and a sparse precision matrix for residual dependence among CpG sites, thereby modeling both genetic regulation and conditional dependence within CpG blocks. Its central methodological feature is the use of unbiased surrogate statistics for masked responses, which avoids imputation while preserving convexity at each optimization stage and enabling non-asymptotic error guarantees under stated assumptions (Alam et al., 8 Jul 2025).

1. Statistical setting and inferential target

missoNet is formulated for regional mQTL analysis in which responses are transformed methylation levels at multiple CpGs within a genomic region and predictors are nearby cis SNP genotypes. The multivariate response formulation is motivated by the observation that treating a block of CpGs as a joint outcome can increase power and reveal coordinated regulation across CpGs. With $n$ samples, $p$ SNPs, and $q$ CpGs, the model is

$$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$

where $\mathbf{X}\in\mathbb{R}^{n\times p}$ has standardized columns, $\mathbf{Y}\in\mathbb{R}^{n\times q}$, and $\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$ satisfies $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$ and $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$. The residual precision matrix is $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$.

The two primary estimands are the sparse coefficient matrix $p$0, which encodes SNP-to-CpG effects, and the sparse precision matrix $p$1, whose zero off-diagonal entries encode conditional independences among CpGs under multivariate normality. This joint estimation problem places missoNet in direct continuity with multivariate regression and Gaussian graphical modeling, but specialized to the response-missingness regime that arises in bisulfite sequencing data.

A key feature of the setup is that missingness affects the responses rather than the predictors. The observed matrix is

$p$2

where $p$3 independently across rows, with column-specific missing rates $p$4. The framework therefore assumes missing completely at random (MCAR) for responses. This assumption is structurally important because the unbiased surrogate constructions and the associated theory are derived under precisely this masking model (Alam et al., 8 Jul 2025).

2. Surrogate covariance construction under response missingness

The central obstacle is that the empirical response covariance $p$5 and cross-covariance $p$6 are unavailable when $p$7 is partially masked. Pairwise deletion and imputation-based strategies can yield biased or indefinite covariance estimates, which is particularly problematic when the downstream objective includes a log-determinant term.

missoNet addresses this by replacing unavailable complete-data quantities with unbiased surrogates corrected by the first and second moments of the mask. The mask moments are

$p$8

and

$p$9

Using these quantities, missoNet defines

$q$0

Under MCAR, these are unbiased for $q$1 and $q$2. The predictor covariance remains the standard

$q$3

because predictors are fully observed.

These surrogate statistics are then substituted into a missing-data analogue of the MRCE objective. For complete data, MRCE minimizes a penalized Gaussian likelihood in $q$4. missoNet replaces $q$5 and $q$6 with $q$7 and $q$8, yielding the surrogate objective

$q$9

This objective is convex in either $$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$0 or $$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$1 when the other is fixed, provided the covariance input is positive semidefinite. The positive-semidefiniteness issue is therefore not peripheral but fundamental to the design of the full algorithm (Alam et al., 8 Jul 2025).

3. Three-stage convex estimation procedure

missoNet solves the surrogate problem through three sequential convex programs. The staged design separates sparse regression initialization, residual-network estimation, and coefficient refinement using the learned residual precision structure.

Stage	Optimization target	Computational strategy
I	Initial $$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$2 with $$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$3	Column-wise Lasso via coordinate descent or proximal gradient
II	$$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$4 from PSD-corrected residual covariance	PSD projection, then graphical Lasso with GLASSO or QUIC
III	Refined $$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$5 with fixed $$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$6	FISTA with Nesterov acceleration and backtracking

Stage I fixes $$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$7 and solves

$$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$8

This decomposes column-wise. For each response $$\mathbf{Y} = \mathbf{X}\mathbf{B}^{*} + \boldsymbol{\varepsilon},$$9,

$\mathbf{X}\in\mathbb{R}^{n\times p}$0

Stage II constructs a plug-in residual covariance

$\mathbf{X}\in\mathbb{R}^{n\times p}$1

Because $\mathbf{X}\in\mathbb{R}^{n\times p}$2 may be indefinite under missingness, directly solving a graphical Lasso objective can be non-convex or unbounded. missoNet therefore projects it to the nearest positive semidefinite matrix in max norm: $\mathbf{X}\in\mathbb{R}^{n\times p}$3 and then estimates the precision matrix by

$\mathbf{X}\in\mathbb{R}^{n\times p}$4

Stage III refines the coefficient matrix with $\mathbf{X}\in\mathbb{R}^{n\times p}$5 fixed: $\mathbf{X}\in\mathbb{R}^{n\times p}$6 The implementation uses FISTA with Nesterov acceleration and backtracking line search. The iterative updates are

$\mathbf{X}\in\mathbb{R}^{n\times p}$7

with soft-thresholding

$\mathbf{X}\in\mathbb{R}^{n\times p}$8

The Armijo condition is used for backtracking, and the convergence rate improves from $\mathbf{X}\in\mathbb{R}^{n\times p}$9 to $\mathbf{Y}\in\mathbb{R}^{n\times q}$0. This staged decomposition suggests that missoNet is designed less as a monolithic block-coordinate routine than as a sequence of convex subproblems with carefully controlled inputs (Alam et al., 8 Jul 2025).

4. Tuning, assumptions, and theoretical guarantees

missoNet supports both BIC-style tuning and cross-validation. The joint BIC for selecting $\mathbf{Y}\in\mathbb{R}^{n\times q}$1 is

$\mathbf{Y}\in\mathbb{R}^{n\times q}$2

The paper also specifies stage-specific BICs corresponding to the three estimation stages. For practical selection, a one-standard-error rule, denoted cv.1se, is recommended for parsimony and better false-discovery-rate control.

The theoretical analysis is non-asymptotic and quantifies the effects of missingness, dimensionality, sparsity, and noise. Stage I relies on a restricted eigenvalue condition for each response-specific support set $\mathbf{Y}\in\mathbb{R}^{n\times q}$3, with cone

$\mathbf{Y}\in\mathbb{R}^{n\times q}$4

and

$\mathbf{Y}\in\mathbb{R}^{n\times q}$5

Under the stated deviation bounds and tuning condition $\mathbf{Y}\in\mathbb{R}^{n\times q}$6, Proposition 1 yields

$\mathbf{Y}\in\mathbb{R}^{n\times q}$7

Stage II adds an irrepresentability condition for graph selection, expressed through the Kronecker matrix $\mathbf{Y}\in\mathbb{R}^{n\times q}$8. After PSD projection, the residual covariance obeys

$\mathbf{Y}\in\mathbb{R}^{n\times q}$9

With $\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$0, Proposition 2 gives

$\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$1

and

$\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$2

Stage III assumes a restricted eigenvalue condition for the refined matrix loss,

$\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$3

and defines

$\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$4

Theorem 1 then gives

$\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$5

The theory is explicitly tied to the assumptions that response missingness is MCAR, errors are sub-Gaussian, columns of $\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$6 are standardized, and the missingness rates $\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$7 are estimated from data to compute $\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$8 and $\boldsymbol{\varepsilon}\in\mathbb{R}^{n\times q}$9. This suggests that the guarantees are strong within the paper’s formal regime, but are not advertised as covering more general missing-data mechanisms (Alam et al., 8 Jul 2025).

5. Empirical behavior in simulations and CARTaGENE mQTL analysis

The simulation studies examine both scaling behavior and direct method comparison. Under AR(1) errors and MCAR masks, Stage I error $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$0 decreases at $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$1, increases with $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$2, and grows with missingness $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$3. Stage II precision error $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$4 scales like $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$5 and grows with missingness and response dimension $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$6. Stage III coefficient error $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$7 improves over Stage I uniformly and scales like $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$8 and $$\mathbb{E}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\mathbf{0}$$9, with clearer benefits when residual correlation $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$0 is high. For $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$1, at $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$2 and $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$3, the Stage I error is $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$4 versus Stage III $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$5; at $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$6, the corresponding values are $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$7 versus $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$8. The paper also notes that at $$\mathrm{Cov}(\boldsymbol{\varepsilon}\mid \mathbf{X})=\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*}$$9 missingness with small $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$0, Stage III may underperform Stage I, but improves with larger $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$1.

Across comparative experiments involving MRCE, Lasso, cglasso, and missoNet variants, missoNet [cv.min] consistently achieves the lowest prediction error for $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$2. In Model 1 with $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$3 missingness, the reported prediction errors are missoNet $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$4, MRCE $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$5, and Lasso $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$6. For network estimation, missoNet [BIC] often attains the lowest KLL; for example, in Model 2 with $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$7 missingness, the reported KLL is $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$8. For sparsity recovery of $$\boldsymbol{\Theta}_{\varepsilon\varepsilon}^{*}=(\boldsymbol{\Sigma}_{\varepsilon\varepsilon}^{*})^{-1}$$9, missoNet [cv.1se] yields high TPR, excellent TNR, and top MCC; in Model 1 with $p$00 missingness, the reported TNR is $p$01 and MCC is $p$02. Computation is also favorable: missoNet [BIC] is reported as fast due to convex programs, warm starts, local search, and parallelization, whereas MRCE is slow in all settings and cglasso is generally slower, becoming faster only in some high-dimensional instances.

The real-data analysis uses 202 adults from CARTaGENE with targeted bisulfite sequencing methylomes. Preprocessing regressed CpG methylation residuals on cohort, ACPA status, age, sex, and a cell composition proxy. Residuals were grouped into 62 candidate CpG regions, and nearby SNPs within $p$03 kb were considered after an initial univariate screen at $p$04. Training used 115 high-depth samples and independent validation used 87 earlier samples. Average missingness was approximately $p$05 in training and $p$06 in validation, with an additional $p$07 MCAR masking applied on complete CpGs for robustness.

On validation, missoNet [cv.1se] achieved the best mean squared prediction error:

• missoNet [cv.1se]: $p$08
• Lasso [cv.1se]: $p$09
• MRCE [cv.min]: $p$10
• cglasso [BIC]: $p$11

For out-of-sample log-likelihood, cglasso [BIC] had the highest value at $p$12 but the worst MSE, while missoNet [BIC] achieved $p$13 with over-dense networks and missoNet [cv.1se] achieved $p$14, which the study characterizes as balancing fit and parsimony. False discovery checks using 10 synthetic null SNPs per region showed that missoNet [cv.1se] attained

$p$15

with $p$16 and $p$17. By contrast, Lasso [cv.min] retained many dummy SNPs, with $p$18. Across 9,238 CpG-SNP pairs tested in 16 illustrative regions, missoNet [cv.1se] selected 957 nonzero associations, of which 841 ($p$19) matched known mQTLs from Gaunt et al. 2016 or involved CpGs within 100 bp of known mQTL sites. Some regions had no associations, including chr20:32,255,222–32,255,294, while others, including chr1:8,958,863–8,959,033, yielded novel candidates for follow-up (Alam et al., 8 Jul 2025).

6. Software implementation, operational guidance, and scope

missoNet is available as an R package on CRAN via install.packages("missoNet"). The package takes fully observed $p$20 with standardized columns and $p$21 with missing responses set to zero, together with a mask or inferred missingness probabilities $p$22. It internally computes $p$23, $p$24, and $p$25. The principal outputs are $p$26, $p$27, $p$28, and $p$29, along with fit diagnostics such as BIC grids, cross-validation errors, and plotting utilities including a BIC heatmap. Tuning can be specified by "cv" or "BIC", with options including nlambda, lambda.min.ratio, warm starts, fast=TRUE for local search, and parallel=TRUE for parallelization (Alam et al., 8 Jul 2025).

The package defines maximal regularization values as

$p$30

and then constructs logarithmic sequences to smaller values. Operationally, Stage I consists of $p$31 independent Lasso fits, Stage II uses graphical Lasso on the PSD-corrected residual covariance, and Stage III uses FISTA with backtracking and stopping criteria based on the Armijo condition and relative change thresholds. The paper emphasizes that warm starts, local search, and parallelization reduce runtime substantially, with missoNet [BIC] being the fastest setting and cross-validation trading additional cost for generalizability.

Several limitations are explicit. The missingness model is MCAR for responses; the theory assumes sub-Gaussian errors, restricted eigenvalue conditions, and irrepresentability for graph selection; and estimation of the column-specific missing rates $p$32 is required. The empirical studies further indicate that very high missing rates, such as at least $p$33, combined with small sample size can limit the gains of Stage III refinement. A plausible implication is that missoNet is best matched to moderate-to-high-dimensional genomic settings where missingness is substantial but not overwhelming, and where conditional dependence among responses is biologically meaningful.

The name has a separate usage in collider-physics literature: the DeepMET study refers to a “missoNet”-style per-particle weighting design for missing transverse momentum estimation, which is distinct from the genomic regression framework described here (Collaboration, 15 Sep 2025).