Plotting ancestral states or posterior probabilities at nodes in which one state has a PP>0.9 (or some other value)

A phytools user asks:

“I am trying to plot the ancestral reconstruction of discrete chatacter in a tree. I fit a ARD model and stochastic mapping with make.simmap function. I would like to know how to plot the posterior probabilities of the ancestral reconstruction of a qualitative trait in a tree. But only those nodes with high support (e.g. 0.9).”

This is not too difficult. First load packages:

library(phytools)

For this demo I will use the Caribbean anole tree & ecomorph state. Note that I think this is probably not enough data to fit the "ARD" model as it has a lot of parameters, but just for demonstrative purposes:

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 
## Levels: CG GB TC TG Tr Tw

trees<-make.simmap(tree,x,model="ARD",nsim=100)

## Warning in log(comp[1:M + N]): NaNs produced

## make.simmap is sampling character histories conditioned on the transition matrix
## 
## Q =
##               CG            GB            TC            TG            Tr
## CG -4.253868e-01  6.697818e-02  0.0930139740  1.046513e-01  0.0266286973
## GB  6.324657e-06 -6.237397e-04  0.0005920676  2.033148e-07  0.0000000000
## TC  3.144188e-02  1.031776e-05 -0.0705403837  8.211499e-06  0.0389796571
## TG  4.123322e-05  9.432208e-02  0.0000000000 -9.564207e-02  0.0005619329
## Tr  2.743656e-06  3.235536e-05  0.0360525975  4.092340e-02 -0.0995735586
## Tw  0.000000e+00  3.845547e-02  0.0437195607  0.000000e+00  0.0000000000
##               Tw
## CG  1.341147e-01
## GB  2.514404e-05
## TC  1.003162e-04
## TG  7.168285e-04
## Tr  2.256246e-02
## Tw -8.217503e-02
## (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

## Warning in rstate(L[as.character(root), ] * pi/sum(L[as.character(root), :
## Some probabilities (slightly?) < 0. Setting p < 0 to zero.

## Warning in rstate(p/sum(p)): Some probabilities (slightly?) < 0. Setting p
## < 0 to zero.

## Done.

obj<-summary(trees)
plotTree(tree,ftype="i",lwd=1,fsize=0.6,
    ylim=c(-4,Ntip(tree)))
nodes<-(1:tree$Nnode+Ntip(tree))[apply(obj$ace,1,
    function(x,alpha) any(x>(1-alpha)),alpha=0.1)]
colors<-setNames(palette()[1:length(levels(x))],
    colnames(obj$ace))
nodelabels(node=nodes,pie=obj$ace[as.character(nodes),],
    cex=0.5,piecol=colors)
add.simmap.legend(x=0,y=8,colors=colors,prompt=FALSE,
    horizontal=FALSE,shape="circle",fsize=0.9)

Now, just to double check what we've done, let's plot all the posterior probabilities:

plot(obj,ftype="i",lwd=1,fsize=0.6,ylim=c(-4,Ntip(tree)),
    cex=c(0.6,0.3),colors=colors)
add.simmap.legend(x=0,y=8,colors=colors,prompt=FALSE,
    horizontal=FALSE,shape="circle",fsize=0.9)

This code would also work for fan trees, and could be modified to work with class "ace" objects or with other reconstructions. In the case of ace, there is one important caveat. We can't rely on the matrix containing marginal likelihoods having row names corresponding to node numbers. Instead, we must remember that the rows of lik.anc are in the numerical order of the nodes.

That's it.