Package 'DeCovarT'

Additive logistic transform (unconstrained coordinates to the simplex)

Description

Inverse additive log-ratio map $\boldsymbol{\psi}:\boldsymbol{\rho}\mapsto\boldsymbol{p}$ of Aitchison ($\mathrm{alr}^{-1}$): a softmax with the last category $J$ pinned as reference ($\rho_J\equiv 0$). Solvers and vcov.decovart_fit() use the isometric log-ratio chart isometric_logistic() instead. This ALR helper is retained for the vignette appendix.

Recovers the unconstrained additive log-ratio coordinates $\rho_j=\ln(p_j/p_J)$ for $j=1,\ldots,J-1$, with the last part $p_J$ as reference. This is Aitchison's additive log-ratio ($\mathrm{alr}$) transform, equivalently the multinomial-logit link with reference category $J$ (see compositions::alr() and vignette("theory-decovart-generative-model", package = "DeCovarT")).

Usage

additive_logistic(rho)

additive_log_ratio(p)

Arguments

rho	Numeric vector $\boldsymbol{\rho}\in\mathbb{R}^{J-1}$ of unconstrained additive log-ratio coordinates (reference cell type $J$).
p	Numeric vector $\boldsymbol{p}$ on the open simplex.

Details

With $A=\sum_{k=1}^{J-1}\mathrm{e}^{\rho_k}+1$,

$p_j=\frac{\mathrm{e}^{\rho_j}}{A}\quad(j<J),\qquad p_J=\frac{1}{A}.$

Equivalently, $\boldsymbol{p}=\boldsymbol{\psi}(\boldsymbol{\rho})$ with $\boldsymbol{\psi}(\boldsymbol{\rho})\propto (\mathrm{e}^{\rho_1},\ldots,\mathrm{e}^{\rho_{J-1}},1)^{\mathsf{T}}$. Jacobians and Hessians of both maps are derived in the package vignette vignette("theory-decovart-generative-model", package = "DeCovarT"). See also compositions::alrInv().

Value

Numeric vector $\boldsymbol{p}\in\mathbb{R}^{J}$ on the unit simplex ($\mathbf{1}^{\mathsf{T}}\boldsymbol{p}=1$, $\boldsymbol{p}\ge\mathbf{0}$).

Numeric vector $\boldsymbol{\rho}\in\mathbb{R}^{J-1}$.

Numerical stability

Evaluated through the log-sum-exp shift $p_i\propto\exp(\tilde{\rho}_i-\max_k\tilde{\rho}_k)$ with $\tilde{\boldsymbol{\rho}}=(\rho_1,\ldots,\rho_{J-1},0)$, which is algebraically identical to the ratio above but never forms $\mathrm{e}^{\rho_i}$ directly. The naive quotient overflows to NaN for $\rho_i\gtrsim 710$, which unconstrained ascent can reach when the MLE approaches a simplex face; the shifted form returns an exact zero instead and keeps $\boldsymbol{\Sigma}(\boldsymbol{p})$ positive definite.

See Also

The inverse map (additive log-ratio) is documented as additive_log_ratio() on this help page.

Examples

rho <- c(0.2, -0.5)
p <- additive_logistic(rho)
sum(p)
additive_log_ratio(p)
# Stable far out in ALR space, where exp(rho) would overflow.
additive_logistic(c(800, 1200))

---

Algorithm-similarity correlation from a Monte Carlo benchmark

Description

Pearson correlations $r_{ab}=\mathrm{cor}(\hat p_a,\hat p_b)$ across Monte Carlo replicates (and cell types). This is behavioural similarity: two solvers can correlate near 1 while remaining systematically biased. Hierarchical order uses Ward D2 on $d_{ab} = 1 - r_{ab}$.

Usage

algorithm_similarity(benchmark, facet_rows = NULL, facet_cols = NULL)

Arguments

benchmark	List from run_simulation_benchmark().
facet_rows, facet_cols	Optional scenario columns. When supplied, correlations are computed inside each facet cell.

Value

A tibble with algorithm_x, algorithm_y, correlation, and any facet columns.

See Also

plot_algorithm_similarity(), pivot_mc_estimates()

Examples

set.seed(1)
genes <- paste0("g", 1:2)
cts <- paste0("ct", 1:2)
mu <- matrix(c(20, 22, 22, 20), nrow = 2, dimnames = list(genes, cts))
Sigma <- array(
  c(1, 0, 0, 1, 1, 0, 0, 1),
  dim = c(2, 2, 2),
  dimnames = list(genes, genes, cts)
)
out <- run_simulation_benchmark(
  tibble::tibble(
    true_theta = list(list(p = c(0.5, 0.5), mu = mu, sigma = Sigma))
  ),
  deconvolution_functions = list(
    "nnls" = list(FUN = deconvolute_ratios_nnls),
    "rlm" = list(FUN = deconvolute_ratios_rlm)
  ),
  n = 4,
  cores = 1
)
algorithm_similarity(out)

---

Parametric bootstrap for DeCovarT proportions and boundary tests

Description

The conditional law of the bulk profile is fully specified, so the asymptotic calibrations of lrt_decovart() can be replaced by simulation. Given plug-in moments and a composition $\boldsymbol{p}$, the function draws $\boldsymbol{Y}^{(b)}\sim \mathcal{N}_{G}(\boldsymbol{\mu}\boldsymbol{p}, \boldsymbol{\Sigma}(\boldsymbol{p}))$, refits DeCovarT, and returns the empirical distribution of $\hat{\boldsymbol{p}}^{(b)}$ together with percentile intervals.

Usage

bootstrap_decovart(
  y = NULL,
  mean_signature_matrix,
  Sigma,
  bulk_expression = NULL,
  p = NULL,
  null_value = NULL,
  n_boot = 500L,
  level = 0.95,
  n_replicates = NULL,
  epsilon = 10^-8,
  itmax = 500L
)

Arguments

y	Numeric vector (or one-column matrix) $\boldsymbol{y}\in\mathbb{R}^{G}$.
mean_signature_matrix	Numeric matrix $\boldsymbol{\mu}\in\mathcal{M}_{G\times J}$ (plug-in means).
Sigma	Array of cell-type covariances in $\mathcal{M}_{G\times G\times J}$.
bulk_expression	Numeric matrix $\boldsymbol{Y}\in \mathcal{M}_{G\times N}$ of replicate bulk profiles sharing one composition, used instead of y for the pooled (population) test.
p	Composition to simulate from. Defaults to the MLE computed from y / bulk_expression.
null_value	Named numeric vector of hypothesised ratios, e.g. c(celltype_2 = 0).
n_boot	Number of bootstrap replicates.
level	Confidence level for percentile intervals.
n_replicates	Number of bulk columns simulated per replicate; defaults to the number of observed columns.
epsilon, itmax	Relative convergence tolerance and iteration budget passed to stats::optim().

Details

Supplying null_value switches to a restricted bootstrap: replicates are generated under the restricted MLE and the likelihood-ratio statistic is recomputed on each, giving a Monte Carlo p-value that does not rely on the $50{:}50$ asymptotic weights. For a composite null the remaining nuisance proportions are profiled at their restricted estimate, so the procedure is a plug-in (approximate) rather than an exact bootstrap.

Reference uncertainty is not propagated here: $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}_j$ are treated as known, exactly as in the frequentist likelihood. To resample the purified observations that produced those plug-in moments, use reference_bootstrap_decovart().

Value

A list with estimates ($J\times$ n_boot matrix), interval (percentile bounds), p_simulated, and, when null_value is supplied, statistic (observed $D$), null_statistics, and p_value.

See Also

lrt_decovart(), chi_bar_square_pvalue(), reference_bootstrap_decovart()

Other decovart_inference: boundary_diagnostics(), chi_bar_square_pvalue(), confint_profile_decovart(), equivariance_check_decovart(), lrt_decovart(), multistart_decovart(), profile_loglik_decovart(), reference_bootstrap_decovart(), restricted_mle_decovart()

Examples

mu <- matrix(c(20, 40, 15, 40, 20, 25), nrow = 3)
colnames(mu) <- paste0("ct", 1:2)
Sigma <- array(c(diag(3), diag(3)), dim = c(3, 3, 2))
y <- drop(mu %*% c(0.6, 0.4))
set.seed(1)
bootstrap_decovart(y, mu, Sigma, n_boot = 20)$interval

---

