Recently, a colleague contacted me about using the *phytools* function `fitmultiMk`

(described in an ‘in review’ pre-print here) to test a hypothesis that the *state* of one discrete character influenced the *rate* of evolution in a second, where the character history of the former is uncertain, rather than known.

Much as I’d like him to use our new method, I had to inform him that I thought a more *appropriate* analysis for his question would be to set up a “custom” Pagel ’94 - type model. What do I mean by that? Well, the Pagel (1994) approach is *most often* thought to be a test of correlated binary trait evolution; however, what it *actually* is is a test about whether the *state* of trait 1 influences the *rate* of trait 2 and/or *vice versa*. Depending on the specific manner in which the traits influence each others evolution, this might lead to an accumulation of certain combinations of character levels over time – in other words, an evolutionary correlation between the two characters.

On the other hand, if all we’re interested in doing (as my colleague indicated he was) is testing whether the condition of trait 1 influenced the back-and-forth *rate* of character 2, this too can be tested in much the same way as Pagel (1994) describes – just with a slightly different design matrix of our model.

For his benefit, and for anyone else who might be interested, I’m going to review how to do this here. Note that if both our traits are *binary* in nature, we can also using `phytools::fitPagel`

, which will do us the favor of automating all the following steps.

To start, and to make sure we’re doing what I’m going to claim we’re doing, let’s begin by simulating one discrete character (`x`

), along with a second (`y`

) whose *rate* of evolution depends on the state of `x`

. To do all this, I’m going to use *phytools* functions `pbtree`

(to obtain the tree), `sim.history`

(to simulate a discrete character history for `x`

), and `sim.multiMk`

(to simulate our second character, `y`

– the one whose rate is meant to depend on the level of `x`

). I’m going to call the levels of my trait `x`

`"slow"`

and `"fast"`

– though, of course, in practice we needn’t know which levels of `x`

correspond to sedentary or speedy rates of evolution for `y`

!

```
## load packages
library(phytools)
```

```
## simulate 400-taxon tree
tree<-pbtree(n=400,scale=1)
tree
```

```
##
## Phylogenetic tree with 400 tips and 399 internal nodes.
##
## Tip labels:
## t4, t5, t21, t315, t316, t67, ...
##
## Rooted; includes branch lengths.
```

```
## simulate history of x
Qx<-matrix(c(-0.5,0.5,0.5,-0.5),2,2,
dimnames=list(c("slow","fast"),c("slow","fast")))
Qx
```

```
## slow fast
## slow -0.5 0.5
## fast 0.5 -0.5
```

```
sim_tree<-sim.history(tree,Qx,anc="slow")
```

```
## Done simulation(s).
```

```
plot(sim_tree,setNames(palette()[c(2,4)],c("fast","slow")),
ftype="off",lwd=1)
```

(Remember, this is just the generating history of our rate regimes for `y`

!)

```
## get tip states of x
x<-getStates(sim_tree,"tips")
head(x)
```

```
## t4 t5 t21 t315 t316 t67
## "slow" "fast" "slow" "slow" "slow" "slow"
```

```
## now simulate second character, y,
## whose rate depends on x
Qy<-matrix(c(-1,1,0,1,-2,1,0,1,-2),3,3,
dimnames=list(letters[1:3],letters[1:3]))
QQ<-setNames(list(Qy,10*Qy),c("slow","fast"))
## this is a list of Q matrices for simulation
QQ
```

```
## $slow
## a b c
## a -1 1 0
## b 1 -2 1
## c 0 1 -2
##
## $fast
## a b c
## a -10 10 0
## b 10 -20 10
## c 0 10 -20
```

We should be able to see here that I’m about to simulate an ordered, three-state character with levels `"a"`

, `"b"`

, and `"c"`

– that will evolve with a rate of evolution that’s 10 \(\times\) higher on the red edges of our tree, from above, than it does on the blue edges.

```
y<-sim.multiMk(sim_tree,QQ)
head(y,20)
```

```
## t4 t5 t21 t315 t316 t67 t130 t131 t106 t85 t391 t392 t93 t364 t365 t33 t3 t263 t379 t380
## c b b a a b b a a a a a a b b b b c c c
## Levels: a b c
```

