R/yule.R

   1 ## yule.R (2011-11-03)
   2
   3 ##     Fits Yule Model to a Phylogenetic Tree
   4
   5 ## yule: standard Yule model (constant birth rate)
   6 ## yule.cov: Yule model with covariates
   7
   8 ## Copyright 2003-2011 Emmanuel Paradis
   9
  10 ## This file is part of the R-package `ape'.
  11 ## See the file ../COPYING for licensing issues.
  12
  13 yule <- function(phy, use.root.edge = FALSE)
  14 {
  15     if (!is.binary.tree(phy))
  16         stop("tree must be dichotomous to fit the Yule model.")
  17
  18     X <- sum(phy$edge.length)
  19     nb.node <- phy$Nnode
  20
  21     if (!is.null(phy$root.edge) && use.root.edge) X <- X + phy$root.edge
  22     else nb.node <- nb.node - 1
  23
  24     lambda <- nb.node/X
  25     se <- lambda/sqrt(nb.node)
  26     loglik <- -lambda * X + lfactorial(phy$Nnode) + nb.node * log(lambda)
  27     obj <- list(lambda = lambda, se = se, loglik = loglik)
  28     class(obj) <- "yule"
  29     obj
  30 }
  31
  32 yule.cov <- function(phy, formula, data = NULL)
  33 {
  34     if (is.null(data)) data <- parent.frame()
  35     n <- length(phy$tip.label)
  36     nb.node <- phy$Nnode
  37     if (!is.null(phy$node.label)) phy$node.label <- NULL
  38     bt <- sort(branching.times(phy)) # branching times (from present to past)
  39     bt <- rev(bt) # branching times from past to present
  40     ni <- cumsum(rev(table(bt))) + 1
  41     X <- model.matrix(formula, data)
  42     Xi <- X[phy$edge[, 1], , drop = FALSE]
  43     Xj <- X[phy$edge[, 2], , drop = FALSE]
  44     dev <- function(b) {
  45         2 * sum(((1/(1 + exp(-(Xi %*% b)))) +
  46                  (1/(1 + exp(-(Xj %*% b)))))
  47                 * phy$edge.length/2) -
  48          2 * (sum(log(ni[-length(ni)])) +
  49               sum(log((1/(1 + exp(-(X[-(1:(n + 1)), , drop = FALSE] %*% b)))))))
  50     }
  51     out <- nlm(function(p) dev(p), p = c(rep(0, ncol(X) - 1), -1),
  52                hessian = TRUE)
  53     Dev <- out$minimum
  54     para <- matrix(NA, ncol(X), 2)
  55     para[, 1] <- out$estimate
  56     if (any(out$gradient == 0))
  57       warning("The likelihood gradient seems flat in at least one dimension (null gradient):\ncannot compute the standard-errors of the parameters.\n")
  58     else para[, 2] <- sqrt(diag(solve(out$hessian)))
  59     rownames(para) <- colnames(X)
  60     colnames(para) <- c("Estimate", "StdErr")
  61     ## fit the intercept-only model:
  62     X <- model.matrix(~ 1, data = data.frame(X))
  63     Xi <- X[phy$edge[, 1], , drop = FALSE]
  64     Xj <- X[phy$edge[, 2], , drop = FALSE]
  65     Dev.null <- nlm(function(p) dev(p), p = -1)$minimum
  66     cat("\n---- Yule Model with Covariates ----\n\n")
  67     cat("    Phylogenetic tree:", deparse(substitute(phy)), "\n")
  68     cat("       Number of tips:", n, "\n")
  69     cat("      Number of nodes:", nb.node, "\n")
  70     cat("             Deviance:", Dev, "\n")
  71     cat("       Log-likelihood:", -Dev/2, "\n\n")
  72     cat("  Parameter estimates:\n")
  73     print(para)
  74     cat("\n")
  75     cat("Null Deviance:", Dev.null, "\n")
  76     cat("  Test of the fitted model: ")
  77     chi <- Dev.null - Dev
  78     df <- nrow(para) - 1
  79     cat("chi^2 =", round(chi, 3), "  df =", df,
  80         "  P =", round(1 - pchisq(chi, df), 3), "\n")
  81 }