Validate generative-model parameters $\theta$

Description

Checks that true_theta stores the DeCovarT mixture parameters:

• p: mixing proportions of length $J$, or a matrix $J\times N$ of sample-wise proportions (each column on the simplex);

• mu: mean signature $\boldsymbol{\mu}\in\mathcal{M}_{G\times J}$ (also accepts $J\times G$);

• sigma and/or Theta: array $G\times G\times J$ of covariances and/or precision matrices.

Usage

check_true_theta(
  true_theta,
  require_p = TRUE,
  J = NULL,
  second_moment = c("either", "sigma", "Theta")
)

Arguments

true_theta	Named list with the fields above.
require_p	Logical; if TRUE (default), p must be present and valid. If FALSE, a missing p is allowed (callers may default it).
J	Optional expected number of cell types. When NULL, inferred from mu / the third dimension of sigma or Theta.
second_moment	Which second-moment field to require: "sigma", "Theta", or "either" (default: at least one of sigma / Theta). Matching is case-insensitive.

Value

TRUE invisibly if the structure is valid; otherwise stops with an informative error.

See Also

compute_average_overlap(), compute_average_jeffreys()

Examples

theta <- list(
  p = c(0.5, 0.5),
  mu = cbind(c(0, 0), c(3, 0)),
  sigma = array(c(diag(2), diag(2)), dim = c(2, 2, 2))
)
check_true_theta(theta)

---

Chi-bar-square tail probability for boundary likelihood-ratio tests

Description

Upper-tail probability of the mixture

$\bar{\chi}^{2} = \sum_{i=0}^{q} w_i\,\chi^{2}_{s+i}, \qquad w_i=\binom{q}{i}2^{-q},$

with $q$ constraints active on the boundary of the parameter space and $s$ interior constraints. This is Case 9 of (Self and Liang 1987), which generalises the $\tfrac{1}{2}\chi^{2}_{0}+\tfrac{1}{2}\chi^{2}_{1}$ result of (Chernoff 1954). For $q=0$ the mixture collapses to the ordinary $\chi^{2}_{s}$ of (Wilks 1938).

Usage

chi_bar_square_pvalue(
  statistic,
  n_boundary = 1L,
  df_interior = 0L,
  weights = NULL
)

Arguments

statistic	Numeric vector of observed likelihood-ratio statistics $D=2(\hat{\ell}_1-\hat{\ell}_0)$.
n_boundary	Integer $q$: number of constraints active on the boundary (e.g. null ratios set to zero).
df_interior	Integer $s$: number of interior (two-sided) constraints tested simultaneously.
weights	Optional numeric vector of length n_boundary + 1 giving the mixing weights $w_0,\ldots,w_q$. Defaults to the binomial weights.

Details

The binomial weights hold when the block of the Fisher information corresponding to the active constraints is diagonal after orthogonalisation against the nuisance block. Under correlated constraints the mixing weights depend on the geometry of the approximating tangent cone and the binomial choice is only an approximation; Self and Liang also exhibit a nuisance-parameter configuration whose limit is not a chi-square mixture. Supply weights explicitly, or calibrate by simulation with bootstrap_decovart(), whenever $q>1$ and the constrained cell types share signal.

Value

Numeric vector of upper-tail probabilities.

Atom at zero

$\chi^{2}_{0}$ is a point mass at the origin, so its strict upper tail $\Pr(\chi^{2}_{0}>D)$ is zero for every $D\ge 0$. The component therefore only removes weight: with one active constraint and no interior constraint the p-value is $\tfrac{1}{2}\Pr(\chi^{2}_{1}>D)$, capped at $0.5$ when $D=0$. Note stats::pchisq() switches to the closed tail at exactly q = 0 for df = 0; the strict convention is used here so the p-value stays continuous in statistic.

References

Self SG, Liang K (1987). “Asymptotic Properties of Maximum Likelihood Estimators and Likelihood Ratio Tests under Nonstandard Conditions.” Journal of the American Statistical Association, 82(398), 605–610. ISSN 0162-1459. doi:10.1080/01621459.1987.10478472.

Chernoff H (1954). “On the Distribution of the Likelihood Ratio.” The Annals of Mathematical Statistics, 25(3), 573–578. ISSN 0003-4851. doi:10.1214/aoms/1177728725.

Wilks SS (1938). “The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses.” The Annals of Mathematical Statistics, 9(1), 60–62. ISSN 0003-4851. doi:10.1214/aoms/1177732360.

See Also

lrt_decovart()

Other decovart_inference: bootstrap_decovart(), boundary_diagnostics(), confint_profile_decovart(), equivariance_check_decovart(), lrt_decovart(), multistart_decovart(), profile_loglik_decovart(), reference_bootstrap_decovart(), restricted_mle_decovart()

Examples

# One active boundary constraint: 50:50 mixture of chi2_0 and chi2_1.
chi_bar_square_pvalue(2.71, n_boundary = 1L)
# Interior (Wilks) calibration.
chi_bar_square_pvalue(3.84, n_boundary = 0L, df_interior = 1L)

---

One-dominant composition with a target normalised Shannon entropy

Description

Returns a length-$J$ simplex vector of the form $(1-(J-1)q,\,q,\ldots,q)$ whose normalised Shannon entropy compute_shannon_entropy() equals h_star. The uniform composition (h_star = 1) and a Dirac mass (h_star = 0) are returned exactly.

Usage

composition_from_entropy(h_star, n_celltypes = 3L, nms = NULL)

Arguments

h_star	Target $H^{\star}\in[0,1]$.
n_celltypes	Integer $J\ge 2$.
nms	Optional names for the returned vector.

Value

Numeric simplex vector of length $J$.

See Also

compute_shannon_entropy()

Examples

p <- composition_from_entropy(0.5, n_celltypes = 3L)
compute_shannon_entropy(p)

---

Average pairwise Jeffreys divergence of a Gaussian mixture

Description

Returns the $p_j p_{\ell}$-weighted mean of pairwise Jeffreys (symmetrised KL) divergences between purified Gaussians

$\overline{J} = \frac{ \sum_{1\le j<\ell\le J} p_j p_{\ell}\, J(f_{j},f_{\ell}) }{ \sum_{1\le j<\ell\le J} p_j p_{\ell} } \in[0,\infty),$

with $f_{j}=\mathcal{N}_{G}( \boldsymbol{\mu}_{\cdot j},\boldsymbol{\Sigma}_{j} )$. If p is omitted it defaults to the equi-balanced vector $1/J$, which recovers the uniform pairwise average $2/(J(J-1))\sum_{j<\ell}J(f_{j},f_{\ell})$.

Usage

compute_average_jeffreys(true_theta, J = NULL)

Arguments

true_theta	List validated by check_true_theta(): mu ($G\times J$), sigma ($G\times G\times J$), and optionally p (length $J$ or $J\times N$). If p is missing it is set to $(1/J,\ldots,1/J)$.
J	Number of cell types. Defaults to the third dimension of sigma.

Value

Scalar average pairwise Jeffreys divergence.

See Also

check_true_theta(), compute_average_overlap(), compute_glmnet_gene_scores()

Examples

set.seed(1)
theta <- list(
  mu = cbind(c(0, 0), c(3, 0)),
  sigma = array(c(diag(2), diag(2)), dim = c(2, 2, 2))
)
compute_average_jeffreys(theta)

---

Average pairwise overlap of a Gaussian mixture

Description

Returns MixSim's BarOmega: the unweighted mean of pairwise overlaps

$\overline{\omega} = \frac{2}{J(J-1)} \sum_{1\le j<\ell\le J} \bigl(\Omega_{j\ell}+\Omega_{\ell j}\bigr) \in[0,1]$

(up to the MixSim numerical convention), where $\Omega_{j\ell}=\Pr_{X\sim f_j}(X\text{ classified as }\ell)$ already uses the mixture weights $\boldsymbol{p}$ inside the Bayes / MAP rule of MixSim::overlap(). Do not multiply the directional masses by $p_j$ again.

For $G\ge 4$ MixSim's Davies quadrature becomes expensive. The helper then switches to overlap_gaussian_mc(): stratified Sobol draws, inverse-transform sampling through precomputed Cholesky factors, and log-density MAP comparisons (n_mc draws per component; default 10,000).

Usage

compute_average_overlap(
  true_theta,
  J = NULL,
  n_mc = 10000L,
  seed = NULL,
  verbose = TRUE
)

