First bash at IQ-TREE impl

Maths not right yet

R/Concordance.R

@@ -34,6 +34,9 @@
 #' Briefly, this interprets each edge as defining *four* clades, and counts
 #' the status of quartets that contain exactly one leaf from each of these
 #' clades.  This is a subset of the quartets considered by other methods.
+#' `method = "iqtree"` uses the \insertCite{Minh2020;textual}{TreeSearch}
+#' approach as [implemented in IQ-TREE](
+#' https://github.com/iqtree/iqtree2/issues/415).
 #' 
 #' 
 #' `QuartetConcordance()` is computed exactly, using all quartets.
@@ -113,7 +116,7 @@ QuartetConcordance <- function (tree, dataset = NULL, method = "split") {
   dataset <- dataset[tipLabels, drop = FALSE]
   characters <- PhyDatToMatrix(dataset, ambigNA = TRUE)
 
-  if (method == "minh") {
+  if (method %in% c("minh", "iqtree")) {
     if (attr(tree, "order") != "preorder") {
       stop("Tree must be in preorder; try `tree <- Preorder(tree)`")
     }
@@ -165,11 +168,9 @@ QuartetConcordance <- function (tree, dataset = NULL, method = "split") {
         tab <- table(q, char)
         nCol <- dim(tab)[[2]]
         if (nCol > 1L) {
-
           # RULES:
           # We must pick one leaf per row
-          # To be parsimony informative, we must pick two leaves from each of
-          #   two columns
+          # Only parsimony informative are decisive
           # To be concordant, leaves 1 and 2 must come from the same column
 
           # Example table:
@@ -189,19 +190,54 @@ QuartetConcordance <- function (tree, dataset = NULL, method = "split") {
               prod(tab[1:2, j], tab[3:4, i])
           }))
 
-          decisive <- sum(apply(combn(seq_len(nCol), 2), 2, function(ij) {
-            # With three columns, we'll encounter i = 1, 1, 2; j = 2, 3, 3
-            i <- ij[[1]]
-            j <- ij[[2]]
-            sum(apply(combn(4, 2), 2, function(kl) {
-              # We can pick taxa AB, AC, AD, BC, BD, CD from column i,
-              # and the other pair from col j
-              pair1 <- kl
-              pair2 <- (1:4)[-kl]
-              prod(tab[pair1, i], tab[pair2, j])
+          if (method == "minh") {
+            # To be parsimony informative, we must pick two leaves from each of
+            #   two columns
+
+
+            decisive <- sum(apply(combn(seq_len(nCol), 2), 2, function(ij) {
+              # With three columns, we'll encounter i = 1, 1, 2; j = 2, 3, 3
+              i <- ij[[1]]
+              j <- ij[[2]]
+              sum(apply(combn(4, 2), 2, function(kl) {
+                # We can pick taxa AB, AC, AD, BC, BD, CD from column i,
+                # and the other pair from col j
+                pair1 <- kl
+                pair2 <- (1:4)[-kl]
+                prod(tab[pair1, i], tab[pair2, j])
+              }))
             }))
-          }))
-
+          } else {
+            # Under the IQ-definition, to be parsimony informative,
+            # we must pick two leaves from one column.
+            # Hence 0012 is parsimony informative (!)
+
+            parsInf <- sum(vapply(seq_len(nCol), function(i) {
+              sum(apply(combn(4, 2), 2, function(kl) {
+                # We can pick taxa AB, AC, AD, BC, BD, CD from column i,
+                # and the other pair from any other column
+                pair1 <- kl
+                prod(tab[pair1, i], rowSums(tab[-pair1, -i, drop = FALSE]))
+              }))
+            }, double(1)))
+
+            # inclusion-exclusion principle:
+            countedTwice <- sum(apply(combn(seq_len(nCol), 2), 2, function(ij) {
+              # With three columns, we'll encounter i = 1, 1, 2; j = 2, 3, 3
+              i <- ij[[1]]
+              j <- ij[[2]]
+              sum(apply(combn(4, 2), 2, function(kl) {
+                # We can pick taxa AB, AC, AD, BC, BD, CD from column i,
+                # and the other pair from col j
+                pair1 <- kl
+                pair2 <- (1:4)[-kl]
+                prod(tab[pair1, i], tab[pair2, j])
+              }))
+            }))
+
+            decisive <- parsInf - countedTwice
+
+          }
           c(concordant, decisive)
         } else {
           c(0L, 0L)

tests/testthat/test-Concordance.R

@@ -34,9 +34,11 @@ test_that("QuartetConcordance(method = minh)", {
   tree <- Preorder(tree)
 
   # plot(tree); nodelabels();
-  expect_concordance <- function(i, expectation) {
-    expect_equal(QuartetConcordance(tree, dat[, i], method = "minh"),
-                 setNames(expectation, 11:15))
+  expect_concordance <- function(i, expectation, iq = "minh") {
+    expect_equal(
+      QuartetConcordance(tree, dat[, i],
+                         method = ifelse(iq == "iq", "iqtree", "minh")),
+      setNames(expectation, 11:15))
   }
   # Expectations computed by working through tables manually
   expect_concordance(1, c(NaN, NaN, 1, NaN, NaN))
@@ -49,7 +51,6 @@ test_that("QuartetConcordance(method = minh)", {
   expect_concordance(8, c(NaN, NaN, 0, NaN, NaN))
   expect_concordance(9, c(NaN, 0, 1 / 2, 0, NaN))
 
-  expect_concordance(6, rep(NaN, 5))
   # Values calculated from summing results above
   expect_equal(unname(QuartetConcordance(tree, dat, method = "minh")), 
                c(5 + 3 + 1,
@@ -63,8 +64,19 @@ test_that("QuartetConcordance(method = minh)", {
                   6 + 2,
                   4 + 3))
 
-  expect_equal(unname(QuartetConcordance(tree, dat, method = "minh")), 
+  # Expectations computed by iq-tree
+  expect_concordance(iq = "iq", 1, c(NaN, NaN, 1, NaN, NaN))
+  expect_concordance(iq = "iq", 2, c(1, NaN, NaN, NaN, NaN))
+  expect_concordance(iq = "iq", 3, c(0.6, NaN, NaN, NaN, 0.5))
+  expect_concordance(iq = "iq", 4, c(0, 0, 2 / 3, NaN, NaN))
+  expect_concordance(iq = "iq", 5, c(0, 0, 1 / 3, 0, 3 / 8))
+  expect_concordance(iq = "iq", 6, rep(NaN, 5))
+  expect_concordance(iq = "iq", 7, c(1 / 5, NaN, NaN, NaN, NaN))
+  expect_concordance(iq = "iq", 8, c(NaN, NaN, 0, NaN, NaN))
+  expect_concordance(iq = "iq", 9, c(NaN, 0, 8.29, 0, NaN))
+  expect_equal(unname(QuartetConcordance(tree, dat, method = "iqtree")), 
                c(56.7, 0, 85.4, 0, 62.5), tolerance = 0.01)
+
 })
 
 test_that("QuartetConcordance() calculates correct values - weighting", {