Package 'DeCovarT'

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

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)

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,
  mean_signature_matrix,
  Sigma,
  p = NULL,
  null_value = NULL,
  n_boot = 500L,
  level = 0.95,
  n_replicates = NULL,
  epsilon = 10^-8,
  itmax = 500L
)

Arguments

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), n_failed (refits that errored, left as NA columns), 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: bootstrap_single_cell(), boundary_diagnostics(), chi_bar_square_pvalue(), confint_decovart(), confint_ilr_backtransform(), confint_profile_decovart(), confint_profile_single_cell(), confint_single_cell(), equivariance_check_decovart(), fit_firth_decovart(), fit_firth_single_cell(), lrt_decovart(), lrt_single_cell(), multistart_decovart(), reference_bootstrap_decovart(), reference_bootstrap_single_cell(), restricted_mle_decovart(), restricted_mle_single_cell()

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

Description

Draws $\boldsymbol{Y}^{(b)}\sim\mathcal{N}(\hat c\tilde{\boldsymbol{\mu}} \hat{\boldsymbol{p}},\hat\phi\boldsymbol{S}_W(\hat{\boldsymbol{p}}))$ and refits the profile model. Percentile intervals are returned for the composition and for the two scales. The reference moments stay fixed.

Usage

bootstrap_single_cell(
  y,
  mean_signature_matrix,
  Sigma,
  n_boot = 200L,
  level = 0.95,
  epsilon = 1e-06,
  itmax = 40L
)

Arguments

Value

A list with estimates, interval, n_failed (refits that errored, left as NA columns), scales, scale_interval and fit.

See Also

confint_single_cell(), bootstrap_decovart()

Other decovart_inference: bootstrap_decovart(), boundary_diagnostics(), chi_bar_square_pvalue(), confint_decovart(), confint_ilr_backtransform(), confint_profile_decovart(), confint_profile_single_cell(), confint_single_cell(), equivariance_check_decovart(), fit_firth_decovart(), fit_firth_single_cell(), lrt_decovart(), lrt_single_cell(), multistart_decovart(), reference_bootstrap_decovart(), reference_bootstrap_single_cell(), restricted_mle_decovart(), restricted_mle_single_cell()

Examples

mu <- matrix(c(2, 0.2, 0.2, 2, 0.2, 3, 0.2, 3), nrow = 4)
colnames(mu) <- c("ct1", "ct2")
sigma <- array(rep(diag(0.05, 4), 2), c(4, 4, 2))
y <- 20 * drop(mu %*% c(0.7, 0.3)) + c(0.3, -0.2, 0.1, -0.1)
set.seed(1)
bootstrap_single_cell(y, mu, sigma, n_boot = 5)$interval

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

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(), bootstrap_single_cell(), boundary_diagnostics(), confint_decovart(), confint_ilr_backtransform(), confint_profile_decovart(), confint_profile_single_cell(), confint_single_cell(), equivariance_check_decovart(), fit_firth_decovart(), fit_firth_single_cell(), lrt_decovart(), lrt_single_cell(), multistart_decovart(), reference_bootstrap_decovart(), reference_bootstrap_single_cell(), restricted_mle_decovart(), restricted_mle_single_cell()

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)

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

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)

Description

Returns the $p_j p_{\ell}$-weighted mean of pairwise Jeffreys (symmetrised KL) divergences between purified Gaussians, with the closed form of the multivariate-normal Kullback–Leibler divergence (https://statproofbook.github.io/P/mvn-kl.html)

$\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

Value

Scalar average pairwise Jeffreys divergence.

See Also

compute_average_overlap(), compute_average_riemannian()

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)

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

Value

Scalar average pairwise overlap (MixSim BarOmega).

See Also

overlap_gaussian_mc(), compute_average_jeffreys()

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)

Description

Average over $j<\ell$ of the affine-invariant Riemannian metric on the cone of symmetric positive-definite matrices, $d_R(\Sigma_j,\Sigma_\ell)=\lVert\log(\Sigma_j^{-1/2}\Sigma_\ell \Sigma_j^{-1/2})\rVert_F$, computed from the eigenvalues of $R^{-\mathsf{T}}\Sigma_\ell R^{-1}$ with $\Sigma_j=R^{\mathsf{T}}R$. The Frobenius chord $\lVert\Sigma_j-\Sigma_\ell\rVert_F$ is not used: it leaves the cone, is not inversion-invariant, and suffers the swelling effect. With $J=2$ the average is the pairwise distance.

Usage

compute_average_riemannian(true_theta, J = NULL)

Arguments

Value

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

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)

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,
  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

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))

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

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

Description