Arguments

true_theta	List validated by check_true_theta(): p (length $J$ or $J\times N$), mu ($G\times J$), sigma ($G\times G\times J$).
J	Number of cell types (components). Defaults to the third dimension of sigma.
n_mc	Monte Carlo draws per component when $G\ge 4$.
seed	Optional seed forwarded to overlap_gaussian_mc().
verbose	If TRUE (default), announce the MixSim-to-MC switch with cli when it is installed.

Value

Scalar average pairwise overlap (MixSim BarOmega).

See Also

check_true_theta(), overlap_gaussian_mc()

Examples

set.seed(1)
theta <- list(
  p = c(0.5, 0.5),
  mu = cbind(c(0, 0), c(3, 0)),
  sigma = array(c(diag(2), diag(2)), dim = c(2, 2, 2))
)
compute_average_overlap(theta)

---

Mean pairwise AIRM distance of component covariances

Description

Mean pairwise AIRM distance of component covariances

Usage

compute_average_riemannian(true_theta, J = NULL)

Arguments

true_theta	List validated by check_true_theta(): mu ($G\times J$), sigma ($G\times G\times J$), and optionally p (length $J$ or $J\times N$). If p is missing it is set to $(1/J,\ldots,1/J)$.
J	Number of cell types. Defaults to the third dimension of sigma.

Value

Scalar average of $d_R(\Sigma_j,\Sigma_\ell)$ over $j<\ell$.

See Also

spd_affine_invariant_distance()

Examples

a <- diag(2)
b <- matrix(c(2, 0.5, 0.5, 1), nrow = 2)
theta <- list(
  mu = cbind(c(0, 0), c(1, 0)),
  sigma = array(c(a, b), dim = c(2, 2, 2))
)
compute_average_riemannian(theta)
# J = 2 has a single pair, so the average equals the pairwise AIRM.
all.equal(
  compute_average_riemannian(theta),
  spd_affine_invariant_distance(a, b)
)

---

Compute deconvolution benchmark metrics

Description

Returns three blocks: composition / regression scores, Monte Carlo summaries (ADEMP-style), and optimisation / runtime diagnostics. When $\boldsymbol{p}^{\star}$ is a matrix $J\times N$ (or $\hat{\boldsymbol{p}}$ is), cell-type Pearson correlation and Monte Carlo summaries are computed across samples, not within a single composition.

Usage

compute_benchmark_metrics(
  y,
  mean_signature_matrix,
  estimated_p,
  true_ratios = NULL,
  se = NULL,
  lower = NULL,
  upper = NULL,
  elapsed = NULL,
  memory_bytes = NULL,
  kkt_residual = NULL,
  numerical_converged = NULL,
  loglik_hat = NULL,
  loglik_true = NULL,
  converged = NULL,
  local_maximum = NULL,
  min_eigenvalue = NULL,
  max_eigenvalue = NULL,
  presence_threshold = 1e-04,
  level = 0.95,
  coverage_interval = "wilson",
  algorithm = NA_character_
)

Arguments

y	Numeric vector (or one-column matrix) $\boldsymbol{y}\in\mathbb{R}^{G}$.
mean_signature_matrix	Numeric matrix $\boldsymbol{\mu}\in\mathcal{M}_{G\times J}$ (plug-in means).
estimated_p	Estimated proportions: a length-$J$ vector or a numeric matrix $J\times N$.
true_ratios	Optional ground truth with the same layout. When omitted, reconstitution scores compare $\hat{\boldsymbol{y}}=\boldsymbol{\mu}\hat{\boldsymbol{p}}$ with $\boldsymbol{y}$.
se	Optional standard errors matching estimated_p. Convolution-likelihood solvers fill these from vcov_ilr_delta() (expected Fisher information pulled back through the ILR chart). Mean-only solvers (NNLS, LSEI) leave them missing.
lower, upper	Optional interval bounds matching estimated_p (used for coverage and mean width). When omitted but se is supplied, Wald intervals at level are used.
elapsed, memory_bytes, kkt_residual	Optional per-sample runtime diagnostics (recycled to $N$).
numerical_converged	Optional logical per sample: the solver returned a finite simplex vector.
loglik_hat, loglik_true	Optional per-sample log-likelihoods used for theoretical convergence (regret $\ell(\boldsymbol{p}^{\star})-\ell(\hat{\boldsymbol{p}})$).
converged	Optional list of length $N$ with homogenised optimiser fields (iterations, criterion, code, rdm). Closed-form solvers store NA entries.
local_maximum	Optional logical per sample: ILR score small and Hessian negative definite ($\lambda_{\min}<0$ and $\lambda_{\max}<0$).
min_eigenvalue, max_eigenvalue	Optional ILR Hessian extrema.
presence_threshold	Threshold $\varepsilon$ for presence / F1 / false-positive mass (default 1e-4).
level	Wald coverage level when lower / upper are omitted (default 0.95).
coverage_interval	Interval for the Monte Carlo coverage rate $\hat\pi$: "wilson" (default), "wald", or "agresti_coull". This is not the interval for $p_j$ itself.
algorithm	Character label stored on each metrics row (typically the solver name). Recycled to length one.

Value

A named list with:

• regression: global (one row per sample) and cell_type (one row per cell type);

• monte_carlo: one row per cell type, including theoretical_se (filled later from expected Fisher at $p^{\star}$);

• optimisation: one row per sample.

Examples

mu <- matrix(c(20, 22, 22, 20), nrow = 2,
             dimnames = list(paste0("g", 1:2), paste0("ct", 1:2)))
y <- drop(mu %*% c(0.4, 0.6))
compute_benchmark_metrics(y, mu, estimated_p = c(0.45, 0.55),
                          true_ratios = c(0.4, 0.6))

---

Gene scores from multinomial elastic-net cell-type classification

Description

Fits a multinomial (or binomial) elastic net (glmnet::glmnet()) that predicts cell type from expression features. Inputs are purified expression profiles $\boldsymbol{X}\in\mathcal{M}_{G\times J\times N}$ (genes $\times$ cell types $\times$ samples) and length-$J$ cell-type labels. Variability across samples replaces synthetic isotropic noise. Gene scores are the sum over classes of absolute coefficients at a chosen $\lambda$ (intercept excluded). For nested / CV selection of $\lambda$, see the experimental compute_glmnet_gene_scores_cv() helper (not shipped in the package build).

Usage

compute_glmnet_gene_scores(
  expression_profiles,
  celltype_labels,
  alpha = 0.5,
  lambda = NULL,
  ...
)

Arguments

expression_profiles	Numeric array $G\times J\times N$ of purified profiles.
celltype_labels	Character or factor labels of length $J$ (one per cell-type slice).
alpha	Elastic-net mixing parameter in $[0,1]$ (default 0.5).
lambda	Optional penalty value at which coefficients are extracted. When NULL, uses the smallest $\lambda$ on the fitted path.
...	Additional arguments forwarded to glmnet::glmnet().

Value

Named numeric vector of length $G$ (gene scores; larger means stronger multinomial signal).

See Also

compute_average_jeffreys(), compute_average_overlap()

Examples

## Not run: 
set.seed(1)
G <- 4L
J <- 3L
N <- 12L
profiles <- array(0, dim = c(G, J, N))
for (j in seq_len(J)) {
  profiles[j, j, ] <- 5 + stats::rnorm(N, sd = 0.2)
}
labels <- paste0("ct", seq_len(J))
scores <- compute_glmnet_gene_scores(profiles, labels)
names(scores)[which.max(scores)]

## End(Not run)

---

Normalised Shannon entropy of a discrete distribution

Description

For a probability vector $\boldsymbol{p}\in\Delta^{J-1}$ over $J$ classes (here: cell types), the Shannon entropy is

$H(\boldsymbol{p}) = -\sum_{j=1}^{J} p_j \log p_j,$

with the convention $0\log 0 = 0$. Dividing by the maximum entropy $\log J$ (uniform over all $J$ classes) yields normalised Shannon entropy

$H^{\star}(\boldsymbol{p}) = \frac{H(\boldsymbol{p})}{\log J} \in[0,1],$

so $H^{\star}=0$ for a Dirac mass on one type and $H^{\star}=1$ for the uniform distribution over the $J$ cell types. Zero masses are dropped only inside the sum; the normaliser still uses the original class count $J$.

