ALDEx2 from a behavioural spec; R oracle agrees within MC noise (v0.2.0)

.github/workflows/ci.yml

@@ -43,10 +43,50 @@ jobs:
     steps:
       - name: Install the pinned Debian toolchain
         run: |
-          apt-get update
-          apt-get install --no-install-recommends -y ca-certificates git agda-bin=2.6.4.3-1+b2 agda-stdlib=2.1-4
+          set -euo pipefail
+          curl -fsSL -o "$RUNNER_TEMP/agda.tar.xz" \
+            https://github.com/agda/agda/releases/download/v2.7.0.1/Agda-v2.7.0.1-linux.tar.xz
+          echo "${AGDA_SHA256}  $RUNNER_TEMP/agda.tar.xz" | sha256sum --check
+          tar -xJf "$RUNNER_TEMP/agda.tar.xz" -C "$RUNNER_TEMP"
+          echo "Agda_datadir=$RUNNER_TEMP/Agda-v2.7.0.1/data" >> "$GITHUB_ENV"
+          echo "AGDA_BIN=$RUNNER_TEMP/Agda-v2.7.0.1/bin/agda" >> "$GITHUB_ENV"
+      - name: Type-check and axiom audit
+        run: proofs/bootstrap.sh
+      - name: Gate self-test (the gate must reject planted defects)
+        run: proofs/tests/gate-selftest.sh
+  oracle:
+    name: R reference oracle (R 4.5.0, Bioconductor 3.22)
+    runs-on: ubuntu-24.04
+    timeout-minutes: 60
+    env:
+      R_LIBS_USER: ${{ github.workspace }}/.r-lib
+    steps:
       - uses: actions/checkout@11d5960a326750d5838078e36cf38b85af677262 # v4
+      - uses: r-lib/actions/setup-r@f9a764fea8d5c63df6ef9a5c7795bf7deb5d7e05 # v2
         with:
-          persist-credentials: false
-      - name: Proof gate (axiom audit, type-check, reject controls, gate self-test)
-        run: proofs/check.sh
+          r-version: '4.5.0'
+          use-public-rspm: true
+      - uses: actions/cache@0057852bfaa89a56745cba8c7296529d2fc39830 # v4
+        with:
+          path: .r-lib
+          key: r-4.5.0-bioc-3.22-oracle-v1
+      - name: Install the Bioconductor 3.22 references
+        run: |
+          mkdir -p "$R_LIBS_USER"
+          Rscript -e 'options(Ncpus = 4, 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")'
+      - uses: julia-actions/setup-julia@4c0cb0fce8556fdb04a90347310e5db8b1f98fb9 # v2
+        with:
+          version: '1.12'
+      - uses: julia-actions/cache@d10a6fd8f31b12404a54613ebad242900567f2b9 # v2
+      - name: Native vs reference agreement
+        run: |
+          export R_HOME="$(R RHOME)"
+          # Skip Pkg's bulk precompile: it reports a failed extension without its error.
+          # Loading directly below surfaces the real error if there is one.
+          JULIA_PKG_PRECOMPILE_AUTO=0 julia --project=test/oracle -e 'using Pkg; Pkg.develop(path="."); Pkg.instantiate()'
+          julia --project=test/oracle -e 'using CompositionalDA, RCall; @assert Base.get_extension(CompositionalDA, :CompositionalDARCallExt) !== nothing'
+          julia --project=test/oracle test/oracle/runoracle.jl

.gitignore

@@ -1,3 +1,4 @@
 /Manifest.toml
 *.agdai
 _build/
+/test/oracle/Manifest.toml

Project.toml

@@ -1,28 +1,39 @@
 name = "CompositionalDA"
 uuid = "a77efed1-127c-5b9d-bf19-e0eccfa50304"
+version = "0.2.0"
 authors = ["Jonathan D.A. Jewell <j.d.a.jewell@open.ac.uk>"]
-version = "0.1.0"
 
 [deps]
 DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+HypothesisTests = "09f84164-cd44-5f33-b23f-e6b0d136a0d5"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
+[weakdeps]
+RCall = "6f49c342-dc21-5d91-9882-a32aef131414"
+
+[extensions]
+CompositionalDARCallExt = "RCall"
+
 [compat]
 Aqua = "0.8"