One entry point for the interior intervals of the population convolution (covariance weights $p_j^2$, no technical scales). The default, method = "ilr", is confint_ilr_backtransform().

ilr

Wald interval for $\boldsymbol{\rho}$, mapped by the inverse ILR chart.

wald

Delta method on $\boldsymbol{p}$. Endpoints can leave $[0,1]$.

firth

Wald interval at the Jeffreys-penalised estimate (fit_firth_decovart()). interval_ilr is the same estimate mapped through the inverse chart.

profile

Likelihood-ratio interval (confint_profile_decovart()). This is the interior calibration. A null on a face is lrt_decovart(), not an interval.

bootstrap

Parametric percentile interval (bootstrap_decovart()).

reference

Resample the reference, then refit (reference_bootstrap_decovart()). Two or more donors in every cell type use the donor bootstrap and keep the bulk fixed. Otherwise columns are resampled within each cell type and the bulk is regenerated. Cell-type labels are not shuffled; that check is equivariance_check_decovart().

Columns of y share one composition. fit_decovart() instead estimates a separate composition per column, and confint.decovart_fit() is the per-column delta method.

Usage

confint_decovart(
  y,
  mean_signature_matrix = NULL,
  Sigma = NULL,
  method = c("ilr", "wald", "firth", "profile", "bootstrap", "reference"),
  level = 0.95,
  p = NULL,
  n_draws = 2000L,
  n_boot = 500L,
  seed = NULL,
  reference_profiles = NULL,
  donor_ids = NULL,
  resample = NULL,
  regenerate_bulk = NULL,
  ...
)

Arguments

Value

A list with method, estimate and interval. Other elements depend on method.

See Also

lrt_decovart(), confint_single_cell()

Other decovart_inference: bootstrap_decovart(), bootstrap_single_cell(), boundary_diagnostics(), chi_bar_square_pvalue(), confint_ilr_backtransform(), confint_profile_decovart(), confint_profile_single_cell(), confint_single_cell(), equivariance_check_decovart(), fit_firth_decovart(), fit_firth_single_cell(), lrt_decovart(), lrt_single_cell(), multistart_decovart(), reference_bootstrap_decovart(), reference_bootstrap_single_cell(), restricted_mle_decovart(), restricted_mle_single_cell()

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_decovart(y, mu, Sigma)$interval
confint_decovart(y, mu, Sigma, method = "wald")$interval

Description

Forms a Wald interval for the unconstrained ILR coordinate $\boldsymbol{\rho}$ from the expected Fisher information, then maps it to the simplex by the inverse chart $\boldsymbol{\psi}$. For two cell types the map is monotone, so the endpoints of the interval for $\rho$ are transformed. For more than two cell types the joint normal draw $\boldsymbol{\rho}^{(b)}$ is transformed and the interval is the marginal quantile of $\boldsymbol{\psi}(\boldsymbol{\rho}^{(b)})$. Every transformed draw is a composition, so the interval lies in $(0,1)$.

This is the default of confint_decovart(). The symmetric delta-method interval on $\boldsymbol{p}$ remains confint.decovart_fit(). Columns of y are replicates that share one composition: their number multiplies the Fisher information.

Usage

confint_ilr_backtransform(
  y,
  mean_signature_matrix,
  Sigma,
  p = NULL,
  level = 0.95,
  n_draws = 2000L,
  seed = NULL,
  epsilon = 10^-8,
  itmax = 500L
)

Arguments

Value

A list with estimate, interval on the simplex, rho, rho_interval (component-wise Wald bounds for $\boldsymbol{\rho}$), vcov_rho, vcov, mapping ("endpoint" or "monte_carlo") and n_obs. On a simplex face the interval is NA.

See Also

confint_decovart(), vcov_ilr_delta(), confint.decovart_fit()

Other decovart_inference: bootstrap_decovart(), bootstrap_single_cell(), boundary_diagnostics(), chi_bar_square_pvalue(), confint_decovart(), confint_profile_decovart(), confint_profile_single_cell(), confint_single_cell(), equivariance_check_decovart(), fit_firth_decovart(), fit_firth_single_cell(), lrt_decovart(), lrt_single_cell(), multistart_decovart(), reference_bootstrap_decovart(), reference_bootstrap_single_cell(), restricted_mle_decovart(), restricted_mle_single_cell()

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_ilr_backtransform(y, mu, Sigma)$interval

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,
  mean_signature_matrix,
  Sigma,
  level = 0.95,
  celltypes = colnames(mean_signature_matrix),
  epsilon = 10^-8,
  itmax = 500L
)

