Dear Lior. > I don't know what property of my data creates this phenomenon. Could > this be related to the very broad range of trait values (spanning two > orders of magnitude)? Based on the trait distribution, you might consider transforming your data to a log-scale.
> 1. Is there a way to let the software determine the optimal number of > rate shifts? Something like what the SURFACE package does with OU models? /rateshift/ is a likelihood method, so you can fit models with 0, 1, 2, 3, shifts (and so on), and then compare models using AIC or likelihood ratio tests. As we show in Chapter 5 of our book (https://press.princeton.edu/books/phylogenetic-comparative-methods-in-r), it's even possible to compare the /rateshift/ model to other models, such as the "EB" model of /geiger::fitContinuous/, in which the rate of evolution changes as a continuous (rather than step) function of time. > 2. Is it possible to test specific hypotheses, i.e. choose the > positions of rate shifts a-priori and calculate the likelihood? Yes, this corresponds to the method of O'Meara et al. (2006; /Evolution/) and is implemented in the /phytools/ function /brownie.lite/. In the case of a specific set of hypothesized rate shifts, one simply needs to first "paint" these on the tree using the /phytools/ function/make.era.map/, and then fit your model of interest with /brownie.lite/. This model can *also* be compared to /rateshift/ using standard methods. We assume that each /fitted/ (rather than set) rate shift would consume one degree of freedom. > 3. Can I forbid the software from placing shifts on branches leading > to clades smaller than X species? One *could* constrain /rateshift/ to not place a rate shift (say) any higher (above the root) than a certain total depth below the tips. (/rateshift/ doesn't give this control to the user, but it would be an easy feature to add.) On the other hand, I'm not sure I would advise doing this simply because a rate shift very close to the present is probably indicative of unexplained variance in the tip taxa, rather than a genuine temporal rate shift. This could be due, for example, to measurement error / uncertainty in the species means for the trait we're studying. I hope this is helpful. All the best, Liam Liam J. Revell Professor of Biology, University of Massachusetts Boston Web: http://faculty.umb.edu/liam.revell/ Book: Phylogenetic Comparative Methods in R <https://press.princeton.edu/books/phylogenetic-comparative-methods-in-r> (/Princeton University Press/, 2022) On 5/17/2023 10:37 AM, Lior Glick wrote: > CAUTION: EXTERNAL SENDER > Hi Liam, and thanks a lot for the quick fix! I highly appreciate it. > I was able to install from github and the commands now completes > successfully. Unfortunately, convergence was not reached with 100 > iterations. This is probably not surprising if the likelihood surface > is flat. I don't know what property of my data creates this > phenomenon. Could this be related to the very broad range of trait > values (spanning two orders of magnitude)? > > Some other questions related to this type of analysis: > 1. Is there a way to let the software determine the optimal number of > rate shifts? Something like what the SURFACE package does with OU models? > 2. Is it possible to test specific hypotheses, i.e. choose the > positions of rate shifts a-priori and calculate the likelihood? > 3. Can I forbid the software from placing shifts on branches leading > to clades smaller than X species? > > As far as I can tell, these features were not implemented in > rateshift(), so my questions are, whether such analyses have been > implemented somewhere else, and do you see any reason why this > shouldn't work (if I make the relevant adjustments to your code). > > Thanks again. > Lior > > On Tue, May 16, 2023 at 9:48 PM Liam J. Revell <liam.rev...@umb.edu> > wrote: > > Hi Lior. > > I just updated /rateshift/ in /phytools/ on GitHub with these > adjustments. It seems to work well and will run your dataset > without failing. > > I did observe an odd thing with the two-rate model for your data > which is that the likelihood surface for the position of the > rate-shift seems to be totally flat. I don't know what would give > your data this property. (In fact, I don't know if I could > simulate data with this attribute if I tried.) I found this result > to be so worrying that I took your original tree & simulated data > containing a single rate shift -- but I found that I was able to > recover the temporal position of this shift without an enormous > amount of trouble (though optimization /is/ very slow). This > suggests that I am not seeing the result of a software bug or an > idiosyncratic problem with the tree. > > To update /phytools/ from GitHub you just need to have the > /devtools/ CRAN package installed and then run > /devtools:::install_github("liamrevell/phytools")/ from the > command line of an interactive session in R. > > If you want to try to visualize the likelihood surface for your > rate shift, you make want to try this /phytools/ function > /likeSurface.rateshift/. See more information here: > > http://blog.phytools.org/2016/08/function-to-plot-likelihood-surface-for.html > > <http://blog.phytools.org/2016/08/function-to-plot-likelihood-surface-for.html>. > Unfortunately, it is also bound to be pretty slow for larger > phylogenetic trees. > > Let us know if this is of help. > > Sincerely, Liam > > Liam J. Revell > Professor of Biology, University of Massachusetts Boston > Web: http://faculty.umb.edu/liam.revell/ > <http://faculty.umb.edu/liam.revell/> > Book: Phylogenetic Comparative Methods in R > <https://press.princeton.edu/books/phylogenetic-comparative-methods-in-r> > (/Princeton University Press/, 2022) > > > On 5/16/2023 11:09 AM, Liam J. Revell wrote: >> CAUTION: EXTERNAL SENDER >> >> Hi Lior. >> >> I'm working on this for you now. >> >> Basically a lot of optimization iterations with /rateshift/ tend to >> fail, so I'm basically doing two things: (1) preventing a failed >> optimization from causing the whole run to fail using /try/; and (2) >> allowing the user to parallelize the multiple optimization iterations >> using /foreach/. >> >> When I get both of these fixes working I'll push the update to the >> /phytools/ GitHub and you can install the package from there. >> >> All the best, Liam >> >> Liam J. Revell >> Professor of Biology, University of Massachusetts Boston >> Web:http://faculty.umb.edu/liam.revell/ >> <http://faculty.umb.edu/liam.revell/> >> Book: Phylogenetic Comparative Methods in R >> >> <https://press.princeton.edu/books/phylogenetic-comparative-methods-in-r> >> <https://press.princeton.edu/books/phylogenetic-comparative-methods-in-r> >> (/Princeton University Press/, 2022) >> >> >> On 5/16/2023 9:57 AM, Lior Glick wrote: >>> CAUTION: EXTERNAL SENDER >>> Hello, >>> I am trying to run the rateshift function from phytools on a rather >>> large phylogeny (500 species). >>> With nrates=1 (null model) everything works fine. Using nrates=2 also >>> completed successfully, but did not reach convergence: >>>> bm_multi = rateshift(tree, data, nrates=2) >>>> bm_multi$convergence >>> 0 >>> >>> I therefore tried to increase to niter=100, but then I got the >>> following error: >>> Optimization progress: >>> |........Error in solve.default(model2$hessian) : >>> Lapack routine dgesv: system is exactly singular: U[3,3] = 0 >>> Calls: rateshift >>> Execution halted >>> >>> I also tried changing to method='REML' and got a slightly different >>> error: >>> Optimization progress: >>> >>> |............................................................................Error >>> in optimHess(sig2, lik2, tree = tree, x = x) : >>> non-finite finite-difference value [2] >>> Calls: rateshift >>> Execution halted >>> >>> I saw in an old post >>> >>> <http://blog.phytools.org/2017/10/fix-to-optimization-routine-for-non.html> >>> >>> <http://blog.phytools.org/2017/10/fix-to-optimization-routine-for-non.html> >>> that >>> a similar issue was fixed infit.bd <http://fit.bd/> >>> <http://fit.bd/> <http://fit.bd/>. >>> Is there any fix for rateshit? Alternatively, has anything similar >>> been implemented in another function/package? Or any other suggestion >>> on how to go around this? >>> >>> Thanks! >>> >>> _______________________________________________ >>> R-sig-phylo mailing list -R-sig-phylo@r-project.org >>> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo >>> <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo> >>> Searchable archive >>> athttp://www.mail-archive.com/r-sig-phylo@r-project.org/ >>> <http://www.mail-archive.com/r-sig-phylo@r-project.org/> >> [[alternative HTML version deleted]] >> >> _______________________________________________ >> R-sig-phylo mailing list -R-sig-phylo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo >> <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo> >> Searchable archive >> athttp://www.mail-archive.com/r-sig-phylo@r-project.org/ >> <http://www.mail-archive.com/r-sig-phylo@r-project.org/> > [[alternative HTML version deleted]] _______________________________________________ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
