8 \title{Ancestral Character Estimation}
10 This function estimates ancestral character states, and the associated
11 uncertainty, for continuous and discrete characters.
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.
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
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, ...)
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
48 \item{model}{a character specifying the model (ignored if \code{method
49 = "GLS"}), or a numeric matrix if \code{type = "discrete"} (see
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.}
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
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
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.
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.
113 an object of class \code{"ace"} with the following elements:
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
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
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.}
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},
139 Felsenstein, J. (1985) Phylogenies and the comparative
140 method. \emph{American Naturalist}, \bold{125}, 1--15.
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.
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.
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.
156 \author{Emmanuel Paradis, Ben Bolker}
158 \code{\link{corBrownian}}, \code{\link{corGrafen}},
159 \code{\link{corMartins}}, \code{\link{compar.ou}},
160 \code{\link[stats]{anova}}
163 ### Just some random data...
166 ### Compare the three methods for continuous characters:
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':
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)