Arguments

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(), bootstrap_single_cell(), boundary_diagnostics(), chi_bar_square_pvalue(), confint_decovart(), confint_ilr_backtransform(), confint_profile_single_cell(), confint_single_cell(), equivariance_check_decovart(), fit_firth_decovart(), fit_firth_single_cell(), lrt_decovart(), lrt_single_cell(), multistart_decovart(), reference_bootstrap_decovart(), reference_bootstrap_single_cell(), restricted_mle_decovart(), restricted_mle_single_cell()

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)

Description

Inverts the likelihood-ratio test for one coordinate of the profile log-likelihood. The two scales are concentrated out at every trial value. An endpoint that lands on a face is not a Wilks interval; use lrt_single_cell() for a boundary null.

Usage

confint_profile_single_cell(
  y,
  mean_signature_matrix,
  Sigma,
  level = 0.95,
  celltypes = colnames(mean_signature_matrix),
  variance = c("ml", "reml"),
  epsilon = 1e-08,
  itmax = 100L
)

Arguments

Value

A matrix with columns estimate, lower and upper.

See Also

confint_single_cell(), lrt_single_cell()

Other decovart_inference: bootstrap_decovart(), bootstrap_single_cell(), boundary_diagnostics(), chi_bar_square_pvalue(), confint_decovart(), confint_ilr_backtransform(), confint_profile_decovart(), confint_single_cell(), equivariance_check_decovart(), fit_firth_decovart(), fit_firth_single_cell(), lrt_decovart(), lrt_single_cell(), multistart_decovart(), reference_bootstrap_decovart(), reference_bootstrap_single_cell(), restricted_mle_decovart(), restricted_mle_single_cell()

Examples

mu <- matrix(c(2, 0.2, 0.2, 2, 0.2, 3, 0.2, 3), nrow = 4)
colnames(mu) <- c("ct1", "ct2")
sigma <- array(rep(diag(0.05, 4), 2), c(4, 4, 2))
y <- 20 * drop(mu %*% c(0.7, 0.3)) + c(0.3, -0.2, 0.1, -0.1)
confint_profile_single_cell(y, mu, sigma)

Description

The composition and the two technical scales are separated by the Schur complement of the observed Hessian. method chooses the interval for the composition. The default is the inverse ILR chart applied to a Wald interval for $\boldsymbol{\rho}$, using the profile information $\{-\nabla^{2}_{\boldsymbol{\rho}}\ell^{\ast}\}^{-1}$.

Intervals for $(c,\phi)$ in nuisance$interval are Wald intervals on the log scale, mapped back by the exponential, from the Schur complement that removes $\boldsymbol{\rho}$. nuisance$interval_wald is the symmetric normal interval on the original scale. Bootstrap methods replace both with percentile intervals.

ilr

Default. Inverse ILR chart.

wald

Delta method on the composition. Endpoints can leave $[0,1]$.

firth

Wald interval at the penalised profile estimate. interval_ilr maps that estimate through the inverse chart.

profile

Interior likelihood-ratio interval. A face null is lrt_single_cell().

bootstrap

Parametric percentile interval.

reference

Donor bootstrap, or a within-type cell resample that regenerates the bulk. See reference_bootstrap_single_cell().

These intervals treat the reference covariance as known. They do not add a sandwich or higher-order correction for that plug-in.

variance = "reml" corrects for the $\nu=G-J$ residual degrees of freedom of the variance scale: covariances are inflated by $G/\nu$, Wald and ILR quantiles come from $t_{\nu}$ (a multivariate $t$ for the ILR map when $J>2$), $\phi$ gets the chi-square interval $[G\hat\phi/\chi^2_{\nu;1-\alpha/2}, G\hat\phi/\chi^2_{\nu;\alpha/2}]$, and the profile interval uses the cut-off of confint_profile_single_cell(). The maximum-likelihood $\hat\phi$ is biased by about $\nu/G$, so with few genes ($G=10$, $J=2$) the "ml" intervals cover about 0.85 to 0.92 at the 0.95 level; the bootstrap methods are unaffected by variance.

Usage

confint_single_cell(
  y,
  mean_signature_matrix = NULL,
  Sigma = NULL,
  method = c("ilr", "wald", "firth", "profile", "bootstrap", "reference"),
  level = 0.95,
  variance = c("ml", "reml"),
  n_draws = 2000L,
  n_boot = 200L,
  seed = NULL,
  reference_profiles = NULL,
  donor_ids = NULL,
  resample = NULL,
  regenerate_bulk = NULL,
  ...
)

Arguments

Value

A list with method, estimate (the fitted composition), interval and nuisance. Bootstrap results also include the replicate matrices.

