What this vignette is for

A marginal reconstruction asks, node by node, what state was this ancestor in? A joint reconstruction asks a different question: what sequence of states across the whole tree best explains the data? If your biology is about trajectories — when a trait arose, how many times it was lost, whether two clades acquired it independently — the joint reconstruction is the estimate you want, because stringing together the best marginal state at each node can produce a history that is improbable as a path.

The problem is that classical joint algorithms return one history and no uncertainty. This vignette walks through the fix described in Boyko et al. (2025): sample the joint posterior, cluster the samples into a handful of distinct narratives (“reconstruction axes”), and read off where on the tree those narratives actually disagree.

Everything here runs in under a minute on the example data. You will need Rtsne, dbscan and proxy alongside corHMM.

library(corHMM)
library(ape)

state_cols <- c("#F2F2F2", "#7EA6E0", "#B85450")  # one colour per character state

1. Fit a model and get the best joint reconstruction

We use the primates data shipped with corHMM: two binary traits, which corHMM collapses to the three observed combinations 0|0, 0|1 and 1|1.

data(primates)
phy <- ladderize(multi2di(primates[[1]]))
cor_dat <- primates[[2]]

focal <- corHMM(phy, cor_dat,
                rate.cat = 1,
                model = "ER",
                node.states = "joint",
                root.p = "maddfitz")
phy <- ladderize(focal$phy)   # the tree corHMM actually used

focal$states holds the single best joint history. Plotting it gives the estimate most papers stop at.

plot(phy, cex = 0.35, label.offset = 1, no.margin = TRUE)
tiplabels(pch = 21, bg = state_cols[focal$tip.states], cex = 0.6)
nodelabels(pch = 21, bg = state_cols[focal$states], cex = 0.6)
legend("topleft", legend = c("0|0", "0|1", "1|1"), pch = 21,
       pt.bg = state_cols, bty = "n", title = "State")

2. Sample the joint posterior

compute_joint_ci() samples complete histories in proportion to their joint posterior probability, using the two-pass conditional sampling of stochastic character mapping. It keeps sampling in batches until max_samples is reached, discards duplicates, and drops likelihood outliers.

set.seed(0104)
joint_ci <- compute_joint_ci(focal, batch_size = 100, max_samples = 1000)

state_df <- joint_ci$state_df   # rows = reconstructions, columns = internal nodes
n_recon <- nrow(state_df)
n_recon
## [1] 758

That is 758 distinct histories that are all consistent with the same fitted model. Their likelihoods span a narrow band, which is the point: many of them are nearly as good as the best.

par(mar = c(4, 4, 2, 1))
hist(joint_ci$lnliks, breaks = 40, col = "grey80", border = "white",
     main = "", xlab = "log-likelihood of sampled history", ylab = "count")
abline(v = joint_ci$best_lnlik, lty = 2, lwd = 2, col = "#B85450")
text(joint_ci$best_lnlik, par("usr")[4] * 0.9, "best joint  ", adj = 1, col = "#B85450")

Before any clustering, it is worth looking at the raw sample. plot_raw_recon() draws one row per reconstruction and one column per internal node.

plot_raw_recon(as.matrix(state_df),
               state_colors = state_cols,
               state_labels = c("0|0", "0|1", "1|1"))

Most nodes are the same colour in every row — those are settled. A handful of columns are mottled, and those are where the histories genuinely diverge. There is visible vertical structure, but rows arrive in sampling order, so we need to sort them.

3. A distance between reconstructions

Two histories that differ at five scattered nodes are not really telling different stories; two that differ at five adjacent nodes are. phylo_aware_dist_local() encodes exactly that: it sums exp(1 / d) over all pairs of disagreeing nodes, where d is the node distance on a cladogram, so clustered disagreements are penalised much more heavily than isolated ones.

Register it with proxy and build the full distance matrix.

cladogram <- phy
cladogram$edge.length <- rep(1, nrow(cladogram$edge))
phylo_dist <- ape::dist.nodes(cladogram)[-(1:Ntip(phy)), -(1:Ntip(phy))]

try(proxy::pr_DB$set_entry(FUN = phylo_aware_dist_local, names = "PhyloSpreadDistLocal"),
    silent = TRUE)

D <- proxy::dist(as.matrix(state_df),
                 method = "PhyloSpreadDistLocal",
                 phylo_dist = phylo_dist)

