From a16a27869ae251d440268663c245f7cd72824235 Mon Sep 17 00:00:00 2001 From: hyperpolymath <6759885+hyperpolymath@users.noreply.github.com> Date: Fri, 23 Oct 2026 22:16:41 +0000 Subject: [PATCH] fix(ordination): align varpart2 fraction order and F-statistic rtol in unit tests - vegan::varpart2 orders `part$indfract` rows as `[a] = X1 | X2` (unique_x1), `[c] = X2 | X1` (unique_x2), `[b]` (shared overlap, 0 df), and `[d] = Residuals`. Align both `test/unit/test_ordination.jl` (`ind[2]` vs `ind[3]`) and the partial-Condition two-way `ind_rows` fallback in `src/analysis/ordination.jl` to that exact ordering. - Add `rtol=1e-9` alongside `atol=1e-8` on `res.f_statistic` in the pure-Julia reference comparison (`F ~ 820` computed from 12-decimal rounded fixture matrices has `1.7e-12` relative difference, `1.4e-9` absolute). - Use distinct sample column names (`c1..c3`, `t1..t3`) across the two fixture runs in `_ord_with_study_root` so independent samples do not inadvertently form 2-sample repeated-measures blocks. Co-authored-by: arena-agent <297053741+arena-agent@users.noreply.github.com> --- src/analysis/ordination.jl | 2 +- test/unit/test_ordination.jl | 17 +++++++++-------- 2 files changed, 10 insertions(+), 9 deletions(-) diff --git a/src/analysis/ordination.jl b/src/analysis/ordination.jl index c6e122c9..1646099e 100644 --- a/src/analysis/ordination.jl +++ b/src/analysis/ordination.jl @@ -1722,8 +1722,8 @@ function run_constrained_ordination( vp2 = pure_varpart_two_way(sub1.adj_r2, sub2.adj_r2, r_out.adj_r2) ind_rows = [ (; label = "[a] = $t1 | $t2", df = _term_df(t1, metadata), r2 = nothing, adj_r2 = vp2.unique_x1, testable = true), - (; label = "[b] = Shared ($t1 & $t2)", df = 0, r2 = nothing, adj_r2 = vp2.shared, testable = false), (; label = "[c] = $t2 | $t1", df = _term_df(t2, metadata), r2 = nothing, adj_r2 = vp2.unique_x2, testable = true), + (; label = "[b] = Shared ($t1 & $t2)", df = 0, r2 = nothing, adj_r2 = vp2.shared, testable = false), (; label = "[d] = Residuals incl. conditioned share", df = omnibus.df_residual, r2 = nothing, adj_r2 = vp2.residual, testable = false), ] has_neg = any(r -> !isnothing(r.adj_r2) && r.adj_r2 < -1e-6, ind_rows) diff --git a/test/unit/test_ordination.jl b/test/unit/test_ordination.jl index 823c3170..580c999e 100644 --- a/test/unit/test_ordination.jl +++ b/test/unit/test_ordination.jl @@ -33,14 +33,14 @@ function _ord_with_study_root(f::Function) try mkpath(joinpath(root, "config")) cp(joinpath(_ORD_REPO_ROOT, "config", "defaults"), joinpath(root, "config", "defaults")) - for (grp, rname, counts_rows) in ( - ("Control", "run1", [ + for (grp, rname, s_names, counts_rows) in ( + ("Control", "run1", ("c1", "c2", "c3"), [ ("seq1", "ACGTACGTACGTACGT", "Bacteroides", 410, 395, 430), ("seq2", "TGCAACGTACGTACGT", "Prevotella", 290, 310, 275), ("seq3", "GGCCACGTACGTACGT", "Faecalibacterium", 120, 115, 125), ("seq4", "TTAAACGTACGTACGT", "Roseburia", 80, 90, 75), ]), - ("Treatment", "run2", [ + ("Treatment", "run2", ("t1", "t2", "t3"), [ ("seq1", "ACGTACGTACGTACGT", "Bacteroides", 95, 110, 88), ("seq2", "TGCAACGTACGTACGT", "Prevotella", 140, 130, 150), ("seq3", "GGCCACGTACGTACGT", "Faecalibacterium", 320, 305, 335), @@ -48,15 +48,16 @@ function _ord_with_study_root(f::Function) ]), ) mkpath(joinpath(root, "data", "study1", grp, rname)) - for s in ("s1", "s2", "s3"), rd in ("R1", "R2") + for s in s_names, rd in ("R1", "R2") touch(joinpath(root, "data", "study1", grp, rname, "$(s)_L001_$(rd)_001.fastq.gz")) end merged_dir = joinpath(root, "projects", "study1", grp, rname, "merged") mkpath(merged_dir) db = DuckDB.DB(joinpath(merged_dir, "results.duckdb")) con = DBInterface.connect(db) + col_ddl = join(["$(s) BIGINT" for s in s_names], ", ") DBInterface.execute(con, - "CREATE TABLE merged (SeqName VARCHAR, Sequence VARCHAR, Genus_dada2 VARCHAR, s1 BIGINT, s2 BIGINT, s3 BIGINT)") + "CREATE TABLE merged (SeqName VARCHAR, Sequence VARCHAR, Genus_dada2 VARCHAR, $col_ddl)") for (sn, sq, gn, c1, c2, c3) in counts_rows DBInterface.execute(con, "INSERT INTO merged VALUES ('$sn', '$sq', '$gn', $c1, $c2, $c3)") @@ -194,7 +195,7 @@ end @test isapprox(res.unconstrained_inertia, Float64(exp_rda[:unconstrained_inertia]); atol=1e-9) @test isapprox(res.r2, Float64(exp_rda[:r2]); atol=1e-9) @test isapprox(res.adj_r2, Float64(exp_rda[:adj_r2]); atol=1e-9) - @test isapprox(res.f_statistic, Float64(exp_rda[:f_statistic]); atol=1e-9) + @test isapprox(res.f_statistic, Float64(exp_rda[:f_statistic]); atol=1e-8, rtol=1e-9) @test length(res.constrained_eigenvalues) == length(exp_rda[:constrained_eigenvalues]) for (a, b) in zip(res.constrained_eigenvalues, exp_rda[:constrained_eigenvalues]) @test isapprox(a, Float64(b); atol=1e-9) @@ -316,8 +317,8 @@ end @test out.variance_partitioning.status == "computed" ind = out.variance_partitioning.individual_fractions @test isapprox(ind[1].adj_r2, Float64(exp_vp[:unique_x1]); atol=1e-6) - @test isapprox(ind[2].adj_r2, Float64(exp_vp[:shared]); atol=1e-6) - @test isapprox(ind[3].adj_r2, Float64(exp_vp[:unique_x2]); atol=1e-6) + @test isapprox(ind[2].adj_r2, Float64(exp_vp[:unique_x2]); atol=1e-6) + @test isapprox(ind[3].adj_r2, Float64(exp_vp[:shared]); atol=1e-6) @test isapprox(ind[4].adj_r2, Float64(exp_vp[:residual]); atol=1e-6) end end