From fea949b74f4bb6b0b14e8ddc76286e7c168b02d8 Mon Sep 17 00:00:00 2001 From: "Jonathan D.A. Jewell" <6759885+hyperpolymath@users.noreply.github.com> Date: Sat, 7 Nov 2026 15:50:01 +0000 Subject: [PATCH] feat: implement two-group ANCOM-BC and strengthen DA verification Correct ALDEx2 log2 units and input boundaries. Add an independent paper-derived ANCOM-BC estimator, reference comparisons, scoped algebraic proofs and explicit release gates for MetaManifold issue 23. Draft pending full conformance and registration. --- .github/workflows/ci.yml | 36 +++- Project.toml | 4 +- README.md | 37 ++-- docs/ALDEx2-SPEC.md | 10 + docs/ANCOMBC-SPEC.md | 76 ++++++++ docs/ISSUE-23-ACCEPTANCE.md | 51 +++++ docs/R-SUNSET.md | 12 +- ext/CompositionalDARCallExt.jl | 61 +++++- proofs/PROOF-STATUS.md | 22 +++ proofs/agda/CompositionalDA/All.agda | 2 + proofs/agda/CompositionalDA/BiasModel.agda | 44 +++++ proofs/agda/CompositionalDA/PosteriorP.agda | 36 ++++ proofs/tests/gate-selftest.sh | 6 + src/CompositionalDA.jl | 6 +- src/Core/composition.jl | 30 ++- src/Methods/ALDEx2.jl | 40 ++-- src/Methods/ANCOMBC.jl | 197 ++++++++++++++++++++ test/boundary_aldex2.jl | 35 ++++ test/method_ancombc.jl | 51 +++++ test/oracle/ancombc_oracle.jl | 23 +++ test/oracle/runoracle.jl | 43 ++++- test/posterior_aggregation.jl | 18 ++ test/runtests.jl | 3 + 23 files changed, 769 insertions(+), 74 deletions(-) create mode 100644 docs/ANCOMBC-SPEC.md create mode 100644 docs/ISSUE-23-ACCEPTANCE.md create mode 100644 proofs/agda/CompositionalDA/BiasModel.agda create mode 100644 proofs/agda/CompositionalDA/PosteriorP.agda create mode 100644 src/Methods/ANCOMBC.jl create mode 100644 test/boundary_aldex2.jl create mode 100644 test/method_ancombc.jl create mode 100644 test/oracle/ancombc_oracle.jl create mode 100644 test/posterior_aggregation.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index e3d1b02..f80b6a6 100644 --- a/.github/workflows/ci.yml +++ b/.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 }} diff --git a/Project.toml b/Project.toml index ccc1d78..94c8554 100644 --- a/Project.toml +++ b/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 "] [deps] diff --git a/README.md b/README.md index 6834607..9210403 100644 --- a/README.md +++ b/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. diff --git a/docs/ALDEx2-SPEC.md b/docs/ALDEx2-SPEC.md index f8bbda6..dbe4332 100644 --- a/docs/ALDEx2-SPEC.md +++ b/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. diff --git a/docs/ANCOMBC-SPEC.md b/docs/ANCOMBC-SPEC.md new file mode 100644 index 0000000..9fbbe31 --- /dev/null +++ b/docs/ANCOMBC-SPEC.md @@ -0,0 +1,76 @@ + + + +# 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.** diff --git a/docs/ISSUE-23-ACCEPTANCE.md b/docs/ISSUE-23-ACCEPTANCE.md new file mode 100644 index 0000000..24a9459 --- /dev/null +++ b/docs/ISSUE-23-ACCEPTANCE.md @@ -0,0 +1,51 @@ + + + +# 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. diff --git a/docs/R-SUNSET.md b/docs/R-SUNSET.md index 552df3c..3d96dfe 100644 --- a/docs/R-SUNSET.md +++ b/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 diff --git a/ext/CompositionalDARCallExt.jl b/ext/CompositionalDARCallExt.jl index eb8b597..a6248c2 100644 --- a/ext/CompositionalDARCallExt.jl +++ b/ext/CompositionalDARCallExt.jl @@ -10,7 +10,8 @@ module CompositionalDARCallExt using CompositionalDA -using CompositionalDA: Provenance, DAResult, validate_counts, CompositionError +using CompositionalDA: Provenance, DAResult, validate_counts, CompositionError, + validate_two_groups, validate_mc using DataFrames using RCall @@ -24,16 +25,18 @@ the same column names. Requires R with the ALDEx2 package installed. function CompositionalDA.aldex2_reference(counts::AbstractMatrix{<:Integer}, features::AbstractVector, groups::AbstractVector; mc_samples::Integer = 128, seed::Integer = 1, runner = f -> f()) - validate_counts(counts) - size(counts, 2) == length(features) || throw(CompositionError("feature names do not match columns")) - levels = sort(unique(string.(groups))) - length(levels) == 2 || throw(CompositionError("aldex2 needs exactly 2 groups")) + design = validate_two_groups(counts, features, groups) + validate_mc(mc_samples, seed) + seed <= typemax(Int32) || throw(ArgumentError("R reference seed must fit Int32")) + levels = design.levels + keep = vec(any(>(0), counts; dims = 1)) + count(keep) >= 2 || throw(CompositionError("ALDEx2 needs at least two nonzero features")) df = runner() do - reads = permutedims(Matrix{Int}(counts)) # ALDEx2 wants features × samples + reads = permutedims(Float64.(counts[:, keep])) # R integers are only 32-bit @rput reads - conds = string.(groups); lv = levels + conds = design.groups; lv = levels @rput conds lv - fnames = string.(features); @rput fnames + fnames = design.features[keep]; @rput fnames R""" set.seed($seed) reads <- as.data.frame(reads); rownames(reads) <- fnames @@ -56,8 +59,48 @@ function CompositionalDA.aldex2_reference(counts::AbstractMatrix{<:Integer}, fea qvalue = res[!, "we_eBH"], wi_pvalue = res[!, "wi_ep"], wi_qvalue = res[!, "wi_eBH"], diff_btw = res[!, "diff_btw"], diff_win = res[!, "diff_win"], rab_all = res[!, "rab_all"]) prov = Provenance(:aldex2, :r_reference; seed, mc_samples, r_package = "ALDEx2 $ver", - contrast = "$(levels[2]) vs $(levels[1])") + contrast = "$(levels[2]) vs $(levels[1])", log_base = 2, + excluded_all_zero = design.features[.!keep]) return DAResult(table, prov) end +function CompositionalDA.ancombc_reference(counts::AbstractMatrix{<:Integer}, features::AbstractVector, + groups::AbstractVector; prevalence::Real=0.1, structural_zeros::Bool=false, + negative_lower_bound::Bool=false, conserve::Bool=true, tolerance::Real=1e-7, + max_iterations::Integer=1000, runner=f -> f()) + design = validate_two_groups(counts, features, groups) + 0 <= prevalence < 1 || throw(ArgumentError("prevalence must be in [0,1)")) + 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 + R""" + rownames(reads) <- fnames + colnames(reads) <- paste0("s", seq_len(ncol(reads))) + meta <- data.frame(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", + 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 + ids <- if ("taxon" %in% colnames(ref$lfc)) ref$lfc$taxon else rownames(ref$lfc) + native_shape <- data.frame(feature=ids, effect=ref$lfc[["conditioncontrast"]], + se=ref$se[["conditioncontrast"]], statistic=ref$W[["conditioncontrast"]], + pvalue=ref$p_val[["conditioncontrast"]], qvalue=ref$q_val[["conditioncontrast"]]) + reference_version <- as.character(packageVersion("ANCOMBC")) + """ + 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, + conserve, tolerance, max_iterations, + contrast="$(design.levels[2]) vs $(design.levels[1])", log_base="e")) +end + end # module diff --git a/proofs/PROOF-STATUS.md b/proofs/PROOF-STATUS.md index f599bb8..cfbcca4 100644 --- a/proofs/PROOF-STATUS.md +++ b/proofs/PROOF-STATUS.md @@ -41,6 +41,28 @@ The exact commands, exit codes and counts of the most recent run are in ## What is proved +### Additional method identities (v0.3 development) + +`CompositionalDA.PosteriorP` proves over exact rationals that the sum of +elementwise minima is bounded above by each separate sum, hence their minimum +(`consistent-sign-bound`). This is finite algebra behind ALDEx2's directional +aggregation, not p-value calibration. `test/posterior_aggregation.jl` mirrors +it and supplies a strict counterexample to averaging per-instance minima. +Division by a common positive Monte-Carlo count is not separately formalised. + +`CompositionalDA.BiasModel`, parameterised by a commutative ring, proves +`observational-shift` (opposite abundance/sampling shifts leave their sum +unchanged), `bias-cancellation` (subtracting the supplied bias removes it) and +`residual-bias` (an estimate leaves exactly the bias estimation error). + +These statements have no postulates. They do not discharge statistical +identification assumptions, verify a mixture optimiser, establish its global +convergence, prove the floating-point implementation equivalent or establish +calibration/FDR control. New mutation controls require rejection of an invalid +aggregation bound and addition of bias where subtraction was required. See +`docs/ANCOMBC-SPEC.md`. Historical receipts below cover their named source +versions only; a new release requires a fresh complete gate run. + Every module is reachable from `proofs/agda/CompositionalDA/All.agda`. ### P1 — the centred log-ratio transform diff --git a/proofs/agda/CompositionalDA/All.agda b/proofs/agda/CompositionalDA/All.agda index 2411273..0515b43 100644 --- a/proofs/agda/CompositionalDA/All.agda +++ b/proofs/agda/CompositionalDA/All.agda @@ -16,3 +16,5 @@ import CompositionalDA.CLR import CompositionalDA.CLR.Instance import CompositionalDA.BenjaminiHochberg import CompositionalDA.BenjaminiHochberg.StepUp +import CompositionalDA.PosteriorP +import CompositionalDA.BiasModel diff --git a/proofs/agda/CompositionalDA/BiasModel.agda b/proofs/agda/CompositionalDA/BiasModel.agda new file mode 100644 index 0000000..342ea15 --- /dev/null +++ b/proofs/agda/CompositionalDA/BiasModel.agda @@ -0,0 +1,44 @@ +-- SPDX-License-Identifier: MPL-2.0 +-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +{-# OPTIONS --safe --without-K #-} + +open import Algebra.Bundles using (CommutativeRing) + +module CompositionalDA.BiasModel {c ℓ} (R : CommutativeRing c ℓ) where + +open CommutativeRing R +open import Relation.Binary.Reasoning.Setoid setoid + +-- The observable linear predictor cannot distinguish an abundance shift +-- from an opposite sampling-fraction shift. An identification assumption is +-- necessary; count data alone cannot discharge it. +observational-shift : ∀ μ d s → (μ + s) + (d - s) ≈ μ + d +observational-shift μ d s = begin + (μ + s) + (d - s) + ≈⟨ +-cong refl (+-comm d (- s)) ⟩ + (μ + s) + ((- s) + d) + ≈⟨ sym (+-assoc (μ + s) (- s) d) ⟩ + ((μ + s) + (- s)) + d + ≈⟨ +-cong (+-assoc μ s (- s)) refl ⟩ + (μ + (s + (- s))) + d + ≈⟨ +-cong (+-cong refl (-‿inverseʳ s)) refl ⟩ + (μ + 0#) + d + ≈⟨ +-cong (+-identityʳ μ) refl ⟩ + μ + d ∎ + +-- Subtraction removes EXACTLY the supplied bias. No claim that the estimated +-- bias equals the population bias is assumed to have been proved here. +bias-cancellation : ∀ μ ν δ → ((μ - ν) + δ) - δ ≈ μ - ν +bias-cancellation μ ν δ = begin + ((μ - ν) + δ) - δ + ≈⟨ +-assoc (μ - ν) δ (- δ) ⟩ + (μ - ν) + (δ + (- δ)) + ≈⟨ +-cong refl (-‿inverseʳ δ) ⟩ + (μ - ν) + 0# + ≈⟨ +-identityʳ (μ - ν) ⟩ + μ - ν ∎ + +-- Residual bias is precisely estimation error, with the opposite sign. +residual-bias : ∀ μ ν δ estimate → + ((μ - ν) + δ) - estimate ≈ (μ - ν) + (δ - estimate) +residual-bias μ ν δ estimate = +-assoc (μ - ν) δ (- estimate) diff --git a/proofs/agda/CompositionalDA/PosteriorP.agda b/proofs/agda/CompositionalDA/PosteriorP.agda new file mode 100644 index 0000000..c715aa2 --- /dev/null +++ b/proofs/agda/CompositionalDA/PosteriorP.agda @@ -0,0 +1,36 @@ +-- SPDX-License-Identifier: MPL-2.0 +-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) +{-# OPTIONS --safe --without-K #-} + +-- The finite aggregation inequality behind ALDEx2's consistent-sign rule. +-- This is algebra over exact rationals, not a calibration/FDR theorem. +module CompositionalDA.PosteriorP where + +open import Data.Nat.Base using (ℕ; zero; suc) +open import Data.Vec.Base using (Vec; []; _∷_; zipWith) +open import Data.Rational.Unnormalised.Base using (ℚᵘ; 0ℚᵘ; _+_; _⊓_; _≤_) +open import Data.Rational.Unnormalised.Properties + using (≤-refl; +-mono-≤; p⊓q≤p; p⊓q≤q; ⊓-glb) + +sum : ∀ {n} → Vec ℚᵘ n → ℚᵘ +sum [] = 0ℚᵘ +sum (x ∷ xs) = x + sum xs + +sum-min≤left : ∀ {n} (xs ys : Vec ℚᵘ n) → + sum (zipWith _⊓_ xs ys) ≤ sum xs +sum-min≤left [] [] = ≤-refl +sum-min≤left (x ∷ xs) (y ∷ ys) = + +-mono-≤ (p⊓q≤p x y) (sum-min≤left xs ys) + +sum-min≤right : ∀ {n} (xs ys : Vec ℚᵘ n) → + sum (zipWith _⊓_ xs ys) ≤ sum ys +sum-min≤right [] [] = ≤-refl +sum-min≤right (x ∷ xs) (y ∷ ys) = + +-mono-≤ (p⊓q≤q x y) (sum-min≤right xs ys) + +-- Multiplication by the common positive 1/mc gives the corresponding +-- inequality of means. This theorem does not assume either vector is random. +consistent-sign-bound : ∀ {n} (xs ys : Vec ℚᵘ n) → + sum (zipWith _⊓_ xs ys) ≤ (sum xs ⊓ sum ys) +consistent-sign-bound xs ys = + ⊓-glb (sum-min≤left xs ys) (sum-min≤right xs ys) diff --git a/proofs/tests/gate-selftest.sh b/proofs/tests/gate-selftest.sh index b2bc322..a63a74d 100755 --- a/proofs/tests/gate-selftest.sh +++ b/proofs/tests/gate-selftest.sh @@ -185,6 +185,12 @@ expect_reject "stepup: monotonicity reversed" run_gate "$NAMESPACE/BenjaminiHoch expect_reject "stepup: scale written as rank/m" run_gate "$NAMESPACE/BenjaminiHochberg/StepUp.agda" \ 'sed -i "s|^scale m k p = p \* mkℚᵘ (+ m) k|scale m k p = p * mkℚᵘ (+ suc k) m|" "$1"' +expect_reject "posterior: minimum replaced by maximum" run_gate "$NAMESPACE/PosteriorP.agda" \ + 'sed -i "s|sum (zipWith _⊓_ xs ys) ≤ (sum xs ⊓ sum ys)|sum xs ≤ (sum xs ⊓ sum ys)|" "$1"' + +expect_reject "bias: correction adds bias instead of subtracting" run_gate "$NAMESPACE/BiasModel.agda" \ + 'sed -i "s|((μ - ν) + δ) - δ ≈ μ - ν|((μ - ν) + δ) + δ ≈ μ - ν|" "$1"' + # --- the axiom audit must reject an unchecked or unsound module ------------- expect_reject "audit: module dropped from the gate entry" run_audit "$NAMESPACE/All.agda" \ diff --git a/src/CompositionalDA.jl b/src/CompositionalDA.jl index 08b1fbb..0073838 100644 --- a/src/CompositionalDA.jl +++ b/src/CompositionalDA.jl @@ -5,8 +5,8 @@ CompositionalDA Compositional differential-abundance methods for count tables, implemented -natively in Julia. Implemented: ALDEx2 (native `aldex2`). Planned: ANCOM-BC -(v0.3, not yet in `src/`). Input is always a table of raw integer counts with +natively in Julia: ALDEx2 (`aldex2`) and original two-group ANCOM-BC (`ancombc`). +See the method specifications and acceptance gates. Input is raw integer counts with samples as rows and features (taxa) as columns. The core layer is exported: closure, the centred log-ratio, zero policies, @@ -30,6 +30,7 @@ include("Core/composition.jl") include("Core/bh.jl") include("Core/result.jl") include("Methods/ALDEx2.jl") +include("Methods/ANCOMBC.jl") export CompositionError export ZeroPolicy, RefuseZeros, Pseudocount, MultiplicativeReplacement, replace_zeros @@ -37,5 +38,6 @@ export validate_counts, closure, clr export bh_adjust export Provenance, DAResult export aldex2, aldex2_reference, dirichlet_clr_instances +export ancombc, ancombc_reference end # module diff --git a/src/Core/composition.jl b/src/Core/composition.jl index 591ab57..b23c9da 100644 --- a/src/Core/composition.jl +++ b/src/Core/composition.jl @@ -15,7 +15,8 @@ function validate_counts(counts::AbstractMatrix{<:Integer}) d > 1 || throw(CompositionError("a composition needs at least 2 features, got $d")) any(<(0), counts) && throw(CompositionError("negative counts present")) for i in 1:n - sum(@view counts[i, :]) > 0 || throw(CompositionError("sample $i has zero total count")) + # Do not sum machine integers: a valid library can overflow Int/UInt. + any(>(0), @view counts[i, :]) || throw(CompositionError("sample $i has zero total count")) end return counts end @@ -23,6 +24,33 @@ end validate_counts(::AbstractMatrix) = throw(CompositionError("counts must be an integer matrix (raw counts, not normalised values)")) +# Both native methods and their reference adapters share this boundary. In +# particular, converting missing labels to the string "missing" is forbidden. +function validate_two_groups(counts, features, groups) + validate_counts(counts) + n, d = size(counts) + length(features) == d || throw(CompositionError("feature names do not match columns")) + length(groups) == n || throw(CompositionError("group labels do not match samples")) + any(ismissing, features) && throw(CompositionError("feature names cannot be missing")) + any(ismissing, groups) && throw(CompositionError("group labels cannot be missing")) + f = string.(features) + g = string.(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)) + length(levels) == 2 || throw(CompositionError("exactly two groups are required")) + i1, i2 = (findall(==(level), g) for level in 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 + +function validate_mc(mc_samples, seed) + mc_samples >= 2 || throw(ArgumentError("mc_samples must be at least 2")) + 0 <= seed <= typemax(UInt64) || throw(ArgumentError("seed must fit UInt64")) + return nothing +end + """ closure(x) -> Matrix{Float64} diff --git a/src/Methods/ALDEx2.jl b/src/Methods/ALDEx2.jl index f63c177..db61994 100644 --- a/src/Methods/ALDEx2.jl +++ b/src/Methods/ALDEx2.jl @@ -19,19 +19,22 @@ ALDEx2 carries sampling uncertainty into the tests. """ function dirichlet_clr_instances(rng::AbstractRNG, counts::AbstractMatrix{<:Integer}; mc_samples::Integer = 128, prior::Real = 0.5) - mc_samples > 0 || throw(ArgumentError("mc_samples must be positive")) - prior > 0 || throw(ArgumentError("prior must be positive")) + validate_counts(counts) + mc_samples >= 2 || throw(ArgumentError("mc_samples must be at least 2")) + isfinite(prior) && prior > 0 || throw(ArgumentError("prior must be finite and positive")) n, d = size(counts) out = Array{Float64,3}(undef, n, d, mc_samples) g = Vector{Float64}(undef, d) for i in 1:n, k in 1:mc_samples for j in 1:d - g[j] = rand(rng, Gamma(counts[i, j] + prior, 1.0)) + shape = Float64(counts[i, j]) + Float64(prior) + isfinite(shape) || throw(CompositionError("Dirichlet shape is not representable in Float64")) + g[j] = rand(rng, Gamma(shape, 1.0)) end - # Gamma draws with tiny shape can underflow to 0; floor them at the - # smallest positive normal so the log stays finite. - g .= max.(g, floatmin(Float64)) - l = log.(g) + # Flooring a zero draw would silently change the posterior. Refuse an + # unrepresentable draw instead. ALDEx2 reports CLR values in log2 units. + all(v -> isfinite(v) && v > 0, g) || throw(CompositionError("Dirichlet draw underflowed or overflowed; no result produced")) + l = log2.(g) out[i, :, k] .= l .- mean(l) end return out @@ -191,23 +194,20 @@ Bioconductor reference, not proved. """ function aldex2(counts::AbstractMatrix{<:Integer}, features::AbstractVector, groups::AbstractVector; mc_samples::Integer = 128, seed::Integer = 1) - validate_counts(counts) - n, d = size(counts) - length(features) == d || throw(CompositionError("$(length(features)) feature names for $d columns")) - length(groups) == n || throw(CompositionError("$(length(groups)) group labels for $n samples")) - levels = sort(unique(groups)) - length(levels) == 2 || throw(CompositionError("aldex2 needs exactly 2 groups, got $(length(levels))")) - i1 = findall(==(levels[1]), groups) - i2 = findall(==(levels[2]), groups) - (length(i1) >= 2 && length(i2) >= 2) || throw(CompositionError("each group needs at least 2 samples")) - first_is_i1 = groups[1] == levels[1] + design = validate_two_groups(counts, features, groups) + validate_mc(mc_samples, seed) + (; levels, i1, i2) = design + first_is_i1 = design.groups[1] == levels[1] + keep = vec(any(>(0), counts; dims = 1)) + count(keep) >= 2 || throw(CompositionError("ALDEx2 needs at least two nonzero features")) rng = StableRNG(seed) - inst = dirichlet_clr_instances(rng, counts; mc_samples) + inst = dirichlet_clr_instances(rng, counts[:, keep]; mc_samples) table = aldex2_from_instances(inst, i1, i2; rng, first_is_i1) - insertcols!(table, 1, :feature => collect(features)) + insertcols!(table, 1, :feature => design.features[keep]) prov = Provenance(:aldex2, :native; seed, mc_samples, test = "welch+wilcoxon", - contrast = "$(levels[2]) vs $(levels[1])", prior = 0.5) + contrast = "$(levels[2]) vs $(levels[1])", prior = 0.5, + log_base = 2, excluded_all_zero = design.features[.!keep]) return DAResult(table, prov) end diff --git a/src/Methods/ANCOMBC.jl b/src/Methods/ANCOMBC.jl new file mode 100644 index 0000000..1bf2db1 --- /dev/null +++ b/src/Methods/ANCOMBC.jl @@ -0,0 +1,197 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +# Independent implementation of the two-group model in Lin & Peddada (2020), +# doi:10.1038/s41467-020-17041-7, equations 6, 14, 19--22, 31, 33 and 35. +# This is not ANCOM-BC2 and does not implement covariates or repeated measures. + +function ancombc_moments(y, i1, i2) + d = size(y, 2) + means = zeros(2, d) + variances = zeros(2, d) + offsets = zeros(size(y, 1)) + for (g, rows) in enumerate((i1, i2)) + block = y[rows, :] + means[g, :] = vec(mean(block; dims = 1)) + offsets[rows] = vec(mean(block; dims = 2)) .- mean(block) + residual = block .- offsets[rows] .- means[g:g, :] + variances[g, :] = vec(sum(abs2, residual; dims = 1)) ./ length(rows)^2 + end + return (; difference = vec(means[2, :] .- means[1, :]), + variance = vec(sum(variances; dims = 1)), offsets) +end + +# Bounded one-dimensional minimisation. The variance objective need not be +# unimodal, so bracket every sampled local minimum, retaining the best result. +function variance_component(x, v, weights, centre) + residual2 = (x .- centre).^2 + upper = maximum(residual2) + upper == 0 && return 0.0 + objective(t) = sum(weights .* (log.(v .+ t) .+ residual2 ./ (v .+ t))) + grid = vcat(0.0, exp.(range(log(minimum(v) * 1e-6), log(max(upper, minimum(v) * 1e-6)); length = 65))) + values = objective.(grid) + best = argmin(values) + answer, score = grid[best], values[best] + golden = (sqrt(5.0) - 1) / 2 + for j in 2:length(grid)-1 + values[j] <= values[j-1] && values[j] <= values[j+1] || continue + lo, hi = grid[j-1], grid[j+1] + a, b = hi - golden * (hi-lo), lo + golden * (hi-lo) + fa, fb = objective(a), objective(b) + for _ in 1:80 + if fa < fb + hi, b, fb = b, a, fa + a = hi - golden * (hi-lo); fa = objective(a) + else + lo, a, fa = a, b, fb + b = lo + golden * (hi-lo); fb = objective(b) + end + end + candidate = (lo + hi) / 2 + value = objective(candidate) + if value < score + answer, score = candidate, value + end + end + return answer +end + +function mixture_responsibilities(x, v, centres, extra, mixing) + logweights = hcat([log(mixing[k]) .- 0.5 .* (log.(2π .* (v .+ extra[k])) .+ + (x .- centres[k]).^2 ./ (v .+ extra[k])) for k in 1:3]...) + maxima = maximum(logweights; dims = 2) + weights = exp.(logweights .- maxima) + total = sum(weights; dims = 2) + return weights ./ total, sum(maxima .+ log.(total)) +end + +function ancombc_bias(x, v; tolerance = 1e-7, max_iterations = 1000) + length(x) >= 4 || throw(CompositionError("ANCOM-BC requires at least four estimable taxa")) + all(z -> isfinite(z) && z > 0, v) || throw(CompositionError("ANCOM-BC has a zero or non-finite variance; no result produced")) + centres = quantile(x, [0.1, 0.5, 0.9]) + extra = [var(x), 0.0, var(x)] + mixing = [0.1, 0.8, 0.1] + responsibilities, likelihood = mixture_responsibilities(x, v, centres, extra, mixing) + for iteration in 1:max_iterations + oldcentres = copy(centres) + oldextra = copy(extra) + oldmixing = copy(mixing) + mixing = vec(mean(responsibilities; dims = 1)) + all(>(eps(Float64)), mixing) || throw(CompositionError("ANCOM-BC mixture component collapsed; no result produced")) + precision = [responsibilities[:, k] ./ (v .+ extra[k]) for k in 1:3] + masses = sum.(precision) + centres = [sum(precision[k] .* x) / masses[k] for k in 1:3] + # Weighted isotonic projection is the constrained maximiser in the + # means for fixed variances: negative <= null <= positive. + if centres[1] > centres[2] + centres[1:2] .= sum(masses[1:2] .* centres[1:2]) / sum(masses[1:2]) + end + if centres[2] > centres[3] + centres[2:3] .= sum(masses[2:3] .* centres[2:3]) / sum(masses[2:3]) + if centres[1] > centres[2] + centres .= sum(masses .* centres) / sum(masses) + end + end + for k in (1, 3) + candidate = variance_component(x, v, responsibilities[:, k], centres[k]) + objective(t) = sum(responsibilities[:, k] .* (log.(v .+ t) .+ (x .- centres[k]).^2 ./ (v .+ t))) + objective(candidate) <= objective(extra[k]) && (extra[k] = candidate) + end + next_resp, next_likelihood = mixture_responsibilities(x, v, centres, extra, mixing) + isfinite(next_likelihood) || throw(CompositionError("ANCOM-BC mixture is non-finite")) + next_likelihood + 1e-9 * (1 + abs(likelihood)) >= likelihood || + throw(CompositionError("ANCOM-BC likelihood decreased; no result produced")) + change = maximum(abs.(vcat(centres-oldcentres, extra-oldextra, mixing-oldmixing))) + if abs(next_likelihood - likelihood) <= tolerance * (1 + abs(likelihood)) && + change <= sqrt(tolerance) * (1 + maximum(abs, centres)) + classes = [argmax(@view next_resp[i, :]) for i in eachindex(x)] + weights = [1 / (v[i] + extra[classes[i]]) for i in eachindex(x)] + shifts = centres .- centres[2] + wls = sum(weights .* (x .- shifts[classes])) / sum(weights) + return (; bias = centres[2], wls, variance = 1 / sum(weights), + iterations = iteration, likelihood = next_likelihood, + mixing, centres, extra) + end + responsibilities, likelihood = next_resp, next_likelihood + end + throw(CompositionError("ANCOM-BC mixture did not converge in $max_iterations iterations; no result produced")) +end + +""" + ancombc(counts, features, groups; prevalence = 0.1, pseudocount = 1.0, + structural_zeros = true, negative_lower_bound = false, + conserve = true, tolerance = 1e-7, max_iterations = 1000) + +Two-group ANCOM-BC (the original 2020 model), on samples × taxa raw counts. +The effect is a natural-log difference, sorted group 2 minus sorted group 1. +BH adjusts the Wald p-values of tested taxa only. Structural-zero decisions +are reported separately with `status` and undefined (`NaN`) numerical tests. +All-zero and below-prevalence taxa are omitted and named in provenance. + +This estimator assumes enough non-differential taxa to identify sampling bias. +It does not make absolute abundance identifiable from counts alone. Calibration, +floating-point execution and equivalence to the R package are not formal proofs. +See `docs/ANCOMBC-SPEC.md` for the exact model and verification boundaries. +""" +function ancombc(counts::AbstractMatrix{<:Integer}, features::AbstractVector, + groups::AbstractVector; prevalence::Real = 0.1, + pseudocount::Real = 1.0, structural_zeros::Bool = true, + negative_lower_bound::Bool = false, conserve::Bool = true, + tolerance::Real = 1e-7, max_iterations::Integer = 1000) + design = validate_two_groups(counts, features, groups) + isfinite(prevalence) && 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")) + negative_lower_bound && !structural_zeros && throw(ArgumentError("negative_lower_bound requires structural_zeros")) + present = counts .> 0 + keep = vec(any(present; dims = 1)) .& (vec(mean(present; dims = 1)) .>= prevalence) + kept = findall(keep) + status = fill("tested", length(kept)) + structural = falses(length(kept), 2) + if structural_zeros + for (g, rows) in enumerate((design.i1, design.i2)), (j, col) in enumerate(kept) + fraction = mean(present[rows, col]) + structural[j, g] = fraction == 0 || (negative_lower_bound && + fraction - 1.959963984540054 * sqrt(fraction * (1-fraction) / length(rows)) <= 0) + end + for j in eachindex(kept) + if all(structural[j, :]) + status[j] = "structural_zero_both" + elseif any(structural[j, :]) + status[j] = "structural_zero" + end + end + end + tested = findall(==("tested"), status) + length(tested) >= 4 || throw(CompositionError("fewer than four taxa remain for ANCOM-BC bias estimation")) + y = log.(Float64.(counts[:, kept[tested]]) .+ Float64(pseudocount)) + all(isfinite, y) || throw(CompositionError("log counts are not representable in Float64")) + moments = ancombc_moments(y, design.i1, design.i2) + fit = ancombc_bias(moments.difference, moments.variance; tolerance, max_iterations) + d = length(kept) + effect, se, statistic, p = (fill(NaN, d) for _ in 1:4) + effect[tested] = moments.difference .- (conserve ? fit.wls : fit.bias) + se[tested] = sqrt.(moments.variance) .+ (conserve ? sqrt(fit.variance) : 0.0) + statistic[tested] = effect[tested] ./ se[tested] + p[tested] = 2 .* ccdf.(Normal(), abs.(statistic[tested])) + table = DataFrame(; feature = design.features[kept], effect, pvalue = p, + qvalue = bh_adjust(p), se, statistic, status, + structural_zero_group1 = structural[:, 1], structural_zero_group2 = structural[:, 2]) + provenance = Provenance(:ancombc, :native; model = "Lin-Peddada-2020-two-group", + contrast = "$(design.levels[2]) vs $(design.levels[1])", log_base = "e", + prevalence, pseudocount, structural_zeros, negative_lower_bound, conserve, + tolerance, max_iterations, excluded_features = design.features[.!keep], + bias = fit.bias, bias_wls = fit.wls, bias_variance = fit.variance, + iterations = fit.iterations, log_likelihood = fit.likelihood, + mixture_weights = fit.mixing, p_adjust = "BH") + return DAResult(table, provenance) +end + +""" + ancombc_reference(counts, features, groups; runner = f -> f(), kwargs...) + +Bioconductor ANCOMBC's original `ancombc` function, through the optional RCall +extension. A numerical reference, never a fallback for the native method. +""" +function ancombc_reference end diff --git a/test/boundary_aldex2.jl b/test/boundary_aldex2.jl new file mode 100644 index 0000000..66f6c48 --- /dev/null +++ b/test/boundary_aldex2.jl @@ -0,0 +1,35 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +@testset "ALDEx2 data boundary and units" begin + x = [2 4 8; 3 6 9; 5 8 3; 6 9 4] + f = ["x", "y", "z"] + g = ["a", "a", "b", "b"] + run(x, f, g) = aldex2(x, f, g; mc_samples = 8) + @test_throws CompositionError run(x, ["x", "x", "z"], g) + @test_throws CompositionError run(x, ["x", missing, "z"], g) + @test_throws CompositionError run(x, f, ["a", missing, "b", "b"]) + @test_throws CompositionError run(x, f, ["a", "a", "b"]) + @test_throws CompositionError run(x, f, ["a", "a", " ", " "]) + @test_throws ArgumentError aldex2(x, f, g; mc_samples = 1) + @test_throws ArgumentError aldex2(x, f, g; seed = -1) + @test_throws ArgumentError dirichlet_clr_instances(StableRNG(1), x; prior = Inf) + @test_throws CompositionError dirichlet_clr_instances(StableRNG(1), -x) + @test_throws CompositionError dirichlet_clr_instances(StableRNG(1), zeros(Int, 4, 3)) + @test validate_counts(fill(typemax(Int), 2, 2)) == fill(typemax(Int), 2, 2) + + baseline = run(x, f, g) + withzero = run(hcat(zeros(Int, 4), x), vcat(["absent"], f), g) + @test withzero.table == baseline.table + @test withzero.provenance.parameters[:excluded_all_zero] == ["absent"] + @test_throws CompositionError run(hcat(x[:, 1], zeros(Int, 4, 2)), f, g) + + # Independent reconstruction of the first draw catches natural-log output: + # effect and test statistics alone cannot detect the wrong log base. + r = StableRNG(12) + gamma = [rand(r, CompositionalDA.Gamma(x[1, j] + 0.5, 1.0)) for j in 1:3] + expected = log2.(gamma) .- mean(log2.(gamma)) + inst = dirichlet_clr_instances(StableRNG(12), x; mc_samples = 2) + @test inst[1, :, 1] ≈ expected atol = 1e-14 + @test baseline.provenance.parameters[:log_base] == 2 +end diff --git a/test/method_ancombc.jl b/test/method_ancombc.jl new file mode 100644 index 0000000..f30d337 --- /dev/null +++ b/test/method_ancombc.jl @@ -0,0 +1,51 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +@testset "ANCOM-BC original two-group model" begin + @testset "equations 6 and 14, known answer and nuisance offsets" begin + y = [1.0 3 5 7; 3 1 7 5; 2 5 4 9; 4 3 8 5] + moments = CompositionalDA.ancombc_moments(y, [1, 2], [3, 4]) + @test moments.difference == [1, 2, 0, 1] + @test moments.variance == [1, 1, 2.5, 2.5] + @test moments.offsets == zeros(4) + shifted = CompositionalDA.ancombc_moments(y .+ [10, -4, 12, -2], [1, 2], [3, 4]) + @test shifted.difference == moments.difference .+ 2 + @test shifted.variance == moments.variance + @test shifted.offsets == [7, -7, 7, -7] + end + rng = StableRNG(2026) + n, d = 20, 40 + base = collect(range(1000, 5000; length = d)) + x = round.(Int, permutedims(base) .* exp.(0.25 .* randn(rng, n, d))) + x[11:20, 1:3] .*= 8 + x[11:20, 4:5] .= max.(1, x[11:20, 4:5] .÷ 4) + features = ["f$j" for j in 1:d] + groups = vcat(fill("a", 10), fill("b", 10)) + result = ancombc(x, features, groups) + @test all(result.table.qvalue[1:5] .< 0.05) + @test all(result.table.effect[1:3] .> 0) + @test all(result.table.effect[4:5] .< 0) + @test result.provenance.parameters[:iterations] < 1000 + @test result.table == ancombc(x, features, groups).table + + perm = reverse(1:d) + permuted = ancombc(x[:, perm], features[perm], groups) + @test permuted.table.effect ≈ result.table.effect[perm] atol = 1e-8 + @test permuted.table.qvalue ≈ result.table.qvalue[perm] atol = 1e-8 + + # Structural absence is a separate decision, not a fabricated p-value. + absent = vcat(zeros(Int, 10), fill(200, 10)) + withzero = ancombc(hcat(x, absent, zeros(Int, n)), vcat(features, ["absent_a", "absent_all"]), groups) + row = withzero.table[end, :] + @test row.status == "structural_zero" + @test row.structural_zero_group1 && !row.structural_zero_group2 + @test isnan(row.pvalue) && isnan(row.qvalue) && isnan(row.effect) + @test withzero.provenance.parameters[:excluded_features] == ["absent_all"] + @test withzero.table.effect[1:d] == result.table.effect + + @test_throws CompositionError ancombc(x, features, groups; max_iterations = 1) + @test_throws ArgumentError ancombc(x, features, groups; pseudocount = NaN) + @test_throws ArgumentError ancombc(x, features, groups; tolerance = Inf) + @test_throws CompositionError ancombc(fill(1, 4, 4), ["a", "b", "c", "d"], [1, 1, 2, 2]) + @test_throws CompositionError ancombc(x[:, 1:3], features[1:3], groups) +end diff --git a/test/oracle/ancombc_oracle.jl b/test/oracle/ancombc_oracle.jl new file mode 100644 index 0000000..255b1f2 --- /dev/null +++ b/test/oracle/ancombc_oracle.jl @@ -0,0 +1,23 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +using CompositionalDA, DataFrames, RCall, Test +include("datasets.jl") + +# Fixed BEFORE examining the comparison results. These are numerical +# conformance thresholds, not confidence bounds or proofs of calibration. +const ANCOMBC_TOLERANCE = (effect=0.05, se=0.02, pvalue=0.02, qvalue=0.02) +@testset "ANCOM-BC Bioconductor conformance" begin + for (name, (x, f, g, _)) in [("mock", mock_dataset()), ("gut", fixture_dataset("gut")), ("soil", fixture_dataset("soil"))] + native = ancombc(x, f, g; structural_zeros=false) + reference = ancombc_reference(x, f, g; structural_zeros=false) + @test Set(native.table.feature) == Set(reference.table.feature) + positions = indexin(native.table.feature, reference.table.feature) + @test all(!isnothing, positions) + for column in keys(ANCOMBC_TOLERANCE) + gap = maximum(abs.(native.table[!, column] .- reference.table[Int.(positions), column])) + println(name, " ANCOM-BC ", column, " maximum error = ", gap) + @test gap <= ANCOMBC_TOLERANCE[column] + end + end +end diff --git a/test/oracle/runoracle.jl b/test/oracle/runoracle.jl index efb8198..f9c6c51 100644 --- a/test/oracle/runoracle.jl +++ b/test/oracle/runoracle.jl @@ -15,7 +15,9 @@ include("datasets.jl") const MC = 512 const SLACK = 2.0 -const FLOOR = (effect = 0.05, pvalue = 0.01, wi_pvalue = 0.01) +const FLOOR = (effect = 0.05, pvalue = 0.01, qvalue = 0.01, + wi_pvalue = 0.01, wi_qvalue = 0.01, + diff_btw = 0.05, diff_win = 0.05, rab_all = 0.05) """Largest absolute difference of column `c` between two result tables.""" maxdiff(a, b, c) = maximum(abs.(a.table[!, c] .- b.table[!, c])) @@ -28,7 +30,7 @@ expected Wilcoxon p, and whether native is within tolerance of R for each. """ function agree(native, r1, r2) - out = map((:effect, :pvalue, :wi_pvalue)) do c + out = map(keys(FLOOR)) do c dn = maxdiff(native, r1, c); dr = maxdiff(r2, r1, c) c => (native = dn, r_vs_r = dr, ok = dn <= SLACK * dr + FLOOR[c]) end @@ -43,9 +45,9 @@ end r2 = aldex2_reference(x, f, g; mc_samples = MC, seed = 2) a = agree(nat, r1, r2) println(name, ": ", a, " R=", r1.provenance.parameters) - @test a.effect.ok - @test a.pvalue.ok - @test a.wi_pvalue.ok + for column in keys(FLOOR) + @test a[column].ok + end @test cor(nat.table.effect, r1.table.effect) > 0.99 if !isempty(truth) hit(t) = Set(t.table.feature[t.table.qvalue .< 0.05]) @@ -54,3 +56,34 @@ end end end end + +# The p-value and BH aggregations are deterministic conditional on the draws. +# Check identical R draws independently of both random-number generators and +# of the deliberately different effect-size pairing correlations. +@testset "ALDEx2 identical-instance oracle" begin + for (name, (x, f, g, _)) in [("mock", mock_dataset()), ("gut", fixture_dataset("gut")), ("soil", fixture_dataset("soil"))] + reads = permutedims(Float64.(x)); fnames = string.(f); conds = string.(g) + @rput reads fnames conds + R""" + set.seed(73) + rownames(reads) <- fnames + colnames(reads) <- paste0("s", seq_len(ncol(reads))) + clr_reference <- ALDEx2::aldex.clr(reads, conds, mc.samples=64, denom="all", verbose=FALSE) + test_reference <- ALDEx2::aldex.ttest(clr_reference, verbose=FALSE) + draws <- ALDEx2::getMonteCarloInstances(clr_reference) + shared_draws <- array(unlist(draws), dim=c(nrow(reads), 64L, ncol(reads))) + """ + inst = permutedims(rcopy(Array{Float64,3}, R"shared_draws"), (3, 1, 2)) + ref = rcopy(DataFrame, R"test_reference") + levels = sort(unique(g)) + i1, i2 = (findall(==(level), g) for level in levels) + native = CompositionalDA.aldex2_from_instances(inst, i1, i2; + rng=StableRNG(73), first_is_i1=g[1] == levels[1]) + for (julia_col, r_col) in [(:pvalue, :we_ep), (:qvalue, :we_eBH), + (:wi_pvalue, :wi_ep), (:wi_qvalue, :wi_eBH)] + gap = maximum(abs.(native[!, julia_col] .- ref[!, r_col])) + println(name, " shared draws ", julia_col, " maximum error = ", gap) + @test gap <= 1e-12 + end + end +end diff --git a/test/posterior_aggregation.jl b/test/posterior_aggregation.jl new file mode 100644 index 0000000..9fcc05d --- /dev/null +++ b/test/posterior_aggregation.jl @@ -0,0 +1,18 @@ +# SPDX-License-Identifier: MPL-2.0 +# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) + +@testset "Consistent-sign aggregation (exact rational mirror)" begin + # Opposite posterior directions are the load-bearing counterexample: + # averaging each instance's two-sided p-value is not ALDEx2's rule. + g = [1//10, 1//1] + l = [1//1, 1//10] + @test min(sum(g), sum(l)) / 2 == 11//20 + @test sum(min.(g, l)) / 2 == 1//10 + @test sum(min.(g, l)) < min(sum(g), sum(l)) + rng = StableRNG(98) + for _ in 1:100 + left = rand(rng, -20:20, 12) .// 20 + right = rand(rng, -20:20, 12) .// 20 + @test sum(min.(left, right)) <= min(sum(left), sum(right)) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 449032b..901425f 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -18,4 +18,7 @@ include("mutants.jl") include("property_clr.jl") include("property_bh.jl") include("property_aldex2.jl") + include("boundary_aldex2.jl") + include("method_ancombc.jl") + include("posterior_aggregation.jl") end