# ape comparative analysis query

3 messages
Open this post in threaded view
|
Report Content as Inappropriate

## ape comparative analysis query

 I've been comparing variables among objects (taxa) related by known trees, using phylogentically independent contrasts in the ape package, and want to move on to more complex models e.g. by using gls with appropriate correlation terms. My trees contain lots of (hard) polytomies and information about ancestors, which I've been including- creating fully dichotomous trees by using zero branch lengths. For instance: library(ape) tree<-read.tree(text="((B1:3,(((D1:5,C1:0):5,B2:0):4,(B3:3,B4:4):0):0):0,A1:0):0;") plot(tree) Where all the B taxa arise in a true polytomy from A1 and B2 is the ancestor of C1 which is the ancestor of D1 I have information on all the taxa and can create independent contrasts for this: x<-c(B1=47,D1=43,C1=45,B2=50,B3=47,B4=48,A1=48) y<-c(B1=2.9,D1=5.4,C1=2.8,B2=3.5,B3=3.2,B4=3.5,A1=1.8) picx<-pic(x,phy=tree) picy<-pic(y,phy=tree) Which seems to make sense. However, if I try to use anything more sophisticated I run into problems: >dat<-as.data.frame(cbind(x,y)) >model1<-gls(y~x,data=dat, correlation=corBrownian(phy=tree)) Error in corFactor.corStruct(object) : NA/NaN/Inf in foreign function call (arg 1) or equally: > model2<-compar.gee(y~x,data=dat, phy=tree) [1] "Beginning Cgee S-function, @(#) geeformula.q 4.13 98/01/27" [1] "running glm to get initial regression estimate" [1] 15.6000 -0.2625 Error in gee(y ~ x, c(1, 1, 1, 1, 1, 1, 1), data = list(x = c(47, 43,  :         NA/NaN/Inf in foreign function call (arg 14) (NB using the 'data=' argument seems to be necessary in the gls case- if I don't, there is a further problem: 'Row names in dataframe do not match tree tip names. data taken to be in the same order as in tree. in: Initialize.corPhyl(X[[1]], ...)') This seems to go away if I remove the root taxon (A1): tree2<-drop.tip(tree, "A1") x2<-c(B1=47,D1=43,C1=45,B2=50,B3=47,B4=48) y2<-c(B1=2.9,D1=5.4,C1=2.8,B2=3.5,B3=3.2,B4=3.5) dat2<-as.data.frame(cbind(x2,y2)) model2<-gls(y2~x2,data=dat2, correlation=corBrownian(phy=tree2)) This raises various questions: 1) Was I misleading myself that my independent contrasts were valid in the first place? 2) What is it, if anything, about the root taxon that causes this issue, given that other taxa also have zero branch lengths? 3) Is there any way of getting around this and including data on the root taxon, or am I better off just dropping it (ultimately I want to work with much larger trees (up to tens of thousands of taxa) where that one piece of information will become relatively less important) Any help very much appreciated, Chris I'm working with ape 1.8-2 in R 2.1.1 under ubuntu 'Breezy' linux (unfortunately 2.1.1 is the latest easily available in breezy) -- -------------------------------------------------------------- Dr Christopher Knight                      School of Chemistry [hidden email]     The University of Manchester Tel: extn 64414                               Faraday Building      +44 (0)161 3064414                              PO Box 88 Fax: +44 (0)161 3064556                       Sackville Street                                             Manchester M60 1QD  ` · . ,,><(((°>                                            UK ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide! http://www.R-project.org/posting-guide.html
Open this post in threaded view
|
Report Content as Inappropriate

## Re: ape comparative analysis query

 Chris Knight wrote: > > This raises various questions: > > 1) Was I misleading myself that my independent contrasts were valid in > the first place? >   I think partly. Independent contrasts were designed for data for tip taxa only, not ancestral state data. One of the steps in the algorithm requires a correction to account for the fact that ordinarily the values for traits at nodes are estimated and not known, basically lengthening the branch length to that node by a certain amount (See Felsenstein's original paper). If you have known ancestral states, then you should not make this correction. > 2) What is it, if anything, about the root taxon that causes this issue, > given that other taxa also have zero branch lengths? >   When converting a tree to a correlation matrix for use in GLS or GEE, the Brownian motion assumption implies that the covariances among taxa are represented by their shared branch length from the root, and the height of each terminal taxon is equal to the variance of the trait for that taxon. If the tree is ultrametric, then the variances are equal. In your case, the root has zero covariance with each taxon, and zero variance. Now, the error in the gls call is caused by a division by zero in the calculation of the correlation matrix, because the root has zero variance. try: vcv.phylo(tree) vcv.phylo(tree, cor=TRUE) > 3) Is there any way of getting around this and including data on the > root taxon, or am I better off just dropping it (ultimately I want to > work with much larger trees (up to tens of thousands of taxa) where that > one piece of information will become relatively less important) >   I'm surprised that you have known values for the root data. Are you sure that the root taxon is not actually an outgroup? If so, it will have some branch length (variance), and the model should fit OK. HTH, Simon. > Any help very much appreciated, > > Chris > > I'm working with ape 1.8-2 in R 2.1.1 under ubuntu 'Breezy' linux > (unfortunately 2.1.1 is the latest easily available in breezy) > > ______________________________________________ [hidden email] mailing list https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide! http://www.R-project.org/posting-guide.html