theta.tree | R Documentation |
These functions estimate the population parameter \Theta
from a genealogy (coded a as phylogenetic tree) under the coalescent.
theta.tree(phy, theta, fixed = FALSE, analytical = TRUE, log = TRUE)
theta.tree.hetero(phy, theta, fixed = FALSE, log = TRUE)
phy |
an object of class |
theta |
a numeric vector. |
fixed |
a logical specifying whether to estimate |
analytical |
a logical specifying whether to use analytical
formulae to estimate |
.
log |
a logical specifying whether to return the likelihoods on a
log scale (the default); ignored if |
With theta.tree
, the tree phy
is considered as a
genealogy with contemporaneous samples, and therefore should be
ultrametric. With theta.tree.hetero
, the samples may be
heterochronous so phy
can be non-ultrametric. If phy
is
ultrametric, both functions return the same results.
By default, \theta
is estimated by maximum likelihood and
the value given in theta
is used as starting value for the
minimisation function (if several values are given as a vector the
first one is used). If fixed = TRUE
, then the [log-]likelihood
values are returned corresponding to each value in theta
.
The present implementation does a numerical optimisation of the
log-likelihood function (with nlminb
) with the
first partial derivative as gradient. It is possible to solve the
latter and have a direct analytical MLE of \theta
(and
its standard-error), but this does not seem to be faster.
If fixed = FALSE
, a list with two elements:
theta |
the maximum likelihood estimate of |
logLik |
the log-likelihood at its maximum. |
If fixed = TRUE
, a numeric vector with the [log-]likelihood
values.
Emmanuel Paradis
Kingman, J. F. C. (1982) The coalescent. Stochastic Processes and their Applications, 13, 235–248.
Kingman, J. F. C. (1982) On the genealogy of large populations. Journal of Applied Probability, 19A, 27–43.
Wakeley, J. (2009) Coalescent Theory: An Introduction. Greenwood Village, CO: Roberts and Company Publishers.
theta.h
, theta.s
, theta.k
tr <- rcoal(50)
(o <- theta.tree(tr))
theta.tree(tr, 10, analytical = FALSE) # uses nlminb()
## profile log-likelihood:
THETA <- seq(0.5, 2, 0.01)
logLikelihood <- theta.tree(tr, THETA, fixed = TRUE)
plot(THETA, logLikelihood, type = "l")
xx <- seq(o$theta - 1.96 * o$se, o$theta + 1.96 * o$se, 0.01)
yy <- theta.tree(tr, xx, fixed = TRUE)
polygon(c(xx, rev(xx)), c(yy, rep(0, length(xx))),
border = NA, col = "lightblue")
segments(o$theta, 0, o$theta, o$logLik, col = "blue")
abline(v = 1, lty = 3)
legend("topright", legend = expression("log-likelihood",
"True " * theta, hat(theta) * " (MLE)", "95%\ conf. interv."),
lty = c(1, 3, 1, 1), lwd = c(1, 1, 1, 15),
col = c("black", "black", "blue", "lightblue"))
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.