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