Over the past couple of days I’ve posted (1, 2) about fitting an ordered discrete character evolution (M*k*) model in which the transition rates between adjacent states vary one from the other according to a functional form. Last night I added R code for a function (along with some class methods) to fit this model, for the condition in which our functional form is a quadratic polynomial that can be different for our forward & backward transition processes.

I’ve now upgraded & added this model to *phytools* in the form of the function `fitfnMk`

(see code here).

Rather than merely being limited to simple quadratic functions, `fitfnMk`

can fit *n*th degree polynomials in which *n* is arbitrary and can be different for forward and backward transition rates. \(n = 2\) corresponds to the quadratic functions we saw before, while \(n = 1\) is a straight line and \(n = 0\) is a constant. (Likewise, \(n = 3\) is a cubic polynomial, as we’ll see below.)

OK, let’s try it.

First, load *phytools* (which will need to be updated from it’s GitHub page in order to follow along).

```
library(phytools)
packageVersion("phytools")
```

```
## [1] '2.1.16'
```

Now, I’m going to simulate ordered discrete character data on a stochastic phylogeny according to two different 3rd degree polynomial functional forms: one, the forward rates, with *high* variation (i.e., high amplitude in my cubic polynomial function) in rates among adjacent states; and a second, the backward rates, with *low* variation. (If you’re curious about where these cubic functions came from, I used this site to generate them in standard form based on a set of four points I wanted on each curve!)

```
k<-21
x<-0:(k-2)+0.5
q1<-0.0105*x^3-0.31491*x^2+2.43862*x-0.64189
q2<--0.00233*x^3+0.06998*x^2-0.54192*x+3.25375
par(mar=c(5.1,4.1,2.1,1.1))
plot(x,q1,type="b",col="blue",bty="n",las=1,
axes=FALSE,xlab="",ylab="transition rate (q)",
ylim=c(0,max(c(q1,q2))))
lines(x,q2,type="b",col="red")
labs<-mapply(function(x,y) bquote(.(x) %<->% .(y)),
x=0:(k-2),y=1:(k-1))
axis(1,at=seq(0.5,k-1.5,by=1),labels=rep("",k-1))
nulo<-mapply(mtext,text=labs,at=seq(0.5,k-1.5,by=1),
MoreArgs=list(side=1,line=1,las=3,cex=0.7))
axis(2,las=1,cex.axis=0.8)
grid()
legend("bottomleft",c("forward","backward"),
col=c("blue","red"),
lty="solid",pch=1,cex=0.8,title="generating rates")
```

As in my prior posts, I’ll simulate on a *large* tree (1,000 tips in this case), just to give us the best chance of recovering our generating process. On the other hand, these models have lower parameter complexity (specifically, \(n_{forward}+n_{backward}+2\) for our *n*th-degree polynomial in each trait direction) than do many other M*k* model designs. For example, an ordered model with idiosyncratically different substitution rates between adjacent states would have \(2k-2 = 40\) parameters to be estimated for \(k = 21\) character levels, compared to 8 for a cubic polynomial functional form in each of the backward & forward directions.

```
Q<-matrix(0,k,k,dimnames=list(0:(k-1),0:(k-1)))
Q[cbind(1:(k-1),2:k)]<-q1
Q[cbind(2:k,1:(k-1))]<-q2
diag(Q)<--rowSums(Q)
tree<-pbtree(n=1000,scale=10)
y<-sim.Mk(tree,Q,anc=as.character(20))
head(y,30)
```

```
## t707 t989 t990 t460 t461 t309 t310 t517 t977 t978 t865 t703 t704 t928 t929 t482 t378 t379 t100 t74 t404 t405
## 8 5 5 6 9 13 10 12 10 10 8 10 10 10 10 10 11 11 11 8 9 9
## t46 t47 t39 t202 t547 t961 t962 t603
## 9 13 11 13 11 9 8 8
## Levels: 10 11 12 13 14 15 16 17 18 19 20 5 6 7 8 9
```

(I started my simulation with the ancestral condition of trait level `20`

because some inspection of the **Q** matrix told me that this condition had the lowest expected relative frequency at equilibrium under our generating model.)

To start, I’m going to go ahead and fit the generating model – which uses a 3rd degree polynomial in each forward & backward direction. To do that I’ll use `fitfnMk`

and set the argument `degree = 3`

. Note that here I have supplied my character in binary matrix format (using `phytools::to.matrix`

) to ensure that all *possible* levels (between 0 and \(k = 21\)) that the trait can assume are represented in my fitted model.

```
fit_deg3<-fitfnMk(tree,to.matrix(y,0:20),trace=0,
model="polynomial",degree=3,start="smart")
```

```
## Numerically optimizing simple equal-rates ordered model
## to get better random starting values....
```

Cool. Let’s plot our fitted model & overlay the generating process as follows.

