This article is the bivariate Gaussian-convolution toy of the methods paper (G = 2, J = 2). It isolates gene–gene correlation from mean separation. The factorial grid and ADEMP pipeline live in a single script, scripts/fig02_bivariate_toy.R. Scenario builders (build_bivariate_scenario_config(), bivariate_toy_deconvolution_functions()) sit at the top of that file. Sourcing the script from a vignette or from a temp copy in tests defines those functions without launching the 972-scenario pipeline.

run_simulation_benchmark() wraps simulate_bulk_mixture(), deconvolute_ratios(), and compute_benchmark_metrics(). Scenario rows are sequential; sample-level workers live only in deconvolute_ratios(). Reporting conventions (ADEMP, Nature Methods, raincloud / forest plots) are in how to build synthetic scenarios.

%%{init: {"theme": "sandstone"}}%%
flowchart TD
  P["3 compositions<br/>Ba / Mo / Hi"] --> S["972 scenarios"]
  V["2 variance structures<br/>Ho homoscedastic / He heteroscedastic"] --> S
  C["2 CLD mean separations<br/>Sm small / Lg large"] --> S
  R["9 x 9 gene-gene correlations<br/>rho in -0.8 to 0.8, step 0.2"] --> S
  S --> A["7 deconvolution solvers"]

%%{init: {"theme": "sandstone"}}%%
flowchart TD
  P["3 compositions<br/>Ba / Mo / Hi"] --> S["972 scenarios"]
  V["2 variance structures<br/>Ho homoscedastic / He heteroscedastic"] --> S
  C["2 CLD mean separations<br/>Sm small / Lg large"] --> S
  R["9 x 9 gene-gene correlations<br/>rho in -0.8 to 0.8, step 0.2"] --> S
  S --> A["7 deconvolution solvers"]

---
config:
  layout: elk
  theme: sandstone
---
flowchart TB
  ROOT["Bivariate toy"]
  ROOT --> Sm["Small CLD"]
  ROOT --> Lg["Large CLD"]
  Sm --> SmHo["Homoscedastic"]
  Sm --> SmHe["Heteroscedastic"]
  Lg --> LgHo["Homoscedastic"]
  Lg --> LgHe["Heteroscedastic"]
  SmHo --> SmHoBa["Balanced"]
  SmHo --> SmHoMo["Moderately unbalanced"]
  SmHo --> SmHoHi["Highly unbalanced"]
  SmHe --> SmHeBa["Balanced"]
  SmHe --> SmHeMo["Moderately unbalanced"]
  SmHe --> SmHeHi["Highly unbalanced"]
  LgHo --> LgHoBa["Balanced"]
  LgHo --> LgHoMo["Moderately unbalanced"]
  LgHo --> LgHoHi["Highly unbalanced"]
  LgHe --> LgHeBa["Balanced"]
  LgHe --> LgHeMo["Moderately unbalanced"]
  LgHe --> LgHeHi["Highly unbalanced"]

---
config:
  layout: elk
  theme: sandstone
---
flowchart TB
  ROOT["Bivariate toy"]
  ROOT --> Sm["Small CLD"]
  ROOT --> Lg["Large CLD"]
  Sm --> SmHo["Homoscedastic"]
  Sm --> SmHe["Heteroscedastic"]
  Lg --> LgHo["Homoscedastic"]
  Lg --> LgHe["Heteroscedastic"]
  SmHo --> SmHoBa["Balanced"]
  SmHo --> SmHoMo["Moderately unbalanced"]
  SmHo --> SmHoHi["Highly unbalanced"]
  SmHe --> SmHeBa["Balanced"]
  SmHe --> SmHeMo["Moderately unbalanced"]
  SmHe --> SmHeHi["Highly unbalanced"]
  LgHo --> LgHoBa["Balanced"]
  LgHo --> LgHoMo["Moderately unbalanced"]
  LgHo --> LgHoHi["Highly unbalanced"]
  LgHe --> LgHeBa["Balanced"]
  LgHe --> LgHeMo["Moderately unbalanced"]
  LgHe --> LgHeHi["Highly unbalanced"]

