out <- nlm(function(p) dev(p), p = c(rep(0, ncol(X) - 1), -1),
hessian = TRUE)
}
out <- nlm(function(p) dev(p), p = c(rep(0, ncol(X) - 1), -1),
hessian = TRUE)
@@ -55,10+54,15 @@ yule.cov <- function(phy, formula, data = NULL)
para <- matrix(NA, ncol(X), 2)
para[, 1] <- out$estimate
if (any(out$gradient == 0))
para <- matrix(NA, ncol(X), 2)
para[, 1] <- out$estimate
if (any(out$gradient == 0))
- warning("The likelihood gradient seems flat in at least one dimension (null gradient):\ncannot compute the standard-errors of the transition rates.\n")
+ warning("The likelihood gradient seems flat in at least one dimension (null gradient):\ncannot compute the standard-errors of the parameters.\n")
else para[, 2] <- sqrt(diag(solve(out$hessian)))
rownames(para) <- colnames(X)
colnames(para) <- c("Estimate", "StdErr")
else para[, 2] <- sqrt(diag(solve(out$hessian)))
rownames(para) <- colnames(X)
colnames(para) <- c("Estimate", "StdErr")
+ ## fit the intercept-only model:
+ X <- model.matrix(~ 1, data = data.frame(X))
+ Xi <- X[phy$edge[, 1], , drop = FALSE]
+ Xj <- X[phy$edge[, 2], , drop = FALSE]
+ Dev.null <- nlm(function(p) dev(p), p = -1)$minimum