CRAN Package Check Results for Package phylocurve

Last updated on 2019-11-26 00:52:01 CET.

Check Details

Version: 2.0.9
Check: examples
Result: ERROR
    Running examples in 'phylocurve-Ex.R' failed
    The error most likely occurred in:
    
    > base::assign(".ptime", proc.time(), pos = "CheckExEnv")
    > ### Name: K.mult
    > ### Title: Test phylogenetic signal (Kmult) using phylogenetic simulation
    > ### Aliases: K.mult
    >
    > ### ** Examples
    >
    > rand.data <- sim.traits()
    > null.model <- evo.model(tree = rand.data$tree,
    + Y = rand.data$trait_data,method = "Pairwise ML")
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
    :
     --- package (from environment) ---
    phylocurve
     --- call from context ---
    evo.model(tree = rand.data$tree, Y = rand.data$trait_data, method = "Pairwise ML")
     --- call from argument ---
    if (class(Y) == "data.frame") {
     species_col <- match(species.id, colnames(Y))
     if (is.na(species_col))
     stop("Name of species.id much match a column in Y (if supplying a data frame).")
     trait_names <- colnames(Y)
     if (trait_names[1] != species.id) {
     Y <- data.frame(Y[, species_col], Y[, -species_col, drop = FALSE])
     }
     colnames(Y) <- c("species", trait_names[-species_col])
     if (!intraspecific) {
     new_Y <- as.matrix(Y[, 2:ncol(Y), drop = FALSE])
     rownames(new_Y) <- as.character(Y[, 1])
     Y <- new_Y
     }
    } else if (class(Y) == "matrix") {
     if (is.null(rownames(Y))) {
     rownames(Y) <- tree$tip.label
     warning("Row names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     nms <- name.check(phy = tree, data.names = rownames(Y))
     if (class(nms) == "list") {
     rownames(Y) <- tree$tip.label
     warning("Row names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     if (is.null(colnames(Y)))
     colnames(Y) <- paste("V", 1:ncol(Y), sep = "")
    } else if (class(Y) == "numeric") {
     if (is.null(names(Y))) {
     warning("No species names provided for Y. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(Y) <- tree$tip.label
     }
     else {
     nms <- name.check(phy = tree, data.names = names(Y))
     if (class(nms) == "list") {
     warning("Names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(Y) <- tree$tip.label
     }
     Y <- matrix(data = Y, nrow = nspecies, ncol = 1, dimnames = list(names(Y),
     "V1"))
     }
    }
     --- R stacktrace ---
    where 1: evo.model(tree = rand.data$tree, Y = rand.data$trait_data, method = "Pairwise ML")
    
     --- value of length: 2 type: logical ---
    [1] FALSE FALSE
     --- function from context ---
    function (tree, Y, fixed.effects = NA, species.groups, trait.groups,
     model = "BM", diag.phylocov = FALSE, method = "Pairwise REML",
     force.zero.phylocov = character(), species.id = "species",
     max.combn = 10000, painted.edges, ret.level = 2, plot.LL.surface = FALSE,
     par.init.iters = 50, fixed.par = numeric(), multirate = FALSE,
     subset = TRUE, bounds)
    {
     tree <- reorder(tree, "postorder")
     if (length(tree$edge.length) > nrow(tree$edge))
     tree$edge.length <- tree$edge.length[1:nrow(tree$edge)]
     if (model == "OU")
     model <- "OUfixedRoot"
     if (missing(species.groups)) {
     nrates.species <- 1
     species.group.levels <- NA
     }
     else {
     species.groups <- species.groups[match(tree$tip.label,
     names(species.groups))]
     if (missing(painted.edges)) {
     painted.edges <- paint.edges(tree, species.groups)
     }
     nrates.species <- nlevels(species.groups)
     species.group.levels <- levels(species.groups)
     }
     if (missing(trait.groups)) {
     nrates.traits <- 1
     trait.group.levels <- NA
     }
     else {
     nrates.traits <- nlevels(trait.groups)
     trait.group.levels <- levels(trait.groups)
     }
     if (multirate == FALSE & nrates.traits > 1) {
     multirate <- TRUE
     warning("Setting multirate to TRUE because trait.groups was specified.",
     immediate. = TRUE)
     }
     if (method == "Pairwise REML" | method == "Full REML")
     REML <- 1
     else REML <- 0
     evo.model.args <- as.list(environment())
     evo.model.args$REML <- evo.model.args$trait.group.levels <- evo.model.args$nrates.traits <- evo.model.args$species.group.levels <- evo.model.args$nrates.species <- NULL
     nspecies <- length(tree$tip.label)
     nedge <- nrow(tree$edge)
     intraspecific <- FALSE
     if (class(fixed.effects) == "logical") {
     X <- numeric(0)
     X1 <- matrix(1, nspecies)
     rownames(X1) <- tree$tip.label
     }
     else {
     X <- fixed.effects
     if (class(X) == "data.frame") {
     species_col <- match(species.id, colnames(X))
     if (is.na(species_col))
     stop("Name of species.id much match a column in X (if supplying a data frame).")
     trait_names <- colnames(X)
     if (trait_names[1] != species.id) {
     X <- data.frame(X[, species_col], X[, -species_col,
     drop = FALSE])
     }
     colnames(X) <- c("species", trait_names[-species_col])
     if (!intraspecific) {
     new_X <- as.matrix(X[, 2:ncol(X), drop = FALSE])
     rownames(new_X) <- as.character(X[, 1])
     X <- new_X
     }
     }
     else if (class(X) == "matrix") {
     if (is.null(rownames(X))) {
     rownames(X) <- tree$tip.label
     warning("Row names of X do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     nms <- name.check(phy = tree, data.names = rownames(X))
     if (class(nms) == "list") {
     rownames(X) <- tree$tip.label
     warning("Row names of X do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     if (is.null(colnames(X)))
     colnames(X) <- paste("X", 1:ncol(X), sep = "")
     }
     else if (class(X) == "numeric") {
     if (is.null(names(X))) {
     warning("No species names provided for X. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(X) <- tree$tip.label
     }
     else {
     nms <- name.check(phy = tree, data.names = names(X))
     if (class(nms) == "list") {
     warning("Names of X do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(X) <- tree$tip.label
     }
     X <- matrix(data = X, nrow = nspecies, ncol = 1,
     dimnames = list(names(X), "X1"))
     }
     }
     X <- X[tree$tip.label, , drop = FALSE]
     evo.model.args$fixed.effects <- X
     X1 <- cbind(1, X)
     }
     if (class(Y) == "data.frame") {
     species_col <- match(species.id, colnames(Y))
     if (is.na(species_col))
     stop("Name of species.id much match a column in Y (if supplying a data frame).")
     trait_names <- colnames(Y)
     if (trait_names[1] != species.id) {
     Y <- data.frame(Y[, species_col], Y[, -species_col,
     drop = FALSE])
     }
     colnames(Y) <- c("species", trait_names[-species_col])
     if (!intraspecific) {
     new_Y <- as.matrix(Y[, 2:ncol(Y), drop = FALSE])
     rownames(new_Y) <- as.character(Y[, 1])
     Y <- new_Y
     }
     }
     else if (class(Y) == "matrix") {
     if (is.null(rownames(Y))) {
     rownames(Y) <- tree$tip.label
     warning("Row names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     nms <- name.check(phy = tree, data.names = rownames(Y))
     if (class(nms) == "list") {
     rownames(Y) <- tree$tip.label
     warning("Row names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     if (is.null(colnames(Y)))
     colnames(Y) <- paste("V", 1:ncol(Y), sep = "")
     }
     else if (class(Y) == "numeric") {
     if (is.null(names(Y))) {
     warning("No species names provided for Y. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(Y) <- tree$tip.label
     }
     else {
     nms <- name.check(phy = tree, data.names = names(Y))
     if (class(nms) == "list") {
     warning("Names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(Y) <- tree$tip.label
     }
     Y <- matrix(data = Y, nrow = nspecies, ncol = 1,
     dimnames = list(names(Y), "V1"))
     }
     }
     if (all(complete.cases(Y)))
     missing_data <- FALSE
     if (ncol(Y) == 1)
     if (method == "Pairwise ML")
     method <- "Full ML"
     else if (method == "Pairwise REML")
     method <- "Full REML"
     if (!intraspecific & !missing_data) {
     Y <- Y[tree$tip.label, , drop = FALSE]
     evo.model.args$Y <- Y
     ndims <- ncol(Y)
     diagndims <- diag(ndims)
     diag2 <- diag(2)
     if (method == "Full ML" | method == "Full REML") {
     MC <- FALSE
     }
     else {
     ncombn <- round(exp(lfactorial(ndims) - (log(2) +
     lfactorial(ndims - 2))))
     if (ncombn <= max.combn) {
     MC <- FALSE
     subsets <- matrix(0, 2, ncombn)
     counter <- 1
     for (r in 1:(ndims - 1)) {
     for (s in (r + 1):ndims) {
     subsets[, counter] <- c(r, s)
     counter <- counter + 1
     }
     }
     }
     else {
     MC <- TRUE
     subsets <- apply(matrix(0, max.combn, 2), 1,
     function(X) sample(x = ndims, size = 2, replace = FALSE))
     }
     }
     if (nrates.species == 1) {
     simple_BM <- function(temp_tree, ditree, ret_pars = FALSE) {
     pY$edge_len <- ditree$edge.length
     a_results <- do.call(multipic, pY)
     suminvV_and_logdV <- c(a_results[[2]], a_results[[3]])
     a <- a_results[[1]]
     if (length(X) > 0) {
     pX$edge_len <- ditree$edge.length
     x_results <- do.call(multipic, pX)
     x <- x_results[[1]]
     XANC <- x_results$root[1, ]
     XX <- crossprod(cbind(0, x)) + tcrossprod(c(1,
     XANC)) * suminvV_and_logdV[1]
     betas <- solve(crossprod(x), crossprod(x, a))
     C <- crossprod(a - x %*% betas)/(nspecies -
     REML)
     }
     else C <- crossprod(a)/(nspecies - REML)
     colnames(C) <- rownames(C) <- colnames(Y)
     if (ret_pars)
     if (length(X > 0)) {
     XY <- crossprod(cbind(0, x), a) + tcrossprod(c(1,
     XANC), a_results$root[1, ]) * suminvV_and_logdV[1]
     BETAS <- solve(XX, XY)
     predicted <- X1 %*% BETAS
     multipic.results <- list(X.multipic = x_results,
     Y.multipic = a_results)
     }
     else {
     predicted <- matrix(1, nspecies) %*% a_results$root[1,
     ]
     multipic.results <- list(Y.multipic = a_results)
     }
     if (diag.phylocov)
     C[upper.tri(C)] <- C[lower.tri(C)] <- 0
     diagC <- diag(C)
     if (nrates.traits > 1) {
     temp_C <- C
     for (i in 1:nrates.traits) {
     group_i <- which(trait.groups == trait.group.levels[i])
     diag(temp_C)[group_i] <- mean(diagC[group_i]) *
     if (subset)
     1
     else length(group_i)
     diagC[group_i] <- diag(temp_C)[group_i]
     }
     temp_C <- matrix(nearPD(x = temp_C, corr = FALSE,
     keepDiag = TRUE)$mat, ncol(C), ncol(C))
     C <- temp_C
     }
     else if (multirate) {
     diag(C) <- mean(diagC) * if (subset)
     1
     else ndims
     C <- matrix(nearPD(x = C, corr = FALSE, keepDiag = TRUE)$mat,
     ncol(C), ncol(C))
     diagC <- diag(C)
     }
     if (length(force.zero.phylocov) > 0) {
     other <- (1:ncol(C))[-match(force.zero.phylocov,
     colnames(C))]
     C[force.zero.phylocov, other] <- C[other, force.zero.phylocov] <- 0
     }
     if (method == "Full REML" | method == "Full ML") {
     const <- (nspecies * ndims - ndims * REML) *
     (log(2 * pi))
     ndims_logd <- ndims * suminvV_and_logdV[2]
     ndims_invV <- ndims * log(suminvV_and_logdV[1])
     log_Csub <- determinant(C, logarithm = TRUE)$modulus[[1]]
     if (length(X) > 0) {
     YY <- sum(tcrossprod(a - x %*% betas, backsolve(chol(C),
     diagndims, transpose = TRUE))^2)
     if (class(YY) == "try-error")
     YY <- sum(diag(solve(C) %*% crossprod(a -
     x %*% betas)))
     logdetXWX <- determinant(XX)$modulus[[1]] *
     ncol(C) - determinant(C)$modulus[[1]] *
     (ncol(X) + 1)
     }
     else {
     YY <- try(sum(tcrossprod(a, backsolve(chol(C),
     diagndims, transpose = TRUE))^2), silent = TRUE)
     if (class(YY) == "try-error")
     YY <- sum(diag(solve(C) %*% crossprod(a)))
     logdetXWX <- (-log_Csub + ndims_invV)
     }
     logL <- -(const + (nspecies * log_Csub + ndims_logd) +
     REML * logdetXWX + YY)/2
     }
     else if (method == "Pairwise REML" | method ==
     "Pairwise ML") {
     const <- (nspecies * 2 - 2 * REML) * (log(2 *
     pi))
     ndims_logd <- 2 * suminvV_and_logdV[2]
     ndims_invV <- 2 * log(suminvV_and_logdV[1])
     logL <- 0
     if (MC)
     subsets <- apply(matrix(0, max.combn, 2),
     1, function(X) sample(x = ndims, size = 2,
     replace = FALSE))
     for (i in 1:ncol(subsets)) {
     r <- subsets[1, i]
     s <- subsets[2, i]
     Csub <- C[c(r, s), c(r, s)]
     log_Csub <- log(Csub[1, 1] * Csub[2, 2] -
     Csub[1, 2]^2)
     if (length(X) > 0) {
     try(YY <- sum(tcrossprod(a[, c(r, s)] -
     x %*% betas[, c(r, s)], backsolve(chol(Csub),
     diag2, transpose = TRUE))^2), silent = TRUE)
     if (class(YY) == "try-error")
     YY <- sum(diag(solve(Csub) %*% crossprod(a[,
     c(r, s)] - x %*% betas[, c(r, s)])))
     logdetXWX <- determinant(XX)$modulus[[1]] *
     2 - log_Csub * (ncol(X) + 1)
     }
     else {
     try(YY <- sum(tcrossprod(a[, c(r, s)],
     backsolve(chol(Csub), diag2, transpose = TRUE))^2),
     silent = TRUE)
     if (class(YY) == "try-error")
     YY <- sum(diag(solve(Csub) %*% crossprod(a[,
     c(r, s)])))
     logdetXWX <- (-log_Csub + ndims_invV)
     }
     logL <- logL - (const + (nspecies * log_Csub +
     ndims_logd) + REML * logdetXWX + YY)
     }
     logL <- logL/i * ncombn/2
     }
     if (ret_pars)
     return(list(phylocov = C, logL = logL, predicted = predicted,
     transf_tree = temp_tree, multipic.results = multipic.results))
     else return(logL)
     }
     }
     else {
     ind <- 0L:1L
     tempY <- matrix(0, nspecies * 2, 1)
     tempY_span <- 1:(nspecies * 2)
     simple_BM <- function(temp_tree, ditree, ret_pars = FALSE) {
     pY$edge_len <- ditree$edge.length
     a_results <- do.call(multipic, pY)
     suminvV_and_logdV <- c(a_results[[2]], a_results[[3]])
     a <- a_results[[1]]
     Y_anc <- a_results$root[1, ]
     if (length(X) > 0) {
     pX$edge_len <- ditree$edge.length
     x_results <- do.call(multipic, pX)
     x <- x_results[[1]]
     XANC <- x_results$root[1, ]
     YANC <- a_results$root[1, ]
     XX <- crossprod(cbind(0, x)) + tcrossprod(c(1,
     XANC)) * suminvV_and_logdV[1]
     betas <- solve(crossprod(x), crossprod(x, a))
     XY <- crossprod(cbind(0, x), a) + tcrossprod(c(1,
     XANC), YANC) * suminvV_and_logdV[1]
     BETAS <- solve(XX, XY)
     C <- crossprod(a - x %*% betas)/(nspecies -
     REML)
     anc <- X1 %*% BETAS
     }
     else {
     anc <- matrix(1, nspecies) %*% a_results$root[1,
     ]
     C <- crossprod(a)/(nspecies - REML)
     }
     colnames(C) <- rownames(C) <- colnames(Y)
     if (ret_pars)
     if (length(X > 0)) {
     predicted <- anc
     multipic.results <- list(X.multipic = x_results,
     Y.multipic = a_results)
     }
     else {
     predicted <- anc
     multipic.results <- list(Y.multipic = a_results)
     }
     resid <- fast_transform(vcv(temp_tree)) %*% (Y -
     anc)
     rownames(resid) <- rownames(Y)
     species.subsets <- Clist <- diagClist <- vector("list",
     nrates.traits)
     for (j in 1:nrates.species) {
     species.subsets[[j]] <- which(!is.na(match(species.groups,
     species.group.levels[j])))
     Clist[[j]] <- t(resid[species.subsets[[j]],
     ]) %*% resid[species.subsets[[j]], ]/(length(species.subsets[[j]]) -
     REML)
     if (diag.phylocov)
     Clist[[j]][upper.tri(Clist[[j]])] <- Clist[[j]][lower.tri(Clist[[j]])] <- 0
     diagClist[[j]] <- diag(Clist[[j]])
     if (nrates.traits > 1) {
     temp_C <- Clist[[j]]
     for (i in 1:nrates.traits) {
     group_i <- which(trait.groups == trait.group.levels[i])
     diag(temp_C)[group_i] <- mean(diagClist[[j]][group_i]) *
     if (subset)
     1
     else length(group_i)
     }
     temp_C <- matrix(nearPD(x = temp_C, corr = FALSE,
     keepDiag = TRUE)$mat, ncol(Clist[[j]]),
     ncol(Clist[[j]]))
     Clist[[j]] <- temp_C
     diagClist[[j]] <- diag(Clist[[j]])
     }
     else if (multirate) {
     diag(Clist[[j]]) <- mean(diagClist[[j]]) *
     if (subset)
     1
     else length(diagClist[[j]])
     Clist[[j]] <- matrix(nearPD(x = Clist[[j]],
     corr = FALSE, keepDiag = TRUE)$mat, ncol(Clist[[j]]),
     ncol(Clist[[j]]))
     diagClist[[j]] <- diag(Clist[[j]])
     }
     if (length(force.zero.phylocov) > 0) {
     other <- which(is.na(match(colnames(Clist[[j]]),
     force.zero.phylocov)))
     Clist[[j]][force.zero.phylocov, other] <- Clist[[j]][other,
     force.zero.phylocov] <- 0
     }
     if (is.null(temp_tree$root.edge))
     temp_tree$edge.length <- c(temp_tree$edge.length[1:nedge],
     0)
     else temp_tree$edge.length <- c(temp_tree$edge.length[1:nedge],
     temp_tree$root.edge)
     }
     if (method == "Full REML" | method == "Full ML") {
     if (length(X) > 0) {
     ret <- Rinv5(nspecies = nspecies, nedge = nedge,
     edge.length = temp_tree$edge.length, anc = tree$edge[,
     1], des = tree$edge[, 2], L = X1, R = Y,
     painted.edges = painted.edges, phylocov_list = Clist)
     }
     else {
     ret <- Rinv3(nspecies = nspecies, nedge = nedge,
     edge.length = temp_tree$edge.length, anc = tree$edge[,
     1], des = tree$edge[, 2], R = Y, painted.edges = painted.edges,
     phylocov_list = Clist)
     }
     if (method == "Full REML")
     logL <- ret$ReLogLik
     else logL <- ret$LogLik
     }
     else if (method == "Pairwise REML" | method ==
     "Pairwise ML") {
     logL <- 0
     if (MC)
     subsets <- apply(matrix(0, max.combn, 2),
     1, function(X) sample(x = ndims, size = 2,
     replace = FALSE))
     anc_nodes <- temp_tree$edge[, 1] - 1
     des_nodes <- temp_tree$edge[, 2] - 1
     len_vec <- temp_tree$edge.length
     for (i in 1:ncol(subsets)) {
     r <- subsets[1, i]
     s <- subsets[2, i]
     tempC <- lapply(Clist, function(X) X[c(r,
     s), c(r, s)])
     tempY[tempY_span] <- Y[, c(r, s)]
     tempCmat <- matrix(0, 2 * nrates.species,
     2)
     for (zz in 1:length(tempC)) {
     tempCmat[1:2 + 2 * (zz - 1), ] <- tempC[[zz]]
     }
     if (length(X) > 0) {
     ret <- Rinv6(nspecies = nspecies, nedge = nedge,
     ngroups = nrates.species, ind = ind,
     len_vec = len_vec, anc = anc_nodes, des = des_nodes,
     L = X1, R = tempY, painted_edges = painted.edges,
     phylocovs = tempCmat, REML = REML)
     }
     else {
     ret <- Rinv4(nspecies = nspecies, nedge = nedge,
     ngroups = nrates.species, ind = ind,
     len_vec = len_vec, anc = anc_nodes, des = des_nodes,
     R = tempY, painted_edges = painted.edges,
     phylocovs = tempCmat, REML = REML)
     }
     logL <- logL + ret
     }
     logL <- logL/i * ncombn
     }
     if (ret_pars) {
     theta <- matrix(0, ncol(X1), ndims)
     for (i in 1:ndims) theta[, i] <- theta_Rinv6(nspecies = nspecies,
     nedge = nedge, ngroups = nrates.species,
     ind = 0L, len_vec = temp_tree$edge.length,
     anc = tree$edge[, 1] - 1, des = tree$edge[,
     2] - 1, L = X1, R = Y[, i, drop = FALSE],
     painted_edges = painted.edges, phylocovs = matrix(sapply(Clist,
     function(X) X[i, i])), REML = REML)
     predicted <- X1 %*% theta
     }
     if (ret_pars)
     return(list(phylocov = Clist, logL = logL,
     predicted = predicted, transf_tree = temp_tree))
     else return(logL)
     }
     }
     }
     else {
     stop("Missing data and intraspecific observations are not yet supported.")
     }
     if (!is.binary.tree(tree)) {
     binary <- TRUE
     ditree <- multi2di(tree, random = FALSE)
     }
     else {
     binary <- FALSE
     ditree <- tree
     }
     pY <- prep_multipic(Y, tree)
     if (length(X) > 0)
     pX <- prep_multipic(X, tree)
     if (model != "BM") {
     if (model != "BM")
     if ((model == "OUfixedRoot" | model == "OUrandomRoot") &
     !is.ultrametric(tree) & nrates.species > 1)
     stop("Non-ultrametric trees not currently supported for OU model with multiple species groups.")
     dis <- pruningwise.distFromRoot(reorder(tree, "pruningwise"))[1:nspecies]
     Tmax <- mean(dis)
     bounds.default <- matrix(c(1e-07/Tmax, max(50/Tmax, 5),
     1e-07/Tmax, max(50/Tmax, 5), 1e-07, 1, 1e-06, 1,
     1e-05, 3, min(-4, -3/Tmax), 0), ncol = 2, byrow = TRUE)
     rownames(bounds.default) <- c("alpha", "alpha", "lambda",
     "kappa", "delta", "rate")
     colnames(bounds.default) <- c("min", "max")
     starting.values.default <- c(0.5/Tmax, 0.5/Tmax, 0.5,
     0.5, 0.5, -1/Tmax)
     if (!missing(bounds)) {
     bounds.default[, 1] <- bounds[1]
     bounds.default[, 2] <- bounds[2]
     if (any(starting.values.default < bounds[1]) | any(starting.values.default >
     bounds[1])) {
     starting.values.default[1:length(starting.values.default)] <- mean(bounds[1:2])
     }
     }
     names(starting.values.default) <- c("OUrandomRoot", "OUfixedRoot",
     "lambda", "kappa", "delta", "EB")
     model_i <- match(model, names(starting.values.default))
     bounds.default <- bounds.default[model_i, , drop = FALSE]
     starting.values.default <- starting.values.default[model_i]
     names(starting.values.default) <- rownames(bounds.default)
     if (length(fixed.par) > 0) {
     plot.LL.surface <- FALSE
     bounds.default[1:2] <- rep(fixed.par, 2)
     starting.values.default[1] <- fixed.par
     }
     if (model != "EB" & model != "BM") {
     starting.values.default <- log(starting.values.default)
     bounds.default <- log(bounds.default)
     }
     if (model[1] == "BM")
     starting.values.default <- bounds.default <- NULL
     geiger_args <- list(x = tree, model = if (model == "OUfixedRoot" |
     model == "OUrandomRoot") "OU" else model, par = 5)
     names(geiger_args)[3] <- names(starting.values.default)
     if (model == "EB")
     names(geiger_args)[3] <- "a"
     if (model == "OUrandomRoot" | model == "OUfixedRoot" |
     model == "OU") {
     temp_args <- list(phy = tree, model = model, parameters = list(pars = 5))
     names(temp_args$parameters) <- names(starting.values.default)
     }
     model_LL <- function(par, ret_pars = FALSE) {
     if (model != "EB")
     par <- exp(par)
     if (model == "OUrandomRoot" | model == "OUfixedRoot" |
     model == "OU") {
     temp_args$parameters <- list(alpha = par)
     temp_tree <- do.call(transf.branch.lengths, temp_args)$tree
     }
     else {
     geiger_args[[3]] <- par
     temp_tree <- do.call(rescale, geiger_args)
     }
     temp_tree <- reorder(temp_tree, "postorder")
     if (!binary)
     temp_ditree <- multi2di(temp_tree, random = FALSE)
     else temp_ditree <- temp_tree
     ret <- try(simple_BM(temp_tree = temp_tree, ditree = temp_ditree,
     ret_pars = ret_pars), silent = TRUE)
     if (class(ret) == "try-error")
     return(-.Machine$double.xmax)
     else return(ret)
     }
     if (MC)
     try_length <- max(50, par.init.iters)
     else try_length <- par.init.iters
     if (length(fixed.par) > 0)
     try_length <- 1
     if (model != "EB")
     try_theta <- log(seq(exp(min(bounds.default)), exp(max(bounds.default)),
     length = try_length))
     else try_theta <- (seq((min(bounds.default)), (max(bounds.default)),
     length = try_length))
     try_LL <- double(length(try_theta))
     for (i in 1:length(try_theta)) {
     try_LL[i] <- model_LL(try_theta[i])
     }
     if (plot.LL.surface & !MC)
     try({
     plot(if (model != "EB")
     exp(try_theta)
     else try_theta, try_LL)
     }, silent = TRUE)
     starting.values.default <- try_theta[which.max(try_LL)]
     if (MC & length(fixed.par) == 0) {
     if (model != "EB")
     try_theta <- exp(try_theta)
     if (model != "EB")
     bounds.default <- exp(bounds.default)
     best <- double(round(length(try_theta)/2))
     for (i in 1:length(best)) {
     best[i] <- BIC(lm(try_LL ~ poly(try_theta, degree = i)))
     }
     degree <- which.min(best)
     LL_surface <- lm(try_LL ~ poly(try_theta, degree = degree))
     theta_start <- try_theta[which(predict(LL_surface) ==
     max(predict(LL_surface)))][1]
     o <- optim(par = theta_start, fn = function(X) predict.lm(LL_surface,
     newdata = data.frame(try_theta = X)), control = list(fnscale = -1),
     method = "L-BFGS-B", lower = bounds.default[1],
     upper = bounds.default[2])
     if (plot.LL.surface)
     try({
     plot(try_theta, try_LL)
     points(try_theta, predict.lm(LL_surface), col = "blue")
     points(o$par, o$value, col = "red", pch = 16)
     }, silent = TRUE)
     ret <- model_LL(par = if (model != "EB")
     log(o$par)
     else o$par, ret_pars = ret.level > 1)
     ret$model.par[[rownames(bounds.default)]] <- o$par
     if (ret.level > 1 & (model == "OUfixedRoot" | model ==
     "OUrandomRoot" | model == "OU")) {
     if (nrates.species > 1) {
     for (i in 1:nrates.species) ret$phylocov[[i]] <- 2 *
     ret$phylocov[[i]] * ret$alpha
     }
     else ret$phylocov <- 2 * ret$phylocov * ret$alpha
     }
     }
     else {
     if (length(fixed.par) > 0) {
     o <- list(par = starting.values.default, value = model_LL(starting.values.default))
     }
     else o <- optim(starting.values.default, model_LL,
     control = list(fnscale = -1), method = "L-BFGS-B",
     lower = bounds.default[1], upper = bounds.default[2])
     ret <- model_LL(par = o$par, ret_pars = ret.level >
     1)
     if (model != "EB")
     o$par <- exp(o$par)
     if (ret.level > 1) {
     ret$model.par[[rownames(bounds.default)]] <- o$par
     if (ret.level > 1 & (model == "OUfixedRoot" |
     model == "OUrandomRoot" | model == "OU")) {
     if (nrates.species > 1) {
     for (i in 1:nrates.species) ret$phylocov[[i]] <- 2 *
     ret$phylocov[[i]] * ret$model.par[[1]]
     }
     else ret$phylocov <- 2 * ret$phylocov * ret$model.par[[1]]
     }
     }
     }
     }
     else ret <- simple_BM(temp_tree = tree, ditree = ditree,
     ret_pars = ret.level > 1)
     if (ret.level <= 1)
     return(ret)
     if (ret.level > 1) {
     ret$method <- method
     sigma2.mult <- vector("list", nrates.species)
     if (nrates.species > 1) {
     for (i in 1:nrates.species) {
     if (nrates.traits < 2)
     sigma2.mult[[i]] <- mean(diag(ret$phylocov[[i]]))
     else {
     for (j in 1:nrates.traits) {
     group_j <- which(trait.groups == trait.group.levels[j])
     sigma2.mult[[i]] <- c(sigma2.mult[[i]], mean(diag(ret$phylocov[[i]])[group_j]) *
     if (subset) 1 else length(group_j))
     }
     names(sigma2.mult[[i]]) <- levels(trait.groups)
     }
     }
     names(sigma2.mult) <- levels(species.groups)
     }
     else {
     if (nrates.traits < 2)
     sigma2.mult <- mean(diag(ret$phylocov))
     else {
     sigma2.mult <- double()
     for (j in 1:nrates.traits) {
     group_j <- which(trait.groups == trait.group.levels[j])
     sigma2.mult <- c(sigma2.mult, mean(diag(ret$phylocov)[group_j]) *
     if (subset) 1 else length(group_j))
     }
     names(sigma2.mult) <- levels(trait.groups)
     }
     }
     if (ret.level < 3)
     ret$phylocov <- NULL
     ret <- c(ret, list(sigma2.mult = sigma2.mult, evo.model.args = evo.model.args))
     }
     class(ret) <- "evo.model"
     return(ret)
    }
    <bytecode: 0xcad4ef8>
    <environment: namespace:phylocurve>
     --- function search by body ---
    Function evo.model in namespace phylocurve has this body.
     ----------- END OF FAILURE REPORT --------------
    Fatal error: the condition has length > 1
Flavor: r-devel-linux-x86_64-debian-clang

Version: 2.0.9
Check: examples
Result: ERROR
    Running examples in ‘phylocurve-Ex.R’ failed
    The error most likely occurred in:
    
    > base::assign(".ptime", proc.time(), pos = "CheckExEnv")
    > ### Name: K.mult
    > ### Title: Test phylogenetic signal (Kmult) using phylogenetic simulation
    > ### Aliases: K.mult
    >
    > ### ** Examples
    >
    > rand.data <- sim.traits()
    > null.model <- evo.model(tree = rand.data$tree,
    + Y = rand.data$trait_data,method = "Pairwise ML")
     ----------- FAILURE REPORT --------------
     --- failure: the condition has length > 1 ---
     --- srcref ---
    :
     --- package (from environment) ---
    phylocurve
     --- call from context ---
    evo.model(tree = rand.data$tree, Y = rand.data$trait_data, method = "Pairwise ML")
     --- call from argument ---
    if (class(Y) == "data.frame") {
     species_col <- match(species.id, colnames(Y))
     if (is.na(species_col))
     stop("Name of species.id much match a column in Y (if supplying a data frame).")
     trait_names <- colnames(Y)
     if (trait_names[1] != species.id) {
     Y <- data.frame(Y[, species_col], Y[, -species_col, drop = FALSE])
     }
     colnames(Y) <- c("species", trait_names[-species_col])
     if (!intraspecific) {
     new_Y <- as.matrix(Y[, 2:ncol(Y), drop = FALSE])
     rownames(new_Y) <- as.character(Y[, 1])
     Y <- new_Y
     }
    } else if (class(Y) == "matrix") {
     if (is.null(rownames(Y))) {
     rownames(Y) <- tree$tip.label
     warning("Row names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     nms <- name.check(phy = tree, data.names = rownames(Y))
     if (class(nms) == "list") {
     rownames(Y) <- tree$tip.label
     warning("Row names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     if (is.null(colnames(Y)))
     colnames(Y) <- paste("V", 1:ncol(Y), sep = "")
    } else if (class(Y) == "numeric") {
     if (is.null(names(Y))) {
     warning("No species names provided for Y. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(Y) <- tree$tip.label
     }
     else {
     nms <- name.check(phy = tree, data.names = names(Y))
     if (class(nms) == "list") {
     warning("Names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(Y) <- tree$tip.label
     }
     Y <- matrix(data = Y, nrow = nspecies, ncol = 1, dimnames = list(names(Y),
     "V1"))
     }
    }
     --- R stacktrace ---
    where 1: evo.model(tree = rand.data$tree, Y = rand.data$trait_data, method = "Pairwise ML")
    
     --- value of length: 2 type: logical ---
    [1] FALSE FALSE
     --- function from context ---
    function (tree, Y, fixed.effects = NA, species.groups, trait.groups,
     model = "BM", diag.phylocov = FALSE, method = "Pairwise REML",
     force.zero.phylocov = character(), species.id = "species",
     max.combn = 10000, painted.edges, ret.level = 2, plot.LL.surface = FALSE,
     par.init.iters = 50, fixed.par = numeric(), multirate = FALSE,
     subset = TRUE, bounds)
    {
     tree <- reorder(tree, "postorder")
     if (length(tree$edge.length) > nrow(tree$edge))
     tree$edge.length <- tree$edge.length[1:nrow(tree$edge)]
     if (model == "OU")
     model <- "OUfixedRoot"
     if (missing(species.groups)) {
     nrates.species <- 1
     species.group.levels <- NA
     }
     else {
     species.groups <- species.groups[match(tree$tip.label,
     names(species.groups))]
     if (missing(painted.edges)) {
     painted.edges <- paint.edges(tree, species.groups)
     }
     nrates.species <- nlevels(species.groups)
     species.group.levels <- levels(species.groups)
     }
     if (missing(trait.groups)) {
     nrates.traits <- 1
     trait.group.levels <- NA
     }
     else {
     nrates.traits <- nlevels(trait.groups)
     trait.group.levels <- levels(trait.groups)
     }
     if (multirate == FALSE & nrates.traits > 1) {
     multirate <- TRUE
     warning("Setting multirate to TRUE because trait.groups was specified.",
     immediate. = TRUE)
     }
     if (method == "Pairwise REML" | method == "Full REML")
     REML <- 1
     else REML <- 0
     evo.model.args <- as.list(environment())
     evo.model.args$REML <- evo.model.args$trait.group.levels <- evo.model.args$nrates.traits <- evo.model.args$species.group.levels <- evo.model.args$nrates.species <- NULL
     nspecies <- length(tree$tip.label)
     nedge <- nrow(tree$edge)
     intraspecific <- FALSE
     if (class(fixed.effects) == "logical") {
     X <- numeric(0)
     X1 <- matrix(1, nspecies)
     rownames(X1) <- tree$tip.label
     }
     else {
     X <- fixed.effects
     if (class(X) == "data.frame") {
     species_col <- match(species.id, colnames(X))
     if (is.na(species_col))
     stop("Name of species.id much match a column in X (if supplying a data frame).")
     trait_names <- colnames(X)
     if (trait_names[1] != species.id) {
     X <- data.frame(X[, species_col], X[, -species_col,
     drop = FALSE])
     }
     colnames(X) <- c("species", trait_names[-species_col])
     if (!intraspecific) {
     new_X <- as.matrix(X[, 2:ncol(X), drop = FALSE])
     rownames(new_X) <- as.character(X[, 1])
     X <- new_X
     }
     }
     else if (class(X) == "matrix") {
     if (is.null(rownames(X))) {
     rownames(X) <- tree$tip.label
     warning("Row names of X do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     nms <- name.check(phy = tree, data.names = rownames(X))
     if (class(nms) == "list") {
     rownames(X) <- tree$tip.label
     warning("Row names of X do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     if (is.null(colnames(X)))
     colnames(X) <- paste("X", 1:ncol(X), sep = "")
     }
     else if (class(X) == "numeric") {
     if (is.null(names(X))) {
     warning("No species names provided for X. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(X) <- tree$tip.label
     }
     else {
     nms <- name.check(phy = tree, data.names = names(X))
     if (class(nms) == "list") {
     warning("Names of X do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(X) <- tree$tip.label
     }
     X <- matrix(data = X, nrow = nspecies, ncol = 1,
     dimnames = list(names(X), "X1"))
     }
     }
     X <- X[tree$tip.label, , drop = FALSE]
     evo.model.args$fixed.effects <- X
     X1 <- cbind(1, X)
     }
     if (class(Y) == "data.frame") {
     species_col <- match(species.id, colnames(Y))
     if (is.na(species_col))
     stop("Name of species.id much match a column in Y (if supplying a data frame).")
     trait_names <- colnames(Y)
     if (trait_names[1] != species.id) {
     Y <- data.frame(Y[, species_col], Y[, -species_col,
     drop = FALSE])
     }
     colnames(Y) <- c("species", trait_names[-species_col])
     if (!intraspecific) {
     new_Y <- as.matrix(Y[, 2:ncol(Y), drop = FALSE])
     rownames(new_Y) <- as.character(Y[, 1])
     Y <- new_Y
     }
     }
     else if (class(Y) == "matrix") {
     if (is.null(rownames(Y))) {
     rownames(Y) <- tree$tip.label
     warning("Row names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     nms <- name.check(phy = tree, data.names = rownames(Y))
     if (class(nms) == "list") {
     rownames(Y) <- tree$tip.label
     warning("Row names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     }
     if (is.null(colnames(Y)))
     colnames(Y) <- paste("V", 1:ncol(Y), sep = "")
     }
     else if (class(Y) == "numeric") {
     if (is.null(names(Y))) {
     warning("No species names provided for Y. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(Y) <- tree$tip.label
     }
     else {
     nms <- name.check(phy = tree, data.names = names(Y))
     if (class(nms) == "list") {
     warning("Names of Y do not match tip labels of tree. Assuming data are in order of tips.",
     immediate. = TRUE)
     names(Y) <- tree$tip.label
     }
     Y <- matrix(data = Y, nrow = nspecies, ncol = 1,
     dimnames = list(names(Y), "V1"))
     }
     }
     if (all(complete.cases(Y)))
     missing_data <- FALSE
     if (ncol(Y) == 1)
     if (method == "Pairwise ML")
     method <- "Full ML"
     else if (method == "Pairwise REML")
     method <- "Full REML"
     if (!intraspecific & !missing_data) {
     Y <- Y[tree$tip.label, , drop = FALSE]
     evo.model.args$Y <- Y
     ndims <- ncol(Y)
     diagndims <- diag(ndims)
     diag2 <- diag(2)
     if (method == "Full ML" | method == "Full REML") {
     MC <- FALSE
     }
     else {
     ncombn <- round(exp(lfactorial(ndims) - (log(2) +
     lfactorial(ndims - 2))))
     if (ncombn <= max.combn) {
     MC <- FALSE
     subsets <- matrix(0, 2, ncombn)
     counter <- 1
     for (r in 1:(ndims - 1)) {
     for (s in (r + 1):ndims) {
     subsets[, counter] <- c(r, s)
     counter <- counter + 1
     }
     }
     }
     else {
     MC <- TRUE
     subsets <- apply(matrix(0, max.combn, 2), 1,
     function(X) sample(x = ndims, size = 2, replace = FALSE))
     }
     }
     if (nrates.species == 1) {
     simple_BM <- function(temp_tree, ditree, ret_pars = FALSE) {
     pY$edge_len <- ditree$edge.length
     a_results <- do.call(multipic, pY)
     suminvV_and_logdV <- c(a_results[[2]], a_results[[3]])
     a <- a_results[[1]]
     if (length(X) > 0) {
     pX$edge_len <- ditree$edge.length
     x_results <- do.call(multipic, pX)
     x <- x_results[[1]]
     XANC <- x_results$root[1, ]
     XX <- crossprod(cbind(0, x)) + tcrossprod(c(1,
     XANC)) * suminvV_and_logdV[1]
     betas <- solve(crossprod(x), crossprod(x, a))
     C <- crossprod(a - x %*% betas)/(nspecies -
     REML)
     }
     else C <- crossprod(a)/(nspecies - REML)
     colnames(C) <- rownames(C) <- colnames(Y)
     if (ret_pars)
     if (length(X > 0)) {
     XY <- crossprod(cbind(0, x), a) + tcrossprod(c(1,
     XANC), a_results$root[1, ]) * suminvV_and_logdV[1]
     BETAS <- solve(XX, XY)
     predicted <- X1 %*% BETAS
     multipic.results <- list(X.multipic = x_results,
     Y.multipic = a_results)
     }
     else {
     predicted <- matrix(1, nspecies) %*% a_results$root[1,
     ]
     multipic.results <- list(Y.multipic = a_results)
     }
     if (diag.phylocov)
     C[upper.tri(C)] <- C[lower.tri(C)] <- 0
     diagC <- diag(C)
     if (nrates.traits > 1) {
     temp_C <- C
     for (i in 1:nrates.traits) {
     group_i <- which(trait.groups == trait.group.levels[i])
     diag(temp_C)[group_i] <- mean(diagC[group_i]) *
     if (subset)
     1
     else length(group_i)
     diagC[group_i] <- diag(temp_C)[group_i]
     }
     temp_C <- matrix(nearPD(x = temp_C, corr = FALSE,
     keepDiag = TRUE)$mat, ncol(C), ncol(C))
     C <- temp_C
     }
     else if (multirate) {
     diag(C) <- mean(diagC) * if (subset)
     1
     else ndims
     C <- matrix(nearPD(x = C, corr = FALSE, keepDiag = TRUE)$mat,
     ncol(C), ncol(C))
     diagC <- diag(C)
     }
     if (length(force.zero.phylocov) > 0) {
     other <- (1:ncol(C))[-match(force.zero.phylocov,
     colnames(C))]
     C[force.zero.phylocov, other] <- C[other, force.zero.phylocov] <- 0
     }
     if (method == "Full REML" | method == "Full ML") {
     const <- (nspecies * ndims - ndims * REML) *
     (log(2 * pi))
     ndims_logd <- ndims * suminvV_and_logdV[2]
     ndims_invV <- ndims * log(suminvV_and_logdV[1])
     log_Csub <- determinant(C, logarithm = TRUE)$modulus[[1]]
     if (length(X) > 0) {
     YY <- sum(tcrossprod(a - x %*% betas, backsolve(chol(C),
     diagndims, transpose = TRUE))^2)
     if (class(YY) == "try-error")
     YY <- sum(diag(solve(C) %*% crossprod(a -
     x %*% betas)))
     logdetXWX <- determinant(XX)$modulus[[1]] *
     ncol(C) - determinant(C)$modulus[[1]] *
     (ncol(X) + 1)
     }
     else {
     YY <- try(sum(tcrossprod(a, backsolve(chol(C),
     diagndims, transpose = TRUE))^2), silent = TRUE)
     if (class(YY) == "try-error")
     YY <- sum(diag(solve(C) %*% crossprod(a)))
     logdetXWX <- (-log_Csub + ndims_invV)
     }
     logL <- -(const + (nspecies * log_Csub + ndims_logd) +
     REML * logdetXWX + YY)/2
     }
     else if (method == "Pairwise REML" | method ==
     "Pairwise ML") {
     const <- (nspecies * 2 - 2 * REML) * (log(2 *
     pi))
     ndims_logd <- 2 * suminvV_and_logdV[2]
     ndims_invV <- 2 * log(suminvV_and_logdV[1])
     logL <- 0
     if (MC)
     subsets <- apply(matrix(0, max.combn, 2),
     1, function(X) sample(x = ndims, size = 2,
     replace = FALSE))
     for (i in 1:ncol(subsets)) {
     r <- subsets[1, i]
     s <- subsets[2, i]
     Csub <- C[c(r, s), c(r, s)]
     log_Csub <- log(Csub[1, 1] * Csub[2, 2] -
     Csub[1, 2]^2)
     if (length(X) > 0) {
     try(YY <- sum(tcrossprod(a[, c(r, s)] -
     x %*% betas[, c(r, s)], backsolve(chol(Csub),
     diag2, transpose = TRUE))^2), silent = TRUE)
     if (class(YY) == "try-error")
     YY <- sum(diag(solve(Csub) %*% crossprod(a[,
     c(r, s)] - x %*% betas[, c(r, s)])))
     logdetXWX <- determinant(XX)$modulus[[1]] *
     2 - log_Csub * (ncol(X) + 1)
     }
     else {
     try(YY <- sum(tcrossprod(a[, c(r, s)],
     backsolve(chol(Csub), diag2, transpose = TRUE))^2),
     silent = TRUE)
     if (class(YY) == "try-error")
     YY <- sum(diag(solve(Csub) %*% crossprod(a[,
     c(r, s)])))
     logdetXWX <- (-log_Csub + ndims_invV)
     }
     logL <- logL - (const + (nspecies * log_Csub +
     ndims_logd) + REML * logdetXWX + YY)
     }
     logL <- logL/i * ncombn/2
     }
     if (ret_pars)
     return(list(phylocov = C, logL = logL, predicted = predicted,
     transf_tree = temp_tree, multipic.results = multipic.results))
     else return(logL)
     }
     }
     else {
     ind <- 0L:1L
     tempY <- matrix(0, nspecies * 2, 1)
     tempY_span <- 1:(nspecies * 2)
     simple_BM <- function(temp_tree, ditree, ret_pars = FALSE) {
     pY$edge_len <- ditree$edge.length
     a_results <- do.call(multipic, pY)
     suminvV_and_logdV <- c(a_results[[2]], a_results[[3]])
     a <- a_results[[1]]
     Y_anc <- a_results$root[1, ]
     if (length(X) > 0) {
     pX$edge_len <- ditree$edge.length
     x_results <- do.call(multipic, pX)
     x <- x_results[[1]]
     XANC <- x_results$root[1, ]
     YANC <- a_results$root[1, ]
     XX <- crossprod(cbind(0, x)) + tcrossprod(c(1,
     XANC)) * suminvV_and_logdV[1]
     betas <- solve(crossprod(x), crossprod(x, a))
     XY <- crossprod(cbind(0, x), a) + tcrossprod(c(1,
     XANC), YANC) * suminvV_and_logdV[1]
     BETAS <- solve(XX, XY)
     C <- crossprod(a - x %*% betas)/(nspecies -
     REML)
     anc <- X1 %*% BETAS
     }
     else {
     anc <- matrix(1, nspecies) %*% a_results$root[1,
     ]
     C <- crossprod(a)/(nspecies - REML)
     }
     colnames(C) <- rownames(C) <- colnames(Y)
     if (ret_pars)
     if (length(X > 0)) {
     predicted <- anc
     multipic.results <- list(X.multipic = x_results,
     Y.multipic = a_results)
     }
     else {
     predicted <- anc
     multipic.results <- list(Y.multipic = a_results)
     }
     resid <- fast_transform(vcv(temp_tree)) %*% (Y -
     anc)
     rownames(resid) <- rownames(Y)
     species.subsets <- Clist <- diagClist <- vector("list",
     nrates.traits)
     for (j in 1:nrates.species) {
     species.subsets[[j]] <- which(!is.na(match(species.groups,
     species.group.levels[j])))
     Clist[[j]] <- t(resid[species.subsets[[j]],
     ]) %*% resid[species.subsets[[j]], ]/(length(species.subsets[[j]]) -
     REML)
     if (diag.phylocov)
     Clist[[j]][upper.tri(Clist[[j]])] <- Clist[[j]][lower.tri(Clist[[j]])] <- 0
     diagClist[[j]] <- diag(Clist[[j]])
     if (nrates.traits > 1) {
     temp_C <- Clist[[j]]
     for (i in 1:nrates.traits) {
     group_i <- which(trait.groups == trait.group.levels[i])
     diag(temp_C)[group_i] <- mean(diagClist[[j]][group_i]) *
     if (subset)
     1
     else length(group_i)
     }
     temp_C <- matrix(nearPD(x = temp_C, corr = FALSE,
     keepDiag = TRUE)$mat, ncol(Clist[[j]]),
     ncol(Clist[[j]]))
     Clist[[j]] <- temp_C
     diagClist[[j]] <- diag(Clist[[j]])
     }
     else if (multirate) {
     diag(Clist[[j]]) <- mean(diagClist[[j]]) *
     if (subset)
     1
     else length(diagClist[[j]])
     Clist[[j]] <- matrix(nearPD(x = Clist[[j]],
     corr = FALSE, keepDiag = TRUE)$mat, ncol(Clist[[j]]),
     ncol(Clist[[j]]))
     diagClist[[j]] <- diag(Clist[[j]])
     }
     if (length(force.zero.phylocov) > 0) {
     other <- which(is.na(match(colnames(Clist[[j]]),
     force.zero.phylocov)))
     Clist[[j]][force.zero.phylocov, other] <- Clist[[j]][other,
     force.zero.phylocov] <- 0
     }
     if (is.null(temp_tree$root.edge))
     temp_tree$edge.length <- c(temp_tree$edge.length[1:nedge],
     0)
     else temp_tree$edge.length <- c(temp_tree$edge.length[1:nedge],
     temp_tree$root.edge)
     }
     if (method == "Full REML" | method == "Full ML") {
     if (length(X) > 0) {
     ret <- Rinv5(nspecies = nspecies, nedge = nedge,
     edge.length = temp_tree$edge.length, anc = tree$edge[,
     1], des = tree$edge[, 2], L = X1, R = Y,
     painted.edges = painted.edges, phylocov_list = Clist)
     }
     else {
     ret <- Rinv3(nspecies = nspecies, nedge = nedge,
     edge.length = temp_tree$edge.length, anc = tree$edge[,
     1], des = tree$edge[, 2], R = Y, painted.edges = painted.edges,
     phylocov_list = Clist)
     }
     if (method == "Full REML")
     logL <- ret$ReLogLik
     else logL <- ret$LogLik
     }
     else if (method == "Pairwise REML" | method ==
     "Pairwise ML") {
     logL <- 0
     if (MC)
     subsets <- apply(matrix(0, max.combn, 2),
     1, function(X) sample(x = ndims, size = 2,
     replace = FALSE))
     anc_nodes <- temp_tree$edge[, 1] - 1
     des_nodes <- temp_tree$edge[, 2] - 1
     len_vec <- temp_tree$edge.length
     for (i in 1:ncol(subsets)) {
     r <- subsets[1, i]
     s <- subsets[2, i]
     tempC <- lapply(Clist, function(X) X[c(r,
     s), c(r, s)])
     tempY[tempY_span] <- Y[, c(r, s)]
     tempCmat <- matrix(0, 2 * nrates.species,
     2)
     for (zz in 1:length(tempC)) {
     tempCmat[1:2 + 2 * (zz - 1), ] <- tempC[[zz]]
     }
     if (length(X) > 0) {
     ret <- Rinv6(nspecies = nspecies, nedge = nedge,
     ngroups = nrates.species, ind = ind,
     len_vec = len_vec, anc = anc_nodes, des = des_nodes,
     L = X1, R = tempY, painted_edges = painted.edges,
     phylocovs = tempCmat, REML = REML)
     }
     else {
     ret <- Rinv4(nspecies = nspecies, nedge = nedge,
     ngroups = nrates.species, ind = ind,
     len_vec = len_vec, anc = anc_nodes, des = des_nodes,
     R = tempY, painted_edges = painted.edges,
     phylocovs = tempCmat, REML = REML)
     }
     logL <- logL + ret
     }
     logL <- logL/i * ncombn
     }
     if (ret_pars) {
     theta <- matrix(0, ncol(X1), ndims)
     for (i in 1:ndims) theta[, i] <- theta_Rinv6(nspecies = nspecies,
     nedge = nedge, ngroups = nrates.species,
     ind = 0L, len_vec = temp_tree$edge.length,
     anc = tree$edge[, 1] - 1, des = tree$edge[,
     2] - 1, L = X1, R = Y[, i, drop = FALSE],
     painted_edges = painted.edges, phylocovs = matrix(sapply(Clist,
     function(X) X[i, i])), REML = REML)
     predicted <- X1 %*% theta
     }
     if (ret_pars)
     return(list(phylocov = Clist, logL = logL,
     predicted = predicted, transf_tree = temp_tree))
     else return(logL)
     }
     }
     }
     else {
     stop("Missing data and intraspecific observations are not yet supported.")
     }
     if (!is.binary.tree(tree)) {
     binary <- TRUE
     ditree <- multi2di(tree, random = FALSE)
     }
     else {
     binary <- FALSE
     ditree <- tree
     }
     pY <- prep_multipic(Y, tree)
     if (length(X) > 0)
     pX <- prep_multipic(X, tree)
     if (model != "BM") {
     if (model != "BM")
     if ((model == "OUfixedRoot" | model == "OUrandomRoot") &
     !is.ultrametric(tree) & nrates.species > 1)
     stop("Non-ultrametric trees not currently supported for OU model with multiple species groups.")
     dis <- pruningwise.distFromRoot(reorder(tree, "pruningwise"))[1:nspecies]
     Tmax <- mean(dis)
     bounds.default <- matrix(c(1e-07/Tmax, max(50/Tmax, 5),
     1e-07/Tmax, max(50/Tmax, 5), 1e-07, 1, 1e-06, 1,
     1e-05, 3, min(-4, -3/Tmax), 0), ncol = 2, byrow = TRUE)
     rownames(bounds.default) <- c("alpha", "alpha", "lambda",
     "kappa", "delta", "rate")
     colnames(bounds.default) <- c("min", "max")
     starting.values.default <- c(0.5/Tmax, 0.5/Tmax, 0.5,
     0.5, 0.5, -1/Tmax)
     if (!missing(bounds)) {
     bounds.default[, 1] <- bounds[1]
     bounds.default[, 2] <- bounds[2]
     if (any(starting.values.default < bounds[1]) | any(starting.values.default >
     bounds[1])) {
     starting.values.default[1:length(starting.values.default)] <- mean(bounds[1:2])
     }
     }
     names(starting.values.default) <- c("OUrandomRoot", "OUfixedRoot",
     "lambda", "kappa", "delta", "EB")
     model_i <- match(model, names(starting.values.default))
     bounds.default <- bounds.default[model_i, , drop = FALSE]
     starting.values.default <- starting.values.default[model_i]
     names(starting.values.default) <- rownames(bounds.default)
     if (length(fixed.par) > 0) {
     plot.LL.surface <- FALSE
     bounds.default[1:2] <- rep(fixed.par, 2)
     starting.values.default[1] <- fixed.par
     }
     if (model != "EB" & model != "BM") {
     starting.values.default <- log(starting.values.default)
     bounds.default <- log(bounds.default)
     }
     if (model[1] == "BM")
     starting.values.default <- bounds.default <- NULL
     geiger_args <- list(x = tree, model = if (model == "OUfixedRoot" |
     model == "OUrandomRoot") "OU" else model, par = 5)
     names(geiger_args)[3] <- names(starting.values.default)
     if (model == "EB")
     names(geiger_args)[3] <- "a"
     if (model == "OUrandomRoot" | model == "OUfixedRoot" |
     model == "OU") {
     temp_args <- list(phy = tree, model = model, parameters = list(pars = 5))
     names(temp_args$parameters) <- names(starting.values.default)
     }
     model_LL <- function(par, ret_pars = FALSE) {
     if (model != "EB")
     par <- exp(par)
     if (model == "OUrandomRoot" | model == "OUfixedRoot" |
     model == "OU") {
     temp_args$parameters <- list(alpha = par)
     temp_tree <- do.call(transf.branch.lengths, temp_args)$tree
     }
     else {
     geiger_args[[3]] <- par
     temp_tree <- do.call(rescale, geiger_args)
     }
     temp_tree <- reorder(temp_tree, "postorder")
     if (!binary)
     temp_ditree <- multi2di(temp_tree, random = FALSE)
     else temp_ditree <- temp_tree
     ret <- try(simple_BM(temp_tree = temp_tree, ditree = temp_ditree,
     ret_pars = ret_pars), silent = TRUE)
     if (class(ret) == "try-error")
     return(-.Machine$double.xmax)
     else return(ret)
     }
     if (MC)
     try_length <- max(50, par.init.iters)
     else try_length <- par.init.iters
     if (length(fixed.par) > 0)
     try_length <- 1
     if (model != "EB")
     try_theta <- log(seq(exp(min(bounds.default)), exp(max(bounds.default)),
     length = try_length))
     else try_theta <- (seq((min(bounds.default)), (max(bounds.default)),
     length = try_length))
     try_LL <- double(length(try_theta))
     for (i in 1:length(try_theta)) {
     try_LL[i] <- model_LL(try_theta[i])
     }
     if (plot.LL.surface & !MC)
     try({
     plot(if (model != "EB")
     exp(try_theta)
     else try_theta, try_LL)
     }, silent = TRUE)
     starting.values.default <- try_theta[which.max(try_LL)]
     if (MC & length(fixed.par) == 0) {
     if (model != "EB")
     try_theta <- exp(try_theta)
     if (model != "EB")
     bounds.default <- exp(bounds.default)
     best <- double(round(length(try_theta)/2))
     for (i in 1:length(best)) {
     best[i] <- BIC(lm(try_LL ~ poly(try_theta, degree = i)))
     }
     degree <- which.min(best)
     LL_surface <- lm(try_LL ~ poly(try_theta, degree = degree))
     theta_start <- try_theta[which(predict(LL_surface) ==
     max(predict(LL_surface)))][1]
     o <- optim(par = theta_start, fn = function(X) predict.lm(LL_surface,
     newdata = data.frame(try_theta = X)), control = list(fnscale = -1),
     method = "L-BFGS-B", lower = bounds.default[1],
     upper = bounds.default[2])
     if (plot.LL.surface)
     try({
     plot(try_theta, try_LL)
     points(try_theta, predict.lm(LL_surface), col = "blue")
     points(o$par, o$value, col = "red", pch = 16)
     }, silent = TRUE)
     ret <- model_LL(par = if (model != "EB")
     log(o$par)
     else o$par, ret_pars = ret.level > 1)
     ret$model.par[[rownames(bounds.default)]] <- o$par
     if (ret.level > 1 & (model == "OUfixedRoot" | model ==
     "OUrandomRoot" | model == "OU")) {
     if (nrates.species > 1) {
     for (i in 1:nrates.species) ret$phylocov[[i]] <- 2 *
     ret$phylocov[[i]] * ret$alpha
     }
     else ret$phylocov <- 2 * ret$phylocov * ret$alpha
     }
     }
     else {
     if (length(fixed.par) > 0) {
     o <- list(par = starting.values.default, value = model_LL(starting.values.default))
     }
     else o <- optim(starting.values.default, model_LL,
     control = list(fnscale = -1), method = "L-BFGS-B",
     lower = bounds.default[1], upper = bounds.default[2])
     ret <- model_LL(par = o$par, ret_pars = ret.level >
     1)
     if (model != "EB")
     o$par <- exp(o$par)
     if (ret.level > 1) {
     ret$model.par[[rownames(bounds.default)]] <- o$par
     if (ret.level > 1 & (model == "OUfixedRoot" |
     model == "OUrandomRoot" | model == "OU")) {
     if (nrates.species > 1) {
     for (i in 1:nrates.species) ret$phylocov[[i]] <- 2 *
     ret$phylocov[[i]] * ret$model.par[[1]]
     }
     else ret$phylocov <- 2 * ret$phylocov * ret$model.par[[1]]
     }
     }
     }
     }
     else ret <- simple_BM(temp_tree = tree, ditree = ditree,
     ret_pars = ret.level > 1)
     if (ret.level <= 1)
     return(ret)
     if (ret.level > 1) {
     ret$method <- method
     sigma2.mult <- vector("list", nrates.species)
     if (nrates.species > 1) {
     for (i in 1:nrates.species) {
     if (nrates.traits < 2)
     sigma2.mult[[i]] <- mean(diag(ret$phylocov[[i]]))
     else {
     for (j in 1:nrates.traits) {
     group_j <- which(trait.groups == trait.group.levels[j])
     sigma2.mult[[i]] <- c(sigma2.mult[[i]], mean(diag(ret$phylocov[[i]])[group_j]) *
     if (subset) 1 else length(group_j))
     }
     names(sigma2.mult[[i]]) <- levels(trait.groups)
     }
     }
     names(sigma2.mult) <- levels(species.groups)
     }
     else {
     if (nrates.traits < 2)
     sigma2.mult <- mean(diag(ret$phylocov))
     else {
     sigma2.mult <- double()
     for (j in 1:nrates.traits) {
     group_j <- which(trait.groups == trait.group.levels[j])
     sigma2.mult <- c(sigma2.mult, mean(diag(ret$phylocov)[group_j]) *
     if (subset) 1 else length(group_j))
     }
     names(sigma2.mult) <- levels(trait.groups)
     }
     }
     if (ret.level < 3)
     ret$phylocov <- NULL
     ret <- c(ret, list(sigma2.mult = sigma2.mult, evo.model.args = evo.model.args))
     }
     class(ret) <- "evo.model"
     return(ret)
    }
    <bytecode: 0x5582e4b92e90>
    <environment: namespace:phylocurve>
     --- function search by body ---
    Function evo.model in namespace phylocurve has this body.
     ----------- END OF FAILURE REPORT --------------
    Fatal error: the condition has length > 1
Flavor: r-devel-linux-x86_64-debian-gcc

Version: 2.0.9
Check: whether package can be installed
Result: ERROR
    Installation failed.
Flavor: r-devel-windows-ix86+x86_64-gcc8

Version: 2.0.9
Check: whether package can be installed
Result: WARN
    Found the following significant warnings:
     Warning: 'rgl_init' failed, running with rgl.useNULL = TRUE
Flavor: r-release-osx-x86_64