This is preferable to changing the logarithm base to the number of positive masses $J'$, which would renormalise only on the support and hide sparsity relative to the full panel of cell types.

Usage

compute_shannon_entropy(ratios)

Arguments

ratios

Numeric vector of cellular proportions $\boldsymbol{p}\in\Delta^{J-1}$ (length $J$).

Value

Scalar normalised entropy $H^{\star}\in[0,1]$.

Illustration

Panel B of the figure contrasts a peaked (type-specific, low entropy) share vector with a flat (multi-type, high entropy) one for $J=5$ classes.

See Also

composition_from_entropy()

Examples

compute_shannon_entropy(c(1, 0, 0)) # Dirac -> 0
compute_shannon_entropy(rep(1 / 3, 3)) # uniform -> 1

---

Profile-likelihood confidence intervals for cellular ratios

Description

Inverts the likelihood-ratio test of lrt_decovart() coordinate by coordinate: the interval collects every $c$ with $2\{\ell(\hat{\boldsymbol{p}})-\ell^{\mathrm{prof}}_j(c)\} \le\chi^{2}_{1,1-\alpha}$. Because likelihood ratios are invariant under reparametrisation, the interval respects $0\le p_j\le 1$ without the delta-method linearisation that makes Wald intervals unreliable near a simplex face.

Usage

confint_profile_decovart(
  y = NULL,
  mean_signature_matrix,
  Sigma,
  bulk_expression = NULL,
  level = 0.95,
  celltypes = NULL,
  epsilon = 10^-8,
  itmax = 500L
)

Arguments

y	Numeric vector (or one-column matrix) $\boldsymbol{y}\in\mathbb{R}^{G}$.
mean_signature_matrix	Numeric matrix $\boldsymbol{\mu}\in\mathcal{M}_{G\times J}$ (plug-in means).
Sigma	Array of cell-type covariances in $\mathcal{M}_{G\times G\times J}$.
bulk_expression	Numeric matrix $\boldsymbol{Y}\in \mathcal{M}_{G\times N}$ of replicate bulk profiles sharing one composition, used instead of y for the pooled (population) test.
level	Confidence level $1-\alpha$.
celltypes	Cell types to profile; defaults to all.
epsilon, itmax	Relative convergence tolerance and iteration budget passed to stats::optim().

Value

A matrix with one row per cell type and columns estimate, lower, upper.

See Also

confint.decovart_fit() for the Wald analogue, lrt_decovart()

Other decovart_inference: bootstrap_decovart(), boundary_diagnostics(), chi_bar_square_pvalue(), equivariance_check_decovart(), lrt_decovart(), multistart_decovart(), profile_loglik_decovart(), reference_bootstrap_decovart(), restricted_mle_decovart()

Examples

mu <- matrix(c(20, 40, 15, 40, 20, 25), nrow = 3)
colnames(mu) <- paste0("ct", 1:2)
Sigma <- array(c(diag(3), diag(3)), dim = c(3, 3, 2))
y <- drop(mu %*% c(0.6, 0.4))
confint_profile_decovart(y, mu, Sigma)

---

Monte Carlo interval for an empirical coverage probability

Description

Interval for a binomial coverage rate $\hat\pi=X/N$ from independent coverage indicators. Wilson is the default because the parameter is bounded in $[0,1]$ and Wald intervals behave poorly near 0 or 1.

Usage

coverage_mc_interval(
  covered,
  conf_level = 0.95,
  method = c("wilson", "wald", "agresti_coull")
)

Arguments

covered	Logical or 0/1 coverage indicators (one per replicate).
conf_level	Confidence level for the interval around $\hat\pi$ (default 0.95).
method	"wilson" (default), "wald", or "agresti_coull".

Value

A list with n, successes, coverage, mcse, lower, upper, and method.

Examples

coverage_mc_interval(c(TRUE, TRUE, TRUE, FALSE))

---

Parallel deconvolution of a bulk expression matrix

Description

For each column $\boldsymbol{y}_{\cdot i}$ of the bulk matrix $\boldsymbol{Y}\in\mathcal{M}_{G\times N}$, estimates $\hat{\boldsymbol{p}}_{\cdot i}$ with every supplied algorithm. Samples are iterated with furrr when cores > 1; algorithms are sequential, so workers are never nested. When covariance information is provided, DeCovarT methods maximise $\ell_{\boldsymbol{y}\,|\,\boldsymbol{\zeta}}(\boldsymbol{p})$ under $\boldsymbol{y}\,|\,(\boldsymbol{\zeta},\boldsymbol{p})\sim \mathcal{N}_{G}( \boldsymbol{\mu}\boldsymbol{p}, \boldsymbol{\Sigma}(\boldsymbol{p}) )$.

Usage

deconvolute_ratios(
  signature_matrix,
  bulk_expression,
  true_ratios = NULL,
  Sigma = NULL,
  deconvolution_functions = NULL,
  standardise = FALSE,
  scaled = FALSE,
  cores = getOption("mc.cores", 1L),
  coverage_interval = "wilson",
  verbose = FALSE,
  progress_every = 10L
)

Arguments

signature_matrix	Mean signature $\boldsymbol{\mu}\in\mathcal{M}_{G\times J}$ (rownames = genes, colnames = cell types).
bulk_expression	Bulk matrix $\boldsymbol{Y}\in\mathcal{M}_{G\times N}$ (rownames = genes, colnames = samples).
true_ratios	Optional ground-truth proportions: a length-$J$ numeric vector (recycled over samples) or a numeric matrix $J\times N$ (or $N\times J$).
Sigma	Optional array $(\boldsymbol{\Sigma}_j)_{j}\in\mathcal{M}_{G\times G\times J}$ of cell-type covariances (numeric; off-diagonals may be negative).
deconvolution_functions	Named list; each element has FUN and optional additional_parameters.
standardise	If TRUE, apply a gene-wise affine z-score computed once from $\boldsymbol{\mu}$ to bulk, means and covariances (see Details). Cell-type-wise or sample-wise transforms are not supported.
scaled	Deprecated. TRUE (log2 mixing) always errors.
cores	Number of furrr workers for the sample loop. Defaults to getOption("mc.cores", 1L). Use cores = 1 to stay sequential (CRAN examples and nested callers). With cores > 1, .map_samples() uses furrr_options(seed = TRUE), which assigns an independent L'Ecuyer-CMRG stream per worker (R 4.1 / future streams). Do not offset a shared seed by the sample index.
coverage_interval	Passed to compute_benchmark_metrics() ("wilson", "wald", or "agresti_coull").
verbose	If TRUE, emit progress for sequential sample loops (every progress_every samples). Parallel furrr runs log a one-line notice instead of per-sample ticks.
progress_every	Sample-progress interval when verbose is TRUE and cores = 1. Defaults to 10L.

Value

A named list from compute_benchmark_metrics() with regression (global and cell-type subtables), monte_carlo, and optimisation (per-sample elapsed time, memory, KKT residual, and $\hat{\boldsymbol{p}}$). First-generation solvers still call repair_simplex(); the three native DeCovarT maps (ILR or $p/\sum p$) already lie on the simplex. When Sigma is supplied, Newton-Raphson, Marquardt-Levenberg, L-BFGS-B, BFGS (gradient), and simulated annealing attach ILR Wald standard errors from vcov_ilr_delta() so monte_carlo can report coverage, mean model SE, and the SE/SD ratio. Those SEs are the expected-Fisher delta-method bounds used by confint.decovart_fit(), not diagonal entries of the observed Hessian. Mean-only solvers (NNLS, LSEI) leave Wald columns missing: those estimators do not use the convolution likelihood.

References