Now, of course if we *knew* the true history of `x`

on the tree we could go ahead and *use* `fitmultiMk`

(the very thing I started off by saying was not appropriate).

Just for fun, and only because we *know* the true history in this simulation, let’s do precisely that.

```
MODEL<-Qy
diag(MODEL)<-0
fitmultiMk(sim_tree,y,model=MODEL)
```

```
## Object of class "fitmultiMk".
##
## Fitted value of Q[fast]:
## a b c
## a -9.326254 9.326254 0.000000
## b 9.326254 -18.652509 9.326254
## c 0.000000 9.326254 -9.326254
##
## Fitted value of Q[slow]:
## a b c
## a -1.042671 1.042671 0.000000
## b 1.042671 -2.085342 1.042671
## c 0.000000 1.042671 -1.042671
##
## Fitted (or set) value of pi:
## a b c
## 0.3333333 0.3333333 0.3333333
##
## Log-likelihood: -280.170925
##
## Optimization method used was "nlminb"
```

Most likely, we can see here that `fitmultiMk`

gives very satisfying estimates of our generating values of \(\mathbf{Q}_{fast}\) and \(\mathbf{Q}_{slow}\)!

On the other hand, if we *don’t* know the real history of `x`

, as would be the case for most realistic datasets, a better approach (IMO) is to use the Pagel ’94 style approach.

To fit this kind of model we first need to create a new composite character of `x`

and `y`

as follows. Though it may not be necessary here, in general it’s a good idea to make sure that our trait factor includes *all* possible levels (i.e., combinations) of our two traits – not just the combinations that are observed in our data. The easiest way to do that in R is using the function `interaction`

as follows.

```
xy<-setNames(interaction(y,x),names(x))
head(xy)
```

```
## t4 t5 t21 t315 t316 t67
## c.slow b.fast b.slow a.slow a.slow b.slow
## Levels: a.fast b.fast c.fast a.slow b.slow c.slow
```

Only because I want to, let’s update the levels to `a+fast`

, `b+fast`

, etc. (This just looks nicer.)

```
levels(xy)<-gsub(".","+",levels(xy),fixed=TRUE)
head(xy)
```

```
## t4 t5 t21 t315 t316 t67
## c+slow b+fast b+slow a+slow a+slow b+slow
## Levels: a+fast b+fast c+fast a+slow b+slow c+slow
```

Finally, we’re ready to build our design matrix for our Pagel-type model.

This is going to look like a normal M*k* model design matrix, but in which we only permit changes (say) between `"a+fast"`

\(\leftrightarrow\) `"b+fast"`

, etc., *and* `"a+slow"`

\(\leftrightarrow\) `"b+slow"`

, etc., *and* `"a+fast"`

\(\leftrightarrow\) `"a+slow"`

, etc.

According to our hypothesis, each of these types of change, of course, will occur at different rates.

```
MODEL<-matrix(0,6,6,dimnames=list(levels(xy),levels(xy)))
MODEL
```

```
## a+fast b+fast c+fast a+slow b+slow c+slow
## a+fast 0 0 0 0 0 0
## b+fast 0 0 0 0 0 0
## c+fast 0 0 0 0 0 0
## a+slow 0 0 0 0 0 0
## b+slow 0 0 0 0 0 0
## c+slow 0 0 0 0 0 0
```

```
MODEL[1,2]<-MODEL[2,1]<-MODEL[2,3]<-MODEL[3,2]<-1
MODEL
```

```
## a+fast b+fast c+fast a+slow b+slow c+slow
## a+fast 0 1 0 0 0 0
## b+fast 1 0 1 0 0 0
## c+fast 0 1 0 0 0 0
## a+slow 0 0 0 0 0 0
## b+slow 0 0 0 0 0 0
## c+slow 0 0 0 0 0 0
```

```
MODEL[4,5]<-MODEL[5,4]<-MODEL[5,6]<-MODEL[6,5]<-2
MODEL
```

