Fix ANCOM-BC estimation and qualify pinned reference comparisons

.github/workflows/ci.yml

@@ -3,7 +3,7 @@
 name: CI
 on:
   push:
-    branches: [main]
+    branches: [main, 'fix/issue-23-reference-pins']
     tags: ['v*']
   pull_request:
 permissions:
@@ -77,18 +77,24 @@ jobs:
       - uses: actions/cache@0057852bfaa89a56745cba8c7296529d2fc39830 # v4
         with:
           path: .r-lib
-          key: r-4.5.0-bioc-3.22-${{ matrix.reference }}-oracle-v2
+          key: r-4.5.0-bioc-3.22-cran-2026-02-13-${{ matrix.reference }}-oracle-v3
       - 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: |
+          export R_LIBS_USER="$GITHUB_WORKSPACE/.r-lib"
+          echo "R_LIBS_USER=$R_LIBS_USER" >> "$GITHUB_ENV"
           mkdir -p "$R_LIBS_USER"
           Rscript -e 'options(Ncpus = 2, timeout = 3600)
+            # ANCOMBC 2.12.0 imports the pre-S7 CVXR API. Freeze CRAN as well
+            # as Bioconductor; latest CRAN no longer supplies that API.
+            options(repos = c(CRAN = "https://packagemanager.posit.co/cran/__linux__/noble/2026-02-13"))
             if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
             BiocManager::install(version = "3.22", ask = FALSE, update = FALSE)
             p <- Sys.getenv("REFERENCE_PACKAGE")
             if (!requireNamespace(p, quietly = TRUE)) BiocManager::install(p, ask = FALSE, update = FALSE)
             stopifnot(packageVersion(p) == Sys.getenv("REFERENCE_VERSION"))
+            if (p == "ANCOMBC") stopifnot(packageVersion("CVXR") == "1.0.15")
             cat(p, as.character(packageVersion(p)), "\n")'
       - uses: julia-actions/setup-julia@4c0cb0fce8556fdb04a90347310e5db8b1f98fb9 # v2
         with:
@@ -109,6 +115,18 @@ jobs:
           julia --project=test/oracle -e 'using RCall; println(rcopy(R"R.version.string"))'
           julia --project=test/oracle -e 'using CompositionalDA, RCall; @assert Base.get_extension(CompositionalDA, :CompositionalDARCallExt) !== nothing; println("extension loaded")'
       - name: Native vs reference agreement
+        if: matrix.reference == 'ALDEx2'
         run: |
           export R_HOME="$(R RHOME)"
           julia --project=test/oracle test/oracle/${{ matrix.script }}
+      - name: Pristine ANCOMBC known-defect control (three expected failures)
+        if: matrix.reference == 'ANCOMBC'
+        run: |
+          export R_HOME="$(R RHOME)"
+          ANCOMBC_PRISTINE_DEFECT_CONTROL=1 julia --project=test/oracle test/oracle/ancombc_oracle.jl
+      - name: ANCOMBC conformance with the pinned upstream variance repair
+        if: always() && matrix.reference == 'ANCOMBC'
+        run: |
+          curl --fail --location --retry 3 https://raw.githubusercontent.com/hyperpolymath/ANCOMBC/6d094731ed85e9b0458053ad3510e9bdc238e152/R/ancombc_bias_correct.R -o /tmp/ancombc-variance-repair.R
+          export R_HOME="$(R RHOME)"
+          ANCOMBC_ESTIMATOR_SOURCE=/tmp/ancombc-variance-repair.R julia --project=test/oracle test/oracle/ancombc_oracle.jl

README.md

@@ -17,6 +17,12 @@ backend.
 > gate now checks all statistics and shared Monte-Carlo instances. ANCOM-BC's
 > precise model and limitations are described in `docs/ANCOMBC-SPEC.md`.
 
+The v0.3 implementation, tests, proofs and integration work were developed
+with OpenAI Codex assistance. General registration requires a human maintainer
+to review and understand these contributions; automated checks and code review
+do not replace that review. Verification scope is recorded in
+`proofs/PROOF-STATUS.md` and `docs/ISSUE-23-ACCEPTANCE.md`.
+
 ## Input contract
 
 Raw integer counts, as a samples × features matrix. Normalised or transformed

docs/ANCOMBC-SPEC.md

