1 ## birthdeath.R (2007-10-30)
3 ## Estimation of Speciation and Extinction Rates
4 ## with Birth-Death Models
6 ## birthdeath: standard model
7 ## bd.ext: extended version
9 ## Copyright 2002-2007 Emmanuel Paradis
11 ## This file is part of the R-package `ape'.
12 ## See the file ../COPYING for licensing issues.
14 birthdeath <- function(phy)
16 if (class(phy) != "phylo") stop('object "phy" is not of class "phylo"')
17 N <- length(phy$tip.label)
18 x <- c(NA, branching.times(phy))
19 dev <- function(a, r) {
20 -2 * (lfactorial(N - 1)
24 - 2 * sum(log(exp(r * x[2:N]) - a)))
26 out <- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE)
27 if (out$estimate[1] < 0) {
28 out <- nlm(function(p) dev(0, p), 0.2, hessian = TRUE)
29 para <- c(0, out$estimate)
30 se <- c(0, sqrt(diag(solve(out$hessian))))
34 se <- sqrt(diag(solve(out$hessian)))
37 ## compute the 95 % profile likelihood CIs
38 ## (not very clean... but seems to work -- EP 2003-03-29)
39 CI <- matrix(NA, 2, 2)
40 foo <- function(p) dev(p, para[2]) - 3.84 - Dev
44 while (foo(lo) < 0) lo <- lo - inc
45 while (foo(up) < 0) up <- up + inc
46 CI[1, 1] <- uniroot(foo, lower = lo, upper = para[1])$root
47 if (CI[1, 1] < 0) CI[1, 1] <- 0
48 CI[1, 2] <- uniroot(foo, lower = para[1], upper = up)$root
49 foo <- function(p) dev(para[1], p) - 3.84 - Dev
52 while (foo(lo) < 0) lo <- lo - inc
53 while (foo(up) < 0) up <- up + inc
54 CI[2, 1] <- uniroot(foo, lower = lo, upper = para[2])$root
55 CI[2, 2] <- uniroot(foo, lower = para[2], upper = up)$root
56 names(para) <- names(se) <- rownames(CI) <- c("d/b", "b-d")
57 colnames(CI) <- c("lo", "up")
58 obj <- list(tree = deparse(substitute(phy)), N = N,
59 dev = Dev, para = para, se = se, CI = CI)
60 class(obj) <- "birthdeath"
64 print.birthdeath <- function(x, ...)
66 cat("\nEstimation of Speciation and Extinction Rates\n")
67 cat(" with Birth-Death Models\n\n")
68 cat(" Phylogenetic tree:", x$tree, "\n")
69 cat(" Number of tips:", x$N, "\n")
70 cat(" Deviance:", x$dev, "\n")
71 cat(" Log-likelihood:", -(x$dev)/2, "\n")
72 cat(" Parameter estimates:\n")
73 cat(" d / b =", x$para[1], " StdErr =", x$se[1], "\n")
74 cat(" b - d =", x$para[2], " StdErr =", x$se[2], "\n")
75 cat(" (b: speciation rate, d: extinction rate)\n")
76 cat(" Profile likelihood 95% confidence intervals:\n")
77 cat(" d / b: [", x$CI[1, 1], ", ", x$CI[1, 2], "]", "\n", sep = "")
78 cat(" b - d: [", x$CI[2, 1], ", ", x$CI[2, 2], "]", "\n\n", sep = "")
81 bd.ext <- function(phy, S)
83 if (class(phy) != "phylo") stop('object "phy" is not of class "phylo"')
84 if (!is.null(names(S))) {
85 if (all(names(S) %in% phy$tip.label)) S <- S[phy$tip.label]
86 else warning('the names of argument "S" and the names of the tip labels
87 did not match: the former were ignored in the analysis.')
90 x <- branching.times(phy)
92 trm.br <- phy$edge.length[phy$edge[, 2] <= N]
95 -2 * (lfactorial(N - 1)
97 + (3 * N) * log(1 - a)
99 - 2 * sum(log(exp(r * x[2:N]) - a))
101 + sum((S - 1) * log(exp(r * trm.br) - 1))
102 - sum((S + 1) * log(exp(r * trm.br) - a)))
104 out <- nlm(function(p) dev(p[1], p[2]), c(0, 0.2), hessian = TRUE)
105 if (out$estimate[1] < 0) {
106 out <- nlm(function(p) dev(0, p), 0.2, hessian = TRUE)
107 para <- c(0, out$estimate)
108 se <- c(0, sqrt(diag(solve(out$hessian))))
112 se <- sqrt(diag(solve(out$hessian)))
115 cat("\nExtended Version of the Birth-Death Models to\n")
116 cat(" Estimate Speciation and Extinction Rates\n\n")
117 cat(" Data: phylogenetic:", deparse(substitute(phy)), "\n")
118 cat(" taxonomic:", deparse(substitute(S)), "\n")
119 cat(" Number of tips:", N, "\n")
120 cat(" Deviance:", Dev, "\n")
121 cat(" Log-likelihood:", -Dev/2, "\n")
122 cat(" Parameter estimates:\n")
123 cat(" d / b =", para[1], " StdErr =", se[1], "\n")
124 cat(" b - d =", para[2], " StdErr =", se[2], "\n")
125 cat(" (b: speciation rate, d: extinction rate)\n")