-StableRNGs = "1"
-Test = "1"
 DataFrames = "1.6"
+Distributions = "0.25"
+HypothesisTests = "0.11, 0.12"
 LinearAlgebra = "1"
 Random = "1"
+RCall = "0.14"
+StableRNGs = "1"
 Statistics = "1"
+Test = "1"
 julia = "1.10"
 
 [extras]
 Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "StableRNGs", "Test"]
+test = ["Aqua", "Test"]

README.md

@@ -10,9 +10,11 @@ live in a separate package,
 It needs no Python and no TensorFlow. R is needed only for the optional reference
 backend.
 
-> **Status: v0.1.0, core layer only.** Closure, the centred log-ratio, zero
-> policies, Benjamini–Hochberg and the result/provenance types are implemented.
-> The methods land one per release: ALDEx2 (v0.2), then ANCOM-BC (v0.3).
+> **Status: v0.2.0.** The core layer (closure, the centred log-ratio, zero
+> policies, Benjamini–Hochberg, result/provenance types) and two-group ALDEx2
+> (`aldex2`) are implemented. 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. ANCOM-BC lands in v0.3.
 > Nothing here claims a method that is not yet in `src/`.
 
 ## Input contract

docs/ALDEx2-SPEC.md

@@ -0,0 +1,156 @@
+<!-- SPDX-License-Identifier: CC-BY-SA-4.0 -->
+<!-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) -->
+
+# ALDEx2: behavioural specification
+
+This is the specification the native `aldex2` is written from. It describes
+**what** the reported numbers are, as measurable behaviour, so that an
+implementation can be written from it alone.
+
+## Provenance
+
+The method is from Fernandes et al. 2013 (PLoS ONE 8:e67019) and 2014
+(Microbiome 2:15). The papers do not pin down two details that decide the
+reported numbers:
+
+- how the per-instance p-values are combined into an expected p-value;
+- how the random pairs behind the effect size are drawn.
+
+Bioconductor ALDEx2's public documentation describes the first only loosely.
+NEWS for 1.33 says "all p-values calculated are now posterior p-values with
+consistent sign", and the vignette calls `we.ep` "a posterior predictive
+p-value".
+
+On 2026-10-02 the native and Bioconductor 1.42.0 results disagreed beyond
+Monte-Carlo noise. To find out why, the coding agent (Claude, session
+https://claude.ai/code/session_019Y6e175rh1LbrsEhQvCrJ1) read the reference
+package's R source (GPL ≥ 3): `aldex.ttest`, `t.fast`, `wilcox.fast` and
+`aldex.effect`. No code was copied. The behaviour below was then confirmed
+numerically on the reference package's own Monte-Carlo instances: on the mock,
+gut and soil datasets, the p-value rules reproduce its `we.ep`, `we.eBH` and
+`wi.ep` to within 1e-14.
+
+The owner ruled that the native code is written from this document by someone
+who has **not** read the reference source. That is what happened: a separate
+agent that was barred from the reference source, R and web search wrote
+`src/Methods/ALDEx2.jl` from this document (spec committed alone as `87a6316`,
+before the code).
+
+## Notation
+
+- 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.
+- `inst[i, j, k]` is the CLR value of sample `i`, feature `j`, in Monte-Carlo
+  instance `k = 1..mc`. These are produced by the existing
+  `dirichlet_clr_instances`, which is unchanged by this document.
+- **Test group 1** is the condition of the *first sample in input order*.
+  **Test group 2** is the other condition. With `n1` and `n2` samples
+  respectively, `n1 + n2 = n`.
+
+**Calling convention.** The function is `aldex2_from_instances(inst, i1, i2; rng, first_is_i1 = true)`.
+
+- `i1` and `i2` are the sample rows of `L1` and `L2`.
+- `first_is_i1` says whether the first input sample belongs to `L1`. If it does,
+  `i1` is test group 1; if not, `i2` is.
+- `aldex2` computes `first_is_i1` from `groups[1]`.
+- The property test P6 calls the function without this flag and gets the
+  default.
+
+Test group 1 only matters for the exact Wilcoxon p-value, because that
+distribution is discrete. The rule is part of the reference behaviour: on the
+soil dataset, whose first sample is in `L2`, using `L1` instead moves `wi_pvalue`
+by up to 0.0097.
+
+## Per-instance one-sided p-values
+
+For each instance `k` and feature `j`, let `a` be the values of test group 1 and
+`b` the values of test group 2.
+
+- **Welch.** `t = (mean(a) − mean(b)) / sqrt(s_a²/n1 + s_b²/n2)`, using sample
+  variances. The degrees of freedom are Welch–Satterthwaite. The upper-tail
+  p-value is `pW = P(T_df ≥ t)`.
+- **Wilcoxon rank sum.** `U = #{(x, y) : x ∈ a, y ∈ b, x > y}`.
+  - If `n1 < 50` **and** `n2 < 50`: `pR = P(U' ≥ U)` under the exact null
+    distribution of the Mann–Whitney statistic for sizes `(n1, n2)`.
+  - Otherwise use the normal approximation **without continuity correction**:
+    `pR = 1 − Φ((U − n1·n2/2) / sqrt(n1·n2·(n1+n2+1)/12))`.
+  - CLR instances are continuous, so ties have probability zero. If any do
+    occur, count a tie as ½ in `U`; the result in that case is not
+    oracle-checked.
+
+## Consistent-sign expected p-values
+
+For a per-instance upper-tail p-value `p` (either `pW` or `pR`), define the two
+doubled directional values:
+
+```
+g = min(1, 2p)          # "group 1 higher"
+l = min(1, 2(1 − p))    # "group 1 lower"   (note: 1 − p, not the lower tail)
+```
+
+Then, per test:
+
+- **Expected p:** `min( mean_k g_k , mean_k l_k )`. This is reported as
+  `pvalue` for Welch and `wi_pvalue` for Wilcoxon.
+- **Expected BH:** within each instance, BH-adjust the vector of `g` over
+  features, and separately the vector of `l` over features, using the existing
+  `bh_adjust`. Average each over instances, then take the elementwise minimum of
+  the two averages. This is reported as `qvalue` and `wi_qvalue`.
+
+This is not the mean of two-sided p-values: by Jensen, the mean of the minimum
+is at most the minimum of the means. The two agree when a feature's direction is
+the same in every instance.
+
+## Effect size
+
+The effect is computed on the **sorted** levels, so its sign is `L2 − L1`. The
+test-group rule does not apply here.
+
+For each feature `j`:
+
+- `x1` pools all instance values of `L1` samples, `n_L1·mc` values in total.
+  `x2` pools those of `L2` samples.
+- `m = min(length(x1), length(x2))`.
+- **Within-group spread.** For each group `g`, take two independent uniformly
+  random permutations of `x_g` and keep the first `m` elements of each. `win_g`
+  is their elementwise absolute difference.
+- **Between-group difference.** Take one further independent random
+  permutation of `x1` and one of `x2`, each truncated to `m`. Then
+  `btw = perm(x2) − perm(x1)`.
+  - These draws are **independent** of the within-group draws. Reusing the same
+    elements in `btw` and `win` biases the effect: on soil, the gap was 0.18
+    against noise of 0.03.
+- `winmax = max.(win_1, win_2)`.
+- `ratio = btw ./ winmax`. Any `0/0` becomes `0`; a nonzero value divided by 0
+  stays `±Inf`.
+- The reported values are:
+  - `effect = median(ratio)`
+  - `diff_btw = median(btw)`
+  - `diff_win = median(winmax)`
+  - `rab_all = median` of every value of feature `j` over all samples and
+    instances.
+
+**Deliberate difference from the reference, required by P6.** The permutations
+are drawn as **index** permutations once and shared by every feature. Relabelling
+the features then permutes every output column exactly. The reference draws them
+per feature. The marginal distribution of each feature's output is the same
+under both; only the correlation across features differs. Sharing gives exact
+permutation-equivariance.
+
+## Determinism
+
+All randomness comes from one `StableRNG(seed)`. The draws happen in this order:
+
+1. the Dirichlet instances;
+2. the effect-size index permutations.
+
+## What is checked, and how
+
+- **Exact:** given the same instances, the p-value columns must match the
+  reference to floating-point precision. A diagnostic run on the reference's own
+  instances reproduced this.
+- **Statistical:** with independent Monte-Carlo draws, the oracle requires
+  `|native − R| ≤ 2 × |R(seed 1) − R(seed 2)| + floor` per column
+  (`test/oracle/runoracle.jl`).
+- **Not proved:** calibration; see `PROOF-STATUS.md` and the README.

docs/DATASETS.md

@@ -0,0 +1,50 @@
+<!-- SPDX-License-Identifier: CC-BY-SA-4.0 -->
+<!-- SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) -->
+
+# Oracle datasets
+
+The R-reference oracle (`test/oracle/runoracle.jl`) compares each native method
+with its Bioconductor implementation on three two-group count tables.
+
+| Name | Kind | Samples | Features | Groups | Known answer |
+|---|---|---|---|---|---|
+| `mock` | synthetic, generated in `test/oracle/datasets.jl` | 20 | 40 | a 10 / b 10 | yes: features 1–3 up 8×, 4–5 down 4× |
+| `gut` | stool 16S, Parkinson's disease vs healthy controls | 108 | 202 | H 55 / PAR 53 | no |
+| `soil` | blueberry soil vs rhizosphere 16S | 63 | 300 | Soil 39 / Rhizosphere 24 | no |
+
+`mock` is not a MetaManifold fixture; that repository holds no count-table
+fixtures. Real tables have no ground truth, so for `gut` and `soil` the oracle
+checks agreement with R only.
+
+## Source and licence of `gut` and `soil`
+
+Both come from the deposit accompanying Nearing et al. (2022), which compared
+differential-abundance methods on 38 two-group 16S datasets:
+
+- Nearing J (2021). *16S_rRNA_Microbiome_Datasets*. figshare.
+  doi:[10.6084/m9.figshare.14531724.v1](https://doi.org/10.6084/m9.figshare.14531724.v1).
+  Licence **CC BY 4.0** (as stated by the figshare record).
+- Nearing JT et al. (2022). Microbiome differential abundance methods produce
+  different results across 38 datasets. *Nat Commun* 13:342.
+  doi:10.1038/s41467-022-28034-z.
+
+File: `DA_COMPARE_DATA_21_03_08.tar.gz`, 253,151,868 bytes, md5
+`3bf04686142ba000b3b6f7825cc1a606` (as published by figshare), sha256
+`53337642e26735fc138ed17089f6a467135719413bafbd1f9e5bff048c164b84`.
+Members used: `Hackathon/Studies/par_scheperjans/` and
+`Hackathon/Studies/Blueberry/`, the raw `*_ASVs_table.tsv` (not the rarefied
+`*_rare.tsv`) and the metadata TSV (`comparison` column).
+
+**Changes made** (CC BY 4.0 §3(a)(1)(B)): samples intersected with the
+metadata; features kept if present in ≥ 25% of samples, then at most the 300
+most prevalent; transposed to the package's samples × features convention at
+load time. `test/oracle/make_fixtures.jl` reproduces the fixtures exactly from
+the checksummed tarball. Each fixture carries a REUSE `.license` sidecar.
+
+The deposit's CC BY 4.0 covers these processed tables. The licences of the
+original studies behind each table were not separately verified, and the
+original-study citations are not recorded here because they were not confirmed
+from the deposit itself.
+
+The fixture data are CC BY 4.0, not MPL-2.0; they are test data only and are
+not part of the package's loaded code.

docs/R-SUNSET.md

@@ -3,15 +3,17 @@
 
 # Sunset criterion for the R reference backend
 
-> **Status:** the extension described here is planned and not yet in the tree:
-> there is no `ext/` directory and no `test/oracle_*.jl` yet.
+> **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.
 
 The R backend (`ext/CompositionalDARCallExt.jl`, loaded only when RCall is
-present) is a bridge, not a product. It never owns R state. Callers pass an
-`evaluator` function, so a host application keeps its own R session and lock.
+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 method, once the native implementation agrees with the R
-reference within the tolerances stated in `test/oracle_*.jl` on all three
+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
 further maintenance. When no method still needs it, the extension is deleted.