```
par(mar=c(5.1,4.1,2.1,1.1))
plot(fit_deg3)
x<-seq(0.5,(k-2)+0.5,length.out=100)
q1<-0.0105*x^3-0.31491*x^2+2.43862*x-0.64189
q2<--0.00233*x^3+0.06998*x^2-0.54192*x+3.25375
lines(x,q1,col="blue",lty="dotted")
lines(x,q2,col="red",lty="dotted")
legend("bottomright","generating",lty="dotted",cex=0.8)
```

That’s pretty neat.

For comparison, and also to illustrate how it’s done, let’s consider a 3rd degree (cubic) polynomial for our *forward* rates, and a 0th degree polynomial (i.e., a constant value) for our *backward* rates.

For some reason I don’t yet understand, this model tends to get stuck on local optima – so I’m going to try to avoid that by starting it out with good values for our model parameters based on our first, more parameter-rich fitted model, as follows.

```
start<-c(fit_deg3$par[1:4],
mean(as.Qmatrix(fit_deg3)[cbind(2:k,1:(k-1))]))
start
```

```
## [1] 0.01536061 -0.47942863 4.09660798 -5.13130610 2.03754710
```

```
fit_deg3_0<-fitfnMk(tree,to.matrix(y,0:20),trace=1,
model="polynomial",degree=c(3,0),start=start)
```

```
## 0: 1955.1073: 0.0153606 -0.479429 4.09661 -5.13131 2.03755
## 1: 1932.0547: 0.0152060 -0.479441 4.09661 -5.13131 2.03755
## 2: 1931.9107: 0.0151942 -0.479441 4.09661 -5.13131 2.03755
## 3: 1931.6715: 0.0152072 -0.479827 4.09653 -5.13132 2.03758
## 4: 1931.1978: 0.0152245 -0.479811 4.09624 -5.13137 2.03784
## 5: 1917.8966: 0.0153249 -0.475968 4.02347 -5.14463 2.10657
## 6: 1897.9140: 0.0150427 -0.471537 4.02840 -5.15449 2.39150
## 7: 1897.6184: 0.0150662 -0.471535 4.02840 -5.15449 2.39150
## 8: 1897.5859: 0.0150613 -0.471533 4.02838 -5.15447 2.39150
## 9: 1897.5834: 0.0150595 -0.471531 4.02836 -5.15446 2.39150
## 10: 1895.8596: 0.0137907 -0.432399 3.71798 -4.91417 2.35956
## 11: 1892.9882: 0.0150594 -0.465541 3.88053 -4.10995 2.66633
## 12: 1892.4594: 0.0148992 -0.456926 3.73718 -3.31105 2.72103
## 13: 1891.3925: 0.0139966 -0.417587 3.18560 -0.732851 2.80694
## 14: 1890.6960: 0.0128802 -0.371870 2.58506 1.83552 2.83813
## 15: 1890.1392: 0.0115832 -0.320794 1.94321 4.39389 2.81774
## 16: 1889.7662: 0.0102876 -0.270886 1.33144 6.71900 2.75470
## 17: 1889.7067: 0.00999686 -0.259920 1.19955 7.18899 2.71919
## 18: 1889.7018: 0.0100780 -0.263227 1.24260 7.00224 2.71232
## 19: 1889.7012: 0.0101390 -0.265739 1.27585 6.85872 2.71089
## 20: 1889.7012: 0.0101460 -0.266052 1.28037 6.83751 2.71076
## 21: 1889.7012: 0.0101458 -0.266053 1.28052 6.83622 2.71069
## 22: 1889.7012: 0.0101456 -0.266046 1.28046 6.83639 2.71067
```

(`trace = 1`

told R to print out the progress of our numerical optimization as it ran.)

We can plot this model and overlay the generating conditions as well. The backward rates (in red) should form a straight line, while the forward rates (in blue) of our fitted model should take the form of a cubic polynomial.

```
par(mar=c(5.1,4.1,2.1,1.1))
plot(fit_deg3_0)
x<-seq(0.5,(k-2)+0.5,length.out=100)
q1<-0.0105*x^3-0.31491*x^2+2.43862*x-0.64189
q2<--0.00233*x^3+0.06998*x^2-0.54192*x+3.25375
lines(x,q1,col="blue",lty="dotted")
lines(x,q2,col="red",lty="dotted")
legend("bottomright","generating",lty="dotted",cex=0.8)
```

Awesome.

Lastly, comparing between these different models is easy.

The first model we fit (`fit_deg3`

) has the other one as a special case – but consumes three additional parameters. As such we can compare them either using AIC or a likelihood-ratio test.

Let’s do the former using a generic `anova`

call.

```
anova(fit_deg3_0,fit_deg3)
```

```
## log(L) d.f. AIC weight
## fit_deg3_0 -1889.701 5 3789.402 0.91538721
## fit_deg3 -1889.082 8 3794.165 0.08461279
```

Most likely, this tells us that although our more complex model has a better (less negative) log-likelihood, the simpler model with a cubic polynomial functional form only for the forward (and not the backward) rates of evolution is a better explanation of our data, accounting for parameter complexity.

That’s all folks!

## No comments:

## Post a Comment

Note: due to the very large amount of spam, all comments are now automatically submitted for moderation.