diff --git a/NEWS.md b/NEWS.md index 8fa93f015..9cdab6971 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. + - `MaximizeParsimony()` now contracts zero-length (unsupported) branches into polytomies by default (`collapse = TRUE`), deduplicating the returned trees on the resulting collapsed topologies, à la TNT's "collapse zero-length diff --git a/R/Concordance.R b/R/Concordance.R index 2ba797245..7d25e4964 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,15 +87,29 @@ 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 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()]. #' #' @returns #' `ClusteringConcordance(return = "all")` returns a 3D array where each @@ -116,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. #' @@ -158,7 +172,7 @@ ClusteringConcordance <- function( tree, dataset, return = "edge", - normalize = TRUE + chanceCorrect = TRUE ) { # Check inputs if (is.null(dataset)) { @@ -243,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)) @@ -257,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(, )` @@ -271,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) } } @@ -295,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 @@ -308,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) @@ -344,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) @@ -528,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) @@ -729,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 @@ -741,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%. @@ -765,8 +789,20 @@ 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"`, 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 -#' 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: +#' - `"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 #' @importFrom utils combn @@ -776,7 +812,9 @@ QuartetConcordance <- function( tree, dataset = NULL, weight = TRUE, - return = "edge" + return = "edge", + unit = c("nrqs", "quartet"), + chanceCorrect = TRUE ) { if (is.null(dataset)) { warning("Cannot calculate concordance without `dataset`.") @@ -785,6 +823,15 @@ QuartetConcordance <- function( if (!inherits(dataset, "phyDat")) { stop("`dataset` must be a phyDat object.") } + unit <- match.arg(unit) + 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.") + } + } + } tipLabels <- intersect(TipLabels(tree), names(dataset)) if (!length(tipLabels)) { warning("No overlap between tree labels and dataset.") @@ -816,14 +863,37 @@ QuartetConcordance <- function( dimnames = dimnames(characters) ) + options <- c("character", "site", "default") + return <- options[[pmatch(tolower(trimws(return)), options, + nomatch = length(options))]] + + if (unit == "nrqs") { + # Return: + return(.NrqsConcordance(logiSplits, charInt, weight, return, splits, + chanceCorrect)) + } + 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))]] - + + # Chance correction (option): re-zero the observed concordant/decisive ratio + # against the ratio expected under the same fixed-marginal null used for the + # 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 NRQS) + # or by Monte-Carlo tip-shuffle, then re-zero the pooled ratio. + doNorm <- !isFALSE(chanceCorrect) + if (doNorm) { + base <- if (isTRUE(chanceCorrect)) { + .QuartetExpect(charInt, logiSplits) + } else { + .QuartetMC(charInt, logiSplits, chanceCorrect) + } + eNum <- base[["concordant"]] + eDen <- base[["decisive"]] + } if (return == "default") { if (isTRUE(weight)) { @@ -836,37 +906,358 @@ 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) + ) + } + } + } +} + +# 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) { + quartet_expect(logiSplits, charInt) +} + +.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 ("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 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 NRQS add over state pairs: +# W_char = sum_{i 0` guard. + be <- nrqs_expect(logiSplits, charInt) + baseNumEdge <- be[["numEdge"]] + baseNumChar <- be[["numChar"]] + baseDenM <- be[["denM"]] + } + } + + for (ci in seq_len(nChar)) { + col <- charInt[, ci] + obs <- .CharNrqsContrib(col, logiSplits, pos) + numEdge[, ci] <- obs[["numEdge"]] + numChar[, ci] <- obs[["numChar"]] + denM[, ci] <- obs[["denM"]] + wcTot[ci] <- obs[["wcTot"]] + 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"]] + } + } + + 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) + 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 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_) + # 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)) + } else { + # char: one value per character + if (isTRUE(weight)) { + denom <- colSums(denM) + 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_) + if (doNorm) { + sBase <- ifelse(baseDenM > 0, baseNumChar / baseDenM, NA_real_) + # 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 + } + } +} + +# 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. +.CharNrqsContrib <- 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 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 + 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) +} + +# 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 `.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. +.CharNrqsMC <- 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 <- .CharNrqsContrib(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()) diff --git a/R/RcppExports.R b/R/RcppExports.R index 100dd4d15..a66ea624d 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) +} + +nrqs_expect <- function(splits, characters) { + .Call(`_TreeSearch_nrqs_expect`, splits, characters) +} + #' Expected mutual information between two partitions #' #' Computes the mutual information expected purely by chance between two @@ -27,7 +35,7 @@ MaddisonSlatkin_clear_cache <- function() { #' 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 new file mode 100644 index 000000000..38187da82 --- /dev/null +++ b/dev/benchmarks/frac-quart/chance_baseline.R @@ -0,0 +1,86 @@ +# Validate the exact hypergeometric chance baseline used by +# 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 +# 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/benchmarks/frac-quart/cpp_expect_parity.R b/dev/benchmarks/frac-quart/cpp_expect_parity.R new file mode 100644 index 000000000..10972edd3 --- /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 +# `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. +# +# 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 <- 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)) + 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/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..b29e053c6 --- /dev/null +++ b/dev/benchmarks/frac-quart/validate_trit.R @@ -0,0 +1,229 @@ +# 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). +# 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 = "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)) +} + +# 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 = "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)) +} + +# 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 = "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 = "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", + 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", 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", + 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..7fe19b3d4 --- /dev/null +++ b/dev/plans/frac-quart-weight-agreement.md @@ -0,0 +1,272 @@ +# 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` 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: 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 + pair by the *edge's* content, which is incoherent as shared information. + `min(W_c,W_k)` also drops into `ConcordanceTable`'s `hBest·n` slot (§5). diff --git a/inst/REFERENCES.bib b/inst/REFERENCES.bib index 0c16caba7..8fd1e150c 100644 --- a/inst/REFERENCES.bib +++ b/inst/REFERENCES.bib @@ -1,24 +1,4 @@ -@article{LapointeCucumel1997, - title = {The average consensus procedure: combination of weighted trees containing identical or overlapping sets of taxa}, - author = {Lapointe, Fran{\c c}ois-Joseph and Cucumel, Guy}, - year = {1997}, - journal = {Systematic Biology}, - volume = {46}, - number = {2}, - pages = {306--312}, - doi = {10.1093/sysbio/46.2.306} -} - -@article{Arias2004, - title = {Profile parsimony ({{PP}}): an analysis under implied weights ({{IW}})}, - author = {Arias, J. Salvador and {Miranda-Esquivel}, Daniel Rafael}, - year = {2004}, - volume = {20}, - pages = {56--63}, - doi = {10.1111/j.1096-0031.2003.00001.x}, - journal = {Cladistics}, - number = {1} -} +@comment{List entries alphabetically by key} @article{Aberer2013, title = {Pruning rogue taxa improves phylogenetic accuracy: an efficient algorithm and webservice}, @@ -42,6 +22,17 @@ @article{Archie1989 doi = {10.2307/2992286} } +@article{Arias2004, + title = {Profile parsimony ({{PP}}): an analysis under implied weights ({{IW}})}, + author = {Arias, J. Salvador and {Miranda-Esquivel}, Daniel Rafael}, + year = {2004}, + volume = {20}, + pages = {56--63}, + doi = {10.1111/j.1096-0031.2003.00001.x}, + journal = {Cladistics}, + number = {1} +} + @article{Bien2011, title = {Hierarchical Clustering with Prototypes via Minimax Linkage}, author = {Bien, Jacob and Tibshirani, Robert}, @@ -93,6 +84,27 @@ @article{Carter1990 doi = {10.1137/0403005} } +@article{Chao1984, + author = {Chao, Anne}, + title = {Nonparametric estimation of the number of classes in a population}, + journal = {Scandinavian Journal of Statistics}, + year = {1984}, + volume = {11}, + number = {4}, + pages = {265--270} +} + +@article{Chao2012, + author = {Chao, Anne and Jost, Lou}, + title = {Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size}, + journal = {Ecology}, + year = {2012}, + volume = {93}, + number = {12}, + pages = {2533--2547}, + doi = {10.1890/11-1952.1} +} + @article{Congreve2016, author = {Congreve, Curtis R. and Lamsdell, James C.}, doi = {10.1111/pala.12236}, @@ -357,6 +369,17 @@ @article{Gonzalez1985 doi = {10.1016/0304-3975(85)90224-5} } +@article{Good1953, + author = {Good, I. J.}, + title = {The population frequencies of species and the estimation of population parameters}, + journal = {Biometrika}, + year = {1953}, + volume = {40}, + number = {3/4}, + pages = {237--264}, + doi = {10.2307/2333344} +} + @article{Hartigan1979, title = {Algorithm {{AS}} 136: A {{{\emph{K}}}}-Means Clustering Algorithm}, shorttitle = {Algorithm {{AS}} 136}, @@ -400,6 +423,17 @@ @article{Kumar2012 number = {2} } +@article{LapointeCucumel1997, + title = {The average consensus procedure: combination of weighted trees containing identical or overlapping sets of taxa}, + author = {Lapointe, Fran{\c c}ois-Joseph and Cucumel, Guy}, + year = {1997}, + journal = {Systematic Biology}, + volume = {46}, + number = {2}, + pages = {306--312}, + doi = {10.1093/sysbio/46.2.306} +} + @article{Maddison1991, title = {Null models for the number of evolutionary steps in a character on a phylogenetic tree}, author = {Maddison, Wayne P. and Slatkin, M}, @@ -474,6 +508,28 @@ @article{Muller2005 number = {1} } +@article{Nelson1992, + title = {Information content and fractional weight of three-taxon statements}, + author = {Nelson, Gareth and Ladiges, Pauline Y.}, + year = {1992}, + journal = {Systematic Biology}, + volume = {41}, + number = {4}, + pages = {490--494}, + doi = {10.1093/sysbio/41.4.490} +} + +@article{Nguyen2015, + title = {{IQ-TREE}: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies}, + author = {Nguyen, Lam-Tung and Schmidt, Heiko A. and von Haeseler, Arndt and Minh, Bui Quang}, + year = {2015}, + volume = {32}, + pages = {268--274}, + doi = {10.1093/molbev/msu300}, + journal = {Molecular Biology and Evolution}, + number = {1} +} + @article{Nixon1999, author = {Nixon, Kevin C.}, doi = {10.1111/j.1096-0031.1999.tb00277.x}, @@ -705,6 +761,17 @@ @article{Sun2018 doi = {10.1098/rspb.2018.1780} } +@article{Thompson1933, + title = {On the Likelihood that One Unknown Probability Exceeds Another in View of the Evidence of Two Samples}, + author = {Thompson, William R.}, + year = {1933}, + volume = {25}, + pages = {285--294}, + doi = {10.2307/2332286}, + journal = {Biometrika}, + number = {3/4} +} + @article{Vinther2008, author = {Vinther, Jakob and {Van Roy}, Peter and Briggs, Derek E. G.}, doi = {10.1038/nature06474}, @@ -753,57 +820,3 @@ @article{Wilkinson2017 number = {4} } -@article{Nguyen2015, - title = {{IQ-TREE}: A Fast and Effective Stochastic Algorithm for Estimating Maximum-Likelihood Phylogenies}, - author = {Nguyen, Lam-Tung and Schmidt, Heiko A. and von Haeseler, Arndt and Minh, Bui Quang}, - year = {2015}, - volume = {32}, - pages = {268--274}, - doi = {10.1093/molbev/msu300}, - journal = {Molecular Biology and Evolution}, - number = {1} -} - -@article{Thompson1933, - title = {On the Likelihood that One Unknown Probability Exceeds Another in View of the Evidence of Two Samples}, - author = {Thompson, William R.}, - year = {1933}, - volume = {25}, - pages = {285--294}, - doi = {10.2307/2332286}, - journal = {Biometrika}, - number = {3/4} -} - -@article{Chao1984, - author = {Chao, Anne}, - title = {Nonparametric estimation of the number of classes in a population}, - journal = {Scandinavian Journal of Statistics}, - year = {1984}, - volume = {11}, - number = {4}, - pages = {265--270} -} - -@article{Good1953, - author = {Good, I. J.}, - title = {The population frequencies of species and the estimation of population parameters}, - journal = {Biometrika}, - year = {1953}, - volume = {40}, - number = {3/4}, - pages = {237--264}, - doi = {10.2307/2333344} -} - -@article{Chao2012, - author = {Chao, Anne and Jost, Lou}, - title = {Coverage-based rarefaction and extrapolation: standardizing samples by completeness rather than size}, - journal = {Ecology}, - year = {2012}, - volume = {93}, - number = {12}, - pages = {2533--2547}, - doi = {10.1890/11-1952.1} -} - diff --git a/man/ConcordanceTable.Rd b/man/ConcordanceTable.Rd index 9d6a0831b..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,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{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; 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()} 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 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()} (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 7dc208b54..c5f46b630 100644 --- a/man/SiteConcordance.Rd +++ b/man/SiteConcordance.Rd @@ -9,11 +9,18 @@ \alias{SharedPhylogeneticConcordance} \title{Concordance factors} \usage{ -ClusteringConcordance(tree, dataset, return = "edge", normalize = TRUE) +ClusteringConcordance(tree, dataset, return = "edge", chanceCorrect = TRUE) MutualClusteringConcordance(tree, dataset) -QuartetConcordance(tree, dataset = NULL, weight = TRUE, return = "edge") +QuartetConcordance( + tree, + dataset = NULL, + weight = TRUE, + return = "edge", + unit = c("nrqs", "quartet"), + chanceCorrect = TRUE +) PhylogeneticConcordance(tree, dataset) @@ -37,20 +44,62 @@ 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{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; 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()} 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 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()} (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.} +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}. +}} } \value{ \code{ClusteringConcordance(return = "all")} returns a 3D array where each @@ -71,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. @@ -99,6 +148,23 @@ 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 = "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 split's internal numbering. @@ -122,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 @@ -145,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 @@ -157,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\%. @@ -236,6 +311,8 @@ myMatrix <- congreveLamsdellMatrices[[10]] ClusteringConcordance(TreeTools::NJTree(myMatrix), myMatrix) } \references{ +\insertAllCited{} + \insertAllCited{} } \seealso{ 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 091d463c4..c0d66c83f 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 +} +// 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(nrqs_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 46e70baf8..a0b440028 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_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); @@ -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_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 new file mode 100644 index 000000000..b6b201920 --- /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(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 +// (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 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_nrqs_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 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); + 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_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]; + numEdge(s, c) += eWk[M]; + numChar(s, c) += eWc[M]; + } + } + } + } + return List::create(_["numEdge"] = numEdge, _["numChar"] = numChar, + _["denM"] = denM); +} diff --git a/src/expected_mi.cpp b/src/expected_mi.cpp index f8a3cdef8..6886cb7e9 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 05a75f53f..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,6 +138,269 @@ test_that("QuartetConcordance() handles non-integer data", { QuartetConcordance(tree, MatrixToPhyDat(intSet))) }) +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 = "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 = "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 = "nrqs", + chanceCorrect = TRUE)) + # Coverage normalisation makes NRQS no laxer than quartet. + expect_true(mean(qt, na.rm = TRUE) <= + mean(QuartetConcordance(tree, identChar, unit = "quartet", + chanceCorrect = FALSE), + na.rm = TRUE)) + + # 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 = "nrqs")), + rep(NA_real_, 5)) + + # `unit` is validated. + expect_error(QuartetConcordance(tree, identChar, unit = "trits"), + "should be one of") + + # `return` aliases mirror the quartet path. + 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 = '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 = "nrqs", chanceCorrect = FALSE) + + # 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 = '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 = "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 = '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 + # 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 = "nrqs", chanceCorrect = FALSE) + 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("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))) + + # 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 = "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", chanceCorrect = 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 + # 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 = "nrqs", + chanceCorrect = TRUE) + expect_equal(unname(corrected[col789]), 1) +}) + +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} + 0, 0, 0, 0, 0, 0, 0, 1, 1, # nested: {t8,t9} + 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 -- 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 = "nrqs", return = ret, + weight = w, chanceCorrect = TRUE) + set.seed(1) + 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)) + expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04) + } + } +}) + +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, chanceCorrect = TRUE) + set.seed(1) + mc <- QuartetConcordance(tree, dat, unit = "quartet", return = ret, + weight = w, chanceCorrect = 8000) + expect_identical(is.na(exact), is.na(mc)) + expect_lt(max(abs(exact - mc), na.rm = TRUE), 0.04) + } + } + + # 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)) { + 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() '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, + 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 = "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) + 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) }) @@ -218,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));") @@ -238,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) @@ -246,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 846f8b972..3535ffb23 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) ```