I have been working on the phytools function `rateshift`

, for
locating one or multiple temporal shifts in the rate of evolution for a
continuous character on the tree.

This function has nice `print`

and `plot`

methods, and
would seem to be quite nice - however it suffers from one fairly critical
problem which is that it doesn't really work. That is to say, it frequently
fails to find the MLE solution.

I have put a little effort into overhauling the guts of this function now;
with the main and most critical update being that the function now uses
the phytools function `brownie.lite`

to find the MLE rates
conditioned on a given set of shift points; rather than simultaneously
optimizing rates & shift points. This effectively halves the dimensionality
of the optimization, and this seems to have had a very nice effect on
performance (although there is still the chance of being stuck on local
optima, as we'll see below).

```
library(devtools)
install_github("liamrevell/phytools",quiet=TRUE)
library(phytools)
```

```
## Loading required package: ape
## Loading required package: maps
```

Now, we can start by simulating a tree and data in which there are two shift points, and thus three Brownian rates:

```
set.seed(99)
tree<-pbtree(n=50,scale=100)
x<-sim.rates(tree<-make.era.map(tree,c(0,60,85)),
setNames(c(5,20,1),1:3))
## these are our simulated regimes:
plot(tree,setNames(c("purple","red","blue"),1:3),fsize=0.8,ftype="i",
mar=c(3.1,0.1,0.1,0.1))
axis(1)
```

Let's fit one, two, three, and four rate models to these data:

```
fit1<-rateshift(tree,x,nrates=1,niter=10)
```

```
## Optimization progress:
## |..........|
## Done.
```

```
fit1
```

```
## ML 1-rate model:
## s^2(1) se(1) k logL
## value 7.4779 1.4956 2 -206.1193
##
## This is a one-rate model.
##
## Frequency of best fit: 1
##
## R thinks it has found the ML solution.
```

```
fit2<-rateshift(tree,x,nrates=2,niter=10)
```

```
## Optimization progress:
## |..........|
## Done.
```

```
fit2
```

```
## ML 2-rate model:
## s^2(1) se(1) s^2(2) se(2) k logL
## value 13.8711 3.3579 0.5263 0.2124 4 -193.311
##
## Shift point(s) between regimes (height above root):
## 1|2 se(1|2)
## value 88.4271 0.0436
##
## Frequency of best fit: 0.9
##
## R thinks it has found the ML solution.
```

```
fit3<-rateshift(tree,x,nrates=3,niter=10)
```

```
## Optimization progress:
## |..........|
## Done.
```

```
fit3
```

```
## ML 3-rate model:
## s^2(1) se(1) s^2(2) se(2) s^2(3) se(3) k logL
## value 7e-04 NaN 16.5352 4.7994 0.5158 0.2057 6 -190.4577
##
## Shift point(s) between regimes (height above root):
## 1|2 se(1|2) 2|3 se(2|3)
## value 48.8951 10.6356 88.4271 0.0387
##
## Frequency of best fit: 0.7
##
## R thinks it has found the ML solution.
```

```
fit4<-rateshift(tree,x,nrates=4,niter=10)
```

```
## Optimization progress:
## |..........|
## Done.
```

```
fit4
```

```
## ML 4-rate model:
## s^2(1) se(1) s^2(2) se(2) s^2(3) se(3) s^2(4) se(4) k logL
## value 7e-04 NaN 833.039 NaN 5.0879 2.2729 0.5437 0.2229 8 -188.5406
##
## Shift point(s) between regimes (height above root):
## 1|2 se(1|2) 2|3 se(2|3) 3|4 se(3|4)
## value 68.2645 NaN 69.0033 NaN 90.5243 0.07
##
## Frequency of best fit: 0.2
##
## R thinks it has found the ML solution.
```

```
AIC(fit1,fit2,fit3,fit4)
```

```
## df AIC
## fit1 2 416.2387
## fit2 4 394.6221
## fit3 6 392.9153
## fit4 8 393.0813
```

We can plot our best model, in this case the generating, three-rate model:

```
plot(fit3,fsize=0.8)
```

Finally, I mentioned that it is possible to get trapped on local optima. This is true, somewhat surprisingly, even for the relatively simple, two-rate model. We can see why by visualizing the likelihood surface across the length of the tree:

```
shift<-1:99/100*max(nodeHeights(tree))
logL<-sapply(shift,function(s,tree,x) logLik(rateshift(tree,x,fixed.shift=s,
quiet=TRUE)),tree=tree,x=x)
plot(shift,logL,type="l",lwd=2)
plotTree(tree,color=rgb(0,0,1,0.25),ftype="off",mar=c(5.1,4.1,4.1,2.1),
add=TRUE)
```

We can also do this in two dimensions for the three-rate model:

```
shift<-seq(2,98,by=2)/100*max(nodeHeights(tree))
logL<-sapply(shift,function(s1,s2,tree,x)
sapply(s2,function(s2,s1,tree,x)
logLik(rateshift(tree,x,fixed.shift=if(s1!=s2) c(s1,s2) else s1,
quiet=TRUE)),s1=s1,tree=tree,x=x),s2=shift,tree=tree,x=x)
image(shift,shift,logL,col=gray.colors(60,0,1),xlab="shift 1 (or 2)",
ylab="shift 2 (or 1)")
points(fit3$shift[1],fit3$shift[2],pch=4) ## MLE
points(fit3$shift[2],fit3$shift[1],pch=4)
```

Or, alternatively:

```
library(lattice)
wireframe(logL,row.values=shift,column.values=shift,xlab="shift 1 (or 2)",
ylab="shift 2 (or 1)",colorkey=TRUE,drape=TRUE)
```

Interesting!

## No comments:

## Post a Comment