R/birthdeath.R

   1 ## birthdeath.R (2012-04-20)
   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-2012 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 || a > 1) return(1e100)
  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         inv.hessian <- try(solve(out$hessian))
  32         se <-
  33             if (class(inv.hessian) == "try-error") NA
  34             else sqrt(diag(inv.hessian))
  35         se <- c(0, se)
  36     }
  37     else {
  38         para <- out$estimate
  39         inv.hessian <- try(solve(out$hessian))
  40         se <-
  41             if (class(inv.hessian) == "try-error") c(NA, NA)
  42             else sqrt(diag(inv.hessian))
  43     }
  44     Dev <- out$minimum
  45
  46     ## 95% profile likelihood CIs
  47
  48     ## which: index of the parameter (1 or 2)
  49     ## s: sign of the increment (-1 or +1)
  50     foo <- function(which, s) {
  51         i <- 0.1
  52
  53         if (which == 1) {
  54             p <- para[1] + s * i
  55             bar <- function() dev(p, para[2])
  56         } else { # which == 2
  57             p <- para[2] + s * i
  58             bar <- function() dev(para[1], p)
  59         }
  60
  61         while (i > 1e-9) {
  62             while (bar() < Dev + 3.84) p <- p + s * i
  63             p <- p - s * i
  64             i <- i / 10
  65         }
  66         p
  67     }
  68
  69     CI <- mapply(foo, c(1, 2, 1, 2), c(-1, -1, 1, 1))
  70     dim(CI) <- c(2, 2)
  71
  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"
  77     obj
  78 }
  79
  80 print.birthdeath <- function(x, ...)
  81 {
  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 = "")
  95 }
  96
  97 bd.ext <- function(phy, S, conditional = TRUE)
  98 {
  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.')
 105     }
 106     N <- length(S)
 107     x <- branching.times(phy)
 108     x <- c(x[1], x)
 109     trm.br <- phy$edge.length[phy$edge[, 2] <= N]
 110     if (conditional) {
 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)
 116                   + (N - 2) * log(r)
 117                   + N * log(1 - a)
 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)))
 121         }
 122     } else {
 123         dev <- function(a, r) {
 124             if (a >= 1 || a < 0 || r <= 0) return(1e50)
 125             -2 * (lfactorial(N - 1)
 126                   + (N - 2) * log(r)
 127                   + (3 * N) * log(1 - a)
 128                   + 2 * r * sum(x[2:N])
 129                   - 2 * sum(log(exp(r * x[2:N]) - a))
 130                   + r * sum(trm.br)
 131                   + sum((S - 1) * log(exp(r * trm.br) - 1))
 132                   - sum((S + 1) * log(exp(r * trm.br) - a)))
 133         }
 134     }
 135     out <- nlm(function(p) dev(p[1], p[2]), c(0.1, 0.2), hessian = TRUE)
 136     para <- out$estimate
 137     se <- sqrt(diag(solve(out$hessian)))
 138     Dev <- out$minimum
 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")
 150 }