R/compar.ou.R

   1 ## compar.ou.R (2009-05-10)
   2
   3 ##   Ornstein--Uhlenbeck Model for Continuous Characters
   4
   5 ## Copyright 2005-2009 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 compar.ou <- function(x, phy, node = NULL, alpha = NULL)
  11 {
  12     if (!inherits(phy, "phylo"))
  13       stop('object "phy" is not of class "phylo".')
  14     if (!is.numeric(x)) stop("'x' must be numeric.")
  15     if (!is.null(names(x))) {
  16         if (all(names(x) %in% phy$tip.label)) x <- x[phy$tip.label]
  17         else warning('the names of argument "x" and the tip labels of the tree did not match: the former were ignored in the analysis.')
  18     }
  19     nb.tip <- length(phy$tip.label)
  20     root <- nb.tip + 1
  21     if (is.null(node)) node <- numeric(0)
  22     if (root %in% node) node <- node[node != root]
  23     bt <- branching.times(phy)
  24     Tmax <- bt[1]
  25     Wend <- matrix(0, nb.tip, length(node) + 1)
  26     colnames(Wend) <- c(names(sort(bt[node])), as.character(root))
  27     Wstart <- Wend
  28     Wstart[, ncol(Wstart)] <- Tmax
  29     root2tip <- .Call("seq_root2tip", phy$edge, nb.tip,
  30                       phy$Nnode, PACKAGE = "ape")
  31     for (i in 1:nb.tip) {
  32         last.change <- names(Tmax)
  33         for (j in root2tip[[i]]) {#[-1]) {# don't need to look at the root
  34             if (j %in% node) {
  35                 jb <- as.character(j)
  36                 Wend[i, last.change] <- Wstart[i, jb] <- bt[jb]
  37                 last.change <- jb
  38             }
  39         }
  40     }
  41     W <- cophenetic.phylo(phy)
  42     dev <- function(p) {
  43         ##M <- rowSums(exp(-p[1] * Wstart) - exp(-p[1] * Wend) * p[-(1:2)])
  44         ##M <- -rowSums(exp(-p[1] * Wstart) - exp(-p[1] * Wend) * p[-(1:2)])
  45         M <- rowSums((exp(-p[1] * Wend) - exp(-p[1] * Wstart)) * p[-(1:2)])
  46         V <- exp(-p[1]*W) * (1 - exp(-2*p[1]*(Tmax - W/2)))
  47         nb.tip*log(2*pi*p[2]) + log(det(V)) +
  48           (t(x - M) %*% chol2inv(V) %*% (x - M)) / p[2]
  49     }
  50     if (is.null(alpha))
  51       out <- nlm(function(p) dev(p),
  52                  p = c(0.1, 1, rep(mean(x), ncol(Wstart))),
  53                  hessian = TRUE)
  54     else
  55       out <- nlm(function(p) dev(c(alpha, p)),
  56                  p = c(1, rep(mean(x), ncol(Wstart))),
  57                  hessian = TRUE)
  58     para <- cbind(out$estimate, sqrt(diag(solve(out$hessian))))
  59     nms <- c("sigma2", paste("theta", 1:ncol(Wstart), sep = ""))
  60     if (is.null(alpha)) nms <- c("alpha", nms)
  61     dimnames(para) <- list(nms, c("estimate", "stderr"))
  62     obj <- list(deviance = out$minimum, para = para, call = match.call())
  63     class(obj) <- "compar.ou"
  64     obj
  65 }