fit_tree {bnpsd}R Documentation

Fit a tree structure to a coancestry matrix


Implements a heuristic algorithm to find the optimal tree topology based on joining pairs of subpopulations with the highest between-coancestry values, and averaging parent coancestries for the merged nodes (a variant of WPGMA hierarchical clustering stats::hclust()). Branch lengths are optimized using the non-negative least squares approach nnls::nnls(), which minimize the residual square error between the given coancestry matrix and the model coancestry implied by the tree. This algorithm recovers the true tree when the input coancestry truly belongs to this tree and there is no estimation error (i.e. it is the inverse function of coanc_tree()), and performs well for coancestry estimates (i.e. the result of simulating genotypes from the true tree, i.e. from draw_all_admix(), and estimating the coancestry using popkin::popkin() followed by popkin::inbr_diag()).


fit_tree(coancestry, method = "mcquitty")



The coancestry matrix to fit a tree to.


The hierarchical clustering method (passed to stats::hclust()). While all stats::hclust() methods are supported, only two really make sense for this application: "mcquitty" (i.e. WPGMA, default) and "average" (UPGMA).


The tree is bifurcating by construction, but edges may equal zero if needed, effectively resulting in multifurcation, although this code makes no attempt to merge nodes with zero edges. For best fit to arbitrary data, the root edge is always fit to the data (may also be zero). Data fit may be poor if the coancestry does not correspond to a tree, particularly if there is admixture between subpopulations.


A phylo object from the ape package (see ape::read.tree()), with two additional list elements at the end:

See Also

coanc_tree()) for the inverse function (when the coancestry corresponds exactly to a tree).

tree_reorder() for reordering tree structure so that tips appear in a particular desired order.


# create a random tree
k <- 5
tree <- rtree( k )
# true coancestry matrix for this tree
coanc <- coanc_tree( tree )

# now recover tree from this coancestry matrix:
tree2 <- fit_tree( coanc )
# tree2 equals tree!

[Package bnpsd version 1.3.13 Index]