Generative model

With p₁ + p₂ = 1, only one free ILR coordinate is estimated (see derivatives under simplex transforms). Even when mean profiles are similar, gene–gene correlation can degrade mean-only solvers, and is partly recovered once Σ(p) = ∑_jp_j²Σ_j enters the likelihood.

Factorial design (972 scenarios)

N = 500 Monte Carlo replicates per scenario.

CLD is centroid (mean-profile) Euclidean separation ‖μ_⋅1 − μ_⋅2‖₂, not an Aitchison distance; see the descriptor notes. Small CLD uses μ_⋅1 = (20, 22), μ_⋅2 = (22, 20) ($d=2\sqrt{2}\approx 2.83$, cosine  ≈ 0.995). Large CLD uses (20, 40) and (40, 20) ($d=20\sqrt{2}\approx 28.3$, cosine 0.8).

Scenario ID

Each row is labelled B{index}_{variance}_{composition}_{CLD}. The index is the factorial row number (1–972). Example: B1_Ho_Ba_Sm is the first row, homoscedastic, balanced proportions, small CLD.

The slim design grid is output/fig02/bivariate_config.rds (no true_theta). Geometry, MixSim BarOmega, Hellinger, and Jeffreys live in bivariate_descriptors.rds. Convolution parameters are bivariate_theta.rds. ADEMP metrics (regression, monte_carlo, optimisation, call) are bivariate_benchmark.rds. All four files share the ID key. read_simulation_artefacts(dir, "bivariate", assemble = TRUE) rejoins them for plotting.

The Monte Carlo $\hat{\boldsymbol{p}}$ draws themselves are the cell-type columns of optimisation (one row per replicate × algorithm × scenario). pivot_mc_estimates() stacks those columns into long form. Interior Wald intervals for convolution-likelihood solvers use vcov_ilr_delta() (expected Fisher information through the ILR chart), matching confint.decovart_fit(). Coverage, mean_model_se, and se_sd_ratio (mean model SE / empirical SD) are then finite. Mean-only solvers (NNLS, LSEI) leave those Wald columns missing: they do not maximise the convolution likelihood, so an ILR interval would not apply.

library(DeCovarT)
source("scripts/fig02_bivariate_toy.R")

bivariate_config <- build_bivariate_scenario_config()
nrow(bivariate_config)

Log-likelihood along the simplex