Erkkilä T, Lehmusvaara S, Ruusuvuori P, Visakorpi T, Shmulevich I, Lähdesmäki H (2010). “Probabilistic Analysis of Gene Expression Measurements from Heterogeneous Tissues.” Bioinformatics, 26. doi:10.1093/bioinformatics/btq406. Ahn J, Yuan Y, Parmigiani G, Suraokar MB, Diao L, Wistuba II, Wang W (2013). “DeMix: Deconvolution for Mixed Cancer Transcriptomes Using Raw Measured Data.” Bioinformatics (Oxford, England), 29. doi:10.1093/bioinformatics/btt301. Wang Z, Cao S, Morris JS, Ahn J, Liu R, Tyekucheva S, Gao F, Li B, Lu W, Tang X, Wistuba II, Bowden M, Mucci L, Loda M, Parmigiani G, Holmes CC, Wang W (2018). “Transcriptome Deconvolution of Heterogeneous Tumor Samples with Immune Infiltration.” iScience, 9. doi:10.1016/j.isci.2018.10.028. Anghel CV, Quon G, Haider S, Nguyen F, Deshwar AG, Morris QD, Boutros PC (2015). “ISOpureR: An R Implementation of a Computational Purification Algorithm of Mixed Tumour Profiles.” BMC Bioinformatics, 16. doi:10.1186/s12859-015-0597-x. Ogundijo OE, Wang X (2017). “A Sequential Monte Carlo Approach to Gene Expression Deconvolution.” PLOS ONE, 12. doi:10.1371/journal.pone.0186167. Hafemeister C, Satija R (2019). “Normalization and variance stabilization of single-cell RNA-seq data using regularized negative binomial regression.” Genome Biology, 20(1), 296. doi:10.1186/s13059-019-1874-1. Chion M, Leroy A (2023). “A Bayesian Framework for Multivariate Differential Analysis Accounting for Missing Data.” arXiv. doi:10.48550/arxiv.2307.08975. Newman A, Liu C, Green M, Gentles A, Feng W, Xu Y, Hoang C, Diehn M, Alizadeh A (2015). “Robust Enumeration of Cell Subsets from Tissue Expression Profiles.” Nature methods, 12. doi:10.1038/nmeth.3337.

See Also

deconvolute_ratios_Marquardt_Levenberg()

Examples

set.seed(1)
genes <- paste0("g", 1:2)
cts <- paste0("ct", 1:2)
mu <- matrix(c(20, 22, 22, 20), nrow = 2, dimnames = list(genes, cts))
Sigma <- array(
  c(1, 0, 0, 1, 1, 0, 0, 1),
  dim = c(2, 2, 2),
  dimnames = list(genes, genes, cts)
)
Y <- simulate_bulk_mixture(mu, Sigma, p = c(0.5, 0.5), n = 2)$Y
deconvolute_ratios(
  signature_matrix = mu,
  bulk_expression = Y,
  true_ratios = c(0.5, 0.5),
  Sigma = Sigma,
  deconvolution_functions = list(
    "nnls" = list(FUN = deconvolute_ratios_nnls)
  ),
  cores = 1
)

---

DeCovarT MLE of cellular proportions for one bulk sample

Description

Estimates $\hat{\boldsymbol{p}}=\arg\max_{\boldsymbol{p}} \ell_{\boldsymbol{y}\,|\,\boldsymbol{\zeta}}(\boldsymbol{p})$ under the Gaussian convolution model

$\boldsymbol{y}\,|\,(\boldsymbol{\zeta},\boldsymbol{p}) \sim\mathcal{N}_{G}\!\bigl(\boldsymbol{\mu}\boldsymbol{p},\, \boldsymbol{\Sigma}(\boldsymbol{p})\bigr), \qquad \boldsymbol{\Sigma}(\boldsymbol{p})=\sum_{j=1}^{J}p_j^{2}\boldsymbol{\Sigma}_j,$

subject to the simplex constraint $\mathbf{1}^{\mathsf{T}}\boldsymbol{p}=1$, $\boldsymbol{p}\ge\mathbf{0}$. Optimisation is performed in unconstrained ILR coordinates $\boldsymbol{z}\in\mathbb{R}^{J-1}$ via $\boldsymbol{p}=\operatorname{softmax}(\mathbf{V}\boldsymbol{z})$ (Marquardt–Levenberg default; see other methods below and vignette("theory-decovart-generative-model", package = "DeCovarT")).

Usage

deconvolute_ratios_cibersort(y, mean_signature_matrix)

deconvolute_ratios_rlm(y, mean_signature_matrix)

deconvolute_ratios_nnls(y, mean_signature_matrix)

deconvolute_ratios_deconrnaseq(y, mean_signature_matrix)

deconvolute_ratios_Marquardt_Levenberg(
  y,
  mean_signature_matrix,
  Sigma,
  epsilon = 10^-4,
  itmax = 200,
  return_model = FALSE,
  initial_p = NULL,
  dirichlet_alpha = 1
)

deconvolute_ratios_simulated_annealing(
  y,
  mean_signature_matrix,
  Sigma,
  epsilon = 10^-4,
  itmax = 200,
  return_model = FALSE,
  initial_p = NULL,
  dirichlet_alpha = 1
)

deconvolute_ratios_L_BFGS_B(
  y,
  mean_signature_matrix,
  Sigma,
  epsilon = 10^-4,
  itmax = 200,
  return_model = FALSE,
  initial_p = NULL,
  dirichlet_alpha = 1
)

deconvolute_ratios_Newton_Raphson(
  y,
  mean_signature_matrix,
  Sigma,
  epsilon = 10^-4,
  itmax = 200,
  return_model = FALSE,
  initial_p = NULL,
  dirichlet_alpha = 1
)

deconvolute_ratios_gradient_descent(
  y,
  mean_signature_matrix,
  Sigma,
  epsilon = 10^-4,
  itmax = 200,
  return_model = FALSE,
  initial_p = NULL,
  dirichlet_alpha = 1
)

Arguments

y	Numeric vector (or one-column matrix) $\boldsymbol{y}\in\mathbb{R}^{G}$.
mean_signature_matrix	Numeric matrix $\boldsymbol{\mu}\in\mathcal{M}_{G\times J}$ (plug-in means).
Sigma	Array of cell-type covariances in $\mathcal{M}_{G\times G\times J}$.
epsilon, itmax	Absolute convergence tolerance and maximum number of iterations for the optimiser (same roles as reltol / maxit in stats::optim()).
return_model	If TRUE, return a named list with coefficients, ILR coordinates, log-likelihood and optimiser diagnostics instead of the proportion vector.
initial_p	Optional start: NULL / "barycentre" (default equi-balanced), a length-$J$ numeric vector, "dirichlet", or "qp" (mean-only simplex QP). See starting_simplex(). Starts on a simplex face are nudged into the interior so the ILR map is defined (ILR methods) and so L-BFGS-B does not start on a degenerate $\boldsymbol{\Sigma}(\boldsymbol{p})$.
dirichlet_alpha	Dirichlet concentration when initial_p = "dirichlet" (default 1: uniform on the simplex; $\alpha>1$ centre-biased; $\alpha<1$ face-biased). Independent draws are independent restarts; use multistart_decovart() to sequence several.

Details

Plug-in signature. Argument mean_signature_matrix is the mean $\boldsymbol{\mu}$, used as a proxy for the unobserved cell-type profiles $\boldsymbol{x}_{\cdot j}$. This is the frequentist plug-in; recovering sample-specific latents is a Bayesian / MAP problem (.map_gaussian_convolution()).

Value

Named numeric vector $\hat{\boldsymbol{p}}$ on the simplex (ILR methods), or that list when return_model = TRUE. Benchmark metrics are computed by deconvolute_ratios().

Functions

• deconvolute_ratios_cibersort(): Linear baseline $\hat{\boldsymbol{y}}=\boldsymbol{\mu}\hat{\boldsymbol{p}}$ via nu-SVR (CIBERSORT-style); no covariance prior is used.

• deconvolute_ratios_rlm(): Robust linear model $\boldsymbol{y}\approx\boldsymbol{\mu}\boldsymbol{p}$ (MASS::rlm()), as in Monaco et al. (2019).

• deconvolute_ratios_nnls(): Non-negative least squares for $\boldsymbol{y}\approx\boldsymbol{\mu}\boldsymbol{p}$ (nnls::nnls()), then simplex projection.

• deconvolute_ratios_deconrnaseq(): Equality- and inequality-constrained least squares on the simplex (limSolve::lsei()), in the spirit of deconRNASeq.

• deconvolute_ratios_simulated_annealing(): Simulated annealing on ILR coordinates $\boldsymbol{z}$ (stats::optim() with method = "SANN").

• deconvolute_ratios_L_BFGS_B(): Box-constrained L-BFGS-B directly in $\boldsymbol{p}$ (stats::optim() method = "L-BFGS-B"). The box keeps each coordinate in $[0,1]$; the returned vector is closed by $p/\sum p$ (no repair_simplex() clipping). Wald standard errors still use the ILR expected-Fisher map vcov_ilr_delta(), not a $\boldsymbol{p}$-space Hessian.

