mcmc2anc: Ancestral character reconstruction from an MCMC sample

In dosreislab/mcmc3r: Tools to work with MCMCtree and BPP

Ancestral character reconstruction from an MCMC sample

Description

Obtain the ancestral reconstruction of characters from an MCMC sample using Felsenstein (1973) model of continuous character evolution.

Usage

mcmc2anc(
  tree,
  M,
  mcmc,
  time.name,
  rate.name,
  tip.ages = NULL,
  method = c("quick", "proper"),
  thin
)

Arguments

tree	a rooted, strictly bifurcating phylogeny
M	s x k matrix of k continuous morphological measurements for s species
mcmc	data frame with MCMC output from MCMCtree
time.name	character vector of length one
rate.name	character vector of length one
tip.ages	numeric, the ages of the terminal taxa in the tree
method	character, the reconstruction method, see details
thin	numeric, the fraction of mcmc samples to keep, ignored if method == "quick"

Details

If method == "quick", the function first calculates the mean of divergence times and morphological rates in the MCMC sample. These are used to reconstruct the branch lengths in units of morphological evolution, and then Eq. (7) in Felsenstein (1973) is used to calculate the ancestral reconstruction. This results in a "mean" reconstruction along each internal node of the phylogeny. This method is meant to be quick for exploratory data analysis.

If method == "proper", the function applies Eq. (7) in Felsenstein (1973) to each observation from the thinned MCMC. The result is a posterior sample of reconstructions at each internal node. This method should be preferred.

Note time.name is the name format used for the node ages in the MCMC dataframe, usually of the form time.name = "t_". Similarly rate.name is the name format used in the MCMC sample for the rates, of the form rate.name = "r_g1_" where the subscript in "g" must be the partition number containing the morphological rates. Note tree must be the same used by MCMCtree when generating the MCMC sample, right down to its Newick representation. Taxon names in tree and M must match.

Value

If method == "quick", the function returns a n x k matrix with the "mean" ancestral reconstruction for the k characters at the n internal nodes of the phylogeny. If method == "proper", it returns a n x k x N array, where N is the number of observations in the thinned MCMC.

Author(s)

Mario dos Reis

References

Felsenstein J (1973) Maximum-likelihood estimation of evolutionary trees from continuous characters. Am J Hum Genet, 25: 471–492.

See Also

write.morpho

Examples

data(carnivores) 
# calculate tip ages correctly, as they're needed by the function:
tips_info <- strsplit(carnivores$tree$tip.label, "\\^" )
back.ages <- as.numeric(unlist(tips_info)[seq(from=2, to=2*19, by=2)])
back.ages <- max(back.ages) - back.ages 
C <- carnivores$C.proc
rownames(C) <- rownames(carnivores$M)
recM = mcmc2anc(carnivores$tree, C, mcmc=carnivores$mcmc, time.name="t_", 
       rate.name="r_g1_", tip.ages=back.ages)

x <- seq(1, 87, by=3); y <- seq(2, 87, by=3)
# Plot landmarks for the 19 carnivores:
plot(carnivores$C.proc[,x], carnivores$C.proc[,y], pch='+', cex=.5)
# plot ancestral reconstruction at the root (node 20):
points(recM["20",x], recM["20",y], pch=19, col="red")
# mean shape
mS <- apply(carnivores$C.proc, 2, mean)
points(mS[x], mS[y], pch=19, col="blue")

# full (proper) reconstruction
recF <- mcmc2anc(carnivores$tree, C, mcmc=carnivores$mcmc, time.name="t_", 
       rate.name="r_g1_", tip.ages=back.ages, method="proper", thin=.0125)

# calculate mean across full mcmc reconstruction
mSF <- apply(recF, c(1,2), mean)

# prepare area for plotting reconstruction:
plot(mSF["20",x], mSF["20",y], ty='n')
# add posterior sample for the root
for (i in 1:dim(recF)[3]) 
   points(recF["20",x,i], recF["20",y,i], cex=.7, pch='+', col="darkgray")
# add proper mean reconstruction at the root
points(mSF["20",x], mSF["20",y], pch=19, col="orange")
# add quick mean reconstruction at the root (from previous step)
points(recM["20",x], recM["20",y], pch='+', col="red")

## Not run: 
# Convert the quick reconstruction to an array, as is the standard in
# morphometrics software
recA <- matrix2array(recM, 3)
options(rgl.printRglwidget = TRUE)
rgl::plot3d(recA[,,"20"], ty='s', size=2, col="red", aspect=FALSE)

## End(Not run)

dosreislab/mcmc3r documentation built on April 5, 2025, 4:06 p.m.