X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Fbirthdeath.R;h=7d53fd29218613b8a3047a8182ceee2c8556b46b;hb=refs%2Fheads%2Fmaster;hp=4c65ecd9248e1b1ca497bdc16957f4bcbc208e4e;hpb=a8d86ed3981c64b6d20316d54d3ed4d02c2b887e;p=ape.git diff --git a/R/birthdeath.R b/R/birthdeath.R index 4c65ecd..7d53fd2 100644 --- a/R/birthdeath.R +++ b/R/birthdeath.R @@ -1,4 +1,4 @@ -## birthdeath.R (2010-10-17) +## birthdeath.R (2012-04-20) ## Estimation of Speciation and Extinction Rates ## with Birth-Death Models @@ -6,7 +6,7 @@ ## birthdeath: standard model ## bd.ext: extended version -## Copyright 2002-2010 Emmanuel Paradis +## Copyright 2002-2012 Emmanuel Paradis ## This file is part of the R-package `ape'. ## See the file ../COPYING for licensing issues. @@ -17,7 +17,7 @@ birthdeath <- function(phy) N <- length(phy$tip.label) x <- c(NA, branching.times(phy)) dev <- function(a, r) { - if (r <= 0) return(1e50) + if (r < 0 || a > 1) return(1e100) -2 * (lfactorial(N - 1) + (N - 2) * log(r) + r * sum(x[3:N]) @@ -28,32 +28,47 @@ birthdeath <- function(phy) 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, @@ -79,39 +94,47 @@ print.birthdeath <- function(x, ...) 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")