From 871f125838a9e3d1da442a8be9ce36151a76d7e6 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 12:47:51 +0100 Subject: [PATCH 01/16] feat(concordance): add unit="trit" (Nelson-Ladiges) to QuartetConcordance 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 --- R/Concordance.R | 147 ++++++++++- dev/benchmarks/frac-quart/multistate_trit.R | 88 +++++++ dev/benchmarks/frac-quart/validate_trit.R | 229 ++++++++++++++++++ dev/benchmarks/frac-quart/verify_conc_disc.R | 73 ++++++ .../frac-quart/verify_direction_dep.R | 47 ++++ dev/benchmarks/frac-quart/verify_trit_ranks.R | 59 +++++ dev/plans/frac-quart-weight-agreement.md | 222 +++++++++++++++++ inst/REFERENCES.bib | 13 + man/SiteConcordance.Rd | 35 ++- tests/testthat/test-Concordance.R | 100 ++++++++ 10 files changed, 1007 insertions(+), 6 deletions(-) create mode 100644 dev/benchmarks/frac-quart/multistate_trit.R create mode 100644 dev/benchmarks/frac-quart/validate_trit.R create mode 100644 dev/benchmarks/frac-quart/verify_conc_disc.R create mode 100644 dev/benchmarks/frac-quart/verify_direction_dep.R create mode 100644 dev/benchmarks/frac-quart/verify_trit_ranks.R create mode 100644 dev/plans/frac-quart-weight-agreement.md diff --git a/R/Concordance.R b/R/Concordance.R index 2ba797245..7c755c9f9 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -765,8 +765,31 @@ MutualClusteringConcordance <- function(tree, dataset) { #' `QuartetConcordance(return = "char")` returns a numeric vector giving the #' concordance index calculated at each site, averaged across all splits. #' +#' With `unit = "trit"` (see below) the same vectors are returned, but scored in +#' redundancy-corrected currency: values are typically lower, and reach 1 only +#' where a split is *displayed* by the character rather than merely concordant +#' with many of its quartets. +#' #' @param weight Logical specifying whether to weight sites according to the #' number of quartets they are decisive for. +#' @param unit Character specifying the currency in which quartets are counted: +#' - `"quartet"` (default): each resolved quartet counts once, so a character +#' is credited with the full combinatorial volume of quartets it resolves; +#' - `"trit"`: quartets are counted as the *independent* information they carry +#' \insertCite{Nelson1992}{TreeSearch}. Because the quartets resolved by a +#' split of sizes \eqn{(k, t - k)} are logically redundant -- +#' \eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} -- only +#' \eqn{(k - 1)(t - k - 1)} of the \eqn{\binom{k}{2}\binom{t - k}{2}} +#' resolved quartets are independent. Concordance is then measured against +#' this reduced (redundancy-corrected) content, so that only a character +#' whose split is *identical* to the tree split scores full marks; nested +#' (compatible) characters receive genuine partial support, and crossing +#' (incompatible) characters score lower still. A multistate character is +#' scored in the same currency: trits are counted independently within each +#' pair of states (sizes \eqn{n_i}, \eqn{n_j} give a per-state-pair weight of +#' \eqn{4 / (n_i n_j)}) and summed across pairs, so multistate and binary +#' characters remain directly comparable. +#' @references \insertAllCited{} #' @importFrom ape keep.tip #' @importFrom cli cli_progress_bar cli_progress_update #' @importFrom utils combn @@ -776,7 +799,8 @@ QuartetConcordance <- function( tree, dataset = NULL, weight = TRUE, - return = "edge" + return = "edge", + unit = c("quartet", "trit") ) { if (is.null(dataset)) { warning("Cannot calculate concordance without `dataset`.") @@ -785,6 +809,7 @@ QuartetConcordance <- function( if (!inherits(dataset, "phyDat")) { stop("`dataset` must be a phyDat object.") } + unit <- match.arg(unit) tipLabels <- intersect(TipLabels(tree), names(dataset)) if (!length(tipLabels)) { warning("No overlap between tree labels and dataset.") @@ -816,14 +841,19 @@ QuartetConcordance <- function( dimnames = dimnames(characters) ) + options <- c("character", "site", "default") + return <- options[[pmatch(tolower(trimws(return)), options, + nomatch = length(options))]] + + if (unit == "trit") { + # Return: + return(.TritConcordance(logiSplits, charInt, weight, return, splits)) + } + raw_counts <- quartet_concordance(logiSplits, charInt) num <- raw_counts$concordant den <- raw_counts$decisive - options <- c("character", "site", "default") - return <- options[[pmatch(tolower(trimws(return)), options, - nomatch = length(options))]] - if (return == "default") { if (isTRUE(weight)) { @@ -869,6 +899,113 @@ QuartetConcordance <- function( } } +# Nelson-Ladiges fractional ("trit") currency for QuartetConcordance(). +# +# The quartets a bipartition of sizes (k, t - k) resolves are the 4-cycles of +# the complete bipartite graph K_{k, t-k}; the Nelson-Ladiges entailment +# (ab,cd) + (ab,ce) -> (ab,de) is GF(2) cycle addition, so only the cyclomatic +# number (k - 1)(t - k - 1) of them are independent ("trits"). +# +# A multistate character is the disjoint union of its state-pairs: a quartet is +# decisive only when two taxa share one state and two share another, so every +# decisive quartet lives in the K_{n_i, n_j} block between two states. Those +# blocks are edge-disjoint and the entailment never crosses them (it forces the +# shared "c, d, e" taxa into a single state), so trits add over state pairs: +# W_char = sum_{i 0, m * aij / wk, 0) + numChar[, ci] <- numChar[, ci] + ifelse(wc > 0, m * aij / wc, 0) + wcTot[ci] <- wcTot[ci] + wc[1] # wc is constant across splits + } + } + } + + informative <- wcTot > 0 + + if (return == "default") { + # edge: one value per split + if (isTRUE(weight)) { + denom <- rowSums(denM) + ret <- ifelse(denom == 0, NA_real_, rowSums(numEdge) / denom) + } else { + # Mean per-site quality over informative characters; uninformative + # characters carry no trits and are dropped, as in the quartet path. + sEdge <- ifelse(denM > 0, numEdge / denM, NA_real_) + ret <- if (any(informative)) { + rowMeans(sEdge[, informative, drop = FALSE], na.rm = TRUE) + } else { + rep(NA_real_, nSplit) + } + ret[is.nan(ret)] <- NA_real_ + } + setNames(ret, names(splits)) + } else { + # char: one value per character + if (isTRUE(weight)) { + denom <- colSums(denM) + ifelse(denom == 0, NA_real_, colSums(numChar) / denom) + } else { + sChar <- ifelse(denM > 0, numChar / denM, NA_real_) + vapply(seq_len(nChar), function(ci) { + if (informative[ci]) { + m <- mean(sChar[, ci], na.rm = TRUE) + if (is.nan(m)) NA_real_ else m + } else { + NA_real_ + } + }, double(1)) + } + } +} + .ExpectedMICache <- new.env(hash = TRUE, parent = emptyenv()) # @param a must be a vector of length <= 2 diff --git a/dev/benchmarks/frac-quart/multistate_trit.R b/dev/benchmarks/frac-quart/multistate_trit.R new file mode 100644 index 000000000..d11e38e1d --- /dev/null +++ b/dev/benchmarks/frac-quart/multistate_trit.R @@ -0,0 +1,88 @@ +# Verify: for a multistate character, GF(2) rank of its decisive quartets +# equals Sum_{i trits add across pairs. + +gf2rank <- function(M) { # rank over GF(2) + if (nrow(M) == 0 || ncol(M) == 0) return(0L) + M <- M %% 2L; r <- 0L; nc <- ncol(M) + for (col in seq_len(nc)) { + piv <- which(M[, col] == 1L) + piv <- piv[piv > r] + if (!length(piv)) next + p <- piv[1]; r <- r + 1L + M[c(r, p), ] <- M[c(p, r), ] + others <- which(M[, col] == 1L); others <- others[others != r] + if (length(others)) M[others, ] <- sweep(M[others, , drop = FALSE], 2, M[r, ], `+`) %% 2L + if (r == nrow(M)) break + } + r +} + +# edge id for unordered pair (a,b) among t taxa +eid <- function(a, b, t) { if (a > b) { tmp <- a; a <- b; b <- tmp }; (a - 1) * t + b } + +# All decisive quartets of `char` (state vector), each as its char-induced +# 4-cycle edge vector over C(t,2)... encoded on t*t (sparse-safe) columns. +quartet_rows <- function(char, split = NULL, concordant_only = FALSE) { + t <- length(char) + qs <- combn(t, 4) + rows <- list() + for (k in seq_len(ncol(qs))) { + q <- qs[, k]; st <- char[q] + tb <- table(st) + if (length(tb) == 2 && all(tb == 2)) { # "xx yy": decisive + states <- as.integer(names(tb)) + g1 <- q[st == states[1]]; g2 <- q[st == states[2]] # the two same-state pairs + if (concordant_only) { + s1 <- split[g1]; s2 <- split[g2] + # concordant iff each same-state pair lies wholly on one split side, + # and the two pairs are on opposite sides + ok <- length(unique(s1)) == 1 && length(unique(s2)) == 1 && s1[1] != s2[1] + if (!ok) next + } + v <- integer(t * t) + for (a in g1) for (b in g2) v[eid(a, b, t)] <- 1L # cross edges = 4-cycle + rows[[length(rows) + 1]] <- v + } + } + if (!length(rows)) return(matrix(0L, 0, t * t)) + do.call(rbind, rows) +} + +Aij <- function(p, q, r, s) max(p - 1, 0) * max(s - 1, 0) + max(q - 1, 0) * max(r - 1, 0) + +set.seed(1) +cat(sprintf("%-28s %6s %6s %6s\n", "case", "rankD", "SumW", "match")) +for (trial in 1:6) { + t <- sample(6:9, 1) + r <- sample(2:4, 1) + char <- sample(seq_len(r), t, replace = TRUE) + while (length(unique(char)) < 2 || any(tabulate(char) == 1 & tabulate(char) > 0 & FALSE)) char <- sample(seq_len(r), t, replace = TRUE) + ns <- as.integer(tabulate(char)); ns <- ns[ns > 0] + sumW <- sum(outer(seq_along(ns), seq_along(ns), Vectorize(function(i, j) + if (i < j) (ns[i] - 1) * (ns[j] - 1) else 0))) + rD <- gf2rank(quartet_rows(char)) + cat(sprintf("t=%d states=%-14s %6d %6d %6s\n", + t, paste(ns, collapse = ","), rD, sumW, rD == sumW)) +} + +cat("\n-- concordant vs split --\n") +cat(sprintf("%-34s %6s %6s %6s\n", "case", "rankC", "SumAij", "match")) +for (trial in 1:6) { + t <- sample(6:9, 1); r <- sample(2:4, 1) + char <- sample(seq_len(r), t, replace = TRUE) + while (length(unique(char)) < 2) char <- sample(seq_len(r), t, replace = TRUE) + split <- sample(c(TRUE, FALSE), t, replace = TRUE) + states <- sort(unique(char)) + sumA <- 0 + for (ii in seq_along(states)) for (jj in seq_along(states)) if (ii < jj) { + i <- states[ii]; j <- states[jj] + p <- sum(char == i & split); qc <- sum(char == i & !split) + rr <- sum(char == j & split); ss <- sum(char == j & !split) + sumA <- sumA + Aij(p, qc, rr, ss) + } + rC <- gf2rank(quartet_rows(char, split, concordant_only = TRUE)) + cat(sprintf("t=%d states=%-12s split=%-2d/%-2d %6d %6d %6s\n", + t, paste(as.integer(tabulate(char))[tabulate(char) > 0], collapse = ","), + sum(split), sum(!split), rC, sumA, rC == sumA)) +} diff --git a/dev/benchmarks/frac-quart/validate_trit.R b/dev/benchmarks/frac-quart/validate_trit.R new file mode 100644 index 000000000..b0f9142e9 --- /dev/null +++ b/dev/benchmarks/frac-quart/validate_trit.R @@ -0,0 +1,229 @@ +# Validation of the unit = "trit" path in QuartetConcordance(). +# +# Ordering follows the advisor's priority: +# 1. Cell extraction cross-checked against the verified C++ kernel (conc/dec). +# 2. A_ij <= min(Wc, Wk) asserted on every real (split, char, pair). +# 3. Spec ladder (t = 8) reproduced from the per-pair formula. +# 4. An INDEPENDENT re-implementation of the measure matches the package on +# binary (congreveLamsdell), multistate, and missing-data inputs. +# 5. `return` parsing quirk mirrored between quartet and trit paths. +# 6. congreveLamsdell quartet -> trit edge-score movement + regime census. +# +# Run: Rscript dev/benchmarks/frac-quart/validate_trit.R + +suppressMessages(pkgload::load_all(".", quiet = TRUE)) +library("TreeTools", quietly = TRUE) + +ok <- function(label, cond) { + cat(sprintf("[%s] %s\n", if (isTRUE(cond)) "PASS" else "FAIL", label)) + if (!isTRUE(cond)) stop("FAILED: ", label) +} + +# ---- Rebuild the (logiSplits, charInt) the function feeds to its kernels ---- +build_inputs <- function(tree, dataset) { + tipLabels <- intersect(TipLabels(tree), names(dataset)) + dataset <- dataset[tipLabels, drop = FALSE] + splits <- as.Splits(tree, dataset) + logiSplits <- vapply(seq_along(splits), function(i) as.logical(splits[[i]]), + logical(NTip(dataset))) + contrast <- attr(dataset, "contrast") + charLevels <- attr(dataset, "allLevels") + isInapp <- charLevels == "-" + isAmbig <- rowSums(contrast[, colnames(contrast) != "-", drop = FALSE]) > 1 + isGrouping <- !isAmbig & !isInapp + groupingCols <- apply(contrast[isGrouping, , drop = FALSE] > 0, 1, which) + levelToInt <- rep(NA_integer_, length(charLevels)) + levelToInt[isGrouping] <- as.integer(groupingCols) + characters <- PhyDatToMatrix(dataset) + charInt <- array(levelToInt[match(characters, charLevels)], + dim = dim(characters), dimnames = dimnames(characters)) + list(logiSplits = logiSplits, charInt = charInt, splits = splits) +} + +ch2 <- function(x) x * (x - 1) / 2 +posm <- function(z) { z[z < 0] <- 0; z } + +# Per (split, char) cells summed over state-pairs, returned as component lists +# so the same crosstab feeds both the kernel cross-check and the trit measure. +pair_cells <- function(logiSplits, charInt) { + nSplit <- ncol(logiSplits); nChar <- ncol(charInt) + out <- vector("list", nChar) + for (ci in seq_len(nChar)) { + col <- charInt[, ci]; scored <- !is.na(col) + states <- sort(unique(col[scored])) + prs <- list() + if (length(states) >= 2) { + for (a in seq_len(length(states) - 1L)) for (b in seq(a + 1L, length(states))) { + inI <- scored & col == states[a]; inJ <- scored & col == states[b] + prs[[length(prs) + 1L]] <- list( + aI = colSums(logiSplits & inI), bI = colSums(!logiSplits & inI), + aJ = colSums(logiSplits & inJ), bJ = colSums(!logiSplits & inJ)) + } + } + out[[ci]] <- prs + } + out +} + +# ---------- 1. Cell extraction vs the verified kernel ---------- +{ + data("congreveLamsdellMatrices", package = "TreeSearch") + dataset <- congreveLamsdellMatrices[[1]] + tree <- TreeSearch::referenceTree + inp <- build_inputs(tree, dataset) + cells <- pair_cells(inp$logiSplits, inp$charInt) + nSplit <- ncol(inp$logiSplits); nChar <- ncol(inp$charInt) + reconC <- reconD <- matrix(0, nSplit, nChar) + for (ci in seq_len(nChar)) for (p in cells[[ci]]) { + reconC[, ci] <- reconC[, ci] + ch2(p$aI) * ch2(p$bJ) + ch2(p$bI) * ch2(p$aJ) + reconD[, ci] <- reconD[, ci] + p$aI * p$bI * p$aJ * p$bJ + } + reconD <- reconD + reconC + kern <- TreeSearch:::quartet_concordance(inp$logiSplits, inp$charInt) + ok("cell -> concordant matches kernel", all(abs(reconC - kern$concordant) < 1e-9)) + ok("cell -> decisive matches kernel", all(abs(reconD - kern$decisive) < 1e-9)) + + # ---------- 2. A_ij <= min(Wc, Wk) on every pair ---------- + worst <- -Inf + for (ci in seq_len(nChar)) for (p in cells[[ci]]) { + aij <- posm(p$aI - 1) * posm(p$bJ - 1) + posm(p$bI - 1) * posm(p$aJ - 1) + nI <- p$aI + p$bI; nJ <- p$aJ + p$bJ; mA <- p$aI + p$aJ; tP <- nI + nJ + wc <- posm(nI - 1) * posm(nJ - 1); wk <- posm(mA - 1) * posm(tP - mA - 1) + worst <- max(worst, max(aij - pmin(wc, wk))) + } + ok("A_ij <= min(Wc, Wk) everywhere", worst <= 1e-9) +} + +# ---------- 3. Spec ladder (t = 8) from the per-pair formula ---------- +pairQ <- function(p, q, r, s) { # edge quality A / Wk for one 2x2 table + A <- posm(p - 1) * posm(s - 1) + posm(q - 1) * posm(r - 1) + m <- p + r; t <- p + q + r + s + Wk <- posm(m - 1) * posm(t - m - 1) + A / Wk +} +ladder <- rbind( + identical = c(4, 0, 0, 4), + nested = c(3, 0, 2, 3), + mildCrossing = c(3, 1, 1, 3), + maxCrossing = c(2, 2, 2, 2)) +Qlad <- apply(ladder, 1, function(x) pairQ(x[1], x[2], x[3], x[4])) +expected <- c(identical = 1, nested = 0.5, mildCrossing = 4 / 9, maxCrossing = 2 / 9) +ok("ladder reproduces spec table", all(abs(Qlad - expected) < 1e-9)) +ok("only identical scores 1", sum(abs(Qlad - 1) < 1e-9) == 1L) +ok("crossing < nested", Qlad["mildCrossing"] < Qlad["nested"] && + Qlad["maxCrossing"] < Qlad["mildCrossing"]) + +# ---------- 4. Independent re-implementation vs the package ---------- +# Deliberately different code shape: per-pair matrices, explicit pooling. +ref_trit <- function(tree, dataset, weight = TRUE, return = "edge") { + inp <- build_inputs(tree, dataset) + cells <- pair_cells(inp$logiSplits, inp$charInt) + nSplit <- ncol(inp$logiSplits); nChar <- ncol(inp$charInt) + numE <- numC <- den <- matrix(0, nSplit, nChar); wcTot <- numeric(nChar) + for (ci in seq_len(nChar)) for (p in cells[[ci]]) { + aij <- posm(p$aI - 1) * posm(p$bJ - 1) + posm(p$bI - 1) * posm(p$aJ - 1) + nI <- p$aI + p$bI; nJ <- p$aJ + p$bJ; mA <- p$aI + p$aJ; tP <- nI + nJ + wc <- posm(nI - 1) * posm(nJ - 1); wk <- posm(mA - 1) * posm(tP - mA - 1) + m <- pmin(wc, wk) + den[, ci] <- den[, ci] + m + numE[, ci] <- numE[, ci] + ifelse(wk > 0, m * aij / wk, 0) + numC[, ci] <- numC[, ci] + ifelse(wc > 0, m * aij / wc, 0) + wcTot[ci] <- wcTot[ci] + wc[1] + } + info <- wcTot > 0 + ret <- pmatch(tolower(trimws(return)), c("character", "site", "default"), + nomatch = 3L) + if (ret == 3L) { # edge + if (isTRUE(weight)) { + d <- rowSums(den); v <- ifelse(d == 0, NA_real_, rowSums(numE) / d) + } else { + sE <- ifelse(den > 0, numE / den, NA_real_) + v <- if (any(info)) rowMeans(sE[, info, drop = FALSE], na.rm = TRUE) + else rep(NA_real_, nSplit) + v[is.nan(v)] <- NA_real_ + } + setNames(v, names(inp$splits)) + } else { # char + if (isTRUE(weight)) { + d <- colSums(den); ifelse(d == 0, NA_real_, colSums(numC) / d) + } else { + sC <- ifelse(den > 0, numC / den, NA_real_) + vapply(seq_len(nChar), function(ci) if (info[ci]) { + mm <- mean(sC[, ci], na.rm = TRUE); if (is.nan(mm)) NA_real_ else mm + } else NA_real_, double(1)) + } + } +} + +same <- function(a, b) all((is.na(a) & is.na(b)) | (abs(a - b) < 1e-9), na.rm = FALSE) + +data("congreveLamsdellMatrices", package = "TreeSearch") +binDat <- congreveLamsdellMatrices[[1]] +refTree <- TreeSearch::referenceTree +for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) { + pkg <- QuartetConcordance(refTree, binDat, weight = w, return = r, unit = "trit") + rf <- ref_trit(refTree, binDat, weight = w, return = r) + ok(sprintf("binary pkg==ref weight=%s return=%s", w, r), same(pkg, rf)) +} + +# Multistate toy (3- and 4-state characters, complete) +msMat <- matrix(c(0, 0, 1, 1, 2, 2, 2, 0, # 3 states + 0, 0, 0, 1, 1, 2, 3, 3, # 4 states + 0, 0, 1, 1, 2, 2, 0, 1), 8, # 3 states + dimnames = list(paste0("t", 1:8), NULL)) +msDat <- MatrixToPhyDat(msMat) +msTree <- BalancedTree(8) +for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) { + pkg <- QuartetConcordance(msTree, msDat, weight = w, return = r, unit = "trit") + rf <- ref_trit(msTree, msDat, weight = w, return = r) + ok(sprintf("multi pkg==ref weight=%s return=%s", w, r), same(pkg, rf)) +} + +# Missing / ambiguous / inapplicable toy +naMat <- matrix(c(0, 0, 1, 1, "?", "?", 0, 1, + 0, 1, "-", 1, 0, "(01)", 1, 0, + 0, 0, 0, 1, 1, 2, "?", 2), 8, + dimnames = list(paste0("t", 1:8), NULL)) +naDat <- MatrixToPhyDat(naMat) +for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) { + pkg <- QuartetConcordance(msTree, naDat, weight = w, return = r, unit = "trit") + rf <- ref_trit(msTree, naDat, weight = w, return = r) + ok(sprintf("missing pkg==ref weight=%s return=%s", w, r), same(pkg, rf)) +} + +# ---------- 5. `return` parsing mirrors the quartet path ---------- +e1 <- QuartetConcordance(refTree, binDat, return = "edge", unit = "trit") +e2 <- QuartetConcordance(refTree, binDat, return = "default", unit = "trit") +c1 <- QuartetConcordance(refTree, binDat, return = "char", unit = "trit") +c2 <- QuartetConcordance(refTree, binDat, return = "character", unit = "trit") +c3 <- QuartetConcordance(refTree, binDat, return = "site", unit = "trit") +ok("return edge == default", same(e1, e2)) +ok("return char == character == site", same(c1, c2) && same(c1, c3)) +ok("edge length == n splits, char length == n chars", + length(e1) == length(inp <- as.Splits(refTree)) && + length(c1) == sum(attr(binDat, "weight"))) + +# ---------- 6. congreveLamsdell quartet -> trit movement + census ---------- +qEdge <- QuartetConcordance(refTree, binDat, unit = "quartet") +tEdge <- QuartetConcordance(refTree, binDat, unit = "trit") +cat("\nEdge concordance, quartet vs trit (congreveLamsdell[[1]]):\n") +print(round(rbind(quartet = qEdge, trit = tEdge, delta = tEdge - qEdge), 3)) +cat(sprintf("\nmean quartet = %.3f mean trit = %.3f (trit is stricter)\n", + mean(qEdge, na.rm = TRUE), mean(tEdge, na.rm = TRUE))) + +# Regime census per edge: identical / nested / crossing state-pairs +inpc <- build_inputs(refTree, binDat) +cells <- pair_cells(inpc$logiSplits, inpc$charInt) +regime <- matrix(0L, ncol(inpc$logiSplits), 3, + dimnames = list(names(inpc$splits), + c("identical", "nested", "crossing"))) +for (ci in seq_along(cells)) for (p in cells[[ci]]) { + empt <- (p$aI == 0) + (p$bI == 0) + (p$aJ == 0) + (p$bJ == 0) + regime[, "identical"] <- regime[, "identical"] + (empt == 2) + regime[, "nested"] <- regime[, "nested"] + (empt == 1) + regime[, "crossing"] <- regime[, "crossing"] + (empt == 0) +} +cat("\nState-pair regime census per edge:\n") +print(regime) + +cat("\nAll checks passed.\n") diff --git a/dev/benchmarks/frac-quart/verify_conc_disc.R b/dev/benchmarks/frac-quart/verify_conc_disc.R new file mode 100644 index 000000000..c60fd14c9 --- /dev/null +++ b/dev/benchmarks/frac-quart/verify_conc_disc.R @@ -0,0 +1,73 @@ +# Brute-force check of concordant/discordant/independent-quartet formulae +# for two binary bipartitions with 2x2 taxon table {p,q,r,s}. +choose2 <- function(x) x * (x - 1) / 2 + +# quartet resolution induced by a 0/1 labelling of 4 taxa: returns the +# unordered pair-of-pairs if 2-2, else NA +res <- function(lab) { + if (sum(lab) != 2) return(NA_character_) # not 2-2 + paste(sort(which(lab == lab[order(lab)][1]))[1:2], collapse = "") # placeholder +} +# cleaner: canonical partition key for a 2-2 split of positions 1:4 +key22 <- function(side) { + if (length(unique(side)) != 2 || sum(side == side[1]) != 2) return(NA_character_) + a <- sort(which(side == side[1])) + paste(a, collapse = "") # the pair on 'first' side identifies the 2|2 split up to complement +} + +set.seed(1) +ok <- TRUE +for (rep in 1:2000) { + t <- sample(6:11, 1) + charState <- sample(0:1, t, replace = TRUE) # I = state1 + splitSide <- sample(0:1, t, replace = TRUE) # A = side1 + p <- sum(charState == 1 & splitSide == 1) + q <- sum(charState == 1 & splitSide == 0) + r <- sum(charState == 0 & splitSide == 1) + s <- sum(charState == 0 & splitSide == 0) + + conc <- 0; disc <- 0; decisive <- 0 + combs <- combn(t, 4) + for (j in seq_len(ncol(combs))) { + idx <- combs[, j] + ck <- key22(charState[idx]); sk <- key22(splitSide[idx]) + if (!is.na(ck)) decisive <- decisive + 1 + if (!is.na(ck) && !is.na(sk)) { + # same 2|2 split iff the pair-key matches OR is complementary + same <- (ck == sk) || (ck == paste(setdiff(1:4, as.integer(strsplit(ck,"")[[1]])), collapse="")) + if (same) conc <- conc + 1 else disc <- disc + 1 + } + } + fConc <- choose2(p)*choose2(s) + choose2(q)*choose2(r) + fDisc <- p*q*r*s + n <- p + q + fDecisive <- choose2(n) * choose2(t - n) + if (conc != fConc || disc != fDisc || decisive != fDecisive) { + cat("MISMATCH pqrs=", p,q,r,s, + " conc", conc, fConc, " disc", disc, fDisc, + " dec", decisive, fDecisive, "\n"); ok <- FALSE; break + } + # fractional (independent) agreement with floors + pf <- function(x) pmax(x - 1, 0) + indConc <- pf(p)*pf(s) + pf(q)*pf(r) + indDec <- pf(n)*pf(t - n) # char's own independent quartets + if (indDec > 0) { + FQ <- indConc / indDec + if (FQ < -1e-9 || FQ > 1 + 1e-9) { cat("FQ out of range", FQ, p,q,r,s, "\n"); ok <- FALSE; break } + } + # self-agreement: char == split => q=r=0 => FQ==1 + if (q == 0 && r == 0 && indDec > 0 && abs(FQ - 1) > 1e-9) { + cat("self-agree != 1", FQ, p,q,r,s, "\n"); ok <- FALSE; break + } +} +cat(if (ok) "ALL CHECKS PASS\n" else "FAILED\n") + +# ind_decisive - ind_concordant identity (advisor: = pr + qs - 1 without floors) +set.seed(2); id_ok <- TRUE +for (i in 1:5000) { + p<-sample(0:6,1);q<-sample(0:6,1);r<-sample(0:6,1);s<-sample(0:6,1) + n<-p+q; tt<-p+q+r+s + lhs <- (n-1)*(tt-n-1) - ((p-1)*(s-1) + (q-1)*(r-1)) + if (lhs != p*r + q*s - 1) { cat("identity fail\n"); id_ok<-FALSE; break } +} +cat(if (id_ok) "IDENTITY (no-floor) ind_dec - ind_conc == pr+qs-1: CONFIRMED\n" else "IDENTITY FAIL\n") diff --git a/dev/benchmarks/frac-quart/verify_direction_dep.R b/dev/benchmarks/frac-quart/verify_direction_dep.R new file mode 100644 index 000000000..cf3fcad4d --- /dev/null +++ b/dev/benchmarks/frac-quart/verify_direction_dep.R @@ -0,0 +1,47 @@ +source("gfrank.R") # reuse gf2rank (prints table; ignore) +# Represent the SAME decisive quartet set two ways: +# - by the CHARACTER's pairing (4-cycle in K_{n,l}, I x O) -> Delta_char +# - by the SPLIT's pairing (4-cycle in K_{m,t-m}, A x B) -> Delta_split +# Concordant quartets pair the same either way; discordant differ. +delta_both <- function(p, q, r, s) { + n <- p + q; l <- r + s; m <- p + r; b <- q + s # split sizes + # taxon coords: state (I=1..n grouped p then q), side; build 4 cells + # cell membership: give each taxon (charState in {I,O}, side in {A,B}) + taxa <- rbind( + cbind(I=1, A=1, id=seq_len(p)), # p: I&A + cbind(I=1, A=0, id=seq_len(q)), # q: I&B + cbind(I=0, A=1, id=seq_len(r)), # r: O&A + cbind(I=0, A=0, id=seq_len(s))) # s: O&B + N <- nrow(taxa) + # index taxa within char-graph: I-vertices and O-vertices + Iidx <- which(taxa[,"I"]==1); Oidx <- which(taxa[,"I"]==0) + Aidx <- which(taxa[,"A"]==1); Bidx <- which(taxa[,"A"]==0) + vc <- match(seq_len(N), Iidx); vo <- match(seq_len(N), Oidx) # char graph coord + va <- match(seq_len(N), Aidx); vb <- match(seq_len(N), Bidx) # split graph coord + edgeC <- function(iA, iB) (iA-1)*length(Oidx) + iB # char edge id + edgeS <- function(a, b) (a-1)*length(Bidx) + b # split edge id + rowsC <- list(); rowsS <- list() + quads <- combn(N, 4) + for (j in seq_len(ncol(quads))) { + t4 <- quads[, j] + isI <- taxa[t4,"I"]==1; isA <- taxa[t4,"A"]==1 + if (sum(isI)!=2 || sum(isA)!=2) next # decisive = both resolve 2-2 + Is <- t4[isI]; Os <- t4[!isI]; As <- t4[isA]; Bs <- t4[!isA] + # char pairing: (Is)|(Os) + ec <- c(edgeC(vc[Is[1]],vo[Os[1]]), edgeC(vc[Is[1]],vo[Os[2]]), + edgeC(vc[Is[2]],vo[Os[1]]), edgeC(vc[Is[2]],vo[Os[2]])) + es <- c(edgeS(va[As[1]],vb[Bs[1]]), edgeS(va[As[1]],vb[Bs[2]]), + edgeS(va[As[2]],vb[Bs[1]]), edgeS(va[As[2]],vb[Bs[2]])) + vC <- integer(length(Iidx)*length(Oidx)); vC[ec] <- 1; rowsC[[length(rowsC)+1]] <- vC + vS <- integer(length(Aidx)*length(Bidx)); vS[es] <- 1; rowsS[[length(rowsS)+1]] <- vS + } + mk <- function(L, w) if (length(L)) do.call(rbind, L) else matrix(0,0,w) + c(Dchar = gf2rank(mk(rowsC, length(Iidx)*length(Oidx))), + Dsplit= gf2rank(mk(rowsS, length(Aidx)*length(Bidx))), + char_wt = (n-1)*(l-1), split_wt = (m-1)*(b-1)) +} +for (cfg in list(c(3,1,2,2), c(2,3,4,1), c(5,2,3,4), c(4,1,2,3))) { + d <- delta_both(cfg[1],cfg[2],cfg[3],cfg[4]) + cat(sprintf("p,q,r,s=%-10s Delta_char=%d (charWt=%d) Delta_split=%d (splitWt=%d)\n", + paste(cfg,collapse=","), d["Dchar"], d["char_wt"], d["Dsplit"], d["split_wt"])) +} diff --git a/dev/benchmarks/frac-quart/verify_trit_ranks.R b/dev/benchmarks/frac-quart/verify_trit_ranks.R new file mode 100644 index 000000000..f1830671c --- /dev/null +++ b/dev/benchmarks/frac-quart/verify_trit_ranks.R @@ -0,0 +1,59 @@ +# Independent-quartet counts as GF(2) ranks in the character's bipartite +# cycle space. A character quartet {i1,i2 | o1,o2} = a 4-cycle of K_{n,l} +# (n = I-taxa, l = O-taxa); entailment = cycle addition (Nelson-Ladiges). +gf2rank <- function(M) { # rank over GF(2) + if (!nrow(M) || !ncol(M)) return(0L) + M <- M %% 2; r <- 0L; ncol <- ncol(M); row <- 1L + for (col in seq_len(ncol)) { + piv <- which(M[row:nrow(M), col] == 1) + if (!length(piv)) next + piv <- piv[1] + row - 1L + tmp <- M[row, ]; M[row, ] <- M[piv, ]; M[piv, ] <- tmp + others <- which(M[, col] == 1); others <- others[others != row] + for (o in others) M[o, ] <- (M[o, ] + M[row, ]) %% 2 + row <- row + 1L; r <- r + 1L + if (row > nrow(M)) break + } + r +} + +# taxa: I-vertices 1..n (in A iff <=p), O-vertices 1..l (in A iff <=r) +quartet_vec <- function(i1, i2, o1, o2, l) { + v <- integer(0) # edge id = (i-1)*l + o + e <- c((i1-1)*l+o1, (i1-1)*l+o2, (i2-1)*l+o1, (i2-1)*l+o2) + e +} +counts <- function(p, q, r, s) { + n <- p + q; l <- r + s; nEdge <- n * l + inA_I <- function(i) i <= p + inA_O <- function(o) o <= r + rows_conc <- list(); rows_disc <- list(); rows_dec <- list() + for (i1 in 1:(n-1)) for (i2 in (i1+1):n) for (o1 in 1:(l-1)) for (o2 in (o1+1):l) { + nA <- inA_I(i1) + inA_I(i2) + inA_O(o1) + inA_O(o2) + if (nA != 2) next # split doesn't resolve -> not decisive + e <- c((i1-1)*l+o1, (i1-1)*l+o2, (i2-1)*l+o1, (i2-1)*l+o2) + v <- integer(nEdge); v[e] <- 1 + # concordant iff the two I-taxa are on the same split side + conc <- (inA_I(i1) == inA_I(i2)) + rows_dec[[length(rows_dec)+1]] <- v + if (conc) rows_conc[[length(rows_conc)+1]] <- v + else rows_disc[[length(rows_disc)+1]] <- v + } + mk <- function(L) if (length(L)) do.call(rbind, L) else matrix(0, 0, nEdge) + list(conc = gf2rank(mk(rows_conc)), + disc = gf2rank(mk(rows_disc)), + dec = gf2rank(mk(rows_dec)), + n_conc = length(rows_conc), n_disc = length(rows_disc)) +} +pf <- function(x) pmax(x - 1, 0) +cat(sprintf("%-14s %5s %5s %5s | %-22s %-10s\n", + "p,q,r,s", "conc", "disc", "dec", "A=(p1)(s1)+(q1)(r1)", "n-1 l-1")) +for (cfg in list(c(2,1,1,2), c(3,2,2,3), c(3,1,2,2), c(2,2,2,2), + c(4,1,1,4), c(3,3,1,1), c(2,3,4,1), c(5,2,3,4), + c(1,2,3,4), c(4,0,0,4), c(3,0,2,3), c(2,1,0,3))) { + p<-cfg[1];q<-cfg[2];r<-cfg[3];s<-cfg[4]; n<-p+q; l<-r+s + x <- counts(p,q,r,s) + A <- pf(p)*pf(s) + pf(q)*pf(r) + cat(sprintf("%-14s %5d %5d %5d | A=%-3d disc? dec vs %-4d| char=(%d)(%d)=%d\n", + paste(cfg,collapse=","), x$conc, x$disc, x$dec, A, A + x$disc, n-1, l-1, (n-1)*(l-1))) +} diff --git a/dev/plans/frac-quart-weight-agreement.md b/dev/plans/frac-quart-weight-agreement.md new file mode 100644 index 000000000..a7ff1b578 --- /dev/null +++ b/dev/plans/frac-quart-weight-agreement.md @@ -0,0 +1,222 @@ +# Fractional (Nelson–Ladiges) quartet weighting for `QuartetConcordance()` + +Branch: `frac-quart`. Goal: port Nelson & Ladiges (1992, *Syst. Biol.* 41:490–494) +fractional weighting of three-item statements into the **unrooted quartet** +concordance measure, so that logically redundant quartets are counted as the +information they actually carry (trits) rather than as raw combinatorial volume. + +This file is the settled design. Everything below was derived and, where noted, +verified by brute-force enumeration / GF(2) rank (scratch scripts, this session). + +--- + +## 1. Decisions locked in + +- **Unrooted quartets only.** No polarity / no "informative state" asymmetry. + N&L's rooted `2/n` symmetrises to the quartet weight below. +- **New arg `unit = c("quartet", "trit")`** on `QuartetConcordance()`. It + *replaces* any `fractional` flag. `"quartet"` = current behaviour (raw counts); + `"trit"` = N&L fractional (independent-quartet) currency. +- **Normalization = coverage by the *reported* unit** (option "b"): only a + character-split **identical** to the tree split scores full marks; nesting and + crossing get partial credit. See §4. +- **`weight` (TRUE/FALSE) is retained and orthogonal to `unit`.** It controls the + *pooling amount*, not the per-pair quality. See §5. +- **`return` (edge/char/tree) retained;** it selects which unit is "reported" and + therefore which weight normalises. See §4–5. +- This slots into the package's existing **quality × amount** decomposition + (`ConcordanceTable()` / `QACol()`): quality = per-pair coverage ratio, amount = + N&L fractional weight (trit content). + +--- + +## 2. Currency: trits, quartets = 4-cycles + +For `t` taxa, a binary character (no missing data) is a bipartition of sizes +`(k, t−k)`. The quartets it resolves are the **4-cycles of the complete bipartite +graph `K_{k, t−k}`** (two taxa each side = one 4-cycle). Nelson–Ladiges +entailment `(ab,cd) + (ab,ce) → (ab,de)` is exactly GF(2) cycle addition, so the +independent-quartet count is the cyclomatic number: + +``` +independent quartets of a split of sizes (k, t−k) = (k−1)(t−k−1) # "trits" +raw quartets it resolves = C(k,2)·C(t−k,2) +per-quartet fractional weight = 4 / (k(t−k)) # = indep/raw +``` + +`(k−1)(t−k−1)` is the "fractional weight" `W` of the whole character/split. It is +`0` for a constant (`k=0`) or autapomorphy (`k=1`) — so **uninformative characters +self-zero in trit currency**, no special-casing needed. + +Smallest N&L case (t=5, k=2): `4/(2·3)·3 = 2` trits from 3 raw quartets. ✓ + +**Multistate (solved, not deferred).** A quartet is decisive only when two taxa +share one state and two share another, so every decisive quartet lives in the +`K_{n_i, n_j}` edge-block between two states `i, j`. Those blocks are +edge-disjoint, and the N&L entailment never crosses them (it forces the shared +`c,d,e` taxa into one state), so **trits add over state pairs**: + +``` +W_char = Σ_{i 0` | genuine test; conflict lowers score | + +Nesting is a *passed test*, not homoplasy: no evidence against ⇒ evidence for +(cf. viviparity → Mammalia despite the platypus). It must get partial positive +credit — not exclusion. + +--- + +## 4. The measure (option b): coverage by the reported unit + +Per pair, quality = concordant trits over the **reported unit's** own weight: + +``` +edge (report split k): Q(k,c) = A(k,c) / W_k, W_k = (m−1)(t−m−1) +char (report char c): Q(k,c) = A(k,c) / W_c, W_c = (n−1)(t−n−1) +``` + +`A ≤ W_reported` always (⇒ `Q ∈ [0,1]`), and **per state-pair** +`Q_ij = 1 ⟺ A_ij = W_reported ⟺ identical sub-bipartition`. So for a **binary** +character, only an identical character-split scores full marks; nesting and +crossing are partial, with conflict strictly lowering the score. Verified ladder +(t=8, edge): + +| (p,q,r,s) | relationship | A | W_k | Q | +|---|---|---|---|---| +| 4,0,0,4 | identical | 9 | 9 | 1.00 | +| 3,0,2,3 | nested | 4 | 8 | 0.50 | +| 3,1,1,3 | mild crossing | 4 | 9 | 0.44 | +| 2,2,2,2 | max crossing | 2 | 9 | 0.22 | + +**Multistate: `Q = 1 ⟺ the split is displayed by the character`** (§5 pools the +per-pair `Q_ij` by `M = min(W_c, W_k)`). When every state block lies wholly on +one side of the split, each split-relevant pair scores `Q_ij = 1` and the one +split-orthogonal pair has `W_k^{ij} = 0 ⇒ M = 0` and drops out, so the pooled +`Q = 1` — even though `A_tot < W_c,tot` and the multistate character is *not* +"identical" to the (binary) split. This is the correct generalisation of +"identical" (for a binary character, displayed = identical), and the intended +behaviour: a reproductive character coded `{oviparous | viviparous | ovovivip.}` +that cleanly refines Mammalia fully supports the Mammalia split. A state block +that *straddles* the split makes its pairs crossing (`Q_ij < 1`), lowering the +pooled score. Verified: `t7,t8,t9`-refining char scores `1.0` on every displayed +split, `0.667` on splits that cut a block (see `test-Concordance.R`). + +Raw-currency (`unit="quartet"`) analogue: `Q = conc / (raw quartets of the +reported unit)` — same shape, `conc` and `C(m,2)C(t−m,2)` instead of `A`, `W_k`. + +--- + +## 5. Pooling (`weight` acts on the amount, not the quality) + +`W_k` is constant across characters for a fixed edge, so `weight` must act on the +**amount** (pooling weight), else it goes inert: + +``` +edge FQ(k) = Σ_c M(k,c)·Q(k,c) / Σ_c M(k,c) +char FQ(c) = Σ_k M(k,c)·Q(k,c) / Σ_k M(k,c) +``` + +- `M(k,c)` = shared information = **`min(W_c, W_k)`** (the hBest analogue used by + `ClusteringConcordance`). Preserves "only identical = full marks" (ceiling needs + `A = W_c = W_k`). +- `weight = FALSE` ⇒ `M ≡ 1` (per-pair-equal mean of `Q`). +- `unit = "quartet"` ⇒ `M` = raw shared-quartet count; `Q` uses raw conc. + +**Quality × amount map to the existing table:** `Q` is `QACol`'s *quality* axis; +`W_c` (trit content) is the *amount* axis. The trit port reuses +`ConcordanceTable()` rather than parallelling it. + +--- + +## 6. Build plan + +1. **R prototype first** (binary, no missing). Compute cells `{p,q,r,s}` per + (split, char) from the split logical vectors + character integer vectors; + derive `A`, `W_c`, `W_k`, `Q`, `M`; pool per §5. Keep the existing C++ kernel + for the raw path; the trit path can start in R since cells are a cheap crosstab. +2. **Validate**: (i) brute-force check `A` and the ladder on toy cells; + (ii) run on `congreveLamsdellMatrices` — report edge-score movement + `quartet → trit`, the identical/nested/crossing census per edge, and a + quality/amount table; (iii) confirm `unit="quartet"` reproduces current + `QuartetConcordance` bit-for-bit. +3. **Kernel** (later, perf only): the R per-state-pair path is already correct + and validated for binary + multistate + missing data, so a C++ port of `A` + is now a *performance* optimisation (large matrices), not a correctness + prerequisite. +4. **Tests**: extend `tests/testthat/test-QuartetConcordance*` — `unit` switch, + identical→1, nested partial, crossing < nested, uninformative→dropped, + `unit="quartet"` == legacy. + +## 7. Deferred / open + +- **Multistate**: ~~deferred~~ **DONE.** Trits add over state pairs + (`W_char = Σ_{i= 0 & qt <= 1, na.rm = TRUE)) + + # unit = "quartet" is the default and is byte-identical to omitting `unit`. + expect_identical(QuartetConcordance(tree, identChar), + QuartetConcordance(tree, identChar, unit = "quartet")) + # Coverage normalisation makes trit no laxer than quartet. + expect_true(mean(qt, na.rm = TRUE) <= + mean(QuartetConcordance(tree, identChar, unit = "quartet"), + na.rm = TRUE)) + + # Uninformative characters (constant / autapomorphy) carry no trits -> NA. + autap <- MatrixToPhyDat(matrix( + c(0, 0, 0, 0, 0, 0, 0, 1), 8, dimnames = list(paste0("t", 1:8), NULL))) + expect_equal(unname(QuartetConcordance(tree, autap, unit = "trit")), + rep(NA_real_, 5)) + + # `unit` is validated. + expect_error(QuartetConcordance(tree, identChar, unit = "trits"), + "should be one of") + + # `return` aliases mirror the quartet path. + expect_equal(QuartetConcordance(tree, identChar, return = "edge", unit = "trit"), + QuartetConcordance(tree, identChar, return = "default", unit = "trit")) + cA <- QuartetConcordance(tree, identChar, return = "char", unit = "trit") + expect_equal(cA, + QuartetConcordance(tree, identChar, return = "character", unit = "trit")) + expect_equal(cA, + QuartetConcordance(tree, identChar, return = "site", unit = "trit")) +}) + +test_that("QuartetConcordance() unit = 'trit' gives nested partial credit", { + # {t1..t5 | t6,t7,t8} split; state {t1,t2,t4} nested within side A but not a + # clade, giving cells (p,q,r,s) = (3,0,2,3) -> A/Wk = 4/8 = 0.5 exactly. + tree <- ape::read.tree(text = "(((((t1,t2),t3),t4),t5),(t6,(t7,t8)));") + char <- MatrixToPhyDat(matrix( + c(0, 0, 1, 0, 1, 1, 1, 1), 8, dimnames = list(paste0("t", 1:8), NULL))) + qt <- QuartetConcordance(tree, char, unit = "trit") + + # The split isolating {t6,t7,t8} scores exactly the nested value 0.5. + sp <- as.Splits(tree) + tips <- TipLabels(sp) + member <- vapply(seq_along(sp), function(i) { + side <- tips[as.logical(sp[[i]])] + setequal(side, c("t6", "t7", "t8")) || + setequal(setdiff(tips, side), c("t6", "t7", "t8")) + }, logical(1)) + expect_equal(unname(qt[member]), 0.5) +}) + +test_that("QuartetConcordance() unit = 'trit' supports multistate", { + # Multistate no longer errors: trits sum over state-pairs (same currency). + tree <- BalancedTree(8) + ms <- MatrixToPhyDat(matrix( + c(0, 0, 1, 1, 2, 2, 2, 0, # 3 states + 0, 0, 0, 1, 1, 2, 3, 3), 8, # 4 states + dimnames = list(paste0("t", 1:8), NULL))) + qt <- QuartetConcordance(tree, ms, unit = "trit") + expect_length(qt, 5L) + expect_true(all(qt >= 0 & qt <= 1, na.rm = TRUE)) + expect_false(anyNA(qt)) # both characters are informative on this tree +}) + +test_that("QuartetConcordance() unit = 'trit' scores 1 iff split displayed", { + # A multistate character need not be *identical* to a split to score 1: it + # scores full marks for every split its own tree displays (each state block + # wholly on one side; the split-orthogonal state-pair drops via M = 0). This + # generalises "only identical scores 1" from binary to multistate -- a cleanly + # refining reproductive character fully supports the clade it refines. + tree <- ape::read.tree(text = "((((t1,t2),t3),((t4,t5),t6)),((t7,t8),t9));") + char <- MatrixToPhyDat(matrix( + c(0, 0, 0, 1, 1, 1, 2, 2, 2), 9, + dimnames = list(paste0("t", 1:9), NULL))) + qt <- QuartetConcordance(tree, char, unit = "trit") + sp <- as.Splits(tree) + tips <- TipLabels(sp) + atSplit <- function(members) { + i <- which(vapply(seq_along(sp), function(j) { + side <- tips[as.logical(sp[[j]])] + setequal(side, members) || setequal(setdiff(tips, side), members) + }, logical(1))) + unname(qt[i]) + } + # Splits the character displays -> exactly 1 + expect_equal(atSplit(c("t1", "t2", "t3")), 1) + expect_equal(atSplit(c("t4", "t5", "t6")), 1) + expect_equal(atSplit(c("t1", "t2", "t3", "t4", "t5", "t6")), 1) + # A split cutting through a state block is only partially supported + expect_true(atSplit(c("t1", "t2")) < 1) +}) + test_that(".Rezero() works", { expect_equal(TreeSearch:::.Rezero(seq(0, 1, by = 0.1), 0.1), -1:9 / 9) }) From b5226a06fbaf67a2f752ae40d4897fd09f993fa4 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 13:19:18 +0100 Subject: [PATCH 02/16] feat(concordance): add normalize= chance correction to unit="trit" 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 = 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 --- R/Concordance.R | 307 ++++++++++++++++---- dev/benchmarks/frac-quart/chance_baseline.R | 86 ++++++ dev/plans/frac-quart-weight-agreement.md | 51 +++- man/ConcordanceTable.Rd | 30 +- man/SiteConcordance.Rd | 33 ++- tests/testthat/test-Concordance.R | 82 ++++++ 6 files changed, 516 insertions(+), 73 deletions(-) create mode 100644 dev/benchmarks/frac-quart/chance_baseline.R diff --git a/R/Concordance.R b/R/Concordance.R index 7c755c9f9..7db8c95cd 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -87,15 +87,30 @@ NULL #' #' Matching is case‑insensitive and partial. #' -#' @param normalize Controls how the *expected* mutual information (the zero -#' point of the scale) is determined. -#' - `FALSE`: no chance correction; MI is scaled only by its maximum. -#' - `TRUE`: subtract the analytical expected MI for random association. -#' - ``: subtract an empirical expected MI estimated from that -#' number of random trees. -#' -#' In all cases, 1 corresponds to the maximal attainable MI for the pair -#' (`hBest`), and 0 corresponds to the chosen expectation. +#' @param normalize Controls the zero point of the concordance scale by +#' subtracting the value expected under a chance (fixed-marginal) null, in which +#' each character's tokens are reassigned at random across the leaves while its +#' state frequencies and the split sizes are held fixed. +#' - `FALSE`: no chance correction; the measure is scaled only by its maximum, +#' so 1 marks a perfect match and 0 the measure's own (typically positive) +#' floor. +#' - `TRUE`: subtract the expected value under the null, so that 0 marks random +#' expectation and negative values fall below it. `ClusteringConcordance()` +#' uses an analytical approximation to the expected mutual information (fast +#' and generally accurate for large trees, ~200+ taxa, but neglecting +#' correlation between splits); `QuartetConcordance(unit = "trit")` uses the +#' exact hypergeometric expectation. +#' - a positive integer `n`: estimate that expectation by Monte Carlo. +#' `ClusteringConcordance()` fits each character to `n` random trees (and +#' returns Monte-Carlo standard errors, more accurate for small trees where +#' the analytical approximation is biased); `QuartetConcordance(unit = +#' "trit")` averages over `n` random reassignments of each character's +#' tokens. +#' +#' In all cases 1 corresponds to the maximum attainable value. For +#' `QuartetConcordance()`, chance correction is implemented for `unit = "trit"` +#' only, and the returned values are unclamped (they may fall below \eqn{-1}); +#' clamp to \eqn{[-1, 1]} before plotting with [QCol()] / [QACol()]. #' #' @returns #' `ClusteringConcordance(return = "all")` returns a 3D array where each @@ -800,7 +815,8 @@ QuartetConcordance <- function( dataset = NULL, weight = TRUE, return = "edge", - unit = c("quartet", "trit") + unit = c("quartet", "trit"), + normalize = FALSE ) { if (is.null(dataset)) { warning("Cannot calculate concordance without `dataset`.") @@ -810,6 +826,21 @@ QuartetConcordance <- function( stop("`dataset` must be a phyDat object.") } unit <- match.arg(unit) + if (!isFALSE(normalize)) { + if (!isTRUE(normalize)) { + if (!is.numeric(normalize) || length(normalize) != 1L || + is.na(normalize) || normalize < 1) { + stop("`normalize` must be FALSE, TRUE, or a positive integer.") + } + } + if (unit != "trit") { + # The chance baseline needs the per-state-pair cell counts, which the raw + # quartet path computes in the C++ kernel but does not expose; re-baselining + # the raw currency (and choosing its null; see the 1/3 note in the design + # spec) is deferred to the manuscript-figure review gate. + stop("`normalize` chance-correction is implemented for `unit = \"trit\"` only.") + } + } tipLabels <- intersect(TipLabels(tree), names(dataset)) if (!length(tipLabels)) { warning("No overlap between tree labels and dataset.") @@ -847,7 +878,8 @@ QuartetConcordance <- function( if (unit == "trit") { # Return: - return(.TritConcordance(logiSplits, charInt, weight, return, splits)) + return(.TritConcordance(logiSplits, charInt, weight, return, splits, + normalize)) } raw_counts <- quartet_concordance(logiSplits, charInt) @@ -924,7 +956,8 @@ QuartetConcordance <- function( # split so both `return`s pool by the same amount. Ambiguous / inapplicable / # absent tokens drop out per character (treated as "?", as in the quartet path), # giving each character its own effective taxon count. -.TritConcordance <- function(logiSplits, charInt, weight, return, splits) { +.TritConcordance <- function(logiSplits, charInt, weight, return, splits, + normalize = FALSE) { nSplit <- ncol(logiSplits) nChar <- ncol(charInt) pos <- function(z) { @@ -932,39 +965,32 @@ QuartetConcordance <- function( z } - # Contributions to the pools, summed over each character's state-pairs. + # Observed contributions to the pools, summed over each character's + # state-pairs; the chance baseline (if requested) uses the same accumulators + # for the expected pools, so the re-zero divides like against like. numEdge <- numChar <- denM <- matrix(0, nSplit, nChar) wcTot <- numeric(nChar) # character trit content (split-free) + doNorm <- !isFALSE(normalize) + if (doNorm) { + baseNumEdge <- baseNumChar <- baseDenM <- matrix(0, nSplit, nChar) + } for (ci in seq_len(nChar)) { col <- charInt[, ci] - scored <- !is.na(col) - states <- sort(unique(col[scored])) - if (length(states) < 2L) { - next # constant / autapomorphic / empty - } - for (a in seq_len(length(states) - 1L)) { - for (b in seq(a + 1L, length(states))) { - inI <- scored & col == states[a] - inJ <- scored & col == states[b] - aI <- colSums(logiSplits & inI) # state i, side A - bI <- colSums(!logiSplits & inI) # state i, side B - aJ <- colSums(logiSplits & inJ) # state j, side A - bJ <- colSums(!logiSplits & inJ) # state j, side B - nI <- aI + bI # n_i (constant across splits) - nJ <- aJ + bJ # n_j - mA <- aI + aJ # taxa of this pair on side A - tP <- nI + nJ # taxa scored in this pair - # Concordant trits; the (x - 1)_+ floors stop self-agreement exceeding 1. - aij <- pos(aI - 1) * pos(bJ - 1) + pos(bI - 1) * pos(aJ - 1) - wc <- pos(nI - 1) * pos(nJ - 1) # pair's character content - wk <- pos(mA - 1) * pos(tP - mA - 1) # pair's split content - m <- pmin(wc, wk) # shared information - denM[, ci] <- denM[, ci] + m - numEdge[, ci] <- numEdge[, ci] + ifelse(wk > 0, m * aij / wk, 0) - numChar[, ci] <- numChar[, ci] + ifelse(wc > 0, m * aij / wc, 0) - wcTot[ci] <- wcTot[ci] + wc[1] # wc is constant across splits + obs <- .CharTritContrib(col, logiSplits, pos) + numEdge[, ci] <- obs[["numEdge"]] + numChar[, ci] <- obs[["numChar"]] + denM[, ci] <- obs[["denM"]] + wcTot[ci] <- obs[["wcTot"]] + if (doNorm && obs[["wcTot"]] > 0) { + base <- if (isTRUE(normalize)) { + .CharTritExpect(col, logiSplits, pos) # exact hypergeometric + } else { + .CharTritMC(col, logiSplits, pos, normalize) # Monte-Carlo tip-shuffle } + baseNumEdge[, ci] <- base[["numEdge"]] + baseNumChar[, ci] <- base[["numChar"]] + baseDenM[, ci] <- base[["denM"]] } } @@ -975,37 +1001,212 @@ QuartetConcordance <- function( if (isTRUE(weight)) { denom <- rowSums(denM) ret <- ifelse(denom == 0, NA_real_, rowSums(numEdge) / denom) + if (doNorm) { + bDen <- rowSums(baseDenM) + base <- ifelse(bDen == 0, NA_real_, rowSums(baseNumEdge) / bDen) + ret <- .RezeroGuarded(ret, base) + } } else { # Mean per-site quality over informative characters; uninformative # characters carry no trits and are dropped, as in the quartet path. sEdge <- ifelse(denM > 0, numEdge / denM, NA_real_) - ret <- if (any(informative)) { - rowMeans(sEdge[, informative, drop = FALSE], na.rm = TRUE) - } else { - rep(NA_real_, nSplit) + ret <- .RowMeanInformative(sEdge, informative, nSplit) + if (doNorm) { + sBase <- ifelse(baseDenM > 0, baseNumEdge / baseDenM, NA_real_) + base <- .RowMeanInformative(sBase, informative, nSplit) + ret <- .RezeroGuarded(ret, base) } - ret[is.nan(ret)] <- NA_real_ } setNames(ret, names(splits)) } else { # char: one value per character if (isTRUE(weight)) { denom <- colSums(denM) - ifelse(denom == 0, NA_real_, colSums(numChar) / denom) + ret <- ifelse(denom == 0, NA_real_, colSums(numChar) / denom) + if (doNorm) { + bDen <- colSums(baseDenM) + base <- ifelse(bDen == 0, NA_real_, colSums(baseNumChar) / bDen) + ret <- .RezeroGuarded(ret, base) + } + ret } else { sChar <- ifelse(denM > 0, numChar / denM, NA_real_) - vapply(seq_len(nChar), function(ci) { - if (informative[ci]) { - m <- mean(sChar[, ci], na.rm = TRUE) - if (is.nan(m)) NA_real_ else m - } else { - NA_real_ - } - }, double(1)) + ret <- .ColMeanInformative(sChar, informative, nChar) + if (doNorm) { + sBase <- ifelse(baseDenM > 0, baseNumChar / baseDenM, NA_real_) + base <- .ColMeanInformative(sBase, informative, nChar) + ret <- .RezeroGuarded(ret, base) + } + ret } } } +# Observed per-character trit contributions, summed over the character's +# state-pairs, returned as per-split vectors. Shared by the observed pass and +# the Monte-Carlo baseline so the two use byte-identical arithmetic. +.CharTritContrib <- function(col, logiSplits, pos) { + nSplit <- ncol(logiSplits) + numEdge <- numChar <- denM <- numeric(nSplit) + wcTot <- 0 + scored <- !is.na(col) + states <- sort(unique(col[scored])) + nStates <- length(states) + if (nStates >= 2L) { + for (a in seq_len(nStates - 1L)) { + for (b in seq(a + 1L, nStates)) { + inI <- scored & col == states[a] + inJ <- scored & col == states[b] + aI <- colSums(logiSplits & inI) # state i, side A + bI <- colSums(!logiSplits & inI) # state i, side B + aJ <- colSums(logiSplits & inJ) # state j, side A + bJ <- colSums(!logiSplits & inJ) # state j, side B + nI <- aI + bI # n_i (constant across splits) + nJ <- aJ + bJ # n_j + mA <- aI + aJ # taxa of this pair on side A + tP <- nI + nJ # taxa scored in this pair + # Concordant trits; the (x - 1)_+ floors stop self-agreement exceeding 1. + aij <- pos(aI - 1) * pos(bJ - 1) + pos(bI - 1) * pos(aJ - 1) + wc <- pos(nI - 1) * pos(nJ - 1) # pair's character content + wk <- pos(mA - 1) * pos(tP - mA - 1) # pair's split content + m <- pmin(wc, wk) # shared information + denM <- denM + m + numEdge <- numEdge + ifelse(wk > 0, m * aij / wk, 0) + numChar <- numChar + ifelse(wc > 0, m * aij / wc, 0) + wcTot <- wcTot + wc[1] # wc is constant across splits + } + } + } + list(numEdge = numEdge, numChar = numChar, denM = denM, wcTot = wcTot) +} + +# Exact expected trit contributions under the fixed-marginal (hypergeometric) +# null: each split's scored side-A size `M` and each pair's state counts +# `n_i, n_j` are held fixed while the character's tokens are reassigned at +# random across the `t` scored leaves. Only `A`, and (for multistate pairs) +# `wk` and `m`, vary; so we accumulate E[m], E[m * A / wk] and E[m * A / wc] +# from the exact pmf, matching `.CharTritContrib`'s pool definitions. +.CharTritExpect <- function(col, logiSplits, pos) { + nSplit <- ncol(logiSplits) + numEdge <- numChar <- denM <- numeric(nSplit) + scored <- !is.na(col) + tc <- sum(scored) + mSideA <- colSums(logiSplits & scored) # M per split + uM <- unique(mSideA) + idx <- match(mSideA, uM) + states <- sort(unique(col[scored])) + nStates <- length(states) + if (nStates >= 2L) { + cnt <- tabulate(match(col[scored], states), nStates) # state counts + for (a in seq_len(nStates - 1L)) { + for (b in seq(a + 1L, nStates)) { + et <- vapply(uM, function(M) .ExpectedTrit(cnt[a], cnt[b], M, tc), + double(3)) + denM <- denM + et[1, idx] + numEdge <- numEdge + et[2, idx] + numChar <- numChar + et[3, idx] + } + } + } + list(numEdge = numEdge, numChar = numChar, denM = denM) +} + +.ExpectedTritCache <- new.env(hash = TRUE, parent = emptyenv()) + +# Exact E[m], E[m * A / wk], E[m * A / wc] for one state-pair, over the +# trivariate hypergeometric placement of `nI` state-i and `nJ` state-j leaves +# (and `t - nI - nJ` other leaves) among the `M` leaves on side A of a split. +.ExpectedTrit <- function(nI, nJ, M, t) { + key <- paste(nI, nJ, M, t, sep = ",") + cached <- .ExpectedTritCache[[key]] + if (!is.null(cached)) { + return(cached) + } + pos1 <- function(z) if (z < 0) 0 else z + nOther <- t - nI - nJ + tP <- nI + nJ # taxa in this pair + wc <- pos1(nI - 1) * pos1(nJ - 1) # pair's character content (fixed) + logDen <- lchoose(t, M) + acc <- c(0, 0, 0) # E[m], E[m*A/wk], E[m*A/wc] + for (p in 0:min(nI, M)) { + for (r in 0:min(nJ, M - p)) { + oA <- M - p - r # other-state leaves on side A + if (oA < 0 || oA > nOther) { + next + } + prob <- exp(lchoose(nI, p) + lchoose(nJ, r) + lchoose(nOther, oA) - logDen) + if (prob <= 0) { + next + } + q <- nI - p + s <- nJ - r + mA <- p + r # this pair's taxa on side A (random) + A <- pos1(p - 1) * pos1(s - 1) + pos1(q - 1) * pos1(r - 1) + wk <- pos1(mA - 1) * pos1(tP - mA - 1) # pair's split content (random) + m <- min(wc, wk) + acc[1] <- acc[1] + prob * m + if (wk > 0) acc[2] <- acc[2] + prob * m * A / wk + if (wc > 0) acc[3] <- acc[3] + prob * m * A / wc + } + } + .ExpectedTritCache[[key]] <- acc + acc +} + +# Monte-Carlo baseline: average `.CharTritContrib` over `nRelabel` random +# reassignments of the character's tokens across its scored leaves (the split, +# the scored set and the state counts are held fixed). Averaging the pools +# (rather than per-split ratios) mirrors the exact ratio-of-expectations path. +.CharTritMC <- function(col, logiSplits, pos, nRelabel) { + nSplit <- ncol(logiSplits) + scored <- !is.na(col) + tokens <- col[scored] + accEdge <- accChar <- accDen <- numeric(nSplit) + for (i in seq_len(nRelabel)) { + shuffled <- col + shuffled[scored] <- sample(tokens) + cc <- .CharTritContrib(shuffled, logiSplits, pos) + accEdge <- accEdge + cc[["numEdge"]] + accChar <- accChar + cc[["numChar"]] + accDen <- accDen + cc[["denM"]] + } + list(numEdge = accEdge / nRelabel, + numChar = accChar / nRelabel, + denM = accDen / nRelabel) +} + +# Mean of per-split quality across informative characters, matching the +# quartet path's treatment of uninformative characters (dropped). +.RowMeanInformative <- function(mat, informative, nSplit) { + ret <- if (any(informative)) { + rowMeans(mat[, informative, drop = FALSE], na.rm = TRUE) + } else { + rep(NA_real_, nSplit) + } + ret[is.nan(ret)] <- NA_real_ + ret +} + +.ColMeanInformative <- function(mat, informative, nChar) { + vapply(seq_len(nChar), function(ci) { + if (informative[ci]) { + m <- mean(mat[, ci], na.rm = TRUE) + if (is.nan(m)) NA_real_ else m + } else { + NA_real_ + } + }, double(1)) +} + +# Re-zero to a chance baseline, guarding the degenerate `zero -> 1` case +# (returns NA rather than dividing by ~0). Values are returned UNCLAMPED, as +# in `ClusteringConcordance`; clamping to [-1, 1] is deferred to plotting. +.RezeroGuarded <- function(value, zero) { + denom <- 1 - zero + ifelse(is.na(value) | is.na(zero) | denom < sqrt(.Machine[["double.eps"]]), + NA_real_, (value - zero) / denom) +} + .ExpectedMICache <- new.env(hash = TRUE, parent = emptyenv()) # @param a must be a vector of length <= 2 diff --git a/dev/benchmarks/frac-quart/chance_baseline.R b/dev/benchmarks/frac-quart/chance_baseline.R new file mode 100644 index 000000000..93b9d9690 --- /dev/null +++ b/dev/benchmarks/frac-quart/chance_baseline.R @@ -0,0 +1,86 @@ +# Validate the exact hypergeometric chance baseline used by +# QuartetConcordance(unit = "trit", normalize = TRUE) against a Monte-Carlo +# tip-shuffle oracle, at the level of a single (split, character, state-pair). +# +# Under the fixed-marginal null the character's tokens are reassigned across the +# t scored leaves at random, holding the pair's state counts (n_i, n_j) and the +# split's scored side-A size M fixed. The cell vector (p,q,r,s) is then +# trivariate-hypergeometric. We need the EXPECTED per-pair pool contributions +# m = min(w_c, w_k) (shared-information pooling weight) +# m*A/w_k (edge-numerator contribution) +# m*A/w_c (char-numerator contribution) +# because w_k = (mA-1)+(tP-mA-1)+ -- and hence m -- are THEMSELVES random for +# multistate pairs (mA = p + r varies). Binary characters are the special case +# mA == M, tP == t, so w_k is fixed and only A varies (the clean +# ClusteringConcordance case). This mirrors `.ExpectedTrit()` in R/Concordance.R. + +pos <- function(z) { z[z < 0] <- 0; z } + +# Observed per-pair triple for one realised cell vector. +pairTriple <- function(p, q, r, s) { + nI <- p + q; nJ <- r + s + mA <- p + r; tP <- nI + nJ + A <- pos(p - 1) * pos(s - 1) + pos(q - 1) * pos(r - 1) + wc <- pos(nI - 1) * pos(nJ - 1) + wk <- pos(mA - 1) * pos(tP - mA - 1) + m <- min(wc, wk) + c(m = m, + mAwk = if (wk > 0) m * A / wk else 0, + mAwc = if (wc > 0) m * A / wc else 0) +} + +# Exact expectation over the trivariate hypergeometric (the R kernel's method). +expectedTriple <- function(nI, nJ, M, t) { + nOther <- t - nI - nJ + stopifnot(nOther >= 0, M >= 0, M <= t) + logC <- function(n, k) if (k < 0 || k > n) -Inf else lchoose(n, k) + denom <- lchoose(t, M) + acc <- c(m = 0, mAwk = 0, mAwc = 0) + for (p in 0:min(nI, M)) { + for (r in 0:min(nJ, M - p)) { + oA <- M - p - r + if (oA < 0 || oA > nOther) next + prob <- exp(logC(nI, p) + logC(nJ, r) + logC(nOther, oA) - denom) + if (prob <= 0) next + acc <- acc + prob * pairTriple(p, nI - p, r, nJ - r) + } + } + acc +} + +# Monte-Carlo oracle: shuffle state labels across t leaves, M fixed on side A. +mcTriple <- function(nI, nJ, M, t, nShuf = 2e5) { + sideA <- c(rep(TRUE, M), rep(FALSE, t - M)) + acc <- c(m = 0, mAwk = 0, mAwc = 0) + for (i in seq_len(nShuf)) { + lab <- sample(rep(1:3, c(nI, nJ, t - nI - nJ))) # 1 = i, 2 = j, 3 = other + p <- sum(lab == 1 & sideA); q <- sum(lab == 1 & !sideA) + r <- sum(lab == 2 & sideA); s <- sum(lab == 2 & !sideA) + acc <- acc + pairTriple(p, q, r, s) + } + acc / nShuf +} + +set.seed(1) +cases <- list( + c(nI = 4, nJ = 4, M = 4, t = 8), # binary-like (nOther = 0): w_k fixed + c(nI = 3, nJ = 5, M = 4, t = 8), + c(nI = 3, nJ = 3, M = 4, t = 9), # multistate (nOther > 0): w_k random + c(nI = 4, nJ = 2, M = 5, t = 10), + c(nI = 5, nJ = 4, M = 6, t = 12), + c(nI = 2, nJ = 2, M = 3, t = 7), + c(nI = 6, nJ = 3, M = 5, t = 15) +) + +worst <- 0 +for (cs in cases) { + ex <- expectedTriple(cs["nI"], cs["nJ"], cs["M"], cs["t"]) + mc <- mcTriple(cs["nI"], cs["nJ"], cs["M"], cs["t"]) + worst <- max(worst, max(abs(ex - mc))) + cat(sprintf("nI=%d nJ=%d M=%2d t=%2d exact=(%.4f,%.4f,%.4f) mc=(%.4f,%.4f,%.4f)\n", + cs["nI"], cs["nJ"], cs["M"], cs["t"], + ex[1], ex[2], ex[3], mc[1], mc[2], mc[3])) +} +cat(sprintf("\nWorst exact-vs-MC discrepancy: %.5f\n", worst)) +stopifnot(worst < 0.05) +cat("PASS: exact hypergeometric baseline matches the MC tip-shuffle oracle.\n") diff --git a/dev/plans/frac-quart-weight-agreement.md b/dev/plans/frac-quart-weight-agreement.md index a7ff1b578..6840a9142 100644 --- a/dev/plans/frac-quart-weight-agreement.md +++ b/dev/plans/frac-quart-weight-agreement.md @@ -212,9 +212,54 @@ char FQ(c) = Σ_k M(k,c)·Q(k,c) / Σ_k M(k,c) toy with `?`/`-`/ambiguous tokens (`validate_trit.R`), but the validation corpus (congreveLamsdell) has no missing data, so the semantics are proven self-consistent (`Q∈[0,1]`, pkg==independent ref) rather than benchmarked. -- **Random-expectation baseline**: a `.Rezero`/`miRand` analogue (crossing floor - ≈ 0.22 in §4 is not zero) — decide whether trit concordance needs chance - correction. Deferred. +- **Random-expectation baseline**: ~~deferred~~ **DONE** (`unit = "trit"`). + New arg `normalize = FALSE / TRUE / ` on `QuartetConcordance()`, mirroring + `ClusteringConcordance()`. The crossing floor (≈ 0.22 in §4) is not zero, so + without correction a maximally conflicting character still scores positive and + the floor is split-size dependent; `normalize` re-zeros against the + concordance *expected by chance*, exactly as `ClusteringConcordance()` re-zeros + MI against `miRand`. + - **Null model (chosen): fixed-marginal randomization** — reassign each + character's tokens across the leaves at random, holding its state counts and + the split sizes fixed (the multivariate-hypergeometric confusion table; the + same null as `ClusteringConcordance`'s `miRand` and `Consistency`'s + `ExpectedLength` tip-shuffle). 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 and missing data. + - **Rejected alternative: the flat Minh-style `1/3`** (three quartet + resolutions ⇒ expected concordance `1/3`). It is a *raw-quartet* heuristic + that ignores marginals and does **not** extend to floored trit counts (the + trit floor is `0.22`, not `1/3`), so it cannot re-zero the trit measure + coherently. Its only remaining use would be re-baselining the *raw* `quartet` + unit — see the raw-unit note below. + - **Two estimators** (mirroring `ClusteringConcordance`): `normalize = TRUE` + computes the **exact** expected pools from the hypergeometric pmf (a small + double sum, cached on `(n_i, n_j, M, t)`); `normalize = ` estimates them + by **Monte-Carlo** tip-shuffle. Validated to agree within MC error + (`dev/benchmarks/frac-quart/chance_baseline.R` at the cell level; + `test-Concordance.R` end-to-end). + - **What varies under the null.** Only `A` (and, for *multistate* pairs where a + pair's side-A count `mA` is itself random, `wk` and `M = min(w_c, w_k)`) vary; + for *binary* characters `mA = M`, `tP = t` so `wk` is fixed and only `A` + varies (the clean `ClusteringConcordance` case). The exact estimator therefore + accumulates `E[m]`, `E[m·A/w_k]`, `E[m·A/w_c]` and re-zeros the *pooled* + edge/char ratio against the pooled expectation, with the **same** `weight` + aggregation as the observed score. + - **Guarantees / guards.** `.Rezero(1, z) = 1`, so a **displayed** split still + scores exactly `1` (the ceiling invariant of §4 survives chance correction — + only the floor is lifted); conflicting characters go **negative** ("below + random", as in `ClusteringConcordance`). Values are returned **unclamped** + (may fall below `−1`); clamp to `[−1, 1]` only at plot time (`QCol`/`QACol`), + never inside the measure. The `z → 1` blow-up is guarded (→ `NA`). + `normalize` defaults to `FALSE`, so the published measure is unchanged unless + the user opts in. + - **Raw `quartet` unit: deferred.** `normalize` currently errors for + `unit = "quartet"`. Re-baselining the raw currency needs the per-state-pair + cells (computed inside the C++ kernel but not exposed) and a decision on its + null (fixed-marginal vs the flat `1/3`); both are folded into the C++ port + (task 2) and the manuscript-figure review (task 3), which will compare + trit-with vs trit-without correction and decide whether the raw unit needs it + too. - **`M = min(W_c,W_k)`** (settled): the pooling amount is the *shared* information, so it must be symmetric between character and split (both `return`s pool by the same `M`); `M = W_k` would make the char return weight a diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index 9d6a0831b..cc617bca8 100644 --- a/man/ConcordanceTable.Rd +++ b/man/ConcordanceTable.Rd @@ -36,17 +36,31 @@ at edges whose descendants are both contain more than \code{largeClade} leaves.} \item{ylab}{Character giving a label for the y axis.} -\item{normalize}{Controls how the \emph{expected} mutual information (the zero -point of the scale) is determined. +\item{normalize}{Controls the zero point of the concordance scale by +subtracting the value expected under a chance (fixed-marginal) null, in which +each character's tokens are reassigned at random across the leaves while its +state frequencies and the split sizes are held fixed. \itemize{ -\item \code{FALSE}: no chance correction; MI is scaled only by its maximum. -\item \code{TRUE}: subtract the analytical expected MI for random association. -\item \verb{}: subtract an empirical expected MI estimated from that -number of random trees. +\item \code{FALSE}: no chance correction; the measure is scaled only by its maximum, +so 1 marks a perfect match and 0 the measure's own (typically positive) +floor. +\item \code{TRUE}: subtract the expected value under the null, so that 0 marks random +expectation and negative values fall below it. \code{ClusteringConcordance()} +uses an analytical approximation to the expected mutual information (fast +and generally accurate for large trees, ~200+ taxa, but neglecting +correlation between splits); \code{QuartetConcordance(unit = "trit")} uses the +exact hypergeometric expectation. +\item a positive integer \code{n}: estimate that expectation by Monte Carlo. +\code{ClusteringConcordance()} fits each character to \code{n} random trees (and +returns Monte-Carlo standard errors, more accurate for small trees where +the analytical approximation is biased); \code{QuartetConcordance(unit = "trit")} averages over \code{n} random reassignments of each character's +tokens. } -In all cases, 1 corresponds to the maximal attainable MI for the pair -(\code{hBest}), and 0 corresponds to the chosen expectation.} +In all cases 1 corresponds to the maximum attainable value. For +\code{QuartetConcordance()}, chance correction is implemented for \code{unit = "trit"} +only, and the returned values are unclamped (they may fall below \eqn{-1}); +clamp to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{plot}{Logical specifying whether to draw the plot.} diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index f61cddce8..3dd1a5604 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -18,7 +18,8 @@ QuartetConcordance( dataset = NULL, weight = TRUE, return = "edge", - unit = c("quartet", "trit") + unit = c("quartet", "trit"), + normalize = FALSE ) PhylogeneticConcordance(tree, dataset) @@ -43,17 +44,31 @@ split–character pair. Matching is case‑insensitive and partial.} -\item{normalize}{Controls how the \emph{expected} mutual information (the zero -point of the scale) is determined. +\item{normalize}{Controls the zero point of the concordance scale by +subtracting the value expected under a chance (fixed-marginal) null, in which +each character's tokens are reassigned at random across the leaves while its +state frequencies and the split sizes are held fixed. \itemize{ -\item \code{FALSE}: no chance correction; MI is scaled only by its maximum. -\item \code{TRUE}: subtract the analytical expected MI for random association. -\item \verb{}: subtract an empirical expected MI estimated from that -number of random trees. +\item \code{FALSE}: no chance correction; the measure is scaled only by its maximum, +so 1 marks a perfect match and 0 the measure's own (typically positive) +floor. +\item \code{TRUE}: subtract the expected value under the null, so that 0 marks random +expectation and negative values fall below it. \code{ClusteringConcordance()} +uses an analytical approximation to the expected mutual information (fast +and generally accurate for large trees, ~200+ taxa, but neglecting +correlation between splits); \code{QuartetConcordance(unit = "trit")} uses the +exact hypergeometric expectation. +\item a positive integer \code{n}: estimate that expectation by Monte Carlo. +\code{ClusteringConcordance()} fits each character to \code{n} random trees (and +returns Monte-Carlo standard errors, more accurate for small trees where +the analytical approximation is biased); \code{QuartetConcordance(unit = "trit")} averages over \code{n} random reassignments of each character's +tokens. } -In all cases, 1 corresponds to the maximal attainable MI for the pair -(\code{hBest}), and 0 corresponds to the chosen expectation.} +In all cases 1 corresponds to the maximum attainable value. For +\code{QuartetConcordance()}, chance correction is implemented for \code{unit = "trit"} +only, and the returned values are unclamped (they may fall below \eqn{-1}); +clamp to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{weight}{Logical specifying whether to weight sites according to the number of quartets they are decisive for.} diff --git a/tests/testthat/test-Concordance.R b/tests/testthat/test-Concordance.R index 88e32f625..83c7adcd2 100644 --- a/tests/testthat/test-Concordance.R +++ b/tests/testthat/test-Concordance.R @@ -233,6 +233,88 @@ test_that("QuartetConcordance() unit = 'trit' scores 1 iff split displayed", { expect_true(atSplit(c("t1", "t2")) < 1) }) +test_that("QuartetConcordance() unit = 'trit' normalize contract", { + tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") + ident <- MatrixToPhyDat(matrix( + c(0, 0, 0, 0, 0, 0, 1, 1, 1), 9, + dimnames = list(paste0("t", 1:9), NULL))) + + # normalize = FALSE is the default and leaves the measure untouched. + expect_identical( + QuartetConcordance(tree, ident, unit = "trit", normalize = FALSE), + QuartetConcordance(tree, ident, unit = "trit")) + + # Invalid normalize is rejected. + expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = 0), + "positive integer") + expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = -3), + "positive integer") + expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = "x"), + "positive integer") + + # Chance correction is implemented for the trit currency only, for now. + expect_error( + QuartetConcordance(tree, ident, unit = "quartet", normalize = TRUE), + "trit") + + # An identical (displayed) split still scores exactly 1 after re-zeroing: + # .Rezero(1, z) == 1, so chance correction lifts the floor without moving the + # ceiling. + sp <- as.Splits(tree) + tips <- TipLabels(sp) + col789 <- which(vapply(seq_along(sp), function(j) { + side <- tips[as.logical(sp[[j]])] + setequal(side, c("t7", "t8", "t9")) || + setequal(setdiff(tips, side), c("t7", "t8", "t9")) + }, logical(1))) + corrected <- QuartetConcordance(tree, ident, unit = "trit", normalize = TRUE) + expect_equal(unname(corrected[col789]), 1) +}) + +test_that("QuartetConcordance() unit = 'trit' exact baseline matches MC", { + tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") + dat <- MatrixToPhyDat(matrix( + c(0, 0, 0, 0, 0, 0, 1, 1, 1, # displays {t7,t8,t9} + 0, 0, 0, 0, 0, 0, 0, 1, 1, # nested: {t8,t9} + 0, 0, 1, 1, 0, 0, 1, 1, 0), 9, # crossing + dimnames = list(paste0("t", 1:9), NULL))) + + # The exact hypergeometric baseline and the Monte-Carlo tip-shuffle baseline + # target the same null, so they agree to within MC error (edge and char). + exact <- QuartetConcordance(tree, dat, unit = "trit", normalize = TRUE) + set.seed(1) + mc <- QuartetConcordance(tree, dat, unit = "trit", normalize = 6000) + expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.03) + + exactC <- QuartetConcordance(tree, dat, unit = "trit", return = "char", + normalize = TRUE) + set.seed(2) + mcC <- QuartetConcordance(tree, dat, unit = "trit", return = "char", + normalize = 6000) + expect_lt(max(abs(exactC - mcC), na.rm = TRUE), 0.03) +}) + +test_that("QuartetConcordance() unit = 'trit' correction: conflict below random", { + tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") + support <- MatrixToPhyDat(matrix( # agrees with the tree: displays {t7,t8,t9} + c(0, 0, 0, 0, 0, 0, 1, 1, 1), 9, + dimnames = list(paste0("t", 1:9), NULL))) + conflict <- MatrixToPhyDat(matrix( # crosses the two halves of the tree + c(0, 0, 1, 1, 0, 0, 1, 1, 0), 9, + dimnames = list(paste0("t", 1:9), NULL))) + + sChar <- QuartetConcordance(tree, support, return = "char", unit = "trit", + normalize = TRUE) + cNone <- QuartetConcordance(tree, conflict, return = "char", unit = "trit") + cChar <- QuartetConcordance(tree, conflict, return = "char", unit = "trit", + normalize = TRUE) + # A character that agrees with the tree scores above random expectation; + # a crossing character is pulled below it, and below its uncorrected score. + expect_gt(sChar, 0) + expect_lt(cChar, cNone) + expect_lt(cChar, 0) +}) + test_that(".Rezero() works", { expect_equal(TreeSearch:::.Rezero(seq(0, 1, by = 0.1), 0.1), -1:9 / 9) }) From 1318e8f4313200b3c0604d10744833de51c0dc18 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 13:24:55 +0100 Subject: [PATCH 03/16] fix(concordance): match cells in weight=FALSE trit re-zero; test multistate 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 --- R/Concordance.R | 25 +++++++++++++++++++------ tests/testthat/test-Concordance.R | 31 ++++++++++++++++++------------- 2 files changed, 37 insertions(+), 19 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index 7db8c95cd..7e3542ba6 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -1010,11 +1010,19 @@ QuartetConcordance <- function( # Mean per-site quality over informative characters; uninformative # characters carry no trits and are dropped, as in the quartet path. sEdge <- ifelse(denM > 0, numEdge / denM, NA_real_) - ret <- .RowMeanInformative(sEdge, informative, nSplit) if (doNorm) { sBase <- ifelse(baseDenM > 0, baseNumEdge / baseDenM, NA_real_) - base <- .RowMeanInformative(sBase, informative, nSplit) - ret <- .RezeroGuarded(ret, base) + # Average observed and baseline over the SAME (split, char) cells: a + # degenerate multistate pair can be dropped from one but not the other + # (observed wk = 0 yet E[m] > 0), which would otherwise re-zero mismatched + # populations. + naMask <- is.na(sEdge) | is.na(sBase) + sEdge[naMask] <- NA_real_ + sBase[naMask] <- NA_real_ + ret <- .RezeroGuarded(.RowMeanInformative(sEdge, informative, nSplit), + .RowMeanInformative(sBase, informative, nSplit)) + } else { + ret <- .RowMeanInformative(sEdge, informative, nSplit) } } setNames(ret, names(splits)) @@ -1031,11 +1039,16 @@ QuartetConcordance <- function( ret } else { sChar <- ifelse(denM > 0, numChar / denM, NA_real_) - ret <- .ColMeanInformative(sChar, informative, nChar) if (doNorm) { sBase <- ifelse(baseDenM > 0, baseNumChar / baseDenM, NA_real_) - base <- .ColMeanInformative(sBase, informative, nChar) - ret <- .RezeroGuarded(ret, base) + # Match cells, as in the edge return above. + naMask <- is.na(sChar) | is.na(sBase) + sChar[naMask] <- NA_real_ + sBase[naMask] <- NA_real_ + ret <- .RezeroGuarded(.ColMeanInformative(sChar, informative, nChar), + .ColMeanInformative(sBase, informative, nChar)) + } else { + ret <- .ColMeanInformative(sChar, informative, nChar) } ret } diff --git a/tests/testthat/test-Concordance.R b/tests/testthat/test-Concordance.R index 83c7adcd2..c1e100f16 100644 --- a/tests/testthat/test-Concordance.R +++ b/tests/testthat/test-Concordance.R @@ -276,22 +276,27 @@ test_that("QuartetConcordance() unit = 'trit' exact baseline matches MC", { dat <- MatrixToPhyDat(matrix( c(0, 0, 0, 0, 0, 0, 1, 1, 1, # displays {t7,t8,t9} 0, 0, 0, 0, 0, 0, 0, 1, 1, # nested: {t8,t9} - 0, 0, 1, 1, 0, 0, 1, 1, 0), 9, # crossing + 0, 0, 1, 1, 0, 0, 1, 1, 0, # crossing (binary) + 0, 0, 0, 1, 1, 1, 2, 2, 2), 9, # 3-state: exercises the random-wk path dimnames = list(paste0("t", 1:9), NULL))) # The exact hypergeometric baseline and the Monte-Carlo tip-shuffle baseline - # target the same null, so they agree to within MC error (edge and char). - exact <- QuartetConcordance(tree, dat, unit = "trit", normalize = TRUE) - set.seed(1) - mc <- QuartetConcordance(tree, dat, unit = "trit", normalize = 6000) - expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.03) - - exactC <- QuartetConcordance(tree, dat, unit = "trit", return = "char", - normalize = TRUE) - set.seed(2) - mcC <- QuartetConcordance(tree, dat, unit = "trit", return = "char", - normalize = 6000) - expect_lt(max(abs(exactC - mcC), na.rm = TRUE), 0.03) + # target the same null, so they agree to within MC error -- for edge and char, + # under both weightings, and crucially through the MULTISTATE pair path where a + # pair's wk and M = min(w_c, w_k) are themselves random under the null. + for (ret in c("edge", "char")) { + for (w in c(TRUE, FALSE)) { + exact <- QuartetConcordance(tree, dat, unit = "trit", return = ret, + weight = w, normalize = TRUE) + set.seed(1) + mc <- QuartetConcordance(tree, dat, unit = "trit", return = ret, + weight = w, normalize = 8000) + # Same cells resolve to NA under exact and MC (guards the cell-matching in + # the weight = FALSE path), and the finite values agree within MC error. + expect_identical(is.na(exact), is.na(mc)) + expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04) + } + } }) test_that("QuartetConcordance() unit = 'trit' correction: conflict below random", { From 2ac6a6c3c4679f3bda063f2405c9496c44267dc8 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 14:51:08 +0100 Subject: [PATCH 04/16] feat(concordance): enable normalize= chance correction for unit="quartet" 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 --- R/Concordance.R | 198 +++++++++++++++++++---- dev/plans/frac-quart-weight-agreement.md | 19 ++- man/ConcordanceTable.Rd | 14 +- man/SiteConcordance.Rd | 14 +- tests/testthat/test-Concordance.R | 41 ++++- 5 files changed, 230 insertions(+), 56 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index 7e3542ba6..ea2f353b7 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -98,19 +98,18 @@ NULL #' expectation and negative values fall below it. `ClusteringConcordance()` #' uses an analytical approximation to the expected mutual information (fast #' and generally accurate for large trees, ~200+ taxa, but neglecting -#' correlation between splits); `QuartetConcordance(unit = "trit")` uses the -#' exact hypergeometric expectation. +#' correlation between splits); `QuartetConcordance()` uses the exact +#' hypergeometric expectation (for either `unit`). #' - a positive integer `n`: estimate that expectation by Monte Carlo. #' `ClusteringConcordance()` fits each character to `n` random trees (and #' returns Monte-Carlo standard errors, more accurate for small trees where -#' the analytical approximation is biased); `QuartetConcordance(unit = -#' "trit")` averages over `n` random reassignments of each character's -#' tokens. +#' the analytical approximation is biased); `QuartetConcordance()` averages +#' over `n` random reassignments of each character's tokens. #' #' In all cases 1 corresponds to the maximum attainable value. For -#' `QuartetConcordance()`, chance correction is implemented for `unit = "trit"` -#' only, and the returned values are unclamped (they may fall below \eqn{-1}); -#' clamp to \eqn{[-1, 1]} before plotting with [QCol()] / [QACol()]. +#' `QuartetConcordance()` (either `unit`), chance-corrected values are returned +#' unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before +#' plotting with [QCol()] / [QACol()]. #' #' @returns #' `ClusteringConcordance(return = "all")` returns a 3D array where each @@ -833,13 +832,6 @@ QuartetConcordance <- function( stop("`normalize` must be FALSE, TRUE, or a positive integer.") } } - if (unit != "trit") { - # The chance baseline needs the per-state-pair cell counts, which the raw - # quartet path computes in the C++ kernel but does not expose; re-baselining - # the raw currency (and choosing its null; see the 1/3 note in the design - # spec) is deferred to the manuscript-figure review gate. - stop("`normalize` chance-correction is implemented for `unit = \"trit\"` only.") - } } tipLabels <- intersect(TipLabels(tree), names(dataset)) if (!length(tipLabels)) { @@ -887,6 +879,23 @@ QuartetConcordance <- function( num <- raw_counts$concordant den <- raw_counts$decisive + # Chance correction (option): re-zero the observed concordant/decisive ratio + # against the ratio expected under the same fixed-marginal null used for the + # trit currency. Only `conc` and `dec` vary under the null (the split sizes + # and state counts are fixed), so we need E[conc] and E[dec] per (split, char), + # computed exactly from the hypergeometric pmf (no floors here, unlike trits) + # or by Monte-Carlo tip-shuffle, then re-zero the pooled ratio. + doNorm <- !isFALSE(normalize) + if (doNorm) { + base <- if (isTRUE(normalize)) { + .QuartetExpect(charInt, logiSplits) + } else { + .QuartetMC(charInt, logiSplits, normalize) + } + eNum <- base[["concordant"]] + eDen <- base[["decisive"]] + } + if (return == "default") { if (isTRUE(weight)) { # Sum numerator and denominator across sites (columns), then divide @@ -898,37 +907,164 @@ QuartetConcordance <- function( NA_real_, split_sums_num / split_sums_den ) + if (doNorm) { + bDen <- rowSums(eDen) + base <- ifelse(bDen == 0, NA_real_, rowSums(eNum) / bDen) + ret <- .RezeroGuarded(ret, base) + } } else { # Mean of ratios per site # Avoid division by zero (0/0 -> NaN -> NA handled by na.rm) ratios <- num / den # Replace NaN/Inf with NA for rowMeans calculation ratios[!is.finite(ratios)] <- NA - ret <- rowMeans(ratios, na.rm = TRUE) + if (doNorm) { + bRatios <- eNum / eDen + bRatios[!is.finite(bRatios)] <- NA + # Average observed and baseline over the same (split, char) cells. + naMask <- is.na(ratios) | is.na(bRatios) + ratios[naMask] <- NA + bRatios[naMask] <- NA + ret <- .RezeroGuarded(rowMeans(ratios, na.rm = TRUE), + rowMeans(bRatios, na.rm = TRUE)) + } else { + ret <- rowMeans(ratios, na.rm = TRUE) + } } setNames(ret, names(splits)) } else { # return = "char" - p <- num / den if (isTRUE(weight)) { - vapply( - seq_len(dim(num)[[2]]), - function(i) { - weighted.mean(num[, i] / den[, i], den[, i]) - }, - double(1) - ) + if (doNorm) { + obs <- vapply(seq_len(ncol(num)), function(i) { + d <- sum(den[, i]); if (d == 0) NA_real_ else sum(num[, i]) / d + }, double(1)) + b <- vapply(seq_len(ncol(eNum)), function(i) { + d <- sum(eDen[, i]); if (d == 0) NA_real_ else sum(eNum[, i]) / d + }, double(1)) + .RezeroGuarded(obs, b) + } else { + vapply( + seq_len(dim(num)[[2]]), + function(i) { + weighted.mean(num[, i] / den[, i], den[, i]) + }, + double(1) + ) + } } else { - vapply( - seq_len(dim(num)[[2]]), - function(i) { - mean(num[den[, i] > 0, i] / den[den[, i] > 0, i]) - }, - double(1) - ) + if (doNorm) { + sObs <- ifelse(den > 0, num / den, NA_real_) + sBase <- ifelse(eDen > 0, eNum / eDen, NA_real_) + naMask <- is.na(sObs) | is.na(sBase) + sObs[naMask] <- NA_real_ + sBase[naMask] <- NA_real_ + obs <- vapply(seq_len(ncol(sObs)), function(i) { + m <- mean(sObs[, i], na.rm = TRUE); if (is.nan(m)) NA_real_ else m + }, double(1)) + b <- vapply(seq_len(ncol(sBase)), function(i) { + m <- mean(sBase[, i], na.rm = TRUE); if (is.nan(m)) NA_real_ else m + }, double(1)) + .RezeroGuarded(obs, b) + } else { + vapply( + seq_len(dim(num)[[2]]), + function(i) { + mean(num[den[, i] > 0, i] / den[den[, i] > 0, i]) + }, + double(1) + ) + } + } + } +} + +# Exact / Monte-Carlo expected raw concordant + decisive quartet counts under +# the fixed-marginal null (the raw-currency analogue of `.CharTritExpect` / +# `.CharTritMC`). Per state-pair, conc = C(p,2)C(s,2) + C(q,2)C(r,2) and +# dec = conc + p*q*r*s are polynomials in the cells (no floors), so their exact +# expectation is a clean sum over the same trivariate-hypergeometric pmf. +.ExpectedQuartetCache <- new.env(hash = TRUE, parent = emptyenv()) + +.ExpectedQuartet <- function(nI, nJ, M, t) { + key <- paste(nI, nJ, M, t, sep = ",") + cached <- .ExpectedQuartetCache[[key]] + if (!is.null(cached)) { + return(cached) + } + choose2 <- function(z) if (z < 2) 0 else z * (z - 1) / 2 + nOther <- t - nI - nJ + logDen <- lchoose(t, M) + eConc <- 0 + eDec <- 0 + for (p in 0:min(nI, M)) { + for (r in 0:min(nJ, M - p)) { + oA <- M - p - r + if (oA < 0 || oA > nOther) { + next + } + prob <- exp(lchoose(nI, p) + lchoose(nJ, r) + lchoose(nOther, oA) - logDen) + if (prob <= 0) { + next + } + q <- nI - p + s <- nJ - r + conc <- choose2(p) * choose2(s) + choose2(q) * choose2(r) + eConc <- eConc + prob * conc + eDec <- eDec + prob * (conc + p * q * r * s) + } + } + res <- c(eConc, eDec) + .ExpectedQuartetCache[[key]] <- res + res +} + +.QuartetExpect <- function(charInt, logiSplits) { + nSplit <- ncol(logiSplits) + nChar <- ncol(charInt) + eConc <- eDec <- matrix(0, nSplit, nChar) + for (ci in seq_len(nChar)) { + col <- charInt[, ci] + scored <- !is.na(col) + states <- sort(unique(col[scored])) + nStates <- length(states) + if (nStates < 2L) { + next + } + tc <- sum(scored) + mSideA <- colSums(logiSplits & scored) + uM <- unique(mSideA) + idx <- match(mSideA, uM) + cnt <- tabulate(match(col[scored], states), nStates) + for (a in seq_len(nStates - 1L)) { + for (b in seq(a + 1L, nStates)) { + e <- vapply(uM, function(M) .ExpectedQuartet(cnt[a], cnt[b], M, tc), + double(2)) + eConc[, ci] <- eConc[, ci] + e[1, idx] + eDec[, ci] <- eDec[, ci] + e[2, idx] + } + } + } + list(concordant = eConc, decisive = eDec) +} + +.QuartetMC <- function(charInt, logiSplits, nRelabel) { + nSplit <- ncol(logiSplits) + nChar <- ncol(charInt) + accConc <- accDec <- matrix(0, nSplit, nChar) + for (s in seq_len(nRelabel)) { + shuffled <- charInt + for (ci in seq_len(nChar)) { + col <- charInt[, ci] + scored <- !is.na(col) + shuffled[scored, ci] <- sample(col[scored]) } + kc <- quartet_concordance(logiSplits, shuffled) + accConc <- accConc + kc[["concordant"]] + accDec <- accDec + kc[["decisive"]] } + list(concordant = accConc / nRelabel, decisive = accDec / nRelabel) } # Nelson-Ladiges fractional ("trit") currency for QuartetConcordance(). diff --git a/dev/plans/frac-quart-weight-agreement.md b/dev/plans/frac-quart-weight-agreement.md index 6840a9142..7fe19b3d4 100644 --- a/dev/plans/frac-quart-weight-agreement.md +++ b/dev/plans/frac-quart-weight-agreement.md @@ -253,13 +253,18 @@ char FQ(c) = Σ_k M(k,c)·Q(k,c) / Σ_k M(k,c) never inside the measure. The `z → 1` blow-up is guarded (→ `NA`). `normalize` defaults to `FALSE`, so the published measure is unchanged unless the user opts in. - - **Raw `quartet` unit: deferred.** `normalize` currently errors for - `unit = "quartet"`. Re-baselining the raw currency needs the per-state-pair - cells (computed inside the C++ kernel but not exposed) and a decision on its - null (fixed-marginal vs the flat `1/3`); both are folded into the C++ port - (task 2) and the manuscript-figure review (task 3), which will compare - trit-with vs trit-without correction and decide whether the raw unit needs it - too. + - **Raw `quartet` unit: DONE.** `normalize` now works for `unit = "quartet"` + too, using the *same* fixed-marginal null. Per state-pair `conc` and `dec` + are polynomials in the cells (no floors), so the exact expectation + `E[conc]`, `E[dec]` is a clean sum over the same trivariate-hypergeometric + pmf as `.ExpectedTrit` (`.ExpectedQuartet` / `.QuartetExpect`); MC via + `.QuartetMC` reshuffles tokens and re-scores through the kernel. The pooled + `conc/dec` ratio is re-zeroed against `E[conc]/E[dec]` (weighted) or with + the same cell-matching as trits (unweighted). Validated exact≈MC across + edge/char × weight × binary/multistate (`test-Concordance.R`), and + `normalize = FALSE` is byte-identical to the published raw measure. The flat + `1/3` baseline stays rejected (it does not adapt to marginals; the + fixed-marginal null generalises it). - **`M = min(W_c,W_k)`** (settled): the pooling amount is the *shared* information, so it must be symmetric between character and split (both `return`s pool by the same `M`); `M = W_k` would make the char return weight a diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index cc617bca8..b420da41c 100644 --- a/man/ConcordanceTable.Rd +++ b/man/ConcordanceTable.Rd @@ -48,19 +48,19 @@ floor. expectation and negative values fall below it. \code{ClusteringConcordance()} uses an analytical approximation to the expected mutual information (fast and generally accurate for large trees, ~200+ taxa, but neglecting -correlation between splits); \code{QuartetConcordance(unit = "trit")} uses the -exact hypergeometric expectation. +correlation between splits); \code{QuartetConcordance()} uses the exact +hypergeometric expectation (for either \code{unit}). \item a positive integer \code{n}: estimate that expectation by Monte Carlo. \code{ClusteringConcordance()} fits each character to \code{n} random trees (and returns Monte-Carlo standard errors, more accurate for small trees where -the analytical approximation is biased); \code{QuartetConcordance(unit = "trit")} averages over \code{n} random reassignments of each character's -tokens. +the analytical approximation is biased); \code{QuartetConcordance()} averages +over \code{n} random reassignments of each character's tokens. } In all cases 1 corresponds to the maximum attainable value. For -\code{QuartetConcordance()}, chance correction is implemented for \code{unit = "trit"} -only, and the returned values are unclamped (they may fall below \eqn{-1}); -clamp to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} +\code{QuartetConcordance()} (either \code{unit}), chance-corrected values are returned +unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before +plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{plot}{Logical specifying whether to draw the plot.} diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index 3dd1a5604..f8b3c6cc0 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -56,19 +56,19 @@ floor. expectation and negative values fall below it. \code{ClusteringConcordance()} uses an analytical approximation to the expected mutual information (fast and generally accurate for large trees, ~200+ taxa, but neglecting -correlation between splits); \code{QuartetConcordance(unit = "trit")} uses the -exact hypergeometric expectation. +correlation between splits); \code{QuartetConcordance()} uses the exact +hypergeometric expectation (for either \code{unit}). \item a positive integer \code{n}: estimate that expectation by Monte Carlo. \code{ClusteringConcordance()} fits each character to \code{n} random trees (and returns Monte-Carlo standard errors, more accurate for small trees where -the analytical approximation is biased); \code{QuartetConcordance(unit = "trit")} averages over \code{n} random reassignments of each character's -tokens. +the analytical approximation is biased); \code{QuartetConcordance()} averages +over \code{n} random reassignments of each character's tokens. } In all cases 1 corresponds to the maximum attainable value. For -\code{QuartetConcordance()}, chance correction is implemented for \code{unit = "trit"} -only, and the returned values are unclamped (they may fall below \eqn{-1}); -clamp to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} +\code{QuartetConcordance()} (either \code{unit}), chance-corrected values are returned +unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before +plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{weight}{Logical specifying whether to weight sites according to the number of quartets they are decisive for.} diff --git a/tests/testthat/test-Concordance.R b/tests/testthat/test-Concordance.R index c1e100f16..b455f6012 100644 --- a/tests/testthat/test-Concordance.R +++ b/tests/testthat/test-Concordance.R @@ -252,10 +252,10 @@ test_that("QuartetConcordance() unit = 'trit' normalize contract", { expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = "x"), "positive integer") - # Chance correction is implemented for the trit currency only, for now. - expect_error( - QuartetConcordance(tree, ident, unit = "quartet", normalize = TRUE), - "trit") + # Chance correction also works for the raw quartet currency. + qc <- QuartetConcordance(tree, ident, unit = "quartet", normalize = TRUE) + expect_type(qc, "double") + expect_false(anyNA(qc)) # An identical (displayed) split still scores exactly 1 after re-zeroing: # .Rezero(1, z) == 1, so chance correction lifts the floor without moving the @@ -299,6 +299,39 @@ test_that("QuartetConcordance() unit = 'trit' exact baseline matches MC", { } }) +test_that("QuartetConcordance() unit = 'quartet' exact baseline matches MC", { + tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") + dat <- MatrixToPhyDat(matrix( + c(0, 0, 0, 0, 0, 0, 1, 1, 1, # displays {t7,t8,t9} + 0, 0, 0, 0, 0, 0, 0, 1, 1, # nested + 0, 0, 1, 1, 0, 0, 1, 1, 0, # crossing (binary) + 0, 0, 0, 1, 1, 1, 2, 2, 2), 9, # 3-state: multistate decisive path + dimnames = list(paste0("t", 1:9), NULL))) + + # The raw-currency chance baseline (exact E[conc]/E[dec] from the + # hypergeometric pmf) matches the Monte-Carlo tip-shuffle baseline within MC + # error, across edge/char and both weightings, incl. the multistate path. + for (ret in c("edge", "char")) { + for (w in c(TRUE, FALSE)) { + exact <- QuartetConcordance(tree, dat, unit = "quartet", return = ret, + weight = w, normalize = TRUE) + set.seed(1) + mc <- QuartetConcordance(tree, dat, unit = "quartet", return = ret, + weight = w, normalize = 8000) + expect_identical(is.na(exact), is.na(mc)) + expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04) + } + } + + # normalize = FALSE leaves the published raw measure untouched. + for (ret in c("edge", "char")) for (w in c(TRUE, FALSE)) { + expect_identical( + QuartetConcordance(tree, dat, unit = "quartet", return = ret, weight = w), + QuartetConcordance(tree, dat, unit = "quartet", return = ret, weight = w, + normalize = FALSE)) + } +}) + test_that("QuartetConcordance() unit = 'trit' correction: conflict below random", { tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") support <- MatrixToPhyDat(matrix( # agrees with the tree: displays {t7,t8,t9} From 8461fc6bd447ff04c36974b636d7294037439b1d Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 10 Jul 2026 15:14:44 +0100 Subject: [PATCH 05/16] perf(concordance): port chance-baseline expectation to C++ 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 = ) 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 --- R/Concordance.R | 167 ++------------ R/RcppExports.R | 8 + dev/benchmarks/frac-quart/cpp_expect_parity.R | 108 +++++++++ src/RcppExports.cpp | 24 ++ src/TreeSearch-init.c | 4 + src/concordance_expect.cpp | 215 ++++++++++++++++++ 6 files changed, 382 insertions(+), 144 deletions(-) create mode 100644 dev/benchmarks/frac-quart/cpp_expect_parity.R create mode 100644 src/concordance_expect.cpp diff --git a/R/Concordance.R b/R/Concordance.R index ea2f353b7..51a484d22 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -980,73 +980,15 @@ QuartetConcordance <- function( } } -# Exact / Monte-Carlo expected raw concordant + decisive quartet counts under -# the fixed-marginal null (the raw-currency analogue of `.CharTritExpect` / -# `.CharTritMC`). Per state-pair, conc = C(p,2)C(s,2) + C(q,2)C(r,2) and -# dec = conc + p*q*r*s are polynomials in the cells (no floors), so their exact -# expectation is a clean sum over the same trivariate-hypergeometric pmf. -.ExpectedQuartetCache <- new.env(hash = TRUE, parent = emptyenv()) - -.ExpectedQuartet <- function(nI, nJ, M, t) { - key <- paste(nI, nJ, M, t, sep = ",") - cached <- .ExpectedQuartetCache[[key]] - if (!is.null(cached)) { - return(cached) - } - choose2 <- function(z) if (z < 2) 0 else z * (z - 1) / 2 - nOther <- t - nI - nJ - logDen <- lchoose(t, M) - eConc <- 0 - eDec <- 0 - for (p in 0:min(nI, M)) { - for (r in 0:min(nJ, M - p)) { - oA <- M - p - r - if (oA < 0 || oA > nOther) { - next - } - prob <- exp(lchoose(nI, p) + lchoose(nJ, r) + lchoose(nOther, oA) - logDen) - if (prob <= 0) { - next - } - q <- nI - p - s <- nJ - r - conc <- choose2(p) * choose2(s) + choose2(q) * choose2(r) - eConc <- eConc + prob * conc - eDec <- eDec + prob * (conc + p * q * r * s) - } - } - res <- c(eConc, eDec) - .ExpectedQuartetCache[[key]] <- res - res -} - +# Expected raw concordant + decisive quartet counts under the fixed-marginal +# null. The exact expectation (per state-pair, an exact sum over the +# trivariate-hypergeometric pmf) is computed in C++ (`quartet_expect` in +# src/concordance_expect.cpp) for speed; `.QuartetExpect` is a thin R wrapper. +# The `.QuartetMC` Monte-Carlo baseline reshuffles tokens and re-scores through +# the C++ `quartet_concordance` kernel. (The C++ path is validated bit-for-bit +# against the earlier R implementation; see dev/benchmarks/frac-quart.) .QuartetExpect <- function(charInt, logiSplits) { - nSplit <- ncol(logiSplits) - nChar <- ncol(charInt) - eConc <- eDec <- matrix(0, nSplit, nChar) - for (ci in seq_len(nChar)) { - col <- charInt[, ci] - scored <- !is.na(col) - states <- sort(unique(col[scored])) - nStates <- length(states) - if (nStates < 2L) { - next - } - tc <- sum(scored) - mSideA <- colSums(logiSplits & scored) - uM <- unique(mSideA) - idx <- match(mSideA, uM) - cnt <- tabulate(match(col[scored], states), nStates) - for (a in seq_len(nStates - 1L)) { - for (b in seq(a + 1L, nStates)) { - e <- vapply(uM, function(M) .ExpectedQuartet(cnt[a], cnt[b], M, tc), - double(2)) - eConc[, ci] <- eConc[, ci] + e[1, idx] - eDec[, ci] <- eDec[, ci] + e[2, idx] - } - } - } - list(concordant = eConc, decisive = eDec) + quartet_expect(logiSplits, charInt) } .QuartetMC <- function(charInt, logiSplits, nRelabel) { @@ -1109,6 +1051,15 @@ QuartetConcordance <- function( doNorm <- !isFALSE(normalize) if (doNorm) { baseNumEdge <- baseNumChar <- baseDenM <- matrix(0, nSplit, nChar) + if (isTRUE(normalize)) { + # Exact hypergeometric expectation for all characters at once (C++; + # `trit_expect` in src/concordance_expect.cpp). Uninformative characters + # yield zero pools, matching the observed loop's `wcTot > 0` guard. + be <- trit_expect(logiSplits, charInt) + baseNumEdge <- be[["numEdge"]] + baseNumChar <- be[["numChar"]] + baseDenM <- be[["denM"]] + } } for (ci in seq_len(nChar)) { @@ -1118,12 +1069,8 @@ QuartetConcordance <- function( numChar[, ci] <- obs[["numChar"]] denM[, ci] <- obs[["denM"]] wcTot[ci] <- obs[["wcTot"]] - if (doNorm && obs[["wcTot"]] > 0) { - base <- if (isTRUE(normalize)) { - .CharTritExpect(col, logiSplits, pos) # exact hypergeometric - } else { - .CharTritMC(col, logiSplits, pos, normalize) # Monte-Carlo tip-shuffle - } + if (doNorm && !isTRUE(normalize) && obs[["wcTot"]] > 0) { + base <- .CharTritMC(col, logiSplits, pos, normalize) # Monte-Carlo shuffle baseNumEdge[, ci] <- base[["numEdge"]] baseNumChar[, ci] <- base[["numChar"]] baseDenM[, ci] <- base[["denM"]] @@ -1229,78 +1176,10 @@ QuartetConcordance <- function( list(numEdge = numEdge, numChar = numChar, denM = denM, wcTot = wcTot) } -# Exact expected trit contributions under the fixed-marginal (hypergeometric) -# null: each split's scored side-A size `M` and each pair's state counts -# `n_i, n_j` are held fixed while the character's tokens are reassigned at -# random across the `t` scored leaves. Only `A`, and (for multistate pairs) -# `wk` and `m`, vary; so we accumulate E[m], E[m * A / wk] and E[m * A / wc] -# from the exact pmf, matching `.CharTritContrib`'s pool definitions. -.CharTritExpect <- function(col, logiSplits, pos) { - nSplit <- ncol(logiSplits) - numEdge <- numChar <- denM <- numeric(nSplit) - scored <- !is.na(col) - tc <- sum(scored) - mSideA <- colSums(logiSplits & scored) # M per split - uM <- unique(mSideA) - idx <- match(mSideA, uM) - states <- sort(unique(col[scored])) - nStates <- length(states) - if (nStates >= 2L) { - cnt <- tabulate(match(col[scored], states), nStates) # state counts - for (a in seq_len(nStates - 1L)) { - for (b in seq(a + 1L, nStates)) { - et <- vapply(uM, function(M) .ExpectedTrit(cnt[a], cnt[b], M, tc), - double(3)) - denM <- denM + et[1, idx] - numEdge <- numEdge + et[2, idx] - numChar <- numChar + et[3, idx] - } - } - } - list(numEdge = numEdge, numChar = numChar, denM = denM) -} - -.ExpectedTritCache <- new.env(hash = TRUE, parent = emptyenv()) - -# Exact E[m], E[m * A / wk], E[m * A / wc] for one state-pair, over the -# trivariate hypergeometric placement of `nI` state-i and `nJ` state-j leaves -# (and `t - nI - nJ` other leaves) among the `M` leaves on side A of a split. -.ExpectedTrit <- function(nI, nJ, M, t) { - key <- paste(nI, nJ, M, t, sep = ",") - cached <- .ExpectedTritCache[[key]] - if (!is.null(cached)) { - return(cached) - } - pos1 <- function(z) if (z < 0) 0 else z - nOther <- t - nI - nJ - tP <- nI + nJ # taxa in this pair - wc <- pos1(nI - 1) * pos1(nJ - 1) # pair's character content (fixed) - logDen <- lchoose(t, M) - acc <- c(0, 0, 0) # E[m], E[m*A/wk], E[m*A/wc] - for (p in 0:min(nI, M)) { - for (r in 0:min(nJ, M - p)) { - oA <- M - p - r # other-state leaves on side A - if (oA < 0 || oA > nOther) { - next - } - prob <- exp(lchoose(nI, p) + lchoose(nJ, r) + lchoose(nOther, oA) - logDen) - if (prob <= 0) { - next - } - q <- nI - p - s <- nJ - r - mA <- p + r # this pair's taxa on side A (random) - A <- pos1(p - 1) * pos1(s - 1) + pos1(q - 1) * pos1(r - 1) - wk <- pos1(mA - 1) * pos1(tP - mA - 1) # pair's split content (random) - m <- min(wc, wk) - acc[1] <- acc[1] + prob * m - if (wk > 0) acc[2] <- acc[2] + prob * m * A / wk - if (wc > 0) acc[3] <- acc[3] + prob * m * A / wc - } - } - .ExpectedTritCache[[key]] <- acc - acc -} +# The exact trit expectation E[m], E[m*A/wk], E[m*A/wc] is computed in C++ +# (`trit_expect` in src/concordance_expect.cpp) for all characters at once, and +# consumed directly by `.TritConcordance`; there is no per-character R exact +# helper. The Monte-Carlo baseline below is the opt-in `normalize = ` path. # Monte-Carlo baseline: average `.CharTritContrib` over `nRelabel` random # reassignments of the character's tokens across its scored leaves (the split, diff --git a/R/RcppExports.R b/R/RcppExports.R index afce48989..4d29eae1d 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -20,6 +20,14 @@ MaddisonSlatkin_clear_cache <- function() { invisible(.Call(`_TreeSearch_MaddisonSlatkin_clear_cache`)) } +quartet_expect <- function(splits, characters) { + .Call(`_TreeSearch_quartet_expect`, splits, characters) +} + +trit_expect <- function(splits, characters) { + .Call(`_TreeSearch_trit_expect`, splits, characters) +} + #' Expected mutual information between two partitions #' #' Computes the mutual information expected purely by chance between two diff --git a/dev/benchmarks/frac-quart/cpp_expect_parity.R b/dev/benchmarks/frac-quart/cpp_expect_parity.R new file mode 100644 index 000000000..0e8afba7f --- /dev/null +++ b/dev/benchmarks/frac-quart/cpp_expect_parity.R @@ -0,0 +1,108 @@ +# Parity guard for the C++ chance-baseline kernels `quartet_expect` and +# `trit_expect` (src/concordance_expect.cpp), which compute the exact expected +# concordant/decisive quartet counts (raw) and the expected trit pools +# (E[m], E[m*A/wk], E[m*A/wc]) under the fixed-marginal (hypergeometric) null. +# +# The kernels replaced an earlier pure-R implementation; this script keeps a +# self-contained R reference (the same trivariate-hypergeometric double sum) and +# checks the C++ output matches it to machine precision across random binary, +# multistate and missing-data matrices. Run after any change to the kernels. + +suppressMessages(pkgload::load_all( + "C:/Users/pjjg18/GitHub/worktrees/TreeSearch/frac-quart-cpp", + quiet = TRUE, recompile = TRUE)) +suppressMessages(library(TreeTools)) + +# Build logiSplits + charInt exactly as QuartetConcordance() does internally. +build_inputs <- function(tree, dat) { + tl <- intersect(TipLabels(tree), names(dat)) + dat <- dat[tl, drop = FALSE] + splits <- as.Splits(tree, dat) + logiSplits <- vapply(seq_along(splits), function(i) as.logical(splits[[i]]), + logical(NTip(dat))) + contrast <- attr(dat, "contrast"); lv <- attr(dat, "allLevels") + isInapp <- lv == "-" + isAmbig <- rowSums(contrast[, colnames(contrast) != "-", drop = FALSE]) > 1 + isGrouping <- !isAmbig & !isInapp + gCols <- apply(contrast[isGrouping, , drop = FALSE] > 0, 1, which) + l2i <- rep(NA_integer_, length(lv)); l2i[isGrouping] <- as.integer(gCols) + ch <- PhyDatToMatrix(dat) + charInt <- array(l2i[match(ch, lv)], dim = dim(ch), dimnames = dimnames(ch)) + list(logiSplits = logiSplits, charInt = charInt) +} + +# ---- self-contained R reference (trivariate-hypergeometric double sum) ---- +pos1 <- function(z) if (z < 0) 0 else z +choose2 <- function(z) if (z < 2) 0 else z * (z - 1) / 2 + +ref_pair <- function(nI, nJ, M, t) { # returns c(Econc, Edec, Em, Ewk, Ewc) + nOther <- t - nI - nJ; tP <- nI + nJ + wc <- pos1(nI - 1) * pos1(nJ - 1) + logDen <- lchoose(t, M) + acc <- c(0, 0, 0, 0, 0) + for (p in 0:min(nI, M)) for (r in 0:min(nJ, M - p)) { + oA <- M - p - r + if (oA < 0 || oA > nOther) next + prob <- exp(lchoose(nI, p) + lchoose(nJ, r) + lchoose(nOther, oA) - logDen) + if (prob <= 0) next + q <- nI - p; s <- nJ - r; mA <- p + r + conc <- choose2(p) * choose2(s) + choose2(q) * choose2(r) + dec <- conc + p * q * r * s + A <- pos1(p - 1) * pos1(s - 1) + pos1(q - 1) * pos1(r - 1) + wk <- pos1(mA - 1) * pos1(tP - mA - 1); m <- min(wc, wk) + acc <- acc + prob * c(conc, dec, m, if (wk > 0) m * A / wk else 0, + if (wc > 0) m * A / wc else 0) + } + acc +} + +ref_expect <- function(charInt, logiSplits) { + nSplit <- ncol(logiSplits); nChar <- ncol(charInt) + out <- list(eConc = matrix(0, nSplit, nChar), eDec = matrix(0, nSplit, nChar), + eM = matrix(0, nSplit, nChar), eWk = matrix(0, nSplit, nChar), + eWc = matrix(0, nSplit, nChar)) + for (ci in seq_len(nChar)) { + col <- charInt[, ci]; scored <- !is.na(col) + states <- sort(unique(col[scored])); nStates <- length(states) + if (nStates < 2L) next + tc <- sum(scored); mSideA <- colSums(logiSplits & scored) + uM <- unique(mSideA); idx <- match(mSideA, uM) + cnt <- tabulate(match(col[scored], states), nStates) + for (a in seq_len(nStates - 1L)) for (b in seq(a + 1L, nStates)) { + e <- vapply(uM, function(M) ref_pair(cnt[a], cnt[b], M, tc), double(5)) + out$eConc[, ci] <- out$eConc[, ci] + e[1, idx] + out$eDec[, ci] <- out$eDec[, ci] + e[2, idx] + out$eM[, ci] <- out$eM[, ci] + e[3, idx] + out$eWk[, ci] <- out$eWk[, ci] + e[4, idx] + out$eWc[, ci] <- out$eWc[, ci] + e[5, idx] + } + } + out +} + +set.seed(7) +worst <- 0 +for (rep in 1:10) { + nTip <- sample(9:44, 1) + tree <- RandomTree(nTip, root = FALSE) + m <- replicate(sample(5:30, 1), { + ns <- sample(2:4, 1) + x <- as.character(sample(0:(ns - 1), nTip, replace = TRUE)) + if (runif(1) < 0.3) x[sample(nTip, sample(1:3, 1))] <- "?" + x + }) + rownames(m) <- paste0("t", seq_len(nTip)) + io <- build_inputs(tree, MatrixToPhyDat(m)) + ls <- io$logiSplits; ci <- io$charInt + R <- ref_expect(ci, ls) + cq <- quartet_expect(ls, ci); ct <- trit_expect(ls, ci) + d <- max(abs(cq$concordant - R$eConc), abs(cq$decisive - R$eDec), + abs(ct$numEdge - R$eWk), abs(ct$numChar - R$eWc), + abs(ct$denM - R$eM)) + worst <- max(worst, d) + cat(sprintf("rep %2d: nTip=%2d nChar=%2d max|C++ - R| = %.2e\n", + rep, nTip, ncol(ci), d)) +} +cat(sprintf("\nWorst |C++ - R| over all kernels: %.3e\n", worst)) +stopifnot(worst < 1e-9) +cat("PASS: C++ expectation kernels match the R reference to machine precision.\n") diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index dc52e9b59..f46299b2e 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -31,6 +31,30 @@ BEGIN_RCPP return R_NilValue; END_RCPP } +// quartet_expect +List quartet_expect(const LogicalMatrix splits, const IntegerMatrix characters); +RcppExport SEXP _TreeSearch_quartet_expect(SEXP splitsSEXP, SEXP charactersSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const LogicalMatrix >::type splits(splitsSEXP); + Rcpp::traits::input_parameter< const IntegerMatrix >::type characters(charactersSEXP); + rcpp_result_gen = Rcpp::wrap(quartet_expect(splits, characters)); + return rcpp_result_gen; +END_RCPP +} +// trit_expect +List trit_expect(const LogicalMatrix splits, const IntegerMatrix characters); +RcppExport SEXP _TreeSearch_trit_expect(SEXP splitsSEXP, SEXP charactersSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const LogicalMatrix >::type splits(splitsSEXP); + Rcpp::traits::input_parameter< const IntegerMatrix >::type characters(charactersSEXP); + rcpp_result_gen = Rcpp::wrap(trit_expect(splits, characters)); + return rcpp_result_gen; +END_RCPP +} // expected_mi double expected_mi(const IntegerVector& ni, const IntegerVector& nj); RcppExport SEXP _TreeSearch_expected_mi(SEXP niSEXP, SEXP njSEXP) { diff --git a/src/TreeSearch-init.c b/src/TreeSearch-init.c index 7b166a616..890fd4255 100644 --- a/src/TreeSearch-init.c +++ b/src/TreeSearch-init.c @@ -21,6 +21,8 @@ extern SEXP _TreeSearch_mi_key(SEXP, SEXP); // extern SEXP _TreeSearch_astar_search_r(SEXP, SEXP, SEXP); extern SEXP _TreeSearch_quartet_concordance(SEXP, SEXP); +extern SEXP _TreeSearch_quartet_expect(SEXP, SEXP); +extern SEXP _TreeSearch_trit_expect(SEXP, SEXP); extern SEXP _TreeSearch_ts_fitch_score(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _TreeSearch_ts_na_char_steps(SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _TreeSearch_ts_char_steps(SEXP, SEXP, SEXP, SEXP, SEXP); @@ -75,6 +77,8 @@ static const R_CallMethodDef callMethods[] = { // {"_TreeSearch_astar_search_r", (DL_FUNC) &_TreeSearch_astar_search_r, 3}, {"_TreeSearch_quartet_concordance",(DL_FUNC) &_TreeSearch_quartet_concordance, 2}, + {"_TreeSearch_quartet_expect", (DL_FUNC) &_TreeSearch_quartet_expect, 2}, + {"_TreeSearch_trit_expect", (DL_FUNC) &_TreeSearch_trit_expect, 2}, {"_TreeSearch_ts_fitch_score", (DL_FUNC) &_TreeSearch_ts_fitch_score, 12}, {"_TreeSearch_ts_na_char_steps", (DL_FUNC) &_TreeSearch_ts_na_char_steps, 5}, {"_TreeSearch_ts_char_steps", (DL_FUNC) &_TreeSearch_ts_char_steps, 5}, diff --git a/src/concordance_expect.cpp b/src/concordance_expect.cpp new file mode 100644 index 000000000..4faefac8e --- /dev/null +++ b/src/concordance_expect.cpp @@ -0,0 +1,215 @@ +#include +#include +#include +#include + +using namespace Rcpp; + +// Expected concordance under the fixed-marginal (hypergeometric) null used by +// QuartetConcordance(normalize = TRUE). For one character each token is +// reassigned at random across the scored leaves, holding the state counts and +// the split sizes fixed. For a state-pair (i, j) with counts nI, nJ over t +// scored leaves, of which M are on side A of a split, the 2x2 cell vector +// (p, q, r, s) is trivariate-hypergeometric: +// P(p, r) = C(nI,p) C(nJ,r) C(t-nI-nJ, M-p-r) / C(t, M) +// with q = nI - p, s = nJ - r, and mA = p + r on side A. These kernels sum the +// per-pair contributions exactly (double sum over p, r), mirroring the R +// reference `.ExpectedQuartet`/`.QuartetExpect` and `.ExpectedTrit`/ +// `.CharTritExpect` in R/Concordance.R. The p-outer / r-inner accumulation +// order matches the R code so the results agree to machine precision. + +// C(z, 2) with the same z < 2 -> 0 guard as the R `choose2`. +static inline double choose2(double z) { + return z < 2.0 ? 0.0 : z * (z - 1.0) / 2.0; +} + +// (z - 1) floored at 0, matching the R `pos(z - 1)`. +static inline double posm1(double z) { + return z < 1.0 ? 0.0 : z - 1.0; +} + +// Per-character scored state counts + per-split side-A size, shared by both +// kernels. `stateOf[t]` is the compact 0-based state index of a scored leaf, or +// -1 if the leaf is NA for this character. +struct CharCounts { + int tc; // number of scored leaves + std::vector cnt; // count of each compact state + std::vector mSideA; // per split: scored leaves on side A +}; + +static CharCounts char_counts(const LogicalMatrix& splits, + const IntegerMatrix& characters, int c) { + const int n_taxa = splits.nrow(); + const int n_splits = splits.ncol(); + CharCounts out; + out.tc = 0; + out.mSideA.assign(n_splits, 0); + + // Map raw state values to compact 0-based indices (sorted ascending, as R's + // sort(unique(...)) does; the pair loop is symmetric so order only affects + // which is "i" vs "j", not the sums). + std::vector raw; + raw.reserve(n_taxa); + std::vector stateOf(n_taxa, -1); + for (int t = 0; t < n_taxa; ++t) { + int st = characters(t, c); + if (IntegerVector::is_na(st)) continue; + ++out.tc; + raw.push_back(st); + } + std::sort(raw.begin(), raw.end()); + raw.erase(std::unique(raw.begin(), raw.end()), raw.end()); + out.cnt.assign(raw.size(), 0); + + for (int t = 0; t < n_taxa; ++t) { + int st = characters(t, c); + if (IntegerVector::is_na(st)) continue; + int idx = int(std::lower_bound(raw.begin(), raw.end(), st) - raw.begin()); + stateOf[t] = idx; + ++out.cnt[idx]; + for (int s = 0; s < n_splits; ++s) { + if (splits(t, s)) ++out.mSideA[s]; + } + } + return out; +} + +// Exact E[conc], E[dec] for one state-pair over all needed side-A sizes. +// Returns, for each M in 0..tc, the pair's (Econc, Edec); callers pick M[split]. +static void expected_quartet_by_M(int nI, int nJ, int tc, + std::vector& eConc, + std::vector& eDec) { + const int nOther = tc - nI - nJ; + eConc.assign(tc + 1, 0.0); + eDec.assign(tc + 1, 0.0); + // Hoist the log-choose terms out of the M / p / r loops (the dominant cost). + std::vector lcI(nI + 1), lcJ(nJ + 1), lcO(nOther + 1), lcT(tc + 1); + for (int p = 0; p <= nI; ++p) lcI[p] = R::lchoose(nI, p); + for (int r = 0; r <= nJ; ++r) lcJ[r] = R::lchoose(nJ, r); + for (int o = 0; o <= nOther; ++o) lcO[o] = R::lchoose(nOther, o); + for (int M = 0; M <= tc; ++M) lcT[M] = R::lchoose(tc, M); + for (int M = 0; M <= tc; ++M) { + const double lt = lcT[M]; + double aConc = 0.0, aDec = 0.0; + const int pMax = std::min(nI, M); + for (int p = 0; p <= pMax; ++p) { + const int q = nI - p; + const double c_p = choose2(p), c_q = choose2(q), lp = lcI[p]; + const int rMax = std::min(nJ, M - p); + for (int r = 0; r <= rMax; ++r) { + const int oA = M - p - r; + if (oA < 0 || oA > nOther) continue; + const double prob = std::exp(lp + lcJ[r] + lcO[oA] - lt); + if (prob <= 0.0) continue; + const int s = nJ - r; + const double conc = c_p * choose2(s) + c_q * choose2(r); + aConc += prob * conc; + aDec += prob * (conc + double(p) * q * r * s); + } + } + eConc[M] = aConc; + eDec[M] = aDec; + } +} + +// [[Rcpp::export]] +List quartet_expect(const LogicalMatrix splits, const IntegerMatrix characters) { + const int n_splits = splits.ncol(); + const int n_chars = characters.ncol(); + NumericMatrix eConc(n_splits, n_chars), eDec(n_splits, n_chars); + + for (int c = 0; c < n_chars; ++c) { + const CharCounts cc = char_counts(splits, characters, c); + const int nStates = int(cc.cnt.size()); + if (nStates < 2) continue; + std::vector ec, ed; + for (int a = 0; a < nStates - 1; ++a) { + for (int b = a + 1; b < nStates; ++b) { + expected_quartet_by_M(cc.cnt[a], cc.cnt[b], cc.tc, ec, ed); + for (int s = 0; s < n_splits; ++s) { + const int M = cc.mSideA[s]; + eConc(s, c) += ec[M]; + eDec(s, c) += ed[M]; + } + } + } + } + return List::create(_["concordant"] = eConc, _["decisive"] = eDec); +} + +// Exact E[m], E[m*A/wk], E[m*A/wc] for one trit state-pair over all side-A +// sizes. wc = (nI-1)+ (nJ-1)+ is fixed; wk = (mA-1)+ (tP-mA-1)+ and +// m = min(wc, wk) vary with mA = p + r. Mirrors R `.ExpectedTrit`. +static void expected_trit_by_M(int nI, int nJ, int tc, + std::vector& eM, + std::vector& eMAwk, + std::vector& eMAwc) { + const int nOther = tc - nI - nJ; + const int tP = nI + nJ; + const double wc = posm1(nI) * posm1(nJ); + eM.assign(tc + 1, 0.0); + eMAwk.assign(tc + 1, 0.0); + eMAwc.assign(tc + 1, 0.0); + // Hoist the log-choose terms out of the M / p / r loops (the dominant cost). + std::vector lcI(nI + 1), lcJ(nJ + 1), lcO(nOther + 1), lcT(tc + 1); + for (int p = 0; p <= nI; ++p) lcI[p] = R::lchoose(nI, p); + for (int r = 0; r <= nJ; ++r) lcJ[r] = R::lchoose(nJ, r); + for (int o = 0; o <= nOther; ++o) lcO[o] = R::lchoose(nOther, o); + for (int M = 0; M <= tc; ++M) lcT[M] = R::lchoose(tc, M); + for (int M = 0; M <= tc; ++M) { + const double lt = lcT[M]; + double aM = 0.0, aWk = 0.0, aWc = 0.0; + const int pMax = std::min(nI, M); + for (int p = 0; p <= pMax; ++p) { + const int q = nI - p; + const double lp = lcI[p]; + const int rMax = std::min(nJ, M - p); + for (int r = 0; r <= rMax; ++r) { + const int oA = M - p - r; + if (oA < 0 || oA > nOther) continue; + const double prob = std::exp(lp + lcJ[r] + lcO[oA] - lt); + if (prob <= 0.0) continue; + const int s = nJ - r; + const int mA = p + r; + const double A = posm1(p) * posm1(s) + posm1(q) * posm1(r); + const double wk = posm1(mA) * posm1(tP - mA); + const double m = std::min(wc, wk); + aM += prob * m; + if (wk > 0.0) aWk += prob * m * A / wk; + if (wc > 0.0) aWc += prob * m * A / wc; + } + } + eM[M] = aM; + eMAwk[M] = aWk; + eMAwc[M] = aWc; + } +} + +// [[Rcpp::export]] +List trit_expect(const LogicalMatrix splits, const IntegerMatrix characters) { + const int n_splits = splits.ncol(); + const int n_chars = characters.ncol(); + NumericMatrix numEdge(n_splits, n_chars); + NumericMatrix numChar(n_splits, n_chars); + NumericMatrix denM(n_splits, n_chars); + + for (int c = 0; c < n_chars; ++c) { + const CharCounts cc = char_counts(splits, characters, c); + const int nStates = int(cc.cnt.size()); + if (nStates < 2) continue; + std::vector eM, eWk, eWc; + for (int a = 0; a < nStates - 1; ++a) { + for (int b = a + 1; b < nStates; ++b) { + expected_trit_by_M(cc.cnt[a], cc.cnt[b], cc.tc, eM, eWk, eWc); + for (int s = 0; s < n_splits; ++s) { + const int M = cc.mSideA[s]; + denM(s, c) += eM[M]; + numEdge(s, c) += eWk[M]; + numChar(s, c) += eWc[M]; + } + } + } + } + return List::create(_["numEdge"] = numEdge, _["numChar"] = numChar, + _["denM"] = denM); +} From df97c05d4d5c0a3ea29f84f226f9f6697943fa69 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 21 Aug 2026 11:25:39 +0100 Subject: [PATCH 06/16] refactor(concordance)!: name the chance correction and the NRQS currency `normalize` becomes `chanceCorrect` in `ClusteringConcordance()`, `QuartetConcordance()` and `ConcordanceTable()`. The argument moves the zero point to the value expected under a fixed-marginal null; the scaling to a maximum of 1 was never governed by it, so the old name described something the argument does not do. `unit = "trit"` becomes `unit = "nrqs"`, after the quantity it counts: non-redundant quartet statements. Of the quartets a split of sizes (k, t - k) resolves, (k - 1)(t - k - 1) are non-redundant under the Nelson-Ladiges entailment. `QuartetConcordance()` chance-corrects by default, as `ClusteringConcordance()` already does, so values may be negative where agreement falls below chance expectation. `chanceCorrect = FALSE` gives the uncorrected measure. Documentation gains the caveat that the NRQS ceiling of 1 bounds each split individually, so it holds for `return = "edge"` but not for `return = "char"`, which averages over every split in the tree. # Conflicts: # NEWS.md --- NEWS.md | 31 +++ R/Concordance.R | 250 +++++++++++------- R/RcppExports.R | 6 +- dev/benchmarks/frac-quart/chance_baseline.R | 2 +- dev/benchmarks/frac-quart/cpp_expect_parity.R | 4 +- dev/benchmarks/frac-quart/validate_trit.R | 22 +- man/ConcordanceTable.Rd | 4 +- man/SiteConcordance.Rd | 107 +++++--- man/expected_mi.Rd | 2 +- src/RcppExports.cpp | 8 +- src/TreeSearch-init.c | 4 +- src/concordance_expect.cpp | 10 +- src/expected_mi.cpp | 2 +- tests/testthat/test-Concordance.R | 175 +++++++----- vignettes/tree-search.Rmd | 6 +- 15 files changed, 402 insertions(+), 231 deletions(-) diff --git a/NEWS.md b/NEWS.md index 463268405..bb13b61f5 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,5 +1,36 @@ # To integrate into 2.0.0 notes +- `QuartetConcordance()` gains `unit`, selecting the currency in which quartets + are counted. `unit = "nrqs"` (the new default) counts only non-redundant + quartet statements: of the quartets resolved by a split of sizes + (_k_, _t_ − _k_), just (_k_ − 1)(_t_ − _k_ − 1) are non-redundant under the + Nelson–Ladiges entailment, so a character now scores full marks only where it + displays the split, rather than wherever it agrees with a large combinatorial + volume of quartets. Pass `unit = "quartet"` to recover the raw count. + + Note that "full marks where it displays the split" bounds each split + *individually*, and so applies to `return = "edge"`. It does **not** apply to + `return = "char"`, which averages over every split in the tree: a character + agrees exactly with at most one split and is merely compatible with -- silent + about -- the rest, and coverage scores silence as a failure to cover. A + character identical to one of the tree's own splits therefore scores well + below 1 (0.49 for a 22|26 split of a 48-leaf tree, 0.25 for a 2|46 split), and + the attainable maximum depends on the tree. Per-character `unit = "nrqs"` + values rank characters within a tree but must not be read against a ceiling of + 1 or compared across trees; use `unit = "quartet"` where a per-character index + bounded at 1 is wanted. Documentation only -- no numeric behaviour changed. + +- The `normalize` argument of `ClusteringConcordance()`, `ConcordanceTable()` + and `QuartetConcordance()` is renamed `chanceCorrect`, which is what it + controls: the zero point is moved to the value expected under a fixed-marginal + null. (The scaling to a maximum of 1 was never governed by this argument.) + +- `QuartetConcordance()` now applies that chance correction by default + (`chanceCorrect = TRUE`), matching `ClusteringConcordance()`. Values may + therefore be negative, denoting agreement below chance expectation; pass + `chanceCorrect = FALSE` for the uncorrected measure reported by earlier + versions. + - Fixed: `AdditionTree(constraint = )` silently returned a constraint-violating tree for around one addition order in eleven. Taxa are added to a tree seeded from the first three of them, which is built before the constraint is diff --git a/R/Concordance.R b/R/Concordance.R index 51a484d22..90a607e5d 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -60,18 +60,18 @@ NULL #' so that 1 corresponds to the maximum possible mutual information for each #' split–character pair (`hBest`). #' -#' The `normalize` argument specifies how the zero point is defined. +#' The `chanceCorrect` argument specifies how the zero point is defined. #' -#' - If `normalize = FALSE`, zero corresponds to *zero* MI, without correcting -#' for the positive bias that arises because MI is rarely exactly zero in -#' finite samples. +#' - If `chanceCorrect = FALSE`, zero corresponds to *zero* MI, without +#' correcting for the positive bias that arises because MI is rarely exactly +#' zero in finite samples. #' -#' - If `normalize = TRUE`, the expected MI is computed using an analytical +#' - If `chanceCorrect = TRUE`, the expected MI is computed using an analytical #' approximation based on the distribution of character tokens. This is fast #' and generally accurate for large trees (~200+ taxa), but does not account #' for correlation between splits. #' -#' - If `normalize` is a positive integer `n`, the expected MI is estimated +#' - If `chanceCorrect` is a positive integer `n`, the expected MI is estimated #' empirically by fitting each character to `n` uniformly random trees and #' averaging the resulting MI values. This Monte Carlo approach provides a #' more accurate baseline for small trees, for which the analytical @@ -87,7 +87,7 @@ NULL #' #' Matching is case‑insensitive and partial. #' -#' @param normalize Controls the zero point of the concordance scale by +#' @param chanceCorrect Controls the zero point of the concordance scale by #' subtracting the value expected under a chance (fixed-marginal) null, in which #' each character's tokens are reassigned at random across the leaves while its #' state frequencies and the split sizes are held fixed. @@ -130,7 +130,7 @@ NULL #' corresponds to a split (an edge of the tree) and gives the normalized mutual #' information between that split and the character data, averaged across all #' characters. -#' When `normalize = TRUE` (default), values are scaled relative to random +#' When `chanceCorrect = TRUE` (default), values are scaled relative to random #' expectation; when `FALSE`, raw mutual information normalized by `hBest` is #' returned. #' @@ -172,7 +172,7 @@ ClusteringConcordance <- function( tree, dataset, return = "edge", - normalize = TRUE + chanceCorrect = TRUE ) { # Check inputs if (is.null(dataset)) { @@ -257,7 +257,7 @@ ClusteringConcordance <- function( hh["hSplit", , , drop = FALSE] - hh["hJoint", , , drop = FALSE], NULL) miRand <- `rownames<-`(hh["miRand", , , drop = FALSE], NULL) - norm <- if (isFALSE(normalize)) { + norm <- if (isFALSE(chanceCorrect)) { ifelse(hBest == 0, NA, mi / hBest) } else { ifelse(hBest == 0, NA, .Rezero(mi / hBest, miRand / hBest)) @@ -271,8 +271,8 @@ ClusteringConcordance <- function( charMax <- vapply(charSplits, ClusteringEntropy, double(1))[ attr(dataset, "index")] charInfo <- MutualClusteringInfo(tree, charSplits)[at[["index"]]] - if (is.numeric(normalize)) { - rTrees <- replicate(normalize, RandomTree(tree), simplify = FALSE) + if (is.numeric(chanceCorrect)) { + rTrees <- replicate(chanceCorrect, RandomTree(tree), simplify = FALSE) # Score each random tree against `charSplits` separately: characters with # ambiguous tokens yield splits over different tip subsets, and the # vectorised `MutualClusteringInfo(, )` @@ -285,12 +285,12 @@ ClusteringConcordance <- function( double(length(charSplits)) ))[, attr(dataset, "index"), drop = FALSE] randMean <- colMeans(randInfo) - var <- rowSums((t(randInfo) - randMean) ^ 2) / (normalize - 1) - mcse <- sqrt(var / normalize) + var <- rowSums((t(randInfo) - randMean) ^ 2) / (chanceCorrect - 1) + mcse <- sqrt(var / chanceCorrect) randTreeInfo <- rowSums(randInfo) randTreeMean <- mean(randTreeInfo) treeVar <- var(randTreeInfo) - mcseTree <- sqrt(treeVar / normalize) + mcseTree <- sqrt(treeVar / chanceCorrect) } } @@ -309,7 +309,7 @@ ClusteringConcordance <- function( best <- rowSums(hBest[1, , , drop = FALSE], dims = 2) ifelse(!is.na(best) & best == 0, NA_real_, - if (isTRUE(normalize)) { + if (isTRUE(chanceCorrect)) { .Rezero( rowSums(mi[1, , , drop = FALSE], dims = 2) / best, rowSums(miRand[1, , , drop = FALSE], dims = 2) / best @@ -322,27 +322,28 @@ ClusteringConcordance <- function( # one <- hh["hChar", 1, , drop = TRUE] # All rows equal one <- charMax - zero <- if (isFALSE(normalize)) { + zero <- if (isFALSE(chanceCorrect)) { 0 - } else if (isTRUE(normalize)) { + } else if (isTRUE(chanceCorrect)) { apply(hh["miRand", , ], 2, max) } else { randMean } ret <- (charInfo - zero) / (one - zero) - if (is.numeric(normalize)) { + if (is.numeric(chanceCorrect)) { mcseInfo <- ((one - charInfo) / (one - zero) ^ 2) * mcse mcseInfo[mcseInfo < sqrt(.Machine$double.eps)] <- 0 structure(ret, hMax = charMax, mcse = mcseInfo) } else { # The characterwise return is deliberately NOT random-expectation - # normalized for logical `normalize`: `charInfo` is + # normalized for logical `chanceCorrect`: `charInfo` is # MutualClusteringInfo() against the whole tree, whereas the # analytic `zero` baseline above is per-single-split expected MI, so # subtracting it would mix incompatible quantities (and the # entropy-weighted variant was abandoned -- see the note below the - # @return docs). Only the Monte-Carlo path (numeric `normalize`) - # offers a same-scale empirical baseline. So return charInfo scaled + # @return docs). Only the Monte-Carlo path (numeric + # `chanceCorrect`) offers a same-scale empirical baseline, so we + # return charInfo scaled # by its maximum (hBest-like), as shipped since the original # implementation (#205). structure(charInfo / charMax, hMax = charMax) @@ -358,16 +359,16 @@ ClusteringConcordance <- function( idx <- cbind(1, bestMatch, seq_along(bestMatch)) return(weighted.mean(norm[idx], hBest[idx])) - if (isFALSE(normalize)) { + if (isFALSE(chanceCorrect)) { } else { one <- sum(hh["hChar", 1, ]) - zero <- if (isTRUE(normalize)) { + zero <- if (isTRUE(chanceCorrect)) { sum(apply(hh["miRand", , ], 2, max)) } else { randTreeMean } ret <- (sum(charInfo) - zero) / (one - zero) - if (is.numeric(normalize)) { + if (is.numeric(chanceCorrect)) { mcseInfo <- ((one - sum(charInfo)) / (one - zero) ^ 2) * mcseTree mcseInfo[mcseInfo < sqrt(.Machine$double.eps)] <- 0 structure(ret, mcse = mcseInfo) @@ -542,11 +543,11 @@ QALegend <- function(where = c(0.1, 0.3, 0.1, 0.3), n = 5, Col = QACol, #' @export ConcordanceTable <- function(tree, dataset, Col = QACol, largeClade = 0, xlab = "Edge", ylab = "Character", - normalize = TRUE, plot = TRUE, + chanceCorrect = TRUE, plot = TRUE, marginSize = 0L, paintSize = 0L, palette = "default", ...) { cc <- ClusteringConcordance(tree, dataset, return = "all", - normalize = normalize) + chanceCorrect = chanceCorrect) nodes <- seq_len(dim(cc)[[2]]) info <- cc["hBest", , ] * cc["n", , ] amount <- info / max(info, na.rm = TRUE) @@ -743,9 +744,9 @@ MutualClusteringConcordance <- function(tree, dataset) { #' @rdname SiteConcordance #' @details -#' `QuartetConcordance()` is the proportion of quartets (sets of four leaves) -#' that are decisive for a split which are also concordant with it -#' For example, a quartet with the characters `0 0 0 1` is not decisive, as +#' `QuartetConcordance()` measures the agreement between each character and each +#' split in the currency of quartets (sets of four leaves). +#' A quartet with the characters `0 0 0 1` is not decisive, as #' all relationships between those leaves are equally parsimonious. #' But a quartet with characters `0 0 1 1` is decisive, and is concordant #' with any tree that groups the first two leaves together to the exclusion @@ -755,12 +756,21 @@ MutualClusteringConcordance <- function(tree, dataset) { #' quartets that are decisive for a branch. #' Doing so circumvents the criticisms of \insertCite{Goloboff2024;textual}{TreeSearch}. #' -#' By default, the reported value weights each site by the number of quartets -#' it is decisive for. This value can be interpreted as the proportion of -#' all decisive quartets that are concordant with a split. +#' The quartets that a split resolves are not logically independent: because +#' \eqn{(ab, cd)} and \eqn{(ab, ce)} jointly entail \eqn{(ab, de)} +#' \insertCite{Nelson1992}{TreeSearch}, a split of sizes \eqn{(k, t - k)} +#' resolves \eqn{\binom{k}{2}\binom{t - k}{2}} quartets, of which only +#' \eqn{(k - 1)(t - k - 1)} are non-redundant. +#' `unit` selects which of the two is counted. By default, agreement is +#' measured against the non-redundant content, so that a character scores full +#' marks only where it *displays* the split, rather than merely agreeing with a +#' large combinatorial volume of its quartets. +#' +#' By default, the reported value weights each site by the quartet content it +#' shares with the split. #' If `weight = FALSE`, the reported value is the mean of the concordance -#' value for each site. -#' Consider a split associated with two sites: +#' value for each site. +#' Consider a split associated with two sites (counting in `unit = "quartet"`): #' one that is concordant with 25% of 96 decisive quartets, and #' a second that is concordant with 75% of 4 decisive quartets. #' If `weight = TRUE`, the split concordance will be 24 + 3 / 96 + 4 = 27%. @@ -779,30 +789,50 @@ MutualClusteringConcordance <- function(tree, dataset) { #' `QuartetConcordance(return = "char")` returns a numeric vector giving the #' concordance index calculated at each site, averaged across all splits. #' -#' With `unit = "trit"` (see below) the same vectors are returned, but scored in -#' redundancy-corrected currency: values are typically lower, and reach 1 only -#' where a split is *displayed* by the character rather than merely concordant -#' with many of its quartets. +#' With `unit = "nrqs"`, note that this average is taken over *every* split in +#' `tree`. A character can agree exactly with at most one split, and is merely +#' compatible with -- that is, silent about -- the rest; because agreement is +#' measured against each split's non-redundant content, silence counts as a +#' failure to cover. Per-character values are therefore much smaller than +#' per-edge ones and **do not reach 1 even for a character identical to one of +#' the tree's own splits** (0.49 for a 22|26 split of a 48-leaf tree, 0.25 for +#' a 2|46 split). The attainable maximum depends on the tree, so these values +#' rank characters against one another but should not be read against an +#' absolute ceiling of 1, nor compared across trees. Where a per-character +#' index bounded at 1 is wanted, use `unit = "quartet"`, whose per-character +#' path attains 1 for a character identical to any split of the tree. +#' +#' With `unit = "quartet"` the same vectors are returned, but scored in the raw +#' quartet currency: values are typically higher, as a character is credited +#' with the full combinatorial volume of quartets that it resolves. #' #' @param weight Logical specifying whether to weight sites according to the -#' number of quartets they are decisive for. +#' quartet content that they share with each split. #' @param unit Character specifying the currency in which quartets are counted: -#' - `"quartet"` (default): each resolved quartet counts once, so a character -#' is credited with the full combinatorial volume of quartets it resolves; -#' - `"trit"`: quartets are counted as the *independent* information they carry -#' \insertCite{Nelson1992}{TreeSearch}. Because the quartets resolved by a -#' split of sizes \eqn{(k, t - k)} are logically redundant -- -#' \eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} -- only -#' \eqn{(k - 1)(t - k - 1)} of the \eqn{\binom{k}{2}\binom{t - k}{2}} -#' resolved quartets are independent. Concordance is then measured against -#' this reduced (redundancy-corrected) content, so that only a character -#' whose split is *identical* to the tree split scores full marks; nested -#' (compatible) characters receive genuine partial support, and crossing -#' (incompatible) characters score lower still. A multistate character is -#' scored in the same currency: trits are counted independently within each -#' pair of states (sizes \eqn{n_i}, \eqn{n_j} give a per-state-pair weight of -#' \eqn{4 / (n_i n_j)}) and summed across pairs, so multistate and binary -#' characters remain directly comparable. +#' - `"nrqs"` (default): only non-redundant quartet statements (NRQS) are +#' counted. Of the \eqn{\binom{k}{2}\binom{t - k}{2}} quartets resolved by +#' a split of sizes \eqn{(k, t - k)}, just \eqn{(k - 1)(t - k - 1)} are +#' non-redundant under the entailment +#' \eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} +#' \insertCite{Nelson1992}{TreeSearch}. Measuring agreement against this +#' reduced content means that, *at a given split*, only a character whose +#' own split is identical scores full marks; nested (compatible) characters +#' receive genuine partial support, and crossing (incompatible) characters +#' score lower still. This bound applies to each split individually, and so +#' to `return = "edge"`; it does not carry over to `return = "char"`, which +#' averages across all splits -- see Value, above. +#' A multistate character is scored in the same currency: +#' NRQS are counted within each pair of states (sizes \eqn{n_i}, \eqn{n_j} +#' give a per-state-pair weight of \eqn{4 / (n_i n_j)}) and summed across +#' pairs, so multistate and binary characters remain directly comparable. +#' Note that non-redundancy is a logical property: NRQS are not mutually +#' independent in a statistical sense, so their count measures quartet +#' content rather than information in bits. +#' - `"quartet"`: each resolved quartet counts once, so a character is +#' credited with the full combinatorial volume of quartets it resolves. +#' This recovers the measure reported by earlier versions of this package, +#' the direct analogue of the site concordance factor +#' \insertCite{Minh2020}{TreeSearch}. #' @references \insertAllCited{} #' @importFrom ape keep.tip #' @importFrom cli cli_progress_bar cli_progress_update @@ -814,8 +844,8 @@ QuartetConcordance <- function( dataset = NULL, weight = TRUE, return = "edge", - unit = c("quartet", "trit"), - normalize = FALSE + unit = c("nrqs", "quartet"), + chanceCorrect = TRUE ) { if (is.null(dataset)) { warning("Cannot calculate concordance without `dataset`.") @@ -825,11 +855,11 @@ QuartetConcordance <- function( stop("`dataset` must be a phyDat object.") } unit <- match.arg(unit) - if (!isFALSE(normalize)) { - if (!isTRUE(normalize)) { - if (!is.numeric(normalize) || length(normalize) != 1L || - is.na(normalize) || normalize < 1) { - stop("`normalize` must be FALSE, TRUE, or a positive integer.") + if (!isFALSE(chanceCorrect)) { + if (!isTRUE(chanceCorrect)) { + if (!is.numeric(chanceCorrect) || length(chanceCorrect) != 1L || + is.na(chanceCorrect) || chanceCorrect < 1) { + stop("`chanceCorrect` must be FALSE, TRUE, or a positive integer.") } } } @@ -868,10 +898,10 @@ QuartetConcordance <- function( return <- options[[pmatch(tolower(trimws(return)), options, nomatch = length(options))]] - if (unit == "trit") { + if (unit == "nrqs") { # Return: - return(.TritConcordance(logiSplits, charInt, weight, return, splits, - normalize)) + return(.NrqsConcordance(logiSplits, charInt, weight, return, splits, + chanceCorrect)) } raw_counts <- quartet_concordance(logiSplits, charInt) @@ -881,16 +911,16 @@ QuartetConcordance <- function( # Chance correction (option): re-zero the observed concordant/decisive ratio # against the ratio expected under the same fixed-marginal null used for the - # trit currency. Only `conc` and `dec` vary under the null (the split sizes + # NRQS currency. Only `conc` and `dec` vary under the null (the split sizes # and state counts are fixed), so we need E[conc] and E[dec] per (split, char), - # computed exactly from the hypergeometric pmf (no floors here, unlike trits) + # computed exactly from the hypergeometric pmf (no floors here, unlike NRQS) # or by Monte-Carlo tip-shuffle, then re-zero the pooled ratio. - doNorm <- !isFALSE(normalize) + doNorm <- !isFALSE(chanceCorrect) if (doNorm) { - base <- if (isTRUE(normalize)) { + base <- if (isTRUE(chanceCorrect)) { .QuartetExpect(charInt, logiSplits) } else { - .QuartetMC(charInt, logiSplits, normalize) + .QuartetMC(charInt, logiSplits, chanceCorrect) } eNum <- base[["concordant"]] eDen <- base[["decisive"]] @@ -1009,33 +1039,57 @@ QuartetConcordance <- function( list(concordant = accConc / nRelabel, decisive = accDec / nRelabel) } -# Nelson-Ladiges fractional ("trit") currency for QuartetConcordance(). +# Nelson-Ladiges fractional ("nrqs") currency for QuartetConcordance(). # # The quartets a bipartition of sizes (k, t - k) resolves are the 4-cycles of # the complete bipartite graph K_{k, t-k}; the Nelson-Ladiges entailment # (ab,cd) + (ab,ce) -> (ab,de) is GF(2) cycle addition, so only the cyclomatic -# number (k - 1)(t - k - 1) of them are independent ("trits"). +# number (k - 1)(t - k - 1) of them are non-redundant (the NRQS). # # A multistate character is the disjoint union of its state-pairs: a quartet is # decisive only when two taxa share one state and two share another, so every # decisive quartet lives in the K_{n_i, n_j} block between two states. Those # blocks are edge-disjoint and the entailment never crosses them (it forces the -# shared "c, d, e" taxa into a single state), so trits add over state pairs: +# shared "c, d, e" taxa into a single state), so NRQS add over state pairs: # W_char = sum_{i 0` guard. - be <- trit_expect(logiSplits, charInt) + be <- nrqs_expect(logiSplits, charInt) baseNumEdge <- be[["numEdge"]] baseNumChar <- be[["numChar"]] baseDenM <- be[["denM"]] @@ -1064,13 +1118,14 @@ QuartetConcordance <- function( for (ci in seq_len(nChar)) { col <- charInt[, ci] - obs <- .CharTritContrib(col, logiSplits, pos) + obs <- .CharNrqsContrib(col, logiSplits, pos) numEdge[, ci] <- obs[["numEdge"]] numChar[, ci] <- obs[["numChar"]] denM[, ci] <- obs[["denM"]] wcTot[ci] <- obs[["wcTot"]] - if (doNorm && !isTRUE(normalize) && obs[["wcTot"]] > 0) { - base <- .CharTritMC(col, logiSplits, pos, normalize) # Monte-Carlo shuffle + if (doNorm && !isTRUE(chanceCorrect) && obs[["wcTot"]] > 0) { + # Monte-Carlo shuffle + base <- .CharNrqsMC(col, logiSplits, pos, chanceCorrect) baseNumEdge[, ci] <- base[["numEdge"]] baseNumChar[, ci] <- base[["numChar"]] baseDenM[, ci] <- base[["denM"]] @@ -1091,7 +1146,7 @@ QuartetConcordance <- function( } } else { # Mean per-site quality over informative characters; uninformative - # characters carry no trits and are dropped, as in the quartet path. + # characters carry no NRQS and are dropped, as in the quartet path. sEdge <- ifelse(denM > 0, numEdge / denM, NA_real_) if (doNorm) { sBase <- ifelse(baseDenM > 0, baseNumEdge / baseDenM, NA_real_) @@ -1138,10 +1193,10 @@ QuartetConcordance <- function( } } -# Observed per-character trit contributions, summed over the character's +# Observed per-character NRQS contributions, summed over the character's # state-pairs, returned as per-split vectors. Shared by the observed pass and # the Monte-Carlo baseline so the two use byte-identical arithmetic. -.CharTritContrib <- function(col, logiSplits, pos) { +.CharNrqsContrib <- function(col, logiSplits, pos) { nSplit <- ncol(logiSplits) numEdge <- numChar <- denM <- numeric(nSplit) wcTot <- 0 @@ -1161,7 +1216,7 @@ QuartetConcordance <- function( nJ <- aJ + bJ # n_j mA <- aI + aJ # taxa of this pair on side A tP <- nI + nJ # taxa scored in this pair - # Concordant trits; the (x - 1)_+ floors stop self-agreement exceeding 1. + # Concordant NRQS; the (x - 1)_+ floors stop self-agreement exceeding 1. aij <- pos(aI - 1) * pos(bJ - 1) + pos(bI - 1) * pos(aJ - 1) wc <- pos(nI - 1) * pos(nJ - 1) # pair's character content wk <- pos(mA - 1) * pos(tP - mA - 1) # pair's split content @@ -1176,16 +1231,17 @@ QuartetConcordance <- function( list(numEdge = numEdge, numChar = numChar, denM = denM, wcTot = wcTot) } -# The exact trit expectation E[m], E[m*A/wk], E[m*A/wc] is computed in C++ -# (`trit_expect` in src/concordance_expect.cpp) for all characters at once, and -# consumed directly by `.TritConcordance`; there is no per-character R exact -# helper. The Monte-Carlo baseline below is the opt-in `normalize = ` path. +# The exact NRQS expectation E[m], E[m*A/wk], E[m*A/wc] is computed in C++ +# (`nrqs_expect` in src/concordance_expect.cpp) for all characters at once, and +# consumed directly by `.NrqsConcordance`; there is no per-character R exact +# helper. The Monte-Carlo baseline below is the opt-in +# `chanceCorrect = ` path. -# Monte-Carlo baseline: average `.CharTritContrib` over `nRelabel` random +# Monte-Carlo baseline: average `.CharNrqsContrib` over `nRelabel` random # reassignments of the character's tokens across its scored leaves (the split, # the scored set and the state counts are held fixed). Averaging the pools # (rather than per-split ratios) mirrors the exact ratio-of-expectations path. -.CharTritMC <- function(col, logiSplits, pos, nRelabel) { +.CharNrqsMC <- function(col, logiSplits, pos, nRelabel) { nSplit <- ncol(logiSplits) scored <- !is.na(col) tokens <- col[scored] @@ -1193,7 +1249,7 @@ QuartetConcordance <- function( for (i in seq_len(nRelabel)) { shuffled <- col shuffled[scored] <- sample(tokens) - cc <- .CharTritContrib(shuffled, logiSplits, pos) + cc <- .CharNrqsContrib(shuffled, logiSplits, pos) accEdge <- accEdge + cc[["numEdge"]] accChar <- accChar + cc[["numChar"]] accDen <- accDen + cc[["denM"]] diff --git a/R/RcppExports.R b/R/RcppExports.R index 4d29eae1d..da8f9b5cb 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -24,8 +24,8 @@ quartet_expect <- function(splits, characters) { .Call(`_TreeSearch_quartet_expect`, splits, characters) } -trit_expect <- function(splits, characters) { - .Call(`_TreeSearch_trit_expect`, splits, characters) +nrqs_expect <- function(splits, characters) { + .Call(`_TreeSearch_nrqs_expect`, splits, characters) } #' Expected mutual information between two partitions @@ -35,7 +35,7 @@ trit_expect <- function(splits, characters) { #' block sizes (marginals) of each partition are fixed but the items are #' associated at random. Subtracting this baseline from an observed mutual #' information yields a chance-corrected ("adjusted") mutual information, as -#' applied by [`SiteConcordance`]`(normalize = TRUE)`; it is most material for +#' applied by [`SiteConcordance`]`(chanceCorrect = TRUE)`; it is most material for #' small or unbalanced partitions, where raw mutual information is appreciably #' inflated by chance agreement. #' diff --git a/dev/benchmarks/frac-quart/chance_baseline.R b/dev/benchmarks/frac-quart/chance_baseline.R index 93b9d9690..38187da82 100644 --- a/dev/benchmarks/frac-quart/chance_baseline.R +++ b/dev/benchmarks/frac-quart/chance_baseline.R @@ -1,5 +1,5 @@ # Validate the exact hypergeometric chance baseline used by -# QuartetConcordance(unit = "trit", normalize = TRUE) against a Monte-Carlo +# QuartetConcordance(unit = "nrqs", chanceCorrect = TRUE) against a Monte-Carlo # tip-shuffle oracle, at the level of a single (split, character, state-pair). # # Under the fixed-marginal null the character's tokens are reassigned across the diff --git a/dev/benchmarks/frac-quart/cpp_expect_parity.R b/dev/benchmarks/frac-quart/cpp_expect_parity.R index 0e8afba7f..10972edd3 100644 --- a/dev/benchmarks/frac-quart/cpp_expect_parity.R +++ b/dev/benchmarks/frac-quart/cpp_expect_parity.R @@ -1,5 +1,5 @@ # Parity guard for the C++ chance-baseline kernels `quartet_expect` and -# `trit_expect` (src/concordance_expect.cpp), which compute the exact expected +# `nrqs_expect` (src/concordance_expect.cpp), which compute the exact expected # concordant/decisive quartet counts (raw) and the expected trit pools # (E[m], E[m*A/wk], E[m*A/wc]) under the fixed-marginal (hypergeometric) null. # @@ -95,7 +95,7 @@ for (rep in 1:10) { io <- build_inputs(tree, MatrixToPhyDat(m)) ls <- io$logiSplits; ci <- io$charInt R <- ref_expect(ci, ls) - cq <- quartet_expect(ls, ci); ct <- trit_expect(ls, ci) + cq <- quartet_expect(ls, ci); ct <- nrqs_expect(ls, ci) d <- max(abs(cq$concordant - R$eConc), abs(cq$decisive - R$eDec), abs(ct$numEdge - R$eWk), abs(ct$numChar - R$eWc), abs(ct$denM - R$eM)) diff --git a/dev/benchmarks/frac-quart/validate_trit.R b/dev/benchmarks/frac-quart/validate_trit.R index b0f9142e9..b29e053c6 100644 --- a/dev/benchmarks/frac-quart/validate_trit.R +++ b/dev/benchmarks/frac-quart/validate_trit.R @@ -1,4 +1,4 @@ -# Validation of the unit = "trit" path in QuartetConcordance(). +# Validation of the unit = "nrqs" path in QuartetConcordance(). # # Ordering follows the advisor's priority: # 1. Cell extraction cross-checked against the verified C++ kernel (conc/dec). @@ -161,7 +161,7 @@ data("congreveLamsdellMatrices", package = "TreeSearch") binDat <- congreveLamsdellMatrices[[1]] refTree <- TreeSearch::referenceTree for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) { - pkg <- QuartetConcordance(refTree, binDat, weight = w, return = r, unit = "trit") + pkg <- QuartetConcordance(refTree, binDat, weight = w, return = r, unit = "nrqs", chanceCorrect = FALSE) rf <- ref_trit(refTree, binDat, weight = w, return = r) ok(sprintf("binary pkg==ref weight=%s return=%s", w, r), same(pkg, rf)) } @@ -174,7 +174,7 @@ msMat <- matrix(c(0, 0, 1, 1, 2, 2, 2, 0, # 3 states msDat <- MatrixToPhyDat(msMat) msTree <- BalancedTree(8) for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) { - pkg <- QuartetConcordance(msTree, msDat, weight = w, return = r, unit = "trit") + pkg <- QuartetConcordance(msTree, msDat, weight = w, return = r, unit = "nrqs", chanceCorrect = FALSE) rf <- ref_trit(msTree, msDat, weight = w, return = r) ok(sprintf("multi pkg==ref weight=%s return=%s", w, r), same(pkg, rf)) } @@ -186,17 +186,17 @@ naMat <- matrix(c(0, 0, 1, 1, "?", "?", 0, 1, dimnames = list(paste0("t", 1:8), NULL)) naDat <- MatrixToPhyDat(naMat) for (w in c(TRUE, FALSE)) for (r in c("edge", "char")) { - pkg <- QuartetConcordance(msTree, naDat, weight = w, return = r, unit = "trit") + pkg <- QuartetConcordance(msTree, naDat, weight = w, return = r, unit = "nrqs", chanceCorrect = FALSE) rf <- ref_trit(msTree, naDat, weight = w, return = r) ok(sprintf("missing pkg==ref weight=%s return=%s", w, r), same(pkg, rf)) } # ---------- 5. `return` parsing mirrors the quartet path ---------- -e1 <- QuartetConcordance(refTree, binDat, return = "edge", unit = "trit") -e2 <- QuartetConcordance(refTree, binDat, return = "default", unit = "trit") -c1 <- QuartetConcordance(refTree, binDat, return = "char", unit = "trit") -c2 <- QuartetConcordance(refTree, binDat, return = "character", unit = "trit") -c3 <- QuartetConcordance(refTree, binDat, return = "site", unit = "trit") +e1 <- QuartetConcordance(refTree, binDat, return = "edge", unit = "nrqs", chanceCorrect = FALSE) +e2 <- QuartetConcordance(refTree, binDat, return = "default", unit = "nrqs", chanceCorrect = FALSE) +c1 <- QuartetConcordance(refTree, binDat, return = "char", unit = "nrqs", chanceCorrect = FALSE) +c2 <- QuartetConcordance(refTree, binDat, return = "character", unit = "nrqs", chanceCorrect = FALSE) +c3 <- QuartetConcordance(refTree, binDat, return = "site", unit = "nrqs", chanceCorrect = FALSE) ok("return edge == default", same(e1, e2)) ok("return char == character == site", same(c1, c2) && same(c1, c3)) ok("edge length == n splits, char length == n chars", @@ -204,8 +204,8 @@ ok("edge length == n splits, char length == n chars", length(c1) == sum(attr(binDat, "weight"))) # ---------- 6. congreveLamsdell quartet -> trit movement + census ---------- -qEdge <- QuartetConcordance(refTree, binDat, unit = "quartet") -tEdge <- QuartetConcordance(refTree, binDat, unit = "trit") +qEdge <- QuartetConcordance(refTree, binDat, unit = "quartet", chanceCorrect = FALSE) +tEdge <- QuartetConcordance(refTree, binDat, unit = "nrqs", chanceCorrect = FALSE) cat("\nEdge concordance, quartet vs trit (congreveLamsdell[[1]]):\n") print(round(rbind(quartet = qEdge, trit = tEdge, delta = tEdge - qEdge), 3)) cat(sprintf("\nmean quartet = %.3f mean trit = %.3f (trit is stricter)\n", diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index b420da41c..c34e96780 100644 --- a/man/ConcordanceTable.Rd +++ b/man/ConcordanceTable.Rd @@ -11,7 +11,7 @@ ConcordanceTable( largeClade = 0, xlab = "Edge", ylab = "Character", - normalize = TRUE, + chanceCorrect = TRUE, plot = TRUE, marginSize = 0L, paintSize = 0L, @@ -36,7 +36,7 @@ at edges whose descendants are both contain more than \code{largeClade} leaves.} \item{ylab}{Character giving a label for the y axis.} -\item{normalize}{Controls the zero point of the concordance scale by +\item{chanceCorrect}{Controls the zero point of the concordance scale by subtracting the value expected under a chance (fixed-marginal) null, in which each character's tokens are reassigned at random across the leaves while its state frequencies and the split sizes are held fixed. diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index f8b3c6cc0..c5f46b630 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -9,7 +9,7 @@ \alias{SharedPhylogeneticConcordance} \title{Concordance factors} \usage{ -ClusteringConcordance(tree, dataset, return = "edge", normalize = TRUE) +ClusteringConcordance(tree, dataset, return = "edge", chanceCorrect = TRUE) MutualClusteringConcordance(tree, dataset) @@ -18,8 +18,8 @@ QuartetConcordance( dataset = NULL, weight = TRUE, return = "edge", - unit = c("quartet", "trit"), - normalize = FALSE + unit = c("nrqs", "quartet"), + chanceCorrect = TRUE ) PhylogeneticConcordance(tree, dataset) @@ -44,7 +44,7 @@ split–character pair. Matching is case‑insensitive and partial.} -\item{normalize}{Controls the zero point of the concordance scale by +\item{chanceCorrect}{Controls the zero point of the concordance scale by subtracting the value expected under a chance (fixed-marginal) null, in which each character's tokens are reassigned at random across the leaves while its state frequencies and the split sizes are held fixed. @@ -71,26 +71,34 @@ unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{weight}{Logical specifying whether to weight sites according to the -number of quartets they are decisive for.} +quartet content that they share with each split.} \item{unit}{Character specifying the currency in which quartets are counted: \itemize{ -\item \code{"quartet"} (default): each resolved quartet counts once, so a character -is credited with the full combinatorial volume of quartets it resolves; -\item \code{"trit"}: quartets are counted as the \emph{independent} information they carry -\insertCite{Nelson1992}{TreeSearch}. Because the quartets resolved by a -split of sizes \eqn{(k, t - k)} are logically redundant -- -\eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} -- only -\eqn{(k - 1)(t - k - 1)} of the \eqn{\binom{k}{2}\binom{t - k}{2}} -resolved quartets are independent. Concordance is then measured against -this reduced (redundancy-corrected) content, so that only a character -whose split is \emph{identical} to the tree split scores full marks; nested -(compatible) characters receive genuine partial support, and crossing -(incompatible) characters score lower still. A multistate character is -scored in the same currency: trits are counted independently within each -pair of states (sizes \eqn{n_i}, \eqn{n_j} give a per-state-pair weight of -\eqn{4 / (n_i n_j)}) and summed across pairs, so multistate and binary -characters remain directly comparable. +\item \code{"nrqs"} (default): only non-redundant quartet statements (NRQS) are +counted. Of the \eqn{\binom{k}{2}\binom{t - k}{2}} quartets resolved by +a split of sizes \eqn{(k, t - k)}, just \eqn{(k - 1)(t - k - 1)} are +non-redundant under the entailment +\eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} +\insertCite{Nelson1992}{TreeSearch}. Measuring agreement against this +reduced content means that, \emph{at a given split}, only a character whose +own split is identical scores full marks; nested (compatible) characters +receive genuine partial support, and crossing (incompatible) characters +score lower still. This bound applies to each split individually, and so +to \code{return = "edge"}; it does not carry over to \code{return = "char"}, which +averages across all splits -- see Value, above. +A multistate character is scored in the same currency: +NRQS are counted within each pair of states (sizes \eqn{n_i}, \eqn{n_j} +give a per-state-pair weight of \eqn{4 / (n_i n_j)}) and summed across +pairs, so multistate and binary characters remain directly comparable. +Note that non-redundancy is a logical property: NRQS are not mutually +independent in a statistical sense, so their count measures quartet +content rather than information in bits. +\item \code{"quartet"}: each resolved quartet counts once, so a character is +credited with the full combinatorial volume of quartets it resolves. +This recovers the measure reported by earlier versions of this package, +the direct analogue of the site concordance factor +\insertCite{Minh2020}{TreeSearch}. }} } \value{ @@ -112,7 +120,7 @@ expectation. \code{NA} is returned when \code{hBest = 0} (no information potenti corresponds to a split (an edge of the tree) and gives the normalized mutual information between that split and the character data, averaged across all characters. -When \code{normalize = TRUE} (default), values are scaled relative to random +When \code{chanceCorrect = TRUE} (default), values are scaled relative to random expectation; when \code{FALSE}, raw mutual information normalized by \code{hBest} is returned. @@ -140,10 +148,22 @@ each corresponding split in \code{tree}. \code{QuartetConcordance(return = "char")} returns a numeric vector giving the concordance index calculated at each site, averaged across all splits. -With \code{unit = "trit"} (see below) the same vectors are returned, but scored in -redundancy-corrected currency: values are typically lower, and reach 1 only -where a split is \emph{displayed} by the character rather than merely concordant -with many of its quartets. +With \code{unit = "nrqs"}, note that this average is taken over \emph{every} split in +\code{tree}. A character can agree exactly with at most one split, and is merely +compatible with -- that is, silent about -- the rest; because agreement is +measured against each split's non-redundant content, silence counts as a +failure to cover. Per-character values are therefore much smaller than +per-edge ones and \strong{do not reach 1 even for a character identical to one of +the tree's own splits} (0.49 for a 22|26 split of a 48-leaf tree, 0.25 for +a 2|46 split). The attainable maximum depends on the tree, so these values +rank characters against one another but should not be read against an +absolute ceiling of 1, nor compared across trees. Where a per-character +index bounded at 1 is wanted, use \code{unit = "quartet"}, whose per-character +path attains 1 for a character identical to any split of the tree. + +With \code{unit = "quartet"} the same vectors are returned, but scored in the raw +quartet currency: values are typically higher, as a character is credited +with the full combinatorial volume of quartets that it resolves. \code{PhylogeneticConcordance()} returns a numeric vector giving the phylogenetic information of each split in \code{tree}, named according to the @@ -168,16 +188,16 @@ is quantified using mutual information (MI). All reported values are scaled so that 1 corresponds to the maximum possible mutual information for each split–character pair (\code{hBest}). -The \code{normalize} argument specifies how the zero point is defined. +The \code{chanceCorrect} argument specifies how the zero point is defined. \itemize{ -\item If \code{normalize = FALSE}, zero corresponds to \emph{zero} MI, without correcting -for the positive bias that arises because MI is rarely exactly zero in -finite samples. -\item If \code{normalize = TRUE}, the expected MI is computed using an analytical +\item If \code{chanceCorrect = FALSE}, zero corresponds to \emph{zero} MI, without +correcting for the positive bias that arises because MI is rarely exactly +zero in finite samples. +\item If \code{chanceCorrect = TRUE}, the expected MI is computed using an analytical approximation based on the distribution of character tokens. This is fast and generally accurate for large trees (~200+ taxa), but does not account for correlation between splits. -\item If \code{normalize} is a positive integer \code{n}, the expected MI is estimated +\item If \code{chanceCorrect} is a positive integer \code{n}, the expected MI is estimated empirically by fitting each character to \code{n} uniformly random trees and averaging the resulting MI values. This Monte Carlo approach provides a more accurate baseline for small trees, for which the analytical @@ -191,9 +211,9 @@ character‑tree and the supplied phylogeny. High values identify characters whose signal is well represented anywhere in the tree, even if concentrated on a single edge. -\code{QuartetConcordance()} is the proportion of quartets (sets of four leaves) -that are decisive for a split which are also concordant with it -For example, a quartet with the characters \verb{0 0 0 1} is not decisive, as +\code{QuartetConcordance()} measures the agreement between each character and each +split in the currency of quartets (sets of four leaves). +A quartet with the characters \verb{0 0 0 1} is not decisive, as all relationships between those leaves are equally parsimonious. But a quartet with characters \verb{0 0 1 1} is decisive, and is concordant with any tree that groups the first two leaves together to the exclusion @@ -203,12 +223,21 @@ In contrast to the site concordance factor quartets that are decisive for a branch. Doing so circumvents the criticisms of \insertCite{Goloboff2024;textual}{TreeSearch}. -By default, the reported value weights each site by the number of quartets -it is decisive for. This value can be interpreted as the proportion of -all decisive quartets that are concordant with a split. +The quartets that a split resolves are not logically independent: because +\eqn{(ab, cd)} and \eqn{(ab, ce)} jointly entail \eqn{(ab, de)} +\insertCite{Nelson1992}{TreeSearch}, a split of sizes \eqn{(k, t - k)} +resolves \eqn{\binom{k}{2}\binom{t - k}{2}} quartets, of which only +\eqn{(k - 1)(t - k - 1)} are non-redundant. +\code{unit} selects which of the two is counted. By default, agreement is +measured against the non-redundant content, so that a character scores full +marks only where it \emph{displays} the split, rather than merely agreeing with a +large combinatorial volume of its quartets. + +By default, the reported value weights each site by the quartet content it +shares with the split. If \code{weight = FALSE}, the reported value is the mean of the concordance value for each site. -Consider a split associated with two sites: +Consider a split associated with two sites (counting in \code{unit = "quartet"}): one that is concordant with 25\% of 96 decisive quartets, and a second that is concordant with 75\% of 4 decisive quartets. If \code{weight = TRUE}, the split concordance will be 24 + 3 / 96 + 4 = 27\%. diff --git a/man/expected_mi.Rd b/man/expected_mi.Rd index 5002750e4..a87f7abae 100644 --- a/man/expected_mi.Rd +++ b/man/expected_mi.Rd @@ -22,7 +22,7 @@ partitions of the same \code{N} items, under the hypergeometric null in which th block sizes (marginals) of each partition are fixed but the items are associated at random. Subtracting this baseline from an observed mutual information yields a chance-corrected ("adjusted") mutual information, as -applied by \code{\link{SiteConcordance}}\code{(normalize = TRUE)}; it is most material for +applied by \code{\link{SiteConcordance}}\code{(chanceCorrect = TRUE)}; it is most material for small or unbalanced partitions, where raw mutual information is appreciably inflated by chance agreement. } diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index f46299b2e..ff69d5d3f 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -43,15 +43,15 @@ BEGIN_RCPP return rcpp_result_gen; END_RCPP } -// trit_expect -List trit_expect(const LogicalMatrix splits, const IntegerMatrix characters); -RcppExport SEXP _TreeSearch_trit_expect(SEXP splitsSEXP, SEXP charactersSEXP) { +// nrqs_expect +List nrqs_expect(const LogicalMatrix splits, const IntegerMatrix characters); +RcppExport SEXP _TreeSearch_nrqs_expect(SEXP splitsSEXP, SEXP charactersSEXP) { BEGIN_RCPP Rcpp::RObject rcpp_result_gen; Rcpp::RNGScope rcpp_rngScope_gen; Rcpp::traits::input_parameter< const LogicalMatrix >::type splits(splitsSEXP); Rcpp::traits::input_parameter< const IntegerMatrix >::type characters(charactersSEXP); - rcpp_result_gen = Rcpp::wrap(trit_expect(splits, characters)); + rcpp_result_gen = Rcpp::wrap(nrqs_expect(splits, characters)); return rcpp_result_gen; END_RCPP } diff --git a/src/TreeSearch-init.c b/src/TreeSearch-init.c index 890fd4255..c7d16856e 100644 --- a/src/TreeSearch-init.c +++ b/src/TreeSearch-init.c @@ -22,7 +22,7 @@ extern SEXP _TreeSearch_mi_key(SEXP, SEXP); // extern SEXP _TreeSearch_astar_search_r(SEXP, SEXP, SEXP); extern SEXP _TreeSearch_quartet_concordance(SEXP, SEXP); extern SEXP _TreeSearch_quartet_expect(SEXP, SEXP); -extern SEXP _TreeSearch_trit_expect(SEXP, SEXP); +extern SEXP _TreeSearch_nrqs_expect(SEXP, SEXP); extern SEXP _TreeSearch_ts_fitch_score(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _TreeSearch_ts_na_char_steps(SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _TreeSearch_ts_char_steps(SEXP, SEXP, SEXP, SEXP, SEXP); @@ -78,7 +78,7 @@ static const R_CallMethodDef callMethods[] = { // {"_TreeSearch_astar_search_r", (DL_FUNC) &_TreeSearch_astar_search_r, 3}, {"_TreeSearch_quartet_concordance",(DL_FUNC) &_TreeSearch_quartet_concordance, 2}, {"_TreeSearch_quartet_expect", (DL_FUNC) &_TreeSearch_quartet_expect, 2}, - {"_TreeSearch_trit_expect", (DL_FUNC) &_TreeSearch_trit_expect, 2}, + {"_TreeSearch_nrqs_expect", (DL_FUNC) &_TreeSearch_nrqs_expect, 2}, {"_TreeSearch_ts_fitch_score", (DL_FUNC) &_TreeSearch_ts_fitch_score, 12}, {"_TreeSearch_ts_na_char_steps", (DL_FUNC) &_TreeSearch_ts_na_char_steps, 5}, {"_TreeSearch_ts_char_steps", (DL_FUNC) &_TreeSearch_ts_char_steps, 5}, diff --git a/src/concordance_expect.cpp b/src/concordance_expect.cpp index 4faefac8e..b6b201920 100644 --- a/src/concordance_expect.cpp +++ b/src/concordance_expect.cpp @@ -6,7 +6,7 @@ using namespace Rcpp; // Expected concordance under the fixed-marginal (hypergeometric) null used by -// QuartetConcordance(normalize = TRUE). For one character each token is +// QuartetConcordance(chanceCorrect = TRUE). For one character each token is // reassigned at random across the scored leaves, holding the state counts and // the split sizes fixed. For a state-pair (i, j) with counts nI, nJ over t // scored leaves, of which M are on side A of a split, the 2x2 cell vector @@ -137,10 +137,10 @@ List quartet_expect(const LogicalMatrix splits, const IntegerMatrix characters) return List::create(_["concordant"] = eConc, _["decisive"] = eDec); } -// Exact E[m], E[m*A/wk], E[m*A/wc] for one trit state-pair over all side-A +// Exact E[m], E[m*A/wk], E[m*A/wc] for one NRQS state-pair over all side-A // sizes. wc = (nI-1)+ (nJ-1)+ is fixed; wk = (mA-1)+ (tP-mA-1)+ and // m = min(wc, wk) vary with mA = p + r. Mirrors R `.ExpectedTrit`. -static void expected_trit_by_M(int nI, int nJ, int tc, +static void expected_nrqs_by_M(int nI, int nJ, int tc, std::vector& eM, std::vector& eMAwk, std::vector& eMAwc) { @@ -186,7 +186,7 @@ static void expected_trit_by_M(int nI, int nJ, int tc, } // [[Rcpp::export]] -List trit_expect(const LogicalMatrix splits, const IntegerMatrix characters) { +List nrqs_expect(const LogicalMatrix splits, const IntegerMatrix characters) { const int n_splits = splits.ncol(); const int n_chars = characters.ncol(); NumericMatrix numEdge(n_splits, n_chars); @@ -200,7 +200,7 @@ List trit_expect(const LogicalMatrix splits, const IntegerMatrix characters) { std::vector eM, eWk, eWc; for (int a = 0; a < nStates - 1; ++a) { for (int b = a + 1; b < nStates; ++b) { - expected_trit_by_M(cc.cnt[a], cc.cnt[b], cc.tc, eM, eWk, eWc); + expected_nrqs_by_M(cc.cnt[a], cc.cnt[b], cc.tc, eM, eWk, eWc); for (int s = 0; s < n_splits; ++s) { const int M = cc.mSideA[s]; denM(s, c) += eM[M]; diff --git a/src/expected_mi.cpp b/src/expected_mi.cpp index 05f25c6c3..248575d86 100644 --- a/src/expected_mi.cpp +++ b/src/expected_mi.cpp @@ -34,7 +34,7 @@ inline double l2factorial(int n) { //' block sizes (marginals) of each partition are fixed but the items are //' associated at random. Subtracting this baseline from an observed mutual //' information yields a chance-corrected ("adjusted") mutual information, as -//' applied by [`SiteConcordance`]`(normalize = TRUE)`; it is most material for +//' applied by [`SiteConcordance`]`(chanceCorrect = TRUE)`; it is most material for //' small or unbalanced partitions, where raw mutual information is appreciably //' inflated by chance agreement. //' diff --git a/tests/testthat/test-Concordance.R b/tests/testthat/test-Concordance.R index b455f6012..1856a0707 100644 --- a/tests/testthat/test-Concordance.R +++ b/tests/testthat/test-Concordance.R @@ -20,6 +20,11 @@ test_that("_Concordance() handles tip mismatch", { }) test_that("QuartetConcordance() works", { + # This block locks the raw-quartet contract, so it names `unit` and + # `chanceCorrect` explicitly rather than relying on the (NRQS) defaults. + Naive <- function(...) { + QuartetConcordance(..., unit = "quartet", chanceCorrect = FALSE) + } tree <- BalancedTree(8) splits <- as.Splits(tree) mataset <- matrix(c(0, 0, 0, 0, 1, 1, 1, 1, 0, @@ -32,15 +37,15 @@ test_that("QuartetConcordance() works", { expect_error(QuartetConcordance(tree, mataset), "`dataset` must be a phyDat object") dat <- MatrixToPhyDat(mataset) - expect_equal(unname(QuartetConcordance(tree, dat[, 1])), rep(1, 5)) + expect_equal(unname(Naive(tree, dat[, 1])), rep(1, 5)) # plot(tree); nodelabels(); - expect_equal(QuartetConcordance(tree, dat[, 2]), + expect_equal(Naive(tree, dat[, 2]), c("10" = 1/9, "11" = 0, "12" = 0, "13" = 1/9, "14" = 0, "15" = 0)[names(as.Splits(tree))]) allQuartets <- combn(8, 4) for (charI in seq_len(ncol(mataset))) { - qc <- QuartetConcordance(tree, dat[, charI]) + qc <- Naive(tree, dat[, charI]) for (splitI in seq_along(splits)) { split <- splits[[splitI]] logiSplit <- as.logical(split) @@ -68,7 +73,7 @@ test_that("QuartetConcordance() works", { } } - expect_equal(QuartetConcordance(tree, dat[, c(1:4, 6)]), + expect_equal(Naive(tree, dat[, c(1:4, 6)]), c("10" = (36 + 2 + 9 + 12) / (36 + 18 + 18 + 12 + 6), "11" = ( 6 + 0 + 6 + 2) / ( 6 + 9 + 6 + 2 + 1), "12" = ( 6 + 0 + 0 + 2) / ( 6 + 9 + 9 + 2 + 1), @@ -133,30 +138,34 @@ test_that("QuartetConcordance() handles non-integer data", { QuartetConcordance(tree, MatrixToPhyDat(intSet))) }) -test_that("QuartetConcordance() unit = 'trit' locks the contract", { +test_that("QuartetConcordance() unit = 'nrqs' locks the contract", { tree <- BalancedTree(8) # Character identical to the {t1..t4 | t5..t8} split: that split scores 1, # and it is the *only* split scoring 1 (nested/crossing get partial credit). identChar <- MatrixToPhyDat(matrix( c(0, 0, 0, 0, 1, 1, 1, 1), 8, dimnames = list(paste0("t", 1:8), NULL))) - qt <- QuartetConcordance(tree, identChar, unit = "trit") + qt <- QuartetConcordance(tree, identChar, unit = "nrqs", + chanceCorrect = FALSE) expect_equal(max(qt, na.rm = TRUE), 1) expect_equal(sum(abs(qt - 1) < 1e-9, na.rm = TRUE), 1L) expect_true(all(qt >= 0 & qt <= 1, na.rm = TRUE)) - # unit = "quartet" is the default and is byte-identical to omitting `unit`. + # `unit = "nrqs"` and `chanceCorrect = TRUE` are the defaults, so a bare call + # is byte-identical to naming them. expect_identical(QuartetConcordance(tree, identChar), - QuartetConcordance(tree, identChar, unit = "quartet")) - # Coverage normalisation makes trit no laxer than quartet. + QuartetConcordance(tree, identChar, unit = "nrqs", + chanceCorrect = TRUE)) + # Coverage normalisation makes NRQS no laxer than quartet. expect_true(mean(qt, na.rm = TRUE) <= - mean(QuartetConcordance(tree, identChar, unit = "quartet"), + mean(QuartetConcordance(tree, identChar, unit = "quartet", + chanceCorrect = FALSE), na.rm = TRUE)) - # Uninformative characters (constant / autapomorphy) carry no trits -> NA. + # Uninformative characters (constant / autapomorphy) carry no NRQS -> NA. autap <- MatrixToPhyDat(matrix( c(0, 0, 0, 0, 0, 0, 0, 1), 8, dimnames = list(paste0("t", 1:8), NULL))) - expect_equal(unname(QuartetConcordance(tree, autap, unit = "trit")), + expect_equal(unname(QuartetConcordance(tree, autap, unit = "nrqs")), rep(NA_real_, 5)) # `unit` is validated. @@ -164,22 +173,22 @@ test_that("QuartetConcordance() unit = 'trit' locks the contract", { "should be one of") # `return` aliases mirror the quartet path. - expect_equal(QuartetConcordance(tree, identChar, return = "edge", unit = "trit"), - QuartetConcordance(tree, identChar, return = "default", unit = "trit")) - cA <- QuartetConcordance(tree, identChar, return = "char", unit = "trit") - expect_equal(cA, - QuartetConcordance(tree, identChar, return = "character", unit = "trit")) - expect_equal(cA, - QuartetConcordance(tree, identChar, return = "site", unit = "trit")) + ByReturn <- function(r) { + QuartetConcordance(tree, identChar, return = r, unit = "nrqs") + } + expect_equal(ByReturn("edge"), ByReturn("default")) + cA <- ByReturn("char") + expect_equal(cA, ByReturn("character")) + expect_equal(cA, ByReturn("site")) }) -test_that("QuartetConcordance() unit = 'trit' gives nested partial credit", { +test_that("QuartetConcordance() unit = 'nrqs' gives nested partial credit", { # {t1..t5 | t6,t7,t8} split; state {t1,t2,t4} nested within side A but not a # clade, giving cells (p,q,r,s) = (3,0,2,3) -> A/Wk = 4/8 = 0.5 exactly. tree <- ape::read.tree(text = "(((((t1,t2),t3),t4),t5),(t6,(t7,t8)));") char <- MatrixToPhyDat(matrix( c(0, 0, 1, 0, 1, 1, 1, 1), 8, dimnames = list(paste0("t", 1:8), NULL))) - qt <- QuartetConcordance(tree, char, unit = "trit") + qt <- QuartetConcordance(tree, char, unit = "nrqs", chanceCorrect = FALSE) # The split isolating {t6,t7,t8} scores exactly the nested value 0.5. sp <- as.Splits(tree) @@ -192,20 +201,20 @@ test_that("QuartetConcordance() unit = 'trit' gives nested partial credit", { expect_equal(unname(qt[member]), 0.5) }) -test_that("QuartetConcordance() unit = 'trit' supports multistate", { - # Multistate no longer errors: trits sum over state-pairs (same currency). +test_that("QuartetConcordance() unit = 'nrqs' supports multistate", { + # Multistate no longer errors: NRQS sum over state-pairs (same currency). tree <- BalancedTree(8) ms <- MatrixToPhyDat(matrix( c(0, 0, 1, 1, 2, 2, 2, 0, # 3 states 0, 0, 0, 1, 1, 2, 3, 3), 8, # 4 states dimnames = list(paste0("t", 1:8), NULL))) - qt <- QuartetConcordance(tree, ms, unit = "trit") + qt <- QuartetConcordance(tree, ms, unit = "nrqs", chanceCorrect = FALSE) expect_length(qt, 5L) expect_true(all(qt >= 0 & qt <= 1, na.rm = TRUE)) expect_false(anyNA(qt)) # both characters are informative on this tree }) -test_that("QuartetConcordance() unit = 'trit' scores 1 iff split displayed", { +test_that("QuartetConcordance() unit = 'nrqs' scores 1 iff split displayed", { # A multistate character need not be *identical* to a split to score 1: it # scores full marks for every split its own tree displays (each state block # wholly on one side; the split-orthogonal state-pair drops via M = 0). This @@ -215,7 +224,7 @@ test_that("QuartetConcordance() unit = 'trit' scores 1 iff split displayed", { char <- MatrixToPhyDat(matrix( c(0, 0, 0, 1, 1, 1, 2, 2, 2), 9, dimnames = list(paste0("t", 1:9), NULL))) - qt <- QuartetConcordance(tree, char, unit = "trit") + qt <- QuartetConcordance(tree, char, unit = "nrqs", chanceCorrect = FALSE) sp <- as.Splits(tree) tips <- TipLabels(sp) atSplit <- function(members) { @@ -233,27 +242,61 @@ test_that("QuartetConcordance() unit = 'trit' scores 1 iff split displayed", { expect_true(atSplit(c("t1", "t2")) < 1) }) -test_that("QuartetConcordance() unit = 'trit' normalize contract", { +test_that("QuartetConcordance() 'nrqs' return = 'char' has no ceiling of 1", { + # Documented limitation, not an oversight: the character score is the + # M-weighted mean across EVERY split, and a character agrees exactly with at + # most one of them while being merely compatible with -- silent about -- the + # rest. Coverage scores silence as failure to cover, so even a character + # identical to a split of the tree falls short of 1. `unit = "quartet"` is + # the per-character measure that does attain 1. + tree <- BalancedTree(8) + tips <- TipLabels(tree) + sp <- as.Splits(tree) + ByUnit <- function(u, r) { + QuartetConcordance(tree, char, return = r, unit = u, chanceCorrect = FALSE) + } + for (k in seq_along(sp)) { + char <- MatrixToPhyDat(matrix( + as.integer(as.logical(sp[[k]])), length(tips), + dimnames = list(tips, NULL))) + + # Per-split bound holds: the character scores 1 at the split it duplicates + expect_equal(max(ByUnit("nrqs", "edge")), 1) + + # Naive currency reaches 1 on the character path ... + expect_equal(ByUnit("quartet", "char")[[1]], 1) + + # ... the NRQS currency does not, and must not be read against 1 + nrqsChar <- ByUnit("nrqs", "char")[[1]] + expect_gt(nrqsChar, 0) + expect_lt(nrqsChar, 1) + } +}) + +test_that("QuartetConcordance() unit = 'nrqs' chanceCorrect contract", { tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") ident <- MatrixToPhyDat(matrix( c(0, 0, 0, 0, 0, 0, 1, 1, 1), 9, dimnames = list(paste0("t", 1:9), NULL))) - # normalize = FALSE is the default and leaves the measure untouched. + # chanceCorrect = TRUE is the default; re-zeroing can only lower the measure, + # since a non-negative baseline is subtracted. expect_identical( - QuartetConcordance(tree, ident, unit = "trit", normalize = FALSE), - QuartetConcordance(tree, ident, unit = "trit")) - - # Invalid normalize is rejected. - expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = 0), - "positive integer") - expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = -3), - "positive integer") - expect_error(QuartetConcordance(tree, ident, unit = "trit", normalize = "x"), - "positive integer") + QuartetConcordance(tree, ident, unit = "nrqs", chanceCorrect = TRUE), + QuartetConcordance(tree, ident, unit = "nrqs")) + expect_true(all( + QuartetConcordance(tree, ident, unit = "nrqs", chanceCorrect = TRUE) <= + QuartetConcordance(tree, ident, unit = "nrqs", chanceCorrect = FALSE), + na.rm = TRUE)) + + # Invalid chanceCorrect is rejected. + for (bad in list(0, -3, "x")) { + expect_error(QuartetConcordance(tree, ident, chanceCorrect = bad), + "positive integer") + } # Chance correction also works for the raw quartet currency. - qc <- QuartetConcordance(tree, ident, unit = "quartet", normalize = TRUE) + qc <- QuartetConcordance(tree, ident, unit = "quartet", chanceCorrect = TRUE) expect_type(qc, "double") expect_false(anyNA(qc)) @@ -267,11 +310,12 @@ test_that("QuartetConcordance() unit = 'trit' normalize contract", { setequal(side, c("t7", "t8", "t9")) || setequal(setdiff(tips, side), c("t7", "t8", "t9")) }, logical(1))) - corrected <- QuartetConcordance(tree, ident, unit = "trit", normalize = TRUE) + corrected <- QuartetConcordance(tree, ident, unit = "nrqs", + chanceCorrect = TRUE) expect_equal(unname(corrected[col789]), 1) }) -test_that("QuartetConcordance() unit = 'trit' exact baseline matches MC", { +test_that("QuartetConcordance() unit = 'nrqs' exact baseline matches MC", { tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") dat <- MatrixToPhyDat(matrix( c(0, 0, 0, 0, 0, 0, 1, 1, 1, # displays {t7,t8,t9} @@ -286,11 +330,11 @@ test_that("QuartetConcordance() unit = 'trit' exact baseline matches MC", { # pair's wk and M = min(w_c, w_k) are themselves random under the null. for (ret in c("edge", "char")) { for (w in c(TRUE, FALSE)) { - exact <- QuartetConcordance(tree, dat, unit = "trit", return = ret, - weight = w, normalize = TRUE) + exact <- QuartetConcordance(tree, dat, unit = "nrqs", return = ret, + weight = w, chanceCorrect = TRUE) set.seed(1) - mc <- QuartetConcordance(tree, dat, unit = "trit", return = ret, - weight = w, normalize = 8000) + mc <- QuartetConcordance(tree, dat, unit = "nrqs", return = ret, + weight = w, chanceCorrect = 8000) # Same cells resolve to NA under exact and MC (guards the cell-matching in # the weight = FALSE path), and the finite values agree within MC error. expect_identical(is.na(exact), is.na(mc)) @@ -314,25 +358,28 @@ test_that("QuartetConcordance() unit = 'quartet' exact baseline matches MC", { for (ret in c("edge", "char")) { for (w in c(TRUE, FALSE)) { exact <- QuartetConcordance(tree, dat, unit = "quartet", return = ret, - weight = w, normalize = TRUE) + weight = w, chanceCorrect = TRUE) set.seed(1) mc <- QuartetConcordance(tree, dat, unit = "quartet", return = ret, - weight = w, normalize = 8000) + weight = w, chanceCorrect = 8000) expect_identical(is.na(exact), is.na(mc)) expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04) } } - # normalize = FALSE leaves the published raw measure untouched. + # The published raw measure is recovered by naming both non-default settings, + # and the default no longer returns it. for (ret in c("edge", "char")) for (w in c(TRUE, FALSE)) { - expect_identical( - QuartetConcordance(tree, dat, unit = "quartet", return = ret, weight = w), - QuartetConcordance(tree, dat, unit = "quartet", return = ret, weight = w, - normalize = FALSE)) + raw <- QuartetConcordance(tree, dat, unit = "quartet", return = ret, + weight = w, chanceCorrect = FALSE) + expect_true(all(raw >= 0 & raw <= 1, na.rm = TRUE)) + expect_false(identical( + raw, + QuartetConcordance(tree, dat, return = ret, weight = w))) } }) -test_that("QuartetConcordance() unit = 'trit' correction: conflict below random", { +test_that("QuartetConcordance() 'nrqs' correction: conflict below random", { tree <- ape::read.tree(text = "((((t1,t2),t3),(t4,(t5,t6))),(t7,(t8,t9)));") support <- MatrixToPhyDat(matrix( # agrees with the tree: displays {t7,t8,t9} c(0, 0, 0, 0, 0, 0, 1, 1, 1), 9, @@ -341,11 +388,12 @@ test_that("QuartetConcordance() unit = 'trit' correction: conflict below random" c(0, 0, 1, 1, 0, 0, 1, 1, 0), 9, dimnames = list(paste0("t", 1:9), NULL))) - sChar <- QuartetConcordance(tree, support, return = "char", unit = "trit", - normalize = TRUE) - cNone <- QuartetConcordance(tree, conflict, return = "char", unit = "trit") - cChar <- QuartetConcordance(tree, conflict, return = "char", unit = "trit", - normalize = TRUE) + sChar <- QuartetConcordance(tree, support, return = "char", unit = "nrqs", + chanceCorrect = TRUE) + cNone <- QuartetConcordance(tree, conflict, return = "char", unit = "nrqs", + chanceCorrect = FALSE) + cChar <- QuartetConcordance(tree, conflict, return = "char", unit = "nrqs", + chanceCorrect = TRUE) # A character that agrees with the tree scores above random expectation; # a crossing character is pulled below it, and below its uncorrected score. expect_gt(sChar, 0) @@ -438,13 +486,15 @@ test_that("ClusteringConcordance() gives sensible values", { randomset <- matrix(sample(0:1, 8 * 1000, replace = TRUE), 8, dimnames = list(letters[1:8], NULL)) rat <- MatrixToPhyDat(randomset) - expect_equal(ClusteringConcordance(tree, rat, normalize = TRUE), c("10" = 0), + expect_equal(ClusteringConcordance(tree, rat, chanceCorrect = TRUE), + c("10" = 0), tolerance = 0.05) }) test_that("ClusteringConcordance(return = 'char') Monte-Carlo handles ambiguity", { # Regression for T-330: characters whose ambiguous tokens drop different tips - # give `charSplits` over heterogeneous tip sets. The Monte-Carlo `normalize` + # give `charSplits` over heterogeneous tip sets. The Monte-Carlo + # `chanceCorrect` # path scored a *list* of random trees against that list in one call, which # could not reconcile a common label set ("Old and new labels must match"). tree <- ape::read.tree(text = "((a, b, c, d, e), (f, g, h));") @@ -458,7 +508,7 @@ test_that("ClusteringConcordance(return = 'char') Monte-Carlo handles ambiguity" set.seed(1) # Previously errored: "Old and new labels must match" - cc <- ClusteringConcordance(tree, dat, return = "char", normalize = 10L) + cc <- ClusteringConcordance(tree, dat, return = "char", chanceCorrect = 10L) nChar <- length(attr(dat, "index")) expect_length(cc, nChar) expect_length(attr(cc, "mcse"), nChar) @@ -466,7 +516,8 @@ test_that("ClusteringConcordance(return = 'char') Monte-Carlo handles ambiguity" expect_equal(unname(cc[2]), 1) expect_equal(unname(attr(cc, "mcse")[2]), 0) # Monte-Carlo score never exceeds the un-normalized score (subtracts a baseline) - bare <- ClusteringConcordance(tree, dat, return = "char", normalize = FALSE) + bare <- ClusteringConcordance(tree, dat, return = "char", + chanceCorrect = FALSE) finite <- is.finite(cc) & is.finite(bare) expect_true(all(cc[finite] <= bare[finite] + 1e-8)) }) diff --git a/vignettes/tree-search.Rmd b/vignettes/tree-search.Rmd index 7cdd545a3..d09fa9305 100644 --- a/vignettes/tree-search.Rmd +++ b/vignettes/tree-search.Rmd @@ -174,6 +174,10 @@ smaller number of characters contradicting it. whether a dataset is unanimous or divided in support of a grouping, independently of the method of tree reconstruction. +By default, these are corrected for chance: 1 denotes unanimous support, 0 no +more agreement than would be expected if the character's tokens were shuffled +across the leaves, and negative values conflict. + ```{r concordance} concordance <- QuartetConcordance(strict, vinther) @@ -184,7 +188,7 @@ concordance <- QuartetConcordance(strict, vinther) par(mar = rep(0, 4), cex = 0.8) plot(strict) TreeTools::LabelSplits(strict, signif(concordance, 3), - col = TreeTools::SupportColor(concordance / max(concordance)), + col = TreeTools::SupportColor(pmax(0, concordance)), frame = "none", pos = 3L) ``` From 9ec3c0087736f5663d74e61236fd66c5008bcf68 Mon Sep 17 00:00:00 2001 From: "Martin R. Smith" <1695515+ms609@users.noreply.github.com> Date: Fri, 21 Aug 2026 11:00:37 +0100 Subject: [PATCH 07/16] Simplify @return --- R/Concordance.R | 21 ++++----------------- 1 file changed, 4 insertions(+), 17 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index 90a607e5d..51fc88afe 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -763,7 +763,7 @@ MutualClusteringConcordance <- function(tree, dataset) { #' \eqn{(k - 1)(t - k - 1)} are non-redundant. #' `unit` selects which of the two is counted. By default, agreement is #' measured against the non-redundant content, so that a character scores full -#' marks only where it *displays* the split, rather than merely agreeing with a +#' marks only where it displays the split, rather than merely agreeing with a #' large combinatorial volume of its quartets. #' #' By default, the reported value weights each site by the quartet content it @@ -789,22 +789,9 @@ MutualClusteringConcordance <- function(tree, dataset) { #' `QuartetConcordance(return = "char")` returns a numeric vector giving the #' concordance index calculated at each site, averaged across all splits. #' -#' With `unit = "nrqs"`, note that this average is taken over *every* split in -#' `tree`. A character can agree exactly with at most one split, and is merely -#' compatible with -- that is, silent about -- the rest; because agreement is -#' measured against each split's non-redundant content, silence counts as a -#' failure to cover. Per-character values are therefore much smaller than -#' per-edge ones and **do not reach 1 even for a character identical to one of -#' the tree's own splits** (0.49 for a 22|26 split of a 48-leaf tree, 0.25 for -#' a 2|46 split). The attainable maximum depends on the tree, so these values -#' rank characters against one another but should not be read against an -#' absolute ceiling of 1, nor compared across trees. Where a per-character -#' index bounded at 1 is wanted, use `unit = "quartet"`, whose per-character -#' path attains 1 for a character identical to any split of the tree. -#' -#' With `unit = "quartet"` the same vectors are returned, but scored in the raw -#' quartet currency: values are typically higher, as a character is credited -#' with the full combinatorial volume of quartets that it resolves. +#' With `unit = "nrqs"`, per-character values are much smaller than 1 even +#' where a character perfectly matches a split in a tree. Values rank characters +#' on a given tree, but are not comparable between trees. #' #' @param weight Logical specifying whether to weight sites according to the #' quartet content that they share with each split. From bbfc41e57c184ea50e4605a7687431b3aae5f4e5 Mon Sep 17 00:00:00 2001 From: "Martin R. Smith" <1695515+ms609@users.noreply.github.com> Date: Fri, 21 Aug 2026 11:03:03 +0100 Subject: [PATCH 08/16] Simplify param --- R/Concordance.R | 30 ++++++------------------------ 1 file changed, 6 insertions(+), 24 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index 51fc88afe..7d25e4964 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -796,30 +796,12 @@ MutualClusteringConcordance <- function(tree, dataset) { #' @param weight Logical specifying whether to weight sites according to the #' quartet content that they share with each split. #' @param unit Character specifying the currency in which quartets are counted: -#' - `"nrqs"` (default): only non-redundant quartet statements (NRQS) are -#' counted. Of the \eqn{\binom{k}{2}\binom{t - k}{2}} quartets resolved by -#' a split of sizes \eqn{(k, t - k)}, just \eqn{(k - 1)(t - k - 1)} are -#' non-redundant under the entailment -#' \eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} -#' \insertCite{Nelson1992}{TreeSearch}. Measuring agreement against this -#' reduced content means that, *at a given split*, only a character whose -#' own split is identical scores full marks; nested (compatible) characters -#' receive genuine partial support, and crossing (incompatible) characters -#' score lower still. This bound applies to each split individually, and so -#' to `return = "edge"`; it does not carry over to `return = "char"`, which -#' averages across all splits -- see Value, above. -#' A multistate character is scored in the same currency: -#' NRQS are counted within each pair of states (sizes \eqn{n_i}, \eqn{n_j} -#' give a per-state-pair weight of \eqn{4 / (n_i n_j)}) and summed across -#' pairs, so multistate and binary characters remain directly comparable. -#' Note that non-redundancy is a logical property: NRQS are not mutually -#' independent in a statistical sense, so their count measures quartet -#' content rather than information in bits. -#' - `"quartet"`: each resolved quartet counts once, so a character is -#' credited with the full combinatorial volume of quartets it resolves. -#' This recovers the measure reported by earlier versions of this package, -#' the direct analogue of the site concordance factor -#' \insertCite{Minh2020}{TreeSearch}. +#' - `"nrqs"` (default): non-redundant quartet statements (NRQS). +#' At a given split, only a character whose own split is identical scores +#' full marks; nested (compatible) characters receive partial support, and +#' incompatible characters score lower still. +#' - `"quartet"`: each resolved quartet counts once, +#' analogous to the site concordance factor \insertCite{Minh2020}{TreeSearch}. #' @references \insertAllCited{} #' @importFrom ape keep.tip #' @importFrom cli cli_progress_bar cli_progress_update From 5d5f325d172923ea2454f958f1ff6508f3e8d620 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 21 Aug 2026 11:43:13 +0100 Subject: [PATCH 09/16] refactor: follow the MaxMin package's rename to Coreset `ms609/MaxMin` is now `ms609/Coreset`; every exported solver keeps its name, but the progress option moved from `MaxMin.progress` to `Coreset.progress`. `WideSample()`'s dispatch, its dependency guard, the `Suggests`/`Remotes` entries and the Shiny app's option are updated accordingly, as is the `skip_if_not_installed()` guard in `test-WideSample.R` -- which had been skipping silently, since `MaxMin` can no longer be installed under that name. Three references stay as they are: `ExactMaxMin()` is a function name, the "MaxMin optimum" is the objective rather than the package, and `inst/REFERENCES.bib` holds the title of the cited paper. --- DESCRIPTION | 4 ++-- NEWS.md | 2 +- R/WideSample.R | 34 +++++++++++++++---------------- dev/smoke_40k.R | 2 +- inst/Parsimony/global.R | 2 +- inst/WORDLIST | 1 + man/WideSample.Rd | 10 ++++----- man/dot-WideSampleColumnOracle.Rd | 2 +- man/dot-WideSampleMedoid.Rd | 2 +- tests/testthat/test-WideSample.R | 2 +- 10 files changed, 31 insertions(+), 30 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5c3f5f71b..6d3be8e92 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -55,10 +55,10 @@ Imports: utils, Suggests: cluster, + Coreset, future, highs, knitr, - MaxMin, phangorn (>= 2.2.1), PlotTools, promises, @@ -75,7 +75,7 @@ Suggests: vdiffr (>= 1.0.0), zip, Remotes: - ms609/MaxMin, + ms609/Coreset, Config/Needs/check: callr, pkgbuild, diff --git a/NEWS.md b/NEWS.md index bb13b61f5..55544b82d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -303,7 +303,7 @@ are now private and no longer exported. - `WideSample()` now dispatches to the appropriate Max-Min diversity (MMDP) - solver from the `MaxMin` package, choosing the tier automatically + solver from the `Coreset` package, choosing the tier automatically from `length(trees)`. - New functions `LeastSquaresTree()` and `LeastSquaresFit()` search for, and diff --git a/R/WideSample.R b/R/WideSample.R index 7f698faee..74322c55a 100644 --- a/R/WideSample.R +++ b/R/WideSample.R @@ -13,9 +13,9 @@ #' topologies that sit on broad plateaux and under-represents isolated optima. #' `WideSample()` instead selects for topological *spread*, density-blind, by #' dispatching to the appropriate Max-Min Diversity Problem solver from the -#' \pkg{MaxMin} package: +#' \pkg{Coreset} package: #' -# TODO replace {TreeSearch} refs with {MaxMin} once package on CRAN and +# TODO replace {TreeSearch} refs with {Coreset} once package on CRAN and # imported, and remove refs from inst/REFERENCES.bib (DRY) #' \describe{ #' \item{`FarFirst()` (`effort = 1`)}{Greedy farthest-first selection @@ -93,17 +93,17 @@ #' library("TreeTools") #' trees <- as.phylo(0:99, nTip = 8) #' -#' # WideSample() needs the MaxMin package (Max-Min diversity solvers) -#' if (requireNamespace("MaxMin", quietly = TRUE)) { +#' # WideSample() needs the Coreset package (Max-Min diversity solvers) +#' if (requireNamespace("Coreset", quietly = TRUE)) { #' #' # Fast FarFirst subsample (deterministic, matrix-free) #' sub10 <- WideSample(trees, 10, effort = 1) #' length(sub10) # 10 #' -#' # The remaining tiers (DropAdd, Grasp, exact) all dispatch to 'MaxMin'/ +#' # The remaining tiers (DropAdd, Grasp, exact) all dispatch to 'Coreset'/ #' # 'highs' solvers whose runtime is environment-dependent (e.g. a #' # from-source 'highs' build can be far slower than the CRAN binary) and -#' # whose calling convention tracks 'MaxMin' development. Demonstrate them +#' # whose calling convention tracks 'Coreset' development. Demonstrate them #' # interactively only. #' if (interactive()) { #' # Pre-computed distances @@ -146,17 +146,17 @@ WideSample <- function( effort = NULL, maxSeconds = 60 ) { - if (!requireNamespace("MaxMin", quietly = TRUE)) { - stop("`WideSample()` requires the 'MaxMin' package, which provides the ", + if (!requireNamespace("Coreset", quietly = TRUE)) { + stop("`WideSample()` requires the 'Coreset' package, which provides the ", "Max-Min diversity solvers; install it from ", - "https://github.com/ms609/MaxMin", call. = FALSE) + "https://github.com/ms609/Coreset", call. = FALSE) } # Build ceiling: largest N for which we materialize a dense N x N matrix from # a distance function. ~1.1 GB at 12,000; as.matrix.dist overflows near # 46,340 (the dist half-vector exceeds .Machine$integer.max). buildCeiling <- getOption("WideSample.buildCeiling", 12000L) # Exact ceiling: largest N at which auto-selection reaches the exact tier. - # MaxMin::ExactMaxMin() is now a sparse-matrix, heuristic-warm-started solver + # Coreset::ExactMaxMin() is now a sparse-matrix, heuristic-warm-started solver # (~20x faster than the dense form), practical to a few hundred trees at the # small `n` of interest; beyond that the node-packing IP wall bites (the # MaxMin optimum sits near the diameter, where the threshold graph is @@ -256,13 +256,13 @@ WideSample <- function( } else { .WideSampleColumnOracle(dist, trees, nTrees) } - MaxMin::FarFirst(n, colFn, N = nTrees) + Coreset::FarFirst(n, colFn, N = nTrees) }, # Tier 2: DropAdd returns the bare (sorted) index vector; it runs to its # deterministic plateau, with `maxSeconds` as a safety cap. - `2` = MaxMin::DropAdd(n, dmat, maxSeconds = maxSeconds), + `2` = Coreset::DropAdd(n, dmat, maxSeconds = maxSeconds), # Tier 3: Grasp likewise returns the bare index vector (RNG-dependent). - `3` = MaxMin::Grasp(n, dmat, maxSeconds = maxSeconds), + `3` = Coreset::Grasp(n, dmat, maxSeconds = maxSeconds), # Tier 4: exact solver returns the bare (ascending) index vector, like the # other tiers. `4` = { @@ -272,7 +272,7 @@ WideSample <- function( "or 3 (Grasp), or a larger `maxSeconds`.", immediate. = TRUE) } - MaxMin::ExactMaxMin(k = n, dmat, maxSeconds = maxSeconds) + Coreset::ExactMaxMin(k = n, dmat, maxSeconds = maxSeconds) } ) @@ -331,7 +331,7 @@ WideSample <- function( #' distance to all others. Uses the distance matrix when one is available or #' affordable to build; when only a distance function is supplied for a set too #' large to build a matrix, the central medoid is not affordable, so the -#' deterministic peripheral seed ([MaxMin::FarFirst()] with `k = 1`) is returned +#' deterministic peripheral seed ([Coreset::FarFirst()] with `k = 1`) is returned #' as a matrix-free fallback. #' @return Integer index (1-based) of the selected tree. #' @keywords internal @@ -347,7 +347,7 @@ WideSample <- function( } else { colFn <- .WideSampleColumnOracle(dist, trees, nTrees) # Return: - as.integer(MaxMin::FarFirst(colFn, k = 1L, N = nTrees)) + as.integer(Coreset::FarFirst(colFn, k = 1L, N = nTrees)) } } @@ -355,7 +355,7 @@ WideSample <- function( #' #' Returns a function of one 1-based index `i` giving the distances from tree #' `i` to every tree, as required by the distance-column oracle path of -#' [MaxMin::FarFirst()]. Probes +#' [Coreset::FarFirst()]. Probes #' the `(tree, trees)` calling form once up front and fails clearly if the #' supplied `dist` function does not support it. #' @keywords internal diff --git a/dev/smoke_40k.R b/dev/smoke_40k.R index dfa6b019d..33e972888 100644 --- a/dev/smoke_40k.R +++ b/dev/smoke_40k.R @@ -1,5 +1,5 @@ # 40,000-tree matrix-free smoke test for WideSample() over the -# MaxMin::FarFirst() distance-column oracle path. +# Coreset::FarFirst() distance-column oracle path. suppressPackageStartupMessages({ library(TreeTools) # for as.phylo.numeric library(TreeSearch) diff --git a/inst/Parsimony/global.R b/inst/Parsimony/global.R index 8b1094481..adbcd18c9 100644 --- a/inst/Parsimony/global.R +++ b/inst/Parsimony/global.R @@ -2,7 +2,7 @@ # options("TreeSearch.write.code" = TRUE) # Show code as it is written to log logging <- isTRUE(getOption("TreeSearch.logging")) options(shiny.maxRequestSize = 1024 ^ 3) # Allow max 1 GB files -options(MaxMin.progress = FALSE) # Suppress DropAdd progress messages in app +options(Coreset.progress = FALSE) # Suppress DropAdd progress messages in app # Development: prepend .agent-shiny library so library("TreeSearch") finds # the pre-built v2.0.0 install, preventing pkgload from intercepting and diff --git a/inst/WORDLIST b/inst/WORDLIST index e2d9eb350..378d16fc0 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -32,6 +32,7 @@ Chao CharacterHierarchy Chloranthaceae Congreve +Coreset Cucumel Cucumel's Cynipidae diff --git a/man/WideSample.Rd b/man/WideSample.Rd index 045ef6214..3c51dee92 100644 --- a/man/WideSample.Rd +++ b/man/WideSample.Rd @@ -62,7 +62,7 @@ the likelihood or support for that topology. A random draw over-represents topologies that sit on broad plateaux and under-represents isolated optima. \code{WideSample()} instead selects for topological \emph{spread}, density-blind, by dispatching to the appropriate Max-Min Diversity Problem solver from the -\pkg{MaxMin} package: +\pkg{Coreset} package: \describe{ \item{\code{FarFirst()} (\code{effort = 1})}{Greedy farthest-first selection @@ -111,17 +111,17 @@ which automatic selection reaches the exact tier (default \code{200}).} library("TreeTools") trees <- as.phylo(0:99, nTip = 8) -# WideSample() needs the MaxMin package (Max-Min diversity solvers) -if (requireNamespace("MaxMin", quietly = TRUE)) { +# WideSample() needs the Coreset package (Max-Min diversity solvers) +if (requireNamespace("Coreset", quietly = TRUE)) { # Fast FarFirst subsample (deterministic, matrix-free) sub10 <- WideSample(trees, 10, effort = 1) length(sub10) # 10 -# The remaining tiers (DropAdd, Grasp, exact) all dispatch to 'MaxMin'/ +# The remaining tiers (DropAdd, Grasp, exact) all dispatch to 'Coreset'/ # 'highs' solvers whose runtime is environment-dependent (e.g. a # from-source 'highs' build can be far slower than the CRAN binary) and -# whose calling convention tracks 'MaxMin' development. Demonstrate them +# whose calling convention tracks 'Coreset' development. Demonstrate them # interactively only. if (interactive()) { # Pre-computed distances diff --git a/man/dot-WideSampleColumnOracle.Rd b/man/dot-WideSampleColumnOracle.Rd index df6292d6e..56eb6b279 100644 --- a/man/dot-WideSampleColumnOracle.Rd +++ b/man/dot-WideSampleColumnOracle.Rd @@ -9,7 +9,7 @@ \description{ Returns a function of one 1-based index \code{i} giving the distances from tree \code{i} to every tree, as required by the distance-column oracle path of -\code{\link[MaxMin:FarFirst]{MaxMin::FarFirst()}}. Probes +\code{\link[Coreset:FarFirst]{Coreset::FarFirst()}}. Probes the \verb{(tree, trees)} calling form once up front and fails clearly if the supplied \code{dist} function does not support it. } diff --git a/man/dot-WideSampleMedoid.Rd b/man/dot-WideSampleMedoid.Rd index c42347a0a..78eb74b45 100644 --- a/man/dot-WideSampleMedoid.Rd +++ b/man/dot-WideSampleMedoid.Rd @@ -14,7 +14,7 @@ Returns the index of the most central tree -- the medoid, minimizing summed distance to all others. Uses the distance matrix when one is available or affordable to build; when only a distance function is supplied for a set too large to build a matrix, the central medoid is not affordable, so the -deterministic peripheral seed (\code{\link[MaxMin:FarFirst]{MaxMin::FarFirst()}} with \code{k = 1}) is returned +deterministic peripheral seed (\code{\link[Coreset:FarFirst]{Coreset::FarFirst()}} with \code{k = 1}) is returned as a matrix-free fallback. } \keyword{internal} diff --git a/tests/testthat/test-WideSample.R b/tests/testthat/test-WideSample.R index 634d3626f..f70b14ed0 100644 --- a/tests/testthat/test-WideSample.R +++ b/tests/testthat/test-WideSample.R @@ -1,6 +1,6 @@ # Tier 1: runs on CRAN # Tests for WideSample() — Max-Min diversity (MMDP) tree subsampling -skip_if_not_installed("MaxMin") +skip_if_not_installed("Coreset") test_that("n >= length(trees) returns all trees", { trees <- as.phylo(0:9, nTip = 8) From cd0315e7e8b4510b35788daeef3eff50021770f3 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 21 Aug 2026 11:46:29 +0100 Subject: [PATCH 10/16] Add concordance vocabulary to WORDLIST `spell_check_package()` flags four technical terms, all legitimate: the Nelson-Ladiges surname and the NRQS acronym arrive with the concordance rename, while `multistate` and `unclamped` predate it and were already failing the check. --- inst/WORDLIST | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/inst/WORDLIST b/inst/WORDLIST index 378d16fc0..7759d037d 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -79,6 +79,7 @@ JIA Jamoytius KHELLADI KROGMANN +Ladiges LCM LILJEBLAD LOCONTE @@ -104,10 +105,12 @@ Mesonychidae MorphoBank Morphy MorphyLib +multistate NNI Necrolestes Neopilina Novitates +NRQS O'MEARA Odontogriphus OpenMP @@ -179,6 +182,7 @@ TreeDist TreeTools Trueman UHL +unclamped VERBOOM VILHELMSEN VILLAR From f5b7f3a5b5fb84f3f6b8340e7177425056ee5df0 Mon Sep 17 00:00:00 2001 From: "Martin R. Smith" <1695515+ms609@users.noreply.github.com> Date: Fri, 21 Aug 2026 13:03:26 +0100 Subject: [PATCH 11/16] Simplify chanceCorrect --- R/Concordance.R | 22 +++++----------------- 1 file changed, 5 insertions(+), 17 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index 7d25e4964..d5d95d99d 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -61,21 +61,11 @@ NULL #' split–character pair (`hBest`). #' #' The `chanceCorrect` argument specifies how the zero point is defined. -#' -#' - If `chanceCorrect = FALSE`, zero corresponds to *zero* MI, without -#' correcting for the positive bias that arises because MI is rarely exactly -#' zero in finite samples. -#' -#' - If `chanceCorrect = TRUE`, the expected MI is computed using an analytical -#' approximation based on the distribution of character tokens. This is fast -#' and generally accurate for large trees (~200+ taxa), but does not account -#' for correlation between splits. -#' -#' - If `chanceCorrect` is a positive integer `n`, the expected MI is estimated -#' empirically by fitting each character to `n` uniformly random trees and -#' averaging the resulting MI values. This Monte Carlo approach provides a -#' more accurate baseline for small trees, for which the analytical -#' approximation is biased. Monte Carlo standard errors are returned. +#' If `FALSE`, zero corresponds to *zero* MI. If `TRUE`, the expected MI is +#' approximated: this quick estimate is generally accurate for large trees +#' (~200+ taxa), but does not account for correlation between splits. +#' If `chanceCorrect` is a positive integer `n`, the expected MI is estimated +#' by Monte Carlo sampling from `n` uniformly random trees. #' #' @param return Character specifying the summary to return. Options are: #' - `"edge"`: average concordance of each tree split across all characters; @@ -84,8 +74,6 @@ NULL #' - `"tree"`: an overall tree‑level concordance score; #' - `"all"`: a full array of MI components and normalized values for every #' split–character pair. -#' -#' Matching is case‑insensitive and partial. #' #' @param chanceCorrect Controls the zero point of the concordance scale by #' subtracting the value expected under a chance (fixed-marginal) null, in which From c1cede64ad6be913503216d3d11576e66cf3a4c8 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 21 Aug 2026 13:13:01 +0100 Subject: [PATCH 12/16] docs: give chanceCorrect a single definition The argument was described in both `@details` and `@param`, and the two disagreed on what an integer samples: the details said `n` uniformly random trees, the parameter entry said `n` reshuffles of the character's tokens. Both are true, of different functions. The one surviving statement sits under `@param`, and names the null each function actually draws from: `ClusteringConcordance()` scores each character against `n` trees from `RandomTree()`, whereas `QuartetConcordance()` permutes each character's tokens across the scored leaves, holding its state frequencies and the split sizes fixed -- the same fixed-marginal null that its `chanceCorrect = TRUE` expectation evaluates in closed form. `man/SiteConcordance.Rd` is updated to match. --- R/Concordance.R | 44 +++++++++++---------------------- man/SiteConcordance.Rd | 56 +++++++++++------------------------------- 2 files changed, 30 insertions(+), 70 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index d5d95d99d..b775f6ff3 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -60,12 +60,6 @@ NULL #' so that 1 corresponds to the maximum possible mutual information for each #' split–character pair (`hBest`). #' -#' The `chanceCorrect` argument specifies how the zero point is defined. -#' If `FALSE`, zero corresponds to *zero* MI. If `TRUE`, the expected MI is -#' approximated: this quick estimate is generally accurate for large trees -#' (~200+ taxa), but does not account for correlation between splits. -#' If `chanceCorrect` is a positive integer `n`, the expected MI is estimated -#' by Monte Carlo sampling from `n` uniformly random trees. #' #' @param return Character specifying the summary to return. Options are: #' - `"edge"`: average concordance of each tree split across all characters; @@ -75,29 +69,21 @@ NULL #' - `"all"`: a full array of MI components and normalized values for every #' split–character pair. #' -#' @param chanceCorrect Controls the zero point of the concordance scale by -#' subtracting the value expected under a chance (fixed-marginal) null, in which -#' each character's tokens are reassigned at random across the leaves while its -#' state frequencies and the split sizes are held fixed. -#' - `FALSE`: no chance correction; the measure is scaled only by its maximum, -#' so 1 marks a perfect match and 0 the measure's own (typically positive) -#' floor. -#' - `TRUE`: subtract the expected value under the null, so that 0 marks random -#' expectation and negative values fall below it. `ClusteringConcordance()` -#' uses an analytical approximation to the expected mutual information (fast -#' and generally accurate for large trees, ~200+ taxa, but neglecting -#' correlation between splits); `QuartetConcordance()` uses the exact -#' hypergeometric expectation (for either `unit`). -#' - a positive integer `n`: estimate that expectation by Monte Carlo. -#' `ClusteringConcordance()` fits each character to `n` random trees (and -#' returns Monte-Carlo standard errors, more accurate for small trees where -#' the analytical approximation is biased); `QuartetConcordance()` averages -#' over `n` random reassignments of each character's tokens. -#' -#' In all cases 1 corresponds to the maximum attainable value. For -#' `QuartetConcordance()` (either `unit`), chance-corrected values are returned -#' unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before -#' plotting with [QCol()] / [QACol()]. +#' @param chanceCorrect Sets the zero point of the scale; 1 always marks the +#' maximum attainable value. +#' If `FALSE`, zero corresponds to *zero* MI. +#' If `TRUE`, zero is the value expected when each character's tokens are +#' reassigned at random across the leaves, holding its state frequencies and the +#' split sizes fixed: `QuartetConcordance()` computes this expectation exactly, +#' whereas `ClusteringConcordance()` approximates it, accurately for large trees +#' (~200+ taxa) but neglecting correlation between splits. +#' If a positive integer `n`, the expectation is instead sampled -- over `n` +#' random reassignments of each character's tokens (`QuartetConcordance()`), or +#' against `n` uniformly random trees (`ClusteringConcordance()`, which also +#' returns Monte Carlo standard errors). +#' `QuartetConcordance()` returns chance-corrected values unclamped (they may +#' fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before plotting with [QCol()] / +#' [QACol()]. #' #' @returns #' `ClusteringConcordance(return = "all")` returns a 3D array where each diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index c5f46b630..29706453d 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -44,31 +44,21 @@ split–character pair. Matching is case‑insensitive and partial.} -\item{chanceCorrect}{Controls the zero point of the concordance scale by -subtracting the value expected under a chance (fixed-marginal) null, in which -each character's tokens are reassigned at random across the leaves while its -state frequencies and the split sizes are held fixed. -\itemize{ -\item \code{FALSE}: no chance correction; the measure is scaled only by its maximum, -so 1 marks a perfect match and 0 the measure's own (typically positive) -floor. -\item \code{TRUE}: subtract the expected value under the null, so that 0 marks random -expectation and negative values fall below it. \code{ClusteringConcordance()} -uses an analytical approximation to the expected mutual information (fast -and generally accurate for large trees, ~200+ taxa, but neglecting -correlation between splits); \code{QuartetConcordance()} uses the exact -hypergeometric expectation (for either \code{unit}). -\item a positive integer \code{n}: estimate that expectation by Monte Carlo. -\code{ClusteringConcordance()} fits each character to \code{n} random trees (and -returns Monte-Carlo standard errors, more accurate for small trees where -the analytical approximation is biased); \code{QuartetConcordance()} averages -over \code{n} random reassignments of each character's tokens. -} - -In all cases 1 corresponds to the maximum attainable value. For -\code{QuartetConcordance()} (either \code{unit}), chance-corrected values are returned -unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before -plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} +\item{chanceCorrect}{Sets the zero point of the scale; 1 always marks the +maximum attainable value. +If \code{FALSE}, zero corresponds to \emph{zero} MI. +If \code{TRUE}, zero is the value expected when each character's tokens are +reassigned at random across the leaves, holding its state frequencies and the +split sizes fixed: \code{QuartetConcordance()} computes this expectation exactly, +whereas \code{ClusteringConcordance()} approximates it, accurately for large trees +(~200+ taxa) but neglecting correlation between splits. +If a positive integer \code{n}, the expectation is instead sampled -- over \code{n} +random reassignments of each character's tokens (\code{QuartetConcordance()}), or +against \code{n} uniformly random trees (\code{ClusteringConcordance()}, which also +returns Monte Carlo standard errors). +\code{QuartetConcordance()} returns chance-corrected values unclamped (they may +fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before plotting with +\code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} \item{weight}{Logical specifying whether to weight sites according to the quartet content that they share with each split.} @@ -188,22 +178,6 @@ is quantified using mutual information (MI). All reported values are scaled so that 1 corresponds to the maximum possible mutual information for each split–character pair (\code{hBest}). -The \code{chanceCorrect} argument specifies how the zero point is defined. -\itemize{ -\item If \code{chanceCorrect = FALSE}, zero corresponds to \emph{zero} MI, without -correcting for the positive bias that arises because MI is rarely exactly -zero in finite samples. -\item If \code{chanceCorrect = TRUE}, the expected MI is computed using an analytical -approximation based on the distribution of character tokens. This is fast -and generally accurate for large trees (~200+ taxa), but does not account -for correlation between splits. -\item If \code{chanceCorrect} is a positive integer \code{n}, the expected MI is estimated -empirically by fitting each character to \code{n} uniformly random trees and -averaging the resulting MI values. This Monte Carlo approach provides a -more accurate baseline for small trees, for which the analytical -approximation is biased. Monte Carlo standard errors are returned. -} - \code{MutualClusteringConcordance()} provides a character‑wise summary that emphasises each character’s best‑matching split(s). It treats each character as a simple tree and computes the mutual clustering information between this From b56e4d7182c0857ecefbd5685c8b4a82d4f940f1 Mon Sep 17 00:00:00 2001 From: "Martin R. Smith" <1695515+ms609@users.noreply.github.com> Date: Fri, 21 Aug 2026 13:34:33 +0100 Subject: [PATCH 13/16] Simplify --- R/Concordance.R | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index b775f6ff3..90326cbf7 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -69,21 +69,18 @@ NULL #' - `"all"`: a full array of MI components and normalized values for every #' split–character pair. #' -#' @param chanceCorrect Sets the zero point of the scale; 1 always marks the -#' maximum attainable value. -#' If `FALSE`, zero corresponds to *zero* MI. +#' @param chanceCorrect Sets the zero point of the scale. +#' If `FALSE`, zero corresponds to zero MI. #' If `TRUE`, zero is the value expected when each character's tokens are #' reassigned at random across the leaves, holding its state frequencies and the #' split sizes fixed: `QuartetConcordance()` computes this expectation exactly, #' whereas `ClusteringConcordance()` approximates it, accurately for large trees #' (~200+ taxa) but neglecting correlation between splits. -#' If a positive integer `n`, the expectation is instead sampled -- over `n` -#' random reassignments of each character's tokens (`QuartetConcordance()`), or -#' against `n` uniformly random trees (`ClusteringConcordance()`, which also -#' returns Monte Carlo standard errors). -#' `QuartetConcordance()` returns chance-corrected values unclamped (they may -#' fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before plotting with [QCol()] / -#' [QACol()]. +#' If a positive integer `n`, the expectation is sampled over `n` random reassignments +#' of each character's tokens (`QuartetConcordance()`), or against `n` uniformly random +#' trees (`ClusteringConcordance()`). +#' _Hint: Clamp chance-corrected values to \eqn{[-1, 1]} before plotting with [QCol()] / +#' [QACol()]_. #' #' @returns #' `ClusteringConcordance(return = "all")` returns a 3D array where each From 984ba7de327ae91ffcc1af83c77bbda53b96eb4f Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 21 Aug 2026 13:46:01 +0100 Subject: [PATCH 14/16] docs: regenerate SiteConcordance.Rd for the reworded chanceCorrect --- man/SiteConcordance.Rd | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index 1d20ae0e1..7cc3633e1 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -44,21 +44,18 @@ split–character pair. Matching is case‑insensitive and partial.} -\item{chanceCorrect}{Sets the zero point of the scale; 1 always marks the -maximum attainable value. -If \code{FALSE}, zero corresponds to \emph{zero} MI. +\item{chanceCorrect}{Sets the zero point of the scale. +If \code{FALSE}, zero corresponds to zero MI. If \code{TRUE}, zero is the value expected when each character's tokens are reassigned at random across the leaves, holding its state frequencies and the split sizes fixed: \code{QuartetConcordance()} computes this expectation exactly, whereas \code{ClusteringConcordance()} approximates it, accurately for large trees (~200+ taxa) but neglecting correlation between splits. -If a positive integer \code{n}, the expectation is instead sampled -- over \code{n} -random reassignments of each character's tokens (\code{QuartetConcordance()}), or -against \code{n} uniformly random trees (\code{ClusteringConcordance()}, which also -returns Monte Carlo standard errors). -\code{QuartetConcordance()} returns chance-corrected values unclamped (they may -fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before plotting with -\code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} +If a positive integer \code{n}, the expectation is sampled over \code{n} random reassignments +of each character's tokens (\code{QuartetConcordance()}), or against \code{n} uniformly random +trees (\code{ClusteringConcordance()}). +\emph{Hint: Clamp chance-corrected values to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / +\code{\link[=QACol]{QACol()}}}.} \item{weight}{Logical specifying whether to weight sites according to the quartet content that they share with each split.} From 340e816febb234d9237b5f2569f208bbdb981a20 Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 21 Aug 2026 13:50:51 +0100 Subject: [PATCH 15/16] Redoc; no unneeded WORDLIST exceptions --- NEWS.md | 4 +-- R/Concordance.R | 6 ++--- inst/WORDLIST | 4 --- man/ConcordanceTable.Rd | 37 +++++++++------------------ man/SiteConcordance.Rd | 55 +++++++++-------------------------------- 5 files changed, 28 insertions(+), 78 deletions(-) diff --git a/NEWS.md b/NEWS.md index 8121aaa95..f69a23ba1 100644 --- a/NEWS.md +++ b/NEWS.md @@ -3,8 +3,8 @@ - `QuartetConcordance()` gains `unit`, selecting the currency in which quartets are counted. `unit = "nrqs"` (the new default) counts only non-redundant quartet statements: of the quartets resolved by a split of sizes - (_k_, _t_ − _k_), just (_k_ − 1)(_t_ − _k_ − 1) are non-redundant under the - Nelson–Ladiges entailment, so a character now scores full marks only where it + (_k_, _t_ − _k_), just (_k_ − 1)(_t_ − _k_ − 1) are non-redundant, + so a character now scores full marks only where it displays the split, rather than wherever it agrees with a large combinatorial volume of quartets. Pass `unit = "quartet"` to recover the raw count. diff --git a/R/Concordance.R b/R/Concordance.R index a7678798f..e7f0d5fcd 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -786,7 +786,7 @@ MutualClusteringConcordance <- function(tree, dataset) { #' @param weight Logical specifying whether to weight sites according to the #' quartet content that they share with each split. #' @param unit Character specifying the currency in which quartets are counted: -#' - `"nrqs"` (default): non-redundant quartet statements (NRQS). +#' - `"nrqs"` (default): non-redundant quartet statements. #' At a given split, only a character whose own split is identical scores #' full marks; nested (compatible) characters receive partial support, and #' incompatible characters score lower still. @@ -1014,7 +1014,7 @@ QuartetConcordance <- function( # (ab,cd) + (ab,ce) -> (ab,de) is GF(2) cycle addition, so only the cyclomatic # number (k - 1)(t - k - 1) of them are non-redundant (the NRQS). # -# A multistate character is the disjoint union of its state-pairs: a quartet is +# A multi-state character is the disjoint union of its state-pairs: a quartet is # decisive only when two taxa share one state and two share another, so every # decisive quartet lives in the K_{n_i, n_j} block between two states. Those # blocks are edge-disjoint and the entailment never crosses them (it forces the @@ -1119,7 +1119,7 @@ QuartetConcordance <- function( if (doNorm) { sBase <- ifelse(baseDenM > 0, baseNumEdge / baseDenM, NA_real_) # Average observed and baseline over the SAME (split, char) cells: a - # degenerate multistate pair can be dropped from one but not the other + # degenerate multi-state pair can be dropped from one but not the other # (observed wk = 0 yet E[m] > 0), which would otherwise re-zero mismatched # populations. naMask <- is.na(sEdge) | is.na(sBase) diff --git a/inst/WORDLIST b/inst/WORDLIST index 249173c94..733e626c7 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -32,7 +32,6 @@ Chao CharacterHierarchy Chloranthaceae Congreve -Coreset Cucumel Cucumel's Cynipidae @@ -79,7 +78,6 @@ JIA Jamoytius KHELLADI KROGMANN -Ladiges LCM LILJEBLAD LOCONTE @@ -106,12 +104,10 @@ Mesonychidae MorphoBank Morphy MorphyLib -multistate NNI Necrolestes Neopilina Novitates -NRQS O'MEARA Odontogriphus OpenMP diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index c34e96780..49d53a03d 100644 --- a/man/ConcordanceTable.Rd +++ b/man/ConcordanceTable.Rd @@ -36,31 +36,18 @@ at edges whose descendants are both contain more than \code{largeClade} leaves.} \item{ylab}{Character giving a label for the y axis.} -\item{chanceCorrect}{Controls the zero point of the concordance scale by -subtracting the value expected under a chance (fixed-marginal) null, in which -each character's tokens are reassigned at random across the leaves while its -state frequencies and the split sizes are held fixed. -\itemize{ -\item \code{FALSE}: no chance correction; the measure is scaled only by its maximum, -so 1 marks a perfect match and 0 the measure's own (typically positive) -floor. -\item \code{TRUE}: subtract the expected value under the null, so that 0 marks random -expectation and negative values fall below it. \code{ClusteringConcordance()} -uses an analytical approximation to the expected mutual information (fast -and generally accurate for large trees, ~200+ taxa, but neglecting -correlation between splits); \code{QuartetConcordance()} uses the exact -hypergeometric expectation (for either \code{unit}). -\item a positive integer \code{n}: estimate that expectation by Monte Carlo. -\code{ClusteringConcordance()} fits each character to \code{n} random trees (and -returns Monte-Carlo standard errors, more accurate for small trees where -the analytical approximation is biased); \code{QuartetConcordance()} averages -over \code{n} random reassignments of each character's tokens. -} - -In all cases 1 corresponds to the maximum attainable value. For -\code{QuartetConcordance()} (either \code{unit}), chance-corrected values are returned -unclamped (they may fall below \eqn{-1}); clamp to \eqn{[-1, 1]} before -plotting with \code{\link[=QCol]{QCol()}} / \code{\link[=QACol]{QACol()}}.} +\item{chanceCorrect}{Sets the zero point of the scale. +If \code{FALSE}, zero corresponds to zero MI. +If \code{TRUE}, zero is the value expected when each character's tokens are +reassigned at random across the leaves, holding its state frequencies and the +split sizes fixed: \code{QuartetConcordance()} computes this expectation exactly, +whereas \code{ClusteringConcordance()} approximates it, accurately for large trees +(~200+ taxa) but neglecting correlation between splits. +If a positive integer \code{n}, the expectation is sampled over \code{n} random reassignments +of each character's tokens (\code{QuartetConcordance()}), or against \code{n} uniformly random +trees (\code{ClusteringConcordance()}). +\emph{Hint: Clamp chance-corrected values to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / +\code{\link[=QACol]{QACol()}}}.} \item{plot}{Logical specifying whether to draw the plot.} diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index 7cc3633e1..a6112974b 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -40,9 +40,7 @@ by each character's information content; \item \code{"tree"}: an overall tree‑level concordance score; \item \code{"all"}: a full array of MI components and normalized values for every split–character pair. -} - -Matching is case‑insensitive and partial.} +}} \item{chanceCorrect}{Sets the zero point of the scale. If \code{FALSE}, zero corresponds to zero MI. @@ -62,30 +60,12 @@ quartet content that they share with each split.} \item{unit}{Character specifying the currency in which quartets are counted: \itemize{ -\item \code{"nrqs"} (default): only non-redundant quartet statements (NRQS) are -counted. Of the \eqn{\binom{k}{2}\binom{t - k}{2}} quartets resolved by -a split of sizes \eqn{(k, t - k)}, just \eqn{(k - 1)(t - k - 1)} are -non-redundant under the entailment -\eqn{(ab, cd) + (ab, ce) \rightarrow (ab, de)} -\insertCite{Nelson1992}{TreeSearch}. Measuring agreement against this -reduced content means that, \emph{at a given split}, only a character whose -own split is identical scores full marks; nested (compatible) characters -receive genuine partial support, and crossing (incompatible) characters -score lower still. This bound applies to each split individually, and so -to \code{return = "edge"}; it does not carry over to \code{return = "char"}, which -averages across all splits -- see Value, above. -A multistate character is scored in the same currency: -NRQS are counted within each pair of states (sizes \eqn{n_i}, \eqn{n_j} -give a per-state-pair weight of \eqn{4 / (n_i n_j)}) and summed across -pairs, so multistate and binary characters remain directly comparable. -Note that non-redundancy is a logical property: NRQS are not mutually -independent in a statistical sense, so their count measures quartet -content rather than information in bits. -\item \code{"quartet"}: each resolved quartet counts once, so a character is -credited with the full combinatorial volume of quartets it resolves. -This recovers the measure reported by earlier versions of this package, -the direct analogue of the site concordance factor -\insertCite{Minh2020}{TreeSearch}. +\item \code{"nrqs"} (default): non-redundant quartet statements. +At a given split, only a character whose own split is identical scores +full marks; nested (compatible) characters receive partial support, and +incompatible characters score lower still. +\item \code{"quartet"}: each resolved quartet counts once, +analogous to the site concordance factor \insertCite{Minh2020}{TreeSearch}. }} } \value{ @@ -139,22 +119,9 @@ concordance index calculated at each site, averaged across all splits. Unlike the \code{"edge"} result, this vector is unnamed: characters retain no persistent identifier once collapsed to patterns. -With \code{unit = "nrqs"}, note that this average is taken over \emph{every} split in -\code{tree}. A character can agree exactly with at most one split, and is merely -compatible with -- that is, silent about -- the rest; because agreement is -measured against each split's non-redundant content, silence counts as a -failure to cover. Per-character values are therefore much smaller than -per-edge ones and \strong{do not reach 1 even for a character identical to one of -the tree's own splits} (0.49 for a 22|26 split of a 48-leaf tree, 0.25 for -a 2|46 split). The attainable maximum depends on the tree, so these values -rank characters against one another but should not be read against an -absolute ceiling of 1, nor compared across trees. Where a per-character -index bounded at 1 is wanted, use \code{unit = "quartet"}, whose per-character -path attains 1 for a character identical to any split of the tree. - -With \code{unit = "quartet"} the same vectors are returned, but scored in the raw -quartet currency: values are typically higher, as a character is credited -with the full combinatorial volume of quartets that it resolves. +With \code{unit = "nrqs"}, per-character values are much smaller than 1 even +where a character perfectly matches a split in a tree. Values rank characters +on a given tree, but are not comparable between trees. \code{PhylogeneticConcordance()} returns a numeric vector giving the phylogenetic information of each split in \code{tree}, named according to the @@ -209,7 +176,7 @@ resolves \eqn{\binom{k}{2}\binom{t - k}{2}} quartets, of which only \eqn{(k - 1)(t - k - 1)} are non-redundant. \code{unit} selects which of the two is counted. By default, agreement is measured against the non-redundant content, so that a character scores full -marks only where it \emph{displays} the split, rather than merely agreeing with a +marks only where it displays the split, rather than merely agreeing with a large combinatorial volume of its quartets. By default, the reported value weights each site by the quartet content it From c9b5792e0587a896a3f0afe13c5b93ecb523832c Mon Sep 17 00:00:00 2001 From: R script Date: Fri, 21 Aug 2026 13:55:10 +0100 Subject: [PATCH 16/16] Use bullets --- R/Concordance.R | 10 +++++----- man/ConcordanceTable.Rd | 14 ++++++++------ man/SiteConcordance.Rd | 14 ++++++++------ 3 files changed, 21 insertions(+), 17 deletions(-) diff --git a/R/Concordance.R b/R/Concordance.R index e7f0d5fcd..a532443a0 100644 --- a/R/Concordance.R +++ b/R/Concordance.R @@ -70,15 +70,15 @@ NULL #' split–character pair. #' #' @param chanceCorrect Sets the zero point of the scale. -#' If `FALSE`, zero corresponds to zero MI. -#' If `TRUE`, zero is the value expected when each character's tokens are +#' - If `FALSE`, zero corresponds to zero MI. +#' - If `TRUE`, zero is the value expected when each character's tokens are #' reassigned at random across the leaves, holding its state frequencies and the #' split sizes fixed: `QuartetConcordance()` computes this expectation exactly, #' whereas `ClusteringConcordance()` approximates it, accurately for large trees #' (~200+ taxa) but neglecting correlation between splits. -#' If a positive integer `n`, the expectation is sampled over `n` random reassignments -#' of each character's tokens (`QuartetConcordance()`), or against `n` uniformly random -#' trees (`ClusteringConcordance()`). +#' - If a positive integer `n`, the expectation is sampled over `n` random +#' reassignments of each character's tokens (`QuartetConcordance()`), or against +#' `n` uniformly random trees (`ClusteringConcordance()`). #' _Hint: Clamp chance-corrected values to \eqn{[-1, 1]} before plotting with [QCol()] / #' [QACol()]_. #' diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index 49d53a03d..61b010f79 100644 --- a/man/ConcordanceTable.Rd +++ b/man/ConcordanceTable.Rd @@ -37,17 +37,19 @@ at edges whose descendants are both contain more than \code{largeClade} leaves.} \item{ylab}{Character giving a label for the y axis.} \item{chanceCorrect}{Sets the zero point of the scale. -If \code{FALSE}, zero corresponds to zero MI. -If \code{TRUE}, zero is the value expected when each character's tokens are +\itemize{ +\item If \code{FALSE}, zero corresponds to zero MI. +\item If \code{TRUE}, zero is the value expected when each character's tokens are reassigned at random across the leaves, holding its state frequencies and the split sizes fixed: \code{QuartetConcordance()} computes this expectation exactly, whereas \code{ClusteringConcordance()} approximates it, accurately for large trees (~200+ taxa) but neglecting correlation between splits. -If a positive integer \code{n}, the expectation is sampled over \code{n} random reassignments -of each character's tokens (\code{QuartetConcordance()}), or against \code{n} uniformly random -trees (\code{ClusteringConcordance()}). +\item If a positive integer \code{n}, the expectation is sampled over \code{n} random +reassignments of each character's tokens (\code{QuartetConcordance()}), or against +\code{n} uniformly random trees (\code{ClusteringConcordance()}). \emph{Hint: Clamp chance-corrected values to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / -\code{\link[=QACol]{QACol()}}}.} +\code{\link[=QACol]{QACol()}}}. +}} \item{plot}{Logical specifying whether to draw the plot.} diff --git a/man/SiteConcordance.Rd b/man/SiteConcordance.Rd index a6112974b..c29a9880c 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -43,17 +43,19 @@ split–character pair. }} \item{chanceCorrect}{Sets the zero point of the scale. -If \code{FALSE}, zero corresponds to zero MI. -If \code{TRUE}, zero is the value expected when each character's tokens are +\itemize{ +\item If \code{FALSE}, zero corresponds to zero MI. +\item If \code{TRUE}, zero is the value expected when each character's tokens are reassigned at random across the leaves, holding its state frequencies and the split sizes fixed: \code{QuartetConcordance()} computes this expectation exactly, whereas \code{ClusteringConcordance()} approximates it, accurately for large trees (~200+ taxa) but neglecting correlation between splits. -If a positive integer \code{n}, the expectation is sampled over \code{n} random reassignments -of each character's tokens (\code{QuartetConcordance()}), or against \code{n} uniformly random -trees (\code{ClusteringConcordance()}). +\item If a positive integer \code{n}, the expectation is sampled over \code{n} random +reassignments of each character's tokens (\code{QuartetConcordance()}), or against +\code{n} uniformly random trees (\code{ClusteringConcordance()}). \emph{Hint: Clamp chance-corrected values to \eqn{[-1, 1]} before plotting with \code{\link[=QCol]{QCol()}} / -\code{\link[=QACol]{QACol()}}}.} +\code{\link[=QACol]{QACol()}}}. +}} \item{weight}{Logical specifying whether to weight sites according to the quartet content that they share with each split.}