A little while ago I
posted
a function to fit multiple M*k* models to different painted 'regimes' on the tree.

This method differs from others such as the 'covarion' model of Penny et al. (2001) or the related hidden rates model of Beaulieu et al. (2013) in that it treats the regimes as 'known' - rather than integrating over uncertainty in their painting. Thus, the approach is more appropriate for comparing among a limited number of a priori hypotheses - such as a difference in rate between clades or across time periods (such as geological eras), rather than for testing for a relationship between two or more discrete traits or for trying to identify changes in the rate of character evolution for a single discrete trait, absent an a priori hypothesis regarding how the rate may have changed.

I'm trying to write a short note describing this method, but first I need a *simulator*
for the model. In a prior
post, I
simulated under a variable-rate model merely by stretching the branch lengths of the tree.
Unfortunately, this only works if **Q**_{2} is a scalar multiple of
**Q**_{1}. For more complex scenarios in which the various
**Q**_{i} we'd like to simulate are not proportional we'll need a tool that
is a tiny bit more sophisticated.

Using the *phytools* `"simmap"`

object class to encode our paintings, the
way I decided to do this was first by converting our `"simmap"`

object to
an ordinary `"phylo"`

object with singleton nodes instead of mappings along edges.
This way, we can only have one regime per edge. Then I proceeded across all the edges of
the tree with singletons and for each edge computed
**P=**expm(**Q**_{i}×*t*) in which **Q**_{i}
indicates the correct transition matrix for that edge, *t* is the length of the edge,
expm is the matrix exponential, and **P** is the matrix of transition probabilities.
Next, I proceeded to carry out a single preorder tree traversal in which every node is
visited, a character state is assigned using the state at the parent node and the computed
matrix of transition probabilities, **P**, and then the algorithm proceeds to the
daughters of the current node until all nodes have been visited.

Here is what that code looks like:

```
sim.multiMk<-function(tree,Q,anc=NULL,...){
ss<-rownames(Q[[1]])
tt<-map.to.singleton(reorder(tree))
P<-vector(mode="list",length=nrow(tt$edge))
for(i in 1:nrow(tt$edge))
P[[i]]<-expm(Q[[names(tt$edge.length)[i]]]*tt$edge.length[i])
if(is.null(anc)) anc<-sample(ss,1)
STATES<-matrix(NA,nrow(tt$edge),2)
root<-Ntip(tt)+1
STATES[which(tt$edge[,1]==root),1]<-anc
for(i in 1:nrow(tt$edge)){
new<-ss[which(rmultinom(1,1,P[[i]][STATES[i,1],])[,1]==1)]
STATES[i,2]<-new
ii<-which(tt$edge[,1]==tt$edge[i,2])
if(length(ii)>0) STATES[ii,1]<-new
}
x<-as.factor(
setNames(sapply(1:Ntip(tt),function(n,S,E) S[which(E==n)],
S=STATES[,2],E=tt$edge[,2]),tt$tip.label))
x
}
```

For completeness, here's an equivalent simulator for a single **Q** matrix:

```
sim.Mk<-function(tree,Q,anc=NULL,...){
ss<-rownames(Q)
tt<-reorder(tree)
P<-vector(mode="list",length=nrow(tt$edge))
for(i in 1:nrow(tt$edge))
P[[i]]<-expm(Q*tt$edge.length[i])
if(is.null(anc)) anc<-sample(ss,1)
STATES<-matrix(NA,nrow(tt$edge),2)
root<-Ntip(tt)+1
STATES[which(tt$edge[,1]==root),1]<-anc
for(i in 1:nrow(tt$edge)){
new<-ss[which(rmultinom(1,1,P[[i]][STATES[i,1],])[,1]==1)]
STATES[i,2]<-new
ii<-which(tt$edge[,1]==tt$edge[i,2])
if(length(ii)>0) STATES[ii,1]<-new
}
x<-as.factor(
setNames(sapply(1:Ntip(tt),function(n,S,E) S[which(E==n)],
S=STATES[,2],E=tt$edge[,2]),tt$tip.label))
x
}
```

For fun, let's imagine the following regimes:

```
library(phytools)
colors<-setNames(c("blue","red"),0:1)
plot(tree,lwd=1,colors=colors,ftype="off")
```

& then let's simulate two characters, *x* & *y* on this tree in which *x* evolved
via a variable-rate process, whereas *y* evolved under a single set of transition rates:

```
Qx<-setNames(
list(matrix(c(-0.5,0.5,0.5,-0.5),2,2,dimnames=list(letters[1:2],
letters[1:2])),
matrix(c(-5,5,5,-5),2,2,dimnames=list(letters[1:2],
letters[1:2]))),0:1)
Qx
```

```
## $`0`
## a b
## a -0.5 0.5
## b 0.5 -0.5
##
## $`1`
## a b
## a -5 5
## b 5 -5
```

```
Qy<-matrix(c(-1,1,1,-1),2,2,dimnames=list(letters[1:2],letters[1:2]))
Qy
```

```
## a b
## a -1 1
## b 1 -1
```

```
x<-sim.multiMk(tree,Qx)
plot(tree,colors,lwd=1,ftype="off",type="fan",mar=c(1.1,1.1,4.1,1.1))
par(fg="transparent")
tiplabels(pie=to.matrix(x,letters[1:2]),
piecol=RColorBrewer::brewer.pal(n=3,"PRGn")[c(1,3)],
cex=0.3)
par(fg="black")
title(main="Discrete trait simulated via a variable-rate process")
```

```
y<-sim.Mk(tree,Qy)
plot(tree,colors,lwd=1,ftype="off",type="fan",mar=c(1.1,1.1,4.1,1.1))
par(fg="transparent")
tiplabels(pie=to.matrix(y,letters[1:2]),
piecol=RColorBrewer::brewer.pal(n=3,"PRGn")[c(1,3)],
cex=0.3)
par(fg="black")
title(main="Discrete trait simulated via a constant-rate process")
```