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.
> 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/
> 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:
> https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Ffaculty.umb.edu%2Fliam.revell%2F&data=05%7C01%7Cliam.revell%40umb.edu%7Ce62e1382110f40edf43908db561fa1e3%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C638198466088187811%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=7MurBoMRR1PkMdQi3CFEwceft03vicc%2Br1Eyou7MOi8%3D&reserved=0
> <http://faculty.umb.edu/liam.revell/>
> Book: Phylogenetic Comparative Methods in
> R<https://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fpress.princeton.edu%2Fbooks%2Fphylogenetic-comparative-methods-in-r&data=05%7C01%7Cliam.revell%40umb.edu%7Ce62e1382110f40edf43908db561fa1e3%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C638198466088187811%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=ra55BZ92oPq4XN5NdQSGJVBUw1XGfvtifrO8oGQetHQ%3D&reserved=0>
> <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<https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fblog.phytools.org%2F2017%2F10%2Ffix-to-optimization-routine-for-non.html&data=05%7C01%7Cliam.revell%40umb.edu%7Ce62e1382110f40edf43908db561fa1e3%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C638198466088187811%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=pm32AvUIG0XfCVCwENXJZo94kpi39d5tDISaIjCUd78%3D&reserved=0>
> <http://blog.phytools.org/2017/10/fix-to-optimization-routine-for-non.html>
> that
> a similar issue was fixed in
> fit.bd<https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Ffit.bd%2F&data=05%7C01%7Cliam.revell%40umb.edu%7Ce62e1382110f40edf43908db561fa1e3%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C638198466088187811%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=pw6Nxmq3IelxC3SrihllOaKXxKouoaGwz0kD2AlipNg%3D&reserved=0>
> <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.orghttps://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-phylo&data=05%7C01%7Cliam.revell%40umb.edu%7Ce62e1382110f40edf43908db561fa1e3%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C638198466088187811%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=yF1EWG3xcJ8bTyE6m%2B%2BjEI3CyflxvIJ4LrdtWo3m5vQ%3D&reserved=0
> <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
> Searchable archive
> athttps://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.mail-archive.com%2Fr-sig-phylo%40r-project.org%2F&data=05%7C01%7Cliam.revell%40umb.edu%7Ce62e1382110f40edf43908db561fa1e3%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C638198466088187811%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=tQd5MlcqFPCpAzm509QtjhShaoBHz0qbK8yociaviXA%3D&reserved=0
>
> [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-phylo mailing list -
> R-sig-phylo@r-project.orghttps://nam10.safelinks.protection.outlook.com/?url=https%3A%2F%2Fstat.ethz.ch%2Fmailman%2Flistinfo%2Fr-sig-phylo&data=05%7C01%7Cliam.revell%40umb.edu%7Ce62e1382110f40edf43908db561fa1e3%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C638198466088187811%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=yF1EWG3xcJ8bTyE6m%2B%2BjEI3CyflxvIJ4LrdtWo3m5vQ%3D&reserved=0
> <https://stat.ethz.ch/mailman/listinfo/r-sig-phylo>
> Searchable archive at
> https://nam10.safelinks.protection.outlook.com/?url=http%3A%2F%2Fwww.mail-archive.com%2Fr-sig-phylo%40r-project.org%2F&data=05%7C01%7Cliam.revell%40umb.edu%7Ce62e1382110f40edf43908db561fa1e3%7Cb97188711ee94425953c1ace1373eb38%7C0%7C0%7C638198466088187811%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C3000%7C%7C%7C&sdata=tQd5MlcqFPCpAzm509QtjhShaoBHz0qbK8yociaviXA%3D&reserved=0
> <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/