4. Consensus clustering

A single t-SNE plus DBSCAN run is quick but seed-dependent, and it is easy to over-read a clustering that would not survive a different seed. The published approach repeats the embedding over many seeds, records how often each pair of reconstructions lands in the same cluster, and clusters that co-association matrix.

Two parameters to set. minPts is the smallest group you are willing to call a cluster; here we use 5% of the sample, though for the several-thousand-reconstruction datasets in the paper 0.5% is more appropriate. eps is the neighbourhood radius, read off the knee of a k-nearest-neighbour distance plot.

min_pts <- round(0.05 * n_recon)
eps <- 3

set.seed(1)
Y1 <- Rtsne::Rtsne(D)$Y
dbscan::kNNdistplot(Y1, k = min_pts - 1, minPts = min_pts)
abline(h = eps, lty = 2, col = "#B85450")

Where the curve has a sharp bend, take eps there. This one rises smoothly, which is itself informative: the reconstructions do not fall into cleanly separated dense islands. We take a value on the lower shoulder, because anything above about 4 swallows the whole sample into a single cluster.

table(dbscan::dbscan(Y1, eps = eps, minPts = min_pts)$cluster)   # 0 = noise
## 
##   0   1   2   3   4 
## 142 192 151 196  77

This is exactly the fragility the consensus step exists to absorb. Repeat across seeds and let agreement, rather than one lucky embedding, decide the grouping. Twenty-five seeds is enough here; the paper used 1,000.

n_seeds <- 25
co <- matrix(0, n_recon, n_recon)

for (s in seq_len(n_seeds)) {
  set.seed(s)
  Y <- Rtsne::Rtsne(D)$Y
  g <- dbscan::dbscan(Y, eps = eps, minPts = min_pts)$cluster
  co <- co + (outer(g, g, "==") & outer(g != 0, g != 0, "&"))
}
co <- co / n_seeds   # probability that a pair clusters together

hc <- hclust(as.dist(1 - co), method = "ward.D2")

Choosing the number of clusters is a judgement call, made by looking at where the dendrogram has room to cut and whether the co-association matrix shows clean blocks.

k <- 4
clusters <- cutree(hc, k = k)
clust_cols <- hcl.colors(k, "Dark 3")

par(mfrow = c(1, 2), mar = c(2, 4, 3, 1))
plot(hc, labels = FALSE, hang = -1, main = "Consensus dendrogram",
     xlab = "", sub = "", ylab = "Ward.D2 height")
# rect.hclust draws left to right, which is not cluster order
rect.hclust(hc, k = k, border = clust_cols[unique(clusters[hc$order])])

ord <- order(clusters)
par(mar = c(3, 3, 3, 1))
image(co[ord, ord], col = hcl.colors(64, "Blues 3", rev = TRUE), axes = FALSE,
      main = "Co-association, ordered by cluster")
box(col = "grey60")

Dark blocks on the diagonal are sets of reconstructions that cluster together no matter the seed. Here 4 blocks are clearly separated, so we cut there. If your matrix looks like one smear, that is a real result: your reconstructions form a continuum rather than distinct narratives, and a single summary is honest.

5. What the axes look like

Re-draw the raw sample, now sorted by cluster.

plot_raw_recon(as.matrix(state_df), clusters, clust_cols,
               state_colors = state_cols,
               state_labels = c("0|0", "0|1", "1|1"))

The mottled columns from before have resolved into blocks: each axis commits to a particular state over a particular stretch of the tree. Section 6 pins down which nodes carry those commitments.

Crucially, these are not ranked alternatives. Their likelihood distributions overlap almost completely, so the data cannot distinguish between them.

par(mar = c(4, 5, 3, 1))
plot_stacked_densities(joint_ci$lnliks, clusters, cols = clust_cols,
                       cluster_labels = paste("axis", seq_len(k)),
                       best_lnlik = joint_ci$best_lnlik,
                       main = "Log-likelihood by reconstruction axis")

6. Where on the tree do they disagree?

node_importance_by_cluster() scores every node for every cluster. A node scores high when the cluster is internally consistent about it and differs from everything outside the cluster, so high scorers are the nodes that define an axis.

node_imp <- node_importance_by_cluster(state_df, clusters)

