[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 the