• deconvolute_ratios_Newton_Raphson(): Newton–Raphson / nlminb on ILR coordinates $\boldsymbol{z}$ using analytic gradient and Hessian (stats::nlminb()).

• deconvolute_ratios_gradient_descent(): BFGS quasi-Newton ascent on ILR coordinates $\boldsymbol{z}$ (stats::optim() method = "BFGS").

References

Monaco G, Lee B, Xu W, Mustafah S, Hwang YY, Carré C, Burdin N, Visan L, Ceccarelli M, Poidinger M, Zippelius A, Pedro de Magalhães J, Larbi A (2019). “RNA-Seq Signatures Normalized by mRNA Abundance Allow Absolute Deconvolution of Human Immune Cell Types.” Cell Reports, 26. doi:10.1016/j.celrep.2019.01.041.

See Also

deconvolute_ratios(), isometric_logistic(), starting_simplex()

Examples

set.seed(1)
genes <- paste0("g", 1:2)
cts <- paste0("ct", 1:2)
mu <- matrix(c(20, 22, 22, 20), nrow = 2, dimnames = list(genes, cts))
Sigma <- array(
  c(1, 0, 0, 1, 1, 0, 0, 1),
  dim = c(2, 2, 2),
  dimnames = list(genes, genes, cts)
)
y <- drop(mu %*% c(0.6, 0.4) + rnorm(2, sd = 0.1))
deconvolute_ratios_Marquardt_Levenberg(y, mu, Sigma, itmax = 50)

---

Generalised least squares with a fixed residual covariance

Description

Fits $\boldsymbol{y}\approx\boldsymbol{\mu}\boldsymbol{p}$ by MASS::lm.gls() with a known $G\times G$ residual covariance W that does not depend on $\boldsymbol{p}$, then repair_simplex(). Use fixed_gls_covariance() for the global-diagonal competitor $\operatorname{diag}\{\sum_j \bar p_j^2\boldsymbol{\Sigma}_j\}$. This is the mean-only GLS gold standard for covariance-structure benchmarks, not the DeCovarT convolution $\boldsymbol{\Sigma}(\boldsymbol{p})=\sum_j p_j^2\boldsymbol{\Sigma}_j$.

Usage

deconvolute_ratios_gls(y, mean_signature_matrix, W, inverse = TRUE)

Arguments

y	Numeric bulk vector of length $G$.
mean_signature_matrix	Numeric matrix $G\times J$.
W	Residual covariance (when inverse = TRUE) or precision (when inverse = FALSE), of size $G\times G$.
inverse	Passed to MASS::lm.gls(); TRUE means W is a covariance matrix.

Value

Named simplex vector of length $J$.

See Also

fixed_gls_covariance(), deconvolute_ratios_deconrnaseq()

Examples

mu <- cbind(ct1 = c(20, 40), ct2 = c(40, 20))
y <- drop(mu %*% c(0.4, 0.6))
Sigma <- array(c(diag(2), 2 * diag(2)), dim = c(2, 2, 2))
deconvolute_ratios_gls(y, mu, fixed_gls_covariance(Sigma))

---

Describe a Gaussian-convolution simulation scenario

Description

Summarises the generating law $\boldsymbol{y}\mid(\boldsymbol{\zeta},\boldsymbol{p})\sim \mathcal{N}_{G}(\boldsymbol{\mu}\boldsymbol{p}, \sum_j p_j^2\boldsymbol{\Sigma}_j)$ at five layers: composition, mean geometry, covariance / SPD diagnostics, network sparsity, and tangent Fisher information (mean versus covariance split). MixSim BarOmega and average pairwise Hellinger of the component Gaussians are kept in descriptors (Hellinger is the unweighted pairwise mean; hellinger_weighted is optional). Jeffreys / symmetrised KL is returned in supplementary. Mean-geometry columns include the highest pairwise cosine and the mean Euclidean (CLD) gap.

Usage

describe_simulation_scenario(
  true_theta,
  adjacency = NULL,
  active_tol = 1e-08,
  include_mixsim = TRUE
)

Arguments

true_theta	List with p, mu, and sigma (and optionally Theta precision). Parsed by check_true_theta() / .parse_true_theta().
adjacency	Optional $G\times G\times J$ adjacency array, or a length-$J$ list of adjacencies, used for network density when supplied.
active_tol	Threshold for counting an active simplex component.
include_mixsim	Logical. If TRUE (default), compute MixSim BarOmega via compute_average_overlap(). Set FALSE to skip.

Value

A list with:

• theta_true: the convolution parameters p, mu, sigma;

• descriptors: one-row tibble of kept scenario statistics in six families (composition, mean geometry, SPD of $\boldsymbol{\Sigma}(\boldsymbol{p})$, tangent Fisher, network, component overlap). MixSim BarOmega uses compute_average_overlap() (Davies quadrature for $G<4$, Sobol Monte Carlo otherwise). Average pairwise Hellinger and affine-invariant Riemannian distance of the $\Sigma_j$ are kept alongside. hellinger_weighted is optional. SPD columns include both $\kappa\{\boldsymbol{\Sigma}(\boldsymbol{p})\}$ (kappa_sigma_p) and the reciprocal $\lambda_{\min}/\lambda_{\max}$ (kappa_sigma_reciprocal);

• supplementary: one-row tibble of Jeffreys / symmetrised KL, recorded but not treated as a primary score;

• call: the matched call (match.call()).

See Also

run_simulation_benchmark(), expected_fisher_unconstrained(), compute_shannon_entropy(), composition_from_entropy(), helmert_basis(), compute_average_overlap(), compute_average_riemannian()

Examples

theta <- list(
  p = c(0.5, 0.5),
  mu = cbind(c(10, 0), c(0, 10)),
  sigma = array(c(diag(2), diag(2)), dim = c(2, 2, 2))
)
out <- describe_simulation_scenario(theta)
out$descriptors$h_star
out$call

---

Equicorrelation (constant-correlation) Gram matrix

Description

$R=(1-\rho)I+\rho\mathbf{1}\mathbf{1}^{\mathsf{T}}$. Positive semidefinite for $\rho\in[-1/(J-1),1]$.

Usage

equicorrelation_gram(n_celltypes, target_cosine = 0)

Arguments

n_celltypes	Integer $J\ge 2$.
target_cosine	Pairwise cosine $\rho$.

Value

Symmetric $J\times J$ matrix with unit diagonal.

Examples

equicorrelation_gram(3L, 0.4)

---

Expected Fisher information of unconstrained $\boldsymbol{p}$

Description

For the multivariate-normal mean–covariance map of the DeCovarT convolution, $\boldsymbol{y}\sim\mathcal{N}_{G}(\boldsymbol{\mu}\boldsymbol{p}, \boldsymbol{\Sigma}(\boldsymbol{p}))$ with $\boldsymbol{\Sigma}(\boldsymbol{p})=\sum_j p_j^{2}\boldsymbol{\Sigma}_j$ and precision $\boldsymbol{\Omega}(\boldsymbol{p})=\boldsymbol{\Sigma}(\boldsymbol{p})^{-1}$, the expected Fisher information has entries

$I(\boldsymbol{p})_{jk} = \boldsymbol{\mu}_{\cdot j}^{\top} \boldsymbol{\Omega}(\boldsymbol{p}) \boldsymbol{\mu}_{\cdot k} + 2 p_j p_k\, \mathrm{tr}\bigl( \boldsymbol{\Omega}(\boldsymbol{p})\boldsymbol{\Sigma}_j \boldsymbol{\Omega}(\boldsymbol{p})\boldsymbol{\Sigma}_k \bigr).$

The first summand is the mean contribution (an $\boldsymbol{\Omega}$-inner product of signature columns); the second is the covariance contribution of the quadratic map $\boldsymbol{p}\mapsto\boldsymbol{\Sigma}(\boldsymbol{p})$. See the multivariate-normal formula on https://en.wikipedia.org/wiki/Fisher_information#Multivariate_normal_distribution.

Usage

expected_fisher_unconstrained(p, mean_signature_matrix, Sigma)

Arguments

p	Numeric proportions on the open simplex.
mean_signature_matrix	Mean signature $\boldsymbol{\mu}$ ($G\times J$).
Sigma	Cell-type covariances $G\times G\times J$.

Value

Symmetric $J\times J$ expected Fisher information matrix $I(\boldsymbol{p})$.

See Also

vcov_ilr_delta(), vcov.decovart_fit(), confint.decovart_fit(), .inner_product()

