Nothing
## diversi.time.R (2007-09-22)
## Analysis of Diversification with Survival Models
## Copyright 2002-2007 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
diversi.time <- function(x, census = NULL, censoring.codes = c(1, 0),
Tc = NULL)
{
n <- length(x)
if (is.null(census)) {
k <- n
census <- rep(censoring.codes[1], n)
}
else k <- sum(census == censoring.codes[1])
u <- n - k
S <- sum(x)
delta <- k / S
var.delta <- delta^2 / k
loglik.A <- k * log(delta) - delta * S
tk <- x[census == censoring.codes[1]]
tu <- x[census == censoring.codes[2]]
fb <- function(b)
1/b - sum(x^b * log(x))/sum(x^b) + sum(log(tk))/k
beta <- uniroot(fb, interval = c(1e-7, 10))$root
Sp <- sum(x^beta)
alpha <- (k / Sp)^(1/beta)
var.alpha <- 1/ ((k * beta / alpha^2) + beta * (beta - 1) * alpha^(beta - 2) * Sp)
ax <- alpha * x
var.beta <- 1 / (k / beta^2 + sum(ax^beta * log(ax)))
loglik.B <- k*(log(alpha) + log(beta)) +
(beta - 1)*(k*log(alpha) + sum(log(tk)))- Sp * alpha^beta
if (is.null(Tc)) Tc <- median(x)
tk1 <- tk[tk < Tc]
tk2 <- tk[tk >= Tc]
tu1 <- tu[tu < Tc]
tu2 <- tu[tu >= Tc]
k1 <- length(tk1)
k2 <- k - k1
u1 <- length(tu1)
u2 <- u - u1
tmp <- (k2 + u2) * Tc
delta1 <- k1 / (sum(tk1) + sum(tu1) + tmp)
delta2 <- k2 / (sum(tk2) + sum(tu2) - tmp)
var.delta1 <- delta1^2 / k1
var.delta2 <- delta2^2 / k2
tmp <- Tc * (delta2 - delta1)
loglik.C <- k1 * log(delta1) - delta1 * sum(tk1) + k2 * log(delta2) +
k2 * tmp - delta2 * sum(tk2) - delta1 * sum(tu1) +
u2 * tmp - delta2 * sum(tu2)
cat("\nAnalysis of Diversification with Survival Models\n\n")
cat("Data:", deparse(substitute(x)), "\n")
cat("Number of branching times:", n, "\n")
cat(" accurately known:", k, "\n")
cat(" censored:", u, "\n\n")
cat("Model A: constant diversification\n")
cat(" log-likelihood =", round(loglik.A, 3),
" AIC =", round(-2 * loglik.A + 2, 3), "\n")
cat(" delta =", round(delta, 6), " StdErr =", round(sqrt(var.delta), 6), "\n\n")
cat("Model B: diversification follows a Weibull law\n")
cat(" log-likelihood =", round(loglik.B, 3),
" AIC =", round(-2 * loglik.B + 4, 3), "\n")
cat(" alpha =", round(alpha, 6), " StdErr =", round(sqrt(var.alpha), 6), "\n")
cat(" beta =", round(beta, 6), " StdErr =", round(sqrt(var.beta), 6), "\n\n")
cat("Model C: diversification changes with a breakpoint at time =", Tc, "\n")
cat(" log-likelihood =", round(loglik.C, 3),
" AIC =", round(-2 * loglik.C + 4, 3), "\n")
cat(" delta1 =", round(delta1, 6), " StdErr =", round(sqrt(var.delta1), 6), "\n")
cat(" delta2 =", round(delta2, 6), " StdErr =", round(sqrt(var.delta2), 6), "\n\n")
cat("Likelihood ratio tests:\n")
c1 <- 2 * (loglik.B - loglik.A)
p1 <- round(1 - pchisq(c1, 1), 4)
c2 <- 2 * (loglik.C - loglik.A)
p2 <- round(1 - pchisq(c2, 1), 4)
cat(" Model A vs. Model B: chi^2 =", round(c1, 3), " df = 1, P =", p1, "\n")
cat(" Model A vs. Model C: chi^2 =", round(c2, 3), " df = 1, P =", p2, "\n")
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.