## Tuesday, May 16, 2023

### Some updates to the function `rateshift` for detecting temporal shifts in the rate of phenotypic evolution on a phylogeny

A phytools user recently reported persistent optimization errors using the `phytools::rateshift` function for a large phylogenetic tree.

`rateshift` fits a model with one or multiple temporal rate shifts on the phylogeny, but whose positions are not known or hypothesized a priori by the investigator. For more information about this method, check out Chapter 5 of my recent Princeton University Press book with Luke Harmon (which, by chance, is 50% off from May 15-26).

I’ve made two small updates (1, 2) to `rateshift` with the hope of improving optimization performance as follows. (1) I wrapped our call to the optimizer in a `try` function call which (in theory) should cause `rateshift` to repeat any optimization iteration that fails, rather than returning just quitting. This is valuable because I’ve discovered that it’s not uncommon for `rateshift` to get into trouble on the likelihood surface and find parameter space where the likelihood is not defined or can’t be computed. (2) I use `foreach` (from the foreach CRAN package) to parallelize multiple optimization iterations across computer cores.

Just to see how this all works, let’s try simulating a tree & dataset with a single, known rate shift – and then fitting our model to these data to try to identify the temporal position of that shift.

To take advantage of these new options in `rateshift`, we should update the phytools package from its GitHub page using devtools.

``````## load phytools
library(phytools)
packageVersion("phytools")
``````
``````##  '1.9.2'
``````

Starting with simulation, I’ll use `make.era.map` to paint my generating regimes on the tree with a shift and rate 75% of the total tree height above the root.

``````## simulate phylogenetic tree
tree<-pbtree(n=200,scale=100)
## simulate rate-shift 75% of the total
## tree height above the root
tt<-make.era.map(tree,c(0,75))
sig2<-setNames(c(10,1),1:2)
x<-sim.rates(tt,sig2)
``````

Now I’ll run the original optimization iterations (`rateshift` does 10 by default) in serial, but with `try` added to help ensure our optimization doesn’t fail. I’ll wrap my function call in `system.time` so that I create a record of how long this takes.

``````system.time(
fit_serial<-rateshift(tree,x,nrates=2)
)
``````
``````## Optimization progress:
## |..........|
## Done.
``````
``````##    user  system elapsed
##   97.88    1.23  531.17
``````
``````fit_serial
``````
``````## ML 2-rate model:
##       s^2(1) se(1) s^2(2)  se(2) k   logL
## value  11.67 3.324  1.151 0.1375 4 -656.3
##
## Shift point(s) between regimes (height above root):
##         1|2 se(1|2)
## value 73.81   3.381
##
## Model fit using ML.
##
## Frequency of best fit: 0.7
##
## R thinks it has found the ML solution.
``````
``````system.time(
fit_parallel<-rateshift(tree,x,nrates=2,parallel=TRUE)
)
``````
``````## Opened cluster with 10 cores.
## Running optimization iterations in parallel.
## Done.
``````
``````##    user  system elapsed
##    1.00    0.28  154.72
``````
``````fit_parallel
``````
``````## ML 2-rate model:
##       s^2(1) se(1) s^2(2)  se(2) k   logL
## value  11.67 3.324  1.151 0.1375 4 -656.3
##
## Shift point(s) between regimes (height above root):
##         1|2 se(1|2)
## value 73.81   3.383
##
## Model fit using ML.
##
## Frequency of best fit: 0.7
##
## R thinks it has found the ML solution.
``````

Hopefully we see that our shift point is close to the generating value (75) and the same across our two different optimization routines.

Lastly, let’s plot our fitted model and a likelihood surface of the rate shift, just for fun.

``````par(mfrow=c(1,2))
plot(fit_parallel,ftype="off")
par(bty="n",mar=c(5.1,4.1,2.1,1.1))
lik2<-likSurface.rateshift(tree,x,nrates=2,
density=40)
``````
``````## Computing likelihood surface for 2-rate (1 rate-shift) model....
``````
``````## Done.
``````
``````par(fg="blue")
abline(v=75,lty="dotted")
text(x=75,y=mean(par()\$usr[3:4]),labels="generating shift",
srt=90,pos=1)
`````` ``````par(fg="black")
``````

OK. I bet that worked even better than you expected.
{"ft_lang_label":"__label__en","ft_lang_prob":0.71487343,"math_prob":0.9216326,"size":3876,"snap":"2023-14-2023-23","text_gpt3_token_len":1078,"char_repetition_ratio":0.10847107,"word_repetition_ratio":0.17687075,"special_character_ratio":0.30340558,"punctuation_ratio":0.16140777,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.97794133,"pos_list":[0,1,2],"im_url_duplicate_count":[null,null,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2023-06-10T04:38:25Z\",\"WARC-Record-ID\":\"<urn:uuid:f7d59340-199d-4e49-979b-bc8afb31f3af>\",\"Content-Length\":\"555031\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:894dc89a-a748-4eba-ac62-4f5f42bf2604>\",\"WARC-Concurrent-To\":\"<urn:uuid:443dd0e7-e7a8-415a-af67-882eb8e521d3>\",\"WARC-IP-Address\":\"172.253.63.121\",\"WARC-Target-URI\":\"http://blog.phytools.org/2023/05/some-updates-to-function-rateshift-for.html\",\"WARC-Payload-Digest\":\"sha1:3M72LYRTH4MBP3K5FRAMWEESYIZKAF2I\",\"WARC-Block-Digest\":\"sha1:THMC5KLHS4VU4KJJ3Y4F5SUEVPZ6S7CM\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2023/CC-MAIN-2023-23/CC-MAIN-2023-23_segments_1685224656963.83_warc_CC-MAIN-20230610030340-20230610060340-00572.warc.gz\"}"}