R/birthdeath.R

   1 ## birthdeath.R (2010-11-04)
   2
   3 ##   Estimation of Speciation and Extinction Rates
   4 ##             with Birth-Death Models
   5
   6 ## birthdeath: standard model
   7 ## bd.ext: extended version
   8
   9 ## Copyright 2002-2010 Emmanuel Paradis
  10
  11 ## This file is part of the R-package `ape'.
  12 ## See the file ../COPYING for licensing issues.
  13
  14 birthdeath <- function(phy)
  15 {
  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) return(1e50)
  21         -2 * (lfactorial(N - 1)
  22               + (N - 2) * log(r)
  23               + r * sum(x[3:N])
  24               + N * log(1 - a)
  25               - 2 * sum(log(exp(r * x[2:N]) - a)))
  26     }
  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         se <- c(0, sqrt(diag(solve(out$hessian))))
  32     }
  33     else {
  34         para <- out$estimate
  35         se <- sqrt(diag(solve(out$hessian)))
  36     }
  37     Dev <- out$minimum
  38     ## compute the 95 % profile likelihood CIs
  39     ## (not very clean... but seems to work -- EP 2003-03-29)
  40     CI <- matrix(NA, 2, 2)
  41     foo <- function(p) dev(p, para[2]) - 3.84 - Dev
  42     inc <- 1e-2
  43     lo <- para[1] - inc
  44     up <- para[1] + inc
  45     while (foo(lo) < 0) lo <- lo - inc
  46     while (foo(up) < 0) up <- up + inc
  47     CI[1, 1] <- uniroot(foo, lower = lo, upper = para[1])$root
  48     if (CI[1, 1] < 0) CI[1, 1] <- 0
  49     CI[1, 2] <- uniroot(foo, lower = para[1], upper = up)$root
  50     foo <- function(p) dev(para[1], p) - 3.84 - Dev
  51     lo <- para[2] - inc
  52     up <- para[2] + inc
  53     while (foo(lo) < 0) lo <- lo - inc
  54     while (foo(up) < 0) up <- up + inc
  55     CI[2, 1] <- uniroot(foo, lower = lo, upper = para[2])$root
  56     CI[2, 2] <- uniroot(foo, lower = para[2], upper = up)$root
  57     names(para) <- names(se) <- rownames(CI) <- c("d/b", "b-d")
  58     colnames(CI) <- c("lo", "up")
  59     obj <- list(tree = deparse(substitute(phy)), N = N,
  60                 dev = Dev, para = para, se = se, CI = CI)
  61     class(obj) <- "birthdeath"
  62     obj
  63 }
  64
  65 print.birthdeath <- function(x, ...)
  66 {
  67     cat("\nEstimation of Speciation and Extinction Rates\n")
  68     cat("            with Birth-Death Models\n\n")
  69     cat("     Phylogenetic tree:", x$tree, "\n")
  70     cat("        Number of tips:", x$N, "\n")
  71     cat("              Deviance:", x$dev, "\n")
  72     cat("        Log-likelihood:", -(x$dev)/2, "\n")
  73     cat("   Parameter estimates:\n")
  74     cat("      d / b =", x$para[1], "  StdErr =", x$se[1], "\n")
  75     cat("      b - d =", x$para[2], "  StdErr =", x$se[2], "\n")
  76     cat("   (b: speciation rate, d: extinction rate)\n")
  77     cat("   Profile likelihood 95% confidence intervals:\n")
  78     cat("      d / b: [", x$CI[1, 1], ", ", x$CI[1, 2], "]", "\n", sep = "")
  79     cat("      b - d: [", x$CI[2, 1], ", ", x$CI[2, 2], "]", "\n\n", sep = "")
  80 }
  81
  82 bd.ext <- function(phy, S, conditional = TRUE)
  83 {
  84     if (!inherits(phy, "phylo"))
  85         stop('object "phy" is not of class "phylo"')
  86     if (!is.null(names(S))) {
  87         if (all(names(S) %in% phy$tip.label)) S <- S[phy$tip.label]
  88         else warning('the names of argument "S" and the tip labels
  89 did not match: the former were ignored.')
  90     }
  91     N <- length(S)
  92     x <- branching.times(phy)
  93     x <- c(x[1], x)
  94     trm.br <- phy$edge.length[phy$edge[, 2] <= N]
  95     if (conditional) {
  96         dev <- function(a, r) {
  97             if (a >= 1 || a < 0 || r <= 0) return(1e50)
  98             ert <- exp(r * trm.br)
  99             zeta <- (ert - 1)/(ert - a)
 100             -2 * (lfactorial(N - 1)
 101                   + (N - 2) * log(r)
 102                   + N * log(1 - a)
 103                   + 2 * r * sum(x[2:N])
 104                   - 2 * sum(log(exp(r * x[2:N]) - a))
 105                   + sum(log(1 - zeta) + (S - 1)*log(zeta)))
 106         }
 107     } else {
 108         dev <- function(a, r) {
 109             if (a >= 1 || a < 0 || r <= 0) return(1e50)
 110             -2 * (lfactorial(N - 1)
 111                   + (N - 2) * log(r)
 112                   + (3 * N) * log(1 - a)
 113                   + 2 * r * sum(x[2:N])
 114                   - 2 * sum(log(exp(r * x[2:N]) - a))
 115                   + r * sum(trm.br)
 116                   + sum((S - 1) * log(exp(r * trm.br) - 1))
 117                   - sum((S + 1) * log(exp(r * trm.br) - a)))
 118         }
 119     }
 120     out <- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE)
 121     para <- out$estimate
 122     se <- sqrt(diag(solve(out$hessian)))
 123     Dev <- out$minimum
 124     cat("\nExtended Version of the Birth-Death Models to\n")
 125     cat("    Estimate Speciation and Extinction Rates\n\n")
 126     cat("    Data: phylogenetic:", deparse(substitute(phy)), "\n")
 127     cat("             taxonomic:", deparse(substitute(S)), "\n")
 128     cat("        Number of tips:", N, "\n")
 129     cat("              Deviance:", Dev, "\n")
 130     cat("        Log-likelihood:", -Dev/2, "\n")
 131     cat("   Parameter estimates:\n")
 132     cat("      d / b =", para[1], "  StdErr =", se[1], "\n")
 133     cat("      b - d =", para[2], "  StdErr =", se[2], "\n")
 134     cat("   (b: speciation rate, d: extinction rate)\n")
 135 }