If type = "continuous"
, the default model is Brownian motion
where characters evolve randomly following a random walk. This model
can be fitted by residual maximum likelihood (the default), maximum
likelihood (Felsenstein 1973, Schluter et al. 1997), least squares
(method = "pic"
, Felsenstein 1985), or generalized least
squares (method = "GLS"
, Martins and Hansen 1997, Cunningham et
al. 1998). In the last case, the specification of phy
and
model
are actually ignored: it is instead given through a
correlation structure with the option corStruct
. In the setting method = "ML"
and model = "BM"
(this used
to be the default until ape 3.0-7) the maximum likelihood
estimation is done simultaneously on the ancestral values and the
variance of the Brownian motion process; these estimates are then used
to compute the confidence intervals in the standard way. The REML
method first estimates the ancestral value at the root (aka, the
phylogenetic mean), then the variance of the Brownian motion process
is estimated by optimizing the residual log-likelihood. The ancestral
values are finally inferred from the likelihood function giving these
two parameters. If method = "pic"
or "GLS"
, the
confidence intervals are computed using the expected variances under
the model, so they depend only on the tree.
It could be shown that, with a continous character, REML results in
unbiased estimates of the variance of the Brownian motion process
while ML gives a downward bias. Therefore the former is recommanded.
For discrete characters (type = "discrete"
), only maximum
likelihood estimation is available (Pagel 1994) (see MPR
for an alternative method). The model is specified through a numeric
matrix with integer values taken as indices of the parameters. The
numbers of rows and of columns of this matrix must be equal, and are
taken to give the number of states of the character. For instance,
matrix(c(0, 1, 1, 0), 2)
will represent a model with two
character states and equal rates of transition, matrix(c(0, 1,
2, 0), 2)
a model with unequal rates, matrix(c(0, 1, 1, 1, 0,
1, 1, 1, 0), 3)
a model with three states and equal rates of
transition (the diagonal is always ignored). There are short-cuts to
specify these models: "ER"
is an equal-rates model (e.g., the
first and third examples above), "ARD"
is an
all-rates-different model (the second example), and "SYM"
is a
symmetrical model (e.g., matrix(c(0, 1, 2, 1, 0, 3, 2, 3, 0),
3)
). If a short-cut is used, the number of states is determined from
the data.
With discrete characters it is necessary to compute the exponential of
the rate matrix. The only possibility until ape 3.0-7) was the
function matexpo
in ape. If use.expm = TRUE
and use.eigen = FALSE
, the function expm
,
in the package of the same name, is used. matexpo
is faster but
quite inaccurate for large and/or asymmetric matrices. In case of
doubt, use the latter. Since ape 3.0-10, it is possible to use
an eigen decomposition avoiding the need to compute the matrix
exponential; see details in Lebl (2013, sect. 3.8.3). This is much
faster and is now the default.