Re: [R-sig-phylo] Add terminal branches to tree
Hi Dave, Thanks a lot. It was exactly this I was trying to do. Thank you all for your advice. Cheers, Sérgio. - Mensagem original - > De: "David Bapst" <dwba...@gmail.com> > Para: "Sergio Ferreira Cardoso" <sergio.ferreira-card...@umontpellier.fr> > Cc: "Eliot Miller" <et...@cornell.edu>, "r-sig-phylo" > <r-sig-phylo@r-project.org> > Enviadas: Terça-feira, 20 De Junho de 2017 0:11:02 > Assunto: Re: [R-sig-phylo] Add terminal branches to tree > Hi Sergio, > > Sorry for being late to the party, but maybe expandTaxonTree in > paleotree does what you're looking for? I wrote it for turning trees > of genera into trees of species, with the resulting generic polytomies > collapsed or not (in case we know that the genera is paraphyletic). > > Cheers, > -Dave > > On Mon, Jun 19, 2017 at 9:36 AM, Sergio Ferreira Cardoso > <sergio.ferreira-card...@umontpellier.fr> wrote: >> Thank you all for your suggestions. I guess the more direct way is to use >> Megaptera. I tried and it worked (I added Man cra1, Man cra2, etc.). >> However, I wanted it to create a new node that would basically include all >> these taxa. I sthere any way to merge all the added species in a new node? >> Thank you all, once again. >> >> >> >> De: "Eliot Miller" <et...@cornell.edu> >> Para: "Christoph Heibl" <christoph.he...@gmx.net> >> Cc: "Liam Revell" <liam.rev...@umb.edu>, "r-sig-phylo" >> <r-sig-phylo@r-project.org>, "Sergio Ferreira Cardoso" >> <sergio.ferreira-card...@umontpellier.fr> >> Enviadas: Segunda-feira, 19 De Junho de 2017 14:40:15 >> Assunto: Re: [R-sig-phylo] Add terminal branches to tree >> >> Just to add to the slew of other good options, I have a small package up on >> GitHub (https://github.com/eliotmiller/addTaxa/tree/master/R) that is >> basically a big loop around Liam's bind.tip. There are examples included >> with the package. The main function you'd be interested in is >> randomlyAddTaxa(). It works fairly well last I checked, but there are some >> limitations and the package is in need of an overhaul. The main limitation >> is when two clades (A and B) are sister to one another and you intend to add >> a taxon into the clade (A+B); it's currently impossible to generate the >> topology (C,(A,B)), even though that's monophyletic and shouldn't be >> prohibited. I aim to fix this at some point within the next month or two if >> anyone's interested. >> Eliot >> >> On Mon, Jun 19, 2017 at 8:30 AM, Christoph Heibl <christoph.he...@gmx.net> >> wrote: >>> >>> Another possibility would be addSingleTip() and addTips() from the >>> megaptera package (https://github.com/heibl/megaptera). The latter function >>> can add any number of tips; the loop that Liam mentioned is already build >>> into that function. This is potentially very slow, but for your tree size it >>> should work. Note that you have to specify the 'tax' argument; it allows you >>> to add tips also to internal nodes of higher rank than genera. In your case >>> where you have only genera you can create 'tax' easily like this: >>> >>> ## 'species' is a vector of tip labels that you want to add >>> ## will work only if your tip labels are of the form >>> "Genus_epithet-or-any-other-string" >>> library(megaptera) >>> tax <- data.frame(genus = strip.spec(species), >>> species = species, >>> stringsAsFactors = FALSE) >>> >>> >>> >>> Christoph Heibl >>> An der Weiherleite 3 >>> 86633 Neuburg an der Donau >>> 08431-53 96 534 (Festnetz) >>> 0176-23 86 57 92 (Mobil) >>> christoph.he...@gmx.net >>> >>> >>> Am 19.06.2017 um 14:06 schrieb Liam Revell: >>> >>> > The function add.species.to.genus may do what you want. It adds a single >>> > species to the group defined by the MRCA of members of a genus, according >>> > to >>> > multiple criteria (randomly and so on). It can add only one species at a >>> > time, so you will need to write a for loop or something to iterate over >>> > the >>> > species that you‚d like to add. >>> > >>> > -- >>> > Liam J. Revell, Associate Professor of Biology >>> > University of Massachusetts Boston >>> > web: http://faculty.umb.edu/liam.revell >>> > email: liam.rev...@umb.ed
Re: [R-sig-phylo] Add terminal branches to tree
Thank you all for your suggestions. I guess the more direct way is to use Megaptera. I tried and it worked (I added Man cra1, Man cra2, etc.). However, I wanted it to create a new node that would basically include all these taxa. I sthere any way to merge all the added species in a new node? Thank you all, once again. > De: "Eliot Miller" <et...@cornell.edu> > Para: "Christoph Heibl" <christoph.he...@gmx.net> > Cc: "Liam Revell" <liam.rev...@umb.edu>, "r-sig-phylo" > <r-sig-phylo@r-project.org>, "Sergio Ferreira Cardoso" > <sergio.ferreira-card...@umontpellier.fr> > Enviadas: Segunda-feira, 19 De Junho de 2017 14:40:15 > Assunto: Re: [R-sig-phylo] Add terminal branches to tree > Just to add to the slew of other good options, I have a small package up on > GitHub ( https://github.com/eliotmiller/addTaxa/tree/master/R ) that is > basically a big loop around Liam's bind.tip. There are examples included with > the package. The main function you'd be interested in is randomlyAddTaxa(). It > works fairly well last I checked, but there are some limitations and the > package is in need of an overhaul. The main limitation is when two clades (A > and B) are sister to one another and you intend to add a taxon into the clade > (A+B); it's currently impossible to generate the topology (C,(A,B)), even > though that's monophyletic and shouldn't be prohibited. I aim to fix this at > some point within the next month or two if anyone's interested. > Eliot > On Mon, Jun 19, 2017 at 8:30 AM, Christoph Heibl < christoph.he...@gmx.net > > wrote: >> Another possibility would be addSingleTip() and addTips() from the megaptera >> package ( https://github.com/heibl/megaptera ). The latter function can add >> any >> number of tips; the loop that Liam mentioned is already build into that >> function. This is potentially very slow, but for your tree size it should >> work. >> Note that you have to specify the 'tax' argument; it allows you to add tips >> also to internal nodes of higher rank than genera. In your case where you >> have >> only genera you can create 'tax' easily like this: >> ## 'species' is a vector of tip labels that you want to add >> ## will work only if your tip labels are of the form >> "Genus_epithet-or-any-other-string" >> library(megaptera) >> tax <- data.frame(genus = strip.spec(species), >> species = species, >> stringsAsFactors = FALSE) >> Christoph Heibl >> An der Weiherleite 3 >> 86633 Neuburg an der Donau >> 08431-53 96 534 (Festnetz) >> 0176-23 86 57 92 (Mobil) >> christoph.he...@gmx.net >> Am 19.06.2017 um 14:06 schrieb Liam Revell: >>> The function add.species.to.genus may do what you want. It adds a single >>> species >>> to the group defined by the MRCA of members of a genus, according to >>> multiple >>> criteria (randomly and so on). It can add only one species at a time, so you >>> will need to write a for loop or something to iterate over the species that >> > you‚d like to add. >> > -- >> > Liam J. Revell, Associate Professor of Biology >> > University of Massachusetts Boston >> > web: http://faculty.umb.edu/liam.revell >> > email: liam.rev...@umb.edu >> > Sent from my Windows 10 phone >> > From: Sergio Ferreira Cardoso> > sergio.ferreira-card...@umontpellier.fr > >> > Sent: Monday, June 19, 2017 6:58 AM >> > To: r-sig-phylo >> > Subject: [R-sig-phylo] Add terminal branches to tree >> > Hello all, >>> I'm using the package 'phytools' to try to add terminal branches to a tree >>> (attached). I tried to use add.everywhere function to add terminal >>> branches. I >>> have 167 terminal taxa inside each of the 7 genera on my tree. Basically, I >>> just wanted to add some dozens of specimens to the end of each branch, but I >>> don't find a way to do it. Is there any tool that would allow me to do this? >> > Like adding 30 branches to 'Phal', adding 25 to 'Phat', etc. >> > Thanks in advance. >> > -- >> > Institut des Sciences de l'Evolution >> > UMR5554, CNRS, IRD, EPHE >> > Université de Montpellier >> > Place Eugène Bataillon >> > 34095 Montpellier Cedex 05 >> > France >> > Email: sergio.ferreira-card...@umontpellier.fr >> > Tel: +33 (4 ) 67 14 46 52 >> > [[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/ >> [[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/ ___ 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/
[R-sig-phylo] Add terminal branches to tree
Hello all, I'm using the package 'phytools' to try to add terminal branches to a tree (attached). I tried to use add.everywhere function to add terminal branches. I have 167 terminal taxa inside each of the 7 genera on my tree. Basically, I just wanted to add some dozens of specimens to the end of each branch, but I don't find a way to do it. Is there any tool that would allow me to do this? Like adding 30 branches to 'Phal', adding 25 to 'Phat', etc. Thanks in advance. -- Institut des Sciences de l'Evolution UMR5554, CNRS, IRD, EPHE Université de Montpellier Place Eugène Bataillon 34095 Montpellier Cedex 05 France Email: sergio.ferreira-card...@umontpellier.fr Tel: +33 (4 ) 67 14 46 52 #NEXUS BEGIN TREES; Title 'Pango_dist_1'; ID 015c8c71a66c9; LINK Taxa = Untitled_Block_of_Taxa; TRANSLATE [0] 1 Manc, [1] 2 Smug, [2] 3 Smut, [3] 4 Phat, [4] 5 Phal, [5] 6 Manj, [6] 7 Manp; TREE 'Tree_pangos_dist1' = ((1:2,(6:1,7:1):1):1,((2:1,3:1):1,(4:1,5:1):1):1):1[% ] [% ] [% setBetweenBits = selected ]; END; ___ 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/
Re: [R-sig-phylo] ape - corMartins
Dear Emmanuel, Thank you for your answer. I did try your suggestion and the alpha estimate was 0.89. I created models with alpha= 0.1, 1, 50 and a Brownian Motion model, which looks like it is the best fitting model. Model df AIC BIClogLik fitOU0.11 17 49.25612 81.06653 -7.628058 fitOU1 2 17 50.42329 82.23371 -8.211647 fitOU50 3 17 50.42330 82.23371 -8.211648 fitBrownian 4 17 42.82985 74.64027 -4.414926 When I do a stepwise regression, it's even more accentuated: > stepwise<-stepAIC(model,direction="both",trace=0) Model df AIC BIClogLik stepwise0.1 1 6 40.14718 51.37439 -14.07359 stepwise12 6 41.66723 52.89444 -14.83362 stepwise50 3 6 41.66724 52.89444 -14.83362 stepwiseBrownian 4 2 27.97681 31.71921 -11.98841 Also, the model fit does not increase with alpha >1. I tried the following: > fitPagel<-gls(fcl~mass+loc_dim+activity+feeding+loc_typ+agility,correlation=corPagel(value=1,phy=tree,fixed=FALSE),data=df,method="ML",weights=varFixed(~vf)) > fitPagel Generalized least squares fit by maximum likelihood Model: fcl ~ mass + loc_dim + activity + feeding + loc_typ + agility Data: df Log-likelihood: -3.274044 Coefficients: (Intercept) mass -0.3378136590.058271471 loc_dim3D activityDiurnal/Nocturnal 0.519780792 -0.202720741 activityNocturnal feedingOccasional_predator 0.057773877 -0.181687109 feedingPredator loc_typFlyer -0.153543188 -0.007564622 loc_typFossorial loc_typScansorial 0.7442470550.291622930 loc_typSemiaquatic loc_typTerrestrial 0.0292194950.449648236 agilityMedium agilityMedium_fast -0.072365659 -0.246153152 agilityMedium_slow agilityVery_slow -0.042370670 -0.255200809 Correlation Structure: corPagel Formula: ~1 Parameter estimate(s): lambda 0.9388239 Variance function: Structure: fixed weights Formula: ~vf Degrees of freedom: 48 total; 32 residual Residual standard error: 0.02930323 #Comparison between models Model df AIC BIClogLik fitOU0.11 17 49.25612 81.06653 -7.628058 fitOU1 2 17 50.42329 82.23371 -8.211647 fitOU50 3 17 50.42330 82.23371 -8.211648 fitBrownian 4 17 42.82985 74.64027 -4.414926 fitPagel5 18 42.54809 76.22971 -3.274044 Model df AIC BIClogLik stepwise0.1 1 6 40.14718 51.37439 -14.07359 stepwise12 6 41.66723 52.89444 -14.83362 stepwise50 3 6 41.66724 52.89444 -14.83362 stepwiseBrownian 4 2 27.97681 31.71921 -11.98841 stepwisePagel5 4 28.65373 36.13854 -10.32687 Since the Pagel estimates lambda values very close to 1, I believe this means my traits are evolving accoring to a Brownian Motion model. I'll try to search in the literature for a discussion on the estimation of alpha. Thank you, once again. Best regards, Sérgio. - Mensagem original - > De: "Emmanuel Paradis" <emmanuel.para...@ird.fr> > Para: "Sergio Ferreira Cardoso" <sff.card...@campus.fct.unl.pt>, > "r-sig-phylo" <r-sig-phylo@r-project.org> > Enviadas: Terça-feira, 24 De Janeiro de 2017 19:21:46 > Assunto: Re: [R-sig-phylo] ape - corMartins > Hi Sérgio, > > It seems your results make good sense. alpha=0 is too far from the > "optimum" value (i.e., the value that maximises the log-lik), so the > optimisation fails to maximise the log-likelihood with respect to this > parameter. It has been written in the literature that it is difficult to > estimate reliably alpha in an OU model, and your results seem to confirm > this. > > Besides, alpha=0 is equivalent to a Brownian motion model. So if your > traits are not really evolving according to this model, the GLS fitting > procedure may be in difficulty. > > You may try: > > correlation = corMartins(value=0.9, phy=tree, fixed=FALSE) > > If this estimates alpha=0.99, then you may probably rely on this result > (after checking the estimates of the coefficients of course). > > Best, > > Emmanuel > > Le 20/01/2017 à 15:16, Sergio Ferreira Cardoso a écrit : >> Dear all, >> >> I tried to estimate a value for
[R-sig-phylo] ape - corMartins
Dear all, I tried to estimate a value for alpha parameter of corMartins (Ornstein-Uhlenbeck) but the output is a bit puzzling. I do as follows: > fitOU<-gls(fcl~mass+loc_dim+activity+feeding+loc_typ+agility,correlation=corMartins(value=0,phy=tree,fixed=FALSE),data=df,method="ML",weights=varFixed(~vf)) > Error in recalc.corStruct(object[[i]], conLin) : NA/NaN/Inf in foreign function call (arg 4) Enter a frame number, or 0 to exit 1: gls(fcl ~ mass + loc_dim + activity + feeding + loc_typ + agility, correlati 2: nlminb(c(coef(glsSt)), function(glsPars) -logLik(glsSt, glsPars), control = 3: objective(.par, ...) 4: logLik(glsSt, glsPars) 5: logLik.glsStruct(glsSt, glsPars) 6: recalc(object, conLin) 7: recalc.modelStruct(object, conLin) 8: recalc(object[[i]], conLin) 9: recalc.corStruct(object[[i]], conLin) It work when I turn the "value" to 1, but then the estimate for alpha is 0.99, which looks like it isn't actually estimating a value and is instead accepting my input as the alpha value. Do you think there is a reason for this to happen? Thank you in advance. Cheers, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology Museu da Lourinhã, GEAL. LATR/IST/CTN - Campus Tecnológico e Nuclear. Lisboa, Portugal [[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/
Re: [R-sig-phylo] Amniote Phylogenetic Tree
Dear François, Thank you very much. Cheers, Sérgio. Em 13/01/2017 6:03 da tarde, "François Michonneau" < francois.michonn...@gmail.com> escreveu: Hi Sergio, You might find what you want within the Open Tree of Life (http://tree.opentreeoflife.org/) whose content can be accessed from R with the rotl package (https://cran.r-project.org/web/packages/rotl/index.html) Cheers, -- François On Fri, Jan 13, 2017 at 7:04 AM, Sergio Ferreira Cardoso <sff.card...@campus.fct.unl.pt> wrote: > Hello all, > > Do you know of any database where one can obtain a good phylogeny of the > Amniota (in NEXUS or NEWICK file)? > Thank you. > > Best regards, > Sérgio. > > -- > Com os melhores cumprimentos, > Sérgio Ferreira Cardoso. > > > > Best regards, > Sérgio Ferreira Cardoso > > > > > MSc. Paleontology > Museu da Lourinhã, GEAL. > LATR/IST/CTN - Campus Tecnológico e Nuclear. > > Lisboa, Portugal > > [[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-ph...@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/
[R-sig-phylo] Amniote Phylogenetic Tree
Hello all, Do you know of any database where one can obtain a good phylogeny of the Amniota (in NEXUS or NEWICK file)? Thank you. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology Museu da Lourinhã, GEAL. LATR/IST/CTN - Campus Tecnológico e Nuclear. Lisboa, Portugal [[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/
[R-sig-phylo] Phylogenetic signal
Dear all, I was trying to find the pylogenetic signal for a trait and I used phylosig(): phylosig(tree, trait, method="lambda", test=TRUE, nsim=999) If I use, for instance, a lambda correlation in a gls, am I expected to obtain the same value? gls(trait~var1+var2, corPagel(value=1,phy=tree,fixed=F)...) Is the lambda estimated for a model the same as the one calculated for a trait alone? If I have a different result, what's it due to? Thanks in advance. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology Museu da Lourinhã, GEAL. LATR/IST/CTN - Campus Tecnológico e Nuclear. Lisboa, Portugal ᐧ [[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/
Re: [R-sig-phylo] Error in trees[[i]] : subscript out of bounds
Deal Liam and Emmanuel, Thank you very much. I tried to open the file this way: obj<-scan(file="TimetreeOfLife2015.tre",n=1,skip=2,what="character") tree<-read.tree(text=obj) It worked but then there were other errors. Something is wrong with te edge lengths. I can't even plot the tree. I tried with read.tree but the errors persisted. When I used another tree, everything worked fine, so I guess I'll just use that one. Once again, thank you. Sérgio. 2016-05-23 20:55 GMT+01:00 Emmanuel Paradis <emmanuel.para...@ird.fr>: > Liam is probably right: there seems to be something missing in the NEXUS > file. > > The Newick file from the same site can also be read with read.tree. > > Emmanuel > > > Le 23/05/2016 18:21, Liam J. Revell a écrit : > >> It looks like the NEXUS file might be badly conformed. >> >> You could try: >> >> obj<-scan(file="TimetreeOfLife2015.tre",n=1,skip=2,what="character") >> tree<-read.tree(text=obj) >> >> That seemed to work for me. >> >> - Liam >> >> Liam J. Revell, Associate Professor of Biology >> University of Massachusetts Boston >> web: http://faculty.umb.edu/liam.revell/ >> email: liam.rev...@umb.edu >> blog: http://blog.phytools.org >> >> On 5/23/2016 12:11 PM, Sergio Ferreira Cardoso wrote: >> >>> Dear all, >>> >>> I am trying to open a .tre file with ape. I downloaded the tree from >>> http://www.timetree.org/book#treefiles. It's a .tre (not .nex) file, but >>> it's essentially a nexus. I tried to open it: >>> >>> tr <- read.nexus(file="TimetreeOfLife2015.tre") >>>> >>> Error in trees[[i]] : subscript out of bounds >>> >>> I tried to change the file name to .nex, but it wouldn't work. I saw some >>> previous threads but I couldn't solve the problem. >>> Does anyone know how can I solve this? >>> >>> Thanks in advance, >>> Sérgio. >>> >>> ps: I have the lastest 3.4 version of ape. >>> >>> >> ___ >> 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/ >> > -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology Museu da Lourinhã, GEAL. Lisboa, Portugal [[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/
[R-sig-phylo] Weird abline
Dear all, I'm analysing a small data set consisting of values for 27 taxa. I am extracting phylogenetic residuals from a regression (as in Revell, 2009) to be used as size corrected values. Almost all the residuals are negative because the regression line is positioned above what I would expect. Do you think this has something to do with the phylogeny? It probably has because everything looks fine with the ordinary least-squares. Thanks in advance. All the best, Sérgio. Here is the coding I used: ## Used packages library(ape, lib.loc=~/R/win-library/3.1) library(caper, lib.loc=~/R/win-library/3.1) library(nlme, lib.loc=~/R/win-library/3.1) library(phytools, lib.loc=~/R/win-library/3.1) library(phangorn, lib.loc=~/R/win-library/3.1) ## fitting PGLS fit_brain-pgls(perimeter~brain,comparative.data(tree,df,Species)) summary(fit_brain) ## cross checking using phytools phyl.resid(tree,brain,perimeter,method=BM) ## quick plotting plot(brain,perimeter) abline(fit_brain) ols-lm(perimeter~brain,data=df) abline(ols) -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] ANOVA type II with gls()
Dear members, Does anyone know a way to make anova SS II with gls class objects? I've tried Anova() from car package but it does not work with gls(). Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/
[R-sig-phylo] R-squared alternative for gls
Hello members, I'm creating relative values os semicircular canal size by fitting a PGLS and extracting residuals (using phyl.resid). I've heard that R-squared isn't meaningful in gls models. What I'm trying to do is to know is which of the two independent variables is best and on a lm () model I would just check the R-squared. Is the alternative, in the case of gls(), looking at the AIC or logLik? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/
Re: [R-sig-phylo] Different Lambda estimations
Hi Emmanuel, Thanks for your answer. So that means I can still use that Lambda estimated value, I suppose. In the case of Lambda1, it then means that species are more similar than expected, right? Best regards, Sérgio. 2015-07-17 15:57 GMT+01:00 Emmanuel Paradis emmanuel.para...@ird.fr: Hi Sérgio, In phytools, lambda is constrained to be positive, not in ape. A negative value of lambda is the consequence of closely related species being more dissimilar than expected under phylogenetic dependence (lambda = 1) or under independence (lambda = 0). This can be interpreted biologically as close species diverging into distinct niches and/or phenotypes. Best, Emmanuel Le 16/07/2015 20:34, Sergio Ferreira Cardoso a écrit : Hello all, I'm trying to estimate a lambda for my phylogeny but something strange happened. I used 3 different packages to confirm the estimated lambda value. My tree is ultrametric (with branchlengths in million years), it has 59 tips and 58 internal nodes. Here is what I did: df-data.frame(y,x) df$Species-rownames(df) tree-read.nexus(TTOL_tree.nex) vf-diag(vcv(tree)) phyl.resid(tree,x,y,method=lambda) ## phytools lambda = 0.741795 fit-gls(y~x, correlation = corPagel(value=1,phy=tree,fixed=FALSE), data = df, method=ML,weights=varFixed(~vf)) ## ape nlme lambda = -0.6863758 fit-pgls(y~x,comparative.data(tree,df,Species), lambda=ML) ## caper lambda = 0.742 Does anyone know why ape's gls estimates a value as different as that? Lambda varies between 0 and 1. Is there any other way to estimate Lambda (just to make another check). This is not the first time Lambda estimations of my trees deliver values as strange as this one... Cheers, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
Re: [R-sig-phylo] R-squared alternative for gls
Dear Dr. Garland, Then there is no problem because I just need to compare between phylogenetic models. Thank you. Best regards, Sérgio. 2015-07-20 14:48 GMT+01:00 Theodore Garland Jr theodore.garl...@ucr.edu: It's not that r squared isn't meaningful for generalized least squares but rather that it cannot be compared directly with values from OLS models. Cheers, Ted From: R-sig-phylo [r-sig-phylo-boun...@r-project.org] on behalf of Sergio Ferreira Cardoso [sff.card...@campus.fct.unl.pt] Sent: Monday, July 20, 2015 4:11 AM To: R phylo mailing list mailing list Subject: [R-sig-phylo] R-squared alternative for gls Hello members, I'm creating relative values os semicircular canal size by fitting a PGLS and extracting residuals (using phyl.resid). I've heard that R-squared isn't meaningful in gls models. What I'm trying to do is to know is which of the two independent variables is best and on a lm () model I would just check the R-squared. Is the alternative, in the case of gls(), looking at the AIC or logLik? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/ -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] Different Lambda estimations
Hello all, I'm trying to estimate a lambda for my phylogeny but something strange happened. I used 3 different packages to confirm the estimated lambda value. My tree is ultrametric (with branchlengths in million years), it has 59 tips and 58 internal nodes. Here is what I did: df-data.frame(y,x) df$Species-rownames(df) tree-read.nexus(TTOL_tree.nex) vf-diag(vcv(tree)) phyl.resid(tree,x,y,method=lambda) ## phytools lambda = 0.741795 fit-gls(y~x, correlation = corPagel(value=1,phy=tree,fixed=FALSE), data = df, method=ML,weights=varFixed(~vf)) ## ape nlme lambda = -0.6863758 fit-pgls(y~x,comparative.data(tree,df,Species), lambda=ML) ## caper lambda = 0.742 Does anyone know why ape's gls estimates a value as different as that? Lambda varies between 0 and 1. Is there any other way to estimate Lambda (just to make another check). This is not the first time Lambda estimations of my trees deliver values as strange as this one... Cheers, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
Re: [R-sig-phylo] Prune very large tree
Thank you all for the suggestions. Creating a vector with the species works very well. :) 2015-07-15 17:59 GMT+01:00 Sergio Ferreira Cardoso sff.card...@campus.fct.unl.pt: Thank you very much Nicholas. Works perfectly!!! :) 2015-07-15 17:22 GMT+01:00 Nicholas Crouch ncro...@uic.edu: Sergio, You can use ape's 'drop.tip', e.g. # Phylogeny is object 'phy' # Species interested in in object 'my.species' v - phy$tip.label %in% my.species drop - phy$tip.label[v==FALSE] trimmed.phy - drop.tip(phy, drop) On Wed, Jul 15, 2015 at 11:15 AM, Sergio Ferreira Cardoso sff.card...@campus.fct.unl.pt wrote: Dear members, I have a tree with almost 10.000 species and I want to prune it. The problem is I only want a tree with 59 of those species. Is there any possible way to make this in R? Is there any code to prune all taxa except the ones on my data? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/ -- Nicholas Crouch Graduate Student, Igić Lab Department of Biological Sciences University of Illinois at Chicago -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] Build consensus trees
Dear all, I have a nexus file with 1000 trees on it with divergence times as branch lengths. Does anyone know about a package that can make a consensus tree out of my 1000 trees? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
Re: [R-sig-phylo] Build consensus trees
Ok. I just solved the problem with adding file= when using wirte.nexus. 2015-07-15 14:56 GMT+01:00 Sergio Ferreira Cardoso sff.card...@campus.fct.unl.pt: Thank you François, I was trying to use consensus() from ape, but when I try to write a nexus file I get an error message: trees-read.nexus(1000_trees.nex) consensus(trees,p=1,check.labels=TRUE) write.nexus(my_tree,consensus_birds.nex) #NEXUS [R-package APE, Wed Jul 15 14:49:14 2015] BEGIN TAXA; DIMENSIONS NTAX = 59; TAXLABELS Rhynchotus_rufescens Apteryx_haastii Dromaius_novaehollandiae (...) ; END; BEGIN TREES; TRANSLATE Error in y$tip.label : $ operator is invalid for atomic vectors My only problem now is to create the file. I'm going to try the softwares you pointed. Best regards, Sérgio. ᐧ 2015-07-15 14:49 GMT+01:00 François Michonneau francois.michonn...@gmail.com: Hi Sergio, Not in R, but you can use PAUP, MrBayes or Dendropy (my recommendation). Cheers, -- François 2015-07-15 7:53 GMT-04:00 Sergio Ferreira Cardoso sff.card...@campus.fct.unl.pt: Dear all, I have a nexus file with 1000 trees on it with divergence times as branch lengths. Does anyone know about a package that can make a consensus tree out of my 1000 trees? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/ -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
Re: [R-sig-phylo] Build consensus trees
Thank you François, I was trying to use consensus() from ape, but when I try to write a nexus file I get an error message: trees-read.nexus(1000_trees.nex) consensus(trees,p=1,check.labels=TRUE) write.nexus(my_tree,consensus_birds.nex) #NEXUS [R-package APE, Wed Jul 15 14:49:14 2015] BEGIN TAXA; DIMENSIONS NTAX = 59; TAXLABELS Rhynchotus_rufescens Apteryx_haastii Dromaius_novaehollandiae (...) ; END; BEGIN TREES; TRANSLATE Error in y$tip.label : $ operator is invalid for atomic vectors My only problem now is to create the file. I'm going to try the softwares you pointed. Best regards, Sérgio. ᐧ 2015-07-15 14:49 GMT+01:00 François Michonneau francois.michonn...@gmail.com: Hi Sergio, Not in R, but you can use PAUP, MrBayes or Dendropy (my recommendation). Cheers, -- François 2015-07-15 7:53 GMT-04:00 Sergio Ferreira Cardoso sff.card...@campus.fct.unl.pt: Dear all, I have a nexus file with 1000 trees on it with divergence times as branch lengths. Does anyone know about a package that can make a consensus tree out of my 1000 trees? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/ -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] Prune very large tree
Dear members, I have a tree with almost 10.000 species and I want to prune it. The problem is I only want a tree with 59 of those species. Is there any possible way to make this in R? Is there any code to prune all taxa except the ones on my data? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/
[R-sig-phylo] phyres function R package caper
Dear all, When I try to get a list os phylogenetic residuals using phyres function from R package I get this message: Error: could not find function phyres. Does anyone know how to solve this problem? phyres(fit.gls1) Error: could not find function phyres Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
Re: [R-sig-phylo] phyres function R package caper
Thank you both. Best regards, Sérgio. No dia 24 de Jun de 2015 18:30, Theodore Garland Jr theodore.garl...@ucr.edu escreveu: Thank you, Liam! Cheers, Ted From: Liam J. Revell [liam.rev...@umb.edu] Sent: Wednesday, June 24, 2015 10:24 AM To: Theodore Garland Jr; Sergio Ferreira Cardoso; R phylo mailing list mailing list Subject: Re: [R-sig-phylo] phyres function R package caper Hi all. To the original question, you should be able to get these values first using gls(...,correlation=corBrownian(...)) in nlme then applying residuals to the fitted model returned by gls. For instance, for data frame X with variables x y, and ultrametric phylogeny tree, you might compute: library(ape) library(nlme) fit-gls(y~x,data=X,correlation=corBrownian(1,tree)) residuals(fit) (phytools also has a function for this, phyl.resid, but it does exactly the same thing as the code above, and thus there is really no reason to prefer that function - except perhaps to cross-check your result for errors.) With regards to Ted's comment, indeed these are different quantities. Though the fitted coefficients from a contrasts regression should be the same as above, the residuals will be different (and there will be one fewer of them, besides). These residuals, from the contrasts regression, should be phylogenetically independent; however they are not longer associated with species, but with 'contrasts' or nodes in the tree. To obtain these residuals from a contrasts regression you should be able to do something like: X-X[tree$tip.label,] ## precautionary pic.x-pic(X[,x],tree) pic.y-pic(X[,y],tree) fit-lm(pic.y~pic.x-1) residuals(fit) There is no particular reason to prefer one set of quantities over the other - it just depends on what subsequent analyses are intended. In the former case, the residuals are associated with species - but these residuals consequently will be phylogenetically correlated thus the tree needs to be taken into consideration in any subsequent analysis. The latter residuals are phylogenetically independent, but no longer associated with species. (I hate to cite myself, but this is discussed in my paper Revell 2009; Evolution.) I hope this is of some help. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 6/24/2015 12:35 PM, Theodore Garland Jr wrote: Hi All, I am going to suggest that when people want any sort of phylogenetic residuals they do some checking on their own to try to verify what, exactly, they are getting. Here's one check you can do. Compute phylogenetically independent contrasts for two traits. Perform a regression (through the origin, of course) of one trait on the other. Save the residuals. Compare them with the phylogenetic residuals you get from some other program that does a PGLS regression (not with any transformation of the branch lengths). Let us know what you find! Cheers, Ted Theodore Garland, Jr., Professor Department of Biology University of California, Riverside Riverside, CA 92521 Office Phone: (951) 827-3524 Facsimile: (951) 827-4286 (not confidential) Email: tgarl...@ucr.edu http://www.biology.ucr.edu/people/faculty/Garland.html http://scholar.google.com/citations?hl=enuser=iSSbrhwJ Director, UCR Institute for the Development of Educational Applications Editor in Chief, Physiological and Biochemical Zoology Fail Lab: Episode One http://testtube.com/faillab/zoochosis-episode-one-evolution http://www.youtube.com/watch?v=c0msBWyTzU0 From: R-sig-phylo [r-sig-phylo-boun...@r-project.org] on behalf of Sergio Ferreira Cardoso [sff.card...@campus.fct.unl.pt] Sent: Wednesday, June 24, 2015 9:21 AM To: R phylo mailing list mailing list Subject: [R-sig-phylo] phyres function R package caper Dear all, When I try to get a list os phylogenetic residuals using phyres function from R package I get this message: Error: could not find function phyres. Does anyone know how to solve this problem? phyres(fit.gls1) Error: could not find function phyres Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/ ___ R-sig-phylo mailing list - R-sig-phylo
Re: [R-sig-phylo] Non Parametric PGLS
Hi Liam, Again, thank you for the answer. Yes, I'm aware that they are phylogenetically correlated and that they need to be subsequently analysed with methods such as PGLS. In fact, when I apply a normality test to ( chol(solve(vcv(tree)))%*%residuals(fit))the transformed residuals are normally distributed. But I understand that the residuals I need to use as size-corrected traits are the ones without the transformation (that, in my case, aren't normal, and that's why I was a bit worried). Thanks a lot for the help. Cheers, Sérgio. 2015-06-18 15:54 GMT+01:00 Liam J. Revell liam.rev...@umb.edu: Hi Sérgio. What Simon ( I, in my blog post) suggested is that to test the 'normality assumption' you need to first transform the residuals with the inverse Cholesky decomposite matrix. This will give you a vector in which the values should be normal independent (assuming that the correlation structure of the residuals is properly given by the tree). However these values are no longer associated with species(!) so they cannot be used in subsequent among-species analyses. So, basically, yes - you should use the residuals before transformation in cross-species analyses (for instance, as 'size corrected' trait values); however you need to remember that they are still phylogenetically correlated. I hope that's clear enough. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 6/18/2015 6:16 AM, Sergio Ferreira Cardoso wrote: Hello again, Thanks for your help. This kind of solved my problem. I normally use some kind of test (shapiro or komolgorov) to test for normality. I know histograms or qqnorm plots a more helpful, but they are more vulnerable to each others interpretation. So, just to make clear one thing: these kind of analysis of the residual error ( http://blog.phytools.org/2013/02/a-comment-on-distribution-of-residuals.html ) has nothing to do with phylogenetic residuals taken from a regression (in order to obtain relative values - size-correction as in Revell, L. J. (2009). Size‐correction and principal components for interspecific comparative studies. Evolution, 63(12), 3258-3268.). So even if residuals from a PGLS aren't normal, this means I can still use them as size-corrected values for a certain trait, correct? Once again, thank you very much for your advices. Best regards, Sérgio. 2015-06-18 4:56 GMT+01:00 Simon Blomberg s.blombe...@uq.edu.au mailto:s.blombe...@uq.edu.au: Hi Sérgio. Liam is right. But we do expect the normalised residuals to be approximately Normal. You can calculate the normalised residuals by pre-multiplying the raw residuals by the inverse of the Cholesky decomposition of the phylogenetic variance-covariance matrix, and then dividing by the estimate of the residual standard deviation (ie sigma). You may have to plug in a value for any parameters that are co-estimated along with the regression (Pagel's lambda, etc.). If you use the gls function in the nlme package to fit your model, then it's all easy: fit - gls(response~explanatory, correlation=corPagel(1, phy=my.tree), data=dat)) res - residuals(fit, type=normalized) Then you can do some test for normality on those (I don't ordinarily recommend such things, although SnowsPenultimateNormalityTest in the TeachingDemos package is the best I have seen). More usefully, you can do a normal quantile-quantile plot to graphically see whether your normalised residuals are normal enough: qqnorm(fit, form=~residuals(., type=n), abline=c(0,1)) See page 239 in Pinheiro and Bates (2000) Mixed-effects models in S and S-PLUS. Cheers, Simon. On 18/06/15 03:09, Liam J. Revell wrote: Hi Sérgio. It might be worth pointing out that we do not expect that the residuals from a phylogenetic regression to be normal. I described this with respect to the phylogenetic ANOVA on my blog ( http://blog.phytools.org/2013/02/a-comment-on-distribution-of-residuals.html ), but this applies equally to phylogenetic regression. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu mailto:liam.rev...@umb.edu blog: http://blog.phytools.org On 6/17/2015 12:40 PM, Sergio Ferreira Cardoso wrote: Hello all, I'm having a problem with a Multiple Regression PGLS analysis that I'm performing. The residuals are not normal and it's difficult to bring them to normality. In these cases, are there any alternatives to the linear
Re: [R-sig-phylo] Non Parametric PGLS
Hello again, Thanks for your help. This kind of solved my problem. I normally use some kind of test (shapiro or komolgorov) to test for normality. I know histograms or qqnorm plots a more helpful, but they are more vulnerable to each others interpretation. So, just to make clear one thing: these kind of analysis of the residual error ( http://blog.phytools.org/2013/02/a-comment-on-distribution-of-residuals.html) has nothing to do with phylogenetic residuals taken from a regression (in order to obtain relative values - size-correction as in Revell, L. J. (2009). Size‐correction and principal components for interspecific comparative studies. Evolution, 63(12), 3258-3268.). So even if residuals from a PGLS aren't normal, this means I can still use them as size-corrected values for a certain trait, correct? Once again, thank you very much for your advices. Best regards, Sérgio. 2015-06-18 4:56 GMT+01:00 Simon Blomberg s.blombe...@uq.edu.au: Hi Sérgio. Liam is right. But we do expect the normalised residuals to be approximately Normal. You can calculate the normalised residuals by pre-multiplying the raw residuals by the inverse of the Cholesky decomposition of the phylogenetic variance-covariance matrix, and then dividing by the estimate of the residual standard deviation (ie sigma). You may have to plug in a value for any parameters that are co-estimated along with the regression (Pagel's lambda, etc.). If you use the gls function in the nlme package to fit your model, then it's all easy: fit - gls(response~explanatory, correlation=corPagel(1, phy=my.tree), data=dat)) res - residuals(fit, type=normalized) Then you can do some test for normality on those (I don't ordinarily recommend such things, although SnowsPenultimateNormalityTest in the TeachingDemos package is the best I have seen). More usefully, you can do a normal quantile-quantile plot to graphically see whether your normalised residuals are normal enough: qqnorm(fit, form=~residuals(., type=n), abline=c(0,1)) See page 239 in Pinheiro and Bates (2000) Mixed-effects models in S and S-PLUS. Cheers, Simon. On 18/06/15 03:09, Liam J. Revell wrote: Hi Sérgio. It might be worth pointing out that we do not expect that the residuals from a phylogenetic regression to be normal. I described this with respect to the phylogenetic ANOVA on my blog ( http://blog.phytools.org/2013/02/a-comment-on-distribution-of-residuals.html), but this applies equally to phylogenetic regression. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 6/17/2015 12:40 PM, Sergio Ferreira Cardoso wrote: Hello all, I'm having a problem with a Multiple Regression PGLS analysis that I'm performing. The residuals are not normal and it's difficult to bring them to normality. In these cases, are there any alternatives to the linear model? I know that for non-phylogenetic analyses other models exist, but is there any alternative method for phylogenetic analysis? Thanks in advance. Best regards, Sérgio. ___ 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/ -- Simon Blomberg, BSc (Hons), PhD, MAppStat, AStat. Senior Lecturer and Consultant Statistician School of Biological Sciences The University of Queensland St. Lucia Queensland 4072 Australia T: +61 7 3365 2506 email: S.Blomberg1_at_uq.edu.au http://www.evolutionarystatistics.org Policies: 1. I will NOT analyse your data for you. 2. Your deadline is your problem. Basically, I'm not interested in doing research and I never have been. I'm interested in understanding, which is quite a different thing. - David Blackwell ___ 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/ -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] Non Parametric PGLS
Hello all, I'm having a problem with a Multiple Regression PGLS analysis that I'm performing. The residuals are not normal and it's difficult to bring them to normality. In these cases, are there any alternatives to the linear model? I know that for non-phylogenetic analyses other models exist, but is there any alternative method for phylogenetic analysis? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] Phylogenetic size-correction
Dear all, I just read this paper: Revell, L. J. (2009). Size‐correction and principal components for interspecific comparative studies. Evolution, 63(12), 3258-3268. There is no coding for R, but it says it can be made using APE. From what I could understand this is as simples as: -bm.gls-gls(trait~Bodymass, correlation=BM,data=DF,weights=varFixed(~vf)) -residuals(bm.gls) Is this correct? Or do I need to ask for some special kind of residuals? Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/
Re: [R-sig-phylo] Phylogenetic size-correction
Hello, I'm sorry for not specifying that. But yes, that's exactly it. Thanks for the help. Cheers, Sérgio. 2015-06-15 12:47 GMT+01:00 Liam J. Revell liam.rev...@umb.edu: Hi Sergio. You didn't specify what BM vf are, but if: BM-corBrownian(1,phy) ## and vf-diag(vcv(phy)) then, yes, this is correct. There is a function in phytools, phy.resid, that also does this calculation. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 6/15/2015 7:41 AM, Sergio Ferreira Cardoso wrote: Dear all, I just read this paper: Revell, L. J. (2009). Size‐correction and principal components for interspecific comparative studies. Evolution, 63(12), 3258-3268. There is no coding for R, but it says it can be made using APE. From what I could understand this is as simples as: -bm.gls-gls(trait~Bodymass, correlation=BM,data=DF,weights=varFixed(~vf)) -residuals(bm.gls) Is this correct? Or do I need to ask for some special kind of residuals? Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
Re: [R-sig-phylo] Phylogenetic PCA
Thank you all for your suggestions. I'll check them all. Thank you once again. Cheers, Sérgio. 2015-06-04 11:35 GMT+01:00 Thomas Püschel thomaspusc...@gmail.com: Hi Sergio, there you go http://blog.phytools.org/2011/04/new-function-phylpca.html http://cran.r-project.org/web/packages/phytools/index.html best wishes 2015-06-04 11:26 GMT+01:00 Sergio Ferreira Cardoso sff.card...@campus.fct.unl.pt: Dear all, I'm wondering if there is a package in R with which I can run a phylogenetic PCA. Please, if you know something or if you you know another software useful to use this technique, let me know. Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/ -- Thomas Püschel PhD student in Adaptive Organismal Biology Computational and Evolutionary Biology Group Faculty of Life Sciences University of Manchester Michael Smith Building Oxford Road Manchester M13 9PT United Kingdom Email: thomas.pusc...@postgrad.manchester.ac.uk thomas.puschelroul...@postgrad.manchester.ac.uk; thomaspusc...@gmail.com -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] Phylogenetic PCA
Dear all, I'm wondering if there is a package in R with which I can run a phylogenetic PCA. Please, if you know something or if you you know another software useful to use this technique, let me know. Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] Residuals vs prediction intervals
Hello all, This isn't a specific question about phylogenies. Whenever I plot prediction intervals for a linear model or a PGLS I can see what observations fall out of the prediction intervals. My question is: if the residuals of whatever regression model I'm using are normally/aproximatelly normally distributed, then should I drop the observations falling out of the prediction intervals? My question is, should take observations out of the analysis because they fall out of the prediction intervals, because they are severe outliers on a a normality plot or both. Thank you very much, in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/
Re: [R-sig-phylo] Residuals vs prediction intervals
Please, ignore this e-mail. Sent by mistake. Senseless. Regards, Sérgio. No dia 16 de Mai de 2015 16:10, Sergio Ferreira Cardoso sff.card...@campus.fct.unl.pt escreveu: Hello all, This isn't a specific question about phylogenies. Whenever I plot prediction intervals for a linear model or a PGLS I can see what observations fall out of the prediction intervals. My question is: if the residuals of whatever regression model I'm using are normally/aproximatelly normally distributed, then should I drop the observations falling out of the prediction intervals? My question is, should take observations out of the analysis because they fall out of the prediction intervals, because they are severe outliers on a a normality plot or both. Thank you very much, in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/
Re: [R-sig-phylo] Non-ultrametric tree PGLS
Thank you very much for your answers. They were really helpful. Regards, Sérgio. 2015-05-14 17:00 GMT+01:00 Theodore Garland Jr theodore.garl...@ucr.edu: Agreed! Similarly, Pagel's lambda ad Grafen's rho were designed from a purely statistical perspective, whereas OU and ACDC models are motivated by ties to a possible set of biological processes. Cheers, Ted *From:* Alejandro Gonzalez Voyer [ alejandro.gonza...@iecologia.unam.mx] *Sent:* Thursday, May 14, 2015 8:47 AM *To:* Theodore Garland Jr *Cc:* Sergio Ferreira Cardoso; R phylo mailing list mailing list *Subject:* Re: [R-sig-phylo] Non-ultrametric tree PGLS Hi Sergio, I would add to Ted’s reply that you are not only considering alternative statistical models but that the evolutionary assumptions from the models you are fitting also differ, and you need to keep this in mind when comparing models. Comparison of AICs or any other estimate of goodness of fit must also involve careful consideration of the assumptions of the models you are comparing. In your particular case, the tree with branch lengths set to equal values (all branch lengths = 1) implies different amount of time to evolve for each of your species (in other words the expected variances - diagonal terms in the variance-covariance matrix - differ between the species), and thus you should consider whether such an assumption makes biological sense in your system. Cheers Alejandro ___ Dr Alejandro Gonzalez Voyer Laboratorio de Conducta Animal Instituto de Ecología Circuito Exterior S/N Ciudad Universitaria Universidad Nacional Autónoma de México México, D.F. 04510 México Tel: +52 55 5622 9044 E-mail: alejandro.gonza...@iecologia.unam.mx El 14/05/2015, a las 10:39, Theodore Garland Jr theodore.garl...@ucr.edu escribió: Hi Sergio, I am not quite understanding the situation nor why you see a problem. if I understand correctly, you are considering these five (5) alternative models for some sort of simple or multiple regression: OLS = star phylogeny PGLS with real-time branch lengths (ultrametric) Pagel's lambda with real-time branch lengths (ultrametric) PGLS with all branch legnths equal to 1.0 Pagel's lambda with all branch lengths equal to 1.0 To help decide which model best fits your data, you can look at AIC or for some comparisons do a likelihood ratio test. My experience is that any of the transform models (Pagel's lambda, Grafen's rho, OU in various implementations, ACDC) can sometimes yield really bizarre results when you start with a non-ultrametric tree. You need to be careful and check the REML likelihood surface for multiple peaks, etc. Cheers, Ted Theodore Garland, Jr., Professor Department of Biology University of California, Riverside Riverside, CA 92521 Office Phone: (951) 827-3524 Facsimile: (951) 827-4286 (not confidential) Email: tgarl...@ucr.edu http://www.biology.ucr.edu/people/faculty/Garland.html http://scholar.google.com/citations?hl=enuser=iSSbrhwJ Director, UCR Institute for the Development of Educational Applications Editor in Chief, Physiological and Biochemical Zoology Fail Lab: Episode One http://testtube.com/faillab/zoochosis-episode-one-evolution http://www.youtube.com/watch?v=c0msBWyTzU0 From: R-sig-phylo [r-sig-phylo-boun...@r-project.org] on behalf of Sergio Ferreira Cardoso [sff.card...@campus.fct.unl.pt] Sent: Thursday, May 14, 2015 8:32 AM To: r-sig-phylo@r-project.org Subject: [R-sig-phylo] Non-ultrametric tree PGLS Hello all, I have an ultrametric phylogenetic tree with divergence times as branch lengths. To see if there was a big difference between using these branch lengths and equal (=1) branch lengths I set all lengths to 1 and ran a PGLS. I ran with Lambda transformations and the estimation is that Lambda is superior than 1 (both with ML and REML estimation). I suppose this is a consequence of the tree being non ultrametric. Is there a solution for this problem or should I, in this case, just ran a GLS (Brownian Motion) to avoid the over estimation of the phylogenetic signal? Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/ ___ R-sig-phylo mailing list - R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo Searchable
[R-sig-phylo] Non-ultrametric tree PGLS
Hello all, I have an ultrametric phylogenetic tree with divergence times as branch lengths. To see if there was a big difference between using these branch lengths and equal (=1) branch lengths I set all lengths to 1 and ran a PGLS. I ran with Lambda transformations and the estimation is that Lambda is superior than 1 (both with ML and REML estimation). I suppose this is a consequence of the tree being non ultrametric. Is there a solution for this problem or should I, in this case, just ran a GLS (Brownian Motion) to avoid the over estimation of the phylogenetic signal? Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/
[R-sig-phylo] Weird estimated Lambda values (PGLS)
Dear all, I performed a PGLS with a tree with divergence times as branch lengths. I used 1 dependent and 1 independent variable. When I estimate the Lambda (in Pagel's method) I get the Maximum Likelihood Lambda ~0 and the Restricted Maximum Likelihood Lambda 0.878054. This is really strange. I was thinking I might have a tree with incorrect branchlengths but I checked the plot of absolute value of contrast vs. standard deviation for both variables and it's ok (according to Garland et al., 1992). With OU transformation I don't have the same problem. However, the p-value of phylogenetic (both with OU or with Brownian Motion PGLS) aren't very different from the ordinary least-squares. Is there any other ordinary procedure/analysis I should make other than the plots I checked to check my tree? (I understand that there is a chance that the characters I am using have influence and maybe the tree doesn't make any difference). Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal ᐧ [[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/
Re: [R-sig-phylo] Non normal PGLS results
Hello, Thank you very much for the answer. I assumed PGLS residuals had to be normal (I saw this thread: https://stat.ethz.ch/pipermail/r-sig-phylo/2012-May/002064.html). In fact, what I'm analysing are standardized residuals (standardized by setting the determinant of the covariance matrix equal to one). I thought something was wrong with the residuals because when estimating the Lamba using ML I get a completely different result from estimating it with REML (e.g., REML est: 0.945; M est: ~0). Do you have any idea of wat could be causing this? Once again, thank you for answering. Best regards, Sérgio. ᐧ 2015-04-18 18:22 GMT+01:00 Liam J. Revell liam.rev...@umb.edu: Hi Sergio. We don't expect the residuals from PGLS to be independent draws from a normal distribution, but multivariate normal with a correlation structure given by the tree. Here I give some more explanation of this on my blog: http://blog.phytools.org/2013/02/a-comment-on-distribution-of-residuals.html . Let us know if this is helpful. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 4/18/2015 1:15 PM, Sergio Ferreira Cardoso wrote: Dear all, I'm performing PGLS's with an ultrametric phylogenetic tree (divergence time as branchlengths). I tied Pagel's Lambda transformation, OU transformation, regular GLS and OLS, to compare results. There is one problem: the residuals of my analyses are not normal. I tried to remove big outliers but it made things even worse, because without them there are even more outliers. So, what should I do? The results are consistent in both 4 tests. But can I trust the results of the PGLS? Is it particularly bad if PGLS residuals aren't normal? Does it critically afect my results? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] Non normal PGLS results
Dear all, I'm performing PGLS's with an ultrametric phylogenetic tree (divergence time as branchlengths). I tied Pagel's Lambda transformation, OU transformation, regular GLS and OLS, to compare results. There is one problem: the residuals of my analyses are not normal. I tried to remove big outliers but it made things even worse, because without them there are even more outliers. So, what should I do? The results are consistent in both 4 tests. But can I trust the results of the PGLS? Is it particularly bad if PGLS residuals aren't normal? Does it critically afect my results? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
Re: [R-sig-phylo] Non normal PGLS results
Hi Dr. Garland, I just ran a normal OLS and the residuals are normal.Actually, I thing my data has a lot of what you call long singleton branches. I'm sending a pdf. of it attached. Best regards, Sérgio. ᐧ 2015-04-18 19:09 GMT+01:00 Theodore Garland Jr theodore.garl...@ucr.edu: If you get such a large difference between ML and REML estimation in this sort of situation then probably either (1) something is wrong with the code (bad search algorithm?) or (2) you have something pathological in your tip data set and/or the tree (e.g., some really long singleton branches or really short sister branches). If you see a big outlier in the residuals, then this is a problem - I would not trust the results. How do the residuals look from a regular OLS analysis? Cheers, Ted From: R-sig-phylo [r-sig-phylo-boun...@r-project.org] on behalf of Sergio Ferreira Cardoso [sff.card...@campus.fct.unl.pt] Sent: Saturday, April 18, 2015 10:59 AM To: Liam J. Revell Cc: r-sig-phylo@r-project.org Subject: Re: [R-sig-phylo] Non normal PGLS results Hello, Thank you very much for the answer. I assumed PGLS residuals had to be normal (I saw this thread: https://stat.ethz.ch/pipermail/r-sig-phylo/2012-May/002064.html). In fact, what I'm analysing are standardized residuals (standardized by setting the determinant of the covariance matrix equal to one). I thought something was wrong with the residuals because when estimating the Lamba using ML I get a completely different result from estimating it with REML (e.g., REML est: 0.945; M est: ~0). Do you have any idea of wat could be causing this? Once again, thank you for answering. Best regards, Sérgio. ᐧ 2015-04-18 18:22 GMT+01:00 Liam J. Revell liam.rev...@umb.edu: Hi Sergio. We don't expect the residuals from PGLS to be independent draws from a normal distribution, but multivariate normal with a correlation structure given by the tree. Here I give some more explanation of this on my blog: http://blog.phytools.org/2013/02/a-comment-on-distribution-of-residuals.html . Let us know if this is helpful. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 4/18/2015 1:15 PM, Sergio Ferreira Cardoso wrote: Dear all, I'm performing PGLS's with an ultrametric phylogenetic tree (divergence time as branchlengths). I tied Pagel's Lambda transformation, OU transformation, regular GLS and OLS, to compare results. There is one problem: the residuals of my analyses are not normal. I tried to remove big outliers but it made things even worse, because without them there are even more outliers. So, what should I do? The results are consistent in both 4 tests. But can I trust the results of the PGLS? Is it particularly bad if PGLS residuals aren't normal? Does it critically afect my results? Thanks in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/ -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal tree_birds Description: Binary data ___ 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/
Re: [R-sig-phylo] PGLS transformations
Thank you Liam. I'll check it now. Best regards, Sérgio. 2015-04-14 18:27 GMT+01:00 Liam J. Revell liam.rev...@umb.edu: Hi Sergio. Yesterday you asked if it would be reasonable to use the tree with branch lengths simulated assuming a Yule process if the topology is known in phylogenetic regression. I wondered what would happen if, for a given tree topology, we simulated a set of trees with branches as if the tree arose under a Yule process, and then fit our regression model to each tree. We could then add the variance among fitted regression slopes to the mean variance of each slope to get the total variance, from which we could test a hypothesis that (for instance) the slope is not different from zero. Well, I tried this (details here: http://blog.phytools.org/2015/ 04/phylogenetic-regression-when-branch.html), and, assuming I have done so correctly, it seems as though the approach is too conservative - leading to low power and a type I error rate lower than the nominal level. Using arbitrary branch lengths (all branch lengths equal to 1.0, Grafen's branch lengths) results in elevated type I error, but not incredibly badly inflated type I error. I thought you other R-sig-phylo readers might be interested in this result (newly obtained), or could correct it if I have done something wrong. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 4/13/2015 2:20 PM, Liam J. Revell wrote: I wouldn't go as far as to call it a suggestion - but this is one thing that you could do. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 4/13/2015 2:04 PM, Sergio Ferreira Cardoso wrote: Hello, So what you're suggesting is that I do that and then fit-gls(X~Y,corr=corBrownian(1,t.pb),data=df,method=ML, right? Thanks. Best regards, Sérgio. ᐧ 2015-04-13 18:42 GMT+01:00 Liam J. Revell liam.rev...@umb.edu mailto:liam.rev...@umb.edu: Hi all. I just wanted to add that it should be straightforward to sample edge lengths as if they arose under a specific process. For instance, here on my blog I show how to sample branching times edge lengths as if they arose under a pure-birth (i.e., 'Yule') process, given a particular input topology: http://blog.phytools.org/2015/__04/sampling-edge-lengths-__ under-yule-process.html http://blog.phytools.org/2015/04/sampling-edge-lengths- under-yule-process.html. For things like the phylogenetic (contrasts) regression, this may be preferable to Grafen's edge lengths, which tend to make recent edge lengths very short, giving very high weight to the associated contrasts. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.__revell/ http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu mailto:liam.rev...@umb.edu blog: http://blog.phytools.org On 4/13/2015 1:13 PM, Sergio Ferreira Cardoso wrote: Hello, Thank you both for the help. Emmanuel, so is there a way to see contrasts in R? Reading this paper - Garland, T., Harvey, P. H., Ives, A. R. (1992). Procedures for the analysis of comparative data using phylogenetically independent contrasts. Systematic Biology, 41(1), 18-32. - I was aware of the importance of standardizing contrasts and of comparing Absolute value of standard contrat vs Standard deviation of contrast. I've learned to do this in Mesquite but I don't know if R alows this to be done. When I run the GLS I ask R to estimate the rho, so, a priori I never know the rho value. Maybe I'm really really confused, but here is the reason why I tthink something isn't right with my analysis: I built a phylogenetic tree in Mesquite, and based on several works I used million years as branch lengths. It makes sense for me because I'll be using a fossil on my analysis. But I tested the tree before adding the extinct taxa, so to make sure everything was OK when the tree was ultrametric. I noticed that, for example, whatever the independent variable (X) was, the alpha from OU was extremely high (0.999182, for instance). It would always be 0.999. I thought maybe I was using the wrong transformation... That's why I ended up trying to know if I needed to do something prior to transforming and analysing the tree. Thank you very much. Best regards
Re: [R-sig-phylo] PGLS transformations
Than you all. Best regards, Sérgio. ᐧ 2015-04-13 22:02 GMT+01:00 Emmanuel Paradis emmanuel.para...@ird.fr: Le 13/04/2015 19:13, Sergio Ferreira Cardoso a écrit : Hello, Thank you both for the help. Emmanuel, so is there a way to see contrasts in R? Yes! The function pic() has been in ape since its first release (ver. 0.1 in Aug 2002), and there is an option to output the contrasts scaled or not. best, Emmanuel Reading this paper - Garland, T., Harvey, P. H., Ives, A. R. (1992). Procedures for the analysis of comparative data using phylogenetically independent contrasts. Systematic Biology, 41(1), 18-32. - I was aware of the importance of standardizing contrasts and of comparing Absolute value of standard contrat vs Standard deviation of contrast. I've learned to do this in Mesquite but I don't know if R alows this to be done. When I run the GLS I ask R to estimate the rho, so, a priori I never know the rho value. Maybe I'm really really confused, but here is the reason why I tthink something isn't right with my analysis: I built a phylogenetic tree in Mesquite, and based on several works I used million years as branch lengths. It makes sense for me because I'll be using a fossil on my analysis. But I tested the tree before adding the extinct taxa, so to make sure everything was OK when the tree was ultrametric. I noticed that, for example, whatever the independent variable (X) was, the alpha from OU was extremely high (0.999182, for instance). It would always be 0.999. I thought maybe I was using the wrong transformation... That's why I ended up trying to know if I needed to do something prior to transforming and analysing the tree. Thank you very much. Best regards, Sérgio. ᐧ 2015-04-13 17:17 GMT+01:00 Emmanuel Paradis emmanuel.para...@ird.fr mailto:emmanuel.para...@ird.fr: Hi Sérgio, There is indeed generally a relationship between branch length transformations and correlation structures. You may check that with the function vcv2phylo, e.g.: tr - rcoal(20) co - corGrafen(1, phy = tr) ts - vcv2phylo(vcv(co)) all.equal(tr, ts) [1] FALSE all.equal(compute.brlen(tr), ts) [1] TRUE compute.brlen() transforms the branch lengths according to Grafen's model with parameter rho = 1 (by default). Some other transformations of branch lengths are available in package geiger. Best, Emmanuel Le 12/04/2015 21:47, Sergio Ferreira Cardoso a écrit : Hi everyone, I'm relatively new in phylogenetic comparative methods. I'm a little confused about branch length transformations. I'm using a tree with divergence time (My) as branch lengths. When I use corPagel, corGrafen or corMartins in R, the branch lengths, are the branch lengths automatically transformed? e.g., gr.mammals-corGrafen(1,phylo,__fixed=F); fit-gls(FCL~logBodymass,__correlation=gr.mammals,data=__ df,method=ML). My question may sound a bit nonsense but I've seen in some papers (e.g., Spoor, F., Garland, T., Krovitz, G., Ryan, T. M., Silcox, M. T., Walker, A. (2007). The primate semicircular canal system and locomotion. *Proceedings of the National Academy of Sciences*, *104*(26), 10808-10812.) the indication that a PGLS was made without branch transformation, but no reference is made to the model (maybe it's corBrownian). Thank you very much. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
Re: [R-sig-phylo] PGLS transformations
Hello, So what you're suggesting is that I do that and then fit-gls(X~Y,corr=corBrownian(1,t.pb),data=df,method=ML, right? Thanks. Best regards, Sérgio. ᐧ 2015-04-13 18:42 GMT+01:00 Liam J. Revell liam.rev...@umb.edu: Hi all. I just wanted to add that it should be straightforward to sample edge lengths as if they arose under a specific process. For instance, here on my blog I show how to sample branching times edge lengths as if they arose under a pure-birth (i.e., 'Yule') process, given a particular input topology: http://blog.phytools.org/2015/04/sampling-edge-lengths- under-yule-process.html. For things like the phylogenetic (contrasts) regression, this may be preferable to Grafen's edge lengths, which tend to make recent edge lengths very short, giving very high weight to the associated contrasts. All the best, Liam Liam J. Revell, Assistant Professor of Biology University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://blog.phytools.org On 4/13/2015 1:13 PM, Sergio Ferreira Cardoso wrote: Hello, Thank you both for the help. Emmanuel, so is there a way to see contrasts in R? Reading this paper - Garland, T., Harvey, P. H., Ives, A. R. (1992). Procedures for the analysis of comparative data using phylogenetically independent contrasts. Systematic Biology, 41(1), 18-32. - I was aware of the importance of standardizing contrasts and of comparing Absolute value of standard contrat vs Standard deviation of contrast. I've learned to do this in Mesquite but I don't know if R alows this to be done. When I run the GLS I ask R to estimate the rho, so, a priori I never know the rho value. Maybe I'm really really confused, but here is the reason why I tthink something isn't right with my analysis: I built a phylogenetic tree in Mesquite, and based on several works I used million years as branch lengths. It makes sense for me because I'll be using a fossil on my analysis. But I tested the tree before adding the extinct taxa, so to make sure everything was OK when the tree was ultrametric. I noticed that, for example, whatever the independent variable (X) was, the alpha from OU was extremely high (0.999182, for instance). It would always be 0.999. I thought maybe I was using the wrong transformation... That's why I ended up trying to know if I needed to do something prior to transforming and analysing the tree. Thank you very much. Best regards, Sérgio. ᐧ 2015-04-13 17:17 GMT+01:00 Emmanuel Paradis emmanuel.para...@ird.fr: Hi Sérgio, There is indeed generally a relationship between branch length transformations and correlation structures. You may check that with the function vcv2phylo, e.g.: tr - rcoal(20) co - corGrafen(1, phy = tr) ts - vcv2phylo(vcv(co)) all.equal(tr, ts) [1] FALSE all.equal(compute.brlen(tr), ts) [1] TRUE compute.brlen() transforms the branch lengths according to Grafen's model with parameter rho = 1 (by default). Some other transformations of branch lengths are available in package geiger. Best, Emmanuel Le 12/04/2015 21:47, Sergio Ferreira Cardoso a écrit : Hi everyone, I'm relatively new in phylogenetic comparative methods. I'm a little confused about branch length transformations. I'm using a tree with divergence time (My) as branch lengths. When I use corPagel, corGrafen or corMartins in R, the branch lengths, are the branch lengths automatically transformed? e.g., gr.mammals-corGrafen(1,phylo,fixed=F); fit-gls(FCL~logBodymass,correlation=gr.mammals,data=df,method=ML). My question may sound a bit nonsense but I've seen in some papers (e.g., Spoor, F., Garland, T., Krovitz, G., Ryan, T. M., Silcox, M. T., Walker, A. (2007). The primate semicircular canal system and locomotion. *Proceedings of the National Academy of Sciences*, *104*(26), 10808-10812.) the indication that a PGLS was made without branch transformation, but no reference is made to the model (maybe it's corBrownian). Thank you very much. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] PGLS transformations
Hi everyone, I'm relatively new in phylogenetic comparative methods. I'm a little confused about branch length transformations. I'm using a tree with divergence time (My) as branch lengths. When I use corPagel, corGrafen or corMartins in R, the branch lengths, are the branch lengths automatically transformed? e.g., gr.mammals-corGrafen(1,phylo,fixed=F); fit-gls(FCL~logBodymass,correlation=gr.mammals,data=df,method=ML). My question may sound a bit nonsense but I've seen in some papers (e.g., Spoor, F., Garland, T., Krovitz, G., Ryan, T. M., Silcox, M. T., Walker, A. (2007). The primate semicircular canal system and locomotion. *Proceedings of the National Academy of Sciences*, *104*(26), 10808-10812.) the indication that a PGLS was made without branch transformation, but no reference is made to the model (maybe it's corBrownian). Thank you very much. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
[R-sig-phylo] PGLS - branch lengths 1
Dear all, I'm trying to perform a PGLS with arbitrary branch lengths (I used branch lengths = 1). The tree is non-ultrametric. Here is a result with Pagel's Lambda: vf-diag(vcv(tree)) fit-gls(FCLrelative~LogBodymass,correlation=Pagel,data=DF,method=ML,weights=varFixed(-vf)) summary(fit) Model: FCLrelative ~ LogBodymass Data: DF AIC BIClogLik 42.28777 49.77257 -17.14388 Correlation Structure: corPagel Formula: ~1 Parameter estimate(s): lambda 1.039835 Variance function: Structure: fixed weights Formula: ~vf Coefficients: Value Std.Errort-value p-value (Intercept)-0.17662926 0.26166738 -0.6750144 0.5030 LogBodymass 0.058657370.057661711.0172672 0.3143 Correlation: (Intr) LogBodymass -0.651 Standardized residuals: Min Q1 Med Q3 Max -2.36346640 -0.21848447 0.08326376 0.46156642 1.15440168 Residual standard error: 0.2388999 Degrees of freedom: 48 total; 46 residual As you see my lambda is 1, what shouldn't happen. Is this normal when we set all branch lengths to an arbitrary value? Thank you very much in advance. Best regards, Sérgio. -- Com os melhores cumprimentos, Sérgio Ferreira Cardoso. Best regards, Sérgio Ferreira Cardoso MSc. Paleontology candidate Departamento de Ciências da Terra - FCT /Universidade Nova de Lisboa Geociências - Universidade de Évora Lisboa, Portugal [[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/
