1 ## diversi.time.R (2007-09-22)
3 ## Analysis of Diversification with Survival Models
5 ## Copyright 2002-2007 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 diversi.time <- function(x, census = NULL, censoring.codes = c(1, 0),
14 if (is.null(census)) {
16 census <- rep(censoring.codes[1], n)
18 else k <- sum(census == censoring.codes[1])
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]]
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
30 alpha <- (k / Sp)^(1/beta)
31 var.alpha <- 1/ ((k * beta / alpha^2) + beta * (beta - 1) * alpha^(beta - 2) * Sp)
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)
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")