I just added a new option in `fitMk`

(for fitting the extended M*k* model of discrete character evolution on a phylogeny) to use a different internal function to compute the likelihood.

Even though it seems (pretty clearly) to be faster, I haven’t updated the default yet just in case there are any bugs that I haven’t identified or anticipated.

This isn’t a new method or new optimization routine – just a different implementation of *exactly the same* calculation that `fitMk`

was doing before.

We can adjust which likelihood function to use via the hidden function argument `lik.func`

. If `lik.func`

is set to `lik.func="lik"`

(the default), then `fitMk`

will use the previous method. If `lik.func="pruning"`

, then `fitMk`

will use the new likelihood calculation implementation. (Both likelihood calculations use the pruning algorithm of Felsenstein! This is just an indication of the two names of the different internally used functions.)

Let’s see a quick demo that shows us the difference.

I’ll analyze primate diel activity pattern evolution on a phylogeny using data from Kirk & Kay (2004).

Start by loading *phytools* plus our data & tree.

```
library(phytools)
data(primate.tree)
data(primate.data)
diel_pattern<-setNames(
as.factor(primate.data$Activity_pattern),
rownames(primate.data))
```

We’ll use `system.time`

to measure the time that each option takes to run. Remember – this is *not* a change in the optimization routine, and we’ll use the same starting values for optimization, so any difference in time should largely be due to a difference in how long it takes to compute the likelihood.

First with the default function.

```
system.time(
ard_diel.1<-fitMk(primate.tree,diel_pattern,
model="ARD",pi="fitzjohn")
)
```

```
## user system elapsed
## 4.42 0.08 13.88
```

```
print(ard_diel.1)
```

```
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## Cathemeral Diurnal Nocturnal
## Cathemeral 0.00000 0.000000 0.000000
## Diurnal 0.00416 -0.006031 0.001871
## Nocturnal 0.00000 0.006281 -0.006281
##
## Fitted (or set) value of pi:
## Cathemeral Diurnal Nocturnal
## 0.000000 0.013462 0.986538
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -30.691752
##
## Optimization method used was "nlminb"
```

Next, with `lik.func="pruning"`

.

```
system.time(
ard_diel.2<-fitMk(primate.tree,diel_pattern,
model="ARD",pi="fitzjohn",lik.func="pruning")
)
```

```
## user system elapsed
## 2.56 0.03 8.36
```

```
print(ard_diel.2)
```

```
## Object of class "fitMk".
##
## Fitted (or set) value of Q:
## Cathemeral Diurnal Nocturnal
## Cathemeral 0.00000 0.000000 0.000000
## Diurnal 0.00416 -0.006031 0.001871
## Nocturnal 0.00000 0.006281 -0.006281
##
## Fitted (or set) value of pi:
## Cathemeral Diurnal Nocturnal
## 0.000000 0.013463 0.986537
## due to treating the root prior as (a) nuisance.
##
## Log-likelihood: -30.691752
##
## Optimization method used was "nlminb"
```

Note that I make no assertions about this being *fast* in any kind of objective way – and (in fact), I’m *certain* that other R packages are better optimized to compute the likelihood under an M*k* model than *phytools*! I just thought it was interesting that I managed to see such a performance difference out of something that I programmed initially as a demo.

Let’s plot ancestral states and our fitted model.

```
plotTree(primate.tree,type="fan",lwd=1,ftype="i",
fsize=0.7,part=0.7,offset=3)
cols<-setNames(c("#DF536B","#DAA520","black"),
levels(diel_pattern))
nodelabels(pie=ancr(ard_diel.2)$ace,piecol=cols,cex=0.4)
tiplabels(pie=ard_diel.1$data[primate.tree$tip.label,],
piecol=cols,cex=0.3)
legend("topleft",names(cols),pch=21,pt.bg=cols,pt.cex=1.5)
par(cex=0.8)
plot(ard_diel.2,width=TRUE,add=TRUE,show.zeros=FALSE,
text=FALSE,color=TRUE,xlim=c(-4,0.9),ylim=c(-1,3.5),
signif=1,cex.traits=1.2)
```