```
## a+fast b+fast c+fast a+slow b+slow c+slow
## a+fast 0 1 0 0 0 0
## b+fast 1 0 1 0 0 0
## c+fast 0 1 0 0 0 0
## a+slow 0 0 0 0 2 0
## b+slow 0 0 0 2 0 2
## c+slow 0 0 0 0 2 0
```

```
MODEL[1,4]<-MODEL[2,5]<-MODEL[3,6]<-
MODEL[4,1]<-MODEL[5,2]<-MODEL[6,3]<-3
MODEL
```

```
## a+fast b+fast c+fast a+slow b+slow c+slow
## a+fast 0 1 0 3 0 0
## b+fast 1 0 1 0 3 0
## c+fast 0 1 0 0 0 3
## a+slow 3 0 0 0 2 0
## b+slow 0 3 0 2 0 2
## c+slow 0 0 3 0 2 0
```

Perfect! I hope that it’s straightforward to see how we might adapt this design matrix for different hypotheses of evolution of our traits, or to more character levels of `x`

or `y`

.

Now let’s go ahead and fit our model!

```
fit_xy<-fitMk(tree,xy,model=MODEL)
fit_xy
```

```
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## a+fast b+fast c+fast a+slow b+slow c+slow
## a+fast -9.979769 9.473190 0.000000 0.506579 0.000000 0.000000
## b+fast 9.473190 -19.452959 9.473190 0.000000 0.506579 0.000000
## c+fast 0.000000 9.473190 -9.979769 0.000000 0.000000 0.506579
## a+slow 0.506579 0.000000 0.000000 -1.511920 1.005341 0.000000
## b+slow 0.000000 0.506579 0.000000 1.005341 -2.517261 1.005341
## c+slow 0.000000 0.000000 0.506579 0.000000 1.005341 -1.511920
##
## Fitted (or set) value of pi:
## a+fast b+fast c+fast a+slow b+slow c+slow
## 0.166667 0.166667 0.166667 0.166667 0.166667 0.166667
## due to treating the root prior as (a) flat.
##
## Log-likelihood: -393.865149
##
## Optimization method used was "nlminb"
##
## R thinks it has found the ML solution.
```

Just for fun, here’s a cool visualization of our data, tree, & model.

```
mat<-matrix(c(1,2),1,2)
layout(mat,widths=c(0.4,0.6))
plotTree(tree,ftype="off",lwd=1,xlim=c(0,1.1))
col1<-palette()[c(2,4)]
col2<-viridisLite::viridis(n=3)
xx<-as.factor(x[tree$tip.label])
yy<-as.factor(y[tree$tip.label])
for(i in 1:Ntip(tree)){
polygon(x=c(1,1.05,1.05,1),y=i+c(-0.5,-0.5,0.5,0.5),
col=col1[as.numeric(xx)[i]],border="transparent")
polygon(x=c(1.05,1.1,1.1,1.05),y=i+c(-0.5,-0.5,0.5,0.5),
col=col2[as.numeric(yy)[i]],border="transparent")
}
plot(fit_xy,width=TRUE,color=TRUE,show.zeros=FALSE,
mar=rep(0,4),spacer=0.15)
```

Really? Is that all?

How about comparison to a null model in which the rate of `y`

*does not* depend on `x`

?

Well, that’s easy! Just as for the standard Pagel ’94 model, the likelihood of our null (i.e., independent) model is just the *sum* of the log-likelihoods of the corresponding M*k* models for our two traits.

```
lik.x<-logLik(fitMk(tree,x,model="ER"))
lik.y<-logLik(fitMk(tree,y,
model=matrix(c(0,1,0,1,0,1,0,1,0),3,3)))
logL.null<-lik.x+lik.y
attr(logL.null,"df")<-attr(lik.x,"df")+attr(lik.y,"df")
logL.null
```

```
## [1] -420.103
## attr(,"df")
## [1] 2
```

```
## custom function for likelihood-ratio test
lr.test<-function(lik1,lik2){
LR<-2*(lik2-lik1)
as.numeric(
pchisq(LR,df=attr(lik2,"df")-attr(lik1,"df"),
lower.tail=FALSE))
}
P<-lr.test(logL.null,logLik(fit_xy))
P
```

```
## [1] 4.356051e-13
```

That’s it!