pegas (version 0.10)

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

# NOT RUN {
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")
# }