Examples

# Two cell types, bivariate Gaussians: I(p) is 2 x 2.
p <- c(0.6, 0.4)
mu <- cbind(c(0, 0), c(3, 0))
Sigma <- array(c(diag(2), diag(2)), dim = c(2, 2, 2))
I <- expected_fisher_unconstrained(p, mu, Sigma)
I

---

Fit the DeCovarT Gaussian-convolution model

Description

Maximum-likelihood estimator of cellular proportions $\boldsymbol{p}_{\cdot i}$ under the multivariate Gaussian convolution

$\boldsymbol{y}_{\cdot i}\mid(\boldsymbol{\mu},\boldsymbol{\Sigma}_j, \boldsymbol{p}_{\cdot i}) \sim \mathcal{N}_{G}\bigl( \boldsymbol{\mu}\boldsymbol{p}_{\cdot i},\; \boldsymbol{\Sigma}(\boldsymbol{p}_{\cdot i}) \bigr), \qquad \boldsymbol{\Sigma}(\boldsymbol{p}) =\sum_{j=1}^{J}p_j^{2}\boldsymbol{\Sigma}_j.$

This is a variance / likelihood specification, not ordinary least squares: cell-type covariances enter the residual law and cannot be written as extra columns of $\boldsymbol{\mu}$. There is therefore no formula / lm interface, and no predict() method for forecasting bulk expression.

The returned object is an S3 class decovart_fit with accessors coef(), fitted(), residuals(), vcov(), nobs(), confint(), print(), summary() and plot(). Residuals are $\boldsymbol{Y}-\boldsymbol{\mu}\hat{\boldsymbol{P}}$; they are not OLS residuals. Goodness of fit is the convolution log-likelihood, not residual sum of squares.

Usage

fit_decovart(
  signature_matrix,
  bulk_expression,
  true_ratios = NULL,
  Sigma = NULL,
  method = c("Marquardt-Levenberg", "L-BFGS-B", "Newton-Raphson"),
  epsilon = 10^-4,
  itmax = 200,
  standardise = FALSE,
  scaled = FALSE,
  n_starts = 0L,
  boundary_tol = 1e-08
)

## S3 method for class 'decovart_fit'
print(x, ...)

## S3 method for class 'decovart_fit'
summary(object, ...)

## S3 method for class 'summary.decovart_fit'
print(x, ...)

## S3 method for class 'decovart_fit'
coef(object, ...)

## S3 method for class 'decovart_fit'
fitted(object, ...)

## S3 method for class 'decovart_fit'
residuals(object, ...)

## S3 method for class 'decovart_fit'
vcov(object, ...)

## S3 method for class 'decovart_fit'
nobs(object, ...)

## S3 method for class 'decovart_fit'
confint(object, parm, level = 0.95, ...)

## S3 method for class 'decovart_fit'
plot(x, ...)

Arguments

signature_matrix	Mean signature $\boldsymbol{\mu}\in\mathcal{M}_{G\times J}$ (rownames = genes, colnames = cell types).
bulk_expression	Bulk matrix $\boldsymbol{Y}\in\mathcal{M}_{G\times N}$ (rownames = genes, colnames = samples).
true_ratios	Optional ground-truth proportions ($J$ vector or $J\times N$ matrix). Stored on the fit; not used by the MLE.
Sigma	Array $(\boldsymbol{\Sigma}_j)_{j}\in\mathcal{M}_{G\times G\times J}$ of cell-type covariances.
method	Optimiser; one of "Marquardt-Levenberg", "L-BFGS-B", "Newton-Raphson" (case-insensitive). These three maps already land on the simplex (ILR or $p/\sum p$); they do not call repair_simplex().
epsilon, itmax	Absolute convergence tolerance and iteration budget, in the same roles as reltol / maxit in stats::optim().
standardise	If TRUE, apply a gene-wise affine z-score computed once from $\boldsymbol{\mu}$ to bulk, means and covariances (see Details). Cell-type-wise or sample-wise transforms are not supported.
scaled	Deprecated. TRUE (log2 mixing) always errors.
n_starts	Number of additional random Dirichlet starts per sample. With n_starts > 0 the best log-likelihood is kept and multistart_decovart() records the spread of attained optima, the only direct probe of multimodality (the realised log-likelihood is not globally concave).
boundary_tol	Threshold on $\min_j\hat{p}_j$ below which a sample is flagged near_boundary by boundary_diagnostics(). Wald intervals are unreliable there; prefer confint_profile_decovart() or lrt_decovart().
x, object	A decovart_fit.
...	Passed to graphics::plot() (plot method) or ignored.
parm	Kept for compatibility with stats::confint(); all simplex coordinates are always returned (subsetting by parm is not implemented). Removing this formal would break the S3 method contract with the generic.
level	Confidence level $1-\alpha$ (default 0.95). Wald intervals use asymptotic normality of the MLE with standard errors from vcov_ilr_delta() / expected_fisher_unconstrained() (see Details of fit_decovart() and of vcov_ilr_delta()): $\hat{p}_j \pm z_{1-\alpha/2}\,\mathrm{SE}_j$.

Details

Standardisation. CIBERSORT requires non-negative expression, no missing values, and a non-log linear scale (Newman et al. 2015). A logarithm is concave, so Jensen's inequality shifts first and second moments and $\log(\boldsymbol{\mu}\boldsymbol{p})\neq(\log\boldsymbol{\mu})\boldsymbol{p}$. A gene-wise affine map $\boldsymbol{x}\mapsto D^{-1}(\boldsymbol{x}-\boldsymbol{m})$ with the same $D,\boldsymbol{m}$ on $\boldsymbol{Y}$, $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}_j^\star=D^{-1}\boldsymbol{\Sigma}_j D^{-1}$ leaves the MLE of $\boldsymbol{p}$ unchanged (equivariance).

Wald covariance. Let $\boldsymbol{\Omega}(\boldsymbol{p})=\boldsymbol{\Sigma}(\boldsymbol{p})^{-1}$. The expected Fisher information of the unconstrained mean–covariance map (multivariate normal; see e.g. the Wikipedia entry Fisher information, multivariate normal) is

$I(\boldsymbol{p})_{jk} = \boldsymbol{\mu}_{\cdot j}^{\top} \boldsymbol{\Omega}(\boldsymbol{p}) \boldsymbol{\mu}_{\cdot k} + 2 p_j p_k\, \mathrm{tr}\bigl( \boldsymbol{\Omega}(\boldsymbol{p})\boldsymbol{\Sigma}_j \boldsymbol{\Omega}(\boldsymbol{p})\boldsymbol{\Sigma}_k \bigr).$

Cramer–Rao gives $\mathrm{Var}(\hat{\boldsymbol{z}})\succeq I_{\boldsymbol{z}}^{-1}$ in ILR coordinates, with $I_{\boldsymbol{z}} =\mathbf{J}_{\boldsymbol{\psi}}^{\top} I(\boldsymbol{p}) \mathbf{J}_{\boldsymbol{\psi}}$ and $\mathbf{J}_{\boldsymbol{\psi}} =\mathbf{S}(\boldsymbol{p})\mathbf{V}$ (jacobian_isometric_logistic()). The delta method maps the bound back to the simplex:

$\mathrm{Var}(\hat{\boldsymbol{p}}) = \mathbf{J}_{\boldsymbol{\psi}} I_{\boldsymbol{z}}^{-1} \mathbf{J}_{\boldsymbol{\psi}}^{\mathsf{T}}.$

Value

An object of class decovart_fit.

References

Newman A, Liu C, Green M, Gentles A, Feng W, Xu Y, Hoang C, Diehn M, Alizadeh A (2015). “Robust Enumeration of Cell Subsets from Tissue Expression Profiles.” Nature methods, 12. doi:10.1038/nmeth.3337.

See Also

deconvolute_ratios(), deconvolute_ratios_Marquardt_Levenberg(), expected_fisher_unconstrained(), vcov_ilr_delta(), coef.decovart_fit(), vcov.decovart_fit(), confint.decovart_fit()

Examples

toy <- readRDS(system.file(
  "extdata", "toy_deconvolution.rds",
  package = "DeCovarT"
))
fit <- fit_decovart(
  signature_matrix = toy$signature_matrix,
  bulk_expression = toy$bulk_expression,
  Sigma = toy$Sigma,
  itmax = 40
)
coef(fit)
nobs(fit)

