better error-catching with ace()

git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@191 6e262413-ae40-0410-9e79-b911bd7a66b7

diff --git a/NEWS b/NEWS
index ec8918ae97f2c4e7d9bfe7e8a7d19fc2a6599041..1c1ae9cab6573ea0ce243036ada8d6b38b25b11c 100644 (file)
--- a/NEWS
+++ b/NEWS
@@ -1,6 +1,11 @@
                CHANGES IN APE VERSION 3.0-5
 
 
+BUG FIXES
+
+    o ace() should better catch errors when SEs cannot be computed.
+
+
 OTHER CHANGES
 
     o write.dna(format = "fasta") now conforms more closely to the

diff --git a/R/ace.R b/R/ace.R
index ed86fc6d6af481104c7bde2816a83f5e9b9bccab..55e32c35ebc1e3c2580ce0af1be4d04650c53de4 100644 (file)
--- a/R/ace.R
+++ b/R/ace.R
@@ -1,12 +1,24 @@
-## ace.R (2011-07-18)
+## ace.R (2012-06-25)
 
 ##   Ancestral Character Estimation
 
-## Copyright 2005-2011 Emmanuel Paradis and Ben Bolker
+## Copyright 2005-2012 Emmanuel Paradis and Ben Bolker
 
 ## This file is part of the R-package `ape'.
 ## See the file ../COPYING for licensing issues.
 
+.getSEs <- function(out)
+{
+    h <- out$hessian
+    if (any(diag(h) == 0)) {
+        warning("The likelihood gradient seems flat in at least one dimension (gradient null):\ncannot compute the standard-errors of the transition rates.\n")
+        se <- rep(NaN, nrow(h))
+    } else {
+        se <- sqrt(diag(solve(h)))
+    }
+    se
+}
+
 ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE,
                 model = if (type == "continuous") "BM" else "ER",
                 scaled = TRUE, kappa = 1, corStruct = NULL, ip = 0.1)
@@ -35,7 +47,7 @@ ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE,
                 if (sig2 < 0) return(1e100)
                 V <- sig2 * vcv(phy)
                 ## next three lines borrowed from dmvnorm() in 'mvtnorm'
-                distval <- mahalanobis(x, center = mu, cov = V)
+                distval <- stats::mahalanobis(x, center = mu, cov = V)
                 logdet <- sum(log(eigen(V, symmetric = TRUE, only.values = TRUE)$values))
                 (nb.tip * log(2 * pi) + logdet + distval)/2
             }
@@ -56,10 +68,10 @@ ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE,
                        p = rep(mu[1], nb.node), hessian = TRUE)
             obj$resloglik <- -out$minimum
             obj$ace <- out$estimate
-            names(obj$ace) <- (nb.tip + 1):(nb.tip + nb.node)
+            names(obj$ace) <- nb.tip + 1:nb.node
             obj$sigma2 <- c(sigma2, se_sgi2)
             if (CI) {
-                se <- sqrt(diag(solve(out$hessian)))
+                se <- .getSEs(out)
                 tmp <- se * qt(0.025, nb.node)
                 obj$CI95 <- cbind(obj$ace + tmp, obj$ace - tmp)
             }
@@ -78,7 +90,7 @@ ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE,
                       as.integer(CI), as.integer(scaled),
                       PACKAGE = "ape")
             obj$ace <- ans[[6]][-(1:nb.tip)]
-            names(obj$ace) <- (nb.tip + 1):(nb.tip + nb.node)
+            names(obj$ace) <- nb.tip + 1:nb.node
             if (CI) {
                 se <- sqrt(ans[[8]])
                 tmp <- se * qnorm(0.025)
@@ -101,7 +113,7 @@ ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE,
                 obj$loglik <- -out$minimum / 2
                 obj$ace <- out$estimate[-1]
                 names(obj$ace) <- (nb.tip + 1):(nb.tip + nb.node)
-                se <- sqrt(diag(solve(out$hessian)))
+                se <- .getSEs(out)
                 obj$sigma2 <- c(out$estimate[1], se[1])
                 if (CI) {
                     tmp <- se[-1] * qt(0.025, nb.node)
@@ -206,12 +218,10 @@ ace <- function(x, phy, type = "continuous", method = "ML", CI = TRUE,
         obj$rates <- out$par
         oldwarn <- options("warn")
         options(warn = -1)
-        h <- nlm(function(p) dev(p), p = obj$rates, iterlim = 1,
-                 stepmax = 0, hessian = TRUE)$hessian
+        out.nlm <- nlm(function(p) dev(p), p = obj$rates, iterlim = 1,
+                       stepmax = 0, hessian = TRUE)
         options(oldwarn)
-        if (any(h == 0))
-          warning("The likelihood gradient seems flat in at least one dimension (gradient null):\ncannot compute the standard-errors of the transition rates.\n")
-        else obj$se <- sqrt(diag(solve(h)))
+        obj$se <- .getSEs(out.nlm)
         obj$index.matrix <- index.matrix
         if (CI) {
             obj$lik.anc <- dev(obj$rates, TRUE)