R/compar.ou.R

   1 ## compar.ou.R (2010-11-04)
   2
   3 ##   Ornstein--Uhlenbeck Model for Continuous Characters
   4
   5 ## Copyright 2005-2010 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     n <- length(phy$tip.label)
  20     root <- n + 1L
  21     if (is.null(node)) node <- numeric(0)
  22     if (is.character(node)) {
  23         if (is.null(phy$node.label))
  24             stop("argument 'node' is character but 'phy' has no node labels")
  25         node <- match(node, phy$node.label) + n
  26         phy$node.label <- NULL
  27     }
  28     if (root %in% node) node <- node[-1]
  29     bt <- branching.times(phy)
  30     Tmax <- bt[1]
  31     Wend <- matrix(0, n, length(node) + 1)
  32     colnames(Wend) <- c(names(sort(bt[node - n])), as.character(root))
  33     Wstart <- Wend
  34     Wstart[, ncol(Wstart)] <- Tmax
  35     root2tip <- .Call("seq_root2tip", phy$edge, n,
  36                       phy$Nnode, PACKAGE = "ape")
  37     for (i in 1:n) {
  38         last.change <- names(Tmax)
  39         for (j in root2tip[[i]]) {
  40             if (j %in% node) {
  41                 jb <- as.character(j)
  42                 Wend[i, last.change] <- Wstart[i, jb] <- bt[jb]
  43                 last.change <- jb
  44             }
  45         }
  46     }
  47     W <- cophenetic.phylo(phy)
  48     dev <- function(p) {
  49         alpha <- p[1]
  50         sigma2 <- p[2]
  51         if (sigma2 <= 0) return(1e100)
  52         theta <- p[-(1:2)]
  53         ## fixed a bug below: must be '%*% theta' instead of '* theta' (2010-03-15)
  54         M <- rowSums((exp(-alpha * Wend) - exp(-alpha * Wstart)) %*% theta)
  55         V <- exp(-alpha * W) * (1 - exp(-2 * alpha * (Tmax - W/2)))
  56         R <- chol(V) # correction by Cecile Ane (2010-11-04)
  57         n * log(2 * pi * sigma2) + 2 * sum(log(diag(R))) +
  58             (t(x - M) %*% chol2inv(R) %*% (x - M)) / sigma2
  59     }
  60
  61     out <- if (is.null(alpha))
  62         nlm(function(p) dev(p),
  63             p = c(0.1, 1, rep(mean(x), ncol(Wstart))), hessian = TRUE)
  64     else nlm(function(p) dev(c(alpha, p)),
  65              p = c(1, rep(mean(x), ncol(Wstart))), hessian = TRUE)
  66
  67     ## if alpha is estimated it may be that the Hessian matrix has the
  68     ## corresponding column and row filled with 0, making solve() fail
  69     se <- if (is.null(alpha) && all(out$hessian[1, ] == 0))
  70         c(NA, sqrt(diag(solve(out$hessian[-1, -1]))))
  71     else sqrt(diag(solve(out$hessian)))
  72
  73     para <- cbind(out$estimate, se)
  74     nms <- c("sigma2", paste("theta", 1:ncol(Wstart), sep = ""))
  75     if (is.null(alpha)) nms <- c("alpha", nms)
  76     dimnames(para) <- list(nms, c("estimate", "stderr"))
  77     obj <- list(deviance = out$minimum, para = para, call = match.call())
  78     class(obj) <- "compar.ou"
  79     obj
  80 }