Implement two-group ANCOM-BC and strengthen DA verification

.github/workflows/ci.yml

@@ -19,7 +19,7 @@ jobs:
     strategy:
       fail-fast: false
       matrix:
-        version: ['1.10', '1.12']
+        version: ['1.10', '1.12.5']
     steps:
       - uses: actions/checkout@11d5960a326750d5838078e36cf38b85af677262 # v4
       - uses: julia-actions/setup-julia@4c0cb0fce8556fdb04a90347310e5db8b1f98fb9 # v2
@@ -51,11 +51,23 @@ jobs:
       - name: Proof gate (axiom audit, type-check, reject controls, gate self-test)
         run: proofs/check.sh
   oracle:
-    name: R reference oracle (R 4.5.0, Bioconductor 3.22)
+    name: ${{ matrix.reference }} oracle (R 4.5.0, Bioconductor 3.22)
     runs-on: ubuntu-24.04
     timeout-minutes: 60
+    strategy:
+      fail-fast: false
+      matrix:
+        include:
+          - reference: ALDEx2
+            version: '1.42.0'
+            script: runoracle.jl
+          - reference: ANCOMBC
+            version: '2.12.0'
+            script: ancombc_oracle.jl
     env:
       R_LIBS_USER: ${{ github.workspace }}/.r-lib
+      REFERENCE_PACKAGE: ${{ matrix.reference }}
+      REFERENCE_VERSION: ${{ matrix.version }}
     steps:
       - uses: actions/checkout@11d5960a326750d5838078e36cf38b85af677262 # v4
       - uses: r-lib/actions/setup-r@f9a764fea8d5c63df6ef9a5c7795bf7deb5d7e05 # v2
@@ -65,24 +77,28 @@ jobs:
       - uses: actions/cache@0057852bfaa89a56745cba8c7296529d2fc39830 # v4
         with:
           path: .r-lib
-          key: r-4.5.0-bioc-3.22-oracle-v1
+          key: r-4.5.0-bioc-3.22-${{ matrix.reference }}-oracle-v2
+      - name: Install reference build dependencies
+        run: sudo apt-get update && sudo apt-get install --no-install-recommends -y libcurl4-openssl-dev libssl-dev libxml2-dev libuv1-dev libgsl-dev
       - name: Install the Bioconductor 3.22 references
         run: |
           mkdir -p "$R_LIBS_USER"
-          Rscript -e 'options(Ncpus = 4, timeout = 3600)
+          Rscript -e 'options(Ncpus = 2, timeout = 3600)
             if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
             BiocManager::install(version = "3.22", ask = FALSE, update = FALSE)
-            for (p in c("ALDEx2")) if (!requireNamespace(p, quietly = TRUE)) BiocManager::install(p, ask = FALSE, update = FALSE)
-            for (p in c("ALDEx2")) cat(p, as.character(packageVersion(p)), "\n")'
+            p <- Sys.getenv("REFERENCE_PACKAGE")
+            if (!requireNamespace(p, quietly = TRUE)) BiocManager::install(p, ask = FALSE, update = FALSE)
+            stopifnot(packageVersion(p) == Sys.getenv("REFERENCE_VERSION"))
+            cat(p, as.character(packageVersion(p)), "\n")'
       - uses: julia-actions/setup-julia@4c0cb0fce8556fdb04a90347310e5db8b1f98fb9 # v2
         with:
-          version: '1.12'
+          version: '1.12.5'
       - uses: julia-actions/cache@d10a6fd8f31b12404a54613ebad242900567f2b9 # v2
       # RCall.jl's own CI does the same: libR's directory must be on the loader path.
       - name: Put libR on the loader path
         run: echo "LD_LIBRARY_PATH=$(R RHOME)/lib${LD_LIBRARY_PATH:+:$LD_LIBRARY_PATH}" >> "$GITHUB_ENV"
-      - name: R can load ALDEx2
-        run: Rscript -e 'library(ALDEx2); cat(R.version.string, "ALDEx2", as.character(packageVersion("ALDEx2")), "\n")'
+      - name: R can load the pinned reference
+        run: Rscript -e 'p <- Sys.getenv("REFERENCE_PACKAGE"); library(p, character.only=TRUE); cat(R.version.string, p, as.character(packageVersion(p)), "\n")'
       - name: Julia can start R through RCall
         run: |
           export R_HOME="$(R RHOME)"
@@ -95,4 +111,4 @@ jobs:
       - name: Native vs reference agreement
         run: |
           export R_HOME="$(R RHOME)"
-          julia --project=test/oracle test/oracle/runoracle.jl
+          julia --project=test/oracle test/oracle/${{ matrix.script }}

Project.toml

