Integrating stochastic character maps across multiple character transition models

I recently fielded an interesting question by Oscar Inostroza from the Universidad de Concepción how to choose among alternative models for the substitution process for stochastic mapping when different models have comparable support. Presently, it is possible to condition on a particular model (such as, for example, model="ER", the “equal rates” model), then either fix the transition matrix Q at it's empirical ML value; or sample Q from its posterior distribution conditioned on the selected model. It is not possible, however, to sample character histories conditioned on multiple alternative models for character change between states.

Ideally what we might like to do in this case is design a reversible-jump MCMC to sample among models for character transition in proportion to their posterior probabilities; however unfortunately this is not practical at the current time. What I suggested instead is that we use Akaike weights to generate stochastic maps under all of the alternative models under consideration, in proportion to (or, approximately in proportion to - since we must use a finite number of simulations) the weight of evidence in support of that model.

There are a few tricks to doing this. In the following, I'll illustrate how it can be done using phytools.

First, some preliminaries. Let's start by loading phytools & simulating some data for the present case. I will simulate data under a "SYM" (symmetric) transition model, but in which the rates are not too dissimilar between states - which I hope will create fairly even support for the "SYM" & "ER" models on a modest sized phylogeny.

## load packages
library(phytools)

## simulate some data
tree<-pbtree(n=40,scale=1)
Q<-matrix(c(-1,0.75,0.25,
        0.75,-1.25,0.5,
        0.25,0.5,-0.75),3,3)
rownames(Q)<-colnames(Q)<-letters[1:3]
Q

##       a     b     c
## a -1.00  0.75  0.25
## b  0.75 -1.25  0.50
## c  0.25  0.50 -0.75

x<-sim.history(tree,Q)$states

## Done simulation(s).

x

## t24 t33 t34 t10  t9 t25 t26 t21 t22  t1 t39 t40 t20 t28 t29 t31 t32 t19 
## "c" "c" "c" "c" "c" "a" "a" "b" "c" "b" "b" "b" "b" "c" "b" "b" "b" "b" 
## t13 t23 t30 t35 t36 t27 t14 t15  t3  t6  t7  t2 t17 t18  t8  t4  t5 t37 
## "b" "a" "a" "a" "a" "a" "a" "a" "a" "c" "a" "a" "a" "a" "a" "a" "a" "c" 
## t38 t16 t11 t12 
## "c" "b" "a" "a"

Next, here are some simple functions that we'll use later to compute the AIC scores and the Akaike weights for each model:

aic<-function(logL,k) 2*k-2*logL
aic.w<-function(aic){
    d.aic<-aic-min(aic)
    exp(-1/2*d.aic)/sum(exp(-1/2*d.aic))
}

Now, we need to fit & (more importantly) obtain the log-likelihood of each fitted model. We could actually do this with multiple functions in R (ace, fitDiscrete, & optim.pml, among others); however here I will use make.simmap, which actually computes the likelihood using code modified from ace.

## compute log likelihoods
logL<-sapply(c("ER","SYM","ARD"),
    function(model,tree,x) make.simmap(tree,x,model)$logL,
    tree=tree,x=x)

## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##            a          b          c
## a -0.7389255  0.3694628  0.3694628
## b  0.3694628 -0.7389255  0.3694628
## c  0.3694628  0.3694628 -0.7389255

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##            a         b          c
## a -0.5566061  0.000000  0.5566061
## b  0.0000000 -1.158184  1.1581840
## c  0.5566061  1.158184 -1.7147901

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##            a           b          c
## a -0.5197351  0.03634396  0.4833911
## b  0.0000000 -0.65052948  0.6505295
## c  1.4205225  0.95975112 -2.3802736

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.

logL

##        ER       SYM       ARD 
## -28.55567 -26.84571 -26.38334

## now let's compute AIC values & weights
## (in a real study, we might use AICc)
AIC<-mapply(aic,logL,c(1,3,6))
AIC

##       ER      SYM      ARD 
## 59.11134 59.69141 64.76668

AIC.W<-aic.w(AIC)
AIC.W

##         ER        SYM        ARD 
## 0.55328515 0.41398768 0.03272717

In our next step, we can “normalize” our weights to the number of simulations we want to use for stochastic mapping. The trick here is to remember that if we just round the product of the weights (which are currently normalized to sum to one) × the desired total number of simulations, we might end up with a total number of simulations less than or greater than our desired number. In the code below I (arbitrarily) add (or subtract) the deficit (or surplus) to models chosen at random; but this need not be our strategy.

## now let's normalize this to our number of stochastic
## mapping simulations
nsim<-1000
Nsim<-round(nsim*AIC.W)
d<-if(sum(Nsim)>nsim) -1 else 1
nsim<-Nsim+d*sample(c(rep(1,abs(nsim-sum(Nsim))),
    rep(0,length(Nsim)-abs(nsim-sum(Nsim)))))
nsim

##  ER SYM ARD 
## 553 414  33

Finally, I will perform nsim stochastic mapping simulations for each model & then combine stochastic maps conducted across all of our models into one object of class "multiPhylo":

## remove any with nsim==0
nsim<-nsim[nsim!=0]
trees<-list()
class(trees)<-"multiPhylo"
for(i in 1:length(nsim)){
    obj<-make.simmap(tree,x,model=names(nsim)[i],nsim=nsim[i])
    ## we could also use:
    # obj<-make.simmap(tree,x,model=names(nsim)[i],nsim=nsim[i],Q="mcmc")
    ## but it would be slower
    if(nsim[i]==1){ 
        obj<-list(obj)
        class(obj)<-"multiPhylo"
    }
    trees<-c(trees,obj)
}

## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##            a          b          c
## a -0.7389255  0.3694628  0.3694628
## b  0.3694628 -0.7389255  0.3694628
## c  0.3694628  0.3694628 -0.7389255

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##            a         b          c
## a -0.5566061  0.000000  0.5566061
## b  0.0000000 -1.158184  1.1581840
## c  0.5566061  1.158184 -1.7147901

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.
## make.simmap is sampling character histories conditioned on the transition matrix
## Q =

##            a           b          c
## a -0.5197351  0.03634396  0.4833911
## b  0.0000000 -0.65052948  0.6505295
## c  1.4205225  0.95975112 -2.3802736

## (estimated using likelihood);
## and (mean) root node prior probabilities
## pi =

##         a         b         c 
## 0.3333333 0.3333333 0.3333333

## Done.

With this object, it is possible to do any of the other standard type of things, such as use describe.simmap to summarize the results of our analysis, including posterior probabilities at nodes, etc. For instance:

obj<-describe.simmap(trees)
obj

## 1000 trees with a mapped discrete character with states:
##  a, b, c 
## 
## trees have 11.26 changes between states on average
## 
## changes are of the following types:
##        a,b   a,c   b,a   b,c   c,a   c,b
## x->y 1.238 3.593 0.254 1.816 1.667 2.692
## 
## mean total time spent in each state is:
##              a         b         c    total
## raw  8.5813823 2.6669467 2.6141229 13.86245
## prop 0.6190378 0.1923864 0.1885758  1.00000

plot(obj)

and so on.

That's all for now on this topic.