Re: [R-sig-phylo] fitContinuous-Early Burst Model
Hi Nic and list, When using Maximum Likelihood to estimate model parameters, it's often useful to define bounds on the search, either for computational efficiency (some parameters, for example a BM rate of 1000 for log body mass in kgs, are so unlikely as to not be worth looking at) or for reasons due to the effect of a particular value on our ability to actually evaluate the likelihood. Luke Harmon programmed default bounds for all the trait evolution parameters when he wrote fitContinuous. If you go into the code, you'll be able to find them here: bounds.default - matrix(c(1e-08, 20, 1e-07, 1, 1e-06, 1, 1e-05, 2.99, 1e-07, 50, -3, 0, 1e-10, 100, -100, 100), nrow = 8, ncol = 2, byrow = TRUE) rownames(bounds.default) - c(beta, lambda, kappa, delta, alpha, a, nv, mu) colnames(bounds.default) - c(min, max) bound.default # print them to the screen The EB parameter is a, and you can see that the default bounds are -3 and 0. As I mentioned in an earlier response, the magnitude of the lower bound can affect your ability to fit the model to some datasets and this is especially so for EB: too negative a value can result in a VCV that is effectively or actually singular, meaning the likelihood cannot be evaluated at that parameter and the search will fail. You need to chose an appropriate lower bound if this is the case. One way to do this is to just play around with values for the bound to see if you can get it to work. If you go this route, you'll need to try a number of values and check that your ML estimate of a is not at the bound (i.e. your bound is too high and your estimate is the bound). The other option i suggested is to compute a reasonable maximum lower bound given the age of your tree; by that I mean the depth in time from root node to the present day or, in the case of a paleo tree, the youngest tip (although as we actually want the total amount of evolutionary time elapsed over the phylogeny, you can also think of your youngest tip as existing at the present). If the root - tip distance is T, and we place a lower limit on what we will allow the magnitude of the rate at the end of the EB process to be, expressed as a fraction of the starting rate (I think i explained this incorrectly before and implied that it was the actual rate at the tip) and we call this min.rate, then a reasonable lower bound for the EB model is log(min.rate) / T. I've used 10^-5 as a min.rate as, in my experience, rates for real data almost never decay beyond this amount, although I'm sure there are cases where this happens. At any rate, this almost always results in a reasonable lower bound for the EB process. In terms of running this on sections of your tree, I assume you are pruning out clades and fitting EB to them separately? In that case you would just obtain a reasonable bound for each by computing the depth from root node to tip in each pruned clade. Hope that helps. Graham Graham Slater Peter Buck Post-Doctoral Fellow Department of Paleobiology National Museum of Natural History The Smithsonian Institution [NHB, MRC 121] P.O. Box 37012 Washington DC 20013-7012 slat...@si.edu www.eeb.ucla.edu/gslater From: Nicolas Campione [nicolas.campi...@mail.utoronto.ca] Sent: Monday, November 05, 2012 2:04 AM To: Slater, Graham Cc: David Bapst; Frank Burbrink; r-sig-phylo@r-project.org Subject: Re: [R-sig-phylo] fitContinuous-Early Burst Model Hello All, Great discussion, I'm learning a lot about the inner workings of some of these packages. So for starters, all my taxa are fossil taxa, no extant, and in terms of branch lengths the shortest branches seem to be internal edges, which can be as short as 0.15 Ma (whole tree spans over 160 Ma). I initially tried Dave's suggestion to add a small amount of time to the branches. This initially worked perfectly, but then, once I tried time scaling the tree using other functions, or tried running the analysis with a slightly different topology, the error re-occurred. So, I would like to understand how to use the 'bounds' option a little better, as it seems that this is the most appropriate solution. I tried playing around with it and it works, but it seems to give me results that are different to those that I have previously obtained. I should expand here and say that I am also trying to fit the models at various points within the tree and it seems to me that the bound values will have to be modified depending where I am. Graham, you suggested that I compute the value of X in bounds = list(a=c(X,0)). Could you elaborate on this a bit more? When you say depth, is this node-based or edge length-based? And in terms of determining the minimum rate, I'm guessing you are referring to where the evolutionary rate is the lowest in the tree, so in the case of body size say, the lowest kg/Ma
[R-sig-phylo] fitContinuous-Early Burst Model
Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] fitContinuous-Early Burst Model
Hi Nic, It's not a problem with the branch lengths, at least not the non-ultrametricity of the tree. It's probably that the lower bound on the exponential change parameter is too low, resulting in a singular phylogenetic variance covariance matrix. Try changing the lower bound as follows: eb.fit - fitContinuous(phy, data, model =EB, bounds = list(a=c(X, 0))) # where X is the lower bound you chose. You can compute an appropriate value for X if you're willing to choose a minimum rate for the end of the EB process. In this case, the lower bound would be log(min.rate) / T where T is the depth of your tree and min.rate is the ending rate. I've found a value of 10^-5 is reasonable in most cases. Graham Graham Slater Peter Buck Post-Doctoral Fellow Department of Paleobiology National Museum of Natural History The Smithsonian Institution [NHB, MRC 121] P.O. Box 37012 Washington DC 20013-7012 slat...@si.edu www.eeb.ucla.edu/gslater On Nov 4, 2012, at 11:59 AM, Nicolas Campione wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] fitContinuous-Early Burst Model
Oops...sorry about that Nicolas. I just ran a quick sim and Graham is right on the money...ultrametricity doesn't matter and setting the lower bound on the exp change param is the way to go. I must be getting senile... This mobile message may have typos. Sorry... On Nov 4, 2012, at 1:07 PM, Graham Slater gsla...@ucla.edu wrote: Hi Nic, It's not a problem with the branch lengths, at least not the non-ultrametricity of the tree. It's probably that the lower bound on the exponential change parameter is too low, resulting in a singular phylogenetic variance covariance matrix. Try changing the lower bound as follows: eb.fit - fitContinuous(phy, data, model =EB, bounds = list(a=c(X, 0))) # where X is the lower bound you chose. You can compute an appropriate value for X if you're willing to choose a minimum rate for the end of the EB process. In this case, the lower bound would be log(min.rate) / T where T is the depth of your tree and min.rate is the ending rate. I've found a value of 10^-5 is reasonable in most cases. Graham Graham Slater Peter Buck Post-Doctoral Fellow Department of Paleobiology National Museum of Natural History The Smithsonian Institution [NHB, MRC 121] P.O. Box 37012 Washington DC 20013-7012 slat...@si.edu www.eeb.ucla.edu/gslater On Nov 4, 2012, at 11:59 AM, Nicolas Campione wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] fitContinuous-Early Burst Model
Nicolas, Not certain how you are time-scaling your tree, but I often get warnings of singularity from various comparative analysis functions due to the presence of very short terminal branches. I have found that such warnings of singularity when dealing with paleo-trees can often be avoided by adding a very small amount to all the terminal branch lengths, say 0.001. This is very simple to do with a few lines of code; alternatively, the function addTermBranchLength() in the library paleotree does this automatically if a tree is passed to it. Also, I would advice you to consider multiple methods for time-scaling your tree of fossil taxa and compare results across those trees, if you are not doing so already. In practice, I find one can sometimes give very different results for some comparative analyses with different time-scaling methods on the same paleontological dataset. -Dave On Sun, Nov 4, 2012 at 10:59 AM, Nicolas Campione nicolas.campi...@mail.utoronto.ca wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- David Bapst Dept of Geophysical Sciences University of Chicago 5734 S. Ellis Chicago, IL 60637 http://home.uchicago.edu/~dwbapst/ http://cran.r-project.org/web/packages/paleotree/index.html http://home.uchicago.edu/%7Edwbapst/ [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] fitContinuous-Early Burst Model
Just to confirm that Graham's fix works. Using the following libraries (and their dependencies) I conducted the simulations: library(TreeSim) library(phytools) sim.bd.taxa(30,1,lambda=0.1,mu=0.05)[[1]][[1]]-tree fastBM(tree)-data fitContinuous(tree,data,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular ###This now works with a parameterized alpha (below)...not sure how good the estimates are (they were simulated under ###BM after all) fitContinuous(tree,data,model=EB, bounds=list(alpha=c(0,1))) Fitting EB model: $Trait1 $Trait1$lnl [1] -264.4297 $Trait1$beta [1] 0.8355526 $Trait1$a [1] 0.001813503 $Trait1$aic [1] 534.8593 $Trait1$aicc [1] 535.1042 $Trait1$k [1] 3 On Sun, Nov 4, 2012 at 2:01 PM, David Bapst dwba...@uchicago.edu wrote: Nicolas, Not certain how you are time-scaling your tree, but I often get warnings of singularity from various comparative analysis functions due to the presence of very short terminal branches. I have found that such warnings of singularity when dealing with paleo-trees can often be avoided by adding a very small amount to all the terminal branch lengths, say 0.001. This is very simple to do with a few lines of code; alternatively, the function addTermBranchLength() in the library paleotree does this automatically if a tree is passed to it. Also, I would advice you to consider multiple methods for time-scaling your tree of fossil taxa and compare results across those trees, if you are not doing so already. In practice, I find one can sometimes give very different results for some comparative analyses with different time-scaling methods on the same paleontological dataset. -Dave On Sun, Nov 4, 2012 at 10:59 AM, Nicolas Campione nicolas.campi...@mail.utoronto.ca wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- David Bapst Dept of Geophysical Sciences University of Chicago 5734 S. Ellis Chicago, IL 60637 http://home.uchicago.edu/~dwbapst/ http://cran.r-project.org/web/packages/paleotree/index.html http://home.uchicago.edu/%7Edwbapst/ [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- NEW EMAIL ADDRESS: frank.burbr...@csi.cuny.edu *** Frank T. Burbrink, Ph.D. Professor Biology Department 6S-143 2800 Victory Blvd. College of Staten Island/CUNY Staten Island, New York 10314 E-Mail:frank.burbr...@csi.cuny.edu Phone:718-982-3961 Web Page: http://scholar.library.csi.cuny.edu/~fburbrink/ *** Washington Monthlyhttp://www.washingtonmonthly.com/magazine/septemberoctober_2012/features/americas_bestbangforthebuck_co039461.php magazine ranks the College of Staten Island as one of “America’s Best-Bang-for-the-Buck Colleges” ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] fitContinuous-Early Burst Model
Yeah, I was concerned it might have been related to a zero-branch-length problem, but some simulation shows that only Graham's fix will let fitContinuous run with a simulated paleo-tree. -Dave library(paleotree) library(geiger) taxa-simFossilTaxa_SRCond(avgtaxa=30,p=0.1,q=0.1,r=0.1) rangesCont-sampleRanges(taxa,r=0.1) cladogram-taxa2cladogram(taxa) #get the true tree at the first appearance times trueTree-taxa2phylo(taxa,obs_time=rangesCont[,1]) trait-rTraitCont(trueTree) #no solution fitContinuous(trueTree,trait,model=EB) #nope #my solution fitContinuous(addTermBranchLength(trueTree),trait,model=EB) #nope #Graham's solution fitContinuous(trueTree,trait,model=EB,bounds=list(alpha=c(0,1))) #yep On Sun, Nov 4, 2012 at 1:17 PM, Frank Burbrink frank.burbr...@csi.cuny.eduwrote: Just to confirm that Graham's fix works. Using the following libraries (and their dependencies) I conducted the simulations: library(TreeSim) library(phytools) sim.bd.taxa(30,1,lambda=0.1,mu=0.05)[[1]][[1]]-tree fastBM(tree)-data fitContinuous(tree,data,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular ###This now works with a parameterized alpha (below)...not sure how good the estimates are (they were simulated under ###BM after all) fitContinuous(tree,data,model=EB, bounds=list(alpha=c(0,1))) Fitting EB model: $Trait1 $Trait1$lnl [1] -264.4297 $Trait1$beta [1] 0.8355526 $Trait1$a [1] 0.001813503 $Trait1$aic [1] 534.8593 $Trait1$aicc [1] 535.1042 $Trait1$k [1] 3 On Sun, Nov 4, 2012 at 2:01 PM, David Bapst dwba...@uchicago.edu wrote: Nicolas, Not certain how you are time-scaling your tree, but I often get warnings of singularity from various comparative analysis functions due to the presence of very short terminal branches. I have found that such warnings of singularity when dealing with paleo-trees can often be avoided by adding a very small amount to all the terminal branch lengths, say 0.001. This is very simple to do with a few lines of code; alternatively, the function addTermBranchLength() in the library paleotree does this automatically if a tree is passed to it. Also, I would advice you to consider multiple methods for time-scaling your tree of fossil taxa and compare results across those trees, if you are not doing so already. In practice, I find one can sometimes give very different results for some comparative analyses with different time-scaling methods on the same paleontological dataset. -Dave On Sun, Nov 4, 2012 at 10:59 AM, Nicolas Campione nicolas.campi...@mail.utoronto.ca wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- David Bapst Dept of Geophysical Sciences University of Chicago 5734 S. Ellis Chicago, IL 60637 http://home.uchicago.edu/~dwbapst/ http://cran.r-project.org/web/packages/paleotree/index.html http://home.uchicago.edu/%7Edwbapst/ [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- NEW EMAIL ADDRESS: frank.burbr...@csi.cuny.edu *** Frank T. Burbrink, Ph.D. Professor Biology Department 6S-143 2800 Victory Blvd. College of Staten Island/CUNY Staten Island, New York 10314 E-Mail:frank.burbr...@csi.cuny.edu Phone:718-982-3961 Web Page: http://scholar.library.csi.cuny.edu/~fburbrink/ *** Washington Monthly
Re: [R-sig-phylo] fitContinuous-Early Burst Model
Just to add one point, if you do run this: fitContinuous(trueTree,trait,model=EB,bounds=list(alpha=c(0,1))) you're actually not fitting an early burst at all - bounding the optimization of alpha with 0 at the lower end restricts you to an accelerating, or late burst model of evolution. To be EB, you need fitContinuous(trueTree,trait,model=EB,bounds=list(alpha=c(-1,0))) and if you had fitContinuous(trueTree,trait,model=EB,bounds=list(alpha=c(-1,1))) or some equivalent that worked, this would be the AC/DC (accelerating/decelerating) model of Blomberg et al. Graham Quoting David Bapst dwba...@uchicago.edu: Yeah, I was concerned it might have been related to a zero-branch-length problem, but some simulation shows that only Graham's fix will let fitContinuous run with a simulated paleo-tree. -Dave library(paleotree) library(geiger) taxa-simFossilTaxa_SRCond(avgtaxa=30,p=0.1,q=0.1,r=0.1) rangesCont-sampleRanges(taxa,r=0.1) cladogram-taxa2cladogram(taxa) #get the true tree at the first appearance times trueTree-taxa2phylo(taxa,obs_time=rangesCont[,1]) trait-rTraitCont(trueTree) #no solution fitContinuous(trueTree,trait,model=EB) #nope #my solution fitContinuous(addTermBranchLength(trueTree),trait,model=EB) #nope #Graham's solution fitContinuous(trueTree,trait,model=EB,bounds=list(alpha=c(0,1))) #yep On Sun, Nov 4, 2012 at 1:17 PM, Frank Burbrink frank.burbr...@csi.cuny.eduwrote: Just to confirm that Graham's fix works. Using the following libraries (and their dependencies) I conducted the simulations: library(TreeSim) library(phytools) sim.bd.taxa(30,1,lambda=0.1,mu=0.05)[[1]][[1]]-tree fastBM(tree)-data fitContinuous(tree,data,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular ###This now works with a parameterized alpha (below)...not sure how good the estimates are (they were simulated under ###BM after all) fitContinuous(tree,data,model=EB, bounds=list(alpha=c(0,1))) Fitting EB model: $Trait1 $Trait1$lnl [1] -264.4297 $Trait1$beta [1] 0.8355526 $Trait1$a [1] 0.001813503 $Trait1$aic [1] 534.8593 $Trait1$aicc [1] 535.1042 $Trait1$k [1] 3 On Sun, Nov 4, 2012 at 2:01 PM, David Bapst dwba...@uchicago.edu wrote: Nicolas, Not certain how you are time-scaling your tree, but I often get warnings of singularity from various comparative analysis functions due to the presence of very short terminal branches. I have found that such warnings of singularity when dealing with paleo-trees can often be avoided by adding a very small amount to all the terminal branch lengths, say 0.001. This is very simple to do with a few lines of code; alternatively, the function addTermBranchLength() in the library paleotree does this automatically if a tree is passed to it. Also, I would advice you to consider multiple methods for time-scaling your tree of fossil taxa and compare results across those trees, if you are not doing so already. In practice, I find one can sometimes give very different results for some comparative analyses with different time-scaling methods on the same paleontological dataset. -Dave On Sun, Nov 4, 2012 at 10:59 AM, Nicolas Campione nicolas.campi...@mail.utoronto.ca wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- David Bapst Dept of Geophysical Sciences University of Chicago 5734 S. Ellis Chicago, IL 60637 http://home.uchicago.edu/~dwbapst/ http://cran.r-project.org/web/packages/paleotree/index.html http://home.uchicago.edu/%7Edwbapst/ [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org
Re: [R-sig-phylo] fitContinuous-Early Burst Model
Good point. They should all work with these data though. This mobile message may have typos. Sorry... On Nov 4, 2012, at 3:10 PM, gsla...@ucla.edu wrote: Just to add one point, if you do run this: fitContinuous(trueTree,trait,model=EB,bounds=list(alpha=c(0,1))) you're actually not fitting an early burst at all - bounding the optimization of alpha with 0 at the lower end restricts you to an accelerating, or late burst model of evolution. To be EB, you need fitContinuous(trueTree,trait,model=EB,bounds=list(alpha=c(-1,0))) and if you had fitContinuous(trueTree,trait,model=EB,bounds=list(alpha=c(-1,1))) or some equivalent that worked, this would be the AC/DC (accelerating/decelerating) model of Blomberg et al. Graham Quoting David Bapst dwba...@uchicago.edu: Yeah, I was concerned it might have been related to a zero-branch-length problem, but some simulation shows that only Graham's fix will let fitContinuous run with a simulated paleo-tree. -Dave library(paleotree) library(geiger) taxa-simFossilTaxa_SRCond(avgtaxa=30,p=0.1,q=0.1,r=0.1) rangesCont-sampleRanges(taxa,r=0.1) cladogram-taxa2cladogram(taxa) #get the true tree at the first appearance times trueTree-taxa2phylo(taxa,obs_time=rangesCont[,1]) trait-rTraitCont(trueTree) #no solution fitContinuous(trueTree,trait,model=EB) #nope #my solution fitContinuous(addTermBranchLength(trueTree),trait,model=EB) #nope #Graham's solution fitContinuous(trueTree,trait,model=EB,bounds=list(alpha=c(0,1))) #yep On Sun, Nov 4, 2012 at 1:17 PM, Frank Burbrink frank.burbr...@csi.cuny.eduwrote: Just to confirm that Graham's fix works. Using the following libraries (and their dependencies) I conducted the simulations: library(TreeSim) library(phytools) sim.bd.taxa(30,1,lambda=0.1,mu=0.05)[[1]][[1]]-tree fastBM(tree)-data fitContinuous(tree,data,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular ###This now works with a parameterized alpha (below)...not sure how good the estimates are (they were simulated under ###BM after all) fitContinuous(tree,data,model=EB, bounds=list(alpha=c(0,1))) Fitting EB model: $Trait1 $Trait1$lnl [1] -264.4297 $Trait1$beta [1] 0.8355526 $Trait1$a [1] 0.001813503 $Trait1$aic [1] 534.8593 $Trait1$aicc [1] 535.1042 $Trait1$k [1] 3 On Sun, Nov 4, 2012 at 2:01 PM, David Bapst dwba...@uchicago.edu wrote: Nicolas, Not certain how you are time-scaling your tree, but I often get warnings of singularity from various comparative analysis functions due to the presence of very short terminal branches. I have found that such warnings of singularity when dealing with paleo-trees can often be avoided by adding a very small amount to all the terminal branch lengths, say 0.001. This is very simple to do with a few lines of code; alternatively, the function addTermBranchLength() in the library paleotree does this automatically if a tree is passed to it. Also, I would advice you to consider multiple methods for time-scaling your tree of fossil taxa and compare results across those trees, if you are not doing so already. In practice, I find one can sometimes give very different results for some comparative analyses with different time-scaling methods on the same paleontological dataset. -Dave On Sun, Nov 4, 2012 at 10:59 AM, Nicolas Campione nicolas.campi...@mail.utoronto.ca wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- David Bapst Dept of Geophysical Sciences University of Chicago 5734 S. Ellis
Re: [R-sig-phylo] fitContinuous-Early Burst Model
Just as a followup to my example (which contains numerous extinct taxa). Using these values for alpha fails: fitContinuous(tree,data,model=EB,bounds=list(alpha=c(-1,0))) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular Above the lower bound of -0.4 (e.g., -0.3, -0.2 etc.) Works here: fitContinuous(tree,data,model=EB,bounds=list(alpha=c(-.3,0))) Fitting EB model: $Trait1 $Trait1$lnl [1] -264.4632 $Trait1$beta [1] 0.9233574 $Trait1$a [1] 0 $Trait1$aic [1] 534.9264 $Trait1$aicc [1] 535.1713 $Trait1$k [1] 3 On Sun, Nov 4, 2012 at 1:07 PM, Graham Slater gsla...@ucla.edu wrote: Hi Nic, It's not a problem with the branch lengths, at least not the non-ultrametricity of the tree. It's probably that the lower bound on the exponential change parameter is too low, resulting in a singular phylogenetic variance covariance matrix. Try changing the lower bound as follows: eb.fit - fitContinuous(phy, data, model =EB, bounds = list(a=c(X, 0))) # where X is the lower bound you chose. You can compute an appropriate value for X if you're willing to choose a minimum rate for the end of the EB process. In this case, the lower bound would be log(min.rate) / T where T is the depth of your tree and min.rate is the ending rate. I've found a value of 10^-5 is reasonable in most cases. Graham Graham Slater Peter Buck Post-Doctoral Fellow Department of Paleobiology National Museum of Natural History The Smithsonian Institution [NHB, MRC 121] P.O. Box 37012 Washington DC 20013-7012 slat...@si.edu www.eeb.ucla.edu/gslater On Nov 4, 2012, at 11:59 AM, Nicolas Campione wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- NEW EMAIL ADDRESS: frank.burbr...@csi.cuny.edu *** Frank T. Burbrink, Ph.D. Professor Biology Department 6S-143 2800 Victory Blvd. College of Staten Island/CUNY Staten Island, New York 10314 E-Mail:frank.burbr...@csi.cuny.edu Phone:718-982-3961 Web Page: http://scholar.library.csi.cuny.edu/~fburbrink/ *** Washington Monthlyhttp://www.washingtonmonthly.com/magazine/septemberoctober_2012/features/americas_bestbangforthebuck_co039461.php magazine ranks the College of Staten Island as one of “America’s Best-Bang-for-the-Buck Colleges” ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Re: [R-sig-phylo] fitContinuous-Early Burst Model
Yeah, my simulations with paleotree agrees with Frank's, except no values below a=0 will work at all without using addTermBranchLength; apparently the terminal ZLBs are having some effect on it. When I use addTermBranchLength, then I can get down to -0.4, at which point I get a slightly new warning message: fitContinuous(addTermBranchLength(trueTree),trait,model=EB,bounds=list(alpha=c(-0.4,0))) Fitting EB model: Error in solve.default(phyvcv) : system is computationally singular: reciprocal condition number = 2.04244e-16 On Sun, Nov 4, 2012 at 2:48 PM, Frank Burbrink frank.burbr...@csi.cuny.eduwrote: Just as a followup to my example (which contains numerous extinct taxa). Using these values for alpha fails: fitContinuous(tree,data,model=EB,bounds=list(alpha=c(-1,0))) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular Above the lower bound of -0.4 (e.g., -0.3, -0.2 etc.) Works here: fitContinuous(tree,data,model=EB,bounds=list(alpha=c(-.3,0))) Fitting EB model: $Trait1 $Trait1$lnl [1] -264.4632 $Trait1$beta [1] 0.9233574 $Trait1$a [1] 0 $Trait1$aic [1] 534.9264 $Trait1$aicc [1] 535.1713 $Trait1$k [1] 3 On Sun, Nov 4, 2012 at 1:07 PM, Graham Slater gsla...@ucla.edu wrote: Hi Nic, It's not a problem with the branch lengths, at least not the non-ultrametricity of the tree. It's probably that the lower bound on the exponential change parameter is too low, resulting in a singular phylogenetic variance covariance matrix. Try changing the lower bound as follows: eb.fit - fitContinuous(phy, data, model =EB, bounds = list(a=c(X, 0))) # where X is the lower bound you chose. You can compute an appropriate value for X if you're willing to choose a minimum rate for the end of the EB process. In this case, the lower bound would be log(min.rate) / T where T is the depth of your tree and min.rate is the ending rate. I've found a value of 10^-5 is reasonable in most cases. Graham Graham Slater Peter Buck Post-Doctoral Fellow Department of Paleobiology National Museum of Natural History The Smithsonian Institution [NHB, MRC 121] P.O. Box 37012 Washington DC 20013-7012 slat...@si.edu www.eeb.ucla.edu/gslater On Nov 4, 2012, at 11:59 AM, Nicolas Campione wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get the following error message, and am unsure as to what I can do to fix this: fitContinuous(tree,trait,model=EB) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular: U[4,4] = 0 My tree is based on fossil taxa and my branch lengths are calculated using time (in millions of years). I seem to be able to test every other model, only EB is giving me troubles. I'm pretty sure it has something to do with the branch lengths, cause if I used other measures, such as all=1 or grafen, I don't get an error. Has anybody had this issue before? Any assistance would be greatly appreciated. Cheers, -Nic- Nicolás E. Campione M.Sc., Ph.D. Candidate Dept. Ecology Evolutionary Biology University of Toronto 25 Wilcocks St. Toronto, ON Canada M5S 3B2 Royal Ontario Museum 100 Queen's Park Toronto, ON Canada M5S 2C6 Office: 416-586-5591 Email: nicolas.campi...@mail.utoronto.ca [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- NEW EMAIL ADDRESS: frank.burbr...@csi.cuny.edu *** Frank T. Burbrink, Ph.D. Professor Biology Department 6S-143 2800 Victory Blvd. College of Staten Island/CUNY Staten Island, New York 10314 E-Mail:frank.burbr...@csi.cuny.edu Phone:718-982-3961 Web Page: http://scholar.library.csi.cuny.edu/~fburbrink/ *** Washington Monthly http://www.washingtonmonthly.com/magazine/septemberoctober_2012/features/americas_bestbangforthebuck_co039461.php magazine ranks the College of Staten Island as one of Americas Best-Bang-for-the-Buck Colleges ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo -- David Bapst Dept of Geophysical Sciences University of Chicago 5734 S. Ellis Chicago, IL 60637 http://home.uchicago.edu/~dwbapst/ http://cran.r-project.org/web/packages/paleotree/index.html http://home.uchicago.edu/%7Edwbapst/ [[alternative HTML
Re: [R-sig-phylo] fitContinuous-Early Burst Model
I think it's convenient to think about how fitContinuous actually works here, or at least to think about the impact of the EB parameter on a transformed tree's branch lengths. Using simulations, try this: require(TreeSim); require(geiger); phy - sim.bd.taxa(50, 1, 0.1, 0.05, 1)[[1]][[1]]; ## simulate a 50 (extant) taxon tree with fossils phy$edge.length; # look at the branch lengths ebphy - exponentialchangeTree(phy, a = -0.4); ## now transform this tree under a strong EB process ebphy$edge.length; # take a look at these edge lengths plot(ebphy); # and plot it what you'll probably find is that the exponentialchangeTree has a lot of very very short edges and from plotting it you'll see that the only long edges are those towards the root. This is such a strong burst that it basically requires all the evolutionary change in our trait to occur in the two edges leading from the root. Such a burst is unlikely in real data. I find it helpful to think about rate half lives. An EB parameter of -0.4 gives a rate half life of log(2) / 0.4 = 1.732868 time units - that is it takes ~1.8 million years for the rate to halve from its initial value. To give an idea of what this translates to, the average mammalian order is ~ 50 my old, which would result in approximately 30 half lives elapsing. That, at least to me, has an exceedingly low prior probability! So while I agree with Dave that some paleo trees will have short internal edges or terminals that might mess up the EB model, I think what is more important is to recognize the effect that the EB parameter has on your edge lengths and to chose an appropriate lower bound when fitting the model to avoid this. Probably, for most paleo trees, an exponential change parameter of -0.4 is never going to be reasonable. Graham Graham Slater Peter Buck Post-Doctoral Fellow Department of Paleobiology National Museum of Natural History The Smithsonian Institution [NHB, MRC 121] P.O. Box 37012 Washington DC 20013-7012 slat...@si.edu www.eeb.ucla.edu/gslater From: r-sig-phylo-boun...@r-project.org [r-sig-phylo-boun...@r-project.org] On Behalf Of David Bapst [dwba...@uchicago.edu] Sent: Sunday, November 04, 2012 4:04 PM To: Frank Burbrink Cc: r-sig-phylo@r-project.org; Nicolas Campione Subject: Re: [R-sig-phylo] fitContinuous-Early Burst Model Yeah, my simulations with paleotree agrees with Frank's, except no values below a=0 will work at all without using addTermBranchLength; apparently the terminal ZLBs are having some effect on it. When I use addTermBranchLength, then I can get down to -0.4, at which point I get a slightly new warning message: fitContinuous(addTermBranchLength(trueTree),trait,model=EB,bounds=list(alpha=c(-0.4,0))) Fitting EB model: Error in solve.default(phyvcv) : system is computationally singular: reciprocal condition number = 2.04244e-16 On Sun, Nov 4, 2012 at 2:48 PM, Frank Burbrink frank.burbr...@csi.cuny.eduwrote: Just as a followup to my example (which contains numerous extinct taxa). Using these values for alpha fails: fitContinuous(tree,data,model=EB,bounds=list(alpha=c(-1,0))) Fitting EB model: Error in solve.default(phyvcv) : Lapack routine dgesv: system is exactly singular Above the lower bound of -0.4 (e.g., -0.3, -0.2 etc.) Works here: fitContinuous(tree,data,model=EB,bounds=list(alpha=c(-.3,0))) Fitting EB model: $Trait1 $Trait1$lnl [1] -264.4632 $Trait1$beta [1] 0.9233574 $Trait1$a [1] 0 $Trait1$aic [1] 534.9264 $Trait1$aicc [1] 535.1713 $Trait1$k [1] 3 On Sun, Nov 4, 2012 at 1:07 PM, Graham Slater gsla...@ucla.edu wrote: Hi Nic, It's not a problem with the branch lengths, at least not the non-ultrametricity of the tree. It's probably that the lower bound on the exponential change parameter is too low, resulting in a singular phylogenetic variance covariance matrix. Try changing the lower bound as follows: eb.fit - fitContinuous(phy, data, model =EB, bounds = list(a=c(X, 0))) # where X is the lower bound you chose. You can compute an appropriate value for X if you're willing to choose a minimum rate for the end of the EB process. In this case, the lower bound would be log(min.rate) / T where T is the depth of your tree and min.rate is the ending rate. I've found a value of 10^-5 is reasonable in most cases. Graham Graham Slater Peter Buck Post-Doctoral Fellow Department of Paleobiology National Museum of Natural History The Smithsonian Institution [NHB, MRC 121] P.O. Box 37012 Washington DC 20013-7012 slat...@si.edu www.eeb.ucla.edu/gslater On Nov 4, 2012, at 11:59 AM, Nicolas Campione wrote: Dear R Phylo List, I'm trying to test an Early Burst model of evolution against a tree using the 'fitContinuous' function in 'geiger'. However, I get