@@ -1,8 +1,8 @@
 # SPDX-License-Identifier: MPL-2.0
 name = "CompositionalDA"
 uuid = "bb1787b7-c25e-80cd-8559-27e3f8d35816"
-version = "0.2.0"
-description = "Native Julia compositional differential abundance: implemented ALDEx2 (v0.2); planned ANCOM-BC (v0.3). No Python, no TensorFlow; optional R bridge (temporary)."
+version = "0.3.0"
+description = "Native Julia two-group compositional differential abundance with explicit numerical and proof boundaries."
 authors = ["Jonathan D.A. Jewell <j.d.a.jewell@open.ac.uk>"]
 
 [deps]

README.md

@@ -4,21 +4,18 @@
 # CompositionalDA.jl
 
 Native-Julia compositional differential abundance for microbiome count tables.
-Methods (implemented / planned): ALDEx2 is implemented natively (`src/Methods/ALDEx2.jl`);
-ANCOM-BC is planned for v0.3 (not yet in `src/`). Multinomial and Dirichlet-multinomial regression
+Methods: ALDEx2 (`src/Methods/ALDEx2.jl`) and the original two-group ANCOM-BC
+model (`src/Methods/ANCOMBC.jl`). Multinomial and Dirichlet-multinomial regression
 live in a separate package,
 [`CompositionalCounts.jl`](https://github.com/hyperpolymath/CompositionalCounts.jl).
 It needs no Python and no TensorFlow. R is needed only for the optional reference
 backend.
 
-> **Status: v0.2.0.** Implemented: the core layer (closure, centred log-ratio,
-> zero policies, Benjamini–Hochberg, result/provenance types) and two-group
-> ALDEx2 (`aldex2`). ALDEx2 is defined by `docs/ALDEx2-SPEC.md` and agrees
-> with Bioconductor ALDEx2 1.42.0 within Monte-Carlo noise on the mock, gut
-> and soil datasets. Planned: ANCOM-BC lands in v0.3; the optional R reference
-> backend (`ext/CompositionalDARCallExt.jl`) will call Bioconductor ALDEx2/
-> ANCOMBC through RCall, shipped as a Julia package extension.
-> Nothing here claims a method that is not yet in `src/`.
+> **Status: v0.3 development.** Native implementations and reference adapters
+> exist. Release and integration are gated by `docs/ISSUE-23-ACCEPTANCE.md`;
+> implementation is not a declaration of reference conformance. The ALDEx2
+> gate now checks all statistics and shared Monte-Carlo instances. ANCOM-BC's
+> precise model and limitations are described in `docs/ANCOMBC-SPEC.md`.
 
 ## Input contract
 
@@ -43,15 +40,16 @@ a planted mutant (`test/mutants.jl`). The Benjamini–Hochberg step-up (P5: q 
 q ≤ 1, order-preserving, at the original indices) is proved over exact
 rationals, assuming the input was sorted correctly; the sort itself is not
 verified, and nothing about floating point is proved. Statistical properties such as calibration,
-consistency and FDR control are **not** proved. They are to be tested against
-the R reference implementations; that oracle is not built yet. The full table is in
+consistency and FDR control are **not** proved. R reference comparisons are
+numerical checks, not proofs of those properties. PosteriorP and BiasModel
+add finite aggregation and bias identities. The full table is in
 [`proofs/PROOF-STATUS.md`](proofs/PROOF-STATUS.md).
 
 ## The R bridge is temporary by design
 
-**Planned, not yet in the tree.** The optional reference backend will call
+The optional reference backend calls
 Bioconductor ALDEx2/ANCOMBC through RCall, shipped as a Julia package extension.
-It will exist for two purposes: to serve as
+It exists for two purposes: to serve as
 the test oracle, and to help researchers move off R. It has a written sunset
 criterion; see [`docs/R-SUNSET.md`](docs/R-SUNSET.md).
 
@@ -78,9 +76,8 @@ natively is not speed: it is that `proofs/PROOF-STATUS.md` lets a review answer
 *"is this safe and sound, or can it produce funny results?"* claim by claim --
 each CLR property is either proved under `--safe` or marked tested, so a
 pipeline maintainer adopts exactly as much certainty as the method carries,
-instead of borrowing all of it from a paper. The Bioconductor oracle in CI
-today is ALDEx2 -- the difference between "validated against the
-reference" and "hoped to match" for the method it covers. ANCOM-BC's
-reference oracle remains planned; until it lands, ANCOM-BC carries the
-package's own tests, not a cross-check against the Bioconductor
-implementation.
+instead of borrowing all of it from a paper. CI includes separate
+Bioconductor comparisons for ALDEx2 and the original two-group ANCOM-BC.
+The tests' existence is not evidence that they pass: the release receipt
+must identify the source revision, reference versions and actual outcomes.
+See `docs/ISSUE-23-ACCEPTANCE.md` for the independent release gates.

docs/ALDEx2-SPEC.md

@@ -38,6 +38,16 @@ before the code).
 
 ## Notation
 
