man/ace.Rd

   1 \name{ace}
   2 \alias{ace}
   3 \alias{print.ace}
   4 \alias{logLik.ace}
   5 \alias{deviance.ace}
   6 \alias{AIC.ace}
   7 \alias{anova.ace}
   8 \title{Ancestral Character Estimation}
   9 \usage{
  10 ace(x, phy, type = "continuous", method = "ML", CI = TRUE,
  11     model = if (type == "continuous") "BM" else "ER",
  12     scaled = TRUE, kappa = 1, corStruct = NULL, ip = 0.1)
  13 \method{print}{ace}(x, digits = 4, ...)
  14 \method{logLik}{ace}(object, ...)
  15 \method{deviance}{ace}(object, ...)
  16 \method{AIC}{ace}(object, ..., k = 2)
  17 \method{anova}{ace}(object, ...)
  18 }
  19 \arguments{
  20   \item{x}{a vector or a factor; an object of class \code{"ace"} in the
  21     case of \code{print}.}
  22   \item{phy}{an object of class \code{"phylo"}.}
  23   \item{type}{the variable type; either \code{"continuous"} or
  24     \code{"discrete"} (or an abbreviation of these).}
  25   \item{method}{a character specifying the method used for
  26     estimation. Three choices are possible: \code{"ML"}, \code{"pic"},
  27     or \code{"GLS"}.}
  28   \item{CI}{a logical specifying whether to return the 95\% confidence
  29     intervals of the ancestral state estimates (for continuous
  30     characters) or the likelihood of the different states (for discrete
  31     ones).}
  32   \item{model}{a character specifying the model (ignored if \code{method
  33       = "GLS"}), or a numeric matrix if \code{type = "discrete"} (see
  34     details).}
  35   \item{scaled}{a logical specifying whether to scale the contrast
  36     estimate (used only if \code{method = "pic"}).}
  37   \item{kappa}{a positive value giving the exponent transformation of
  38     the branch lengths (see details).}
  39   \item{corStruct}{if \code{method = "GLS"}, specifies the correlation
  40     structure to be used (this also gives the assumed model).}
  41   \item{ip}{the initial value(s) used for the ML estimation procedure
  42     when \code{type == "discrete"} (possibly recycled).}
  43   \item{digits}{the number of digits to be printed.}
  44   \item{object}{an object of class \code{"ace"}.}
  45   \item{k}{a numeric value giving the penalty per estimated parameter;
  46     the default is \code{k = 2} which is the classical Akaike
  47     information criterion.}
  48   \item{\dots}{further arguments passed to or from other methods.}
  49 }
  50 \description{
  51   This function estimates ancestral character states, and the associated
  52   uncertainty, for continuous and discrete characters.
  53
  54   \code{logLik}, \code{deviance}, and \code{AIC} are generic functions
  55   used to extract the log-likelihood, the deviance (-2*logLik), or the
  56   Akaike information criterion of a tree. If no such values are
  57   available, \code{NULL} is returned.
  58
  59   \code{anova} is another generic function that is used to compare
  60   nested models: the significance of the additional parameter(s) is
  61   tested with likelihood ratio tests. You must ensure that the models
  62   are effectively nested (if they are not, the results will be
  63   meaningless). It is better to list the models from the smallest to the
  64   largest.
  65 }
  66 \details{
  67   If \code{type = "continuous"}, the default model is Brownian motion
  68   where characters evolve randomly following a random walk. This model
  69   can be fitted by maximum likelihood (the default, Schluter et
  70   al. 1997), least squares (\code{method = "pic"}, Felsenstein 1985), or
  71   generalized least squares (\code{method = "GLS"}, Martins and Hansen
  72   1997, Cunningham et al. 1998). In the latter case, the specification
  73   of \code{phy} and \code{model} are actually ignored: it is instead
  74   given through a correlation structure with the option \code{corStruct}.
  75
  76   For discrete characters (\code{type = "discrete"}), only maximum
  77   likelihood estimation is available (Pagel 1994). The model is
  78   specified through a numeric matrix with integer values taken as
  79   indices of the parameters. The numbers of rows and of columns of this
  80   matrix must be equal, and are taken to give the number of states of
  81   the character. For instance, \code{matrix(c(0, 1, 1, 0), 2)} will
  82   represent a model with two character states and equal rates of
  83   transition, \code{matrix(c(0, 1, 2, 0), 2)} a model with unequal
  84   rates, \code{matrix(c(0, 1, 1, 1, 0, 1, 1, 1, 0), 3)} a model with
  85   three states and equal rates of transition (the diagonal is always
  86   ignored). There are short-cuts to specify these models: \code{"ER"} is
  87   an equal-rates model (e.g., the first and third examples above),
  88   \code{"ARD"} is an all-rates-different model (the second example), and
  89   \code{"SYM"} is a symmetrical model (e.g., \code{matrix(c(0, 1, 2, 1,
  90     0, 3, 2, 3, 0), 3)}). If a short-cut is used, the number of states
  91   is determined from the data.
  92 }
  93 \value{
  94   a list with the following elements:
  95
  96   \item{ace}{if \code{type = "continuous"}, the estimates of the
  97     ancestral character values.}
  98   \item{CI95}{if \code{type = "continuous"}, the estimated 95\%
  99     confidence intervals.}
 100   \item{sigma2}{if \code{type = "continuous"}, \code{model = "BM"}, and
 101     \code{method = "ML"}, the maximum likelihood estimate of the
 102     Brownian parameter.}
 103   \item{rates}{if \code{type = "discrete"}, the maximum likelihood
 104     estimates of the transition rates.}
 105   \item{se}{if \code{type = "discrete"}, the standard-errors of
 106     estimated rates.}
 107   \item{index.matrix}{if \code{type = "discrete"}, gives the indices of
 108     the \code{rates} in the rate matrix.}
 109   \item{loglik}{if \code{method = "ML"}, the maximum log-likelihood.}
 110   \item{lik.anc}{if \code{type = "discrete"}, the scaled likelihoods of
 111     each ancestral state.}
 112   \item{call}{the function call.}
 113 }
 114 \references{
 115   Cunningham, C. W., Omland, K. E. and Oakley, T. H. (1998)
 116   Reconstructing ancestral character states: a critical
 117   reappraisal. \emph{Trends in Ecology & Evolution}, \bold{13},
 118   361--366.
 119
 120   Felsenstein, J. (1985) Phylogenies and the comparative
 121   method. \emph{American Naturalist}, \bold{125}, 1--15.
 122
 123   Martins, E. P. and Hansen, T. F. (1997) Phylogenies and the
 124   comparative method: a general approach to incorporating phylogenetic
 125   information into the analysis of interspecific data. \emph{American
 126     Naturalist}, \bold{149}, 646--667.
 127
 128   Pagel, M. (1994) Detecting correlated evolution on phylogenies: a
 129   general method for the comparative analysis of discrete
 130   characters. \emph{Proceedings of the Royal Society of London. Series
 131     B. Biological Sciences}, \bold{255}, 37--45.
 132
 133   Schluter, D., Price, T., Mooers, A. O. and Ludwig, D. (1997)
 134   Likelihood of ancestor states in adaptive radiation. \emph{Evolution},
 135   \bold{51}, 1699--1711.
 136 }
 137 \author{Emmanuel Paradis, Ben Bolker \email{bolker@zoo.ufl.edu}}
 138 \seealso{
 139   \code{\link{corBrownian}}, \code{\link{corGrafen}},
 140   \code{\link{corMartins}}, \code{\link{compar.ou}},
 141   \code{\link[stats]{anova}}
 142 }
 143 \examples{
 144 ### Just some random data...
 145 data(bird.orders)
 146 x <- rnorm(23)
 147 ### Compare the three methods for continuous characters:
 148 ace(x, bird.orders)
 149 ace(x, bird.orders, method = "pic")
 150 ace(x, bird.orders, method = "GLS",
 151     corStruct = corBrownian(1, bird.orders))
 152 ### For discrete characters:
 153 x <- factor(c(rep(0, 5), rep(1, 18)))
 154 ans <- ace(x, bird.orders, type = "d")
 155 #### Showing the likelihoods on each node:
 156 plot(bird.orders, type = "c", FALSE, label.offset = 1)
 157 co <- c("blue", "yellow")
 158 tiplabels(pch = 22, bg = co[as.numeric(x)], cex = 2, adj = 1)
 159 nodelabels(thermo = ans$lik.anc, piecol = co, cex = 0.75)
 160 ### An example of the use of the argument `ip':
 161 tr <- character(4)
 162 tr[1] <- "((((t10:5.03,t2:5.03):2.74,(t9:4.17,"
 163 tr[2] <- "t5:4.17):3.60):2.80,(t3:4.05,t7:"
 164 tr[3] <- "4.05):6.53):2.32,((t6:4.38,t1:4.38):"
 165 tr[4] <- "2.18,(t8:2.17,t4:2.17):4.39):6.33);"
 166 tr <- read.tree(text = paste(tr, collapse = ""))
 167 y <- c(rep(1, 6), rep(2, 4))
 168 ### The default `ip = 0.1' makes ace fails:
 169 ace(y, tr, type = "d")
 170 ace(y, tr, type = "d", ip = 0.01)
 171 ### Surprisingly, using an initial value farther to the
 172 ### MLE than the default one works:
 173 ace(y, tr, type = "d", ip = 0.3)
 174 }
 175 \keyword{models}