Skip to content

PARITY.md

Latest commit

 

History

History
126 lines (92 loc) · 6.05 KB

PARITY.md

File metadata and controls

126 lines (92 loc) · 6.05 KB
Raw

Cross-language parity (R ↔ Python ↔ Web)

mulea ships as three implementations kept in numerical parity: the R package (reference), the mulea Python library, and the mulea web tool (TypeScript + a WebAssembly eFDR core). This document is the single source of truth for how that parity is established and verified.

Gold-standard fixture

One frozen, checksummed reference result generated by the R package is the gold standard that the Python and Web legs validate against.
File python/tests/fixtures/ora_efdr_reference.csv
SHA-256 c118f8918d7650f42325589c05aba5f90c724fce64d38dfcb75c05e7029c7623
Rows 154 terms
Generator python/tests/fixtures/generate_efdr_reference.R (needs R + mulea)
Columns ontology_id, ontology_name, nr_common_with_tested_elements, nr_common_with_background_elements, p_value, eFDR

Inputs / parameters (identical for every leg):

  • Ontology: inst/extdata/Transcription_factor_RegulonDB_Escherichia_coli_GeneSymbol.gmt
  • Target / background: inst/extdata/target_set.txt · inst/extdata/background_set.txt
  • filter_ontology(min_nr_of_elements = 3, max_nr_of_elements = 400)
  • eFDR: number_of_permutations = 100000, random_seed = 42

The fixture is treated as immutable: web/tests/goldStandard.test.ts asserts its SHA-256, so any change must be deliberate (re-generate + update the hash here).

What is compared, and the tolerances

The result splits into a deterministic part (byte-reproducible across languages) and a Monte-Carlo part (the eFDR, which carries resampling noise):

Quantity Comparison Tolerance Verified by
nr_common_* (observed overlaps) exact integer equality 0 Python + Web tests
p_value (hypergeometric) floating-point equality rtol 1e-9 (Py) / 12 dp (Web) test_parity_efdr_with_r.py, efdrMc.test.ts
eFDR numerical agreement ≤ 0.01 (Web), < 0.05 (Py) same
significant set (score < 0.05) set agreement off the 0.05 boundary ±0.005 band test_parity_efdr_with_r.py

Observed agreement (not just the bound): the Web Monte-Carlo eFDR matches the R fixture to max |Δ| ≈ 1.2 × 10⁻³ across the 154 terms (and ≈ 7 × 10⁻⁴ vs the deterministic analytic path). The deterministic columns (p_value, observed counts) match the fixture exactly — that gate must pass before the eFDR tolerance is judged.

Why two eFDR paths agree

The Web/Python eFDR is computed two ways that must converge:

  • exact (analytic) — the n→∞ closed-form limit of the resampling estimator (deterministic), and
  • resampling (Monte-Carlo) — the shared WASM core running number_of_permutations draws.

runAnalysisMc computes both on every run and reports max |eFDR_MC − eFDR_exact| as a convergence diagnostic. The analytic path is the deterministic oracle; the MC path is what reproduces R's resampling method within Monte-Carlo noise.

Derivation: the analytic path is the S→∞ limit of the resampling estimator

The paper (Turek et al. 2024, p. 4) estimates the expected rank by resampling. Write the observed rank of ontology term j as

R_j = Σ_i 1(p_i ≤ p_j)

where p_i is the hypergeometric p-value of term i. In resampling step s, a target set of the same size n = |select| is drawn at random from the background of size N. The overlap of term i (which has m_i = |term_i ∩ background| genes) with that random target is therefore hypergeometrically distributed:

K_i^s ~ Hypergeometric(N, m_i, n),   p_i^s = P_hyp(X ≥ K_i^s),   R_j^s = Σ_i 1(p_i^s ≤ p_j)

mulea's expected rank is the sample mean R̄_j = (1/S) Σ_s R_j^s. By the law of large numbers it converges to its expectation, which has a closed form because each indicator's expectation is a sum over the (finite) support of the hypergeometric law:

R̄_j  ──S→∞──▶  E[R_j^s] = Σ_i P(p_i^s ≤ p_j)
                         = Σ_i  Σ_{k : p_hyp(k; m_i) ≤ p_j}  PMF(k; N, m_i, n)

This double sum is exactly what web/src/efdr.ts:67-94 accumulates as rExp (every (term, k) pair contributes its hypergeometric mass to a cumulative null-mass curve, read off at p_j). Hence the analytic eFDR min(R̄_j / R_j, 1) is the resampling estimator with the Monte-Carlo noise removed — not an approximation of it, but its exact limit.

Convergence rate. Each R_j^s is a sum of Bernoulli indicators, so the resampling mean has Var(R̄_j) = O(1/S); by the CLT the Monte-Carlo estimate deviates from the analytic value by O(1/√S). The benchmark web/bench/efdr-convergence.mjs confirms this empirically: max|Δ| and RMS|Δ| between the WASM Monte-Carlo and the analytic eFDR fall with a fitted slope ≈ −0.5 on a log-log scale (see efdr-convergence.md / .svg).

Note on the clamp: the Web/Python eFDR is clamped to ≤ 1 (min(·, 1)); base R does not clamp. For this fixture all eFDR ≤ 1, so the legs agree; the clamp only differs on terms where the raw expected/observed-rank ratio exceeds 1.

Reproduce

# Web (node) — analytic + MC vs the fixture
cd web && npm test                      # tests/efdrMc.test.ts, tests/goldStandard.test.ts

# Web (real Chromium) — engine + worker in a browser
cd web && npm run test:browser

# Python — exact vs the fixture
cd python && pytest tests/test_parity_efdr_with_r.py

# Regenerate the gold standard (needs R + mulea), then update the SHA-256 above
Rscript python/tests/fixtures/generate_efdr_reference.R

GSEA coverage (ranked-list)

Ranked-list GSEA is currently in R (fgsea) and the Web tool (weighted-KS ES + permutation, validated against fgsea in VALIDATION.md). The Python leg has ORA/eFDR but no GSEA yet — a known parity gap, deferred (a future mulea.gsea would restore full trifecta symmetry). ORA + eFDR remain in all three legs.

Engine performance

Scaling of the WASM Monte-Carlo core (wall-clock vs permutation depth and problem size) is measured by web/bench/efdr-scaling.mjs; the captured table is in web/bench/efdr-scaling.md.