+ALDEx2's CLR instances and descriptive columns (`diff_btw`, `diff_win`,
+`rab_all`) use **base-2** logarithms. The generic core `clr` uses natural
+logarithms. The dimensionless effect and p-values cannot distinguish these
+bases; the reference gate must compare the descriptive columns too.
+All-zero features are excluded before constructing the Dirichlet posterior
+and named in provenance. Labels must be complete, nonempty and unique for
+features; group lengths must match sample rows. Unrepresentable Gamma draws
+raise an explicit error instead of being floored to an arbitrary positive
+number.
+
 - Input: integer counts with `n` samples and `d` features, and a two-level
   condition per sample.
 - `L1 < L2` are the two condition levels in sorted order.

docs/ANCOMBC-SPEC.md

@@ -0,0 +1,76 @@
+<!-- SPDX-License-Identifier: CC-BY-SA-4.0 -->
+<!-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk> -->
+
+# ANCOM-BC: the implemented model and its limits
+
+`ancombc` implements the **original two-group model** of Lin and Peddada,
+*Analysis of compositions of microbiomes with bias correction*, Nature
+Communications 11, 3514 (2020), doi:10.1038/s41467-020-17041-7. It was derived
+from the published equations, not copied or translated from the GPL R source.
+It is not ANCOM-BC2. Covariates, interactions, repeated measures, global tests
+and pairwise testing of more than two groups are outside this API.
+
+## Data and estimates
+
+Rows are independent samples; columns are taxa. Counts must be nonnegative
+integers, sample totals positive, identifiers unique and groups complete.
+Two samples per group is the computational minimum, not an assurance of
+asymptotic calibration. All-zero taxa and taxa below `prevalence` are removed
+and named in provenance. At least four tested taxa must remain; this is a
+computational refusal rule, not a sufficient statistical sample size.
+
+Structural absence is determined separately in each group from zero
+prevalence; `negative_lower_bound=true` additionally uses the published
+normal lower bound. Those decisions have explicit statuses and **no numerical
+p-value**, effect or q-value. They do not enter BH or the bias fit.
+
+For retained taxa, `y = log(count + pseudocount)`, with explicit default
+pseudocount 1. Within each group, the sample offset is its across-taxon mean
+minus the group grand mean. Taxon means are group means of `y`. The contrast
+is sorted group 2 minus group 1. Residual sum of squares divided by the square
+of the group size gives the variance of each group mean (paper equations 6
+and 14); the two variances add.
+
+The contrast distribution is fitted as three normal components. The middle
+component has known taxon-specific variance `v[i]`; the two alternatives
+have variances `v[i] + extra[k]`. Ordered component means are estimated by
+weighted isotonic projection. Extra variances are nonnegative. An EM/ECM
+iteration evaluates densities in the log domain, verifies nondecreasing
+likelihood and requires both likelihood and parameter convergence.
+Initial means are the 0.1/0.5/0.9 quantiles, weights 0.1/0.8/0.1, and
+alternative extra variances the observed contrast variance.
+
+The variance-coordinate search examines a fixed logarithmic grid, refines
+local brackets by 80 golden-section steps, and retains an update only if it
+improves the conditional objective. This finite numerical search does **not**
+certify a global mixture-likelihood maximum. Component collapse, non-finite
+variance, decreasing likelihood or iteration exhaustion raises
+`CompositionError`; no other method is substituted.
+
+The ordinary effect subtracts the fitted middle mean. With `conserve=true`,
+the effect uses the hard-classification weighted bias estimate (equation 22)
+and the standard error is `sqrt(v[i]) + sqrt(bias_variance)` (equation 35).
+Two-sided normal Wald p-values are BH-adjusted over tested taxa. Natural-log
+effects and all fit settings are recorded in provenance.
+
+## What certainty is available
+
+`BiasModel.agda` proves observational shift invariance, exact bias
+cancellation and the expression for residual estimation error over any
+commutative ring. It does not prove that the estimated bias is the true bias.
+The conditional model can shift abundance and sampling bias in opposite
+directions without changing observations. Identification therefore needs
+assumptions about non-differential taxa; counts cannot establish them alone.
+
+The finite data, independence, weak taxon-dependence, asymptotic normality and
+identification assumptions are not formalised. Neither floating-point
+execution, statistical calibration, FDR control nor biological truth is
+proved. Structural-zero classification is also a statistical decision.
+
+`test/method_ancombc.jl` checks published moment equations against a hand
+calculation, nuisance-shift invariance, planted effects, feature permutation,
+structural zeros and refusals. `test/oracle/ancombc_oracle.jl` is a separate
+conformance gate against ANCOMBC 2.12.0's original `ancombc`. Its thresholds
+were fixed before observing results. **A paper-derived implementation is not
+automatically equivalent to that R release; the gate must pass before such a
+claim or downstream integration is released.**

