+ switch(method, "REML" = {
+ minusLogLik <- function(sig2) {
+ if (sig2 < 0) return(1e100)
+ V <- sig2 * vcv(phy)
+ ## next three lines borrowed from dmvnorm() in 'mvtnorm'
+ distval <- 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
+ }
+ mu <- rep(ace(x, phy, method="pic")$ace[1], nb.tip)
+ out <- nlm(minusLogLik, 1, hessian = TRUE)
+ sigma2 <- out$estimate
+ se_sgi2 <- sqrt(1/out$hessian)
+ tip <- phy$edge[, 2] <= nb.tip
+ minus.REML.BM <- function(p) {
+ x1 <- p[phy$edge[, 1] - nb.tip]
+ x2 <- numeric(length(x1))
+ x2[tip] <- x[phy$edge[tip, 2]]
+ x2[!tip] <- p[phy$edge[!tip, 2] - nb.tip]
+ -(-sum((x1 - x2)^2/phy$edge.length)/(2 * sigma2) -
+ nb.node * log(sigma2))
+ }
+ out <- nlm(function(p) minus.REML.BM(p),
+ 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)
+ obj$sigma2 <- c(sigma2, se_sgi2)
+ if (CI) {
+ se <- sqrt(diag(solve(out$hessian)))
+ tmp <- se * qt(0.025, nb.node)
+ obj$CI95 <- cbind(obj$ace + tmp, obj$ace - tmp)
+ }
+ }, "pic" = {