pegas (version 0.8-2)

theta.tree: Population Parameter THETA Using Genealogy

Description

This function estimates the population parameter $\theta$ from a genealogy (coded a as phylogenetic tree) under the coalescent.

Usage

theta.tree(phy, theta, fixed = FALSE, analytical = TRUE, log = TRUE)

Arguments

phy

an object of class "phylo".

theta

a numeric vector.

fixed

a logical specifying whether to estimate theta (the default), or to return the likelihoods for all values in theta.

analytical

a logical specifying whether to use analytical formulae to estimate theta and its standard-error. If FALSE, a numerical optimisation of the likelihood is performed (this option is ignored if fixed = TRUE)

log

a logical specifying whether to return the likelihoods on a log scale (the default); ignored if fixed = FALSE.

Value

• If fixed = FALSE, a list with two elements:
• theta
  the maximum likelihood estimate of $\theta$;
• logLik
  the log-likelihood at its maximum.
• If fixed = TRUE, a numeric vector with the [log-]likelihood values.

Details

The tree phy is considered as a genealogy, and therefore should be ultrametric. 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.

References

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.

See Also

theta.h, theta.s, theta.k

Examples

tr <- rcoal(50) # assumes theta = 1
theta.tree(tr, 10)
theta.tree(tr, 10, analytical = FALSE) # uses nlminb()
## profile log-likelihood:
THETA <- seq(0.5, 1.5, 0.01)
logLikelihood <- theta.tree(tr, THETA, fixed = TRUE)
plot(THETA, logLikelihood, type = "l")