X-Git-Url: https://git.donarmstrong.com/?a=blobdiff_plain;f=R%2Fcompar.ou.R;h=2e1d3a512dbf63ce810542b902468ef04193032d;hb=d5b85d181648e2761cecc91d75c4c66fa05e4508;hp=83e0d1cf5ab5eb0d939caac0a370a9884bc153d6;hpb=c827059eeafc8cbe41c812b26979543ab287803e;p=ape.git

diff --git a/R/compar.ou.R b/R/compar.ou.R
index 83e0d1c..2e1d3a5 100644
--- a/R/compar.ou.R
+++ b/R/compar.ou.R
@@ -1,36 +1,42 @@
-## compar.ou.R (2006-10-05)
+## compar.ou.R (2010-11-04)
 
 ##   Ornstein--Uhlenbeck Model for Continuous Characters
 
-## Copyright 2005-2006 Emmanuel Paradis
+## Copyright 2005-2010 Emmanuel Paradis
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
 
 compar.ou <- function(x, phy, node = NULL, alpha = NULL)
 {
-    if (class(phy) != "phylo")
-      stop('object "phy" is not of class "phylo".')
+    if (!inherits(phy, "phylo"))
+        stop('object "phy" is not of class "phylo".')
     if (!is.numeric(x)) stop("'x' must be numeric.")
     if (!is.null(names(x))) {
         if (all(names(x) %in% phy$tip.label)) x <- x[phy$tip.label]
         else warning('the names of argument "x" and the tip labels of the tree did not match: the former were ignored in the analysis.')
     }
-    nb.tip <- length(phy$tip.label)
-    root <- nb.tip + 1
+    n <- length(phy$tip.label)
+    root <- n + 1L
     if (is.null(node)) node <- numeric(0)
-    if (root %in% node) node <- node[node != root]
+    if (is.character(node)) {
+        if (is.null(phy$node.label))
+            stop("argument 'node' is character but 'phy' has no node labels")
+        node <- match(node, phy$node.label) + n
+        phy$node.label <- NULL
+    }
+    if (root %in% node) node <- node[-1]
     bt <- branching.times(phy)
     Tmax <- bt[1]
-    Wend <- matrix(0, nb.tip, length(node) + 1)
-    colnames(Wend) <- c(names(sort(bt[node])), as.character(root))
+    Wend <- matrix(0, n, length(node) + 1)
+    colnames(Wend) <- c(names(sort(bt[node - n])), as.character(root))
     Wstart <- Wend
     Wstart[, ncol(Wstart)] <- Tmax
-    root2tip <- .Call("seq_root2tip", phy$edge, nb.tip,
+    root2tip <- .Call("seq_root2tip", phy$edge, n,
                       phy$Nnode, PACKAGE = "ape")
-    for (i in 1:nb.tip) {
+    for (i in 1:n) {
         last.change <- names(Tmax)
-        for (j in root2tip[[i]]) {#[-1]) {# don't need to look at the root
+        for (j in root2tip[[i]]) {
             if (j %in% node) {
                 jb <- as.character(j)
                 Wend[i, last.change] <- Wstart[i, jb] <- bt[jb]
@@ -40,20 +46,31 @@ compar.ou <- function(x, phy, node = NULL, alpha = NULL)
     }
     W <- cophenetic.phylo(phy)
     dev <- function(p) {
-        M <- rowSums(exp(-p[1] * Wstart) - exp(-p[1] * Wend) * p[-(1:2)])
-        V <- exp(-p[1]*W) * (1 - exp(-2*p[1]*(Tmax - W/2)))
-        nb.tip*log(2*pi*p[2]) + log(det(V)) +
-          (t(x - M) %*% chol2inv(V) %*% (x - M)) / p[2]
+        alpha <- p[1]
+        sigma2 <- p[2]
+        if (sigma2 <= 0) return(1e100)
+        theta <- p[-(1:2)]
+        ## fixed a bug below: must be '%*% theta' instead of '* theta' (2010-03-15)
+        M <- rowSums((exp(-alpha * Wend) - exp(-alpha * Wstart)) %*% theta)
+        V <- exp(-alpha * W) * (1 - exp(-2 * alpha * (Tmax - W/2)))
+        R <- chol(V) # correction by Cecile Ane (2010-11-04)
+        n * log(2 * pi * sigma2) + 2 * sum(log(diag(R))) +
+            (t(x - M) %*% chol2inv(R) %*% (x - M)) / sigma2
     }
-    if (is.null(alpha))
-      out <- nlm(function(p) dev(p),
-                 p = c(0.1, 1, rep(mean(x), ncol(Wstart))),
-                 hessian = TRUE)
-    else
-      out <- nlm(function(p) dev(c(alpha, p)),
-                 p = c(1, rep(mean(x), ncol(Wstart))),
-                 hessian = TRUE)
-    para <- cbind(out$estimate, sqrt(diag(solve(out$hessian))))
+
+    out <- if (is.null(alpha))
+        nlm(function(p) dev(p),
+            p = c(0.1, 1, rep(mean(x), ncol(Wstart))), hessian = TRUE)
+    else nlm(function(p) dev(c(alpha, p)),
+             p = c(1, rep(mean(x), ncol(Wstart))), hessian = TRUE)
+
+    ## if alpha is estimated it may be that the Hessian matrix has the
+    ## corresponding column and row filled with 0, making solve() fail
+    se <- if (is.null(alpha) && all(out$hessian[1, ] == 0))
+        c(NA, sqrt(diag(solve(out$hessian[-1, -1]))))
+    else sqrt(diag(solve(out$hessian)))
+
+    para <- cbind(out$estimate, se)
     nms <- c("sigma2", paste("theta", 1:ncol(Wstart), sep = ""))
     if (is.null(alpha)) nms <- c("alpha", nms)
     dimnames(para) <- list(nms, c("estimate", "stderr"))