3.2 Cell-type covariance for the DeCovarT plug-in
2026-11-02
Source:vignettes/03-2-covariance-network-estimation.qmd
Notation (shared by all vignettes; aligned with the article)
Gene g=1,\ldots,G; cell type j=1,\ldots,J (a second type is \ell); bulk sample i=1,\ldots,N; reference donor (replicate) r=1,\ldots,R; cell c, with C_j reference cells of type j. The index i is dropped inside a derivation and restored in the model statements.
Symbol Meaning y_{gi}, \boldsymbol{y}_{\cdot i} Bulk expression on a linear scale m_{ji}, n_i=\sum_jm_{ji} Cells of type j in bulk sample i, and their total \tilde{\boldsymbol{y}}_{\cdot i}=\boldsymbol{y}_{\cdot i}/n_i Bulk expression per cell \boldsymbol{p}_{\cdot i}\in\Delta^{J-1} Cell fractions, p_{ji}=m_{ji}/n_i \boldsymbol{\rho}_{i}\in\mathbb{R}^{J-1} Isometric log-ratio (ILR) coordinates of \boldsymbol{p}_{\cdot i}, \boldsymbol{p}_{\cdot i}=\operatorname{softmax}(\mathbf{V}\boldsymbol{\rho}_i) with a Helmert basis \mathbf{V} \boldsymbol{q}_{\cdot i} Transcript fractions, q_{ji}=p_{ji}S_j/\sum_\ell p_{\ell i}S_\ell \boldsymbol{\mu}_{\cdot j}, \boldsymbol{\Sigma}_j Mean and G\times G gene–gene covariance of a type-j population profile \boldsymbol{\Omega}_j=\boldsymbol{\Sigma}_j^{-1} Precision matrix of type j \boldsymbol{\Gamma}_g J\times J covariance across cell types for gene g ( bMIND,EPIC-unmix,Unico)\boldsymbol{x}_{\cdot j,i} Purified profile of type j in sample i x_{gc}, \boldsymbol{x}_{\cdot c}, s_c Raw UMI count of gene g in reference cell c, and the size factor of cell c \boldsymbol{\lambda}_{\cdot c}, \boldsymbol{z}_{\cdot c} Size-free expression of cell c; its latent log-expression in the Poisson log-normal (PLN) model \tilde{\boldsymbol{\mu}}_{\cdot j}=C_j^{-1}\sum_{c\in j}\boldsymbol{x}_{\cdot c}=S_j\boldsymbol{\theta}_{\cdot j} Per-cell mean of the raw UMI counts of the C_j reference cells of type j \boldsymbol{\theta}_{\cdot j}, S_j=\mathbf{1}^{\top}\tilde{\boldsymbol{\mu}}_{\cdot j} Relative abundance (\mathbf{1}^{\top}\boldsymbol{\theta}_{\cdot j}=1) and mean transcriptome size of a type-j cell \boldsymbol{\Sigma}_{W,j}, \boldsymbol{\Sigma}_{B,j} Within-donor cell-to-cell covariance of type j (one single-cell sample); between-donor covariance of type j (needs R\ge2 donors) \boldsymbol{S}_W(\boldsymbol{v})=\sum_jv_j\boldsymbol{\Sigma}_{W,j}, \boldsymbol{S}_B(\boldsymbol{v})=\sum_jv_j^2\boldsymbol{\Sigma}_{B,j} Linear and quadratic covariance mixtures with weights \boldsymbol{v} \boldsymbol{V}_j Empirical covariance of type j, computed from its reference profiles or cells \boldsymbol{R}_j Empirical correlation matrix of type j (bold; the scalar R is the number of donors) \boldsymbol{\Delta}_j=\operatorname{diag}(\boldsymbol{V}_j) Diagonal matrix of the empirical variances of type j \lambda Penalty of a regularised covariance or precision estimator \psi_g, \mathcal{D}_g Negative binomial dispersion of gene g; deviance of gene g in a count model \delta_{gj} Detection rate of gene g in type j (fraction of cells with a non-zero count) \varpi_{\ell\mid j}, \varpi_{j\ell}, \bar\varpi MixSim overlap: probability of misclassifying a draw of j into \ell, pairwise overlap \varpi_{j\ell}=\varpi_{\ell\mid j}+\varpi_{j\mid\ell}, and its average over pairs ( BarOmega)\boldsymbol{T}_{ji}=\sum_{c=1}^{m_{ji}}\boldsymbol{X}_{jic} Total expression contributed by type j to sample i a_i Sample-wide scale between the reference and the bulk library c_i=a_in_i, \phi_i=a_i^2n_i Mean scale and variance scale of the single-cell model b_g, \tilde b_g=b_g^{\mathrm{bulk}}/b_g^{\mathrm{ref}} Capture efficiency of gene g in one protocol, and the gene-wise platform ratio; \tilde{\boldsymbol{B}}=\operatorname{diag}(\tilde{\boldsymbol{b}}) \boldsymbol{D}_i Diagonal, heteroscedastic bulk measurement noise A tilde marks a quantity rescaled to a reference unit: per cell (\tilde{\boldsymbol{y}}, \tilde{\boldsymbol{\mu}}) or relative to the reference protocol (\tilde{\boldsymbol{b}}). \boldsymbol{\Sigma}_j couples genes within one cell type; \boldsymbol{\Gamma}_g couples cell types within one gene. A symbol absent from this table is defined in the vignette that uses it, with one meaning there.
Scope
This vignette states which covariance DeCovarT needs from a single-cell reference, and how to estimate it. The target is the within-type covariance \boldsymbol{\Sigma}_{W,j} that enters the likelihood, not a biologically interpretable gene regulatory network. The default estimator is linear shrinkage. Three sparse-precision alternatives are derived in 3.3 Precision networks and scored by the same cross-validated criterion. The application to the Suppinger gastruloid reference (solver settings, outputs, plots) is in the covariance chapter of the reproducibility book.
Local symbols, defined once here: \omega_{gg'} is an entry of \boldsymbol{\Omega}_j; d is the maximum degree of a graph; K is the number of cross-validation folds; \alpha is an FDR level; \gamma is the EBIC parameter; \circ is the elementwise product and e^{\boldsymbol{A}} an elementwise exponential.
The target: a covariance, not a graph
DeCovarT evaluates a Gaussian likelihood whose covariance is a sum of cell-type covariances. With one single-cell specimen, the identified object is the within-type cell-to-cell covariance \boldsymbol{\Sigma}_{W,j}, which enters linearly in p_{ji} (one-replicate regime). With purified population profiles, \boldsymbol{\Sigma}_j enters through \sum_jp_{ji}^{2}\boldsymbol{\Sigma}_j (convolution). Either way the estimate must live on the linear count scale of the mean signature (normalisation), not on \log(1+x) or on nonparanormal scores.
Two properties make the covariance, not the precision, the natural object.
- Marginalisation. If the cells of type j have covariance \boldsymbol{\Sigma}_j over all genes, the covariance of a panel \mathcal{P} is the sub-matrix (\boldsymbol{\Sigma}_j)_{\mathcal{P}\mathcal{P}}, whatever other genes are measured. The precision of the panel is \{(\boldsymbol{\Sigma}_j)_{\mathcal{P}\mathcal{P}}\}^{-1}, a Schur complement of \boldsymbol{\Omega}_j, not the sub-matrix (\boldsymbol{\Omega}_j)_{\mathcal{P}\mathcal{P}}. Adding or removing a gene changes every partial correlation and every conditional variance 1/\omega_{gg}; it leaves each covariance entry unchanged. A covariance entry is a property of a gene pair; a graph is a property of the panel.
- The mixture has no sparse precision. Sums such as \sum_jp_{ji}\boldsymbol{\Sigma}_{W,j} or \sum_jp_{ji}^{2}\boldsymbol{\Sigma}_j have dense inverses even when every \boldsymbol{\Omega}_j is sparse. The likelihood never uses the zeros of \boldsymbol{\Omega}_j.
Sparsity of \boldsymbol{\Omega}_j therefore buys DeCovarT only a lower variance, through |E_j|+G free parameters instead of G(G+1)/2. Existence for C_j\le G does not require it: any regulariser provides that. Interpretation of the zeros and the recovery theory of Sec. 2 matter for a network, not for the deconvolution. The default plug-in is therefore a dense shrinkage covariance (Sec. 3.1). Sparse precisions are alternatives, judged on the same held-out score (Sec. 4). When the proportions are known, CSNet estimates type-level co-expression from the same convolution (CSNet) (Su et al. 2024).
| Estimator | Data it sees | Tuning parameter | Object exported to DeCovarT | Graph |
|---|---|---|---|---|
| Shrinkage, target D (Sec. 3.1) | linear counts | \lambda\in[0,1]; closed-form \hat\lambda^{\ast} as start | dense \widehat{\boldsymbol{\Sigma}}_j, variances kept | none |
| Graphical lasso, unpenalised diagonal (graphical lasso) | correlation \boldsymbol{R}_j of linear counts | \lambda | \boldsymbol{\Delta}_j^{1/2}\widehat{\boldsymbol{\Sigma}}_{R}\boldsymbol{\Delta}_j^{1/2} | support of \widehat{\boldsymbol{\Omega}}_j |
SILGGM after the nonparanormal map (nonparanormal) |
nonparanormal scores | FDR level \alpha | linear-scale refit under the zeros of the support | FDR-controlled support |
PLNnetwork (PLNnetwork) |
raw UMIs, offset \log s_c | \lambda | linear moments of the latent (\boldsymbol{\mu}_z,\boldsymbol{\Sigma}_z) (linear-moment map) | support of the latent precision |
Sample-size regimes
Report C_j against G for every type (and time point) before reading any estimate. In the Suppinger reference several types have C_j<G, so the sample covariance \boldsymbol{V}_j has rank at most C_j-1 and no unpenalised estimator exists for them. Three thresholds are often conflated.
The penalised estimate exists. With an off-diagonal penalty, the graphical lasso has a unique positive definite solution for any \lambda>0 once every V_{gg}>0, even when C_j<G (Carter 2025). Shrinkage towards a positive definite target exists for any C_j\ge2. No sparsity of the truth is needed.
-
The maximum-likelihood estimate on a fixed graph exists. Dempster’s refit under the zeros of a graph E,
\widehat{\boldsymbol{\Omega}}_{E} =\arg\max_{\boldsymbol{\Omega}\succ0} \log\det\boldsymbol{\Omega}-\operatorname{tr}(\boldsymbol{V}_j\boldsymbol{\Omega}) \quad\text{s.t.}\quad \omega_{gg'}=0\ \ \forall\,(g,g')\notin E, \tag{1}
exists for almost every sample once C_j reaches the maximum-likelihood threshold \operatorname{mlt}(E). That threshold lies between the clique number and the treewidth plus one, and equals the clique number for chordal graphs (Uhler 2012; Bernstein & Outlaw 2023). Bernstein and Outlaw simulate graphs selected by the graphical lasso from independent Gaussian data (G\le9): the probability that C_j\ge\operatorname{mlt} of the selected graph rises with \lambda and with C_j; with a single observation the graphical lasso never returned the empty graph, the only graph whose estimate then exists (pp. 4–5).
The graph is recovered. Support recovery by the graphical lasso needs C_j\gtrsim d^{2}\log G, an irrepresentability condition, bounded eigenvalues and a minimum |\omega_{gg'}| (Ravikumar et al. 2011). Neighbourhood selection needs d\log G (Chiquet 2015, 41), under its own neighbourhood stability condition (Meinshausen & Bühlmann 2006, assumption 6). With G=500 and d=10, d^{2}\log G\approx620 before constants, and hubs raise d. On models derived from gene expression the irrepresentability conditions never held (Heinävaara et al. 2016), and they involve the unknown \boldsymbol{\Sigma}_j: “we can never verify them on real data” (Chiquet 2015, 77).
The consequence for DeCovarT: threshold 1 is all the plug-in needs. Threshold 2 governs the refits used by SILGGM (nonparanormal). Threshold 3 governs any biological reading of the zeros, which is out of scope here; it can be probed by simulating precisions at the observed C_j with simulate_hierarchical_grn_moments() and measuring recovery.
Estimators
Linear shrinkage (default)
Schäfer and Strimmer shrink the sample covariance towards a target \boldsymbol{T}, \widehat{\boldsymbol{\Sigma}}=\lambda\boldsymbol{T}+(1-\lambda)\boldsymbol{V}_j (Schäfer and Strimmer 2005, eq. 3). Their target D (“diagonal, unequal variance”) keeps the variances and shrinks the correlations towards zero (Table 1, eq. 10):
\hat\Sigma_{gg}=V_{gg}, \qquad \hat\Sigma_{gg'}=(1-\lambda)\,R_{j,gg'}\sqrt{V_{gg}V_{g'g'}}, \qquad \hat\lambda^{\ast} =\frac{\sum_{g\neq g'}\widehat{\operatorname{Var}}(R_{j,gg'})} {\sum_{g\neq g'}R_{j,gg'}^{2}}, \tag{2}
with \hat\lambda^{\ast} truncated to [0,1]. The closed form minimises the Frobenius risk (Ledoit–Wolf lemma; Ledoit & Wolf 2004). The estimate is positive definite for any C_j and G, invariant to the scale of each gene, and costs about two sample covariances. Its precision is dense; partial correlations and a local-FDR graph can be read from it if wanted (eqs. 11–15).
The Frobenius risk is not the Gaussian likelihood that DeCovarT evaluates. \hat\lambda^{\ast} is therefore only the starting point of the cross-validated search of Sec. 4. In corpcor::cov.shrink(), lambda.var = 0 reproduces target D; the default also shrinks the variances towards their median.
Sparse precisions are not the default plug-in. The graphical lasso, the nonparanormal SILGGM estimator and PLNnetwork are derived in 3.3 Precision networks. Each still exports a covariance on the linear scale, and each is scored by the same held-out criterion as shrinkage (Sec. 4).
One criterion for all four: cross-validated Gaussian log-score
Each estimator of Table 1 has one tuning parameter: \lambda for shrinkage (searched around \hat\lambda^{\ast}), \lambda for the graphical lasso, \alpha for SILGGM, \lambda for PLNnetwork. All four are tuned by the same score.
Score. Split the C_j cells of type j into K folds. Fit (\widehat{\boldsymbol{\mu}}^{(-k)},\widehat{\boldsymbol{\Sigma}}^{(-k)}) on the other folds, exported on the linear scale, and score the C_{j,k} held-out cells by their mean Gaussian log-density,
\ell^{(k)} =-\frac12\Bigl[\log\det\widehat{\boldsymbol{\Sigma}}^{(-k)} +\operatorname{tr}\bigl\{(\widehat{\boldsymbol{\Sigma}}^{(-k)})^{-1}\boldsymbol{V}_{j,k}\bigr\} +G\log2\pi\Bigr], \tag{3}
with \boldsymbol{V}_{j,k}=C_{j,k}^{-1}\sum_{c\in k}(\boldsymbol{x}_{\cdot c}-\widehat{\boldsymbol{\mu}}^{(-k)})(\boldsymbol{x}_{\cdot c}-\widehat{\boldsymbol{\mu}}^{(-k)})^{\top}. The tuning value maximises K^{-1}\sum_k\ell^{(k)} (or the one standard-error rule). Minus the cross-validated score estimates the Kullback–Leibler risk of the fitted law up to a constant (Commenges and Jacqmin-Gadda 2015, sect. 2.4.2). It applies to maximum-likelihood, penalised and shrinkage fits alike, whereas AIC and BIC are justified for the maximum-likelihood estimate. The folds split cells of one specimen. The score therefore measures prediction of new cells of that specimen, which is the population of \boldsymbol{\Sigma}_{W,j}, and says nothing about new specimens (Zimmerman et al. 2021).
Why a Gaussian score is legitimate for non-Gaussian cells. For a test cell with mean \boldsymbol{\mu}^{\ast} and covariance \boldsymbol{\Sigma}^{\ast}, whatever its law,
\mathbb{E}\bigl\{-\log\phi_G(\boldsymbol{x};\boldsymbol{\mu},\boldsymbol{\Sigma})\bigr\} =\tfrac12\bigl\{\log\det\boldsymbol{\Sigma} +\operatorname{tr}(\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}^{\ast}) +(\boldsymbol{\mu}-\boldsymbol{\mu}^{\ast})^{\top}\boldsymbol{\Sigma}^{-1}(\boldsymbol{\mu}-\boldsymbol{\mu}^{\ast}) +G\log2\pi\bigr\}. \tag{4}
Only the first two moments enter. The excess over its value at (\boldsymbol{\mu}^{\ast},\boldsymbol{\Sigma}^{\ast}) is, at \boldsymbol{\mu}=\boldsymbol{\mu}^{\ast}, half of Stein’s loss \operatorname{tr}(\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}^{\ast})-\log\det(\boldsymbol{\Sigma}^{-1}\boldsymbol{\Sigma}^{\ast})-G\ge0, which is zero only at \boldsymbol{\Sigma}=\boldsymbol{\Sigma}^{\ast}. The Gaussian log-score is thus a proper scoring rule for the mean and covariance, which is all that DeCovarT’s bulk likelihood uses: the bulk law is Gaussian by the central limit theorem over many cells, not because single cells are Gaussian.
Comparing estimators: nested cross-validation. Tuning needs one loop. Comparing the four tuned estimators needs two: an inner loop that chooses each tuning parameter on the training part of an outer fold, and an outer loop that scores the tuned estimators on cells the tuning never saw. Reporting the inner-loop maximum as the estimator’s performance is optimistic, because it is the maximum over a grid. Compare estimators by their paired outer-fold differences, with the diagonal covariance \boldsymbol{\Delta}_j as the baseline.
LOOCV is not nested cross-validation. Leave-one-out is the choice K=C_j in a single loop. It can serve as the inner or the outer loop of a nested scheme but is not one, and it costs C_j fits per grid value.
Fast approximation. The universal approximate cross-validation criterion approximates leave-one-out from one fit,
\mathrm{UACV}=\Psi(\hat{\boldsymbol{\theta}})+\operatorname{tr}\bigl(\boldsymbol{H}_{\Phi}^{-1}\boldsymbol{U}\bigr), \tag{5}
where \Psi is the mean assessment loss (here minus the log-score), \boldsymbol{H}_{\Phi} the Hessian of the estimating function and \boldsymbol{U} (their \boldsymbol{K}) the mean cross-product of the per-cell assessment and estimating scores; the error is O_p(C_j^{-2}) (Commenges et al. 2012, eqs. 6–7). For a maximum-likelihood fit assessed by its own log-likelihood it is the Takeuchi information criterion divided by 2C_j (p. 7). A smooth penalty enters \boldsymbol{H}_{\Phi} and lowers the correction (p. 8). An \ell_1 penalty is not differentiable at zero, so for the graphical lasso and PLNnetwork the formula is applied on the active set, with the support held fixed; that restriction is our reading, not the paper’s.
Diagnostics, not selectors. EBIC (Foygel and Drton 2010) charges \log C_j+4\gamma\log G per edge and relies on the maximum-likelihood degrees of freedom; at small C_j it returns empty graphs (screening threshold). StARS (Liu et al. 2010) chooses the least regularised graph whose edge instability across subsamples stays below a threshold. It costs B\times L fits for B subsamples and L penalties, targets the support rather than the covariance, and tends to select too many edges in Gaussian graphical models. Thresholded adaptive validation targets the support as well (Laszkiewicz et al. 2021). All three are reported beside the cross-validated choice, not instead of it.
The diagonal: conditional variance, not self-regulation
In a Gaussian graphical model, \omega_{gg}=1/\operatorname{Var}(x_g\mid\boldsymbol{x}_{-g}): the part of gene g’s variance that the other panel genes do not explain. On counts it mixes transcriptional bursting and auto-regulation, regulators outside the panel, and Poisson sampling, and \boldsymbol{V}_j alone cannot separate them; at droplet depths sampling is the largest share for most genes (Sarkar and Stephens 2021). By marginalisation (Sec. 1) \omega_{gg} also changes with the panel.
Penalising the diagonal biases the variances. With a penalised diagonal the KKT conditions give \hat\Sigma_{gg}=V_{gg}+\lambda (Friedman et al. 2008; Mazumder and Hastie 2011). On the correlation scale every implied correlation shrinks by about 1/(1+\lambda), in the same direction as Poisson attenuation (count attenuation). The penalty does not move the empty-graph threshold of screening threshold. DeCovarT’s likelihood uses the variances directly, so this bias enters the deconvolution.
Self-regulation belongs to directed dynamic models. CardamomOT estimates a directed self-interaction \theta_{ii} for each gene and gives it no specific penalty beyond its common \ell_1 term towards a prior network; a positive auto-activation with a negative basal term is “the most parsimonious way to produce bistability for an isolated gene” (Maugé & Ventre 2026, p. 28, citing Ferrell 2002). Such a parameter is linked to a precision only through the stationary covariance of a linearised model. For \mathrm{d}\boldsymbol{x}=\boldsymbol{A}(\boldsymbol{x}-\boldsymbol{\mu})\,\mathrm{d}t+\mathrm{d}\boldsymbol{w} with noise covariance \boldsymbol{Q}, the stationary covariance solves the Lyapunov equation \boldsymbol{A}\boldsymbol{\Sigma}+\boldsymbol{\Sigma}\boldsymbol{A}^{\top}+\boldsymbol{Q}=\boldsymbol{0}. For two genes with self-loops -a_1,-a_2, an edge 1\to2 of weight k and \boldsymbol{Q}=\operatorname{diag}(q_1,q_2), writing \kappa=(a_1+a_2)^{2}q_2+k^{2}q_1,
\omega_{11}=\frac{2(a_1+a_2)\{a_1(a_1+a_2)q_2+k^{2}q_1\}}{q_1\,\kappa}, \qquad \omega_{22}=\frac{2a_2(a_1+a_2)^{2}}{\kappa}, \qquad \omega_{12}=-\frac{2a_2(a_1+a_2)k}{\kappa}. \tag{6}
Each diagonal entry depends on both self-loops, the edge and both noise levels; only for k=0 does it reduce to 2a_g/q_g, which still mixes the self-loop with the noise. For a three-gene chain 1\to2\to3 the stationary \omega_{13} is non-zero for generic parameters although genes 1 and 3 share no edge. A GGM diagonal is therefore neither a self-loop nor its sign. Strong positive self-loops also produce bistable, hence bimodal, within-type laws, which a Gaussian \boldsymbol{\Sigma}_{W,j} cannot represent; the DeCovarT counterpart is to split the type into subtypes.
Decision
Never penalise the diagonal of \boldsymbol{\Omega}_j:
glassowithpenalize.diagonal = FALSE,PLNnetworkwith an unpenalised diagonal, shrinkage target D withlambda.var = 0. Schemes that re-estimate the diagonal after the off-diagonal step are compatible with this rule: SPACE estimates conditional variances from residuals (Peng et al. 2009); the partial-correlation graphical lasso penalises only the partial correlations and is scale invariant, computed by alternating the two blocks (Carter and Molinari 2025); the simplest version is a refit without penalty on a selected support (Eq. 1).
A near-diagonal precision is not, by itself, evidence that genes are independent. The five competing explanations, the screening threshold and the count attenuation are in 3.3 Precision networks. Directed Gaussian networks and joint estimation across types are in the same vignette. Graph scores (Hamming distance, structural Hamming distance, BEELINE, stability selection used as a graph rule) do not choose \widehat{\boldsymbol{\Sigma}}_j. The plug-in is the covariance with the best held-out Gaussian log-score (Sec. 4). Those score definitions stay in the reproducibility file docs/05-grn-metrics.qmd, which is not a listed chapter.