3 ## Fits Yule Model to a Phylogenetic Tree
5 ## yule: standard Yule model (constant birth rate)
6 ## yule.cov: Yule model with covariates
8 ## Copyright 2003-2011 Emmanuel Paradis
10 ## This file is part of the R-package `ape'.
11 ## See the file ../COPYING for licensing issues.
13 yule <- function(phy, use.root.edge = FALSE)
15 if (!is.binary.tree(phy))
16 stop("tree must be dichotomous to fit the Yule model.")
18 X <- sum(phy$edge.length)
21 if (!is.null(phy$root.edge) && use.root.edge) X <- X + phy$root.edge
22 else nb.node <- nb.node - 1
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)
32 yule.cov <- function(phy, formula, data = NULL)
34 if (is.null(data)) data <- parent.frame()
35 n <- length(phy$tip.label)
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]
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)))))))
51 out <- nlm(function(p) dev(p), p = c(rep(0, ncol(X) - 1), -1),
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")
75 cat("Null Deviance:", Dev.null, "\n")
76 cat(" Test of the fitted model: ")
79 cat("chi^2 =", round(chi, 3), " df =", df,
80 " P =", round(1 - pchisq(chi, df), 3), "\n")