This tutorial focuses on the estimation of ancestral character states for discretely valued traits using a continuous-time Markov chain model commonly known as the Mk model.
"phytools" has a function for ancestral character estimation
under a couple of different models that uses the re-rooting method of Yang
(1995). For other models of evolutionary change not covered by this function,
users should try Rich Fitzjohn's
diversitree
package, or the ace function in ape.
OK - let's pull our data from the Anolis ecomorph tree to use in these analyses:
require(phytools)
data(anoletree)
## this is just to pull out the tip states from the
## data object - normally we would read this from file
x<-getStates(anoletree,"tips")
tree<-anoletree
rm(anoletree)
tree
##
## Phylogenetic tree with 82 tips and 81 internal nodes.
##
## Tip labels:
## Anolis_ahli, Anolis_allogus, Anolis_rubribarbus, Anolis_imias, Anolis_sagrei, Anolis_bremeri, ...
##
## Rooted; includes branch lengths.
x
## Anolis_ahli Anolis_allogus Anolis_rubribarbus
## "TG" "TG" "TG"
## Anolis_imias Anolis_sagrei Anolis_bremeri
## "TG" "TG" "TG"
## Anolis_quadriocellifer Anolis_ophiolepis Anolis_mestrei
## "TG" "GB" "TG"
## Anolis_jubar Anolis_homolechis Anolis_confusus
## "TG" "TG" "TG"
## Anolis_guafe Anolis_garmani Anolis_opalinus
## "TG" "CG" "TC"
## Anolis_grahami Anolis_valencienni Anolis_lineatopus
## "TC" "Tw" "TG"
## Anolis_evermanni Anolis_stratulus Anolis_krugi
## "TC" "TC" "GB"
## Anolis_pulchellus Anolis_gundlachi Anolis_poncensis
## "GB" "TG" "GB"
## Anolis_cooki Anolis_cristatellus Anolis_brevirostris
## "TG" "TG" "Tr"
## Anolis_caudalis Anolis_marron Anolis_websteri
## "Tr" "Tr" "Tr"
## Anolis_distichus Anolis_alumina Anolis_semilineatus
## "Tr" "GB" "GB"
## Anolis_olssoni Anolis_insolitus Anolis_whitemani
## "GB" "Tw" "TG"
## Anolis_haetianus Anolis_breslini Anolis_armouri
## "TG" "TG" "TG"
## Anolis_cybotes Anolis_shrevei Anolis_longitibialis
## "TG" "TG" "TG"
## Anolis_strahmi Anolis_marcanoi Anolis_baleatus
## "TG" "TG" "CG"
## Anolis_barahonae Anolis_ricordii Anolis_cuvieri
## "CG" "CG" "CG"
## Anolis_altitudinalis Anolis_oporinus Anolis_isolepis
## "TC" "TC" "TC"
## Anolis_allisoni Anolis_porcatus Anolis_loysiana
## "TC" "TC" "Tr"
## Anolis_guazuma Anolis_placidus Anolis_sheplani
## "Tw" "Tw" "Tw"
## Anolis_alayoni Anolis_angusticeps Anolis_paternus
## "Tw" "Tw" "Tw"
## Anolis_alutaceus Anolis_inexpectatus Anolis_clivicola
## "GB" "GB" "GB"
## Anolis_cupeyalensis Anolis_cyanopleurus Anolis_alfaroi
## "GB" "GB" "GB"
## Anolis_macilentus Anolis_vanidicus Anolis_baracoae
## "GB" "GB" "CG"
## Anolis_noblei Anolis_smallwoodi Anolis_luteogularis
## "CG" "CG" "CG"
## Anolis_equestris Anolis_bahorucoensis Anolis_dolichocephalus
## "CG" "GB" "GB"
## Anolis_hendersoni Anolis_darlingtoni Anolis_aliniger
## "GB" "Tw" "TC"
## Anolis_singularis Anolis_chlorocyanus Anolis_coelestinus
## "TC" "TC" "TC"
## Anolis_occultus
## "Tw"
In visual format, these are the data that we have:
plotTree(tree,type="fan",fsize=0.8,ftype="i")
## setEnv=TRUE for this type is experimental. please be patient with bugs
cols<-setNames(palette()[1:length(unique(x))],sort(unique(x)))
tiplabels(pie=to.matrix(x,sort(unique(x))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=-max(nodeHeights(tree)),fsize=0.8)
Next, let's fit a single-rate model & reconstruct ancestral states at internal nodes in the tree.
## estimate ancestral states under a ER model
fitER<-ace(x,tree,model="ER",type="discrete")
fitER
##
## Ancestral Character Estimation
##
## Call: ace(x = x, phy = tree, type = "discrete", model = "ER")
##
## Log-likelihood: -78.04604
##
## Rate index matrix:
## CG GB TC TG Tr Tw
## CG . 1 1 1 1 1
## GB 1 . 1 1 1 1
## TC 1 1 . 1 1 1
## TG 1 1 1 . 1 1
## Tr 1 1 1 1 . 1
## Tw 1 1 1 1 1 .
##
## Parameter estimates:
## rate index estimate std-err
## 1 0.0231 0.004
##
## Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
## CG GB TC TG Tr Tw
## 0.018197565 0.202238628 0.042841575 0.428609607 0.004383532 0.303729093
round(fitER$lik.anc,3)
## CG GB TC TG Tr Tw
## [1,] 0.018 0.202 0.043 0.429 0.004 0.304
## [2,] 0.008 0.197 0.020 0.501 0.001 0.272
## [3,] 0.003 0.134 0.020 0.716 0.005 0.122
## [4,] 0.002 0.042 0.019 0.858 0.002 0.077
## [5,] 0.000 0.001 0.001 0.995 0.000 0.002
## [6,] 0.000 0.000 0.000 0.998 0.000 0.001
## [7,] 0.000 0.000 0.000 1.000 0.000 0.000
## [8,] 0.000 0.000 0.000 1.000 0.000 0.000
## [9,] 0.000 0.000 0.000 1.000 0.000 0.000
## [10,] 0.000 0.000 0.000 1.000 0.000 0.000
## [11,] 0.000 0.003 0.000 0.997 0.000 0.000
## [12,] 0.000 0.000 0.000 1.000 0.000 0.000
## [13,] 0.000 0.000 0.000 1.000 0.000 0.000
## [14,] 0.000 0.000 0.000 1.000 0.000 0.000
## [15,] 0.000 0.000 0.000 1.000 0.000 0.000
## [16,] 0.000 0.000 0.000 1.000 0.000 0.000
## [17,] 0.007 0.019 0.084 0.786 0.005 0.099
## [18,] 0.020 0.018 0.314 0.428 0.011 0.210
## [19,] 0.030 0.001 0.945 0.014 0.001 0.007
## [20,] 0.035 0.001 0.945 0.012 0.001 0.006
## [21,] 0.004 0.148 0.036 0.692 0.071 0.049
## [22,] 0.003 0.179 0.057 0.721 0.024 0.017
## [23,] 0.001 0.004 0.979 0.013 0.001 0.001
## [24,] 0.001 0.214 0.010 0.767 0.005 0.003
## [25,] 0.000 0.353 0.003 0.640 0.002 0.001
## [26,] 0.001 0.921 0.002 0.073 0.002 0.002
## [27,] 0.001 0.356 0.003 0.637 0.002 0.001
## [28,] 0.001 0.046 0.003 0.947 0.002 0.001
## [29,] 0.001 0.003 0.001 0.013 0.980 0.002
## [30,] 0.000 0.000 0.000 0.000 1.000 0.000
## [31,] 0.000 0.000 0.000 0.000 1.000 0.000
## [32,] 0.000 0.000 0.000 0.000 1.000 0.000
## [33,] 0.009 0.203 0.018 0.492 0.001 0.276
## [34,] 0.022 0.213 0.013 0.495 0.002 0.255
## [35,] 0.020 0.359 0.014 0.311 0.008 0.287
## [36,] 0.003 0.954 0.002 0.020 0.002 0.019
## [37,] 0.000 0.997 0.000 0.001 0.000 0.001
## [38,] 0.092 0.141 0.012 0.582 0.005 0.168
## [39,] 0.003 0.005 0.001 0.985 0.001 0.005
## [40,] 0.000 0.000 0.000 0.999 0.000 0.000
## [41,] 0.000 0.000 0.000 1.000 0.000 0.000
## [42,] 0.000 0.000 0.000 1.000 0.000 0.000
## [43,] 0.000 0.000 0.000 1.000 0.000 0.000
## [44,] 0.000 0.000 0.000 1.000 0.000 0.000
## [45,] 0.000 0.000 0.000 1.000 0.000 0.000
## [46,] 0.000 0.000 0.000 1.000 0.000 0.000
## [47,] 0.911 0.015 0.005 0.049 0.004 0.017
## [48,] 1.000 0.000 0.000 0.000 0.000 0.000
## [49,] 1.000 0.000 0.000 0.000 0.000 0.000
## [50,] 0.006 0.291 0.042 0.136 0.006 0.519
## [51,] 0.002 0.055 0.071 0.026 0.006 0.840
## [52,] 0.002 0.039 0.118 0.019 0.009 0.812
## [53,] 0.004 0.016 0.653 0.010 0.047 0.270
## [54,] 0.001 0.004 0.927 0.002 0.010 0.056
## [55,] 0.000 0.000 0.999 0.000 0.000 0.001
## [56,] 0.000 0.000 1.000 0.000 0.000 0.000
## [57,] 0.000 0.000 0.995 0.000 0.001 0.003
## [58,] 0.000 0.005 0.006 0.002 0.001 0.985
## [59,] 0.000 0.000 0.000 0.000 0.000 1.000
## [60,] 0.000 0.000 0.000 0.000 0.000 0.999
## [61,] 0.000 0.000 0.000 0.000 0.000 1.000
## [62,] 0.000 0.988 0.001 0.002 0.000 0.008
## [63,] 0.000 1.000 0.000 0.000 0.000 0.000
## [64,] 0.000 1.000 0.000 0.000 0.000 0.000
## [65,] 0.000 1.000 0.000 0.000 0.000 0.000
## [66,] 0.000 1.000 0.000 0.000 0.000 0.000
## [67,] 0.000 1.000 0.000 0.000 0.000 0.000
## [68,] 0.000 1.000 0.000 0.000 0.000 0.000
## [69,] 0.039 0.212 0.096 0.264 0.007 0.381
## [70,] 0.079 0.253 0.212 0.118 0.008 0.329
## [71,] 1.000 0.000 0.000 0.000 0.000 0.000
## [72,] 1.000 0.000 0.000 0.000 0.000 0.000
## [73,] 1.000 0.000 0.000 0.000 0.000 0.000
## [74,] 1.000 0.000 0.000 0.000 0.000 0.000
## [75,] 0.048 0.280 0.278 0.071 0.006 0.317
## [76,] 0.043 0.322 0.238 0.062 0.007 0.328
## [77,] 0.003 0.959 0.013 0.004 0.002 0.018
## [78,] 0.000 0.999 0.000 0.000 0.000 0.000
## [79,] 0.007 0.034 0.907 0.010 0.002 0.039
## [80,] 0.000 0.001 0.997 0.000 0.000 0.001
## [81,] 0.000 0.000 1.000 0.000 0.000 0.000
The element lik.anc gives us the marginal ancestral
states, also known as the 'empirical Bayesian posterior probabilities.'
It is fairly straightforward to overlay these posterior probabilities on the tree:
plotTree(tree,type="fan",fsize=0.8,ftype="i")
## setEnv=TRUE for this type is experimental. please be patient with bugs
nodelabels(node=1:tree$Nnode+Ntip(tree),
pie=fitER$lik.anc,piecol=cols,cex=0.5)
tiplabels(pie=to.matrix(x,sort(unique(x))),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=-max(nodeHeights(tree)),fsize=0.8)
An alternative approach to the one outline above is to use an MCMC approach to sample character histories from their posterior probability distribution. This is called stochastic character mapping (Huelsenbeck et al. 2003). The model is the same but in this case we get a sample of unambiguous histories for our discrete character's evolution on the tree - rather than a probability distribution for the character at nodes.
For instance, given the data simulated above - we can generate the stochastic character map as follows:
## simulate single stochastic character map using empirical Bayes method
mtree<-make.simmap(tree,x,model="ER")
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
## CG GB TC TG Tr Tw
## CG -0.11570723 0.02314145 0.02314145 0.02314145 0.02314145 0.02314145
## GB 0.02314145 -0.11570723 0.02314145 0.02314145 0.02314145 0.02314145
## TC 0.02314145 0.02314145 -0.11570723 0.02314145 0.02314145 0.02314145
## TG 0.02314145 0.02314145 0.02314145 -0.11570723 0.02314145 0.02314145
## Tr 0.02314145 0.02314145 0.02314145 0.02314145 -0.11570723 0.02314145
## Tw 0.02314145 0.02314145 0.02314145 0.02314145 0.02314145 -0.11570723
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## CG GB TC TG Tr Tw
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## Done.
plotSimmap(mtree,cols,type="fan",fsize=0.8,ftype="i")
## setEnv=TRUE for this type is experimental. please be patient with bugs
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=-max(nodeHeights(tree)),fsize=0.8)
A single stochastic character map does not mean a whole lot in isolation - we need to look at the whole distribution from a sample of stochastic maps. This can be a bit overwhelming. For instance, the following code generates 100 stochastic character maps from our dataset and plots them in a grid:
mtrees<-make.simmap(tree,x,model="ER",nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =
## CG GB TC TG Tr Tw
## CG -0.11570723 0.02314145 0.02314145 0.02314145 0.02314145 0.02314145
## GB 0.02314145 -0.11570723 0.02314145 0.02314145 0.02314145 0.02314145
## TC 0.02314145 0.02314145 -0.11570723 0.02314145 0.02314145 0.02314145
## TG 0.02314145 0.02314145 0.02314145 -0.11570723 0.02314145 0.02314145
## Tr 0.02314145 0.02314145 0.02314145 0.02314145 -0.11570723 0.02314145
## Tw 0.02314145 0.02314145 0.02314145 0.02314145 0.02314145 -0.11570723
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## CG GB TC TG Tr Tw
## 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
## Done.
par(mfrow=c(10,10))
null<-sapply(mtrees,plotSimmap,colors=cols,lwd=1,ftype="off")
It's possible to summarize a set of stochastic maps in a much more meaningful way. For instance, we can estimate the number of changes of each type, the proportion of time spent in each state, and the posterior probabilities that each internal node is in each state, under our model. For example:
pd<-describe.simmap(mtrees,plot=FALSE)
pd
## 100 trees with a mapped discrete character with states:
## CG, GB, TC, TG, Tr, Tw
##
## trees have 23.84 changes between states on average
##
## changes are of the following types:
## CG,GB CG,TC CG,TG CG,Tr CG,Tw GB,CG GB,TC GB,TG GB,Tr GB,Tw TC,CG
## x->y 0.31 0.29 0.25 0.16 0.32 0.56 0.95 1.29 0.37 1.31 1.27
## TC,GB TC,TG TC,Tr TC,Tw TG,CG TG,GB TG,TC TG,Tr TG,Tw Tr,CG Tr,GB
## x->y 0.48 0.34 0.65 1.01 0.9 3.11 1.75 0.96 1.57 0.14 0.21
## Tr,TC Tr,TG Tr,Tw Tw,CG Tw,GB Tw,TC Tw,TG Tw,Tr
## x->y 0.23 0.17 0.19 0.74 1.49 1.48 0.87 0.47
##
## mean total time spent in each state is:
## CG GB TC TG Tr Tw
## raw 13.59430464 46.3032646 32.4980763 68.0843311 12.88762444 32.2993595
## prop 0.06609863 0.2251371 0.1580131 0.3310417 0.06266259 0.1570469
## total
## raw 205.667
## prop 1.000
plot(pd,fsize=0.6,ftype="i")
## now let's plot a random map, and overlay the posterior probabilities
plotSimmap(mtrees[[1]],cols,type="fan",fsize=0.8,ftype="i")
## setEnv=TRUE for this type is experimental. please be patient with bugs
nodelabels(pie=pd$ace,piecol=cols,cex=0.5)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=-max(nodeHeights(tree)),fsize=0.8)
Finally, since we obtained these inferences under exactly the same model, let's compare the posterior probabilities from stochastic mapping with our marginal ancestral states. In the former case, our probabilities were obtained by sampling from the joint (rather than marginal) probability distribution for the ancestral states.
plot(fitER$lik.anc,pd$ace,xlab="marginal ancestral states",
ylab="posterior probabilities from stochastic mapping")
lines(c(0,1),c(0,1),lty="dashed",col="red",lwd=2)
This tells us that although joint & marginal reconstruction are not the same, the marginal probabilities from joint stochastic sampling and the marginal ancestral states are quite highly correlated - at least in this case study.
Written by Liam J. Revell. Last updated July 3, 2015