---
title: "§2.2 Bivariate toy model (G = 2 genes, J = 2 cell types)"
subtitle: "Full script: `scripts/fig02_bivariate_toy.R`"
author: "Bastien Chassagnol"
date: "`r Sys.Date()`"
vignette: >
  %\VignetteIndexEntry{§2.2 Bivariate toy model}
  %\VignetteEngine{quarto::html}
  %\VignetteEncoding{UTF-8}
knitr:
  opts_chunk:
    collapse: true
    comment: "#>"
bibliography:
  - ../inst/packages.bib
  - ../inst/REFERENCES.bib
link-citations: true
crossref:
  eq-prefix: "Eq."
  tbl-prefix: "Table"
  fig-prefix: "Fig."
  sec-prefix: "Sec."
format:
  html:
    embed-resources: true
    toc: true
    toc-depth: 3
execute:
  echo: true
  eval: false
  warning: false
---

```{r}
#| label: setup
#| include: false
#| eval: true
options(tinytable_html_mathjax = TRUE)
mu_small <- matrix(
  c(20, 22, 22, 20),
  nrow = 2,
  dimnames = list(c("gene_1", "gene_2"), c("ct1", "ct2"))
)
p_bal <- c(0.5, 0.5)
```

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()`](../reference/run_simulation_benchmark.html)
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](04-1-synthetic-scenarios.html#sec-ademp).

```{.tikz}
%%| label: fig-bivariate-design
%%| fig-cap: Bivariate toy factorial design: three compositions, two variance structures, two CLD mean separations, a 9-by-9 gene-gene correlation grid, and seven solvers (3 x 2 x 2 x 9 x 9 = 972 scenarios).
%%| filename: fig-bivariate-design
%%| alt: Four design factors, three compositions, two variance structures, two CLD mean separations, and a 9 by 9 grid of gene-gene correlations rho from -0.8 to 0.8, combine into 972 scenarios, each analysed by 7 deconvolution solvers.
\begin{tikzpicture}
\node[box, text width=3.1cm] (P) at (0,0) {3 compositions\\ Ba / Mo / Hi};
\node[box, text width=3.1cm] (V) at (3.6,0) {2 variance structures\\ Ho homoscedastic,\\ He heteroscedastic};
\node[box, text width=3.1cm] (C) at (7.2,0) {2 CLD mean separations\\ Sm small, Lg large};
\node[box, text width=3.1cm] (R) at (10.8,0) {$9\times 9$ gene--gene correlations\\ $\varrho\in[-0.8,0.8]$, step $0.2$};
\node[keybox] (S) at (5.4,-2.2) {972 scenarios\\ $3\times2\times2\times9\times9$};
\node[keybox] (A) at (5.4,-3.9) {7 deconvolution solvers};
\foreach \n in {P,V,C,R} \draw[edge] (\n.south) -- (S.north);
\draw[edge] (S) -- (A);
\end{tikzpicture}
```

```{mermaid}
%%| label: fig-bivariate-design-tree
%%| eval: true
%%| fig-cap: "Hierarchical reading of the same factorial: small versus large CLD, then homoscedastic versus heteroscedastic variances, then balanced / moderately unbalanced / highly unbalanced compositions (2 x 2 x 3 = 12 meta-scenarios). Each leaf is crossed with the 9-by-9 gene-gene correlation grid."
---
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 {#sec-bivariate-toy}

With $p_1+p_2=1$, only one free ILR coordinate is estimated (see
[derivatives under simplex transforms](https://bastienchassagnol.github.io/DeCovarT/articles/01-1-population-generative-optimisation.html)).
Even when mean profiles are similar, gene–gene correlation can degrade
mean-only solvers, and is partly recovered once
$\boldsymbol{\Sigma}(\boldsymbol{p})=\sum_j p_j^2\boldsymbol{\Sigma}_j$
enters the likelihood.

::: {#fig-toy-design layout-ncol=2}

![Hyperparameter grid used in the bivariate manuscript study: two
centroid geometries, balanced versus unbalanced $\boldsymbol{p}$, a
gene–gene correlation sweep, and the first-generation plus DeCovarT
solvers.](figures/fig_toy_model_parameters.svg){#fig-toy-hyperparams}

![AI-generated 2D density of the two-gene mixture while only the
pairwise correlation changes, at fixed means, marginal variances, and
composition. The ggplot panel below recomputes the same contrast and
annotates MixSim overlap.](figures/fig_toy_2d_density.svg){#fig-toy-density-ai}

Factorial design (left) and correlation-only 2D densities (right) for
the bivariate toy. {#fig-toy-design}
:::

## Factorial design (972 scenarios)

| Factor | Levels |
|---|---|
| Composition $\boldsymbol{p}$ | balanced $(1/2,\,1/2)$; moderately unbalanced $(17/20,\,3/20)$; highly unbalanced $(99/100,\,1/100)$ |
| Mean separation (CLD) | small: $\boldsymbol{\mu}_{\cdot 1}=(20,22)$, $\boldsymbol{\mu}_{\cdot 2}=(22,20)$; large: $(20,40)$, $(40,20)$ |
| Gene–gene corr. CT 1 ($\varrho_1$) | $\{-0.8,\,-0.6,\,\ldots,\,0.8\}$ — 9 levels |
| Gene–gene corr. CT 2 ($\varrho_2$) | $\{-0.8,\,-0.6,\,\ldots,\,0.8\}$ — 9 levels |
| Variance structure | homoscedastic $\sigma^2 = (1,1)$; heteroscedastic $\sigma^2 = (1,2)$ |
| **Total scenarios** | $3 \times 2 \times 9 \times 9 \times 2 = \mathbf{972}$ |

$N = 500$ Monte Carlo replicates per scenario.

**CLD** is centroid (mean-profile) Euclidean separation
$\lVert\boldsymbol{\mu}_{\cdot 1}-\boldsymbol{\mu}_{\cdot 2}\rVert_2$,
not an Aitchison distance; see
[the descriptor notes](04-1-synthetic-scenarios.html#nte-desc-cld).
Small CLD uses $\boldsymbol{\mu}_{\cdot 1}=(20,22)$,
$\boldsymbol{\mu}_{\cdot 2}=(22,20)$ ($d=2\sqrt{2}\approx 2.83$,
cosine $\approx 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.

| Token | Meaning |
|---|---|
| `B` | Bivariate toy ($G=2$, $J=2$) |
| `1`–`972` | Factorial row index |
| `Ho` / `He` | Homoscedastic / heteroscedastic marginal variances |
| `Ba` / `Mo` / `Hi` | Balanced / moderately unbalanced / highly unbalanced $\boldsymbol{p}$ |
| `Sm` / `Lg` | Small / large CLD (Euclidean centroid gap) |

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
$\times$ algorithm $\times$ 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.

```{r}
#| label: gen-model
library(DeCovarT)
source("scripts/fig02_bivariate_toy.R")

bivariate_config <- build_bivariate_scenario_config()
nrow(bivariate_config)
```

```{r}
#| label: fig-toy-density-ggplot
#| fig-cap: >-
#|   Two-gene mixture densities at fixed means
#|   $\boldsymbol{\mu}_{\cdot 1}=(20,22)$, $\boldsymbol{\mu}_{\cdot 2}=(22,20)$,
#|   $\boldsymbol{p}=(1/2,1/2)$ and unit diagonals. Only the common
#|   gene–gene correlation $\varrho$ changes. Labels report MixSim
#|   $\bar\varpi$ (`compute_average_overlap()`).
#| fig-width: 8.4
#| fig-height: 4.4
#| echo: false
set.seed(20260828)
rho_grid <- c(-0.6, 0, 0.6)
n_draw <- 400L
dens_df <- list()
ovl_df <- list()
for (rho in rho_grid) {
  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)
  )
  sim <- DeCovarT::simulate_bulk_mixture(
    mu_small,
    Sigma,
    p = p_bal,
    n = n_draw
  )
  dens_df[[length(dens_df) + 1L]] <- data.frame(
    rho = sprintf("rho = %s", format(rho, nsmall = 1)),
    gene_1 = sim$Y[1, ],
    gene_2 = sim$Y[2, ]
  )
  theta <- list(p = p_bal, mu = mu_small, sigma = Sigma)
  ovl_df[[length(ovl_df) + 1L]] <- data.frame(
    rho = sprintf("rho = %s", format(rho, nsmall = 1)),
    ovl = DeCovarT::compute_average_overlap(theta)
  )
}
dens_df <- do.call(rbind, dens_df)
ovl_df <- do.call(rbind, ovl_df)
ovl_df$lab <- sprintf("overlap = %.2f", ovl_df$ovl)
ggplot2::ggplot(dens_df, ggplot2::aes(x = gene_1, y = gene_2)) +
  ggplot2::stat_density_2d_filled(contour_var = "ndensity", bins = 12) +
  ggplot2::geom_label(
    data = ovl_df,
    ggplot2::aes(x = Inf, y = Inf, label = lab),
    inherit.aes = FALSE,
    hjust = 1.05,
    vjust = 1.2,
    size = 3.1,
    fill = "#fffbe9",
    label.size = 0.15
  ) +
  ggplot2::facet_wrap(~rho, nrow = 1) +
  ggplot2::guides(fill = "none") +
  ggplot2::labs(
    x = "Expression of gene 1",
    y = "Expression of gene 2"
  ) +
  ggplot2::theme_bw(base_size = 11)
```

## Log-likelihood along the simplex

For $J=2$ the free coordinate is $p_1\in(0,1)$. The listing below
evaluates
$\ell_{\boldsymbol{y}\mid\boldsymbol{\zeta}}(p_1,1-p_1)$ 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](https://github.com/quarto-dev/quarto-r/issues/217)).
The same likelihood path is checked in `tests/testthat/`.

```{r}
#| label: fig-toy-loglik
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 {#sec-bivariate-inference}

Five methods are evaluated on each scenario:

| Method | Description | Intervals |
|---|---|---|
| `lsei` | DeconRNASeq: least squares on the simplex by `limSolve::lsei` [@gongDeconRNASeqStatisticalFramework2013] | none |
| `NNLS` | Lawson–Hanson non-negative least squares, renormalised [@dessoleLawsonHansonAlgorithmDeviation2023] | none |
| `Newton-Raphson` | Second-order ascent in ILR coordinates | Wald, ILR, profile, oracle |
| `L-BFGS-B` | Limited-memory quasi-Newton on the closure chart $\boldsymbol p=\boldsymbol w/\mathbf 1^{\top}\boldsymbol w$, $\boldsymbol w\ge0$ | Wald, ILR |
| `Firth` | Jeffreys-penalised likelihood [@firthBiasReductionMaximum1993] | Wald, ILR |

: Methods of the bivariate benchmark. {#tbl-bivariate-methods}

- **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 $\nu$-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 ($\boldsymbol\beta=c\boldsymbol p$,
  $c=\mathbf 1^{\top}\boldsymbol\beta$); `lsei` fixes $c=1$. The two agree
  only when $\hat c=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](01-2-population-inference.html): `wald` (delta
method on $\boldsymbol p$), `ilr` (Wald in ILR coordinates mapped back),
`profile` (likelihood-ratio inversion) and `oracle` ($\hat p_j$ plus or
minus the Wald half-width at $\boldsymbol p^{\star}$). 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\times B$ refits per solver dominated the run time.
Performance is summarised with the ADEMP framework
[@morrisUsingSimulationStudies2019]: **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).

```{r}
#| label: inference
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
$\boldsymbol{p}$ and MixSim overlap for tests.

# Visualisations {#sec-bivariate-viz}

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:

```{r}
#| label: redraw
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
$(\varrho_1,\varrho_2)\in\{(0,0),(-0.8,-0.8),(0.8,0.8),(-0.8,0.8)\}$ and the
twelve meta-scenarios (CLD $\times$ variance $\times$ composition).
Every log-likelihood panel uses the same axes and the same fill,
$\log_{10}(L/\max L)$ on $[-10,0]$, with one legend per page: white
diamond at $\boldsymbol p^{\star}$, red point at the MLE (observed
surface) or at the Monte Carlo mean of $\hat{\boldsymbol p}$ (expected
surface, whose maximiser is $\boldsymbol p^{\star}$ exactly by the
information inequality).

| Output | Description |
|---|---|
| `density_visualisations/purified_density.pdf` | The two cell-type Gaussians |
| `density_visualisations/bulk_density.pdf` | The convolution that DeCovarT fits |
| `density_visualisations/loglik_surface_p.pdf`, `loglik_ilr_profile.pdf` | Log-likelihood of one bulk on $p_1$ and on the ILR coordinate |
| `density_visualisations/loglik_expected_surface_p.pdf`, `loglik_expected_ilr_profile.pdf` | Expected log-likelihood and Monte Carlo mean of $\hat{\boldsymbol p}$ |
| `density_visualisations/loglik_rgl.pdf` | 3D relative likelihood (`rgl`) |
| `performance_visualisations/heatmap_rmse.pdf`, `heatmap_aitchison.pdf` | Tiles on $(\varrho_1,\varrho_2)$, one page per meta-scenario; grey = no finite estimate |
| `performance_visualisations/ci_width_forest.pdf` | Mean interval, coverage [Wilson interval] and width per method and construction |
| `performance_visualisations/p_hat_monte_carlo_raincloud.pdf` | Monte Carlo law of $\hat p$, including pile-up at 0 and 1 |
| `performance_visualisations/p_hat_monte_carlo_quantile_dots.pdf` | The same law as quantile dots |
| `performance_visualisations/kkt_residual_raincloud.pdf` | Projected-score KKT residual |
| `performance_visualisations/convergence_theoretical_vs_numerical_failures.pdf` | Theoretical success, numerical failure, theoretical failure |
| `performance_visualisations/wald_normality_qq.pdf` | Whitened ILR estimate against $N(0,1)$ |
| `performance_visualisations/wald_mahalanobis_chi2_qq.pdf` | $D^2$ against $\chi^2_1$ on $\log_{10}$ axes, coloured by convergence |
| `performance_visualisations/estimate_agreement_hellinger_ward.pdf` | Mean Hellinger distance between paired estimates, Ward clustering; dashed line at the mean over method pairs |
| `performance_visualisations/rmse_colour_aitchison_size_dots.pdf` | Colour = mean RMSE, size = mean Aitchison distance |
| `performance_visualisations/runtime.pdf`, `memory.pdf` | $\log_{10}$ seconds and $\log_{10}$ MiB |

: Plot books of the bivariate scenario. {#tbl-fig02-books}

- 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^2=900$, so the $\chi^2$ book uses $\log_{10}$ 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.

::: {#nte-rgl-alternatives .callout-tip}
## Other one-sample log-likelihood slices (not drawn)

The current `loglik_rgl.pdf` holds the two mean profiles fixed and
varies only $\boldsymbol{p}$ (or its ILR coordinate).
White diamonds mark $p^{\star}$ and red circles the numerical MLE. Expected log-likelihood
surfaces
$\mathbb{E}_{\boldsymbol{y}\mid\boldsymbol{p}^{\star}}[\ell(\boldsymbol{p};\boldsymbol{y})]$
are drawn separately as `loglik_expected_surface_p.pdf` and
`loglik_expected_ilr_profile.pdf` (not overlaid on the one-draw
books). Useful alternatives, left as outlook:

- Hold $\boldsymbol{p}^{\star}$ fixed and vary one column of $\boldsymbol{\mu}$
  (a mean-shift slice) to show how centroid collision flattens $\ell$.
- Hold $\boldsymbol{p}^{\star}$ and $\boldsymbol{\mu}$ fixed and vary one
  correlation $\varrho_j$, so the surface is a covariance diagnostic
  rather than a mixing-weight diagnostic.
- Overlay several independent bulk draws as translucent surfaces to
  display sampling variability of $\ell(\boldsymbol{p};\boldsymbol{y})$.
:::

## Likelihood geometry and solver behaviour {#sec-bivariate-findings}

::: {.callout-note}
The numbers in this section come from the archived run of
`output/archives/fig02/` (Wald intervals only), which still included
NNLS, gradient, simulated annealing, Marquardt–Levenberg and the earlier
L-BFGS-B that renormalised a box-constrained fit. They will be replaced
by those of the production run; the mechanisms do not depend on the
solver list.
:::

The factorial heatmaps remain the right *design* summary: large CLD and
$\varrho_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 $\varrho_1=\varrho_2$. That saddle is a statement about
**identifiability of $\boldsymbol{\mu}\boldsymbol{p}$**, not about which
numerical method is used. The pack in `output/fig02/mle_explanation/`
(scenario descriptors, ILR profiles evaluated at
$\boldsymbol{y}=\boldsymbol{\mu}\boldsymbol{p}^{\star}$, expected-Fisher Wald
standard errors, and Newton versus Marquardt forests) is what separates
three mechanisms that the heatmaps mix: the shape of
$\ell(\boldsymbol{p};\boldsymbol{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
$(\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})^{\top}
(\boldsymbol{y}-\boldsymbol{\mu}\boldsymbol{p})$ on the simplex
[@dessoleLawsonHansonAlgorithmDeviation2023;
@gongDeconRNASeqStatisticalFramework2013; @R-limSolve]. For a convex
objective a local minimiser is global
[@boydConvexOptimization2004, Ch. 4], so those fits are unique
solutions of the *least-squares* problem. They are **not** maximisers of
the convolution log-likelihood: $\boldsymbol{\Sigma}(\boldsymbol{p})$ never
enters the quadratic. When $\boldsymbol{\mu}_{\cdot 1}$ and
$\boldsymbol{\mu}_{\cdot 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 $\ell$ 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](https://bastienchassagnol.github.io/DeCovarT/articles/01-2-population-inference.html)
(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 $\boldsymbol p=\boldsymbol w/\mathbf 1^{\top}\boldsymbol w$ with the
  box $\boldsymbol w\ge0$
  [@byrdLimitedMemoryAlgorithm1995; @zhuAlgorithm778LBFGSB1997]
  ([closure chart](https://bastienchassagnol.github.io/DeCovarT/articles/01-1-population-generative-optimisation.html#sec-closure-box)).
  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 $\ell$ (KKT residual about $0.1$).
- `Newton-Raphson` (`stats::nlminb`) and `Marquardt-Levenberg`
  (`marqLevAlg`; [@R-marqLevAlg; @philippsRobustEfficientOptimization2021])
  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 [@marquardtAlgorithmLeastSquaresEstimation1963] even
  when the Hessian of $\ell$ is indefinite, and stops on the
  relative-distance-to-maximum (RDM) criterion: the ratio of
  *optimisation* error to *statistical* error
  [@commengesNewtonLikeAlgorithmLikelihood2006]. 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 $\ell$ 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](https://bastienchassagnol.github.io/DeCovarT/articles/01-2-population-inference.html#sec-finite-sample):
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 $\boldsymbol{p}^{\star}$ is
what would make $\hat{\boldsymbol{p}}_N$ consistent. The ILR profiles in
the pack are evaluated at the *conditional mean*
$\boldsymbol{y}=\boldsymbol{\mu}\boldsymbol{p}^{\star}$, which isolates the
log-determinant remainder; the forests average over random
$\boldsymbol{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 $\ell(\rho)$ on the plotting grid).

**Large CLD, balanced, $\varrho_j=0$ (`B567`).** The mean map is well
conditioned. The ILR profile spans about $101$ nats and
$\mathrm{SE}_{\rho}=0.071$. Newton and Marquardt agree to numerical
noise: bias of $p_1$ 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, $\varrho_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$,
$\mathrm{SE}_{\rho}=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;
$\mathrm{SE}_{\rho}$ is $0.95$ (both $\varrho=-0.8$) or $0.71$
($\varrho=0$). RMSE of $p_1$ 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^{\star}=(0.99,0.01)$).**
The ILR chart sends the rare type to $\rho\to+\infty$. Fisher curvature
in $\rho$ collapses: $\mathrm{SE}_{\rho}$ is $24$ (both $\varrho=-0.8$) or
$20$ ($\varrho=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_1$ 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](https://bastienchassagnol.github.io/DeCovarT/articles/01-2-population-inference.html#sec-boundary)):
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
$\rho$ 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_2$ 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
$|\rho|$ is large enough to represent $p_2=0.01$. A sharp peak in $p$
can still look like a long walk in $\rho$; 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\times 9$ $(\varrho_1,\varrho_2)$ tiles, which are **not** in the
twelve-row explanation pack, remain the place to see mean-only methods
lose first as $|\varrho|$ 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^{\star}$ 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](R2-covariance-driven.html)
- Moment generator and ADEMP reporting:
  [How to build synthetic scenarios](04-1-synthetic-scenarios.html)
- Regular-case MLE checks:
  [Appendix S1](S1-identifiability.html)

## References