-## birthdeath.R (2009-05-10)
+## birthdeath.R (2012-04-20)
## Estimation of Speciation and Extinction Rates
## with Birth-Death Models
## birthdeath: standard model
## bd.ext: extended version
-## Copyright 2002-2009 Emmanuel Paradis
+## Copyright 2002-2012 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
N <- length(phy$tip.label)
x <- c(NA, branching.times(phy))
dev <- function(a, r) {
+ if (r < 0 || a > 1) return(1e100)
-2 * (lfactorial(N - 1)
+ (N - 2) * log(r)
+ r * sum(x[3:N])
if (out$estimate[1] < 0) {
out <- nlm(function(p) dev(0, p), 0.2, hessian = TRUE)
para <- c(0, out$estimate)
- se <- c(0, sqrt(diag(solve(out$hessian))))
+ inv.hessian <- try(solve(out$hessian))
+ se <-
+ if (class(inv.hessian) == "try-error") NA
+ else sqrt(diag(inv.hessian))
+ se <- c(0, se)
}
else {
para <- out$estimate
- se <- sqrt(diag(solve(out$hessian)))
+ inv.hessian <- try(solve(out$hessian))
+ se <-
+ if (class(inv.hessian) == "try-error") c(NA, NA)
+ else sqrt(diag(inv.hessian))
}
Dev <- out$minimum
- ## compute the 95 % profile likelihood CIs
- ## (not very clean... but seems to work -- EP 2003-03-29)
- CI <- matrix(NA, 2, 2)
- foo <- function(p) dev(p, para[2]) - 3.84 - Dev
- inc <- 1e-2
- lo <- para[1] - inc
- up <- para[1] + inc
- while (foo(lo) < 0) lo <- lo - inc
- while (foo(up) < 0) up <- up + inc
- CI[1, 1] <- uniroot(foo, lower = lo, upper = para[1])$root
- if (CI[1, 1] < 0) CI[1, 1] <- 0
- CI[1, 2] <- uniroot(foo, lower = para[1], upper = up)$root
- foo <- function(p) dev(para[1], p) - 3.84 - Dev
- lo <- para[2] - inc
- up <- para[2] + inc
- while (foo(lo) < 0) lo <- lo - inc
- while (foo(up) < 0) up <- up + inc
- CI[2, 1] <- uniroot(foo, lower = lo, upper = para[2])$root
- CI[2, 2] <- uniroot(foo, lower = para[2], upper = up)$root
+
+ ## 95% profile likelihood CIs
+
+ ## which: index of the parameter (1 or 2)
+ ## s: sign of the increment (-1 or +1)
+ foo <- function(which, s) {
+ i <- 0.1
+
+ if (which == 1) {
+ p <- para[1] + s * i
+ bar <- function() dev(p, para[2])
+ } else { # which == 2
+ p <- para[2] + s * i
+ bar <- function() dev(para[1], p)
+ }
+
+ while (i > 1e-9) {
+ while (bar() < Dev + 3.84) p <- p + s * i
+ p <- p - s * i
+ i <- i / 10
+ }
+ p
+ }
+
+ CI <- mapply(foo, c(1, 2, 1, 2), c(-1, -1, 1, 1))
+ dim(CI) <- c(2, 2)
+
names(para) <- names(se) <- rownames(CI) <- c("d/b", "b-d")
colnames(CI) <- c("lo", "up")
obj <- list(tree = deparse(substitute(phy)), N = N,
cat(" b - d: [", x$CI[2, 1], ", ", x$CI[2, 2], "]", "\n\n", sep = "")
}
-bd.ext <- function(phy, S)
+bd.ext <- function(phy, S, conditional = TRUE)
{
- if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"')
+ if (!inherits(phy, "phylo"))
+ stop('object "phy" is not of class "phylo"')
if (!is.null(names(S))) {
if (all(names(S) %in% phy$tip.label)) S <- S[phy$tip.label]
- else warning('the names of argument "S" and the names of the tip labels
-did not match: the former were ignored in the analysis.')
+ else warning('the names of argument "S" and the tip labels
+did not match: the former were ignored.')
}
N <- length(S)
x <- branching.times(phy)
x <- c(x[1], x)
trm.br <- phy$edge.length[phy$edge[, 2] <= N]
- dev <- function(a, r)
- {
- -2 * (lfactorial(N - 1)
- + (N - 2) * log(r)
- + (3 * N) * log(1 - a)
- + 2 * r * sum(x[2:N])
- - 2 * sum(log(exp(r * x[2:N]) - a))
- + r * sum(trm.br)
- + sum((S - 1) * log(exp(r * trm.br) - 1))
- - sum((S + 1) * log(exp(r * trm.br) - a)))
- }
- out <- nlm(function(p) dev(p[1], p[2]), c(0, 0.2), hessian = TRUE)
- if (out$estimate[1] < 0) {
- out <- nlm(function(p) dev(0, p), 0.2, hessian = TRUE)
- para <- c(0, out$estimate)
- se <- c(0, sqrt(diag(solve(out$hessian))))
- }
- else {
- para <- out$estimate
- se <- sqrt(diag(solve(out$hessian)))
+ if (conditional) {
+ dev <- function(a, r) {
+ if (a >= 1 || a < 0 || r <= 0) return(1e50)
+ ert <- exp(r * trm.br)
+ zeta <- (ert - 1)/(ert - a)
+ -2 * (lfactorial(N - 1)
+ + (N - 2) * log(r)
+ + N * log(1 - a)
+ + 2 * r * sum(x[2:N])
+ - 2 * sum(log(exp(r * x[2:N]) - a))
+ + sum(log(1 - zeta) + (S - 1)*log(zeta)))
+ }
+ } else {
+ dev <- function(a, r) {
+ if (a >= 1 || a < 0 || r <= 0) return(1e50)
+ -2 * (lfactorial(N - 1)
+ + (N - 2) * log(r)
+ + (3 * N) * log(1 - a)
+ + 2 * r * sum(x[2:N])
+ - 2 * sum(log(exp(r * x[2:N]) - a))
+ + r * sum(trm.br)
+ + sum((S - 1) * log(exp(r * trm.br) - 1))
+ - sum((S + 1) * log(exp(r * trm.br) - a)))
+ }
}
+ out <- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE)
+ para <- out$estimate
+ se <- sqrt(diag(solve(out$hessian)))
Dev <- out$minimum
cat("\nExtended Version of the Birth-Death Models to\n")
cat(" Estimate Speciation and Extinction Rates\n\n")