dbd {ape}

Probability Density Under Birth–Death Models

Description

These functions compute the probability density under some birth–death models, that is the probability of obtaining x species after a time t giving how speciation and extinction probabilities vary through time (these may be constant, or even equal to zero for extinction).

Usage

dyule(x, lambda = 0.1, t = 1, log = FALSE)
dbd(x, lambda, mu, t, conditional = FALSE, log = FALSE)
dbdTime(x, birth, death, t, conditional = FALSE,
        BIRTH = NULL, DEATH = NULL, fast = FALSE)

Arguments

x

a numeric vector of species numbers (see Details).

lambda

a numerical value giving the probability of speciation; can be a vector with several values for dyule.

mu

id. for extinction.

t

id. for the time(s).

log

a logical value specifying whether the probabilities should be returned log-transformed; the default is FALSE.

conditional

a logical specifying whether the probabilities should be computed conditional under the assumption of no extinction after time t.

birth, death

a (vectorized) function specifying how the speciation or extinction probability changes through time (see yule.time and below).

BIRTH, DEATH

a (vectorized) function giving the primitive of birth or death.

fast

a logical value specifying whether to use faster integration (see bd.time).

Details

These three functions compute the probabilities to observe x species starting from a single one after time t (assumed to be continuous). The first function is a short-cut for the second one with mu = 0 and with default values for the two other arguments. dbdTime is for time-varying lambda and mu specified as R functions.

dyule is vectorized simultaneously on its three arguments x, lambda, and t, according to R's rules of recycling arguments. dbd is vectorized simultaneously x and t (to make likelihood calculations easy), and dbdTime is vectorized only on x; the other arguments are eventually shortened with a warning if necessary.

The returned value is, logically, zero for values of x out of range, i.e., negative or zero for dyule or if conditional = TRUE. However, it is not checked if the values of x are positive non-integers and the probabilities are computed and returned.

The details on the form of the arguments birth, death, BIRTH, DEATH, and fast can be found in the links below.

Value

a numeric vector.

Note

If you use these functions to calculate a likelihood function, it is strongly recommended to compute the log-likelihood with, for instance in the case of a Yule process, sum(dyule( , log = TRUE)) (see examples).

Author(s)

Emmanuel Paradis

References

Kendall, D. G. (1948) On the generalized “birth-and-death” process. Annals of Mathematical Statistics, 19, 1–15.

See Also

bd.time, yule.time

Examples

x <- 0:10
plot(x, dyule(x), type = "h", main = "Density of the Yule process")
text(7, 0.85, expression(list(lambda == 0.1, t == 1)))

y <- dbd(x, 0.1, 0.05, 10)
z <- dbd(x, 0.1, 0.05, 10, conditional = TRUE)
d <- rbind(y, z)
colnames(d) <- x
barplot(d, beside = TRUE, ylab = "Density", xlab = "Number of species",
        legend = c("unconditional", "conditional on\nno extinction"),
        args.legend = list(bty = "n"))
title("Density of the birth-death process")
text(17, 0.4, expression(list(lambda == 0.1, mu == 0.05, t == 10)))

## Not run: 
### generate 1000 values from a Yule process with lambda = 0.05
x <- replicate(1e3, Ntip(rlineage(0.05, 0)))

### the correct way to calculate the log-likelihood...:
sum(dyule(x, 0.05, 50, log = TRUE))

### ... and the wrong way:
log(prod(dyule(x, 0.05, 50)))

### a third, less preferred, way:
sum(log(dyule(x, 0.05, 50)))

## End(Not run)