diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index f80b6a6..09dc0ed 100644 --- a/.github/workflows/ci.yml +++ b/.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 diff --git a/README.md b/README.md index 9210403..bc115d1 100644 --- a/README.md +++ b/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 diff --git a/docs/ANCOMBC-SPEC.md b/docs/ANCOMBC-SPEC.md index 9fbbe31..d97cf94 100644 --- a/docs/ANCOMBC-SPEC.md +++ b/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. diff --git a/docs/ISSUE-23-ACCEPTANCE.md b/docs/ISSUE-23-ACCEPTANCE.md index 24a9459..fee6306 100644 --- a/docs/ISSUE-23-ACCEPTANCE.md +++ b/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. diff --git a/ext/CompositionalDARCallExt.jl b/ext/CompositionalDARCallExt.jl index a6248c2..9e0c63b 100644 --- a/ext/CompositionalDARCallExt.jl +++ b/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 diff --git a/proofs/PROOF-STATUS.md b/proofs/PROOF-STATUS.md index cfbcca4..f90dc51 100644 --- a/proofs/PROOF-STATUS.md +++ b/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 diff --git a/src/Core/composition.jl b/src/Core/composition.jl index b23c9da..e601fb5 100644 --- a/src/Core/composition.jl +++ b/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 diff --git a/src/Methods/ALDEx2.jl b/src/Methods/ALDEx2.jl index db61994..1d9c0be 100644 --- a/src/Methods/ALDEx2.jl +++ b/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) diff --git a/src/Methods/ANCOMBC.jl b/src/Methods/ANCOMBC.jl index 1bf2db1..6138790 100644 --- a/src/Methods/ANCOMBC.jl +++ b/src/Methods/ANCOMBC.jl @@ -65,6 +65,17 @@ function mixture_responsibilities(x, v, centres, extra, mixing) return weights ./ total, sum(maxima .+ log.(total)) end +function ancombc_wls(x, v, centres, extra, mixing) + # The reference WLS approximation assigns lower/upper tails by the fitted + # component proportions. These are not maximum-posterior classifications. + lower, upper = quantile(x, [mixing[1], 1 - mixing[3]]) + classes = [value < lower ? 1 : value >= upper ? 3 : 2 for value in x] + weights = [1 / (v[i] + extra[classes[i]]) for i in eachindex(x)] + shifts = centres .- centres[2] + return (; wls = sum(weights .* (x .- shifts[classes])) / sum(weights), + variance = 1 / sum(weights), classes) +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")) @@ -103,12 +114,9 @@ function ancombc_bias(x, v; tolerance = 1e-7, max_iterations = 1000) 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), + change <= tolerance * (1 + maximum(abs, vcat(centres, extra, mixing))) + uncertainty = ancombc_wls(x, v, centres, extra, mixing) + return (; bias = centres[2], wls = uncertainty.wls, variance = uncertainty.variance, iterations = iteration, likelihood = next_likelihood, mixing, centres, extra) end @@ -119,7 +127,7 @@ end """ ancombc(counts, features, groups; prevalence = 0.1, pseudocount = 1.0, - structural_zeros = true, negative_lower_bound = false, + structural_zeros = false, 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. @@ -135,7 +143,7 @@ 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, + pseudocount::Real = 1.0, structural_zeros::Bool = false, negative_lower_bound::Bool = false, conserve::Bool = true, tolerance::Real = 1e-7, max_iterations::Integer = 1000) design = validate_two_groups(counts, features, groups) @@ -171,7 +179,9 @@ function ancombc(counts::AbstractMatrix{<:Integer}, features::AbstractVector, 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) + # The reference API's conserve option changes the uncertainty estimate, + # not the fitted coefficient. WLS supplies the bias-variance approximation. + effect[tested] = moments.difference .- 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])) @@ -184,7 +194,7 @@ function ancombc(counts::AbstractMatrix{<:Integer}, features::AbstractVector, 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") + mixture_weights = fit.mixing, bias_classification = "mixture-weight-quantiles", p_adjust = "BH") return DAResult(table, provenance) end diff --git a/test/boundary_aldex2.jl b/test/boundary_aldex2.jl index 66f6c48..490d72e 100644 --- a/test/boundary_aldex2.jl +++ b/test/boundary_aldex2.jl @@ -19,6 +19,10 @@ @test validate_counts(fill(typemax(Int), 2, 2)) == fill(typemax(Int), 2, 2) baseline = run(x, f, g) + numeric = run(x, f, [2, 2, 10, 10]) + @test numeric.table == baseline.table + @test numeric.provenance.parameters[:contrast] == "10 vs 2" + @test_throws CompositionError run(x, f, Any[2, 2, "2", "2"]) withzero = run(hcat(zeros(Int, 4), x), vcat(["absent"], f), g) @test withzero.table == baseline.table @test withzero.provenance.parameters[:excluded_all_zero] == ["absent"] diff --git a/test/method_ancombc.jl b/test/method_ancombc.jl index f30d337..2598f28 100644 --- a/test/method_ancombc.jl +++ b/test/method_ancombc.jl @@ -2,6 +2,18 @@ # SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) @testset "ANCOM-BC original two-group model" begin + @testset "WLS mixture-proportion classes and boundary ties" begin + uncertainty = CompositionalDA.ancombc_wls([-2., -1, 0, 1, 2], ones(5), + [-1., 0, 2], [1., 0, 3], [0.2, 0.6, 0.2]) + @test uncertainty.classes == [1, 2, 2, 2, 3] + @test uncertainty.wls ≈ -2/15 + @test uncertainty.variance ≈ 4/15 + tied = CompositionalDA.ancombc_wls([-1., -1, 0, 1, 1], ones(5), + [-1., 0, 2], [1., 0, 3], [0.25, 0.5, 0.25]) + @test tied.classes == [2, 2, 2, 3, 3] + @test tied.wls ≈ -5/7 + @test tied.variance ≈ 2/7 + end @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]) @@ -27,6 +39,9 @@ @test all(result.table.effect[4:5] .< 0) @test result.provenance.parameters[:iterations] < 1000 @test result.table == ancombc(x, features, groups).table + ordinary = ancombc(x, features, groups; conserve=false) + @test ordinary.table.effect == result.table.effect + @test all(result.table.se .>= ordinary.table.se) perm = reverse(1:d) permuted = ancombc(x[:, perm], features[perm], groups) @@ -35,7 +50,7 @@ # 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) + withzero = ancombc(hcat(x, absent, zeros(Int, n)), vcat(features, ["absent_a", "absent_all"]), groups; structural_zeros=true) row = withzero.table[end, :] @test row.status == "structural_zero" @test row.structural_zero_group1 && !row.structural_zero_group2 diff --git a/test/oracle/ancombc_oracle.jl b/test/oracle/ancombc_oracle.jl index 255b1f2..f5686c7 100644 --- a/test/oracle/ancombc_oracle.jl +++ b/test/oracle/ancombc_oracle.jl @@ -1,23 +1,95 @@ # SPDX-License-Identifier: MPL-2.0 # SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) -using CompositionalDA, DataFrames, RCall, Test +using CompositionalDA, DataFrames, RCall, Statistics, Test include("datasets.jl") +# Diagnostic comparison to an explicitly pinned upstream repair. The ordinary +# job leaves this unset and still tests the unmodified Bioconductor release. +if haskey(ENV, "ANCOMBC_ESTIMATOR_SOURCE") + estimator_source = ENV["ANCOMBC_ESTIMATOR_SOURCE"] + @rput estimator_source + R""" + estimator_environment <- new.env(parent=asNamespace("ANCOMBC")) + sys.source(estimator_source, envir=estimator_environment) + assignInNamespace(".bias_em", estimator_environment$.bias_em, ns="ANCOMBC") + """ + println("REFERENCE OVERRIDE: bias estimator loaded from ", estimator_source) +end + # 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) +const PRISTINE_DEFECT_CONTROL = get(ENV, "ANCOMBC_PRISTINE_DEFECT_CONTROL", "0") == "1" +if PRISTINE_DEFECT_CONTROL + haskey(ENV, "ANCOMBC_ESTIMATOR_SOURCE") && error("the pristine control cannot load a repair") + @test rcopy(String, R"as.character(packageVersion('ANCOMBC'))") == "2.12.0" + # Exact equal-variance optimum: max(residual^2 - known_variance, 0). + # This negative control establishes why the original reference cannot + # qualify the native estimator. An unexpected repair must change the gate. + R""" + boundary_objective <- function(t) -dnorm(1, 0, sqrt(0.99+t), log=TRUE) + boundary_result <- nloptr::neldermead(0, boundary_objective, lower=0)$par + """ + @test rcopy(Float64, R"boundary_result") == 0 + @test rcopy(Bool, R"boundary_objective(0.01) < boundary_objective(boundary_result)") + println("KNOWN-DEFECT CONTROL: unmodified ANCOMBC 2.12.0; three soil comparisons are expected failures, not conformance passes") +end +# ANCOMBC 2.12.0 silently returns its current parameters at max_iter. Its +# original 1000-iteration reference run disagreed on soil. Verify convergence +# explicitly; increasing the work budget is not relaxing numerical tolerances. +const REFERENCE_MAX_ITERATIONS = 20000 +function checked_reference(x, f, g) + R""" + loadNamespace("ANCOMBC") + oracle_bias_fits <- list() + trace(".bias_em", where=asNamespace("ANCOMBC"), print=FALSE, + exit=quote({ + .GlobalEnv$oracle_bias_fits[[length(.GlobalEnv$oracle_bias_fits)+1L]] <- + c(iterations=iterNum, residual=epsilon, + bias=delta_new, left=l1_new, right=l2_new, + extra_left=kappa1_new, extra_right=kappa2_new, + weight_left=pi1_new, weight_null=pi0_new, weight_right=pi2_new, + log_likelihood=sum(log( + pi0_new*dnorm(beta, delta_new, sqrt(nu0)) + + pi1_new*dnorm(beta, delta_new+l1_new, sqrt(nu0+kappa1_new)) + + pi2_new*dnorm(beta, delta_new+l2_new, sqrt(nu0+kappa2_new))))) + })) + """ + reference = try + ancombc_reference(x, f, g; structural_zeros=false, + max_iterations=REFERENCE_MAX_ITERATIONS) + finally + R"untrace('.bias_em', where=asNamespace('ANCOMBC'))" + end + diagnostics = rcopy(Matrix{Float64}, R"do.call(rbind, oracle_bias_fits)") + println("Reference fit (iterations, residual, bias, left, right, extra left/right, weights left/null/right, log likelihood): ", diagnostics) + @test size(diagnostics, 1) > 0 + @test all(diagnostics[:, 1] .< REFERENCE_MAX_ITERATIONS) + @test all(diagnostics[:, 2] .<= 1e-7) + return reference +end @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) + reference = checked_reference(x, f, g) @test Set(native.table.feature) == Set(reference.table.feature) positions = indexin(native.table.feature, reference.table.feature) @test all(!isnothing, positions) + moments = CompositionalDA.ancombc_moments(log.(Float64.(x) .+ 1), + findall(==(sort(unique(g))[1]), g), findall(==(sort(unique(g))[2]), g)) + println(name, " bias diagnostic: native EM=", native.provenance.parameters[:bias], + " WLS=", native.provenance.parameters[:bias_wls], + " log likelihood=", native.provenance.parameters[:log_likelihood], + " reference implied bias=", median(moments.difference .- reference.table[Int.(positions), :effect])) 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] + if PRISTINE_DEFECT_CONTROL && name == "soil" && column in (:effect, :pvalue, :qvalue) + @test_broken gap <= ANCOMBC_TOLERANCE[column] + else + @test gap <= ANCOMBC_TOLERANCE[column] + end end end end