top_n <- 2
important_nodes <- lapply(seq_len(k), function(i) {
  cl_nodes <- node_imp[node_imp$cluster == i, ]
  cl_nodes[order(-cl_nodes$importance)[seq_len(top_n)], ]
})

do.call(rbind, important_nodes)[, c("cluster", "node", "modal_state", "importance")]
cluster node modal_state importance
44 1 44 2 0.8429553
43 1 43 2 0.8425952
64 2 5 2 0.5383430
65 2 6 2 0.5383430
123 3 5 1 0.8894231
124 3 6 1 0.8894231
218 4 41 1 0.3627418
219 4 42 1 0.3158086

Taking the top 2 nodes per cluster keeps the figure legible and always returns something. A fixed threshold (importance > 0.5, as in the paper) is the alternative, but it can leave a diffuse cluster with no nodes at all, which will break the plot below. Read the scores as well as the ranks: an axis whose best node scores near 0.2 is not pinned to any particular ancestor, and should be described as diffuse rather than as a specific competing history.

Note that important_nodes[[i]] must be cluster i’s nodes — build the list by looping over seq_len(k), not over unique(clusters), which returns clusters in order of first appearance.

par(mar = c(1, 1, 1, 1))
plot_clustered_phylo(
  phy = phy,
  state_df = state_df,
  cluster_assignments = clusters,
  node_states = focal$states,
  important_nodes = important_nodes,
  cluster_colors = clust_cols,
  state_colors = state_cols,
  state_labels = c("0|0", "0|1", "1|1"),
  legend_title = "State",
  phylogram_params = list(edge.color = "grey60"),
  point_params = list(cex = 0.8),
  segment_params = list(lty = 2, lwd = 1.4),
  outer_cex = 3.2, inner_rad = 1.0, node_offset = 2.6
)
legend("topleft", inset = c(0.01, 0.17), legend = paste("axis", seq_len(k)),
       pch = 21, pt.bg = clust_cols, pt.cex = 1.6, bty = "n")

Small points are the best joint reconstruction. Each large ringed pie is a diagnostic node for one axis, coloured by that axis and offset from its position on the tree so overlapping claims stay readable; the pie shows the states that axis assigns there. Because clustering has already separated the sample into distinct narratives, these pies show variation within an axis, not an average over everything — which is exactly what a marginal reconstruction would give you instead.

7. Turning axes into numbers

Reconstruction axes are only useful if they change downstream quantities. Counting transitions along edges is the simplest check.

tips <- focal$tip.states   # already in tree order

count_events <- function(node_states) {
  st <- c(tips, node_states)
  anc <- st[phy$edge[, 1]]
  dec <- st[phy$edge[, 2]]
  c(gains = sum(dec > anc), losses = sum(dec < anc))
}

events <- t(apply(as.matrix(state_df), 1, count_events))

fmt <- function(x) sprintf("%.0f (%.0f-%.0f)", median(x), min(x), max(x))
do.call(rbind, lapply(seq_len(k), function(i) {
  z <- events[clusters == i, , drop = FALSE]
  data.frame(axis = i,
             n = sum(clusters == i),
             gains = fmt(z[, "gains"]),
             losses = fmt(z[, "losses"]),
             median_lnlik = round(median(joint_ci$lnliks[clusters == i]), 2))
}))
axis n gains losses median_lnlik
1 191 7 (4-12) 10 (7-16) -53.58
2 292 10 (6-16) 7 (3-12) -52.99
3 134 10 (7-15) 6 (2-10) -53.14
4 141 10 (6-16) 6 (2-14) -54.46

The axes imply genuinely different histories: some favour more gains than losses, others the reverse, at essentially the same likelihood. A single point estimate would have reported one of these and hidden the rest. Each axis is a concrete, path-consistent hypothesis you can carry into whatever comes next.

Adapting this to your data

Function reference

Function Purpose
compute_joint_ci() Sample joint histories from their posterior
phylo_aware_dist_local() Phylogenetically weighted distance between two histories
plot_raw_recon() Raster of the sample, optionally sorted by cluster
node_importance_by_cluster() Score nodes by how well they diagnose a cluster
plot_stacked_densities() Likelihood distributions per cluster
plot_clustered_phylo() Competing histories drawn on the tree

Reference

Boyko, J.D., Gontjes, K.J., Snitkin, E.S., and Smith, S.A. (2025). Resolving competing evolutionary histories in joint ancestral state reconstruction.