R/diversi.time.R

   1 ## diversi.time.R (2007-09-22)
   2
   3 ##   Analysis of Diversification with Survival Models
   4
   5 ## Copyright 2002-2007 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 diversi.time <- function(x, census = NULL, censoring.codes = c(1, 0),
  11                          Tc = NULL)
  12 {
  13     n <- length(x)
  14     if (is.null(census)) {
  15         k <- n
  16         census <- rep(censoring.codes[1], n)
  17     }
  18     else k <- sum(census == censoring.codes[1])
  19     u <- n - k
  20     S <- sum(x)
  21     delta <- k / S
  22     var.delta <- delta^2 / k
  23     loglik.A <- k * log(delta) - delta * S
  24     tk <- x[census == censoring.codes[1]]
  25     tu <- x[census == censoring.codes[2]]
  26     fb <- function(b)
  27       1/b - sum(x^b * log(x))/sum(x^b) + sum(log(tk))/k
  28     beta <- uniroot(fb, interval = c(1e-7, 10))$root
  29     Sp <- sum(x^beta)
  30     alpha <- (k / Sp)^(1/beta)
  31     var.alpha <- 1/ ((k * beta / alpha^2) + beta * (beta - 1) * alpha^(beta - 2) * Sp)
  32     ax <- alpha * x
  33     var.beta <- 1 / (k / beta^2 + sum(ax^beta * log(ax)))
  34     loglik.B <- k*(log(alpha) + log(beta)) +
  35       (beta - 1)*(k*log(alpha) + sum(log(tk)))- Sp * alpha^beta
  36     if (is.null(Tc)) Tc <- median(x)
  37     tk1 <- tk[tk < Tc]
  38     tk2 <- tk[tk >= Tc]
  39     tu1 <- tu[tu < Tc]
  40     tu2 <- tu[tu >= Tc]
  41     k1 <- length(tk1)
  42     k2 <- k - k1
  43     u1 <- length(tu1)
  44     u2 <- u - u1
  45     tmp <- (k2 + u2) * Tc
  46     delta1 <- k1 / (sum(tk1) + sum(tu1) + tmp)
  47     delta2 <- k2 / (sum(tk2) + sum(tu2) - tmp)
  48     var.delta1 <- delta1^2 / k1
  49     var.delta2 <- delta2^2 / k2
  50     tmp <- Tc * (delta2 - delta1)
  51     loglik.C <- k1 * log(delta1) - delta1 * sum(tk1) + k2 * log(delta2) +
  52                   k2 * tmp - delta2 * sum(tk2) - delta1 * sum(tu1) +
  53                     u2 * tmp - delta2 * sum(tu2)
  54     cat("\nAnalysis of Diversification with Survival Models\n\n")
  55     cat("Data:", deparse(substitute(x)), "\n")
  56     cat("Number of branching times:", n, "\n")
  57     cat("         accurately known:", k, "\n")
  58     cat("                 censored:", u, "\n\n")
  59     cat("Model A: constant diversification\n")
  60     cat("    log-likelihood =", round(loglik.A, 3),
  61         "   AIC =", round(-2 * loglik.A + 2, 3), "\n")
  62     cat("    delta =", round(delta, 6), "   StdErr =", round(sqrt(var.delta), 6), "\n\n")
  63     cat("Model B: diversification follows a Weibull law\n")
  64     cat("    log-likelihood =", round(loglik.B, 3),
  65         "   AIC =", round(-2 * loglik.B + 4, 3), "\n")
  66     cat("    alpha =", round(alpha, 6), "   StdErr =", round(sqrt(var.alpha), 6), "\n")
  67     cat("     beta =", round(beta, 6), "   StdErr =", round(sqrt(var.beta), 6), "\n\n")
  68     cat("Model C: diversification changes with a breakpoint at time =", Tc, "\n")
  69     cat("    log-likelihood =", round(loglik.C, 3),
  70         "   AIC =", round(-2 * loglik.C + 4, 3), "\n")
  71     cat("    delta1 =", round(delta1, 6), "   StdErr =", round(sqrt(var.delta1), 6), "\n")
  72     cat("    delta2 =", round(delta2, 6), "   StdErr =", round(sqrt(var.delta2), 6), "\n\n")
  73     cat("Likelihood ratio tests:\n")
  74     c1 <- 2 * (loglik.B - loglik.A)
  75     p1 <- round(1 - pchisq(c1, 1), 4)
  76     c2 <- 2 * (loglik.C - loglik.A)
  77     p2 <- round(1 - pchisq(c2, 1), 4)
  78     cat("    Model A vs. Model B: chi^2 =", round(c1, 3), "   df = 1,    P =", p1, "\n")
  79     cat("    Model A vs. Model C: chi^2 =", round(c2, 3), "   df = 1,    P =", p2, "\n")
  80 }