man/rTraitCont.Rd

   1 \name{rTraitCont}
   2 \alias{rTraitCont}
   3 \title{Continuous Character Simulation}
   4 \usage{
   5 rTraitCont(phy, model = "BM", sigma = 0.1, alpha = 1, theta = 0,
   6            ancestor = FALSE, root.value = 0, ...)
   7 }
   8 \arguments{
   9   \item{phy}{an object of class \code{"phylo"}.}
  10   \item{model}{a character (either \code{"BM"} or \code{"OU"}) or a
  11     function specifying the model (see details).}
  12   \item{sigma}{a numeric vector giving the standard-deviation of the
  13     random component for each branch (can be a single value).}
  14   \item{alpha}{if \code{model = "OU"}, a numeric vector giving the
  15     strength of the selective constraint for each branch (can be a
  16     single value).}
  17   \item{theta}{if \code{model = "OU"}, a numeric vector giving the
  18     optimum for each branch (can be a single value).}
  19   \item{ancestor}{a logical value specifying whether to return the
  20     values at the nodes as well (by default, only the values at the tips
  21     are returned).}
  22   \item{root.value}{a numeric giving the value at the root.}
  23   \item{\dots}{further arguments passed to \code{model} if it is a
  24     function.}
  25 }
  26 \description{
  27   This function simulates the evolution of a continuous character along a
  28   phylogeny. The calculation is done recursively from the root. See
  29   Paradis (2012, pp. 232 and 324) for an introduction.
  30 }
  31 \details{
  32   There are three possibilities to specify \code{model}:
  33
  34 \itemize{
  35   \item{\code{"BM"}:}{a Browian motion model is used. If the arguments
  36   \code{sigma} has more than one value, its length must be equal to the
  37   the branches of the tree. This allows to specify a model with variable
  38   rates of evolution. You must be careful that branch numbering is done
  39   with the tree in ``pruningwise'' order: to see the order of the branches
  40   you can use: \code{tr <- reorder(tr, "p"); plor(tr); edgelabels()}.
  41   The arguments \code{alpha} and \code{theta} are ignored.}
  42
  43   \item{\code{"OU"}:}{an Ornstein-Uhlenbeck model is used. The above
  44   indexing rule is used for the three parameters \code{sigma},
  45   \code{alpha}, and \code{theta}. This may be interesting for the last
  46   one to model varying phenotypic optima. The exact updating formula
  47   from Gillespie (1996) are used which are reduced to BM formula if
  48   \code{alpha = 0}.}
  49
  50   \item{A function:}{it must be of the form \code{foo(x, l)} where
  51   \code{x} is the trait of the ancestor and \code{l} is the branch
  52   length. It must return the value of the descendant. The arguments
  53   \code{sigma}, \code{alpha}, and \code{theta} are ignored.}
  54 }}
  55 \value{
  56   A numeric vector with names taken from the tip labels of
  57   \code{phy}. If \code{ancestor = TRUE}, the node labels are used if
  58   present, otherwise, ``Node1'', ``Node2'', etc.
  59 }
  60 \references{
  61   Gillespie, D. T. (1996) Exact numerical simulation of the
  62   Ornstein-Uhlenbeck process and its integral. \emph{Physical Review E},
  63   \bold{54}, 2084--2091.
  64
  65   Paradis, E. (2012) \emph{Analysis of Phylogenetics and Evolution with
  66     R (Second Edition).} New York: Springer.
  67 }
  68 \author{Emmanuel Paradis}
  69 \seealso{
  70   \code{\link{rTraitDisc}}, \code{\link{rTraitMult}}, \code{\link{ace}}
  71 }
  72 \examples{
  73 data(bird.orders)
  74 rTraitCont(bird.orders) # BM with sigma = 0.1
  75 ### OU model with two optima:
  76 tr <- reorder(bird.orders, "p")
  77 plot(tr)
  78 edgelabels()
  79 theta <- rep(0, Nedge(tr))
  80 theta[c(1:4, 15:16, 23:24)] <- 2
  81 ## sensitive to 'alpha' and 'sigma':
  82 rTraitCont(tr, "OU", theta = theta, alpha=.1, sigma=.01)
  83 ### an imaginary model with stasis 0.5 time unit after a node, then
  84 ### BM evolution with sigma = 0.1:
  85 foo <- function(x, l) {
  86     if (l <= 0.5) return(x)
  87     x + (l - 0.5)*rnorm(1, 0, 0.1)
  88 }
  89 tr <- rcoal(20, br = runif)
  90 rTraitCont(tr, foo, ancestor = TRUE)
  91 ### a cumulative Poisson process:
  92 bar <- function(x, l) x + rpois(1, l)
  93 (x <- rTraitCont(tr, bar, ancestor = TRUE))
  94 plot(tr, show.tip.label = FALSE)
  95 Y <- x[1:20]
  96 A <- x[-(1:20)]
  97 nodelabels(A)
  98 tiplabels(Y)
  99 }
 100 \keyword{datagen}