For J = 2 the free coordinate is p₁ ∈ (0, 1). The listing below evaluates ℓ_y ∣ ζ(p₁, 1 − p₁) on one bulk draw. It is not executed during R CMD build: Quarto’s Windows CLI starts a new R process that cannot see the temporary library (quarto-r#217). The same likelihood path is checked in tests/testthat/.

set.seed(20260828)
rho <- 0.6
R <- matrix(c(1, rho, rho, 1), nrow = 2)
Sigma <- array(c(R, R), dim = c(2, 2, 2))
dimnames(Sigma) <- list(
  rownames(mu_small),
  rownames(mu_small),
  colnames(mu_small)
)
y_one <- DeCovarT::simulate_bulk_mixture(
  mu_small,
  Sigma,
  p = p_bal,
  n = 1
)$Y[, 1]
p1_grid <- seq(0.02, 0.98, by = 0.01)
ll_grid <- vapply(
  p1_grid,
  function(p1) {
    DeCovarT::loglik_multivariate(c(p1, 1 - p1), y_one, mu_small, Sigma)
  },
  numeric(1)
)
ll_df <- data.frame(p1 = p1_grid, loglik = ll_grid)
ggplot2::ggplot(ll_df, ggplot2::aes(x = p1, y = loglik)) +
  ggplot2::geom_line(colour = "#1B4F72", linewidth = 0.8) +
  ggplot2::geom_vline(
    xintercept = 0.5,
    linetype = "dashed",
    colour = "grey30"
  ) +
  ggplot2::labs(
    x = "p1 (p2 = 1 - p1)",
    y = "Log-likelihood"
  ) +
  ggplot2::theme_bw(base_size = 11)

Inference

Seven solvers are evaluated on each scenario:

Performance is summarised using the ADEMP framework (Morris et al. 2019): bias, RMSE, coverage of ILR Wald intervals (Wilson intervals on the coverage rate), SE/SD, and optimiser failure rate. Mean-only baselines have no convolution Wald SE.

N_REPLICATES <- as.integer(Sys.getenv("N_REPLICATES", unset = "500"))
deconvolution_functions <- bivariate_toy_deconvolution_functions(
  itmax = 200L,
  epsilon = 1e-4
)
bivariate <- run_simulation_benchmark(
  scenario_config = bivariate_config,
  deconvolution_functions = deconvolution_functions,
  n = N_REPLICATES,
  cores = 1L
)

config on the in-memory benchmark still carries Shannon entropy of p and MixSim overlap for tests. Saved artefacts split that information: optimisation stores per-sample $\hat{\boldsymbol{p}}$, elapsed time, and memory; regression and monte_carlo are the composition and ADEMP blocks from compute_benchmark_metrics().

Visualisations

Figures are written from the split RDS (no ADEMP refit). Density books go to output/fig02/density_visualisations/; RMSE / MAE / Aitchison tiles, raincloud, quantile-dot raincloud, forest, KKT raincloud, stacked convergence, whitened-ILR Q-Q / χ² Q-Q / KS lollipops, similarity, solver-dot, runtime, and memory books go to output/fig02/performance_visualisations/. Tile heatmaps (ggplot2::geom_tile()) show RMSE, MAE, and Aitchison distance on the (ρ₁, ρ₂) plane: one PDF per metric, one page per meta-scenario (CLD () variance () composition), seven solver panels. Scenario factors are relevelled (small then large CLD; homoscedastic then heteroscedastic; balanced, moderately unbalanced, highly unbalanced) when artefacts are read, and solvers follow nnls, lsei, SA, gradient, LBFGS, Newton-Raphson, Marquardt-Levenberg. Density and performance books share the four correlation corners (ρ₁, ρ₂) ∈ {(0, 0), (−0.8, −0.8), (0.8, 0.8), (−0.8, 0.8)} (2 × 2 × 3 = 12 pages). Similarity heatmaps compute Pearson (r) of paired () within each corner, with Ward D2 clustering of (1-r) and Ward merge heights on the right dendrogram. Wald forests drop NNLS, LSEI, and SA; solid whiskers are (1.96) times the expected-Fisher Wald SE at the true composition, dashed whiskers use the empirical Monte Carlo SD, and coloured geom_label boxes (one per cell type, matching fig03) report RMSE and coverage. The native-space log-likelihood surface uses a log10 colour scale of the relative likelihood (L/L) so boundary peaks are visible. The matching rgl book uses the same log10 relative-likelihood z-axis. Rainclouds plot (p) with solver fill and stroke and a bottom legend. A companion dotsinterval.pdf book replaces the bounded [0, 1] KDE slab with quantile dots. Runtime and memory rainclouds use one page per meta-scenario, the four correlation corners on the x-axis, and a log10 y-axis. Solver-dot pages use the same 9-by-9 ((_1,_2)) grid as the RMSE heatmaps: colour is mean RMSE and size is mean Aitchison distance, with horizontal y-axis labels. ggplot data tables are saved under output/fig02/ggplot_rds/.

Likelihood geometry and solver behaviour

The factorial heatmaps remain the right design summary: large CLD and ρ_j = 0 make the mean map injective, so every solver is accurate; small CLD plus aligned correlations flatten the mean residual and raise RMSE off the diagonal ρ₁ = ρ₂. That saddle is a statement about identifiability of μp, not about which numerical method is used. The pack in output/fig02/mle_explanation/ (scenario descriptors, ILR profiles evaluated at y = μp^⋆, expected-Fisher Wald standard errors, and Newton versus Marquardt forests) is what separates three mechanisms that the heatmaps mix: the shape of ℓ(p; y), the surrogate each solver actually minimises, and the chart / stopping rule that walks that surrogate.

What is being optimised

Mean-only NNLS and DeconRNASeq (lsei in limSolve) minimise a convex quadratic (y − μp)^⊤(y − μp) on the simplex (Dessole et al. 2023; Gong and Szustakowski 2013; Soetaert et al. 2026). For a convex objective a local minimiser is global (Boyd et al. 2004, Ch. 4), so those fits are unique solutions of the least-squares problem. They are not maximisers of the convolution log-likelihood: Σ(p) never enters the quadratic. When μ_⋅1 and μ_⋅2 are close, that quadratic is a ridge along the simplex. A strictly feasible QP with bound constraints then has vertices as KKT candidates, which is why mean-only fits pile up at (1, 0) and (0, 1) whenever the residual is uninformative. Simulated annealing explores the same non-concave ℓ as DeCovarT but without curvature, so it is a robustness check, not an MLE.

The convolution solvers maximise the Gaussian-convolution log-likelihood of the MLE properties article (that article is pkgdown-only; this vignette is the CRAN-shipped companion). None of them inherits Boyd’s global certificate, because that likelihood is not concave. They differ in coordinates and damping:

• L-BFGS-B (stats::optim) runs limited-memory BFGS with box constraints 0 ≤ p_j ≤ 1 (Byrd et al. 1995; Zhu et al. 1997), then closes the simplex by p/∑p. The equality 1^⊤p = 1 is not a constraint of the line search, so Σ(p) can approach singularity mid-iteration; the implementation guards that with a finite penalty. Faces of the cube remain admissible, which is how a limited-memory quasi-Newton method sits on a simplex face after renormalisation.
• Newton-Raphson (stats::nlminb) and Marquardt-Levenberg (marqLevAlg; (Philipps et al. 2023, 2021)) work in ILR coordinates, so the walk stays in the relative interior until a coordinate diverges. Newton uses the analytic Hessian. Marquardt adds a Levenberg damping term so that the working matrix stays positive definite (Marquardt 1963) even when the Hessian of ℓ is indefinite, and stops on the relative-distance-to-maximum (RDM) criterion: the ratio of optimisation error to statistical error (Commenges et al. 2006). On a plateau, RDM is a feature — further progress along a flat ridge does not change inferential accuracy — and a bias: the iterate can halt nearer the barycentre than the true face.

The gradient solver is ILR BFGS without a Hessian, so it inherits the chart but not Newton’s local quadratic rate. Boyd’s Newton analysis (quadratic convergence in a neighbourhood of a strict local maximiser with positive-definite Hessian) applies to nlminb only when the realised ℓ is locally strongly concave. The bivariate small-CLD profiles are not.

One bulk, N = 1

Every Monte Carlo replicate here is a single G = 2 bulk. The MLE is therefore the N = 1 estimator of the finite-sample section: biased at O(1), with Wald SEs that describe the Monte Carlo spread of that one-observation maximiser, not a sampling experiment with growing N. Pooling N i.i.d. columns with a shared p^⋆ is what would make $\hat{\boldsymbol{p}}_N$ consistent. The ILR profiles in the pack are evaluated at the conditional mean y = μp^⋆, which isolates the log-determinant remainder; the forests average over random Y.

Reading the twelve meta-scenarios

Numbers below are from forest_newton_marquardt.csv, theoretical_wald_se.csv, and the range of ilr_loglik_profiles.csv (max minus min of ℓ(ρ) on the plotting grid).

Large CLD, balanced, ρ_j = 0 (B567). The mean map is well conditioned. The ILR profile spans about 101 nats and SE_ρ = 0.071. Newton and Marquardt agree to numerical noise: bias of p₁ about 0, RMSE 0.024, coverage 0.96, and the empirical SD matches the expected-Fisher SE (0.025). This is the textbook interior MLE: strongly concave locally, RDM and Newton stop at the same point.

Large CLD, balanced, ρ_j = −0.8 (B487). Equal negative correlation inflates variance along the anti-diagonal of gene space but does not destroy the mean contrast. Profile range 56, SE_ρ = 0.095, RMSE 0.033 for both second-order solvers. Covariance-aware curvature is still enough.

Small CLD, balanced (B1, B81). Means (20, 22) versus (22, 20) almost collinear. Profile range drops to 1.2–1.7 nats; SE_ρ is 0.95 (both ρ = −0.8) or 0.71 (ρ = 0). RMSE of p₁ is 0.19–0.23 and Wald coverage is conservative (0.98–1). Marquardt is slightly less variable than Newton here (RMSE 0.192 versus 0.234 on B1): damping and RDM cut walks along the ridge, which for a balanced truth is a small bias toward 1/2 that happens to reduce RMSE. Mean-only QPs see an even flatter residual and are the ones that snap to vertices.

Highly unbalanced, small CLD (B325, B405, p^⋆ = (0.99, 0.01)). The ILR chart sends the rare type to ρ → +∞. Fisher curvature in ρ collapses: SE_ρ is 24 (both ρ = −0.8) or 20 (ρ = 0), while the profile range on a finite grid is only 2.4–4.1 nats — a plateau, not a peak. Newton bias for p₁ is −0.22 and −0.18; Marquardt bias is −0.40 and −0.38. Both shrink toward the barycentre, Marquardt more so, exactly as RDM predicts on a flat tail. Coverage stays near 0.92–0.94 because the Wald interval is already huge. This is the geometry behind chi-bar-square calibration on a face (boundary LRTs): the tangent cone is a half-line, Wilks does not apply, and a Wald SE in ILR is not a calibrated rare-type interval.

Highly unbalanced, large CLD (B811, B891). The profile in ρ is not flat in the raw log-likelihood (range 218–392), because a 1% contamination still shifts a well-separated mean. The ILR SE remains large (3.3 and 2.5) because a small p₂ is a long way out on the chart. Newton follows that peak: bias −0.017 and −0.010, RMSE 0.035 and 0.026, coverage 0.97. Marquardt does not: bias −0.26 and −0.29, RMSE 0.35 and 0.37, coverage 0.47 and 0.40. The damping / RDM pair stops in ILR before |ρ| is large enough to represent p₂ = 0.01. A sharp peak in p can still look like a long walk in ρ; Marquardt’s statistical stopping rule then under-shoots a rare type even when the convolution likelihood is informative. That failure is an optimiser-chart interaction, not a lack of covariance signal.

Moderately unbalanced rows sit between these extremes: small-CLD RMSE 0.23–0.31 with Marquardt again more barycentric; large-CLD RMSE about 0.03–0.04 with the two second-order solvers aligned.

What the heatmaps still add

The 9 × 9 (ρ₁, ρ₂) tiles, which are not in the twelve-row explanation pack, remain the place to see mean-only methods lose first as |ρ| grows at small CLD, and to see DeCovarT retain the diagonal saddle. The pack’s message is narrower and sharper: with N = 1, a plateau in ILR makes every interior solver noisy; Marquardt’s RDM is the right stop on a balanced ridge and the wrong stop on a rare type; Newton with a full Hessian tracks p^⋆ once the means separate; convex QPs solve a different problem and are entitled only to that problem’s unique minimiser.

See also

• Variance-driven hybrid (G = 20, J = 3): §2.2
• Moment generator and ADEMP reporting: How to build synthetic scenarios
• Regular-case MLE checks: Appendix S1

References