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, linear = TRUE, ...)
   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{linear}{a logical indicating which parameterisation of the OU
  24     model to use (see details).}
  25   \item{\dots}{further arguments passed to \code{model} if it is a
  26     function.}
  27 }
  28 \description{
  29   This function simulates the evolution of a continuous character along a
  30   phylogeny. The calculation is done recursively from the root. See
  31   Paradis (2006, p. 151) for a brief introduction.
  32 }
  33 \details{
  34   There are three possibilities to specify \code{model}:
  35
  36 \itemize{
  37   \item{\code{"BM"}:}{a Browian motion model is used. If the arguments
  38   \code{sigma} has more than one value, its length must be equal to the
  39   the branches of the tree. This allows to specify a model with variable
  40   rates of evolution. You must be careful that branch numbering is done
  41   with the tree in ``pruningwise'' order: to see the order of the branches
  42   you can use: \code{tr <- reorder(tr, "p"); plor(tr); edgelabels()}.
  43   The arguments \code{alpha} and \code{theta} are ignored.}
  44
  45   \item{\code{"OU"}:}{an Ornstein-Uhlenbeck model is used. The above
  46   indexing rule is used for the three parameters \code{sigma},
  47   \code{alpha}, and \code{theta}. This may be more interesting for the
  48   last one to model varying phenotypic optima.
  49
  50   By default the following formula is used:
  51
  52   \deqn{x_{t''} = x_{t'} - \alpha l (x_{t'} - \theta) + \sigma
  53     l \epsilon}{x(t'') = x(t') - alpha l (x(t') - theta) + sigma
  54     l epsilon}
  55
  56   where \eqn{l (= t'' - t')} is the branch length, and \eqn{\epsilon
  57   \sim N(0, 1)}{\epsilon ~ N(0, 1)}. If \eqn{\alpha > 1}{alpha > 1},
  58   this may lead to chaotic oscillations. Thus an alternative
  59   parameterisation is used if \code{linear = FALSE}:
  60
  61   \deqn{x_{t''} = x_{t'} - (1 - exp(-\alpha l)) * (x_{t'} - \theta) +
  62     \sigma l \epsilon}{x(t'') = x(t') - (1 - exp(-alpha l)) * (x(t') -
  63     theta) + sigma l epsilon}}
  64
  65   \item{A function:}{it must be of the form \code{foo(x, l)} where
  66   \code{x} is the trait of the ancestor and \code{l} is the branch
  67   length. It must return the value of the descendant. The arguments
  68   \code{sigma}, \code{alpha}, and \code{theta} are ignored.}
  69 }}
  70 \value{
  71   A numeric vector with names taken from the tip labels of
  72   \code{phy}. If \code{ancestor = TRUE}, the node labels are used if
  73   present, otherwise, ``Node1'', ``Node2'', etc.
  74 }
  75 \references{
  76   Paradis, E. (2006) \emph{Analyses of Phylogenetics and Evolution with
  77     R.} New York: Springer.
  78 }
  79 \author{Emmanuel Paradis}
  80 \seealso{
  81   \code{\link{rTraitDisc}}, \code{\link{rTraitMult}}, \code{\link{ace}}
  82 }
  83 \examples{
  84 data(bird.orders)
  85 rTraitCont(bird.orders) # BM with sigma = 0.1
  86 ### OU model with two optima:
  87 tr <- reorder(bird.orders, "p")
  88 plot(tr)
  89 edgelabels()
  90 theta <- rep(0, Nedge(tr))
  91 theta[c(1:4, 15:16, 23:24)] <- 2
  92 ## sensitive to 'alpha' and 'sigma':
  93 rTraitCont(tr, "OU", theta = theta, alpha=.1, sigma=.01)
  94 ### an imaginary model with stasis 0.5 time unit after a node, then
  95 ### BM evolution with sigma = 0.1:
  96 foo <- function(x, l) {
  97     if (l <= 0.5) return(x)
  98     x + (l - 0.5)*rnorm(1, 0, 0.1)
  99 }
 100 tr <- rcoal(20, br = runif)
 101 rTraitCont(tr, foo, ancestor = TRUE)
 102 ### a cumulative Poisson process:
 103 bar <- function(x, l) x + rpois(1, l)
 104 (x <- rTraitCont(tr, bar, ancestor = TRUE))
 105 plot(tr, show.tip.label = FALSE)
 106 Y <- x[1:20]
 107 A <- x[-(1:20)]
 108 nodelabels(A)
 109 tiplabels(Y)
 110 }
 111 \keyword{datagen}