Impact of model choice on tree topologies
Martin R. Smith martin.smith@durham.ac.uk
2026-03-06
Source:vignettes/compare-trees.Rmd
compare-trees.Rmd
library("neotrans")
hetModels <- c("hg_ki", "hg2_ki", "hg_b_ki", "hg_m_ki", "hg_bm_ki")
rmModels <- c("rm_by_n_ki", "rm_by_t_ki", "rm_by_nt_ki")
statModels <- c("by_ki", "by_n_ki", "by_nn_ki", "by_t_ki", "by_nt_ki")
nsModels <- c("ns_ki", "ns_n_ki", "ns_t_ki", "ns_nt_ki")
homModels <- c(statModels, nsModels)
kiModels <- c(statModels, rmModels, nsModels, hetModels)
nonRm <- setdiff(kiModels, rmModels)
if (length(TreeDist::GetParallel()) == 0) {
TreeDist::StartParallel(ceiling(parallel::detectCores() * 0.75))
}
marginals <- GetMarginals(KiProjects(), nonRm)
stdErr <- attr(marginals, "stdErr")
# Which datasets favour a model other than by_ki?
bestModel <- nonRm[apply(marginals, 2, which.max)] |>
`names<-`(KiProjects())
projects <- KiProjects()[bestModel != "by_ki"]
# Compare by_ki with the globally best model
kiLength <- vapply(projects, function(pID) {
median(TreeLengths(pID, "by_ki")) %||% NA_real_
}, double(1))
cfLength <- vapply(projects, function(pID) {
median(TreeLengths(pID, bestModel[[pID]])) %||% NA_real_
}, double(1))## Warning in ConnectSSH(): sshLogin environment variable not set in ConnectSSH()
## Warning in ConnectSSH(): sshLogin environment variable not set in ConnectSSH()
relLen <- cfLength / kiLength
# Estimated trees using the best available stationary model
treeChange <- sapply(projects, function(pID) {
cli::cli_progress_message("Evaluating trees: {pID}")
d <- TreeDistances(pID, bestModel[[pID]], "by_ki")
disp <- Dispersion(d)
prec <- PrecisionIncrease(disp[["spread"]][, "mad"], bestModel[[pID]], "by_ki")
c(
medMed = disp[["mdmd"]] %||% NA_real_,
sil = disp[["sil"]] %||% NA_real_,
prec = prec %||% NA_real_
)
})## Warning in ConnectSSH(): sshLogin environment variable not set in ConnectSSH()
## Warning in ConnectSSH(): sshLogin environment variable not set in ConnectSSH()
cli::cli_progress_done()
modelBF <- sapply(projects, function(pID) {
marginals[bestModel[[pID]], pID] - marginals["by_ki", pID]
})
modelErr <- sapply(projects, function(pID) {
stdErr[bestModel[[pID]], pID] + marginals["by_ki", pID]
})
dashCol <- "grey70"
bfLab <- "Model Bayes factor"
OutputPlot("treeDiff", 7.2, 7.2, function() {
par(mfrow = c(2, 2), mar = c(4.2, 4.2, 0.4, 0.4))
plot(treeChange["sil", ] ~ modelBF,
xlab = bfLab, ylab = "Silhouette score",
#pch = 16, col = ModelCol(bestModel),
pch = 3,
xpd = NA, frame.plot = FALSE)
abline(h = 0, lty = 2, col = dashCol)
EpsLine(v = TRUE, h = FALSE)
Panel(1)
plot(treeChange["prec", ] ~ modelBF,
xlab = bfLab, ylab = "Change in precision (MAD, %)",
xpd = NA, frame.plot = FALSE,
#pch = 16, col = ModelCol(bestModel)
pch = 3)
abline(h = 0, lty = 2, col = dashCol)
EpsLine(h = FALSE, v = TRUE)
vioplot::vioplot(treeChange["prec", ], at = 144, add = TRUE, width = 100,
frame.plot = FALSE, wex = 42)
Panel(2)
plot(treeChange["medMed", ] ~ modelBF,
xlab = bfLab, ylab = "Normalized CI Distance between median trees",
type = "n",
#pch = 16, col = ModelCol(bestModel),
frame.plot = FALSE)
abline(h = 0, lty = 2, col = dashCol)
points(modelBF, treeChange["medMed", ], pch = 3, xpd = NA)
EpsLine(h = FALSE, vert = TRUE)
Panel(3)
x <- treeChange["medMed", ]
x <- x[!is.na(x)]
dens <- density(x, from = 0, to = max(x), bw = "nrd0")
dens$x <- pmax(dens$x, 0) # force support ≥ 0
xr <- c(x, -x) # reflect data around 0
dens <- density(xr, from = 0, to = max(x), bw = "nrd0")
# scale density for plotting
scale <- 12
xpos <- 144
y <- dens$x
z <- dens$y / max(dens$y) * scale
# draw violin manually
polygon(c(xpos - z, rev(xpos + z)),
c(y, rev(y)), col = "grey50", border = "black")
vioplot::vioplot(x, at = 144, add = TRUE, width = 100,
frame.plot = FALSE, wex = 42, col = "transparent", border = NA)
plot(modelBF, relLen,
xlab = bfLab, ylab = "Relative tree length", xpd = NA,
log = "y",
#pch = 16, col = ModelCol(bestModel),
pch = 3,
frame.plot = FALSE)
Panel(4)
abline(h = 1, lty = 2, col = dashCol)
EpsLine(vert = TRUE, horiz = FALSE)
vioplot::vioplot(relLen, at = 144, add = TRUE, width = 100,
frame.plot = FALSE, wex = 42)
})
summary(relLen)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.6128 0.9104 1.0068 1.0336 1.0821 1.7931 2
parsModels <- c("by_ki", "ns_ki", "by_t_ki", "by_n_ki", "by_nt_ki",
"hg_ki", "hg_b_ki")
# Get parsimony scores of Bayesian trees
pScores <- PosteriorTreeSteps(KiProjects(), parsModels)## Warning in .GitClone(pID, scriptID): Cloning into '/home/runner/work/_temp/Library/revbayes-repos/4291_hg_b_ki'...
## fatal: could not read Username for 'https://github.com': terminal prompts disabled## Warning in ConnectSSH(): sshLogin environment variable not set in ConnectSSH()## Couldn't clone to /home/runner/work/_temp/Library/revbayes-repos/4291_hg_b_ki## ## Warning in ConnectSSH(): sshLogin environment variable not set in ConnectSSH()
pScores$model <- factor(pScores$scriptID, levels = parsModels)
library("dplyr")
pNorm <- pScores |>
filter(pID %in% KiProjects()) |>
rename(scriptID_fac = model) |>
group_by(pID) |>
mutate(ew = ew / min(ew)) |>
mutate(k10 = k10 / min(k10)) |>
mutate(k1 = k1 / min(k1)) |>
ungroup()
complete_pIDs <- pNorm |>
count(pID, scriptID) |>
count(pID) |>
filter(n == n_distinct(pNorm$scriptID)) |>
pull(pID)
if (length(complete_pIDs) < length(KiProjects())) {
message("Missing some results from ",
paste0(setdiff(KiProjects(), complete_pIDs), collapse = ", "))
}
pNorm <- pScores |>
filter(pID %in% complete_pIDs) |>
rename(scriptID_fac = model) |>
group_by(pID) |>
mutate(ew = ew / min(ew)) |>
mutate(k10 = k10 / min(k10)) |>
mutate(k1 = k1 / min(k1)) |>
ungroup()
ew_split <- split(pNorm$ew, pNorm$scriptID_fac)
ew_padded <- sapply(ew_split, `length<-`, max(lengths(ew_split)))
k1_split <- split(pNorm$k1, pNorm$scriptID_fac)
k1_padded <- sapply(k1_split, `length<-`, max(lengths(k1_split)))
k10_split <- split(pNorm$k10, pNorm$scriptID_fac)
k10_padded <- sapply(k10_split, `length<-`, max(lengths(k10_split)))Overall parsimony summary
OutputPlot("parsimony_scores", 7.2, 7.2 / 2, function() {
mar <- c(3, 4, 0.4, 0.9)
par(mfrow = c(1, 2), mar = mar, cex = 0.8)
parsWidth <- 5.25
parsBin <- 3.8
SpindlePlot((ew_padded - 1) * 100, ylab = "Extra EW parsimony length / %",
clip = c(0, 0.9), width = parsWidth, nBin = parsBin)
abline(h = 0)
Panel(1)
par(mar = mar + c(0, mar[4], 0, -mar[4]))
SpindlePlot((k10_padded - 1) * 100,
ylab = "Extra IW parsimony length (k = 10) / %",
clip = c(0, 0.88), width = 1.26 * parsWidth, nBin = parsBin)
abline(h = 0)
Panel(2)
})
Parsimony results for each dataset
OutputPlot("parsimony_scores_by_pID", 7.2, 7.2 / 2, function() {
for (pid in unique(pNorm$pID)) {
cat("Dataset", pid, "\n")
thisP <- pScores |>
filter(pID %in% pid) |>
rename(scriptID_fac = model) |>
group_by(pID) |>
mutate(ew = ew / min(ew)) |>
mutate(k10 = k10 / min(k10)) |>
mutate(k1 = k1 / min(k1)) |>
ungroup()
ew_split <- split(thisP$ew, thisP$scriptID_fac)
ew_padded <- sapply(ew_split, `length<-`, max(lengths(ew_split)))
k10_split <- split(thisP$k10, thisP$scriptID_fac)
k10_padded <- sapply(k10_split, `length<-`, max(lengths(k10_split)))
# one page with two panels
par(mfrow = c(1, 2), mar = c(4, 4, 0.6, 0.4), cex = 0.8)
SpindlePlot((ew_padded - 1) * 100,
ylab = "Extra EW parsimony length / %",
clip = 0.95, width = 3.5, nBin = 38,
main = paste("Dataset:", pid))
abline(h = 0)
SpindlePlot((k10_padded - 1) * 100,
ylab = "Extra IW parsimony length (k = 10) / %",
clip = 0.93, width = 3.6, nBin = 38)
abline(h = 0)
}
})## Dataset 104
## Dataset 157
## Dataset 175
## Dataset 450
## Dataset 493
## Dataset 563
## Dataset 635
## Dataset 675
## Dataset 692
## Dataset 706
## Dataset 748
## Dataset 950
## Dataset 1113
## Dataset 1210
## Dataset 1271
## Dataset 2131
## Dataset 2553
## Dataset 2800
## Dataset 3199
## Dataset 3200
## Dataset 3244
## Dataset 3351
## Dataset 3392
## Dataset 3405
## Dataset 3408
## Dataset 3445
## Dataset 3448
## Dataset 3603
## Dataset 3646
## Dataset 3655
## Dataset 3705
## Dataset 3710
## Dataset 3711
## Dataset 3755
## Dataset 3757
## Dataset 3804
## Dataset 3832
## Dataset 3833
## Dataset 3927
## Dataset 3929
## Dataset 4111
## Dataset 4220
## Dataset 4230
## Dataset 4305
## Dataset 4308
## Dataset 4309
## Dataset 4310
## Dataset 4467
## Dataset 4649
## Dataset 4747
## Dataset 4761
## Dataset 4790
## Dataset 4867
## Dataset 4910
## Dataset 5099
## Dataset 5186
## Dataset 5201
## Dataset 5228
## Dataset 5230
## Dataset 5255
## Dataset 5268
## Dataset 5327
## Dataset 07203
### Now for the more fine-scale analyses
nModels <- c("by_n_ki", "by_nt_ki")
mlAvailable <- ModelBF(marginals, nModels)
nProjects <- colnames(mlAvailable)
whichModel <- nModels[apply(mlAvailable, 2, which.max)]
treeChange <- matrix(NA_real_, nrow = 6, ncol = length(nProjects),
dimnames = list(c("nsMedMed", "nsSil", "nsPrec",
"stMedMed", "stSil", "stPrec"), nProjects))
stLength <- double(length(nProjects))
nsLength <- stLength
cfLength <- stLength
nModelBF <- numeric(length(nProjects))
nsModelBF <- nModelBF
eps <- log(10)
nullParam <- c(`2.5%` = NA_real_, `25%` = NA_real_, `50%` = NA_real_,
`75%` = NA_real_, `97.5%` = NA_real_, mad = NA_real_)
# Compare estimated trees
for (model in nModels) {
usesModel <- whichModel == model
modProjects <- nProjects[usesModel]
nsModel <- sub("by_", "ns_", model, fixed = TRUE)
cfModel <- c("by_n_ki" = "by_ki", "by_nt_ki" = "by_t_ki")[model]
stLength[usesModel] <- vapply(modProjects, function(pID) {
median(TreeLengths(pID, model)) %||% NA_real_
}, double(1))
nsLength[usesModel] <- vapply(modProjects, function(pID) {
median(TreeLengths(pID, nsModel)) %||% NA_real_
}, double(1))
cfLength[usesModel] <- vapply(modProjects, function(pID) {
median(TreeLengths(pID, cfModel)) %||% NA_real_
}, double(1))
# Estimated trees using the best available stationary model
treeChange[, modProjects] <- sapply(modProjects, function(pID) {
cli::cli_progress_message("Evaluating trees: {pID}")
stD <- TreeDistances(pID, cfModel, model)
stDisp <- Dispersion(stD)
stPrec <- PrecisionIncrease(stDisp[["spread"]][, "mad"], cfModel, model)
nsD <- TreeDistances(pID, cfModel, nsModel)
nsDisp <- Dispersion(nsD)
nsPrec <- PrecisionIncrease(nsDisp[["spread"]][, "mad"], cfModel, nsModel)
c(
nsMedMed = nsDisp[["mdmd"]] %||% NA_real_,
nsSil = nsDisp[["sil"]] %||% NA_real_,
nsPrec = nsPrec %||% NA_real_,
stMedMed = stDisp[["mdmd"]] %||% NA_real_,
stSil = stDisp[["sil"]] %||% NA_real_,
stPrec = stPrec %||% NA_real_
)
})
cli::cli_progress_done()
nModelBF[usesModel] <- marginals[model, modProjects] -
marginals[cfModel, modProjects]
# Estimated value of n using the best available non-stationary model
# nonStatN[, usesModel] <- sapply(modProjects, function(pID) {
# ExistingResults(pID, nsModel, checkRemote = FALSE)[["parameters"]][, "rate_loss"] %||% nullParam
# })
nsModelBF[usesModel] <- marginals[nsModel, modProjects] -
marginals[cfModel, modProjects]
rm(modProjects) # for safety
}## Evaluating trees: 104## Evaluating trees: 175## Evaluating trees: 450## Evaluating trees: 493## Evaluating trees: 635## Evaluating trees: 675## Evaluating trees: 950## Evaluating trees: 1113## Evaluating trees: 2131## Evaluating trees: 3199## Evaluating trees: 3351## Evaluating trees: 3405## Evaluating trees: 3448## Evaluating trees: 3603## Evaluating trees: 3646## Evaluating trees: 3655## Evaluating trees: 3755## Evaluating trees: 3757## Evaluating trees: 3804## Evaluating trees: 3832## Evaluating trees: 4111## Evaluating trees: 4220## Evaluating trees: 4230## Evaluating trees: 4291## Evaluating trees: 4305## Evaluating trees: 4310## Evaluating trees: 4467## Evaluating trees: 4747## Evaluating trees: 4790## Evaluating trees: 5099## Evaluating trees: 5186## Evaluating trees: 5201## Evaluating trees: 5255## Evaluating trees: 07203## Evaluating trees: 157## Evaluating trees: 563## Evaluating trees: 692## Evaluating trees: 706## Evaluating trees: 748## Evaluating trees: 1210## Evaluating trees: 1271## Evaluating trees: 2553## Evaluating trees: 2800## Evaluating trees: 3200## Evaluating trees: 3244## Evaluating trees: 3392## Evaluating trees: 3408## Evaluating trees: 3445## Evaluating trees: 3705## Evaluating trees: 3710## Evaluating trees: 3711## Evaluating trees: 3833## Evaluating trees: 3927## Evaluating trees: 3929## Evaluating trees: 4308## Evaluating trees: 4309## Evaluating trees: 4649## Evaluating trees: 4761## Evaluating trees: 4867## Evaluating trees: 4910## Evaluating trees: 5228## Evaluating trees: 5230## Evaluating trees: 5268## Evaluating trees: 5327
stRelLen <- stLength / cfLength
nsRelLen <- nsLength / cfLength
# Estimated trees using the best available stationary model
treeChange <- rbind(treeChange, sapply(colnames(treeChange), function(pID) {
cli::cli_progress_message("Evaluating trees: {pID}")
hgD <- TreeDistances(pID, "by_ki", "hg_ki")
hgDisp <- Dispersion(hgD)
hgPrec <- PrecisionIncrease(hgDisp[["spread"]][, "mad"], "by_ki", "hg_ki")
c(
hgMedMed = hgDisp[["mdmd"]] %||% NA_real_,
hgSil = hgDisp[["sil"]] %||% NA_real_,
hgPrec = hgPrec %||% NA_real_
)
}))## Evaluating trees: 104## Evaluating trees: 157## Evaluating trees: 175## Evaluating trees: 450## Evaluating trees: 493## Evaluating trees: 563## Evaluating trees: 635## Evaluating trees: 675## Evaluating trees: 692## Evaluating trees: 706## Evaluating trees: 748## Evaluating trees: 950## Evaluating trees: 1113## Evaluating trees: 1210## Evaluating trees: 1271## Evaluating trees: 2131## Evaluating trees: 2553## Evaluating trees: 2800## Evaluating trees: 3199## Evaluating trees: 3200## Evaluating trees: 3244## Evaluating trees: 3351## Evaluating trees: 3392## Evaluating trees: 3405## Evaluating trees: 3408## Evaluating trees: 3445## Evaluating trees: 3448## Evaluating trees: 3603## Evaluating trees: 3646## Evaluating trees: 3655## Evaluating trees: 3705## Evaluating trees: 3710## Evaluating trees: 3711## Evaluating trees: 3755## Evaluating trees: 3757## Evaluating trees: 3804## Evaluating trees: 3832## Evaluating trees: 3833## Evaluating trees: 3927## Evaluating trees: 3929## Evaluating trees: 4111## Evaluating trees: 4220## Evaluating trees: 4230## Evaluating trees: 4291## Evaluating trees: 4305## Evaluating trees: 4308## Evaluating trees: 4309## Evaluating trees: 4310## Evaluating trees: 4467## Evaluating trees: 4649## Evaluating trees: 4747## Evaluating trees: 4761## Evaluating trees: 4790## Evaluating trees: 4867## Evaluating trees: 4910## Evaluating trees: 5099## Evaluating trees: 5186## Evaluating trees: 5201## Evaluating trees: 5228## Evaluating trees: 5230## Evaluating trees: 5255## Evaluating trees: 5268## Evaluating trees: 5327## Evaluating trees: 07203
cli::cli_progress_done()
# Compare estimated trees in HG model
hgModelBF <- marginals["hg_ki", nProjects] - marginals["by_ki", nProjects]
hgLength <- vapply(nProjects, function(pID) {
median(TreeLengths(pID, "hg_ki")) %||% NA_real_
}, double(1))
cfLength <- vapply(nProjects, function(pID) {
median(TreeLengths(pID, "by_ki")) %||% NA_real_
}, double(1))
hgRelLen <- hgLength / cfLength
relLenMax <- max(c(stRelLen, nsRelLen, hgRelLen,
1 / stRelLen, 1 / nsRelLen, 1 / hgRelLen),
na.rm = TRUE)
relLenRange <- c(1 / relLenMax, relLenMax)
bfRange <- range(c(nModelBF, nsModelBF, hgModelBF), na.rm = TRUE)
# Plotting symbols
stPch <- 3
nsPch <- 4
hgPch <- 5
stCol <- ModelCol("bestStat")
nsCol <- ModelCol("bestNS")
hgCol <- ModelCol("bestHet")
dashCol <- "grey70"
par(mfrow = c(2, 2))
plot(nModelBF, stRelLen,
xlab = "Model BF", ylab = "Relative tree length", xpd = NA,
xlim = bfRange, ylim = relLenRange, log = "y",
pch = stPch, col = stCol,
frame.plot = FALSE)
points(nsModelBF, nsRelLen, pch = nsPch, col = nsCol, xpd = NA)
points(hgModelBF, hgRelLen, pch = hgPch, col = hgCol, xpd = NA)
abline(h = 1, lty = 2, col = dashCol)
abline(v = c(1, -1) * eps, lty = 2, col = dashCol)
plot(treeChange["stSil", ] ~ nModelBF,
xlab = "Model BF", ylab = "Silhouette score",
xlim = bfRange,
pch = stPch, col = stCol,
xpd = NA, frame.plot = FALSE)
points(treeChange["nsSil", ] ~ nsModelBF, pch = nsPch, col = nsCol, xpd = NA)
points(treeChange["hgSil", ] ~ hgModelBF, pch = hgPch, col = hgCol, xpd = NA)
abline(h = 0, lty = 2, col = dashCol)
abline(v = c(1, -1) * eps, lty = 2, col = dashCol)
legend("bottomright", legend = c("Stationary n", "Non-stationary n", "Heterogeneous"),
pch = c(stPch, nsPch, hgPch), col = c(stCol, nsCol, hgCol), bty = "n")
plot(treeChange["stMedMed", ] ~ nModelBF,
xlab = "Model BF", ylab = "Normalized CI Distance between median trees",
ylim = c(0, max(treeChange[c("stMedMed", "nsMedMed", "hgMedMed"), ], na.rm = TRUE)),
pch = stPch, col = stCol,
xpd = NA, frame.plot = FALSE)
points(treeChange["nsMedMed", ] ~ nsModelBF, pch = nsPch, col = nsCol, xpd = NA)
points(treeChange["hgMedMed", ] ~ hgModelBF, pch = hgPch, col = hgCol, xpd = NA)
abline(h = 0, lty = 2, col = dashCol)
abline(v = c(1, -1) * eps, lty = 2, col = dashCol)
plot(treeChange["stPrec", ] ~ nModelBF,
xlab = "Model BF", ylab = "Change in precision (MAD, %)",
xpd = NA, frame.plot = FALSE,
xlim = bfRange,
pch = stPch, col = stCol,
ylim = range(treeChange[c("stPrec", "nsPrec", "hgPrec"), ], na.rm = TRUE))
points(treeChange["nsPrec", ] ~ nsModelBF, pch = nsPch, col = nsCol, xpd = NA)
points(treeChange["hgPrec", ] ~ hgModelBF, pch = hgPch, col = hgCol, xpd = NA)
abline(h = 0, lty = 2, col = dashCol)
abline(v = c(1, -1) * eps, lty = 2, col = dashCol)
Sample a representative subset of parsimony trees by identifying
k trees that are as distant as possible from one
another.
(Better algorithms likely exist!)
MostDistant <- function(d, k = 12) {
n <- attr(d, "Size")
if (n < k) {
return(seq_len(n))
}
d <- as.matrix(d)
selected <- integer(k)
meanD <- rowMeans(d)
selected[[1]] <- which.max(meanD)
if (k < 2) {
return(selected)
}
for (i in 2:k) {
remaining <- setdiff(1:n, selected[1:(i-1)])
# For each remaining point, find its minimum distance to any selected point
min_distances <- sapply(remaining, function(j) {
min(d[j, selected[1:(i-1)]])
})
# Select the point with the maximum minimum distance
selected[[i]] <- remaining[which.max(min_distances)]
}
selected
}