Exercise 7: Reconstructing ancestral states for continuous & discrete characters
For better or for worse, the estimation of phenotypic trait values for ancestral nodes in the tree continues to be an important goal in phylogenetic comparative biology.
In this tutorial, we will review ancestral state reconstruction for both continuous & discrete character traits.
Continuous characters
Anole body size evolution
Let's start by estimating ancestral states for a a continuous character, in this case overall body size, on a phylogeny of Anolis lizards from the Caribbean.
To estimate the states for a continuously valued character at ancestral nodes, we need find the states that have the maximum probability of having arisen under our model. These will be our maximum likelihood estimates.
Ancestral character esimation is implemented in a variety of different R functions. The most commonly
used is ace. Here, we start with the phytools function fastAnc.
The files you will need are here:
## load libraries
library(phytools)
## read tree from file
tree<-read.tree("Anolis.tre")
Let's plot our tree first:
## plot tree
plotTree(tree,type="fan",ftype="i",lwd=1)
## read data
obj<-read.csv("anole.data.csv",row.names=1)
svl<-setNames(obj$SVL,rownames(obj))
Now we can estimate ancestral states. We will also compute variances & 95% confidence intervals for each node:
fit<-fastAnc(tree,svl,vars=TRUE,CI=TRUE)
print(fit,printlen=10)
## Ancestral character estimates using fastAnc:
## 101 102 103 104 105 106 107 108
## 4.065918 4.0456 4.047993 4.066923 4.036744 4.049515 4.05453 4.04522
## 109 110
## 3.979493 3.95244 ....
##
## Variances on ancestral states:
## 101 102 103 104 105 106 107 108
## 0.011775 0.00841 0.009452 0.01187 0.01446 0.018113 0.011686 0.007269
## 109 110
## 0.014488 0.012814 ....
##
## Lower & upper 95% CIs:
## lower upper
## 101 3.853231 4.278604
## 102 3.86586 4.22534
## 103 3.857438 4.238548
## 104 3.853382 4.280465
## 105 3.801057 4.27243
## 106 3.785732 4.313299
## 107 3.842651 4.266408
## 108 3.878115 4.212326
## 109 3.743578 4.215408
## 110 3.730566 4.174314
## .... ....
## as we discussed in class, 95% CIs can be broad
range(svl) ## compare to root node
## [1] 3.46201 5.11399
We will learn more about this later, but there are also some methods for visualizing the ancestral state reconstructions for continuous traits on the tree. For example:
## projection of the reconstruction onto the edges of the tree
obj<-contMap(tree,svl,plot=FALSE)
obj
## Object of class "contMap" containing:
##
## (1) A phylogenetic tree with 100 tips and 99 internal nodes.
##
## (2) A mapped continuous trait on the range (3.46201, 5.11399).
plot(obj,type="fan",legend=0.7*max(nodeHeights(tree)),
sig=2,fsize=c(0.7,0.9))
Properties of ancestral states of ancestral state reconstruction
Next, we can explore some of the properties of ancestral state reconstruction of continuous traits in general, starting wih our Brownian model.
First, let's simulate some data:
## simulate a tree & some data
tree<-pbtree(n=26,scale=1,tip.label=LETTERS[26:1])
## simulate with ancestral states
x<-fastBM(tree,internal=TRUE)
## ancestral states
a<-x[as.character(1:tree$Nnode+Ntip(tree))]
## tip data
x<-x[tree$tip.label]
Now, let's estimate ancestral states using fastAnc which uses the re-rooting
algorithm discussed in class:
fit<-fastAnc(tree,x,CI=TRUE)
print(fit,printlen=6)
## Ancestral character estimates using fastAnc:
## 27 28 29 30 31 32
## 0.165976 -0.347447 -0.618424 0.493304 0.688705 0.785969 ....
##
## Lower & upper 95% CIs:
## lower upper
## 27 -0.963468 1.295421
## 28 -0.851736 0.156843
## 29 -0.724182 -0.512666
## 30 -0.063819 1.050428
## 31 0.252114 1.125296
## 32 0.531745 1.040193
## .... ....
We can easily compare these estimates to the (normally unknown) generating values for simulated data:
plot(a,fit$ace,xlab="true states",ylab="estimated states",bg="grey",cex=1.4,
pch=21)
lines(range(c(x,a)),range(c(x,a)),lty="dashed",col="red") ## 1:1 line
One of the most common critiques of ancestral state estimation is that the uncertainty around ancestral states can be large; and, furthermore, that this uncertainty is often ignored. Let's address this by obtaining 95% CIs on ancestral values, and then let's show that the 95% CIs include the generating values around 95% of the time:
## first, let's go back to our previous dataset
print(fit)
## Ancestral character estimates using fastAnc:
## 27 28 29 30 31 32 33
## 0.165976 -0.347447 -0.618424 0.493304 0.688705 0.785969 0.822307
## 34 35 36 37 38 39 40
## 0.493130 0.657181 0.678390 0.661762 0.688963 0.247025 -0.073595
## 41 42 43 44 45 46 47
## 0.602520 1.455541 1.456606 1.318106 1.325461 0.408030 -0.495230
## 48 49 50 51
## -0.494141 -0.497957 0.587861 0.704821
##
## Lower & upper 95% CIs:
## lower upper
## 27 -0.963468 1.295421
## 28 -0.851736 0.156843
## 29 -0.724182 -0.512666
## 30 -0.063819 1.050428
## 31 0.252114 1.125296
## 32 0.531745 1.040193
## 33 0.602556 1.042057
## 34 -0.015224 1.001484
## 35 0.226618 1.087745
## 36 0.281925 1.074855
## 37 0.230731 1.092793
## 38 0.289397 1.088530
## 39 -0.151341 0.645392
## 40 -0.422041 0.274851
## 41 0.488405 0.716635
## 42 1.124685 1.786398
## 43 1.127237 1.785975
## 44 1.050908 1.585304
## 45 1.155108 1.495814
## 46 -0.123436 0.939497
## 47 -0.618476 -0.371985
## 48 -0.559153 -0.429129
## 49 -0.542691 -0.453224
## 50 0.091429 1.084293
## 51 0.623005 0.786638
mean(((a>=fit$CI95[,1]) + (a<=fit$CI95[,2]))==2)
## [1] 0.96
One small tree doesn't tell us much, so let's repeat for a sample of trees & reconstructions:
## custom function that conducts a simulation, estimates ancestral
## states, & returns the fraction on 95% CI
foo<-function(){
tree<-pbtree(n=100)
x<-fastBM(tree,internal=TRUE)
fit<-fastAnc(tree,x[1:length(tree$tip.label)],CI=TRUE)
mean(((x[1:tree$Nnode+length(tree$tip.label)]>=fit$CI95[,1]) +
(x[1:tree$Nnode+length(tree$tip.label)]<=fit$CI95[,2]))==2)
}
## conduct 100 simulations
pp<-replicate(100,foo())
mean(pp)
## [1] 0.9479798
This shows us that although CIs can be large, when the model is correct they will include the generating value (1-α)% of the time.
Discrete characters
Next, we can examine ancestral state reconstruction for discrete character traits.
Most commonly, the estimation of ancestral character states for discretely valued traits uses a continuous-time Markov chain model commonly known as the Mk model.
We will use a tree & dataset of feeding mode in elopomorph eels. The data & tree are in two files as follows:
Now, let's read in the data & tree:
## read data
X<-read.csv("elopomorph.csv",row.names=1)
feed.mode<-setNames(X[,1],rownames(X))
feed.mode
## Albula_vulpes Anguilla_anguilla
## suction suction
## Anguilla_bicolor Anguilla_japonica
## suction suction
## Anguilla_rostrata Ariosoma_anago
## suction suction
## Ariosoma_balearicum Ariosoma_shiroanago
## suction suction
## Bathyuroconger_vicinus Brachysomophis_crocodilinus
## bite bite
## Conger_japonicus Conger_myriaster
## suction suction
## Conger_verreaxi Conger_wilsoni
## suction suction
## Congresox_talabonoides Cynoponticus_ferox
## suction bite
## Dysomma_anguillare Elops_saurus
## bite suction
## Facciolella_gilbertii Gavialiceps_taeniola
## bite bite
## Gnathophis_longicauda Gorgasia_taiwanensis
## suction suction
## Gymnothorax_castaneus Gymnothorax_flavimarginatus
## bite bite
## Gymnothorax_kidako Gymnothorax_moringa
## bite bite
## Gymnothorax_pseudothyrsoideus Gymnothorax_reticularis
## bite bite
## Heteroconger_hassi Ichthyapus_ophioneus
## suction bite
## Kaupichthys_hyoproroides Kaupichthys_nuchalis
## bite bite
## Megalops_cyprinoides Moringua_edwardsi
## suction bite
## Moringua_javanica Muraenesox_bagio
## bite bite
## Muraenesox_cinereus Myrichthys_breviceps
## bite suction
## Myrichthys_maculosus Myrichthys_magnificus
## suction suction
## Myrophis_vafer Nemichthys_scolopaceus
## bite bite
## Nettastoma_melanurum Ophichthus_serpentinus
## bite suction
## Ophichthus_zophochir Oxyconger_leptognathus
## suction bite
## Parabathymyrus_macrophthalmus Paraconger_notialis
## suction suction
## Pisodonophis_cancrivorus Poeciloconger_kapala
## bite suction
## Rhinomuraena_quaesita Rhynchoconger_flavus
## bite suction
## Saurenchelys_fierasfer Scolecenchelys_breviceps
## bite suction
## Scuticaria_tigrina Serrivomer_beanii
## bite bite
## Serrivomer_sector Simenchelys_parasitica
## bite suction
## Uroconger_lepturus Uropterygius_micropterus
## suction bite
## Venefica_proboscidea
## bite
## Levels: bite suction
## read tree
eel.tree<-read.tree("elopomorph.tre")
eel.tree
##
## Phylogenetic tree with 61 tips and 60 internal nodes.
##
## Tip labels:
## Moringua_edwardsi, Kaupichthys_nuchalis, Gorgasia_taiwanensis, Heteroconger_hassi, Venefica_proboscidea, Anguilla_rostrata, ...
##
## Rooted; includes branch lengths.
In visual format, these are the data that we have:
plotTree(eel.tree,type="fan",fsize=0.7,ftype="i",lwd=1)
cols<-setNames(c("red","blue"),levels(feed.mode))
tiplabels(pie=to.matrix(feed.mode[eel.tree$tip.label],
levels(feed.mode)),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=0.8*par()$usr[3],fsize=0.8)
Marginal ancestral state reconstruction
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(feed.mode,eel.tree,model="ER",type="discrete")
fitER
##
## Ancestral Character Estimation
##
## Call: ace(x = feed.mode, phy = eel.tree, type = "discrete", model = "ER")
##
## Log-likelihood: -36.33992
##
## Rate index matrix:
## bite suction
## bite . 1
## suction 1 .
##
## Parameter estimates:
## rate index estimate std-err
## 1 0.0158 0.0045
##
## Scaled likelihoods at the root (type '...$lik.anc' to get them for all nodes):
## bite suction
## 0.5480037 0.4519963
fitER$lik.anc
## bite suction
## [1,] 0.5480037210 0.4519962790
## [2,] 0.6242342545 0.3757657455
## [3,] 0.6560898907 0.3439101093
## [4,] 0.6904137188 0.3095862812
## [5,] 0.7066485827 0.2933514173
## [6,] 0.7158097958 0.2841902042
## [7,] 0.7151789283 0.2848210717
## [8,] 0.0005280939 0.9994719061
## [9,] 0.7709017973 0.2290982027
## [10,] 0.7691863614 0.2308136386
## [11,] 0.6440479559 0.3559520441
## [12,] 0.0442832132 0.9557167868
## [13,] 0.0037653340 0.9962346660
## [14,] 0.0172004443 0.9827995557
## [15,] 0.8217944284 0.1782055716
## [16,] 0.8569460714 0.1430539286
## [17,] 0.7393039669 0.2606960331
## [18,] 0.7174480970 0.2825519030
## [19,] 0.6974032138 0.3025967862
## [20,] 0.7051615816 0.2948384184
## [21,] 0.7765900337 0.2234099663
## [22,] 0.9679563345 0.0320436655
## [23,] 0.7763459919 0.2236540081
## [24,] 0.7355606936 0.2644393064
## [25,] 0.5593038528 0.4406961472
## [26,] 0.5493603890 0.4506396110
## [27,] 0.6441467852 0.3558532148
## [28,] 0.5530135261 0.4469864739
## [29,] 0.2983895523 0.7016104477
## [30,] 0.2715130074 0.7284869926
## [31,] 0.0184977763 0.9815022237
## [32,] 0.0425564134 0.9574435866
## [33,] 0.0093509902 0.9906490098
## [34,] 0.5714981738 0.4285018262
## [35,] 0.4984340090 0.5015659910
## [36,] 0.3238194112 0.6761805888
## [37,] 0.0317279739 0.9682720261
## [38,] 0.0115367770 0.9884632230
## [39,] 0.3773398135 0.6226601865
## [40,] 0.5502264549 0.4497735451
## [41,] 0.7532755531 0.2467244469
## [42,] 0.7565228047 0.2434771953
## [43,] 0.8568317500 0.1431682500
## [44,] 0.9665651556 0.0334348444
## [45,] 0.9844171015 0.0155828985
## [46,] 0.9800267702 0.0199732298
## [47,] 0.9939408074 0.0060591926
## [48,] 0.9993040631 0.0006959369
## [49,] 0.9994145342 0.0005854658
## [50,] 0.9990235521 0.0009764479
## [51,] 0.9803292025 0.0196707975
## [52,] 0.8689467273 0.1310532727
## [53,] 0.8126122235 0.1873877765
## [54,] 0.9847413814 0.0152586186
## [55,] 0.2057221815 0.7942778185
## [56,] 0.0438289116 0.9561710884
## [57,] 0.0454895580 0.9545104420
## [58,] 0.0032374280 0.9967625720
## [59,] 0.3028016114 0.6971983886
## [60,] 0.1537905871 0.8462094129
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(eel.tree,type="fan",fsize=0.7,ftype="i",lwd=1)
nodelabels(node=1:eel.tree$Nnode+Ntip(eel.tree),
pie=fitER$lik.anc,piecol=cols,cex=0.4)
tiplabels(pie=to.matrix(feed.mode[eel.tree$tip.label],
levels(feed.mode)),piecol=cols,cex=0.3)
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=0.8*par()$usr[3],fsize=0.8)
We can do the same with any model for the evolution of our traits on the tree.
Stochastic character mapping
An alternative tactic 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(eel.tree,feed.mode,model="ER")
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## bite suction
## bite -0.01582783 0.01582783
## suction 0.01582783 -0.01582783
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## bite suction
## 0.5 0.5
## Done.
plot(mtree,cols,type="fan",fsize=0.7,ftype="i")
add.simmap.legend(colors=cols,prompt=FALSE,x=0.9*par()$usr[1],
y=0.8*par()$usr[3],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(eel.tree,feed.mode,model="ER",nsim=100)
## make.simmap is sampling character histories conditioned on the transition matrix
##
## Q =
## bite suction
## bite -0.01582783 0.01582783
## suction 0.01582783 -0.01582783
## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =
## bite suction
## 0.5 0.5
## 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<-summary(mtrees)
pd
## 100 trees with a mapped discrete character with states:
## bite, suction
##
## trees have 31.36 changes between states on average
##
## changes are of the following types:
## bite,suction suction,bite
## x->y 17.44 13.92
##
## mean total time spent in each state is:
## bite suction total
## raw 1136.8164746 847.2847672 1984.101
## prop 0.5729629 0.4270371 1.000
plot(pd,fsize=0.6,ftype="i",colors=cols,
ylim=c(-2,Ntip(eel.tree)))
add.simmap.legend(colors=cols[2:1],prompt=FALSE,x=0,
y=-4,vertical=FALSE)
For binary discrete characters, we can also use a method in phytools called densityMap to visualize the posterior probability of being in each state across all the edges and nodes of the tree as follows:
obj<-densityMap(mtrees,states=levels(feed.mode)[2:1],plot=FALSE)
## sorry - this might take a while; please be patient
plot(obj,fsize=c(0.6,1))
Written by Liam J. Revell. Last updated 28 July 2018.