docs/ISSUE-23-ACCEPTANCE.md

@@ -0,0 +1,51 @@
+<!-- SPDX-License-Identifier: CC-BY-SA-4.0 -->
+<!-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk> -->
+
+# MetaManifold issue 23: acceptance boundary
+
+The requested contribution comprises native ALDEx2 and ANCOM-BC, independently
+useful as a Julia package, followed by a narrow integration into the existing
+Differential Abundance route and tab. Multinomial/Dirichlet-multinomial models
+belong to issue 22 and CompositionalCounts.jl.
+
+Release requires all of the following, independently:
+
+1. Julia tests on the declared supported versions, including the WebUI's
+   pinned Julia 1.12.5, and package-quality checks.
+2. Pinned Agda/stdlib proof checking, axiom audit and negative controls.
+3. ALDEx2 comparison of **every** returned statistic, plus identical-instance
+   Welch/Wilcoxon p-value and expected-BH comparisons. Comparing only
+   dimensionless effects and p-values cannot detect the wrong logarithm base.
+4. ANCOM-BC conformance to the declared reference version and settings.
+   A failure must remain visible; numerical thresholds cannot be widened
+   after observing results to turn disagreement into success.
+5. A registered Julia release, installable by name/UUID with registry pins.
+   Joshua's issue comment explicitly rejects URL-pinned dependencies.
+6. WebUI integration on its current main, a registry-resolved Manifest,
+   no changes to renv.lock, and all required backend/frontend/server checks.
+7. Signed, GitHub-verified commits, no merge commits, correct per-file SPDX,
+   passing REUSE and a PR describing versions, validation and limitations.
+
+A local development environment may use `Pkg.develop` to test an integration.
+It is not a compliant released dependency or a substitute for registration.
+Do not submit an upstream WebUI PR with a development path or invented
+registry tree hash in its Manifest. Registration/maintainer acceptance are
+external conditions, not facts that local tests can establish.
+
+## Mathematical dependencies inspected
+
+Mathesis's discovery atlas already points to this package's CLR/BH proofs.
+Its README explicitly makes it an index rather than a proof authority; new
+method-specific proofs belong here. CompositionalCounts.jl has a different
+model and does not supply the missing ANCOM-BC estimator or calibration
+proof. JanusKey concerns reversible file operations; Absolute Zero concerns
+null effects/disclosures; Trope Checker verifies composition of declared IR
+grades. None of those inspected interfaces establishes a DA sampling model
+or FDR guarantee. Residual evidence types are useful for distinguishing
+claims, assumptions and missing evidence, but adding them as a runtime
+dependency would not prove the numerical implementation.
+
+The useful mathematical additions are therefore local: PosteriorP's finite
+minimum/sum inequality and BiasModel's cancellation/observability identities.
+The exact rational and ring statements, the floating-point implementation,
+the reference comparison and the statistical assumptions remain separate.

docs/R-SUNSET.md

@@ -3,17 +3,19 @@
 
 # Sunset criterion for the R reference backend
 
-> **Status (v0.2.0):** ALDEx2 is release 1 of 2 within tolerance (oracle run of
-> 2026-10-02, all 13 checks passing). The frozen fixtures under `test/fixtures/`
-> described below do not exist yet.
+> **Status (v0.3 development):** the previous 13-check ALDEx2 oracle did not
+> cover descriptive-column units or identical-instance q-values. The stronger
+> gate must establish a new baseline before a release qualifies for sunset.
+> ANCOM-BC now has a reference adapter and conformance gate; their existence
+> does not mean the gate has passed. Frozen reference-output fixtures under
+> `test/fixtures/` do not yet exist.
 
 The R backend (`ext/CompositionalDARCallExt.jl`, loaded only when RCall is
 present) is a bridge, not a product. It never owns R state. Every R call runs inside a
 caller-supplied `runner` function, so a host application keeps its own R session
 and lock (MetaManifold passes its `with_r_lock`).
 
-**Criterion.** For each implemented method (currently ALDEx2; ANCOM-BC is
-planned for v0.3), once the native implementation agrees with the R reference
+**Criterion.** For each implemented method, once the native implementation agrees with the R reference
 within the tolerances stated in `test/oracle/runoracle.jl` on all three
 reference datasets (mock, gut, soil) for **two consecutive releases**, that
 method's R path is marked deprecated in the next release and receives no