@@ -6,7 +6,9 @@
 `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.
+from the published equations. The R estimator source was subsequently inspected
+to diagnose reference discrepancies and document its initialization, stopping
+and variance conventions; this is not a clean-room conformance exercise.
 It is not ANCOM-BC2. Covariates, interactions, repeated measures, global tests
 and pairwise testing of more than two groups are outside this API.
 
@@ -19,7 +21,9 @@ 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
+Both native and reference APIs default to `structural_zeros=false`, matching
+the reference method. With `structural_zeros=true`, 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.
@@ -37,6 +41,9 @@ 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.
+Both stopping checks use `tolerance`; the parameter check is scaled by the
+largest absolute parameter. A square-root tolerance would stop too early on
+the soil fixture's slowly converging fit.
 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.
 
@@ -47,9 +54,19 @@ 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).
+The effect subtracts the fitted middle mean, independently of `conserve`,
+matching the pinned reference API's distinction between coefficient and
+variance estimation. With `conserve=true`, the standard error is
+`sqrt(v[i]) + sqrt(bias_variance)`, using the WLS variance approximation in
+equation 31 and the conservative bound in equation 35. The paper's equation
+35 also substitutes the WLS point estimate in its numerator; this API retains
+the EM coefficient when changing the uncertainty option. The WLS point
+estimate remains in provenance for comparison, not as the reported effect.
+For this WLS approximation, taxa below the fitted lower-component proportion's
+quantile form the lower class; taxa at or above the upper-tail quantile form
+the upper class. The middle interval is lower-inclusive and upper-exclusive,
+matching the reference's tie conventions. This is recorded as
+`bias_classification = "mixture-weight-quantiles"` in provenance.
 Two-sided normal Wald p-values are BH-adjusted over tested taxa. Natural-log
 effects and all fit settings are recorded in provenance.
 
@@ -74,3 +91,44 @@ 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.**
+
+The reference wrapper allows 20,000 iterations by default. The pinned R
+implementation returns its current parameters at its iteration cap without a
+convergence error. The oracle records its internal iteration count and residual
+and requires actual convergence. Both original fits converged below 1,000
+iterations: increasing the reference budget did not resolve their disagreement.
+
+## Reference defect and the qualified comparison
+
+Unmodified ANCOMBC 2.12.0 fails the soil effect, p-value and q-value thresholds.
+Its variance optimizer can stop at zero despite an improving positive direction.
+For a single squared residual 1 and known variance 0.99, the exact conditional
+optimum is extra variance 0.01. The bounded Nelder-Mead call used by the reference
+returns zero when started there. This defect also occurs in the latest upstream
+original-ANCOM-BC estimator. It is addressed in
+[ANCOMBC PR 344](https://github.com/FrederickHuangLin/ANCOMBC/pull/344).
+
+For fixed residuals r_i, known positive variances v_i and nonnegative weights
+w_i, the conditional negative log likelihood, apart from a constant, is
+Q(t) = (1/2) sum_i w_i [log(v_i+t) + r_i^2/(v_i+t)]. Its derivative is
+(1/2) sum_i w_i [(v_i+t)-r_i^2]/(v_i+t)^2. It is nonnegative once
+t >= max_i r_i^2, supplying a finite search interval. In the equal-variance
+case the optimum is max(sum_i w_i r_i^2 / sum_i w_i - v, 0).
+These are calculus derivations, not additional Agda theorems.
+
+The positive conformance gate loads only the bias estimator and its scalar
+variance helper from signed upstream-repair commit
+`6d094731ed85e9b0458053ad3510e9bdc238e152` into the otherwise pinned ANCOMBC
+2.12.0 reference environment. The native Julia optimizer and the R repair use
+different searches (grid/golden-section versus Brent). All 27 checks passed at
+native revision `480b10bd405c2766f08408105435b33d5dd05050`, including convergence
+checks, with the original thresholds unchanged. The largest soil discrepancies
+were 3.57e-6 in effect, 4.21e-5 in SE, 1.10e-4 in p and 1.54e-4 in q.
+
+A separately named pristine-reference control reproduces the analytic boundary
+failure and marks exactly those three soil comparisons as expected failures.
+An unexpected pass fails that control. Without that explicit control mode,
+the unmodified-reference oracle still fails normally. A green combined CI is
+**not a claim of equivalence to unmodified ANCOMBC 2.12.0**: qualification is
+to the stated estimator repair, pending upstream review. Changing that pin
+requires a fresh comparison and an updated receipt.

docs/ISSUE-23-ACCEPTANCE.md

@@ -19,12 +19,20 @@ Release requires all of the following, independently:
 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.
+   The current qualification uses the explicitly pinned variance-estimator
+   repair documented in `ANCOMBC-SPEC.md`. Unmodified ANCOMBC 2.12.0 has three
+   known soil discrepancies, reproduced as a separate negative control; it
+   is not the passing conformance target. Upstream review remains pending.
 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.
+8. Human maintainer review and understanding of the AI-assisted code before
+   General registration. General requires Registrator for GitHub packages and
+   normally holds new-package registrations for three days. See its
+   [registration and LLM policies](https://github.com/JuliaRegistries/General#registering-a-package-in-general).
 
 A local development environment may use `Pkg.develop` to test an integration.
 It is not a compliant released dependency or a substitute for registration.

ext/CompositionalDARCallExt.jl

@@ -34,7 +34,7 @@ function CompositionalDA.aldex2_reference(counts::AbstractMatrix{<:Integer}, fea
     df = runner() do
         reads = permutedims(Float64.(counts[:, keep]))     # R integers are only 32-bit
         @rput reads
-        conds = design.groups; lv = levels
+        conds = [g == levels[1] ? "0" : "1" for g in design.groups]; lv = ["0", "1"]
         @rput conds lv
         fnames = design.features[keep]; @rput fnames
         R"""
@@ -43,7 +43,8 @@ function CompositionalDA.aldex2_reference(counts::AbstractMatrix{<:Integer}, fea
         colnames(reads) <- paste0("s", seq_len(ncol(reads)))
         # A character vector, not a factor: ALDEx2 1.42 calls round() on conds,
         # which errors for a factor. It builds its own factor, so the level
-        # order is the same sorted order as `lv`.
+        # order is preserved by the explicit 0/1 encoding, even for original
+        # numeric labels such as 2 and 10 whose lexical order differs.
         stopifnot(identical(sort(unique(conds)), lv))
         res <- ALDEx2::aldex(reads, conds, mc.samples = $mc_samples,
                              test = "t", effect = TRUE, denom = "all", verbose = FALSE)
@@ -65,27 +66,31 @@ function CompositionalDA.aldex2_reference(counts::AbstractMatrix{<:Integer}, fea
 end
 
 function CompositionalDA.ancombc_reference(counts::AbstractMatrix{<:Integer}, features::AbstractVector,
-    groups::AbstractVector; prevalence::Real=0.1, structural_zeros::Bool=false,
+    groups::AbstractVector; prevalence::Real=0.1, pseudocount::Real=1.0, structural_zeros::Bool=false,
     negative_lower_bound::Bool=false, conserve::Bool=true, tolerance::Real=1e-7,
-    max_iterations::Integer=1000, runner=f -> f())
+    max_iterations::Integer=20000, runner=f -> f())
     design = validate_two_groups(counts, features, groups)
     0 <= prevalence < 1 || throw(ArgumentError("prevalence must be in [0,1)"))
+    isfinite(pseudocount) && pseudocount > 0 || throw(ArgumentError("pseudocount must be finite and positive"))
     isfinite(tolerance) && 0 < tolerance < 1 || throw(ArgumentError("tolerance must be in (0,1)"))
     max_iterations > 0 || throw(ArgumentError("max_iterations must be positive"))
     table, version = runner() do
         reads = permutedims(Float64.(counts))
         fnames = design.features
         conds = [g == design.levels[1] ? "reference" : "contrast" for g in design.groups]
-        @rput reads fnames conds prevalence structural_zeros negative_lower_bound conserve tolerance max_iterations
+        @rput reads fnames conds prevalence pseudocount structural_zeros negative_lower_bound conserve tolerance max_iterations
         R"""
         rownames(reads) <- fnames
         colnames(reads) <- paste0("s", seq_len(ncol(reads)))
-        meta <- data.frame(condition=factor(conds, levels=c("reference", "contrast")),
+        # The 2.12.0 generic-matrix input path drops a one-column data.frame
+        # to a factor while subsetting samples. An unused identifier column
+        # preserves its dimensionality; it is not part of the model formula.
+        meta <- data.frame(sample_id=colnames(reads),
+                           condition=factor(conds, levels=c("reference", "contrast")),
                            row.names=colnames(reads))
-        input <- phyloseq::phyloseq(phyloseq::otu_table(reads, taxa_are_rows=TRUE),
-                                   phyloseq::sample_data(meta))
-        reference <- ANCOMBC::ancombc(data=input, tax_level=NULL, formula="condition",
-            p_adj_method="BH", prv_cut=prevalence, lib_cut=0, group="condition",
+        reference <- ANCOMBC::ancombc(data=reads, aggregate_data=reads,
+            meta_data=meta, tax_level=NULL, formula="condition",
+            p_adj_method="BH", pseudo=pseudocount, prv_cut=prevalence, lib_cut=0, group="condition",
             struc_zero=structural_zeros, neg_lb=negative_lower_bound, tol=tolerance,
             max_iter=max_iterations, conserve=conserve, global=FALSE, n_cl=1, verbose=FALSE)
         ref <- reference$res
@@ -98,7 +103,7 @@ function CompositionalDA.ancombc_reference(counts::AbstractMatrix{<:Integer}, fe
         rcopy(DataFrame, R"native_shape"), rcopy(String, R"reference_version")
     end
     return DAResult(table, Provenance(:ancombc, :r_reference;
-        r_package="ANCOMBC $version", prevalence, structural_zeros, negative_lower_bound,
+        r_package="ANCOMBC $version", prevalence, pseudocount, structural_zeros, negative_lower_bound,
         conserve, tolerance, max_iterations,
         contrast="$(design.levels[2]) vs $(design.levels[1])", log_base="e"))
 end

proofs/PROOF-STATUS.md

@@ -26,7 +26,7 @@ Machine-checked proofs for CompositionalDA.jl, in **Agda 2.6.4.3** with
 | Gate | Command | Result |
 | --- | --- | --- |
 | Type-check every proof, audit, then self-test | `proofs/check.sh` | PASS — exit 0, no warnings |
-| Axiom audit on its own | `proofs/tests/axiom-audit.sh` | PASS — 7/7 modules reachable from `CompositionalDA.All`, no postulates, no FFI, no unsound flags, no holes, all `--safe` |
+| Axiom audit on its own | `proofs/tests/axiom-audit.sh` | PASS — 9/9 modules reachable from `CompositionalDA.All` at the v0.3 receipt below, no postulates, no FFI, no unsound flags, no holes, all `--safe` |
 | Gate self-test (does the gate reject anything?) | `proofs/tests/gate-selftest.sh` | PASS — see "Last verified" below for the count |
 
 The CI job `proofs` in `.github/workflows/ci.yml` runs `proofs/check.sh` in a
@@ -256,6 +256,17 @@ skip** — in CI and locally alike.
 
 ## Last verified
 
+### v0.3 development, 2026-11-07
+
+At source `bace2e6454dcf4977f8ab3dcf9ab0eb276403455`, the
+[pinned proof job](https://github.com/hyperpolymath/CompositionalDA.jl/actions/runs/42396769228/job/128475707706)
+passed: nine reachable modules, a clean axiom audit, full type checking and
+20/20 gate controls. The Julia 1.10 and 1.12.5 jobs each passed 202 assertions,
+including Aqua. The separate ANCOM-BC reference gate failed; these proof and
+unit-test results do not qualify that revision for release.
+
+### Historical v0.2 receipt
+
 2026-10-02, after rebasing onto PR #2 (`320513d`, which added `StepUp` and
 the 18-control self-test), one run of the CI interface with Debian 13's
 packages (`agda-bin 2.6.4.3-1+b2`, `agda-stdlib 2.1-4`), interface files

src/Core/composition.jl

@@ -38,9 +38,15 @@ function validate_two_groups(counts, features, groups)
     any(isempty ∘ strip, f) && throw(CompositionError("feature names cannot be empty"))
     any(isempty ∘ strip, g) && throw(CompositionError("group labels cannot be empty"))
     allunique(f) || throw(CompositionError("feature names must be unique"))
-    levels = sort(unique(g))
+    original_levels = try
+        sort(unique(groups))
+    catch
+        throw(CompositionError("group labels must be mutually sortable"))
+    end
+    levels = string.(original_levels)
+    allunique(levels) || throw(CompositionError("distinct group labels must have distinct display names"))
     length(levels) == 2 || throw(CompositionError("exactly two groups are required"))
-    i1, i2 = (findall(==(level), g) for level in levels)
+    i1, i2 = (findall(==(level), groups) for level in original_levels)
     min(length(i1), length(i2)) >= 2 || throw(CompositionError("each group needs at least two samples"))
     return (; features = f, groups = g, levels, i1, i2)
 end

src/Methods/ALDEx2.jl

@@ -189,8 +189,8 @@ sample's condition against the other. The result has the columns `feature`,
 `effect`, `pvalue`/`qvalue` (consistent-sign expected Welch p and expected BH),
 `wi_pvalue`/`wi_qvalue` (Wilcoxon), `diff_btw`, `diff_win` and `rab_all`.
 
-Deterministic for a fixed `seed` (StableRNG). Calibration is tested against the
-Bioconductor reference, not proved.
+Deterministic for a fixed `seed` (StableRNG). Numerical agreement is tested
+against the Bioconductor reference. This is not a proof of calibration.
 """
 function aldex2(counts::AbstractMatrix{<:Integer}, features::AbstractVector,
                 groups::AbstractVector; mc_samples::Integer = 128, seed::Integer = 1)