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