Saturday, September 19, 2015

Update to rateshift method that results in much better optimization for models with multiple shifts

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)

plot of chunk unnamed-chunk-2

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)

plot of chunk unnamed-chunk-4

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)

plot of chunk unnamed-chunk-5

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)

plot of chunk unnamed-chunk-6

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)

plot of chunk unnamed-chunk-7

Interesting!

1 comment:

  1. This is one of the best posts i have seen in buy replica watches a long time. thanks for posting it. while there are different grades of replicas,( some better than others) it is best to read the reviews from the owners on these boards to find out replica Omega watches UK the real quality of a watch you are wanting to purchase. never listen to a seller just arbitrarily saying his product is grade 1 or AAA. he only wants your money.

    ReplyDelete