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: 0x55838698efe0>
<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