[R-sig-phylo] R: Re: From ClustalW2 Tree to Heat Map in R
Hi Gents and Lads, I have a very rapid question with a perhaps not-so-obvious reply. I'm in the process of testing a number of evolutionary models and estimating phylogenetic signal on a certain univariate data set. So something very basic and very simple. The point is, I found BM performs poorly as compared to OU (single peak) and lambda. Thus, I transformed the tree by using the fitted lamdba and/or alpha before calculating Blomberg et al's K statistic. Does this make sense? If yes, the competing models have similar Akaike weigths (OU = 0.5, lambda = 0.3) but give very different estimates of K when their fitted parameters are used to transform the tree branch lenghts. How does discriminate which K estimate is best? Translating in R-esque: require(geiger) require(phytools) fitContinuous(tree, x, model="OU") ## gives relatively low alpha (.07) phylosig(ouTree(tree,alpha=alpha),x,method ="K", test =T, nsim =1000) ## gives very low K fitContinuous(tree, x, model="lambda") ## gives very high lambda (.95) phylosig(lambdaTree(tree,lambda=lambda),x,method ="K", test =T, nsim =1000) ## gives very high K thanks for help, Pas ___ 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] mccr test for non ultrametric trees
Hi all, I would enquire the list about a simple issue. Is there any method implemented to test for changes in diversification rate as applied to fossil (non- ultrametric) trees?. As far as I understand, the methods so far available work on ultrametric trees (I've inspected laser's mccrTest.Rd, and Cusimano's et al. CorSiM). I'm tempted to say that an easy workaround would be to get the real age of nodes by substituting the function branching.times with an an hoc script where needed, and then simulate trees of the same size as the orginal (e.g. with TreeSim functions) to get the null distribution of gamma values. Is that feasible? thanks, Pas ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] R: Comparative Methods and Pseudo-Traits
Hi all, David, I'd say yes, PICs or P-GLS or any other? related PCM could be applied to quasi-traits in principle. On the mathematical grounds, I think most of us will be content with Liam's explanation. On the theoretical grounds, I think that for most quasi-traits the correct answer is in your own words: hidden correlations may justify PCMs as applied to quasi traits. Let's take the growth rate/habitat degradation (by humans) relationship. If you assume similar/close habitats are similarly affected by humans, and take niche conservatism seriously, than you will find yourself assuming that closely related species most probably experience similar habitat conditions, or degradation thereof. This means that you'd expect that the quasi trait "habitat degradation by anthropogenic disturbance" takes closer than expected values among closely related species, which is what PCMs do test or account for. One way to look at evolutionary models is that they describe the path of evolution in the simplest yet most reasonable way. Another way to look at them, is that they provide a basic, quantitative hypothesis about character divergence given the (phylogenetic) history of the species you consider. In the latter sense, to appy PICs to a quasi-trait is much easier to swallow, besides being mathematically feasible. Pas >Messaggio originale >Da: dwba...@uchicago.edu >Data: 10/11/2011 21.36 >A: "R Sig Phylo Listserv" >Ogg: [R-sig-phylo] Comparative Methods and Pseudo-Traits > >Hello all, >A recent discussion set my mind thinking on a particular issue and, once >again, I decided to ask for the general opinion of R-Sig-Phylo denizens. It >may be easier to start with an example. > >Let's say that there exists a worker who is measuring several different >traits across a number of species and then testing for correlations among >these traits. The first test is body size versus growth rate and they use >independent contrasts or PGLS to test for a the correlation, accounting for >phylogeny. Both of these traits are inherited, evolving variables. Now >let's say they'd like to test for the relationship between growth rate and >some metric of the anthropogenic degradation of that species' habitat. Now >what? It is even valid to apply PIC to the habitat degradation metric even >though it is not an inherited, evolving trait? It's unclear to me. > >Let's consider a paleontological example, one which I have found myself >both strongly agreeing and disagreeing with at times. Essentially, how >should we test for extinction selectivity on some trait at a mass >extinction event? Let's say we think body size is a predictor of the risk >of extinction during that event and so we want to test for a correlation >between them (please ignore that extinction would be a discrete variable >for the moment). Do we treat these variable with PIC or PGLS? Is it really >proper to refer to the probability of going extinct during a mass >extinction as an evolving trait? Let's say we did and we got different >results than when we used an analysis which did not account for the >phylogenetic covariance. How should we interpret these results? > >One explanation I know of is that when we apply phylogenetic comparative >methods to these quasi-traits to consider their relationship to another >trait, we are assuming that these variables are actually the result of some >underlying, unobserved set of traits which are evolving along the >phylogeny. This makes sense, maybe in the extinction event case, which >would mean that any PCM analysis would be testing for an evolutionary >relationship between body size and these unobserved traits which predict >extinction. Of course, if extinction risk is largely a function of >non-inherited traits, then the initial assumption may be incorrect (that >extinction risk itself is an evolving trait). Regardless, I don't see how >to apply that explanation to the habitat degradation example. > >So, what do people think? How should we test for correlation when >non-evolving quasi-traits are involved? I'm very interested to hear >people's thoughts on this matter. >-Dave Bapst, UChicago > >-- >David Bapst >Dept of Geophysical Sciences >University of Chicago >5734 S. Ellis >Chicago, IL 60637 >http://home.uchicago.edu/~dwbapst/ > > [[alternative HTML version deleted]] > >___ >R-sig-phylo mailing list >R-sig-phylo@r-project.org >https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] R: Re: R: Re: R: Re: R: ancestral state reconstruction for tips
Folks, I was intending my most recent message to be off-list and didn't realize "r-sig-phylo@r-project.org" was in the CC field, which means I'm a fool. All kidding aside, yes Joe, I estimated ancestral (and tip) values for which I have real data via BM assumption to see how good the fit is. Actually, estimated values are very close to real values for some species, barely so for some others, and absolutely not for others still. The good news is that since there is a single mode of evolution tree wise, deviations from real values really mean that evolution is accelerated, or decelerated, either, in these particular lineages for which a significant deviation from the expected value is noticeable. What I', trying to do now is writing a R routine to back-calculate the "expected" branch lengths for the "unusual" critters, given the fitted ancestral values and tip values of the phenotypes, and assuming BM, in order to compare the actual branch lengths to the expected. The ratio of these ! lengths, if I'm not delusional and definitely lucky, is a per-lineage rate of phenotypic evolution. The bottom line is answering the question: how long should the branch leading to that particular species be if it evolved at the same rate of its sister species?Pas Messaggio originale Da: j...@gs.washington.edu Data: 05/08/2011 21.04 A: "pasquale.r...@libero.it" Cc: , , , "r-sig-phylo@r-project.org" Ogg: Re: R: Re: [R-sig-phylo] R: Re: R: ancestral state reconstruction for tips Folks -- I was intending my most recent message to be apologetic --that I was perhaps overreactive. Certainly Pas has notraised unreasonable objections or been obstructive withmy grants! (Others have). Let me raise an issue so I understand him more clearly:Pas, are you saying that you see phenotypes in the fossilsthat seem incompatible with the Brownian Motion assumption? JoeJoe Felsenstein j...@gs.washington.edu Dept of Genome Sciences and Dept of Biology, Univ. of Washington, Box 5065, Seattle Wa 98195-5065 [[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: Re: R: Re: R: ancestral state reconstruction for tips
Hi All, I'm happy I have stimulated some discussion about this subject matter. For some reason I can't imagine it looks this whole thing is going to be somehow personal and I have not posted this last e-mail to the list as a consequence. Joe, unfotunately I never attended a lecture of yours, and didn't raise trivial distinctions and objections to a grant proposal you submitted. My intention was not to be critical about BM or ICs, or whatever. I just wanted to point it out that things are sometimes a bit too complex and some unreliable predictions from our models may slip out unnoticed evey now and then, as I believe it is apparent reading the literature (including my own, of course). Having said that, my guess was that we *may* use the BM and computations at nodes to see where (in which lineages) do phenotypes appear very different from predictions. For instance, I think it could be somehow possible to use estimated ancestral charactes to see how much the inclusion of some fossil (or new) species changes the estimated value (e.g. by creating a polytomy by the inclusion of the new species), or even back-calculate branch lenghts (under BM assumption) for these unusual phenotypes to see how much evolution accelerates in these lineages (by comparison with real branch lengths). I hope I spoke my mind more clearly at this time. Pas >Messaggio originale >Da: dwba...@uchicago.edu >Data: 05/08/2011 20.23 >A: "Joe Felsenstein" >Cc: "r-sig-phylo@r-project.org" >Ogg: Re: [R-sig-phylo] R: Re: R: ancestral state reconstruction for tips > >As the diversity of explicit models of trait evolution grow, it will >be interesting to see if any consensus develops about which models >hold most often in general and whether any insight is gained into >which conditions predict appearance of different models. > >I think Joe is right that realizing a model is an inaccurate or >imprecise description of reality should impel us to develop better >models of the world around us, because this partly how science moves >forward. However, I don't think pointing out that a model is deficient >requires that that person must themselves develop an alternative. >After all, an alternative model that capture a more realistic level of >complexity may not be possible in some situations (it is certainly >possible in trait evolution models, however.) Requiring such a thing >would put too much pressure on scientific whistle-blowers, who play a >very important role in reminding the rest of us that the world is more >than the models we use to understand it and make our predictions. > >-Dave > > > > >On Fri, Aug 5, 2011 at 10:51 AM, Joe Felsenstein wrote: >> >> Pasquale Raia said: >> >>> Of course Ted is right, but my problem with this computation, or >>> with the >>> simple exercise I was proposing is well another: as a >>> paleontologist I often >>> come across pretty exceptional phenotypes (dwarf hippos and >>> elephants, huge >>> flightless birds, to make a few examples). When you use methods >>> like this (I >>> mean Garland and Ives') and compare the output with those >>> phenotypes, as I did, >>> you immediately realize what the the bottom line is: no matter if >>> they are >>> nodes or tips, by using the expected (under BM) covariance the >>> estimated >>> phenotypes are dull, perfectly reasonable but very different from >>> anything >>> exceptional you may find yourself to work with. This is why I feel >>> it is >>> difficult to rely on those (unobserved) values to begin with. >> >> I think that what is being said is that Brownian Motion is too sedate >> a process >> and does not predict some of the large changes actually seen in the >> fossil >> record. >> >> That's a legitimate point but does put the onus on the maker of the >> point to >> propose some other stochastic process that is tractable and has these >> large >> changes (and that fits with known Mendelian and Darwinian mechanisms). >> Just complaining that the Brownian stochastic process is no good is >> insufficient. >> >> If we want to add the fossils to the calculation, then they will of >> course >> pressure the Brownian Motion process to change more in their vicinity, >> which may help some. >> >> Joe >> >> Joe Felsenstein j...@gs.washington.edu >> Dept of Genome Sciences and Dept of Biology, Univ. of Washington, >> Box 5065, Seattle Wa 98195-5065 >> >> >> [[alternative HTML version deleted]] >> >> ___ >> R-sig-phylo mailing list >> R-sig-phylo@r-project.org >> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo >> > > > >-- >David Bapst >Dept of Geophysical Sciences >University of Chicago >5734 S. Ellis >Chicago, IL 60637 >http://home.uchicago.edu/~dwbapst/ > >___ >R-sig-phylo mailing list >R-sig-phylo@r-project.org >https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > ___ R-sig-phylo mailing list R-sig-
[R-sig-phylo] R: Re: R: ancestral state reconstruction for tips
Of course Ted is right, but my problem with this computation, or with the simple exercise I was proposing is well another: as a paleontologist I often come across pretty exceptional phenotypes (dwarf hippos and elephants, huge flightless birds, to make a few examples). When you use methods like this (I mean Garland and Ives') and compare the output with those phenotypes, as I did, you immediately realize what the the bottom line is: no matter if they are nodes or tips, by using the expected (under BM) covariance the estimated phenotypes are dull, perfectly reasonable but very different from anything exceptional you may find yourself to work with. This is why I feel it is difficult to rely on those (unobserved) values to begin with. Any opinion? Pas >Messaggio originale >Da: theodore.garl...@ucr.edu >Data: 05/08/2011 18.24 >A: "Hunt, Gene", "r-sig-phylo@r-project.org" >Ogg: Re: [R-sig-phylo] R: ancestral state reconstruction for tips > >The methods in the Garland and Ives (2000) paper are in our package of DOS PDAP programs, and should also be functional in the PDAP module of Mesquite. > >Cheers, >Ted > >Theodore Garland, Jr. >Professor >Department of Biology >University of California, Riverside >Riverside, CA 92521 >Office Phone: (951) 827-3524 >Home Phone: (951) 328-0820 >Facsimile: (951) 827-4286 = Dept. office (not confidential) >Email: tgarl...@ucr.edu >http://www.biology.ucr.edu/people/faculty/Garland.html > >Experimental Evolution: Concepts, Methods, and Applications of Selection Experiments >Edited by Theodore Garland, Jr. and Michael R. Rose >http://www.ucpress.edu/book.php?isbn=9780520261808 >(PDFs of chapters are available from me or from the individual authors) > > >From: r-sig-phylo-boun...@r-project.org [r-sig-phylo-boun...@r-project.org] on behalf of Hunt, Gene [hu...@si.edu] >Sent: Friday, August 05, 2011 8:35 AM >To: r-sig-phylo@r-project.org >Subject: Re: [R-sig-phylo] R: ancestral state reconstruction for tips > >Also, the issue of predicting values for unknown tips using data from other species in the tree is considered in this reference: > >Garland, T., and A. R. Ives. 2000. Using the past to predict the present: confidence intervals for regression equations in phylogenetic comparative methods. American Naturalist 155(3):346-364. > >Best, >Gene > > > >On 8/5/11 11:31 AM, "pasquale.r...@libero.it" wrote: > > > > >Hi Morgan, > >this is just stuff for thought, and remember, this is wrong anyway. But you >may try something like this: > >1. compute pics, >2. take the pic value at the ancestral node subtending to your unknown tip, >3. pretend one of the two tips the pic was originally computed on is in fact >your unknown species, >4. modify the square of the summed branch lengths of the two species by using >the "new" bl, >5. use the formula for pics (standardized) to derive your unknown tip value by >using the other (real) species tip value and the new square of summed branch >lengthts > > >but again, remember this is wrong, because contrasts were computed without >your unknown species. With ace everything turns out to be much more complicated >because ancestral value estimations are 'optimized' by taking the entire tree >and distribution of values at once, so to speak. > > > > > >>Messaggio originale >>Da: morgan.g.i.langi...@gmail.com >>Data: 05/08/2011 14.15 >>A: >>Ogg: [R-sig-phylo] ancestral state reconstruction for tips >> >>I was wondering if there is a way to get ancestral state >>reconstructions not for nodes within the tree but for tips that I >>don't know the trait of. I could do this somewhat manually, by taking >>the ancestral state resconstruction from the parent and child nodes >>surrounding where my unknown tip branches off from the tree and >>averaging those results (weighted by the branch length). This approach >>seems kind of clunky, so I was hoping there was something better. >> >> >> >>Morgan Langille >> >>___ >>R-sig-phylo mailing list >>R-sig-phylo@r-project.org >>https://stat.ethz.ch/mailman/listinfo/r-sig-phylo >> > >___ >R-sig-phylo mailing list >R-sig-phylo@r-project.org >https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > > > >-- >Gene Hunt >Curator, Department of Paleobiology >National Museum of Natural History >Smithsonian Institution [NHB, MRC 121] >P.O. Box 37012 >Washington DC 20013-7012 >Phone: 202-633-1331 Fax: 202-786-2832 &g
[R-sig-phylo] R: ancestral state reconstruction for tips
Hi Morgan, this is just stuff for thought, and remember, this is wrong anyway. But you may try something like this: 1. compute pics, 2. take the pic value at the ancestral node subtending to your unknown tip, 3. pretend one of the two tips the pic was originally computed on is in fact your unknown species, 4. modify the square of the summed branch lengths of the two species by using the "new" bl, 5. use the formula for pics (standardized) to derive your unknown tip value by using the other (real) species tip value and the new square of summed branch lengthts but again, remember this is wrong, because contrasts were computed without your unknown species. With ace everything turns out to be much more complicated because ancestral value estimations are 'optimized' by taking the entire tree and distribution of values at once, so to speak. >Messaggio originale >Da: morgan.g.i.langi...@gmail.com >Data: 05/08/2011 14.15 >A: >Ogg: [R-sig-phylo] ancestral state reconstruction for tips > >I was wondering if there is a way to get ancestral state >reconstructions not for nodes within the tree but for tips that I >don't know the trait of. I could do this somewhat manually, by taking >the ancestral state resconstruction from the parent and child nodes >surrounding where my unknown tip branches off from the tree and >averaging those results (weighted by the branch length). This approach >seems kind of clunky, so I was hoping there was something better. > > > >Morgan Langille > >___ >R-sig-phylo mailing list >R-sig-phylo@r-project.org >https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] R: ancestral state reconstruction for tips
Hi Morgan, I believe this is theoretically possible but not advisable. You should take your unknown tips off the tree, estimate ancestral characters, and then compute values at (unknown) tips by using their branch lenghts and taking the inverse of the ACE formula. Yet, this would assume AC to be perfect (estimated with no error), which is simply untrue. This would much be like applying simple regression to estimate unknown Y values for Xs outside the range of Xs the slope was originally computed for. Pas >Messaggio originale >Da: morgan.g.i.langi...@gmail.com >Data: 05/08/2011 14.15 >A: >Ogg: [R-sig-phylo] ancestral state reconstruction for tips > >I was wondering if there is a way to get ancestral state >reconstructions not for nodes within the tree but for tips that I >don't know the trait of. I could do this somewhat manually, by taking >the ancestral state resconstruction from the parent and child nodes >surrounding where my unknown tip branches off from the tree and >averaging those results (weighted by the branch length). This approach >seems kind of clunky, so I was hoping there was something better. > > > >Morgan Langille > >___ >R-sig-phylo mailing list >R-sig-phylo@r-project.org >https://stat.ethz.ch/mailman/listinfo/r-sig-phylo > ___ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] R: nodes and taxa depth
Hi Paolo,
I'm sure there is something better (and possibly simpler) around. Anyway, the
script appended below retrieves from the tree info such as node identity, age,
branch lengths and so on. There are two outputs you might be interested into:
"node.ages" gives the above info for all species and nodes in the tree.
"species.ages" is just "node.ages" limited to the tips.
Hope it may help,
Pas
library(ape)
library(picante)
library(geiger)
phy<-read.tree("mytree.txt")
node.age(phy)->phy.age
cbind(phy.age$edge,phy.age$age, phy$edge.length)->BL.position
max(phy.age$age)-BL.position[,3]->dist.tip
cbind(BL.position,dist.tip)->BL.positions
BL.positions[,5]+BL.positions[,4]->ages
cbind(BL.positions,ages)->BL.positions
as.data.frame(BL.positions)->node.ages
names(node.ages)<-c("parental.node","daughter.node","dist.root","BL","dist.
tip","mrca.age")
## node.ages is a data frame listing as variables the identity of parental and
daughter nodes, the distance from the root and from the present of each node,
the branch lenght and the age of the most recent common ancestor
node.ages[node.ages[,2]species.ages
row.names(species.ages)<-phy$tip
## species ages is node.ages data frame reduced to the tips (species)
>Messaggio originale
>Da: ppi...@uniroma3.it
>Data: 20/06/2011 16.01
>A:
>Ogg: [R-sig-phylo] nodes and taxa depth
>
>Hi all,
>
>I write in order to ask if someone knows a smart
>function returning the depth of nodes and taxa (in the
>same order of ape format) relatively to tree scale.
>I just mean a vector of dim= rbind(n°taxa,n°nodes)x1.
>Does anyone knows a smart strategy?
>node.depth() returns depth in function of number of
>descendant; I need it in function of maximum tree
>length, i.e. the heigth -or depth- of node and taxa
>relatively to tree scale, time in my case.
>Thankyou in advance for your help.
>
>
>Paolo Piras
>
>___
>R-sig-phylo mailing list
>R-sig-phylo@r-project.org
>https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
>
___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] R: Re: R: Re: R: Re: Simulation of Continuous Trait with Active Trend
Hi all,
I tried the function written by Liam, it works pretty well! As an indication
for users, do not use an ultrametric tree, you'd get an error like this:
Error in chol.default(x, pivot = FALSE) :
the leading minor of order 140 is not positive definite
Error in pd.solve(varcov) : x appears to be not positive definite
I applied anc.tree to a multichotomous tree, it still works very well, judging
from fitted values at nodes.
Pas
>Messaggio originale
>Da: liam.rev...@umb.edu
>Data: 22/02/2011 4.58
>A: "pasquale.r...@libero.it"
>Cc: "R Sig Phylo Listserv"
>Ogg: Re: [R-sig-phylo] R: Re: R: Re: Simulation of Continuous Trait with
Active Trend
>
>Hi Pasquale,
>
>To follow up on my earlier message, I have written a preliminary version
>of a simple function to simultaneously estimate the evolutionary
>parameters and ancestral states for Brownian evolution with a trend
>using likelihood. I will paste it at the end of this email. I have
>confirmed that the function returns the same trend parameter as
>fitContinuous(model="trend") and the same ancestral character estimates
>as ace() - when mu (the trend parameter) is forced to zero [note that
>this can only be done by modifying the code]. The parameter "sig2," the
>Brownian rate, is biased by a factor n/(2n-2) for n species relative to
>fitContinuous(model="trend)$Trait1$beta.
>
>Note that this model can generally only be fit for a tree with tips that
>are not contemporaneous - for instance, for a phylogeny containing
>fossil species.
>
>Best, Liam
>
>Function code (also
>http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/anc.trend/v0.1/anc.trend.
R):
>
># written by Liam J. Revell 2011
>anc.trend<-function(phy,x,maxit=2000){
>
> # preliminaries
> # require dependencies
> require(ape)
> # set global
> tol<-1e-8
> # compute C
> D<-dist.nodes(phy)
> ntips<-length(phy$tip.label)
> Cii<-D[ntips+1,]
> C<-D; C[,]<-0
> counts<-vector()
> for(i in 1:nrow(D)) for(j in 1:ncol(D))
> C[i,j]<-(Cii[i]+Cii[j]-D[i,j])/2
> dimnames(C)[[1]][1:length(phy$tip)]<-phy$tip.label
> dimnames(C)[[2]][1:length(phy$tip)]<-phy$tip.label
> C<-C[c(1:ntips,(ntips+2):nrow(C)),c(1:ntips,(ntips+2):ncol(C))]
> # sort x by phy$tip.label
> x<-x[phy$tip.label]
>
> # function returns the negative log-likelihood
> likelihood<-function(theta,x,C){
> a<-theta[1]
> u<-theta[2]
> sig2<-theta[3]
> y<-theta[4:length(theta)]
>
> logLik<-dmnorm(x=c(x,y),mean=(a+diag(C)*u),varcov=sig2*C,log=TRUE)
> return(-logLik)
> }
>
> # get reasonable starting values for the optimizer
> a<-mean(x)
> sig2<-var(x)/max(C)
>
> # perform ML optimization
>
> result<-optim(par=c(a,0,sig2,rep(a,phy$Nnode-1)),likelihood,x=x,C=C,method="
L-BFGS-B",lower=c(-Inf,-Inf,tol,rep(-Inf,phy$Nnode-1)),control=list
(maxit=maxit))
>
> # return the result
> ace<-c(result$par[c(1,4:length(result$par))]);
>names(ace)<-c(as.character(phy$edge[1,1]),rownames(C)[(length(phy$tip.label)
+1):nrow(C)])
>
> return(list(ace=ace,mu=result$par[2],sig2=result$par[3],logL=-result$value,
convergence=result$convergence,message=result$message))
>
>}
>
>--
>Liam J. Revell
>University of Massachusetts Boston
>web: http://www.bio.umb.edu/facstaff/faculty_Revell.html
>(new) email: liam.rev...@umb.edu
>(new) blog: http://phytools.blogspot.com
>
>On 2/17/2011 12:45 PM, Liam J. Revell wrote:
>> Hi Pasquale,
>>
>> I think this may be possible if all tips are not contemporaneous.
>>
>> As Gene points out, the values at the tips of the tree given a trend
>> should be distributed as ~MVN(mu=diag(vcv(tree))*tr+ancestor,
>> V=vcv(tree)*sig2). Thus, we should be able to compute the likelihood of
>> any data pattern at the tips and internal nodes of the tree, given tr,
>> ancestor, and sig2. If we start with only the values at the tips and the
>> phylogeny, we can estimate the states at internal nodes and our model by
>> maximizing the likelihood. This would be our solution.
>>
>> I have written down some code for this using dmnorm() (from the mnormt
>> package) to compute the multivariate normal density; a home-made
>> function to compute vcv(tree) but for all the nodes in the tree; and
>> optim() to maximize the likelihood - but it doesn't seem to work very
>> well yet. I will pass it along if I can figure it out.
>>
>> - Liam
>>
>
___
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
[R-sig-phylo] R: Re: R: Re: R: Re: Simulation of Continuous Trait with Active Trend
Hi Emmanuel, On the theoretical grounds, I agree with you that extinction dynamics must be taken in consideration. This is why, in the first place, I drew interval- phylogenies from the whole tree (reasonably short intervals must have experienced few extinctions anyway). Unfortunately, out of my knowledge no method has ever being developed to study continuous character-associated diversification on non-ultrametric trees. Both yule.cov, and QuaSSE requires ultrametric trees. I'm not sure about Ree's method, must check. Besides maths, one additional difficulty with non-ultrametric trees was pointed out, correctly of course, by Dave: are extinct species nodes or tips? should we qualify all of them the same way or not? We paleontologists must focus on this problem very carefully before any phylo method is applied. Pas >Messaggio originale >Da: emmanuel.para...@ird.fr >Data: 18/02/2011 6.05 >A: "pasquale.r...@libero.it" >Cc: , "R Sig Phylo Listserv", "Liam J. Revell", "Gene Hunt" >Ogg: Re: [R-sig-phylo] R: Re: R: Re: Simulation of Continuous Trait with Active Trend > >Hi Pasquale, > >A note on yule.cov: this method requires ancestral states at the node of >the tree whatever the model used to reconstruct them (not necessarily >BM). Also it assumes no extinction, which could be a realistic >assumption within a given time interval. > >Generally, I think that using this kind of methods (also those in >diversitree) to analyse fossil data is not appropriate because they >analyse phylogenies where extinctions are not observed. You'd better >consider using methods that treat extinctions explicitly, which may >allow you to treat all your time intervals together. > >Emmanuel > >pasquale.r...@libero.it wrote on 18/02/2011 00:22: >> >> Hi all, >> to feed Dave's questions; I have body size data (showing high phylogenetic signal) and a trend verified independently (Cope's). Although all species in the tree are in fact fossil species, I derived from it a number of trees which are time-interval phylogenies, restricted to species living in a particular temporal interval. As such, interval phylogenies are ultrametric. I'm trying ape's yule.cov function to test whether body size did influence diversification, yet yule.cov requires aces which are calculated by assuming brownian motion, which is wrong given I know there is a real trend in the data. This is why I asked if it is possible to calculate aces by assuming BM with a positive (mu>0) trend.Any hint?Pas >> >> >> >> Messaggio originale >> >> Da: dwba...@uchicago.edu >> >> Data: 17/02/2011 18.04 >> >> A: "pasquale.r...@libero.it", "R Sig Phylo Listserv", "Liam J. Revell", "Gene Hunt" >> >> Ogg: Re: [R-sig-phylo] R: Re: Simulation of Continuous Trait with Active Trend >> >> >> >> All- >> Thanks for all the alternative ways to simulate trends! >> >> I'd say estimating ancestral traits always depends on the question you are after. If we are dealing with traits that may be evolving under an active trend, than we know reconstruction could be inaccurate because ancestral states are simply a weighted average of the observed values. Whether that invalidates ones approach depends on the question, how one deals with uncertainty in the ancestral trait estimates and the data itself (i.e. How many fossil taxa do you have? Do you have good a priori evidence for a particular root state? Does the trait have a high or low level of phylogenetic signal, as this influences the uncertainty?). >> >> >> >> I think Finarelli and Flynn have a good point that including a large sample of fossil taxa can help quite a bit in constraining our understanding of active trend patterns, but I think we (as paleontologists) need to give a great deal of consideration to how we want to treat lineages in the fossil record (as internal nodes or as terminal taxa) and how we time-scale our trees. >> >> >> -Dave >> >> >> On Thu, Feb 17, 2011 at 6:19 AM, pasquale.r...@libero.it wrote: >> >> >> >> >> Hi all, >> >> >> >> I have a simple question. Is it possible, or even feasible, to run ancestral >> >> state estimation under Brownian Motion with a trend, as modelled by Emmanuel, >> >> Gene and Liam? I wonder whether this is feasible with real data. >> >> Thanks all >> >> Pas >> >> >> >> >> >> >> >> >> >>> Messaggio originale >> >>> Da: emmanuel.para...@ir
[R-sig-phylo] R: Re: R: Re: Simulation of Continuous Trait with Active Trend
Hi all, to feed Dave's questions; I have body size data (showing high phylogenetic signal) and a trend verified independently (Cope's). Although all species in the tree are in fact fossil species, I derived from it a number of trees which are time-interval phylogenies, restricted to species living in a particular temporal interval. As such, interval phylogenies are ultrametric. I'm trying ape's yule.cov function to test whether body size did influence diversification, yet yule.cov requires aces which are calculated by assuming brownian motion, which is wrong given I know there is a real trend in the data. This is why I asked if it is possible to calculate aces by assuming BM with a positive (mu>0) trend.Any hint?Pas Messaggio originale Da: dwba...@uchicago.edu Data: 17/02/2011 18.04 A: "pasquale.r...@libero.it", "R Sig Phylo Listserv", "Liam J. Revell", "Gene Hunt" Ogg: Re: [R-sig-phylo] R: Re: Simulation of Continuous Trait with Active Trend All- Thanks for all the alternative ways to simulate trends! I'd say estimating ancestral traits always depends on the question you are after. If we are dealing with traits that may be evolving under an active trend, than we know reconstruction could be inaccurate because ancestral states are simply a weighted average of the observed values. Whether that invalidates ones approach depends on the question, how one deals with uncertainty in the ancestral trait estimates and the data itself (i.e. How many fossil taxa do you have? Do you have good a priori evidence for a particular root state? Does the trait have a high or low level of phylogenetic signal, as this influences the uncertainty?). I think Finarelli and Flynn have a good point that including a large sample of fossil taxa can help quite a bit in constraining our understanding of active trend patterns, but I think we (as paleontologists) need to give a great deal of consideration to how we want to treat lineages in the fossil record (as internal nodes or as terminal taxa) and how we time-scale our trees. -Dave On Thu, Feb 17, 2011 at 6:19 AM, pasquale.r...@libero.it wrote: Hi all, I have a simple question. Is it possible, or even feasible, to run ancestral state estimation under Brownian Motion with a trend, as modelled by Emmanuel, Gene and Liam? I wonder whether this is feasible with real data. Thanks all Pas >Messaggio originale >Da: emmanuel.para...@ird.fr >Data: 17/02/2011 6.13 >A: "Hunt, Gene" >Cc: "R Sig Phylo Listserv" >Ogg: Re: [R-sig-phylo] Simulation of Continuous Trait with Active Trend > >Hi all, > >In rTraitCont, the argument model can be a function, eg: > >f <- function(x, l) x + rnorm(1, 1, sqrt(l * 0.1)) > >which is a way to model a trend. Here's an example of what you can do: > >tr <- rbdtree(.1, .05) >tr <- makeNodeLabel(tr) >x <- rTraitCont(tr, model = f, ancestor = TRUE) >bt <- branching.times(tr) >plot(-bt, x[names(bt)]) > >You may add the tip values with: > >n <- Ntip(tr) >points(rep(0, n), x[1:n], pch = 2) > >You can also draw lines linking the ancestors with their descendants in >this "morpho-time-space" since they can be identified in the edge matrix. > >If you don't like this model, just change f. As a side note: > >f <- function(x, l) x + rnorm(1, 0, sqrt(l * 0.1)) >x <- rTraitCont(tr, model = f) > >simulates the same model (ie, BM) than: > >x <- rTraitCont(tr) > >HTH, > >Emmanuel > >Hunt, Gene wrote on 17/02/2011 01:25: >> Hi, >> >> To Liam's extremely helpful and practical reply, I'd add that one can also use theory to generate tip data for a trend. For a given tree, all tip values can be considered a single draw from a multivariate normal distribution, with each dimension representing a tip taxon. For this MVN, the vector of means is equal to the directionality term (mu) times the root-to-tip distance for each terminal taxon, and the variance-covariance is that for BM. So the following two lines also generate tip data for a trend model for phylogeny (phy): >> >> V<- vcv(phy) >> x<- mvrnorm(n=1, diag(V)*mu, V) >> >> Setting n to a higher number will generate more replicates. Liam in his blog talks about why the above code can be slow, but I think it is useful to point out theory. BM can be generated the same way with the multivariate mean equal to the ancestral value repeated as many times as there are tips in the tree. >> >> Best, >> Gene >> >> >> On 2/16/11 1:04 PM, "Liam J. Revell" wrote: >> >> Hi Dave. >> >> I believe BM
[R-sig-phylo] R: Re: Simulation of Continuous Trait with Active Trend
Hi all,
I have a simple question. Is it possible, or even feasible, to run ancestral
state estimation under Brownian Motion with a trend, as modelled by Emmanuel,
Gene and Liam? I wonder whether this is feasible with real data.
Thanks all
Pas
>Messaggio originale
>Da: emmanuel.para...@ird.fr
>Data: 17/02/2011 6.13
>A: "Hunt, Gene"
>Cc: "R Sig Phylo Listserv"
>Ogg: Re: [R-sig-phylo] Simulation of Continuous Trait with Active Trend
>
>Hi all,
>
>In rTraitCont, the argument model can be a function, eg:
>
>f <- function(x, l) x + rnorm(1, 1, sqrt(l * 0.1))
>
>which is a way to model a trend. Here's an example of what you can do:
>
>tr <- rbdtree(.1, .05)
>tr <- makeNodeLabel(tr)
>x <- rTraitCont(tr, model = f, ancestor = TRUE)
>bt <- branching.times(tr)
>plot(-bt, x[names(bt)])
>
>You may add the tip values with:
>
>n <- Ntip(tr)
>points(rep(0, n), x[1:n], pch = 2)
>
>You can also draw lines linking the ancestors with their descendants in
>this "morpho-time-space" since they can be identified in the edge matrix.
>
>If you don't like this model, just change f. As a side note:
>
>f <- function(x, l) x + rnorm(1, 0, sqrt(l * 0.1))
>x <- rTraitCont(tr, model = f)
>
>simulates the same model (ie, BM) than:
>
>x <- rTraitCont(tr)
>
>HTH,
>
>Emmanuel
>
>Hunt, Gene wrote on 17/02/2011 01:25:
>> Hi,
>>
>> To Liam's extremely helpful and practical reply, I'd add that one can also
use theory to generate tip data for a trend. For a given tree, all tip values
can be considered a single draw from a multivariate normal distribution, with
each dimension representing a tip taxon. For this MVN, the vector of means is
equal to the directionality term (mu) times the root-to-tip distance for each
terminal taxon, and the variance-covariance is that for BM. So the following
two lines also generate tip data for a trend model for phylogeny (phy):
>>
>> V<- vcv(phy)
>> x<- mvrnorm(n=1, diag(V)*mu, V)
>>
>> Setting n to a higher number will generate more replicates. Liam in his
blog talks about why the above code can be slow, but I think it is useful to
point out theory. BM can be generated the same way with the multivariate mean
equal to the ancestral value repeated as many times as there are tips in the
tree.
>>
>> Best,
>> Gene
>>
>>
>> On 2/16/11 1:04 PM, "Liam J. Revell" wrote:
>>
>> Hi Dave.
>>
>> I believe BM with a trend should just have a non-zero mean for the
>> random normal changes along individual branches of the tree. It was
>> easy to add this to a function that I already have for fast Brownian
>> motion simulation. I will post this update to my website:
>> http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/ but it is so simple
>> that I have also pasted the code below.
>>
>> - Liam
>>
>> # Simulates BM evolution more quickly.
>> # A trend can be simulated by mu!=0.
>> # mu=0 is standard BM; mu<0 downward trend; mu>0 upward trend.
>> # Written by Liam J. Revell 2011
>> fastBM.trend<-function(tree,a=0,mu=0,sig2=1,internal=FALSE){
>>
>> n<-length(tree$tip)
>>
>> # first simulate changes along each branch
>> x<-rnorm(n=length(tree$edge.length),mean=mu,sd=sqrt(sig2*tree$edge.
length))
>>
>> # now add them up
>> y<-matrix(0,nrow(tree$edge),ncol(tree$edge))
>> for(i in 1:length(x)){
>> if(tree$edge[i,1]==(n+1))
>> y[i,1]<-a
>> else
>> y[i,1]<-y[match(tree$edge[i,1],tree$edge[,2]),2]
>>
>> y[i,2]<-y[i,1]+x[i]
>> }
>>
>> rm(x); x<-c(y[1,1],y[,2])
>> names(x)<-c(n+1,tree$edge[,2])
>> x<-x[as.character(1:(n+tree$Nnode))]
>> names(x)[1:n]<-tree$tip.label
>>
>> # return simulated data
>> if(internal==TRUE)
>> return(x) # include internal nodes
>> else
>> return(x[1:length(tree$tip.label)]) # tip nodes only
>>
>> }
>>
>>
>> --
>> Liam J. Revell
>> University of Massachusetts Boston
>> web: http://www.bio.umb.edu/facstaff/faculty_Revell.html
>> (new) email: liam.rev...@umb.edu
>> (new) blog: http://phytools.blogspot.com
>>
>> On 2/16/2011 12:28 PM, David Bapst wrote:
>>> All-
>>> For a class presentation I was working on, I wanted other students to use
>>> model fitting with fitContinuous() on simulated traits, under different
>>> models of trait evolution. I was surprised to discover there seemed to be
no
>>> easy way to simulate active directional trend mechanisms, unlike other
trait
>>> evolution models.
>>>
>>> One idea I had was to kludge it.
>>> tree<-rtree(20)
>>> trait1<-,rTraitCont(tree,sigma=1)+2*diag(vcv.phylo(tree))
>>>
>>> An alternative would be an OU with the root far from the theta value.
Takes
>>> some playing with the parameters to get something that looks like a trend
in
>>> a traitgram.
>>> trait2<-rTraitCont(tree,model="OU",sigma=1,alpha=0.5,root.value=-10)
>>>
>>> Both of these were strongly preferred as 'trend' in fitC