See Also

vcov_single_cell(), confint_decovart(), lrt_single_cell()

Other decovart_inference: bootstrap_decovart(), bootstrap_single_cell(), boundary_diagnostics(), chi_bar_square_pvalue(), confint_decovart(), confint_ilr_backtransform(), confint_profile_decovart(), confint_profile_single_cell(), equivariance_check_decovart(), fit_firth_decovart(), fit_firth_single_cell(), lrt_decovart(), lrt_single_cell(), multistart_decovart(), reference_bootstrap_decovart(), reference_bootstrap_single_cell(), restricted_mle_decovart(), restricted_mle_single_cell()

Examples

mu <- matrix(c(2, 0.2, 0.2, 2, 0.2, 3, 0.2, 3), nrow = 4)
colnames(mu) <- c("ct1", "ct2")
sigma <- array(rep(diag(0.05, 4), 2), c(4, 4, 2))
y <- 20 * drop(mu %*% c(0.7, 0.3)) + c(0.3, -0.2, 0.1, -0.1)
fit <- confint_single_cell(y, mu, sigma)
fit$interval
fit$nuisance$interval

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

Value

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

Examples

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

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

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
)

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

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 on ILR coordinates $\boldsymbol{z}$ with the analytic score and observed Hessian, safeguarded by the PORT trust region of 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)

Description

Fits one bulk sample from the linear single-cell mixture $\boldsymbol{y}\sim\mathcal{N}(c\tilde{\boldsymbol{\mu}}\boldsymbol{p}, \phi\boldsymbol{S}_W(\boldsymbol{p}))$.

The start is single_cell_nuisance_scales(): nonnegative least squares with a free scale (DeconRNASeq's quadratic programme without the simplex equality), then the closed-form scales at that composition.

With profile = TRUE (default) the two scales are profiled out: for every trial composition they are recomputed in closed form, so the weight matrix $\boldsymbol{S}_W(\boldsymbol{p})^{-1}$ is rebuilt at every iterate, and Newton-Raphson runs on the profile log-likelihood of the $J-1$ ILR coordinates (profile_loglik_single_cell_ilr()). Its limit is the joint maximum-likelihood estimate of $(\boldsymbol{p},c,\phi)$, reached in one loop.

With profile = FALSE the scales stay at their start values (a two-stage plug-in fit), which is not the joint maximum. Supplying both mean_scale and variance_scale fixes them at those values.

Every fit uses the analytic score and observed Hessian with the PORT trust region of stats::nlminb(): the step maximises the quadratic model inside a ball whose radius grows or shrinks with the ratio of actual to predicted increase; there is no line search and no fixed step length. The returned vcov is the delta-method image of the inverse negative profile Hessian; intervals are confint_single_cell().

Usage

deconvolute_single_cell(
  y,
  mean_signature_matrix,
  Sigma,
  epsilon = 1e-08,
  itmax = 100L,
  profile = TRUE,
  initial_p = NULL,
  mean_scale = NULL,
  variance_scale = NULL,
  estimand = c("cellular", "transcript"),
  cell_size = NULL
)

Arguments

Value

A list with composition (the fitted simplex coordinate), p, q, rho, estimand, fit_type ("profile", "two_stage" or "fixed_scales"), mean_scale, variance_scale, platform_scale, n_cells, cell_size, ols_scale, loglik, convergence (0 from stats::nlminb() means success), iterations, p_start, message, vcov_rho, vcov and vcov_nuisance (see vcov_single_cell()). A singular information matrix is returned as NA without failing the fit.

See Also

single_cell_nuisance_scales(), profile_loglik_single_cell_ilr(), deconvolute_ratios_Newton_Raphson(), rescale_single_cell_moments()

Examples

genes <- paste0("g", 1:4)
cts <- paste0("ct", 1:2)
mu <- matrix(
  c(2, 2, 0.2, 0.2, 0.2, 0.2, 3, 3),
  nrow = 4,
  dimnames = list(genes, cts)
)
sigma <- array(0, c(4, 4, 2), dimnames = list(genes, genes, cts))
sigma[,, 1] <- diag(0.05, 4)
sigma[,, 2] <- diag(0.05, 4)
set.seed(1)
sim <- simulate_bulk_from_cells(
  mu,
  sigma,
  p = c(0.8, 0.2),
  n_cells = 2000,
  platform_scale = 1
)
fit <- deconvolute_single_cell(sim$Y[, 1], mu, sigma)
round(fit$p, 3)

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

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

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

Value

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

See Also

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

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

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("Newton-Raphson", "Marquardt-Levenberg", "L-BFGS-B"),
  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

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: