1 ## birthdeath.R (2012-04-20)
3 ## Estimation of Speciation and Extinction Rates
4 ## with Birth-Death Models
6 ## birthdeath: standard model
7 ## bd.ext: extended version
9 ## Copyright 2002-2012 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 (!inherits(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 if (r < 0 || a > 1) return(1e100)
21 -2 * (lfactorial(N - 1)
25 - 2 * sum(log(exp(r * x[2:N]) - a)))
27 out <- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE)
28 if (out$estimate[1] < 0) {
29 out <- nlm(function(p) dev(0, p), 0.2, hessian = TRUE)
30 para <- c(0, out$estimate)
31 inv.hessian <- try(solve(out$hessian))
33 if (class(inv.hessian) == "try-error") NA
34 else sqrt(diag(inv.hessian))
39 inv.hessian <- try(solve(out$hessian))
41 if (class(inv.hessian) == "try-error") c(NA, NA)
42 else sqrt(diag(inv.hessian))
46 ## 95% profile likelihood CIs
48 ## which: index of the parameter (1 or 2)
49 ## s: sign of the increment (-1 or +1)
50 foo <- function(which, s) {
55 bar <- function() dev(p, para[2])
58 bar <- function() dev(para[1], p)
62 while (bar() < Dev + 3.84) p <- p + s * i
69 CI <- mapply(foo, c(1, 2, 1, 2), c(-1, -1, 1, 1))
72 names(para) <- names(se) <- rownames(CI) <- c("d/b", "b-d")
73 colnames(CI) <- c("lo", "up")
74 obj <- list(tree = deparse(substitute(phy)), N = N,
75 dev = Dev, para = para, se = se, CI = CI)
76 class(obj) <- "birthdeath"
80 print.birthdeath <- function(x, ...)
82 cat("\nEstimation of Speciation and Extinction Rates\n")
83 cat(" with Birth-Death Models\n\n")
84 cat(" Phylogenetic tree:", x$tree, "\n")
85 cat(" Number of tips:", x$N, "\n")
86 cat(" Deviance:", x$dev, "\n")
87 cat(" Log-likelihood:", -(x$dev)/2, "\n")
88 cat(" Parameter estimates:\n")
89 cat(" d / b =", x$para[1], " StdErr =", x$se[1], "\n")
90 cat(" b - d =", x$para[2], " StdErr =", x$se[2], "\n")
91 cat(" (b: speciation rate, d: extinction rate)\n")
92 cat(" Profile likelihood 95% confidence intervals:\n")
93 cat(" d / b: [", x$CI[1, 1], ", ", x$CI[1, 2], "]", "\n", sep = "")
94 cat(" b - d: [", x$CI[2, 1], ", ", x$CI[2, 2], "]", "\n\n", sep = "")
97 bd.ext <- function(phy, S, conditional = TRUE)
99 if (!inherits(phy, "phylo"))
100 stop('object "phy" is not of class "phylo"')
101 if (!is.null(names(S))) {
102 if (all(names(S) %in% phy$tip.label)) S <- S[phy$tip.label]
103 else warning('the names of argument "S" and the tip labels
104 did not match: the former were ignored.')
107 x <- branching.times(phy)
109 trm.br <- phy$edge.length[phy$edge[, 2] <= N]
111 dev <- function(a, r) {
112 if (a >= 1 || a < 0 || r <= 0) return(1e50)
113 ert <- exp(r * trm.br)
114 zeta <- (ert - 1)/(ert - a)
115 -2 * (lfactorial(N - 1)
118 + 2 * r * sum(x[2:N])
119 - 2 * sum(log(exp(r * x[2:N]) - a))
120 + sum(log(1 - zeta) + (S - 1)*log(zeta)))
123 dev <- function(a, r) {
124 if (a >= 1 || a < 0 || r <= 0) return(1e50)
125 -2 * (lfactorial(N - 1)
127 + (3 * N) * log(1 - a)
128 + 2 * r * sum(x[2:N])
129 - 2 * sum(log(exp(r * x[2:N]) - a))
131 + sum((S - 1) * log(exp(r * trm.br) - 1))
132 - sum((S + 1) * log(exp(r * trm.br) - a)))
135 out <- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE)
137 se <- sqrt(diag(solve(out$hessian)))
139 cat("\nExtended Version of the Birth-Death Models to\n")
140 cat(" Estimate Speciation and Extinction Rates\n\n")
141 cat(" Data: phylogenetic:", deparse(substitute(phy)), "\n")
142 cat(" taxonomic:", deparse(substitute(S)), "\n")
143 cat(" Number of tips:", N, "\n")
144 cat(" Deviance:", Dev, "\n")
145 cat(" Log-likelihood:", -Dev/2, "\n")
146 cat(" Parameter estimates:\n")
147 cat(" d / b =", para[1], " StdErr =", se[1], "\n")
148 cat(" b - d =", para[2], " StdErr =", se[2], "\n")
149 cat(" (b: speciation rate, d: extinction rate)\n")