[R-sig-phylo] Running chronopl on a multiphylo object
Hello
I would like to use the ape function chronopl on a multiphylo object which
contains a subsample of phylogenies from a Bayesian analysis. I've set up a
loop which runs chronopl iteratively on the sample of trees from the multiphylo
object but I seem to be doing something wrong when trying to concatenate the
ultrametricized trees - output from chronopl - into a new multiphylo object and
write a nexus file to export for use in another program.
This is the loop:
trees-read.nexus(MYTREES.nex)
newmulti-list(0)
for(i in 1:500){
newphy-chronopl(trees[[i]], lambda=0.5)
newmulti-c(newphy, newmulti)}
Any help is welcome!!
Cheers
Alejandro
__
Alejandro Gonzalez Voyer
Post-doc
Estación Biológica de Doñana (CSIC)
Avenida Américo Vespucio s/n
41092 Sevilla
Spain
Tel: +34- 954 466700, ext 1749
E-mail: alejandro.gonza...@ebd.csic.es
Group page: http://consevol.org/people.html
Personal web-page: http://consevol.org/members/alejandro_combo.html
[[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] MRM(), AICc, and other ways to compare effects
Hi Roland, There are two functions for Rob's method to combine phylogenetic and spatial signal, although both - as far as I can tell - are based on gls methods, one uses independent contrasts, and thus I'm not sure if it is set up in such a way to allow you to include factors in the model. However, the other method is not based on independent contrasts but rather incorporates both matrices (the phylogenetic variance-covariance matrix and the geographical distance matrix) in the analyses, much as in typical PGLS analysis, therefore this method should - again if I am not mistaken - allow you to include factors. Cheers, Alejandro __ Alejandro Gonzalez Voyer Post-doc Estación Biológica de Doñana (CSIC) Avenida Américo Vespucio s/n 41092 Sevilla Spain Tel: +34- 954 466700, ext 1749 E-mail: alejandro.gonza...@ebd.csic.es Web-site (Under Construction): Group page: http://consevol.org/index.html Personal web-page: http://consevol.org/members/alejandro.html On 18, Sep 2012, at 4:48 PM, Luke Matthews wrote: Hi Roland, I haven't used Rob's method myself, but it might not be set up for factors. The lnam function isn't set up for factors either, but you should be able to use binary (0/1) variables in either method. You can perform the equivalent of adding a factor by coding it into separate 0/1 variables for the presence or absence of each factor level, leaving out one factor level because that will be your 'base' level. Then you add all these binary codings for the levels to the model at once. You can compare to the model without the full set of binary codings for levels to assess the difference in likelihood between including the factor in the model or leaving it out. Basically this is just having you do some of the steps that the machine does itself when you use a factor in a linear model. Don't do any messing about with taking only some binary level codings in and out, as then you would not be using them as a factor. You should only do that if you have a specific apriori hypothesis about a specific factor level, for example, Durkheim's famous hypothesis about reduced suicide among Catholics. You are correct that in lnam the factor variables (recoded into binary variables for #levels-1) should be entered as x. The phylogeny and geography correlation matrices should both be entered as W2 matrices. You use a syntax such as W2=list(phy.mat , geo.mat). The W1 matrix runs a qualitatively different autocorrelation used in some spatial and network research that models the distribution of the y variable itself rather than residuals of y. It is the distribution of residuals that you want to model. Best Luke -Original Message- From: Roland Sookias [mailto:r.sook...@gmail.com] Sent: Tuesday, September 18, 2012 3:28 AM To: Luke Matthews Cc: Alejandro Gonzalez Voyer; r-sig-phylo@r-project.org Subject: Re: [R-sig-phylo] MRM(), AICc, and other ways to compare effects Cheers Luke and Alejandro Yes, I've been trying to use Rob's method. However I can't see how it allows incorporation of the other variables (communism, religion) I want as they are categorical/factors. How would you go about including these? I am interested in trying the lnam function. However I can't quite follow from the documentation how to implement what I need. Would you put the factor variables communism and religion in as x and the phylogeny and geography matrices as W1 and W2? Cheers Roland On Wed, Sep 12, 2012 at 3:21 PM, Luke Matthews lmatth...@activatenetworks.net wrote: Hi Roland and Alejandro, That is good to hear that Rob has the code already for his method to toggle the various matrices in and out to then assess likelihoods. Roland, you can also do this with the 'lnam' function from 'sna' package. You can fit your models with any combination of autocorrelation matrices and regular independent variables. It might be methodologically interesting to see what using both Freckleton's code and the 'lnam' code by Carter Butts produce, as the authors I'm sure don't know each other or each other's work. It would be neat to get some sense of how these implementations compare given that surely there are differences in the implementations even though they are essentially trying to do the same thing. Best Luke -Original Message- From: Alejandro Gonzalez Voyer [mailto:alejandro.gonza...@ebd.csic.es] Sent: Wednesday, September 12, 2012 2:23 AM To: Roland Sookias Cc: Luke Matthews; r-sig-phylo@r-project.org Subject: Re: [R-sig-phylo] MRM(), AICc, and other ways to compare effects Hi, Yes this is possible, look at Freckleton and Jetz 2009. Rob has R code for two functions allowing to combine geographic and phylogenetic matrices. Cheers Alejandro _ Dr. Alejandro Gonzalez Voyer Estación
Re: [R-sig-phylo] Non-parametric alternative to phylogenetic ANOVA?
Hi Karin, I would suggest setting all branch lengths to a unit value and then using corPagel() in ape as the model of trait evolution. Using corPagel the value of lambda is estimated simultaneously with model fit, lambda is an estimate of the phylogenetic signal in your residuals and basically adjusts the covariance among species so that the residuals fit a brownian model (Freckleton et al 2002 and Revell 2010 provide much more detailed explanations of all this), so it should work well even with equal branch lengths for your tree. You could also do as you suggested and transform branch lengths as suggested by Grafen and then repeat the analyses using corPagel as test of how sensitive results are to different branch length transformations. Based on the help files in ape, I understand corGrafen transforms branch lengths first following Grafen's method and then uses the traditional variance-covariance matrix for the phylogeny, which I understand to be a brownian model (ie corBrownian). Hope this helps. Cheers Alejandro __ Alejandro Gonzalez Voyer Post-doc Estación Biológica de Doñana (CSIC) Avenida Américo Vespucio s/n 41092 Sevilla Spain Tel: +34- 954 466700, ext 1749 E-mail: alejandro.gonza...@ebd.csic.es Web-site (Under Construction): Group page: http://consevol.org/index.html Personal web-page: http://consevol.org/members/alejandro.html On 31, May 2012, at 9:46 AM, Karin Schneeberger wrote: Dear Liam, dear Alejandro Thank you so much for your answers, that was very helpful! Just (hopefully) one last question: As branch lengths are not known for the tree I am using, I calculated them as suggested by Grafen (branch lengths being the number of descending taxa minus 1). So, I would replace corBrownian by corGrafen, proceed as suggested by testing the residuals for normality and conduct post-hoc tests if they don't differ from normal distribution, right? Thank you so much for your help. Best regards, Karin Von: Liam J. Revell liam.rev...@umb.edu An: Karin Schneeberger k.schneeber...@yahoo.de CC: Alejandro Gonzalez alejandro.gonza...@ebd.csic.es; r-sig-phylo@r-project.org r-sig-phylo@r-project.org Gesendet: 17:09 Mittwoch, 30.Mai 2012 Betreff: Re: [R-sig-phylo] Non-parametric alternative to phylogenetic ANOVA? Hi Karin. GLS with x as a factor is a generalized ANOVA which assumes [in the case of gls(...,correlation=corBrownian)] that the residual error in the ANOVA model has evolved by Brownian evolution. If you read your data into data frame Z with row names as species names, for instance: Z-read.table(filename,header=T,row.names=1) tree-read.tree(treefile) and your column name for the factor is x the column name for the continuous response variable is y, then you should just be able to do: fit-gls(y~x,data=Z,correlation=corBrownian(1,tree)) You can then perform various posthoc analyses from the gls object that is produced. For instance summary(fit) anova(fit) residuals(fit) As pointed out by Alejandro, you should check for normality of the residuals in residuals(fit) - not the normality of y before analysis. summary(fit) will also give you parameter estimated (fitted means for each factor) and standard errors. These can be used to conduct posthoc comparison of means using t-tests in the standard way. I hope this helps. All the best, Liam -- Liam J. Revell University of Massachusetts Boston web: http://faculty.umb.edu/liam.revell/ email: liam.rev...@umb.edu blog: http://phytools.blogspot.com On 5/30/2012 10:46 AM, Karin Schneeberger wrote: Hi Alejandro Thank you for the very quick answer. I tried PGLS before, but then was told that GLS is not suitable for multistate categorical variables that can not be ranked (otherwise I would treat them as continuous). Also, with GLS it's as far as I understood not possible to state statistically whether certain groups are greater than others. But I am new into this kind of analysis and am very happy for any help and explanation, as I might be totally wrong. Cheers, Karin Von: Alejandro Gonzalezalejandro.gonza...@ebd.csic.es CC: r-sig-phylo@r-project.orgr-sig-phylo@r-project.org Gesendet: 16:26 Mittwoch, 30.Mai 2012 Betreff: Re: [R-sig-phylo] Non-parametric alternative to phylogenetic ANOVA? Hi Karin, You could use a gls method and look at the distribution of your residuals. It is the residuals which must be normally distributed, which can be checked using diagnostic plots such as a histogram or qq-plot of the residuals of your model. Cheers Alejandro On 30, May 2012, at 4:12 PM, Karin Schneeberger wrote: Dear all I'm trying to compare one trait across three (unordered categorical) groups including 25 species (let's say for example basal metabolic rate of aquatic, terrestrial and aerial
[R-sig-phylo] Estimating branch lengths on a fixed topology from molecular data
Hello, I was wondering if it is possible to estimate branch lengths on a fixed topology from a matrix of nuclear or mitochondrial DNA data. I used to do this in PAUP and would very much like to have an alternative in R. In PAUP one loaded a matrix of molecular data and a topology (without branch lengths) with identical species coverage, then branch lengths were estimated using maximum likelihood (could also use parsimony) and given settings for the substitution model. Any tips would be appreciated! Best wishes Alejandro __ Alejandro Gonzalez Voyer Post-doc Estación Biológica de Doñana (CSIC) Avenida Américo Vespucio s/n 41092 Sevilla Spain Tel: +34- 954 466700, ext 1749 E-mail: alejandro.gonza...@ebd.csic.es Web-site (Under Construction): Group page: http://consevol.org/index.html Personal web-page: http://consevol.org/members/alejandro.html [[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] R-sig-phylo list
Here is the Email: R phylo mailing list mailing list r-sig-phylo@r-project.org Cheers Alex __ Alejandro Gonzalez Voyer Post-doc Estación Biológica de Doñana (CSIC) Avenida Américo Vespucio s/n 41092 Sevilla Spain Tel: +34- 954 466700, ext 1749 E-mail: alejandro.gonza...@ebd.csic.es Web-site: https://docs.google.com/View?id=dfs328dh_14gwwqsxcg [[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] multi-state categorical predictor variables in PGLS
Hi Andrew, As I understand it, gls() is doing a multiple generalized LS regression with as many dummy variables as there are factor levels. Is this a correct characterization? I think you'd get one dummy variable less than factor levels in your characterization (at least in regards to the number of levels for which parameters are estimated), as the gls sets one of the levels as the point of comparison with all other levels. Thus you'd get n-1 dummy variables for which the parameters are estimated. Having such a low value of alpha the results of the phylogenetic gls should be similar (if not identical) to results not taking phylogeny into account, as this suggests you don't have phylogenetic signal in the residuals of your relationship. There is a good paper on this issue by Liam Revell in Methods in Ecology and Evolution. There is also a function in geiger which allows you to run phylogenetic ANOVAs, but if I am not mistaken, the p-value is estimated based on simulations assuming traits evolve via Brownian motion (is this correct?). I've also seen lambda values below 0 in ape, theoretically lambda is described as being bounded between 0 and 1, but it could take values outside the bounds. I would be interested in hearing the thoughts of others in the list regarding whether lambda values for the Phylogenetic gls should be forced to be bounded between 0 and 1. This would more closely follow what has been proposed in the literature wouldn't it? Cheers Alejandro __ Alejandro Gonzalez Voyer Post-doc Estación Biológica de Doñana (CSIC) Avenida Américo Vespucio s/n 41092 Sevilla Spain Tel: +34- 954 466700, ext 1749 E-mail: alejandro.gonza...@ebd.csic.es Web-site: https://docs.google.com/View?id=dfs328dh_14gwwqsxcg On 7, Mar 2011, at 5:52 PM, Andrew Barr wrote: Hi everyone, I am trying to piece together the current best-practices for phylogenetic ANOVA with multi-state predictors. In my dataset, my four-level factor is non-random with respect to phylogeny. That is, if I know which higher level clade an species belongs to, I can predict with pretty good success which factor level it will be in. My understanding is that this situation likely overinflates my degrees of freedom and makes traditional F-tests inappropriate. I came across this paper (Garland et al 1993. Phylogenetic Analysis of Covariance by Computer Simulation. Systematic Biology 42:265 -292.) where the authors empirically recalculate critical values for F-ratios using computer simulations, tree topology, and a model of character evolution. I also have found that I can use PGLS (with ape and nlme) and specify my model like this. gls(myVar~myFactor,corr=corPagel(val=1,phy=myTree,fixed=F),data=myDF) As I understand it, gls() is doing a multiple generalized LS regression with as many dummy variables as there are factor levels. Is this a correct characterization? Does this sidestep the degrees of freedom problem discussed by Garland et al.? Can anybody point me to references discussing the mechanics of this process and why this is an appropriate thing to do? Finally, I get a negative value for estimated lambda. Any ideas on what that means? Thanks to everyone for any advice/references/. Andrew Barr PhD Student University of Texas at Austin results from my model Generalized least squares fit by REML Model: LIWI ~ Hab Data: aggast AIC BIC logLik -65.61627 -56.28418 38.80814 Correlation Structure: corPagel Formula: ~1 Parameter estimate(s): lambda -0.1480891 Coefficients: Value Std.Error t-value p-value (Intercept) 1.4492742 0.01876415 77.23635 0. HabH-0.0224975 0.03149986 -0.71421 0.4798 HabL-0.0668761 0.03066232 -2.18105 0.0360 HabO-0.1630386 0.02567505 -6.35008 0. Correlation: (Intr) HabH HabL HabH -0.686 HabL -0.794 0.485 HabO -0.936 0.594 0.542 Standardized residuals: Min Q1 Med Q3 Max -2.17865325 -0.60297897 -0.09760938 0.41995284 2.91201671 Residual standard error: 0.06913702 Degrees of freedom: 39 total; 35 residual ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo [[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] Phylogenetic ANOVA
Hello, Some colleagues and I are running some phylogenetic ANOVAS using the geiger package. In some of the analyses we get the same phylogentic p-value (very small p-value) even though the F-statistic differs between the two analyses, albeit it being relatively high in both instances. We were wondering why this arises, to get better grip on how the analysis works. We thought it may have to do with the randomizations to calculate the phylogenetic p-value. Or that the F-statistics are quite high... Below are two examples : m11-phy.anova(tree1,tmax,biozone,data.names=X,nsim=1000) Standard ANOVA: Analysis of Variance Table Response: td$data Df Sum Sq Mean Sq F valuePr(F) group 1 967.96 967.96 155.88 3.057e-12 *** Residuals 25 155.246.21 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Phylogenetic p-value: 0.000999001 m12-phy.anova(tree1,wt,biozone,data.names=X,nsim=1000) Standard ANOVA: Analysis of Variance Table Response: td$data Df Sum Sq Mean Sq F valuePr(F) group 1 602.88 602.88 109.01 1.333e-10 *** Residuals 25 138.265.53 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Phylogenetic p-value: 0.000999001 Cheers, Alejandro __ Alejandro Gonzalez Voyer Post-doc NEW ADDRESS NEW E-MAIL Estación Biológica de Doñana (CSIC) Avenida Américo Vespucio s/n 41092 Sevilla Spain E-mail: alejandro.gonza...@ebd.csic.es Tel: +34- 954 466700, ext 1749 Website (From my previous position): http://www.iee.uu.se/zooekol/default.php?type=personalpagelang=enid=146 [[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] Estimating parameters of OU model using OUCH
Hello, I am interested in using the package OUCH to estimate parameters of an OU model of evolution for a suite of traits. I've been able to run the fitContinuous function in geiger without problems, but have hit an obstacle when it comes to OUCH. I think it probably has to do with how the data need to be organized in OUCH. I would greatly appreciate any tips to run the hansen function in OUCH, especially regarding data organization. Cheers Alejandro __ Alejandro Gonzalez Voyer Post-doc NEW ADDRESS NEW E-MAIL Estación Biológica de Doñana (CSIC) Avenida Américo Vespucio s/n 41092 Sevilla Spain E-mail: alejandro.gonza...@ebd.csic.es Tel: +34- 954 466700, ext 1749 Website (From my previous position): http://www.iee.uu.se/zooekol/default.php?type=personalpagelang=enid=146 [[alternative HTML version deleted]] ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
