Merge pull request #44 from ScottClaessens/multi-phylo

Including multiPhylo objects

DESCRIPTION

@@ -35,6 +35,7 @@ Imports:
     methods,
     phangorn,
     phaseR,
+    phytools,
     posterior,
     purrr,
     readr,

R/coev_calculate_delta_theta.R

@@ -110,9 +110,9 @@ coev_calculate_delta_theta <- function(object, response, predictor) {
   id_resp <- which(response  == names(object$variables))
   id_pred <- which(predictor == names(object$variables))
   # medians and median absolute deviations for all variables
-  eta  <- draws$eta[1:object$stan_data$N_tips,]
-  med  <- apply(eta, 2, posterior::rvar_median)
-  diff <- apply(eta, 2, posterior::rvar_mad)
+  eta  <- draws$eta[,1:object$stan_data$N_tips,]
+  med  <- apply(eta, 3, posterior::rvar_median)
+  diff <- apply(eta, 3, posterior::rvar_mad)
   # construct intervention values list
   values1 <- list()
   values2 <- list()

R/coev_fit.R

@@ -25,8 +25,9 @@
 #' @param id A character of length one identifying the variable in the data that
 #'   links rows to tips on the phylogeny. Must refer to a valid column name in
 #'   the data. The id column must exactly match the tip labels in the phylogeny.
-#' @param tree A phylogenetic tree object of class \code{phylo}. The tree must
-#'   be rooted and must include branch lengths.
+#' @param tree A phylogenetic tree object of class \code{phylo} or
+#'   \code{multiPhylo}. The tree(s) must be rooted and must include branch
+#'   lengths.
 #' @param effects_mat (optional) A boolean matrix with row and column names
 #'   exactly matching the variables declared for the model. If not specified,
 #'   all cross-lagged effects will be estimated in the model. If specified, the

R/coev_helpers.R

@@ -154,21 +154,31 @@ run_checks <- function(data, variables, id, tree, effects_mat,
   if (!(id %in% colnames(data))) {
     stop2("Argument 'id' is not a valid column name in the data.")
   }
-  # stop if tree is not a phylo object
-  if (!methods::is(tree, "phylo")) {
-    stop2("Argument 'tree' must be an phylogenetic tree object of class phylo.")
-  }
-  # stop if tree does not have branch length information
-  if (is.null(tree$edge.length)) {
-    stop2("Argument 'tree' does not include branch lengths.")
-  }
-  # stop if tree is not rooted
-  if (!ape::is.rooted(tree)) {
-    stop2("Argument 'tree' must be a rooted tree.")
+  # stop if tree is not a phylo or multiPhylo object
+  if (!(methods::is(tree, "phylo") | methods::is(tree, "multiPhylo"))) {
+    stop2(
+      paste0(
+        "Argument 'tree' must be an phylogenetic tree object of class phylo ",
+        "or multiPhylo."
+        )
+      )
   }
-  # stop if id in data does not match tree tip labels exactly
-  if (!identical(sort(unique(data[,id])), sort(tree$tip.label))) {
-    stop2("The id variable in the data does not match tree tip labels exactly.")
+  tree <- phytools::as.multiPhylo(tree)
+  for (t in 1:length(tree)) {
+    # stop if tree does not have branch length information
+    if (is.null(tree[[t]]$edge.length)) {
+      stop2("All trees in 'tree' argument must include branch lengths.")
+    }
+    # stop if tree is not rooted
+    if (!ape::is.rooted(tree[[t]])) {
+      stop2("All trees in 'tree' argument must be rooted.")
+    }
+    # stop if id in data does not match tree tip labels exactly
+    if (!identical(sort(unique(data[,id])), sort(tree[[t]]$tip.label))) {
+      stop2(
+        "The id variable in the data does not match tree tip labels exactly."
+        )
+    }
   }
   # stop if id in data contains missing values
   if (any(is.na(data[,id]))) {