Skip to content

Commit 592ddd7

Browse files
hyperpolymathclaudeJoshuaJewell
authored
fix(permanova): test each term, report permutations, refuse untestable terms (#35)
adonis2 in the pinned vegan 2.7-3 defaults to by = NULL: one test of the whole model. run_permanova read row 1 as the group term and a comment called the rows sequential terms, so the reported R2/F/p mixed group and run and the per-term list held a single "Model" row. adonis2 is now called with by = "terms"; the top-level result names its term (group). With permutations restricted within individuals (how(blocks = ...)), permute enumerates every permutation when there are fewer than 999: six individuals of two samples give 63, so the smallest attainable p is 1/64. That count was never reported. It is now (`permutations`, `min_p_value`), read from the control adonis2 used. A term constant within every block - group, when each individual belongs to one group - never has its labels exchanged by within-block permutation, so the p-value vegan prints for it does not test a group difference (measured: 0.03125 with a real effect, 0.0625 on noise; it is not fixed at 1). Its p-value is withheld with the reason; R2 and F are still reported, and terms that vary within blocks (run) are tested as before. Tests (test_permanova.jl, against direct vegan::adonis2 calls with the same seed): per-term names, R2, F and p match by = "terms"; 999 permutations unblocked; 63 permutations and floor 1/64 blocked; group's p withheld with its reason while run's p and F match vegan. 184/184 with test_determinism. Reverting by = "terms" and the refusal: 8 fail, 4 error. Frontend shows the reason and the permutation count; typecheck clean, 27/27. Co-authored-by: Claude Opus 5.5 <noreply@anthropic.com> Co-authored-by: Joshua Jewell <developer@joshuajewell.dev>
1 parent 578915c commit 592ddd7

5 files changed

Lines changed: 197 additions & 20 deletions

File tree

‎frontend/src/api/types.ts‎

Lines changed: 18 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -243,15 +243,32 @@ export interface ComparisonRequest extends AnalysisRequest {
243243
aggregate?: boolean
244244
}
245245

246+
export interface PermanovaTerm {
247+
term: string
248+
r2: number
249+
f_statistic: number | null
250+
/** Null when the term could not be tested; `untestable_reason` says why. */
251+
p_value: number | null
252+
/** Set when within-block permutation cannot change this term. */
253+
untestable_reason?: string | null
254+
}
255+
246256
export interface PermanovaResult {
247257
text: string
258+
/** The first term (group, when it varies); r2, f_statistic and p_value are its own. */
259+
term?: string
248260
r2: number | null
249261
f_statistic: number | null
250262
p_value: number | null
263+
untestable_reason?: string | null
264+
/** Permutations actually performed: fewer than 999 when the blocks admit fewer. */
265+
permutations?: number
266+
/** Smallest p-value attainable with that many permutations. */
267+
min_p_value?: number
251268
/** True when permutations were restricted within individuals. */
252269
blocked?: boolean
253270
/** One row per model term, tested sequentially in formula order. */
254-
terms?: { term: string; r2: number; f_statistic: number | null; p_value: number | null }[]
271+
terms?: PermanovaTerm[]
255272
/** PERMDISP on the same distances: do the groups differ in spread? */
256273
dispersion?: PermdispResult | { error: string }
257274
}

‎frontend/src/components/DiversityPanel.tsx‎

Lines changed: 8 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -99,19 +99,20 @@ export function DiversityPanel({ study, runs, option, aggregate }: {
9999
<div className="card" style={{ fontFamily: 'monospace', fontSize: '.82rem', whiteSpace: 'pre-wrap' }}>
100100
<div className="card-title">PERMANOVA Results</div>
101101
<pre style={{ margin: 0 }}>{permanova.text}</pre>
102-
{(permanova.terms?.length
103-
? permanova.terms
104-
: permanova.p_value != null
105-
? [{ term: 'model', r2: permanova.r2 ?? NaN, f_statistic: permanova.f_statistic, p_value: permanova.p_value }]
106-
: []
107-
).map(t => (
102+
{(permanova.terms ?? []).map(t => (
108103
<p key={t.term} style={{ marginTop: 8, fontFamily: 'inherit' }}>
109-
{t.term}: R&sup2; = {t.r2.toFixed(3)} · F = {t.f_statistic?.toFixed(2) ?? 'n/a'} · p = {t.p_value?.toFixed(4) ?? 'n/a'}
104+
{t.term}: R&sup2; = {t.r2.toFixed(3)} · F = {t.f_statistic?.toFixed(2) ?? 'n/a'} ·{' '}
105+
{t.untestable_reason ? `p not reported: ${t.untestable_reason}` : `p = ${t.p_value?.toFixed(4) ?? 'n/a'}`}
110106
</p>
111107
))}
112108
{(permanova.terms?.length ?? 0) > 1 && (
113109
<p style={{ marginTop: 4, fontFamily: 'inherit' }}>Terms are tested in order, each after the ones above it.</p>
114110
)}
111+
{permanova.permutations != null && permanova.min_p_value != null && (
112+
<p style={{ marginTop: 4, fontFamily: 'inherit' }}>
113+
{permanova.permutations} permutations; the smallest attainable p is {permanova.min_p_value.toFixed(4)}.
114+
</p>
115+
)}
115116
{permanova.blocked && (
116117
<p style={{ marginTop: 4, fontFamily: 'inherit' }}>Permutations restricted within individuals.</p>
117118
)}

‎src/analysis/analysis.jl‎

Lines changed: 67 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -1423,14 +1423,38 @@ function run_nmds(mat::Matrix{Float64}; seed::Integer=123, transform::String="no
14231423
end
14241424
end
14251425

1426+
"""
1427+
_varies_within_any_block(values, blocks) -> Bool
1428+
1429+
Whether `values` takes at least two distinct values among the samples of some
1430+
block in `blocks` (one label per sample).
1431+
"""
1432+
function _varies_within_any_block(values::AbstractVector, blocks::AbstractVector)
1433+
first_seen = Dict{eltype(blocks),eltype(values)}()
1434+
for (v, b) in zip(values, blocks)
1435+
isequal(get!(first_seen, b, v), v) || return true
1436+
end
1437+
false
1438+
end
1439+
14261440
"""
14271441
run_permanova(mat, metadata; seed=123, transform="none", blocks=nothing) -> Union{NamedTuple, Nothing}
14281442
1429-
PERMANOVA via vegan::adonis2 with 999 permutations.
1443+
PERMANOVA via vegan::adonis2 with up to 999 permutations.
14301444
`metadata` is a DataFrame with one row per sample and covariate columns.
1445+
Each covariate is tested as its own term (`by = "terms"`), sequentially in
1446+
column order, so each row is conditioned on the terms before it. The top-level
1447+
`term`, `r2`, `f_statistic` and `p_value` are those of the first term.
14311448
`blocks`, one label per sample (e.g. the individual a sample came from),
14321449
restricts permutations to within blocks. It is ignored when no label repeats,
1433-
since singleton blocks admit no permutation at all.
1450+
since singleton blocks admit no permutation at all. A term that is constant
1451+
within every block cannot be tested by within-block permutation (each sample
1452+
never leaves its block, so the permutations never exchange that term's labels
1453+
and the p-value does not test it), so its p-value is withheld and
1454+
`untestable_reason` says why.
1455+
`permutations` is the number of permutations actually performed, which is fewer
1456+
than 999 when the blocks admit fewer, and `min_p_value` the smallest p-value
1457+
that number can produce.
14341458
`transform` is applied before `vegdist`; see [`run_nmds`](@ref).
14351459
"""
14361460
function run_permanova(mat::Matrix{Float64}, metadata::DataFrame; seed::Integer=123,
@@ -1443,6 +1467,10 @@ function run_permanova(mat::Matrix{Float64}, metadata::DataFrame; seed::Integer=
14431467
formula_rhs = join(covariates, " + ")
14441468
mat = transform_counts(mat; method=transform)
14451469
blocked = !isnothing(blocks) && length(blocks) == size(mat, 1) && !allunique(blocks)
1470+
# Within-block permutation never exchanges the labels of a term that is
1471+
# constant inside every block, so its p-value does not test that term.
1472+
untestable = blocked ? [c for c in covariates if !_varies_within_any_block(metadata[!, c], blocks)] :
1473+
String[]
14461474

14471475
with_r_lock(; timeout=R_WAIT_SECONDS[]) do
14481476
meta_r = copy(metadata)
@@ -1452,6 +1480,7 @@ function run_permanova(mat::Matrix{Float64}, metadata::DataFrame; seed::Integer=
14521480
RCall.globalEnv[:seed] = Int(seed)
14531481
RCall.globalEnv[:perm_blocks] = blocked ? String.(blocks) : String[]
14541482
RCall.globalEnv[:disp_cols] = String.(covariates)
1483+
RCall.globalEnv[:perm_untestable] = String.(untestable)
14551484
RCall.reval("""
14561485
set.seed(seed, kind = "Mersenne-Twister", normal.kind = "Inversion", sample.kind = "Rejection")
14571486
dist_mat <- vegdist(mat, method = "bray")
@@ -1463,12 +1492,25 @@ function run_permanova(mat::Matrix{Float64}, metadata::DataFrame; seed::Integer=
14631492
999
14641493
}
14651494
perm_err <- NULL
1495+
perm_n <- NA_integer_
14661496
perm_res <- tryCatch(
1467-
adonis2(form, data = meta_r, permutations = perm_ctrl, parallel = 1),
1468-
error = function(e) { perm_err <<- conditionMessage(e); NULL }
1469-
)
1497+
# vegan 2.7 defaults to by = NULL, one test of the whole model.
1498+
adonis2(form, data = meta_r, permutations = perm_ctrl, by = "terms", parallel = 1),
1499+
error = function(e) { perm_err <<- conditionMessage(e); NULL })
14701500
if (!is.null(perm_res)) {
1471-
perm_text <- paste(capture.output(print(perm_res)), collapse = "\\n")
1501+
# permute enumerates every permutation instead of drawing 999 when
1502+
# the blocks admit fewer, which small blocks make likely; the
1503+
# control adonis2 used records how many it ran.
1504+
perm_n <- attr(perm_res, "control")[["nperm"]]
1505+
# Within-block permutation never exchanges such a term's labels
1506+
# between blocks, so its p-value does not test that term.
1507+
untestable_rows <- which(rownames(perm_res) %in% perm_untestable)
1508+
perm_res[untestable_rows, "Pr(>F)"] <- NA
1509+
perm_text <- paste(c(capture.output(print(perm_res)),
1510+
if (length(untestable_rows) > 0)
1511+
paste0("Pr(>F) withheld for ", paste(rownames(perm_res)[untestable_rows], collapse = ", "),
1512+
": constant within every permutation block, so within-block permutation does not test it.")),
1513+
collapse = "\\n")
14721514
perm_r2 <- perm_res\$R2[1]
14731515
perm_f <- perm_res\$F[1]
14741516
perm_p <- perm_res[["Pr(>F)"]][1]
@@ -1506,6 +1548,7 @@ function run_permanova(mat::Matrix{Float64}, metadata::DataFrame; seed::Integer=
15061548
r2 = RCall.rcopy(RCall.reval("perm_r2"))
15071549
f_stat = RCall.rcopy(RCall.reval("perm_f"))
15081550
p_val = RCall.rcopy(RCall.reval("perm_p"))
1551+
n_perm = RCall.rcopy(RCall.reval("perm_n"))
15091552
# adonis2 tests terms sequentially, so each row is conditioned on the terms above it.
15101553
term_names = String.(vcat(RCall.rcopy(RCall.reval("perm_terms"))))
15111554
term_r2 = vcat(RCall.rcopy(RCall.reval("perm_terms_r2")))
@@ -1520,20 +1563,32 @@ function run_permanova(mat::Matrix{Float64}, metadata::DataFrame; seed::Integer=
15201563
(; groups = [(; group = groups[k], mean = vcat(d[:mean])[k], median = vcat(d[:median])[k]) for k in eachindex(groups)],
15211564
f_statistic = d[:f], df = [Int(d[:df1]), Int(d[:df2])], p_value = d[:p], text = String(d[:text]))
15221565
end
1523-
RCall.reval("rm(mat, meta_r, formula_rhs, seed, perm_blocks, disp_cols, perm_ctrl, dist_mat, form, perm_res, perm_err, perm_text, perm_r2, perm_f, perm_p, perm_terms, perm_terms_r2, perm_terms_f, perm_terms_p, disp, disp_err); gc()")
1566+
RCall.reval("rm(mat, meta_r, formula_rhs, seed, perm_blocks, disp_cols, perm_untestable, perm_ctrl, dist_mat, form, perm_res, perm_err, perm_n, perm_text, perm_r2, perm_f, perm_p, perm_terms, perm_terms_r2, perm_terms_f, perm_terms_p, disp, disp_err); if (exists(\"untestable_rows\")) rm(untestable_rows); gc()")
15241567

15251568
ismissing(txt) && return nothing
15261569
# When R's adonis2 threw, txt is the error message and r2/f/p are missing.
15271570
# Return a named tuple with :message so the route can distinguish and surface it.
15281571
ismissing(r2) && return (; message=string(txt))
1572+
reason(t) = t in untestable ?
1573+
"'$t' is constant within every permutation block, and permutations are restricted " *
1574+
"within blocks, so its labels are never exchanged and the p-value would not test it; " *
1575+
"a between-block term needs whole blocks permuted" :
1576+
nothing
1577+
withheld(t, p) = (t in untestable || ismissing(p)) ? nothing : p
1578+
terms = [(; term=term_names[k],
1579+
r2=term_r2[k],
1580+
f_statistic=ismissing(term_f[k]) ? nothing : term_f[k],
1581+
p_value=withheld(term_names[k], term_p[k]),
1582+
untestable_reason=reason(term_names[k])) for k in eachindex(term_names)]
15291583
(; text=txt,
1584+
term=term_names[1],
15301585
r2=r2,
15311586
f_statistic=ismissing(f_stat) ? nothing : f_stat,
1532-
p_value=ismissing(p_val) ? nothing : p_val,
1533-
terms=[(; term=term_names[k],
1534-
r2=term_r2[k],
1535-
f_statistic=ismissing(term_f[k]) ? nothing : term_f[k],
1536-
p_value=ismissing(term_p[k]) ? nothing : term_p[k]) for k in eachindex(term_names)],
1587+
p_value=withheld(term_names[1], p_val),
1588+
untestable_reason=reason(term_names[1]),
1589+
terms,
1590+
permutations=Int(n_perm),
1591+
min_p_value=1 / (Int(n_perm) + 1),
15371592
blocked,
15381593
dispersion)
15391594
end

‎test/runtests.jl‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,7 @@ const UNIT_FILES = [
5050
"test_sample_reads.jl",
5151
"test_determinism.jl",
5252
"test_differential.jl",
53+
"test_permanova.jl",
5354
"test_write_gates.jl",
5455
]
5556

‎test/unit/test_permanova.jl‎

Lines changed: 103 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,103 @@
1+
# SPDX-License-Identifier: AGPL-3.0-only
2+
# SPDX-FileCopyrightText: 2026 Jonathan D.A. Jewell (hyperpolymath) <j.d.a.jewell@open.ac.uk>
3+
#
4+
# PERMANOVA: each covariate is tested as its own sequential term, checked against
5+
# a direct vegan::adonis2 call; the permutation count and p-value floor are
6+
# reported; and a term within-block permutation cannot test is refused.
7+
# Skips loudly when R/vegan is absent.
8+
9+
using MetaManifold
10+
using MetaManifold.Analysis
11+
using RCall
12+
using DataFrames
13+
using Random
14+
using Test
15+
16+
const _PERM_R = MetaManifold.Server.r_available()
17+
_PERM_R || @info "R/vegan unavailable - the PERMANOVA reference tests are SKIPPED"
18+
19+
"""
20+
_adonis2_reference(mat, meta, perm_expr) -> DataFrame
21+
22+
Run vegan::adonis2 directly on Bray-Curtis distances of `mat`, with
23+
`group + run` terms tested sequentially and permutations given by the R
24+
expression `perm_expr` (which may refer to `ref_meta`), seeded as
25+
`run_permanova` seeds it. Returns the term rows: term, r2, f, p.
26+
"""
27+
function _adonis2_reference(mat, meta, perm_expr)
28+
MetaManifold.RRuntime.with_r_lock() do
29+
RCall.globalEnv[:ref_mat] = mat
30+
RCall.globalEnv[:ref_meta] = meta
31+
RCall.reval("""
32+
set.seed(123, kind = "Mersenne-Twister", normal.kind = "Inversion", sample.kind = "Rejection")
33+
ref_d <- vegdist(ref_mat, method = "bray")
34+
ref_res <- adonis2(ref_d ~ group + run, data = ref_meta,
35+
permutations = $perm_expr, by = "terms", parallel = 1)
36+
ref_rows <- which(!(rownames(ref_res) %in% c("Residual", "Total")))
37+
ref_out <- data.frame(term = rownames(ref_res)[ref_rows], r2 = ref_res\$R2[ref_rows],
38+
f = ref_res\$F[ref_rows], p = ref_res[["Pr(>F)"]][ref_rows])
39+
""")
40+
out = rcopy(RCall.reval("ref_out"))
41+
RCall.reval("rm(ref_mat, ref_meta, ref_d, ref_res, ref_rows, ref_out)")
42+
out
43+
end
44+
end
45+
46+
@testset "PERMANOVA" begin
47+
@testset "_varies_within_any_block" begin
48+
blocks = ["i1", "i1", "i2", "i2"]
49+
@test !Analysis._varies_within_any_block(["a", "a", "b", "b"], blocks)
50+
@test Analysis._varies_within_any_block(["r1", "r2", "r1", "r1"], blocks)
51+
end
52+
53+
if _PERM_R
54+
rng = MersenneTwister(7)
55+
# Six individuals, each sampled in two runs; three individuals per group.
56+
mat = Float64.(rand(rng, 0:60, 12, 40))
57+
individuals = ["i$(cld(k, 2))" for k in 1:12]
58+
meta = DataFrame(sample = ["s$k" for k in 1:12],
59+
group = [cld(k, 2) <= 3 ? "A" : "B" for k in 1:12],
60+
run = repeat(["r1", "r2"], 6))
61+
62+
@testset "each term is tested on its own, as vegan by = \"terms\"" begin
63+
res = Analysis.run_permanova(mat, meta; seed=123)
64+
ref = _adonis2_reference(mat, meta, "999")
65+
@test [t.term for t in res.terms] == ["group", "run"] == ref.term
66+
@test [t.r2 for t in res.terms] ≈ ref.r2
67+
@test [t.f_statistic for t in res.terms] ≈ ref.f
68+
@test [t.p_value for t in res.terms] == ref.p
69+
# The headline result is the group term, not the whole model.
70+
@test res.term == "group"
71+
@test res.r2 ≈ ref.r2[1] && res.p_value == ref.p[1]
72+
@test res.permutations == 999
73+
@test res.min_p_value ≈ 1 / 1000
74+
@test all(t -> isnothing(t.untestable_reason), res.terms)
75+
end
76+
77+
@testset "within-individual blocks: permutation count, floor and refusal" begin
78+
res = Analysis.run_permanova(mat, meta; seed=123, blocks=individuals)
79+
ref = _adonis2_reference(mat, meta,
80+
"permute::how(blocks = factor(rep(paste0('i', 1:6), each = 2)), nperm = 999)")
81+
@test res.blocked
82+
# Six blocks of two admit 2^6 orderings, the observed one excluded.
83+
@test res.permutations == 63
84+
@test res.min_p_value ≈ 1 / 64
85+
# group is constant within every individual, so within-block
86+
# permutation never exchanges its labels. vegan still prints an
87+
# ordinary-looking probability for it (the control: that is the
88+
# number withheld here), but it does not test a group difference.
89+
@test 0 < ref.p[1] <= 1
90+
group = res.terms[1]
91+
@test group.term == "group" && isnothing(group.p_value)
92+
@test occursin("constant within every permutation block", group.untestable_reason)
93+
@test isnothing(res.p_value) && res.untestable_reason == group.untestable_reason
94+
@test group.r2 ≈ ref.r2[1]
95+
@test occursin("withheld for group", res.text)
96+
# run varies within individuals and is tested as vegan tests it.
97+
run = res.terms[2]
98+
@test run.term == "run" && isnothing(run.untestable_reason)
99+
@test run.p_value == ref.p[2]
100+
@test run.f_statistic ≈ ref.f[2]
101+
end
102+
end
103+
end

0 commit comments

Comments
 (0)