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.

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

Every result of the run is in output/fig02/fig02_benchmark.rds, a list keyed by the scenario ID: meta (seed, N, B, date, git commit, session), config (the design grid without the true moments), theta_true (convolution parameters per ID), descriptors (geometry, MixSim BarOmega, Hellinger, Jeffreys), regression and monte_carlo (composition and ADEMP blocks of compute_benchmark_metrics()), optimisation (one row per replicate × algorithm × scenario, with the estimates, convergence, KKT residual, elapsed time and memory), intervals (coverage and width of every interval construction) and data (the simulated bulks and their targets).

Mean-only lsei and NNLS do not maximise the convolution likelihood, so no interval is computed for them.

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

Five methods are evaluated on each scenario:

• Dropped from the figure. Marquardt–Levenberg, projected gradient and simulated annealing remain in the package but leave the benchmark. CIBERSORT is not run at G = 2: its ν-SVR regresses the G bulk values on the J signatures, so two genes give two observations for two coefficients and an intercept. The second-generation methods (MEAD, RNA-Sieve, ReDeconv, DECALS) are left to the covariance-driven scenario, since a two-gene reference cannot carry their reference estimation.
• NNLS followed by renormalisation is least squares on the simplex with a free scale (β = cp, c = 1^⊤β); lsei fixes c = 1. The two agree only when ĉ = 1. On 500 bulks of this design the largest coordinate difference has median 0.002 and 95th percentile 0.015, and the fitted scale has standard deviation 0.03 around 1: close, but not identical. Both are run: lsei keeps the sum-to-one equality and the bounds at every iteration, NNLS renormalises a free-scale fit. A robust linear model is not run on two genes.

The interval constructions are those of population inference: wald (delta method on p), ilr (Wald in ILR coordinates mapped back), profile (likelihood-ratio inversion) and oracle (p̂_j plus or minus the Wald half-width at p^⋆). The parametric bootstrap (bootstrap_decovart()) is not run here: the moments of a scenario are fixed, so coverage is measured over its N bulks, and N × B refits per solver dominated the run time. Performance is summarised with the ADEMP framework (Morris et al. 2019): bias, RMSE, coverage of each interval construction (Wilson intervals on the coverage rate), width, SE/SD, and optimiser failure rate.

The run size is hard-coded at the top of the script: N_REPLICATES (N = 200), CORES (workers for the N bulks of a scenario: a third of half the cores, 14 of 88, since the three figure scripts run at the same time; parallel workers load the installed package) and POSTPROCESS_ONLY (TRUE redraws only).

deconvolution_functions <- bivariate_toy_deconvolution_functions(
  itmax = 200L,
  epsilon = 1e-4,
  intervals = TRUE
)
bivariate <- run_simulation_benchmark(
  scenario_config = bivariate_config,
  deconvolution_functions = deconvolution_functions,
  n = 200L,
  cores = 1L,
  return_data = TRUE
)

config on the in-memory benchmark still carries Shannon entropy of p and MixSim overlap for tests.

Visualisations

Figures are drawn from output/fig02/fig02_plot_data.rds, a list keyed by figure file stem that holds the data handed to each plot function, so any figure can be redrawn without the benchmark:

source("scripts/visualisation/load_visualisation.R")
plot_data <- readRDS("output/fig02/fig02_plot_data.rds")
render_figure_from_data("heatmap_rmse", plot_data, "output/fig02")

Density and performance books share the four correlation corners (𝜚₁, 𝜚₂) ∈ {(0, 0), (−0.8, −0.8), (0.8, 0.8), (−0.8, 0.8)} and the twelve meta-scenarios (CLD × variance × composition). Every log-likelihood panel uses the same axes and the same fill, log₁₀(L/max L) on [−10, 0], with one legend per page: white diamond at p^⋆, red point at the MLE (observed surface) or at the Monte Carlo mean of $\hat{\boldsymbol p}$ (expected surface, whose maximiser is p^⋆ exactly by the information inequality).

• The Q-Q books assess the Wald approximation only; the profile intervals do not rely on it. Squaring a whitened value of −30 gives D² = 900, so the χ² book uses log₁₀ axes.
• The Hellinger distance $\lVert\sqrt{\hat{\boldsymbol p}_A}-\sqrt{\hat{\boldsymbol p}_B}\rVert_2/\sqrt2$ replaces 1 − r for the Pearson correlation: it measures agreement rather than co-movement, it is Euclidean in $\sqrt{\boldsymbol p}$ as Ward’s criterion requires, and it is defined at exact zeros, where the Aitchison distance is not.
• The KS normality lollipops and the MAE heatmaps were dropped: the Q-Q books and the RMSE heatmaps carry the same information.

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 on the closure chart p = w/1^⊤w with the box w ≥ 0 (Byrd et al. 1995; Zhu et al. 1997) (closure chart). The objective is the convolution log-likelihood itself, the box KKT points are exactly the simplex KKT points, and exact zeros are reachable. The archived version instead fitted a box-constrained surrogate and renormalised it, so its faces were not stationary points of ℓ (KKT residual about 0.1).
• 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.3
• Moment generator and ADEMP reporting: How to build synthetic scenarios
• Regular-case MLE checks: Appendix S1

References