---

Fixed residual covariance for a GLS competitor

Description

Evaluates the convolution covariance $\boldsymbol{\Sigma}(\boldsymbol{p})=\sum_j p_j^2\boldsymbol{\Sigma}_j$ at a known composition (default the barycentre $p_j=1/J$) and optionally keeps only the diagonal. That matrix does not depend on the unknown $\boldsymbol{p}$ at fit time. Do not copy it into every slice of a DeCovarT tensor: that would yield $\|p\|_2^2\boldsymbol{\Sigma}_{\mathrm{GLS}}$, which still depends on $p$.

Usage

fixed_gls_covariance(Sigma, p = NULL, diagonal = TRUE)

Arguments

Sigma	Array $G\times G\times J$ of cell-type covariances.
p	Optional composition of length $J$; default $(1/J,\ldots,1/J)$.
diagonal	If TRUE (default), return $\operatorname{diag}\{\boldsymbol{\Sigma}(\boldsymbol{p})\}$.

Value

A $G\times G$ covariance matrix.

See Also

deconvolute_ratios_gls()

Examples

Sigma <- array(c(diag(2), 2 * diag(2)), dim = c(2, 2, 2))
fixed_gls_covariance(Sigma)

---

Generate mean profiles from a target Gram matrix

Description

Builds $\boldsymbol{\mu}\in\mathcal{M}_{G\times J}$ whose unit columns realise a prescribed Gram matrix $R$ (pairwise cosines) via the symmetric square root,

$\boldsymbol{\mu} = s\, \boldsymbol{Q} R^{1/2}, \qquad \boldsymbol{Q}^{\mathsf{T}}\boldsymbol{Q}=\boldsymbol{I}_{J}.$

Then $\boldsymbol{\mu}^{\mathsf{T}}\boldsymbol{\mu}=s^{2}R$ whenever $R$ has unit diagonal. The default $R$ is the equicorrelation matrix $(1-\rho)I+\rho\mathbf{1}\mathbf{1}^{\mathsf{T}}$ (equicorrelation_gram()). Supply target_gram to set unequal pairwise cosines (for example two close pairs and a distant background), provided \(R\) is symmetric nonnegative-definite.

Usage

generate_mean_signature_matrix(
  n_genes,
  n_celltypes,
  mean_scale = 10,
  target_cosine = 0,
  target_gram = NULL,
  seed = NULL,
  nonnegative = FALSE,
  gene_names = NULL,
  celltype_names = NULL
)

Arguments

n_genes	Integer $G$; must be at least n_celltypes.
n_celltypes	Integer $J\ge 2$.
mean_scale	Positive scalar $s$ (centroid norms). Default 10. Hold fixed when studying cosine / collinearity alone.
target_cosine	Numeric pairwise cosine $\rho$ for the equicorrelation Gram. Ignored when target_gram is supplied. Must lie in $[-1/(J-1),1]$.
target_gram	Optional symmetric $J\times J$ cosine matrix (unit diagonal). Overrides target_cosine.
seed	Optional integer. When supplied, the orthonormal frame $\boldsymbol{Q}$ is drawn from a Gaussian QR; otherwise it is deterministic. Ignored when nonnegative = TRUE.
nonnegative	If TRUE, use a disjoint-support nonnegative frame so that $\boldsymbol{\mu}\ge 0$ whenever $R^{1/2}\ge 0$ (required by deconvolute_ratios()). The Gram is unchanged.
gene_names	Optional character vector of length $G$.
celltype_names	Optional character vector of length $J$.

Details

Realisations are unique up to an orthogonal transformation of the gene space: replacing $\boldsymbol{Q}$ by $\boldsymbol{Q}U$ with $U^{\mathsf{T}}U=\boldsymbol{I}_{J}$ leaves the Gram unchanged. By default $\boldsymbol{Q}$ is a deterministic thin QR frame; pass seed for a Haar-like Gaussian frame. Cholesky $R=LL^{\mathsf{T}}$ is an alternative square root used in Monte Carlo simulation; see the how-to-build-synthetic-scenarios vignette appendix.

Value

Numeric matrix $\boldsymbol{\mu}$ with dimensions $G\times J$.

Examples

generate_mean_signature_matrix(
  n_genes = 6L,
  n_celltypes = 2L,
  target_cosine = 0.5
)
k <- matrix(c(1, 0.98, 0.2, 0.98, 1, 0.2, 0.2, 0.2, 1), 3, 3)
generate_mean_signature_matrix(
  n_genes = 8L,
  n_celltypes = 3L,
  target_gram = k
)

---

Helmert contrast matrix for isometric log-ratio coordinates

Description

Returns the $J\times(J-1)$ Helmert sub-matrix $\mathbf{V}$ with

$\mathbf{V}^{\mathsf{T}}\mathbf{V}=\mathbf{I}_{J-1},\qquad \mathbf{V}^{\mathsf{T}}\mathbf{1}=\mathbf{0}.$

Column $k$ contrasts the first $k$ parts against part $k+1$. This is the standard ILR basis used with isometric_log_ratio() / isometric_logistic(); any other valid ILR basis is $\mathbf{V}Q$ for an orthogonal $Q$, which rotates coordinates but leaves simplex geometry, eigenvalues of quadratic forms, and $\mathrm{Var}(\hat{\boldsymbol{p}})$ unchanged.

Usage

helmert_basis(n_parts)

Arguments

n_parts

Integer $J\ge 2$ (number of simplex parts).

Value

Numeric matrix $J\times(J-1)$.

Examples

V <- helmert_basis(3L)
crossprod(V)
drop(crossprod(V, rep(1, 3)))

---

Hessian tensor of the isometric logistic map

Description

Slice $i$ is

$\mathbf{H}_{\psi_i} = p_i\bigl(\boldsymbol{q}_i\boldsymbol{q}_i^{\mathsf{T}} -\mathbf{C}_V\bigr),$

where $\bar{\boldsymbol{v}}=\mathbf{V}^{\mathsf{T}}\boldsymbol{p}$, $\boldsymbol{q}_i=\mathbf{V}_{i\cdot}-\bar{\boldsymbol{v}}$ and $\mathbf{C}_V=\mathbf{V}^{\mathsf{T}}\mathbf{S}(\boldsymbol{p}) \mathbf{V}$. Array shape $(J-1)\times(J-1)\times J$.

Usage

hessian_isometric_logistic(z, V = NULL)

Arguments

z	Numeric vector $\boldsymbol{z}\in\mathbb{R}^{J-1}$.
V	Optional ILR basis; see isometric_logistic().

Value

Numeric array used in the constrained Hessian chain rule.

Examples

hessian_isometric_logistic(c(0.2, -0.5))

---

Isometric logistic transform (ILR coordinates to the simplex)

Description

Inverse isometric log-ratio map $\boldsymbol{\psi}:\boldsymbol{z}\mapsto\boldsymbol{p}$ with a Helmert basis $\mathbf{V}$,

$\boldsymbol{p} = \operatorname{softmax}(\mathbf{V}\boldsymbol{z}) = \mathcal{C}\bigl(\exp(\mathbf{V}\boldsymbol{z})\bigr).$

No cell type is pinned as a reference (unlike additive_logistic()). Evaluated with a log-sum-exp shift so large $\|\boldsymbol{z}\|$ cannot overflow.

$\boldsymbol{z} = \mathbf{V}^{\mathsf{T}}\log\boldsymbol{p} = \mathbf{V}^{\mathsf{T}}\operatorname{clr}(\boldsymbol{p}),$

which is well-defined because $\mathbf{V}^{\mathsf{T}}\mathbf{1}=0$. Requires a strictly positive composition (open simplex).

Usage

isometric_logistic(z, V = NULL)

isometric_log_ratio(p, V = NULL)

Arguments

z	Numeric vector $\boldsymbol{z}\in\mathbb{R}^{J-1}$.
V	Optional ILR basis; see isometric_logistic().
p	Numeric vector on the open simplex.

Value

Numeric vector $\boldsymbol{p}$ on the unit simplex.

Numeric vector $\boldsymbol{z}\in\mathbb{R}^{J-1}$.

See Also

isometric_log_ratio(), helmert_basis(), additive_logistic()

Examples

z <- c(0.2, -0.5)
p <- isometric_logistic(z)
sum(p)
isometric_log_ratio(p)
isometric_logistic(c(800, 1200))

---