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

--
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:22 PM, pasquale.r...@libero.it wrote:


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"<pasquale.r...@libero.it>, "R Sig Phylo
    Listserv"<r-sig-phylo@r-project.org>, "Liam J.
    Revell"<liam.rev...@umb.edu>, "Gene Hunt"<hu...@si.edu>
    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
    <mailto:pasquale.r...@libero.it> <pasquale.r...@libero.it
    <mailto: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 <mailto:emmanuel.para...@ird.fr>
         >Data: 17/02/2011 6.13
         >A: "Hunt, Gene"<hu...@si.edu <mailto:hu...@si.edu>>
         >Cc: "R Sig Phylo Listserv"<r-sig-phylo@r-project.org
        <mailto:r-sig-phylo@r-project.org>>
         >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" <liam.rev...@umb.edu
        <mailto:liam.rev...@umb.edu>> 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/
        <http://anolis.oeb.harvard.edu/%7Eliam/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 <mailto: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
        fitContinuous()
        compared
         >>> to BM and OU1, but I was wondering if there were any
        alternatives I was
         >>> missing.
         >>> -Dave
         >>>
         >>
         >> _______________________________________________
         >> R-sig-phylo mailing list
         >> R-sig-phylo@r-project.org <mailto: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
         >> http://paleobiology.si.edu/staff/individuals/hunt.cfm
         >>
         >> _______________________________________________
         >> R-sig-phylo mailing list
         >> R-sig-phylo@r-project.org <mailto:R-sig-phylo@r-project.org>
         >> https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
         >>
         >
         >--
         >Emmanuel Paradis
         >IRD, Jakarta, Indonesia
         >http://ape.mpl.ird.fr/
         >
         >_______________________________________________
         >R-sig-phylo mailing list
         >R-sig-phylo@r-project.org <mailto: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 <mailto:R-sig-phylo@r-project.org>
        https://stat.ethz.ch/mailman/listinfo/r-sig-phylo




    --
    David Bapst
    Dept of Geophysical Sciences
    University of Chicago
    5734 S. Ellis
    Chicago, IL 60637
    http://home.uchicago.edu/~dwbapst/
    <http://home.uchicago.edu/%7Edwbapst/>



_______________________________________________
R-sig-phylo mailing list
R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo

Reply via email to