feat(concordance): Nelson-Ladiges "trit" currency and chance correction for QuartetConcordance - #276
Open
ms609 wants to merge 5 commits into
Open
feat(concordance): Nelson-Ladiges "trit" currency and chance correction for QuartetConcordance#276ms609 wants to merge 5 commits into
ms609 wants to merge 5 commits into
Conversation
…ance
Add `unit = c("quartet", "trit")` to QuartetConcordance(). "quartet"
(default) is unchanged; "trit" scores in Nelson & Ladiges (1992)
redundancy-corrected currency, where a split of sizes (k, t-k) carries
only (k-1)(t-k-1) independent quartets (4-cycles of K_{k,t-k}; the N&L
entailment is GF(2) cycle addition).
Scoring uses coverage ("option b"): per state-pair, concordant trits
A over the reported unit's own content (Wk for edges, Wc for chars),
pooled by shared information M = min(Wc, Wk). Only a split *displayed*
by the character scores 1; nested/crossing get strict partial credit.
Multistate is handled in the same currency, not deferred: a quartet is
decisive only within a state-pair, those K_{n_i,n_j} blocks are
edge-disjoint, and the entailment never crosses them, so trits add over
pairs (W = sum (n_i-1)(n_j-1), weight 4/(n_i n_j)). Verified by GF(2)
rank of the decisive/concordant sets on 2/3/4-state characters. Missing
/ ambiguous / inapplicable tokens drop per character by construction.
Validation (dev/benchmarks/frac-quart/): cell extraction cross-checked
against the C++ kernel; A <= min(Wc, Wk) everywhere; spec ladder
reproduced; package == an independent re-implementation on binary,
multistate and missing-data inputs across all weight x return combos.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add a random-expectation baseline to QuartetConcordance(unit = "trit"), mirroring ClusteringConcordance()'s normalize interface: normalize = FALSE (default) no correction; published measure unchanged normalize = TRUE exact expected concordance under a fixed-marginal null normalize = <int> Monte-Carlo tip-shuffle estimate of that expectation Null model: reassign each character's tokens across the leaves at random, holding state counts and split sizes fixed (the multivariate-hypergeometric confusion table; the same null as ClusteringConcordance's miRand and Consistency's ExpectedLength). This is the only coherent null for the trit currency -- it re-zeros each (split, char, state-pair) against its own combinatorial floor and extends unchanged to multistate + missing data. The flat 1/3 SCF baseline is a raw-quartet heuristic that does not survive the trit floors; it is documented as the rejected alternative. Implementation: only A (and, for multistate pairs where a pair's side-A count is random, wk and M = min(wc,wk)) vary under the null, so the exact estimator accumulates E[m], E[m*A/wk], E[m*A/wc] from the hypergeometric pmf (cached on (n_i,n_j,M,t)) and re-zeros the pooled edge/char ratio against the pooled expectation, with the same weight aggregation as the observed score. .Rezero guards the z->1 blow-up (-> NA) and returns values unclamped; .Rezero(1,z) = 1 so a displayed split still scores exactly 1 (ceiling invariant preserved), while conflicting characters go negative. Chance correction for the raw unit="quartet" currency is deferred (it needs the per-state-pair cells the C++ kernel does not expose, and a decision on its null) -- normalize errors there for now. Validated: exact vs MC tip-shuffle agree within MC error at the cell level (dev/benchmarks/frac-quart/chance_baseline.R) and end-to-end (test-Concordance.R). Design settled in dev/plans/frac-quart-weight-agreement.md section 7. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…istate
Two gaps in the unit="trit" chance correction, both in the multistate path:
- weight=FALSE re-zero averaged observed and baseline over DIFFERENT
(split, char) cells: observed drops a cell when denM==0 (na.rm), but the
baseline kept it whenever baseDenM>0 (a degenerate multistate pair can have
observed wk=0 yet E[m]>0). Mask both to NA on the union of NA cells so the
edge/char row-means cover the same population before .Rezero. weight=TRUE
(ratio-of-sums) is unaffected.
- The exact-vs-MC end-to-end test used only binary characters, so the
random-wk multistate regime -- the reason E[A] was elaborated into
E[m*A/wk] + E[m] -- had no oracle. Add a 3-state column and cover
edge/char x weight in {TRUE, FALSE}, asserting exact and MC agree on which
cells are NA (guards the cell-matching fix) and within MC error elsewhere.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…tet" Extend the chance-correction option to the raw quartet currency, using the same fixed-marginal (hypergeometric) null as the trit path. Per state-pair the raw concordant/decisive counts are polynomials in the 2x2 cells (no floors), so the exact expectation E[conc], E[dec] is a clean sum over the same trivariate-hypergeometric pmf (.ExpectedQuartet / .QuartetExpect); MC (.QuartetMC) reshuffles tokens and re-scores through the C++ kernel. The pooled conc/dec ratio is re-zeroed against E[conc]/E[dec] (weighted) or with the same NA-cell matching as the trit weight=FALSE path (unweighted). Removes the guard that errored on unit="quartet", normalize != FALSE. .Rezero(1, z) = 1, so a displayed split still scores exactly 1; conflicting characters go negative. normalize = FALSE is byte-identical to the published raw measure. Validated exact vs MC tip-shuffle across return x weight x binary/multistate (test-Concordance.R); all Concordance tests green. Design note in dev/plans/frac-quart-weight-agreement.md section 7 updated (raw unit DONE). This is the R reference for the forthcoming C++ port of the same expectation. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Move the exact hypergeometric expectation for QuartetConcordance(normalize = TRUE) from R to C++ (src/concordance_expect.cpp): `quartet_expect` returns E[conc], E[dec]; `trit_expect` returns E[m], E[m*A/wk], E[m*A/wc]. Both sum the per-state-pair trivariate-hypergeometric double sum, hoisting the log-choose terms out of the M/p/r loops (the dominant cost). The R exact helpers (.ExpectedTrit/.CharTritExpect/.ExpectedQuartet and the per-char .QuartetExpect loop) are replaced by thin wrappers / a single all-characters call in .TritConcordance; the Monte-Carlo (normalize = <int>) path stays in R. New Rcpp exports are registered in the manual R_CallMethodDef table (src/TreeSearch-init.c; the package uses R_useDynamicSymbols(FALSE)). Speedup at split-support scale (48 tips, 96 chars): exact normalize=TRUE quartet 2.88s -> 0.14s (~20x), trit 3.48s -> 0.45s (~8x); a full 1000-matrix run drops from ~100 min to ~10 min. Validated bit-for-bit (max |C++ - R| ~1e-14) against a self-contained R reference across binary/multistate/missing data (dev/benchmarks/frac-quart/cpp_expect_parity.R); the exact-vs-MC end-to-end tests in test-Concordance.R now exercise the C++ path and pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Adds a redundancy-corrected currency and an exact chance correction to
QuartetConcordance(), via two new arguments:unit = c("quartet", "trit")andnormalize.What this adds
unit = "trit"scores concordance in non-redundant quartet statements. The quartets a bipartition of sizes (k, t−k) resolves are the 4-cycles of K_{k,t−k}; Nelson–Ladiges entailment is GF(2) cycle addition, so only the cyclomatic number (k−1)(t−k−1) are independent. Multistate characters decompose into edge-disjoint state-pair blocks, so trits sum over pairs.normalize = TRUEsubtracts the exact expectation under a fixed-marginal hypergeometric null (both split sizes and the character's own state frequencies held fixed), computed as a closed double sum rather than simulated.normalize = <integer>gives a Monte-Carlo baseline instead.src/concordance_expect.cpp) for the expectation, mirroring the R reference to machine precision.Why chance correction matters here
Under the raw quartet currency the chance level is ~1/3 largely independently of split size, so correcting it is nearly a monotone transform — rank statistics are unchanged (Spearman between corrected and uncorrected = 1.0000000000 for the weighted score over 354,638 splits). Under the trit currency the
(x−1)+floors make the chance level drift with split size, so the correction reorders substantially (Spearman 0.39–0.56). Pooling uncorrected trit scores across split sizes therefore confounds signal with a size artefact; corrected, it does not.Validation
Rcpp::compileAttributes()produces no further changes.quartetandwQuartetreproduce previously cached values from an independent 999-alignment simulation study to 5.55e-16 across all 354,638 splits.Empirical behaviour
Scored against 999 simulated 48-taxon datasets (C-index vs graded similarity to the generating tree; ROC-AUC vs presence in it), on the parsimony-statistics subset (n = 36,761) the correction lifts unweighted quartet concordance from below chance (C 0.489, AUC 0.492) to 0.658/0.659, and the weighted score from 0.572 to 0.678. On the ML subset (n = 44,955) the weighted score goes 0.633/0.643 to 0.706/0.715.
Base
Targets
cpp-searchrather thanmaindeliberately: the character-recoding change inQuartetConcordance()introduced by 78b7414 (which alters whether{0,-}-style tokens count as grouping information) is a prerequisite for this work and is present oncpp-searchbut not onmain.🤖 Generated with Claude Code