From edb8885161bb3578edbee6529963b45aa72c7da5 Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Fri, 19 Jan 2024 15:56:17 +0000 Subject: [PATCH 01/15] Update multi.ace_fun.R --- R/multi.ace_fun.R | 42 +++++++++++++++++++++++++++++++++--------- 1 file changed, 33 insertions(+), 9 deletions(-) diff --git a/R/multi.ace_fun.R b/R/multi.ace_fun.R index 3ab10e6c..89f8bd1a 100755 --- a/R/multi.ace_fun.R +++ b/R/multi.ace_fun.R @@ -1,3 +1,27 @@ +## Set default arguments for continuous models +set.continuous.args.ace <- function(method, model, scaled, kappa, corStruct) { + continuous_args <- list(type = "continuous") + continuous_args$model <- ifelse(missing(model), "BM", model) + continuous_args$scaled <- ifelse(missing(scaled), TRUE, scaled) + continuous_args$kappa <- ifelse(missing(kappa), 1, kappa) + if(missing(corStruct)) { + continuous_args$corStruct <- NULL + } else { + continuous_args$corStruct <- corStruct + } + return(continuous_args) +} +## Set default arguments for continuous models with "models" as input (= method or model) +set.continuous.args.ace.models <- function(models, n) { + if(models == "BM" || models == "REML") { + ## Set everything default + return(replicate(n, set.continuous.args.ace(), simplify = FALSE)) + } else { + return(replicate(n, set.continuous.args.ace(method = models), simplify = FALSE)) + } +} + + ## Finding or adding node labels get.node.labels <- function(tree) { if(is.null(tree$node.label)) { @@ -105,16 +129,16 @@ castor.ace <- function(castor_args) { details <- NULL } -# stop("DEBUG multi.ace_fun::castor.ace") + # stop("DEBUG multi.ace_fun::castor.ace") -# asr_mk_model( tree = castor_args$tree, -# tip_states = castor_args$tip_states, -# Nstates = castor_args$Nstates, -# tip_priors = castor_args$tip_priors, -# rate_model = castor_args$rate_model, -# Ntrials = castor_args$Ntrials, -# check_input =castor_args$check_input, -# Nthreads = castor_args$Nthreads) + # asr_mk_model( tree = castor_args$tree, + # tip_states = castor_args$tip_states, + # Nstates = castor_args$Nstates, + # tip_priors = castor_args$tip_priors, + # rate_model = castor_args$rate_model, + # Ntrials = castor_args$Ntrials, + # check_input =castor_args$check_input, + # Nthreads = castor_args$Nthreads) ## Increase the number of trials if unsuccessful From 937a9c1c2b122708a0e18c7c5c4b807d1a6c8294 Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Fri, 19 Jan 2024 15:56:30 +0000 Subject: [PATCH 02/15] Multi ace restructure TODO --- DESCRIPTION | 4 +- NEWS.md | 16 +- R/multi.ace.R | 699 +++++++++++++++++++++++++++++++++++--------------- 3 files changed, 514 insertions(+), 205 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 69e7a390..237f5298 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -7,8 +7,8 @@ Authors@R: c(person("Thomas", "Guillerme", role = c("aut", "cre", "cph"), person("Jack", "Hatfield", role = c("aut", "cph")) ) Maintainer: Thomas Guillerme -Version: 1.8.2 -Date: 2024-01-15 +Version: 1.8.3 +Date: 2024-01-18 Description: A modular package for measuring disparity (multidimensional space occupancy). Disparity can be calculated from any matrix defining a multidimensional space. The package provides a set of implemented metrics to measure properties of the space and allows users to provide and test their own metrics. The package also provides functions for looking at disparity in a serial way (e.g. disparity through time) or per groups as well as visualising the results. Finally, this package provides several statistical tests for disparity analysis. Depends: R (>= 3.6.0), diff --git a/NEWS.md b/NEWS.md index 109cd149..ed63592d 100755 --- a/NEWS.md +++ b/NEWS.md @@ -1,9 +1,21 @@ -dispRity v1.8.2 (2024-01-15) +dispRity v1.8.3 (2024-01-18) ========================= +### NEW FEATURES + * Redesigned `multi.ace` to be more modular and handle both continuous and/or discrete characters. + - [ ] `models` arg now takes character names + complex list (for `ape::ace`) + - [ ] `models` selection is now done after sorting the characters (for default: discrete = ER, continuous = BM/REML) + - [ ] add explanations of new `models` argument in the doc + - [ ] change `castor.options` to `options` in general (with a named list: `list(castor = ..., ape = ...)`) + - [ ] new internal function for running `ape::ace` + - [ ] combine results more smartly + - [ ] `estimation.details` can now take arguments from both `castor` and `ape` + - [ ] new output `dispRity` format + - [ ] also update manual + ### MINOR IMPROVEMENTS - * `custom.subsets` can now take a logical vector for + * `custom.subsets` can now take a logical vector for the `group` argument ### BUG FIXES diff --git a/R/multi.ace.R b/R/multi.ace.R index afbe3d22..75a26495 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -2,10 +2,14 @@ #' #' @description Fast ancestral states estimations run on multiple trees using the Mk model from castor::asr_mk_model. #' -#' @param data A \code{matrix} or \code{list} with the characters for each taxa. +#' @param data A \code{matrix}, \code{data.frame} or \code{list} with the characters for each taxa. #' @param tree A \code{phylo} or \code{mutiPhylo} object (if the \code{tree} argument contains node labels, they will be used to name the output). #' @param models A \code{vector} of models to be passed to \code{castor::asr_mk_model}. -#If left empty, the it will use the \code{\link{fit.ace.model}} function to find the best model using the first tree. See details. + + +#TODO: add info about continuous characters + + #' @param threshold either \code{logical} for applying a relative threshold (\code{TRUE} - default) or no threshold (\code{FALSE}) or a \code{numeric} value of the threshold (e.g. 0.95). See details. #' @param special.tokens optional, a named \code{vector} of special tokens to be passed to \code{\link[base]{grep}} (make sure to protect the character with \code{"\\\\"}). By default \code{special.tokens <- c(missing = "\\\\?", inapplicable = "\\\\-", polymorphism = "\\\\&", uncertainty = "\\\\/")}. Note that \code{NA} values are not compared and that the symbol "@" is reserved and cannot be used. #' @param special.behaviours optional, a \code{list} of one or more functions for a special behaviour for \code{special.tokens}. See details. @@ -18,8 +22,16 @@ #' #' @details #' -#' The \code{models} argument can be a single or a list of transition \code{matrix}, a single or a a vector of built-in model(s) (see below) or a list of both matrices and built-in models: -#' The available built-in models in \code{castor::asr_mk_model} are: +#' Depending on the type of characters \code{models} argument can be either: +#' \itemize{ +#' \item the name of a single model to apply to all characters (if all characters are discrete or all are continuous); see below for the list of available names. For example \code{models = "ER"} applies the Equal Rates model to all characters (assuming they are all discrete characters). +#' \item a vector of model names to apply to different type of characters (see below for the list). For example \code{models = c("ER", "ER", "BM")} applies the Equal Rates model to the two first characters (discrete) and the \code{"BM"} model to the third character (continuous). +#' \item a transition \code{"matrix"} to be applied to all characters (if discrete). For example \code{models = matrix(0.2, 2, 2)}. +#' \item an single named list of arguments to be applied to all characters by passing it to \code{ape::ace} (if continuous). For example \code{models = list(method = "GLS", corStruct = corBrownian(1, my_tree))}. +#' \item an un-ambiguous list of arguments to be passed to either \code{castor::asr_mk_model} (discrete characters) or \code{ape::ace} (continuous characters). For example \code{models = list("char1" = list(transition_matrix = matrix(0.2, 2, 2)), "char2" = list(method = "GLS", corStruct = corBrownian(1, my_tree)))} to be specifically passed to the characters named "char1" and "char2" +#'} +#' +#' The available built-in models for discrete characters in \code{castor::asr_mk_model} are: #' \itemize{ #' \item \code{"ER"} for all equal rates #' \item \code{"SYM"} for symmetric rates @@ -28,7 +40,14 @@ #' \item \code{"SRD"} different stepwise transitions #' } #' See directly \code{castor::asr_mk_model} for more models. -# TODO: add note about fit.ace.model +#' +#' The available built-in models and methods for continuous characters in \code{ape::ace} are: +#' \itemize{ +#' \item \code{"BM"} model: for a default Brownian Motion with the "REML" method +#' \item \code{"REML"} method: for a default Brownian Motion with the "REML" method (same as above) +#' \item \code{"ML"} method: for a default Brownian Motion with the "ML" method +#' \item \code{"pic"} method: for a default Brownian Motion with the "pic" (least squared) method +#'} #' #' The \code{threshold} option allows to convert ancestral states likelihoods into discrete states. When \code{threshold = FALSE}, the ancestral state estimated is the one with the highest likelihood (or at random if likelihoods are equal). When \code{threshold = TRUE}, the ancestral state estimated are all the ones that are have a scaled likelihood greater than the maximum observed scaled likelihood minus the inverse number of possible states (i.e. \code{select_state >= (max(likelihood) - 1/n_states)}). This option makes the threshold selection depend on the number of states (i.e. if there are more possible states, a lower scaled likelihood for the best state is expected). Finally using a numerical value for the threshold option (e.g. \code{threshold = 0.95}) will simply select only the ancestral states estimates with a scaled likelihood equal or greater than the designated value. This option makes the threshold selection absolute. Regardless, if more than one value is select, the uncertainty token (\code{special.tokens["uncertainty"]}) will be used to separate the states. If no value is selected, the uncertainty token will be use between all observed characters (\code{special.tokens["uncertainty"]}). #' @@ -47,6 +66,13 @@ #' @return #' Returns a \code{"matrix"} or \code{"list"} of ancestral states. By default, the function returns the ancestral states in the same format as the input \code{matrix}. This can be changed using the option \code{output = "matrix"} or \code{"list"} to force the class of the output. #' To output the combined ancestral states and input, you can use \code{"combined"} (using the input format) or \code{"combined.matrix"} or \code{"combined.list"}. +#' \emph{NOTE} that if the input data had multiple character types (continuous and discrete) and that \code{"matrix"} or \code{"combined.matrix"} output is requested, the function returns a \code{"data.frame"}. + + +#TODO: add dispRity format + + + # To output the light version to be passed to \code{dispRity} functions (a list of two elements: 1) the input \code{matrix} and 2) a list of ancestral states matrices) you can use \code{output = "dispRity"}. #' #' @examples @@ -134,7 +160,7 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token ## matrix matrix <- data - input_class <- check.class(matrix, c("matrix", "list")) + input_class <- check.class(matrix, c("matrix", "list", "data.frame")) ## Convert the matrix if not a list class_matrix <- class(matrix) if(class_matrix[[1]] == "list") { @@ -156,7 +182,6 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token stop(paste0("Some names in the data or the tree(s) are not matching.\nYou can use dispRity::clean.data(", as.expression(match_call$data), ", ", as.expression(match_call$tree), ") to find out more.")) } - ## Find the node labels (and eventually add them to the trees) node_labels <- lapply(tree, get.node.labels) ## Split the trees and the labels @@ -164,61 +189,15 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token class(tree) <- "multiPhylo" node_labels <- lapply(node_labels, `[[`, 2) - ## The available models for castor - available_models <- c("ER", "SYM", "ARD", "SUEDE", "SRD") - ## models - if(missing(models)) { - models <- replicate(n_characters, "ER", simplify = FALSE) - } else { - ## What is the model list - model_class <- check.class(models, c("character", "matrix", "list")) - - ## Models is a list of one matrix or one character - if(model_class == "list") { - ## Check the class of of the list - list_class <- unique(unlist(lapply(models, class))) - if(length(list_class) == 1 || all(list_class %in% c("matrix", "array"))) { - model_class <- ifelse(list_class[1] %in% c("character", "matrix"), model_class[1], stop.call(call = "", msg = "models must be a list containing characters of matrices.")) - } else { - check.length(models, n_characters, msg = paste0(" should be list of characters or/and matrices of length ", ncol(matrix), ".")) - ## Check all models - silent <- lapply(models, check.model.class, available_models) - } - } - ## Models are character or matrix - switch(model_class, - "character" = { - ## Check the model names - available_models <- c("ER", "SYM", "ARD", "SUEDE", "SRD") - silent <- sapply(models, check.method, all_arguments = available_models, msg = "model") - if(length(models) == 1) { - models <- replicate(n_characters, models, simplify = FALSE) - } else { - check.length(models, n_characters, msg = paste0(" should be a single character string or a vector of models for each ", ncol(matrix), " characters.")) - } - }, - "matrix" = { - ## Check the class of the matrix content - check.class(c(models), c("numeric", "integer"), msg = "models must be numerical matrices") - models <- replicate(n_characters, models, simplify = FALSE) - } - ) - } - ## castor.options - if(missing(castor.options)) { - ## No options - castor.options <- NULL - } else { - ## must be list with names - check.class(castor.options, "list") - if(is.null(names(castor.options))) { - stop("castor.options must be a named list of options for castor::asr_mk_model().", call. = FALSE) - } - } + ######### + ## + ## Handle the other options (threshold, brlen, verbose, parallel, output, estimation.details) + ## + ######### ## threshold check.class(threshold, c("logical", "numeric")) @@ -235,67 +214,7 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token threshold.type <- "absolute" } - ## Special tokens - if(missing(special.tokens)) { - special.tokens <- character() - } - check.class(special.tokens, "character") - not.exist <- function(special.tokens, token) { - name_token <- names(special.tokens[token]) - return(is.null(name_token) || is.na(name_token)) - } - if(not.exist(special.tokens, "missing")) { - special.tokens["missing"] <- "\\?" - } - if(not.exist(special.tokens, "inapplicable")) { - special.tokens["inapplicable"] <- "\\-" - } - if(not.exist(special.tokens, "polymorphism")) { - special.tokens["polymorphism"] <- "\\&" - } - if(not.exist(special.tokens, "uncertainty")) { - special.tokens["uncertainty"] <- "\\/" - } - - ## Checking for the reserved character - reserved <- c("\\@", "@") %in% special.tokens - if(any(reserved)) { - stop("special.tokens cannot contain the character '@' since it is reserved for the dispRity::char.diff function.") - } - - ## Checking whether the special.tokens are unique - if(length(unique(special.tokens)) != length(special.tokens)) { - stop("special.tokens cannot contain duplicated tokens.") - } - - ## If any special token is NA, convert them as "N.A" temporarily - if(any(is.na(special.tokens))) { - matrix <- ifelse(is.na(matrix), "N.A", matrix) - special.tokens[is.na(special.tokens)] <- "N.A" - } - - ## Special behaviours - if(missing(special.behaviours)) { - special.behaviours <- list() - } - check.class(special.behaviours, "list") - if(is.null(special.behaviours$missing)) { - special.behaviours$missing <- function(x,y) return(y) - } - if(is.null(special.behaviours$inapplicable)) { - special.behaviours$inapplicable <- function(x,y) return(y) - } - if(is.null(special.behaviours$polymorphism)) { - special.behaviours$polymorphism <- function(x,y) return(strsplit(x, split = "\\&")[[1]]) - } - if(is.null(special.behaviours$uncertainty)) { - special.behaviours$uncertainty <- function(x,y) return(strsplit(x, split = "\\/")[[1]]) - } - - ## Match the behaviours and tokens in the same order - special.behaviours <- special.behaviours[sort(names(special.behaviours))] - special.tokens <- special.tokens[sort(names(special.tokens))] - + ## brlen multiplier if(!missing(brlen.multiplier)) { ## Check class @@ -364,109 +283,483 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token check.method(estimation.details, return_args, msg = "estimation.details") } - ## Convert the potential missing data + ######### + ## Handle the characters + ######### + + ## Preparing the data if(verbose) cat("Preparing the data:.") - ## Translate the characters using the special behaviours - characters <- unlist(apply(do.call(cbind, apply(matrix, 2, convert.bitwise, special.tokens, special.behaviours, bitwise = FALSE)), 2, list), recursive = FALSE) - if(verbose) cat(".") - - ## Get a list of character states - characters_states <- lapply(characters, function(char) sort(unique(na.omit(unlist(char))))) - if(verbose) cat(".") + ## Detecting the continuous or discrete characters + character_is_continuous <- logical() + ## Looping to allow dropping the levels from matrix + for(col in 1:ncol(matrix)) { + character_is_continuous <- c(character_is_continuous, is.numeric(matrix[, col, drop = TRUE])) + } + do_discrete <- do_continuous <- FALSE + continuous_char_ID <- discrete_char_ID <- numeric() + + ## Split the matrices by character types + if(any(character_is_continuous)) { + ## Split the matrix for continuous characters + matrix_continuous <- matrix[, character_is_continuous] + n_characters_continuous <- sum(character_is_continuous) + do_continuous <- TRUE + continuous_char_ID <- which(character_is_continuous) + } + if(any(!character_is_continuous)) { + ## Split the matrix for discrete characters + matrix_discrete <- matrix[, !character_is_continuous] + ## Convert into characters + matrix_discrete <- apply(matrix_discrete, 2, as.character) + rownames(matrix_discrete) <- rownames(matrix) + n_characters_discrete <- sum(!character_is_continuous) + do_discrete <- TRUE + discrete_char_ID <- which(!character_is_continuous) + } - ## Find invariant characters - invariants <- which(lengths(characters_states) < 2) - if(length(invariants) > 0) { + ## Handle the tokens + if(do_discrete) { + ## Special tokens + if(missing(special.tokens)) { + special.tokens <- character() + } + check.class(special.tokens, "character") + not.exist <- function(special.tokens, token) { + name_token <- names(special.tokens[token]) + return(is.null(name_token) || is.na(name_token)) + } + if(not.exist(special.tokens, "missing")) { + special.tokens["missing"] <- "\\?" + } + if(not.exist(special.tokens, "inapplicable")) { + special.tokens["inapplicable"] <- "\\-" + } + if(not.exist(special.tokens, "polymorphism")) { + special.tokens["polymorphism"] <- "\\&" + } + if(not.exist(special.tokens, "uncertainty")) { + special.tokens["uncertainty"] <- "\\/" + } - ## Stop if they are only invariant characters - if(length(invariants) == n_characters) { - stop.call(match_call$data, " contains only invariant characters.") + ## Checking for the reserved character + reserved <- c("\\@", "@") %in% special.tokens + if(any(reserved)) { + stop("special.tokens cannot contain the character '@' since it is reserved for the dispRity::char.diff function.") } - ## Remove the characters - invariant_characters <- characters[invariants] - invariant_characters_states <- characters_states[invariants] - characters <- characters[-invariants] - characters_states <- characters_states[-invariants] + ## Checking whether the special.tokens are unique + if(length(unique(special.tokens)) != length(special.tokens)) { + stop("special.tokens cannot contain duplicated tokens.") + } - ## Remove the models - models <- models[-invariants] + ## If any special token is NA, convert them as "N.A" temporarily + if(any(is.na(special.tokens))) { + matrix_discrete <- ifelse(is.na(matrix_discrete), "N.A", matrix_discrete) + special.tokens[is.na(special.tokens)] <- "N.A" + } - ## Tell the user - warning(paste0("The character", ifelse(length(invariants > 1), "s", "") , " ", paste0(invariants, collapse = ", "), " are invariant (using the current special behaviours for special characters) and are simply duplicated for each node."), call. = FALSE) + ## Special behaviours + if(missing(special.behaviours)) { + special.behaviours <- list() + } + check.class(special.behaviours, "list") + if(is.null(special.behaviours$missing)) { + special.behaviours$missing <- function(x,y) return(y) + } + if(is.null(special.behaviours$inapplicable)) { + special.behaviours$inapplicable <- function(x,y) return(y) + } + if(is.null(special.behaviours$polymorphism)) { + special.behaviours$polymorphism <- function(x,y) return(strsplit(x, split = "\\&")[[1]]) + } + if(is.null(special.behaviours$uncertainty)) { + special.behaviours$uncertainty <- function(x,y) return(strsplit(x, split = "\\/")[[1]]) + } + + ## Match the behaviours and tokens in the same order + special.behaviours <- special.behaviours[sort(names(special.behaviours))] + special.tokens <- special.tokens[sort(names(special.tokens))] + + ## Translate the characters using the special behaviours + characters_discrete <- unlist(apply(do.call(cbind, apply(matrix_discrete, 2, convert.bitwise, special.tokens, special.behaviours, bitwise = FALSE)), 2, list), recursive = FALSE) + if(verbose) cat(".") + + ## Get a list of character states + characters_states <- lapply(characters_discrete, function(char) sort(unique(na.omit(unlist(char))))) + if(verbose) cat(".") + + ## Find invariant characters + invariants <- which(lengths(characters_states) < 2) + + ## Handle invariant characters + if(length(invariants) > 0) { + has_invariants <- TRUE + + ## Stop if they are only invariant characters + if(length(invariants) == n_characters_discrete) { + warning(match_call$data, " contains only invariant discrete characters.") + } + + ## Remove the characters + invariant_characters <- characters_discrete[invariants] + invariant_characters_states <- characters_states[invariants] + characters_discrete <- characters_discrete[-invariants] + characters_states <- characters_states[-invariants] + + ## Tell the user + warning(paste0("The character", ifelse(length(invariants) > 1, "s", "") , " ", paste0(invariants, collapse = ", "), ifelse(length(invariants) > 1, " are", " is"), " invariant (using the current special behaviours for special characters) and", ifelse(length(invariants) > 1, " are", " is"), " simply duplicated for each node."), call. = FALSE) + } else { + has_invariants <- FALSE + invariant_characters_states <- NULL + } + + if(verbose) cat(".") + + ## Get the character tables + characters_tables <- mapply(convert.char.table, characters_discrete, characters_states, SIMPLIFY = FALSE) + if(verbose) cat(".") + } + + ## Handle the continuous characters + if(do_continous) { + ## Make the continuous characters as lists + characters_continuous <- unlist(apply(matrix_continuous, 2, list), recursive = FALSE) + if(verbose) cat(".") + } + + ######### + ## Handle the models for each character + ######### + + ## Default (missing models) + if(missing(models)) { + if(do_discrete) { + models_discrete <- replicate(n_characters_discrete, "ER", simplify = FALSE) + } + if(do_continous) { + models_continuous <- replicate(n_characters_continuous, set.continuous.args.ace(), simplify = FALSE) + } } else { - invariant_characters_states <- NULL + ## Input models + models_class <- check.class(models, c("character", "list", "matrix")) + + ## Models is a vector of models + if(models_class == "character") { + ## Check the different models + available_models_discrete <- c("ER", "SYM", "ARD", "SUEDE", "SRD") + available_models_continuous <- c("BM", "REML", "ML", "pic") + + ## Unique model + if(length(models) == 1) { + if(do_discrete && !do_continuous) { + check.method(models, available_models_discrete, msg = "model applied to all discrete characters") + models_discrete <- replicate(n_characters_discrete, models, simplify = FALSE) + } + if(!do_discrete && do_continuous) { + check.method(models, available_models_continuous, msg = "model applied to all continuous characters") + models_continuous <- set.continuous.args.ace.models(models, n = n_characters_continuous) + } + if(do_discrete && do_continuous) { + stop("Only one model is specified but both discrete and continuous characters are detected.") + } + } else { + ## Vector of models + if(length(models) != n_characters) { + stop(paste0("Incorrect number of models specified: ", length(models), " models for ", n_characters, " characters.")) + } else { + check.method(models, c(available_models_discrete, available_models_continuous), msg = "models applied to characters") + ## Check models per character types + ## Discrete + if(sum(models %in% available_models_discrete) != n_characters_discrete) { + stop(paste0("Incorrect number of models specified: ", sum(models %in% available_models_discrete), " models for ", n_characters, " discrete characters.")) + } else { + ## Discrete models (valid) + models_discrete <- as.list(models[models %in% available_models_discrete]) + } + ## Continuous + if(sum(models %in% available_models_continuous) != n_characters_continuous) { + stop(paste0("Incorrect number of models specified: ", sum(models %in% available_models_continuous), " models for ", n_characters, " continuous characters.")) + } else { + ## Continuous models (valid) + models_continuous <- sapply(models[models %in% available_models_continuous], set.continuous.args.ace.models, n = 1) + } + } + } + + ##TODO + ## Check the length and if they match the available models and type of characters + + ## Check the different models + available_models_discrete <- c("ER", "SYM", "ARD", "SUEDE", "SRD") + available_models_continuous <- c("BM", "REML", "ML", "pic") + + + if(do_discrete) { + models_discrete <- as.list(models[discrete_char_ID]) + } + if(do_continous) { + ## Make the models list using set.continuous.args.ace + models_continuous <- + } + } + + ## Models is a transition matrix (discrete only) + if(models_class == "matrix") { + if(do_continous) { + stop("Transition matrices can be used as models only for discrete characters.") + } else { + models_discrete <- replicate(n_characters_discrete, models, simplify = FALSE) + } + } + + ## Models is a complicated list + if(models_class == "list") { + if(length(models) == 1) { + if(do_discrete && do_continuous) { + stop("Only one model is specified but both discrete and continuous characters are detected.") + } + ## Set the models for discrete + if(do_discrete) { + models_discrete <- replicate(n_characters_discrete, models, simplify = FALSE) + } + ## Set the models for continuous + if(do_continous) { + models_continuous <- replicate(n_characters_continuous, do.call(set.continuous.args.ace, models), simplify = FALSE) + } + } + + ## Models is a list of models + check.length(models, n_characters, msg = paste0(" list must be the same length as the number of characters (", n_characters, ").")) + ## Separate the models per type + if(do_discrete) { + models_discrete <- models[discrete_char_ID] + } + if(do_continous) { + models_continuous <- models[continuous_char_ID] + ## Format correctly + models_continuous <- lapply(models_continuous, function(x) do.call(set.continuous.args.ace, x)) + } + } + + ## Remove invariant characters + if(do_discrete && has_invariants) { + ## Remove the models + models_discrete <- models_discrete[-invariants] + } } - if(verbose) cat(".") - ## Get the character tables - characters_tables <- mapply(convert.char.table, characters, characters_states, SIMPLIFY = FALSE) - if(verbose) cat(".") - ## Set up the arguments for one tree - args_list <- mapply(make.args, characters_tables, characters_states, models, - MoreArgs = list(castor.options, cores, estimation.details), SIMPLIFY = FALSE) - ## Add up the tree arguments - add.tree <- function(tree, args_list) { - return(lapply(args_list, function(arg, tree) c(arg, tree = list(tree)), tree)) + + + + + + + + + ######### + ## + ## Handle the options + ## + ######### + + + + + + + + + + + + + ## castor.options + warning("TODO: change castor options to just options (parsed to castor or ape)") + warning("TODO: multi.ace. ape options are: CI, ") + + + if(missing(castor.options)) { + ## No options + castor.options <- NULL + } else { + ## must be list with names + check.class(castor.options, "list") + if(is.null(names(castor.options))) { + stop("castor.options must be a named list of options for castor::asr_mk_model().", call. = FALSE) + } } - tree_args_list <- lapply(tree, add.tree, args_list) - if(verbose) cat("Done.\n") - if(do_parallel) { - ## Remove verbose - if(verbose) { - cat(paste0("Running the estimation for ", length(tree), " tree", ifelse(length(tree) > 1, "s ", " "), "using ", cores, " core", ifelse(cores == 1, "...", "s..."))) - was_verbose <- TRUE - verbose <- FALSE + + + + + + + + + + ## Setting up discrete characters + tree_args_list_discrete <- tree_args_list_continous <- NULL + + ## Running discrete characters + if(do_discrete) { + + + ## Set up the arguments for one tree + args_list <- mapply(make.args, characters_tables, characters_states, models, + MoreArgs = list(castor.options, cores, estimation.details), SIMPLIFY = FALSE) + + ## Add up the tree arguments + add.tree <- function(tree, args_list) { + return(lapply(args_list, function(arg, tree) c(arg, tree = list(tree)), tree)) + } + tree_args_list <- lapply(tree, add.tree, args_list) + + if(verbose) cat("Done.\n") + + if(do_parallel) { + ## Remove verbose + if(verbose) { + cat(paste0("Running the estimation for ", length(tree), " tree", ifelse(length(tree) > 1, "s ", " "), "using ", cores, " core", ifelse(cores == 1, "...", "s..."))) + was_verbose <- TRUE + verbose <- FALSE + } else { + was_verbose <- FALSE + } + + ## Set up the cluster + cluster <- parallel::makeCluster(cores) + + ## Get the current environment + current_env <- environment() + ## Get the export lists + export_arguments_list <- c("tree_args_list", + "special.tokens", + "invariants", + "threshold.type", + "threshold", + "verbose", + "characters_states", + "invariant_characters_states") + export_functions_list <- c("one.tree.ace", + "castor.ace", + "tree.data.update", + "add.state.names", + "translate.likelihood") + + ## Export from this environment + parallel::clusterExport(cluster, c(export_arguments_list, export_functions_list), envir = current_env) + + ## Call the cluster + results_out <- parLapply(cl = cluster, tree_args_list, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) + + ## Stop the cluster + parallel::stopCluster(cluster) + + ## Reactivate the verbose + if(was_verbose) { + cat("Done.") + verbose <- TRUE + } } else { - was_verbose <- FALSE + ## Running the ACE for each tree + results_out <- lapply(tree_args_list, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) } + ## Separating results and details + details_out_discrete <- lapply(results_out, `[[`, 2) + results_out_discrete <- lapply(results_out, `[[`, 1) + } - ## Set up the cluster - cluster <- parallel::makeCluster(cores) - - ## Get the current environment - current_env <- environment() - ## Get the export lists - export_arguments_list <- c("tree_args_list", - "special.tokens", - "invariants", - "threshold.type", - "threshold", - "verbose", - "characters_states", - "invariant_characters_states") - export_functions_list <- c("one.tree.ace", - "castor.ace", - "tree.data.update", - "add.state.names", - "translate.likelihood") - - ## Export from this environment - parallel::clusterExport(cluster, c(export_arguments_list, export_functions_list), envir = current_env) - - ## Call the cluster - results_out <- parLapply(cl = cluster, tree_args_list, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) - - ## Stop the cluster - parallel::stopCluster(cluster) - - ## Reactivate the verbose - if(was_verbose) { - cat("Done.") - verbose <- TRUE + + + ## Running continuous characters + if(do_continous) { + + models + tree + data + #TODO: set up the code for continuous characters (model and shit) + warning("TODO: multi.ace") + + + ## Set the model list (should be handled above - when selecting the models) + character_models <- list("char1" = list(method = "REML", + model = "BM", + scaled = TRUE), + "char2" = list(method = "REML", + model = "BM", + scaled = TRUE), + "char3" = list(method = "REML", + model = "BM", + scaled = TRUE), + "char4" = list(method = "REML", + model = "BM", + scaled = TRUE), + "char5" = list(method = "REML", + model = "BM", + scaled = TRUE)) + + ## Set the characters as a list (should be handled above - when splitting discrete characters) + characters <- apply(data, 2, list) + + ## Create the character arguments + character_args <- mapply(function(character, ace.args) return(c(x = character, ace.args)), characters, character_models, SIMPLIFY = FALSE) + + ## Create the character and tree arguments + tree_character_args <- list() + for(one_tree in 1:length(tree)) { + tree_character_args[[one_tree]] <- lapply(character_args, function(character, tree) {character$phy <- tree; return(character)}, tree[[one_tree]]) } - } else { - ## Running the ACE for each tree - results_out <- lapply(tree_args_list, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) + + ## Run all the ace + raw_estimates <- lapply(tree_character_args, lapply, function(x) do.call(what = ape::ace, args = x)) + ## Remove the ugly call + raw_estimates <- lapply(raw_estimates, lapply, function(x) {x$call <- "ape::ace"; return(x)}) + + } + + + ## Standardise the whole procedure to: + #1- characters are sorted in two categories + #2- characters in each category get their models sorted + #3- characters from both categories get their trees attributed + #4- characters and tree are run with do.call + + estimates_continuous <- estimates_discrete <- NULL + if(do_continous) { + ## Run all the ace + estimates_continuous <- lapply(tree_character_args_continuous, lapply, function(x) do.call(what = ape::ace, args = x)) + ## Remove the ugly call + estimates_continuous <- lapply(estimates_continuous, lapply, function(x) {x$call <- "ape::ace"; return(x)}) } - ## Separating results and details - details_out <- lapply(results_out, `[[`, 2) - results_out <- lapply(results_out, `[[`, 1) + if(do_discrete) { + ## Run all the castor + estimates_discrete <- lapply(tree_character_args_continuous, lapply, function(x) do.call(what = castor::asr_mk_model, args = x)) + ## Remove the ugly call + estimates_discrete <- lapply(estimates_discrete, lapply, function(x) {x$call <- "ape::ace"; return(x)}) + } + + + + + + + + + + + + + + ## Reorder the results + #TODO: reorder the results using continuous_char_ID and discrete_char_ID: like + details_out <- details_out[c(continuous_char_ID, discrete_char_ID)] + results_out <- results_out[c(continuous_char_ID, discrete_char_ID)] + warning("TODO: multi.ace") ## Output a matrix make.matrix <- function(results) { @@ -494,6 +787,10 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token #dispRity = return(list("tips" = matrix, "nodes" = output_matrix)) ) + #TODO: set up "dispRity" output + #TODO: make sure node.labels available in output + warning("DEBUG: multi.ace") + ## Results out if(is.null(estimation.details)) { if(length(tree) == 1) { From 55b1491dfe864b4fe279708c707cb2c2770d13ba Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Fri, 19 Jan 2024 15:56:37 +0000 Subject: [PATCH 03/15] Update test-multi.ace.R --- tests/testthat/test-multi.ace.R | 30 ++++++++++++++++++++++++++++++ 1 file changed, 30 insertions(+) diff --git a/tests/testthat/test-multi.ace.R b/tests/testthat/test-multi.ace.R index 366b5bba..a4d0189d 100755 --- a/tests/testthat/test-multi.ace.R +++ b/tests/testthat/test-multi.ace.R @@ -412,3 +412,33 @@ test_that("multi.ace works", { # expect_is(results$details[[2]]$loglikelihood[[1]], "numeric") }) + +test_that("multi.ace works with continuous and mix", { + + ## The tree + tree <- rcoal(15) + tree <- makeNodeLabel(tree) + ## The matrix + data <- space.maker(elements = 15, dimensions = 5, distribution = rnorm, elements.name = tree$tip.label) + + ## Run the multi.ace on the continuous data + test <- multi.ace(data = data, tree = tree, output = "combined.matrix") + + ## Works well for continuous + expect_is(test, "matrix") + expect_equal(dim(test), c(15+14, 5)) + expect_equal(sort(rownames(test)), sort(c(tree$tip.label, tree$node.label))) + expect_equal(unique(apply(test, 2, class)), "numeric") + + ## Mixed characters + data <- as.data.frame(data) + data <- cbind(data, "new_char" = as.character(sample(1:2, 15, replace = TRUE))) + data <- cbind(data, "new_char2" = as.character(sample(1:2, 15, replace = TRUE))) + + ## Works well for mixed characters + test <- multi.ace(data = data, tree = tree, output = "combined.matrix") + expect_is(test, "data.frame") + expect_equal(dim(test), c(15+14, 6)) + expect_equal(sort(rownames(test)), sort(c(tree$tip.label, tree$node.label))) + expect_equal(unique(apply(test, 2, class)), c("numeric", "character")) +}) \ No newline at end of file From 03e69035eddd167a0da52dbcec7a2776ce432ced Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Mon, 22 Jan 2024 12:55:22 +0000 Subject: [PATCH 04/15] chipped away at the updates --- NEWS.md | 4 +- R/multi.ace.R | 338 +++++++++++++++----------------- R/multi.ace_fun.R | 14 +- man/multi.ace.Rd | 23 ++- tests/testthat/test-multi.ace.R | 2 +- 5 files changed, 190 insertions(+), 191 deletions(-) diff --git a/NEWS.md b/NEWS.md index ed63592d..336e7f05 100755 --- a/NEWS.md +++ b/NEWS.md @@ -3,8 +3,10 @@ dispRity v1.8.3 (2024-01-18) ### NEW FEATURES * Redesigned `multi.ace` to be more modular and handle both continuous and/or discrete characters. + - inc. changed default `models` argument from `"ER"` to `"ER"` or `"BM"` depending on the character types. + - inc. change `"castor.options"` to the more generic `"options.args"` - [ ] `models` arg now takes character names + complex list (for `ape::ace`) - - [ ] `models` selection is now done after sorting the characters (for default: discrete = ER, continuous = BM/REML) + - [x] `models` selection is now done after sorting the characters (for default: discrete = ER, continuous = BM/REML) - [ ] add explanations of new `models` argument in the doc - [ ] change `castor.options` to `options` in general (with a named list: `list(castor = ..., ape = ...)`) - [ ] new internal function for running `ape::ace` diff --git a/R/multi.ace.R b/R/multi.ace.R index 75a26495..7f7f94db 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -17,7 +17,7 @@ #' @param verbose \code{logical}, whether to be verbose (\code{TRUE}) or not (\code{FALSE} - default). #' @param parallel \code{logical}, whether to use parallel algorithm (\code{TRUE}) or not (\code{FALSE} - default). #' @param output optional, see Value section below. -#' @param castor.options optional, a named list of options to be passed to function called by \code{castor::asr_mk_model}. +#' @param options.args optional, a named list of options to be passed to function called by \code{castor::asr_mk_model}. #' @param estimation.details optional, whether to also return the details for each estimation as returned by \code{castor::asr_mk_model}. This argument can be left \code{NULL} (default) or be any combination of the elements returned by \code{castor::asr_mk_model} (e.g. \code{c("loglikelihood", "transition_matrix")}). #' #' @details @@ -152,7 +152,7 @@ #' @author Thomas Guillerme #' @export -multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.tokens, special.behaviours, brlen.multiplier, verbose = FALSE, parallel = FALSE, output, castor.options, estimation.details = NULL) { +multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, special.behaviours, brlen.multiplier, verbose = FALSE, parallel = FALSE, output, options.args, estimation.details = NULL) { match_call <- match.call() @@ -319,6 +319,7 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token } ## Handle the tokens + # special.tokens <- character(); special.behaviours <- list() ; warning("DEBUG: multi.ace") if(do_discrete) { ## Special tokens if(missing(special.tokens)) { @@ -422,9 +423,9 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token } ## Handle the continuous characters - if(do_continous) { + if(do_continuous) { ## Make the continuous characters as lists - characters_continuous <- unlist(apply(matrix_continuous, 2, list), recursive = FALSE) + characters_continuous <- apply(matrix_continuous, 2, list) if(verbose) cat(".") } @@ -437,7 +438,7 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token if(do_discrete) { models_discrete <- replicate(n_characters_discrete, "ER", simplify = FALSE) } - if(do_continous) { + if(do_continuous) { models_continuous <- replicate(n_characters_continuous, set.continuous.args.ace(), simplify = FALSE) } } else { @@ -486,27 +487,11 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token } } } - - ##TODO - ## Check the length and if they match the available models and type of characters - - ## Check the different models - available_models_discrete <- c("ER", "SYM", "ARD", "SUEDE", "SRD") - available_models_continuous <- c("BM", "REML", "ML", "pic") - - - if(do_discrete) { - models_discrete <- as.list(models[discrete_char_ID]) - } - if(do_continous) { - ## Make the models list using set.continuous.args.ace - models_continuous <- - } } ## Models is a transition matrix (discrete only) if(models_class == "matrix") { - if(do_continous) { + if(do_continuous) { stop("Transition matrices can be used as models only for discrete characters.") } else { models_discrete <- replicate(n_characters_discrete, models, simplify = FALSE) @@ -524,7 +509,7 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token models_discrete <- replicate(n_characters_discrete, models, simplify = FALSE) } ## Set the models for continuous - if(do_continous) { + if(do_continuous) { models_continuous <- replicate(n_characters_continuous, do.call(set.continuous.args.ace, models), simplify = FALSE) } } @@ -535,7 +520,7 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token if(do_discrete) { models_discrete <- models[discrete_char_ID] } - if(do_continous) { + if(do_continuous) { models_continuous <- models[continuous_char_ID] ## Format correctly models_continuous <- lapply(models_continuous, function(x) do.call(set.continuous.args.ace, x)) @@ -552,13 +537,6 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token - - - - - - - ######### ## ## Handle the options @@ -570,25 +548,19 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token - - - - - - - ## castor.options + ## options.args warning("TODO: change castor options to just options (parsed to castor or ape)") warning("TODO: multi.ace. ape options are: CI, ") - if(missing(castor.options)) { + if(missing(options.args)) { ## No options - castor.options <- NULL + options.args <- NULL } else { ## must be list with names - check.class(castor.options, "list") - if(is.null(names(castor.options))) { - stop("castor.options must be a named list of options for castor::asr_mk_model().", call. = FALSE) + check.class(options.args, "list") + if(is.null(names(options.args))) { + stop("options.args must be a named list of options for castor::asr_mk_model().", call. = FALSE) } } @@ -596,169 +568,171 @@ multi.ace <- function(data, tree, models = "ER", threshold = TRUE, special.token + ######### + ## + ## set the arguments for calls + ## + ######### + ## Setting the continuous characters call + if(do_continuous) { + ## Create the character arguments + character_continuous_args <- mapply(function(character, ace.args) return(c(x = character, ace.args)), characters_continuous, models_continuous, SIMPLIFY = FALSE) - - - - - ## Setting up discrete characters - tree_args_list_discrete <- tree_args_list_continous <- NULL - - ## Running discrete characters + ## Create the character and tree arguments + tree_character_continuous_args <- list() + for(one_tree in 1:length(tree)) { + tree_character_continuous_args[[one_tree]] <- lapply(character_continuous_args, function(character, tree) {character$phy <- tree; return(character)}, tree[[one_tree]]) + } + } + ## Setting the discrete characters call if(do_discrete) { - - ## Set up the arguments for one tree - args_list <- mapply(make.args, characters_tables, characters_states, models, - MoreArgs = list(castor.options, cores, estimation.details), SIMPLIFY = FALSE) - - ## Add up the tree arguments - add.tree <- function(tree, args_list) { - return(lapply(args_list, function(arg, tree) c(arg, tree = list(tree)), tree)) - } - tree_args_list <- lapply(tree, add.tree, args_list) - - if(verbose) cat("Done.\n") + character_discrete_args <- mapply(make.args, characters_tables, characters_states, models_discrete, SIMPLIFY = FALSE) - if(do_parallel) { - ## Remove verbose - if(verbose) { - cat(paste0("Running the estimation for ", length(tree), " tree", ifelse(length(tree) > 1, "s ", " "), "using ", cores, " core", ifelse(cores == 1, "...", "s..."))) - was_verbose <- TRUE - verbose <- FALSE - } else { - was_verbose <- FALSE - } - - ## Set up the cluster - cluster <- parallel::makeCluster(cores) - - ## Get the current environment - current_env <- environment() - ## Get the export lists - export_arguments_list <- c("tree_args_list", - "special.tokens", - "invariants", - "threshold.type", - "threshold", - "verbose", - "characters_states", - "invariant_characters_states") - export_functions_list <- c("one.tree.ace", - "castor.ace", - "tree.data.update", - "add.state.names", - "translate.likelihood") - - ## Export from this environment - parallel::clusterExport(cluster, c(export_arguments_list, export_functions_list), envir = current_env) - - ## Call the cluster - results_out <- parLapply(cl = cluster, tree_args_list, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) - - ## Stop the cluster - parallel::stopCluster(cluster) - - ## Reactivate the verbose - if(was_verbose) { - cat("Done.") - verbose <- TRUE - } - } else { - ## Running the ACE for each tree - results_out <- lapply(tree_args_list, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) + ## Create the character and tree arguments + tree_character_discrete_args <- list() + for(one_tree in 1:length(tree)) { + tree_character_discrete_args[[one_tree]] <- lapply(character_discrete_args, function(character, tree) {character$tree <- tree; return(character)}, tree[[one_tree]]) } - ## Separating results and details - details_out_discrete <- lapply(results_out, `[[`, 2) - results_out_discrete <- lapply(results_out, `[[`, 1) } - ## Running continuous characters - if(do_continous) { - - models - tree - data - #TODO: set up the code for continuous characters (model and shit) - warning("TODO: multi.ace") - - - ## Set the model list (should be handled above - when selecting the models) - character_models <- list("char1" = list(method = "REML", - model = "BM", - scaled = TRUE), - "char2" = list(method = "REML", - model = "BM", - scaled = TRUE), - "char3" = list(method = "REML", - model = "BM", - scaled = TRUE), - "char4" = list(method = "REML", - model = "BM", - scaled = TRUE), - "char5" = list(method = "REML", - model = "BM", - scaled = TRUE)) + ######### + ## + ## run the calls + ## + ######### - ## Set the characters as a list (should be handled above - when splitting discrete characters) - characters <- apply(data, 2, list) + if(do_parallel) { + ## Do parallel! + # if(do_parallel) { + # ## Remove verbose + # if(verbose) { + # cat(paste0("Running the estimation for ", length(tree), " tree", ifelse(length(tree) > 1, "s ", " "), "using ", cores, " core", ifelse(cores == 1, "...", "s..."))) + # was_verbose <- TRUE + # verbose <- FALSE + # } else { + # was_verbose <- FALSE + # } + + # ## Set up the cluster + # cluster <- parallel::makeCluster(cores) + + # ## Get the current environment + # current_env <- environment() + # ## Get the export lists + # export_arguments_list <- c("tree_args_list", + # "special.tokens", + # "invariants", + # "threshold.type", + # "threshold", + # "verbose", + # "characters_states", + # "invariant_characters_states") + # export_functions_list <- c("one.tree.ace", + # "castor.ace", + # "tree.data.update", + # "add.state.names", + # "translate.likelihood") + + # ## Export from this environment + # parallel::clusterExport(cluster, c(export_arguments_list, export_functions_list), envir = current_env) + + # ## Call the cluster + # results_out <- parLapply(cl = cluster, tree_args_list, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) + + # ## Stop the cluster + # parallel::stopCluster(cluster) + + # ## Reactivate the verbose + # if(was_verbose) { + # cat("Done.") + # verbose <- TRUE + # } + # } - ## Create the character arguments - character_args <- mapply(function(character, ace.args) return(c(x = character, ace.args)), characters, character_models, SIMPLIFY = FALSE) - ## Create the character and tree arguments - tree_character_args <- list() - for(one_tree in 1:length(tree)) { - tree_character_args[[one_tree]] <- lapply(character_args, function(character, tree) {character$phy <- tree; return(character)}, tree[[one_tree]]) + } else { + ## Make the functions verbose + if(verbose) { + fun_discrete <- function(...) { + cat(".") + return(castor::asr_mk_model(...)) + } + fun_continuous <- function(...) { + cat(".") + return(ape::ace(...)) + } + } else { + fun_discrete <- castor::asr_mk_model + fun_continuous <- ape::ace } - ## Run all the ace - raw_estimates <- lapply(tree_character_args, lapply, function(x) do.call(what = ape::ace, args = x)) - ## Remove the ugly call - raw_estimates <- lapply(raw_estimates, lapply, function(x) {x$call <- "ape::ace"; return(x)}) - + ## Run the continuous characters + if(do_continuous) { + ## Run all the ace + continuous_estimates <- lapply(tree_character_continuous_args, lapply, function(x) do.call(fun_continuous, x)) + ## Remove the ugly call + continuous_estimates <- lapply(continuous_estimates, lapply, function(x) {x$call <- "ape::ace"; return(x)}) + ## Add node labels + # TODO: to be removed if ape update + add.nodes <- function(obj, phy) { + ## Adding node labels to $ace and $CI95 (if available) + options(warn = -1) + names <- as.integer(names(obj$ace)) + options(warn = 0) + if(!is.null(phy$node.label) && !is.na(names[1])) { + ordered_node_labels <- phy$node.label[c(as.integer(names(obj$ace))- Ntip(phy))] + names(obj$ace) <- ordered_node_labels + if(!is.null(obj$CI95)) { + rownames(obj$CI95) <- ordered_node_labels + } + } + return(obj) + } + for(one_tree in 1:length(tree)) { + continuous_estimates[[one_tree]] <- lapply(continuous_estimates[[one_tree]], add.nodes, phy = tree[[one_tree]]) + } + } + ## Run the discrete characters + if(do_discrete) { + ## Run all the ace + discrete_estimates <- lapply(tree_character_discrete_args, lapply, function(x) do.call(fun_discrete, x)) + } + if(verbose) cat("Done.\n") } + ######### + ## + ## handle the outputs + ## + ######### - ## Standardise the whole procedure to: - #1- characters are sorted in two categories - #2- characters in each category get their models sorted - #3- characters from both categories get their trees attributed - #4- characters and tree are run with do.call - - estimates_continuous <- estimates_discrete <- NULL - if(do_continous) { - ## Run all the ace - estimates_continuous <- lapply(tree_character_args_continuous, lapply, function(x) do.call(what = ape::ace, args = x)) - ## Remove the ugly call - estimates_continuous <- lapply(estimates_continuous, lapply, function(x) {x$call <- "ape::ace"; return(x)}) - } - if(do_discrete) { - ## Run all the castor - estimates_discrete <- lapply(tree_character_args_continuous, lapply, function(x) do.call(what = castor::asr_mk_model, args = x)) - ## Remove the ugly call - estimates_discrete <- lapply(estimates_discrete, lapply, function(x) {x$call <- "ape::ace"; return(x)}) + ## Handle the continuous characters + if(do_continuous) { + ## Get the results in a matrix format + results_continuous <- lapply(lapply(continuous_estimates, lapply, `[[`, "ace"), function(x) do.call(cbind, x)) + ## Get the details for continuous + details_continuous <- NULL + #TODO: add details for continuous } + if(do_discrete) { + ## Get the results in a matrix format + results_discrete <- NULL + ## Apply the threshold method + ## Add invariant characters + ## Get the details + } - - - - - - - - - ## Reorder the results - #TODO: reorder the results using continuous_char_ID and discrete_char_ID: like - details_out <- details_out[c(continuous_char_ID, discrete_char_ID)] - results_out <- results_out[c(continuous_char_ID, discrete_char_ID)] + ## Reorder the output results and details + # details_out <- details_out[c(continuous_char_ID, discrete_char_ID)] + # results_out <- results_out[c(continuous_char_ID, discrete_char_ID)] warning("TODO: multi.ace") ## Output a matrix diff --git a/R/multi.ace_fun.R b/R/multi.ace_fun.R index 89f8bd1a..537d9e6a 100755 --- a/R/multi.ace_fun.R +++ b/R/multi.ace_fun.R @@ -71,17 +71,23 @@ convert.char.table <- function(character, character_states) { } ## Set up the characters arguments for one tree -make.args <- function(character, character_states, model, castor.options, cores, estimation.details) { +make.args <- function(character, character_states, model, castor.options = NULL, cores = NULL, estimation.details = NULL) { ## Get the list of arguments castor_args <- list(tip_states = NULL, Nstates = length(character_states), rate_model = model, tip_priors = character, check_input = FALSE, - Ntrials = 1, - Nthreads = cores, - details = estimation.details) + Ntrials = 1) ## Add options + if(!is.null(cores)) { + castor_args$Nthreads <- cores + } else { + castor_args$Nthreads <- 1 + } + if(!is.null(estimation.details)) { + castor_args$details <- estimation.details + } if(!is.null(castor.options)) { castor.options <- c(castor_args, castor.options) } diff --git a/man/multi.ace.Rd b/man/multi.ace.Rd index 3f924611..6ec5f043 100755 --- a/man/multi.ace.Rd +++ b/man/multi.ace.Rd @@ -20,7 +20,7 @@ multi.ace( ) } \arguments{ -\item{data}{A \code{matrix} or \code{list} with the characters for each taxa.} +\item{data}{A \code{matrix}, \code{data.frame} or \code{list} with the characters for each taxa.} \item{tree}{A \code{phylo} or \code{mutiPhylo} object (if the \code{tree} argument contains node labels, they will be used to name the output).} @@ -47,13 +47,22 @@ multi.ace( \value{ Returns a \code{"matrix"} or \code{"list"} of ancestral states. By default, the function returns the ancestral states in the same format as the input \code{matrix}. This can be changed using the option \code{output = "matrix"} or \code{"list"} to force the class of the output. To output the combined ancestral states and input, you can use \code{"combined"} (using the input format) or \code{"combined.matrix"} or \code{"combined.list"}. +\emph{NOTE} that if the input data had multiple character types (continuous and discrete) and that \code{"matrix"} or \code{"combined.matrix"} output is requested, the function returns a \code{"data.frame"}. } \description{ Fast ancestral states estimations run on multiple trees using the Mk model from castor::asr_mk_model. } \details{ -The \code{models} argument can be a single or a list of transition \code{matrix}, a single or a a vector of built-in model(s) (see below) or a list of both matrices and built-in models: -The available built-in models in \code{castor::asr_mk_model} are: +Depending on the type of characters \code{models} argument can be either: +\itemize{ + \item the name of a single model to apply to all characters (if all characters are discrete or all are continuous); see below for the list of available names. For example \code{models = "ER"} applies the Equal Rates model to all characters (assuming they are all discrete characters). + \item a vector of model names to apply to different type of characters (see below for the list). For example \code{models = c("ER", "ER", "BM")} applies the Equal Rates model to the two first characters (discrete) and the \code{"BM"} model to the third character (continuous). + \item a transition \code{"matrix"} to be applied to all characters (if discrete). For example \code{models = matrix(0.2, 2, 2)}. + \item an single named list of arguments to be applied to all characters by passing it to \code{ape::ace} (if continuous). For example \code{models = list(method = "GLS", corStruct = corBrownian(1, my_tree))}. + \item an un-ambiguous list of arguments to be passed to either \code{castor::asr_mk_model} (discrete characters) or \code{ape::ace} (continuous characters). For example \code{models = list("char1" = list(transition_matrix = matrix(0.2, 2, 2)), "char2" = list(method = "GLS", corStruct = corBrownian(1, my_tree)))} to be specifically passed to the characters named "char1" and "char2" +} + +The available built-in models for discrete characters in \code{castor::asr_mk_model} are: \itemize{ \item \code{"ER"} for all equal rates \item \code{"SYM"} for symmetric rates @@ -63,6 +72,14 @@ The available built-in models in \code{castor::asr_mk_model} are: } See directly \code{castor::asr_mk_model} for more models. +The available built-in models and methods for continuous characters in \code{ape::ace} are: +\itemize{ + \item \code{"BM"} model: for a default Brownian Motion with the "REML" method + \item \code{"REML"} method: for a default Brownian Motion with the "REML" method (same as above) + \item \code{"ML"} method: for a default Brownian Motion with the "ML" method + \item \code{"pic"} method: for a default Brownian Motion with the "pic" (least squared) method +} + The \code{threshold} option allows to convert ancestral states likelihoods into discrete states. When \code{threshold = FALSE}, the ancestral state estimated is the one with the highest likelihood (or at random if likelihoods are equal). When \code{threshold = TRUE}, the ancestral state estimated are all the ones that are have a scaled likelihood greater than the maximum observed scaled likelihood minus the inverse number of possible states (i.e. \code{select_state >= (max(likelihood) - 1/n_states)}). This option makes the threshold selection depend on the number of states (i.e. if there are more possible states, a lower scaled likelihood for the best state is expected). Finally using a numerical value for the threshold option (e.g. \code{threshold = 0.95}) will simply select only the ancestral states estimates with a scaled likelihood equal or greater than the designated value. This option makes the threshold selection absolute. Regardless, if more than one value is select, the uncertainty token (\code{special.tokens["uncertainty"]}) will be used to separate the states. If no value is selected, the uncertainty token will be use between all observed characters (\code{special.tokens["uncertainty"]}). \code{special.behaviours} allows to generate a special rule for the \code{special.tokens}. The functions should can take the arguments \code{character, all_states} with \code{character} being the character that contains the special token and \code{all_states} for the character (which is automatically detected by the function). By default, missing data returns and inapplicable returns all states, and polymorphisms and uncertainties return all present states. diff --git a/tests/testthat/test-multi.ace.R b/tests/testthat/test-multi.ace.R index a4d0189d..315dbc73 100755 --- a/tests/testthat/test-multi.ace.R +++ b/tests/testthat/test-multi.ace.R @@ -414,7 +414,7 @@ test_that("multi.ace works", { }) test_that("multi.ace works with continuous and mix", { - + set.seed(1) ## The tree tree <- rcoal(15) tree <- makeNodeLabel(tree) From 25c0481d720bbbb3f7f91a9d0307bf888515634b Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Mon, 22 Jan 2024 14:05:24 +0000 Subject: [PATCH 05/15] Handled castor output and combined data --- R/multi.ace.R | 58 +++++++++++++++++++++++++++++++---------------- R/multi.ace_fun.R | 17 +++++++++++--- 2 files changed, 52 insertions(+), 23 deletions(-) diff --git a/R/multi.ace.R b/R/multi.ace.R index 7f7f94db..e455eb9d 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -656,22 +656,21 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec } else { ## Make the functions verbose - if(verbose) { - fun_discrete <- function(...) { - cat(".") - return(castor::asr_mk_model(...)) - } - fun_continuous <- function(...) { - cat(".") - return(ape::ace(...)) - } - } else { - fun_discrete <- castor::asr_mk_model - fun_continuous <- ape::ace - } + cat("Running ancestral states estimations:\n") ## Run the continuous characters if(do_continuous) { + + ## Set verbose fun + if(verbose) { + fun_continuous <- function(...) { + cat(".") + return(ape::ace(...)) + } + } else { + fun_continuous <- ape::ace + } + ## Run all the ace continuous_estimates <- lapply(tree_character_continuous_args, lapply, function(x) do.call(fun_continuous, x)) ## Remove the ugly call @@ -698,8 +697,8 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec } ## Run the discrete characters if(do_discrete) { - ## Run all the ace - discrete_estimates <- lapply(tree_character_discrete_args, lapply, function(x) do.call(fun_discrete, x)) + ## Run all the ace for discrete + discrete_estimates <- lapply(tree_character_discrete_args, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) } if(verbose) cat("Done.\n") } @@ -721,20 +720,35 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec if(do_discrete) { ## Get the results in a matrix format - results_discrete <- NULL - - ## Apply the threshold method - - ## Add invariant characters + results_discrete <- lapply(lapply(discrete_estimates, `[[`, 1), function(x) do.call(cbind, x)) ## Get the details + details_discrete <- lapply(discrete_estimates, `[[`, 2) + #TODO: add details for discrete + } + ## Reorder the output results and details # details_out <- details_out[c(continuous_char_ID, discrete_char_ID)] # results_out <- results_out[c(continuous_char_ID, discrete_char_ID)] warning("TODO: multi.ace") + + + ## Handle output + ## Combine the traits + if(do_discrete && do_continuous) { + ## Combine the traits + results_out <- mapply(bind.characters, results_continuous, results_discrete, + MoreArgs = list(order = list("continuous" = continuous_char_ID, "discrete" = discrete_char_ID)), + SIMPLIFY = FALSE) + } + + ## Choose the output type: + warning("TODO") + + ## Output a matrix make.matrix <- function(results) { return(lapply(results, function(data) do.call(cbind, data))) @@ -765,6 +779,10 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec #TODO: make sure node.labels available in output warning("DEBUG: multi.ace") + + ## Handle details + + ## Results out if(is.null(estimation.details)) { if(length(tree) == 1) { diff --git a/R/multi.ace_fun.R b/R/multi.ace_fun.R index 537d9e6a..c7c663d2 100755 --- a/R/multi.ace_fun.R +++ b/R/multi.ace_fun.R @@ -215,7 +215,7 @@ translate.likelihood <- function(character, threshold, select.states, special.to one.tree.ace <- function(args_list, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) { if(verbose) body(castor.ace)[[2]] <- substitute(cat(".")) - if(verbose) cat("Running ancestral states estimations:\n") + # if(verbose) cat("Running ancestral states estimations:\n") ancestral_estimations <- lapply(args_list, castor.ace) ancestral_estimations <- mapply(add.state.names, ancestral_estimations, characters_states, SIMPLIFY = FALSE) @@ -269,7 +269,7 @@ one.tree.ace <- function(args_list, special.tokens, invariants, characters_state ## Replace NAs replace.NA <- function(character, characters_states, special.tokens) { - return(unname(sapply(character, function(x) ifelse(x[[1]] == "NA", paste0(characters_states, collapse = sub("\\\\", "", special.tokens["uncertainty"])), x)))) + return(sapply(character, function(x) ifelse(x[[1]] == "NA", paste0(characters_states, collapse = sub("\\\\", "", special.tokens["uncertainty"])), x))) } ancestral_states <- mapply(replace.NA, ancestral_states, characters_states, MoreArgs = list(special.tokens = special.tokens), SIMPLIFY = FALSE) @@ -286,6 +286,17 @@ one.tree.ace <- function(args_list, special.tokens, invariants, characters_state estimations_details_out <- NULL } - if(verbose) cat(" Done.\n") + # if(verbose) cat(" Done.\n") return(list(results = ancestral_states, details = estimations_details_out)) } + +## Bind the continuous and discrete characters and reorder them +bind.characters <- function(continuous, discrete, order) { + bound <- cbind(as.data.frame(continuous), as.data.frame(discrete)) + ## Get the new character IDs + cont_names <- colnames(continuous) + disc_names <- colnames(discrete) + ## Reorder the characters to match the input order + ordering <- matrix(c(1:ncol(bound), c(order$continuous, order$discrete)), ncol = 2, byrow = FALSE, dimnames = list(c(cont_names, disc_names), c("out", "in"))) + return(bound[, names(sort(ordering[, 2, drop = TRUE]))]) +} \ No newline at end of file From d506641a1eed750359c22329f5c9612bc576fd57 Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Mon, 22 Jan 2024 16:08:06 +0000 Subject: [PATCH 06/15] Reorder arguments - ready for testing --- R/multi.ace.R | 108 +++++++++++++++------------------------------- R/multi.ace_fun.R | 15 ++++++- 2 files changed, 48 insertions(+), 75 deletions(-) diff --git a/R/multi.ace.R b/R/multi.ace.R index e455eb9d..eadb295b 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -189,10 +189,6 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec class(tree) <- "multiPhylo" node_labels <- lapply(node_labels, `[[`, 2) - - - - ######### ## ## Handle the other options (threshold, brlen, verbose, parallel, output, estimation.details) @@ -263,26 +259,6 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec cores <- parallel } - ## output - if(missing(output)) { - output <- class(matrix)[1] - } else { - check.class(output, "character") - available_methods <- c("matrix", "list", "combined", "combined.list", "combined.matrix")#, "dispRity") - check.method(output, available_methods, "output option") - ## Combined - if(output == "combined") { - output <- paste(output, class(matrix)[1], sep = ".") - } - } - - ## Check the estimation details - if(!is.null(estimation.details)) { - ## The return args from castor::asr_mk_model (1.6.6) - return_args <- c("success", "Nstates", "transition_matrix", "loglikelihood", "ancestral_likelihoods") - check.method(estimation.details, return_args, msg = "estimation.details") - } - ######### ## Handle the characters ######### @@ -318,6 +294,23 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec discrete_char_ID <- which(!character_is_continuous) } + ## output + if(missing(output)) { + output <- class(matrix)[1] + } else { + check.class(output, "character") + available_methods <- c("matrix", "list", "combined", "combined.list", "combined.matrix", "dispRity") + check.method(output, available_methods, "output option") + ## Combined + if(output == "combined") { + output <- paste(output, class(matrix)[1], sep = ".") + } + ## Check dispRity + if(output == "dispRity" && do_discrete) { + stop("Only ancestral state estimations for continuous characters can be converted into a dispRity object.\nSelect an other output method.") + } + } + ## Handle the tokens # special.tokens <- character(); special.behaviours <- list() ; warning("DEBUG: multi.ace") if(do_discrete) { @@ -534,9 +527,6 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec } } - - - ######### ## ## Handle the options @@ -544,10 +534,6 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ######### - - - - ## options.args warning("TODO: change castor options to just options (parsed to castor or ape)") warning("TODO: multi.ace. ape options are: CI, ") @@ -564,9 +550,12 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec } } - - - + ## Check the estimation details + if(!is.null(estimation.details)) { + ## The return args from castor::asr_mk_model (1.6.6) + return_args <- c("success", "Nstates", "transition_matrix", "loglikelihood", "ancestral_likelihoods") + check.method(estimation.details, return_args, msg = "estimation.details") + } ######### ## @@ -597,8 +586,6 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec } } - - ######### ## ## run the calls @@ -606,6 +593,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ######### if(do_parallel) { + stop("TODO: do parallel multi.ace") ## Do parallel! # if(do_parallel) { # ## Remove verbose @@ -725,17 +713,8 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Get the details details_discrete <- lapply(discrete_estimates, `[[`, 2) #TODO: add details for discrete - } - - ## Reorder the output results and details - # details_out <- details_out[c(continuous_char_ID, discrete_char_ID)] - # results_out <- results_out[c(continuous_char_ID, discrete_char_ID)] - warning("TODO: multi.ace") - - - ## Handle output ## Combine the traits if(do_discrete && do_continuous) { @@ -744,43 +723,24 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec MoreArgs = list(order = list("continuous" = continuous_char_ID, "discrete" = discrete_char_ID)), SIMPLIFY = FALSE) } - - ## Choose the output type: - warning("TODO") - - - ## Output a matrix - make.matrix <- function(results) { - return(lapply(results, function(data) do.call(cbind, data))) - } - ## Combine the ace matrix with the tips - add.tips <- function(ace, matrix) { - return(rbind(matrix, ace)) + if(do_discrete && !do_continuous) { + results_out <- results_discrete } - ## Output a list from a matrix - make.list <- function(results) { - ## Make into a list - return(unlist(apply(results, 1, list), recursive = FALSE)) + if(do_continuous && !do_discrete) { + results_out <- results_continuous } - - ## Make the basic output matrix - output_matrix <- make.matrix(results_out) ## Handle output output_return <- switch(output, - matrix = output_matrix, - list = lapply(output_matrix, make.list), - combined.matrix = lapply(output_matrix, add.tips, matrix = matrix), - combined.list = lapply(lapply(output_matrix, add.tips, matrix = matrix), make.list)# - #dispRity = return(list("tips" = matrix, "nodes" = output_matrix)) + matrix = results_out, + list = lapply(results_out, make.list), + combined.matrix = lapply(results_out, add.tips, matrix = matrix), + combined.list = lapply(lapply(results_out, add.tips, matrix = matrix), make.list), + dispRity = make.dispRity(data = lapply(results_out, add.tips, matrix = matrix), tree = tree) ) - #TODO: set up "dispRity" output - #TODO: make sure node.labels available in output - warning("DEBUG: multi.ace") - - ## Handle details + warning("TODO: handle details") ## Results out diff --git a/R/multi.ace_fun.R b/R/multi.ace_fun.R index c7c663d2..dd2c130c 100755 --- a/R/multi.ace_fun.R +++ b/R/multi.ace_fun.R @@ -299,4 +299,17 @@ bind.characters <- function(continuous, discrete, order) { ## Reorder the characters to match the input order ordering <- matrix(c(1:ncol(bound), c(order$continuous, order$discrete)), ncol = 2, byrow = FALSE, dimnames = list(c(cont_names, disc_names), c("out", "in"))) return(bound[, names(sort(ordering[, 2, drop = TRUE]))]) -} \ No newline at end of file +} + + + +## Combine the ace matrix with the tips +add.tips <- function(ace, matrix) { + return(rbind(matrix, ace)) +} +## Output a list from a matrix +make.list <- function(results) { + ## Make into a list + return(unlist(apply(results, 1, list), recursive = FALSE)) +} +## Make the dispRity object out of the results \ No newline at end of file From eb781bed73486f6f2c208df6315190ffb2ac073b Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Mon, 22 Jan 2024 16:27:06 +0000 Subject: [PATCH 07/15] multi works --- R/multi.ace.R | 9 +++------ R/multi.ace_fun.R | 2 +- man/multi.ace.Rd | 6 +++--- tests/testthat/test-multi.ace.R | 11 ++++++++--- 4 files changed, 15 insertions(+), 13 deletions(-) diff --git a/R/multi.ace.R b/R/multi.ace.R index eadb295b..ae87a093 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -209,7 +209,6 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Use an absolute threshold threshold.type <- "absolute" } - ## brlen multiplier if(!missing(brlen.multiplier)) { @@ -421,6 +420,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec characters_continuous <- apply(matrix_continuous, 2, list) if(verbose) cat(".") } + if(verbose) cat("Done.\n") ######### ## Handle the models for each character @@ -536,7 +536,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## options.args warning("TODO: change castor options to just options (parsed to castor or ape)") - warning("TODO: multi.ace. ape options are: CI, ") + warning("TODO: multi.ace. ape options are: CI, ") if(missing(options.args)) { @@ -640,11 +640,9 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec # verbose <- TRUE # } # } - - } else { ## Make the functions verbose - cat("Running ancestral states estimations:\n") + if(verbose) cat("Running ancestral states estimations:\n") ## Run the continuous characters if(do_continuous) { @@ -742,7 +740,6 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Handle details warning("TODO: handle details") - ## Results out if(is.null(estimation.details)) { if(length(tree) == 1) { diff --git a/R/multi.ace_fun.R b/R/multi.ace_fun.R index dd2c130c..a35029e9 100755 --- a/R/multi.ace_fun.R +++ b/R/multi.ace_fun.R @@ -312,4 +312,4 @@ make.list <- function(results) { ## Make into a list return(unlist(apply(results, 1, list), recursive = FALSE)) } -## Make the dispRity object out of the results \ No newline at end of file +## Make the dispRity object out of the results (only for continuous) diff --git a/man/multi.ace.Rd b/man/multi.ace.Rd index 6ec5f043..c58e56a0 100755 --- a/man/multi.ace.Rd +++ b/man/multi.ace.Rd @@ -7,7 +7,7 @@ multi.ace( data, tree, - models = "ER", + models, threshold = TRUE, special.tokens, special.behaviours, @@ -15,7 +15,7 @@ multi.ace( verbose = FALSE, parallel = FALSE, output, - castor.options, + options.args, estimation.details = NULL ) } @@ -40,7 +40,7 @@ multi.ace( \item{output}{optional, see Value section below.} -\item{castor.options}{optional, a named list of options to be passed to function called by \code{castor::asr_mk_model}.} +\item{options.args}{optional, a named list of options to be passed to function called by \code{castor::asr_mk_model}.} \item{estimation.details}{optional, whether to also return the details for each estimation as returned by \code{castor::asr_mk_model}. This argument can be left \code{NULL} (default) or be any combination of the elements returned by \code{castor::asr_mk_model} (e.g. \code{c("loglikelihood", "transition_matrix")}).} } diff --git a/tests/testthat/test-multi.ace.R b/tests/testthat/test-multi.ace.R index 315dbc73..89b1f9b1 100755 --- a/tests/testthat/test-multi.ace.R +++ b/tests/testthat/test-multi.ace.R @@ -422,7 +422,7 @@ test_that("multi.ace works with continuous and mix", { data <- space.maker(elements = 15, dimensions = 5, distribution = rnorm, elements.name = tree$tip.label) ## Run the multi.ace on the continuous data - test <- multi.ace(data = data, tree = tree, output = "combined.matrix") + test <- multi.ace(data = data, tree = tree, output = "combined.matrix", verbose = TRUE) ## Works well for continuous expect_is(test, "matrix") @@ -438,7 +438,12 @@ test_that("multi.ace works with continuous and mix", { ## Works well for mixed characters test <- multi.ace(data = data, tree = tree, output = "combined.matrix") expect_is(test, "data.frame") - expect_equal(dim(test), c(15+14, 6)) + expect_equal(dim(test), c(15+14, 7)) expect_equal(sort(rownames(test)), sort(c(tree$tip.label, tree$node.label))) - expect_equal(unique(apply(test, 2, class)), c("numeric", "character")) + + classes <- character() + for(i in 1:ncol(test)) { + classes[i] <- class(test[, i]) + } + expect_equal(unique(classes), c("numeric", "character")) }) \ No newline at end of file From 020618d08bafc6d42a9a9690e6c761306f094707 Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Tue, 23 Jan 2024 12:37:10 +0000 Subject: [PATCH 08/15] sorted details handeling missing options and then done --- R/multi.ace.R | 56 +++++++++++++++++++++++++++------ R/multi.ace_fun.R | 16 ++++++++++ tests/testthat/test-multi.ace.R | 6 ++-- 3 files changed, 66 insertions(+), 12 deletions(-) diff --git a/R/multi.ace.R b/R/multi.ace.R index ae87a093..8e31c090 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -18,7 +18,7 @@ #' @param parallel \code{logical}, whether to use parallel algorithm (\code{TRUE}) or not (\code{FALSE} - default). #' @param output optional, see Value section below. #' @param options.args optional, a named list of options to be passed to function called by \code{castor::asr_mk_model}. -#' @param estimation.details optional, whether to also return the details for each estimation as returned by \code{castor::asr_mk_model}. This argument can be left \code{NULL} (default) or be any combination of the elements returned by \code{castor::asr_mk_model} (e.g. \code{c("loglikelihood", "transition_matrix")}). +#' @param estimation.details optional, whether to also return the details for each estimation as returned by \code{castor::asr_mk_model} or \code{ape::ace}. This argument can be left \code{NULL} (default) or be any combination of the elements returned by \code{castor::asr_mk_model} or \code{ape::ace} (e.g. \code{c("loglikelihood", "transition_matrix", "CI95")}). #' #' @details #' @@ -533,11 +533,8 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## ######### - ## options.args warning("TODO: change castor options to just options (parsed to castor or ape)") - warning("TODO: multi.ace. ape options are: CI, ") - if(missing(options.args)) { ## No options @@ -553,7 +550,19 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Check the estimation details if(!is.null(estimation.details)) { ## The return args from castor::asr_mk_model (1.6.6) - return_args <- c("success", "Nstates", "transition_matrix", "loglikelihood", "ancestral_likelihoods") + return_args_discrete <- c("success", "Nstates", "transition_matrix", "loglikelihood", "ancestral_likelihoods") + return_args_continuous <- c("CI95", "sigma2", "loglik") + if(do_discrete && do_continuous) { + return_args <- c(return_args_discrete, return_args_continuous) + } else { + if(do_discrete) { + return_args <- return_args_discrete + } + if(do_continuous) { + return_args <- return_args_continuous + } + } + ## Check the requested details check.method(estimation.details, return_args, msg = "estimation.details") } @@ -576,8 +585,16 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec } ## Setting the discrete characters call if(do_discrete) { + + ## Set the details to return (if any) + if(any(return_args_discrete %in% estimation.details)) { + details_out <- return_args_discrete[return_args_discrete %in% estimation.details] + } else { + details_out <- NULL + } + ## Set up the arguments for one tree - character_discrete_args <- mapply(make.args, characters_tables, characters_states, models_discrete, SIMPLIFY = FALSE) + character_discrete_args <- mapply(make.args, characters_tables, characters_states, models_discrete, MoreArgs = list(estimation.details = details_out), SIMPLIFY = FALSE) ## Create the character and tree arguments tree_character_discrete_args <- list() @@ -699,9 +716,17 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec if(do_continuous) { ## Get the results in a matrix format results_continuous <- lapply(lapply(continuous_estimates, lapply, `[[`, "ace"), function(x) do.call(cbind, x)) + ## Get the details for continuous - details_continuous <- NULL - #TODO: add details for continuous + if(any(return_args_continuous %in% estimation.details)) { + ## Get which details to grep + details_out <- return_args_continuous[return_args_continuous %in% estimation.details] + + ## Get the details + details_continuous <- lapply(continuous_estimates, lapply, function(x, details_out) return(x[details_out]), details_out) + } else { + details_continuous <- NULL + } } if(do_discrete) { @@ -710,7 +735,6 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Get the details details_discrete <- lapply(discrete_estimates, `[[`, 2) - #TODO: add details for discrete } ## Handle output @@ -720,12 +744,26 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec results_out <- mapply(bind.characters, results_continuous, results_discrete, MoreArgs = list(order = list("continuous" = continuous_char_ID, "discrete" = discrete_char_ID)), SIMPLIFY = FALSE) + ## Return the details per characters + if(is.null(details_continuous)) { + ## Make a list of nulls + details_continuous <- replicate(length(tree), lapply(as.list(1:n_characters_continuous), function(x) return(NULL)), simplify = FALSE) + } + if(is.null(details_discrete[[1]])) { + ## Make a list of nulls + details_discrete <- replicate(length(tree), lapply(as.list(1:n_characters_discrete), function(x) return(NULL)), simplify = FALSE) + } + details_out <- mapply(bind.details, details_continuous, details_discrete, + MoreArgs = list(order = list("continuous" = continuous_char_ID, "discrete" = discrete_char_ID)), + SIMPLIFY = FALSE) } if(do_discrete && !do_continuous) { results_out <- results_discrete + details_out <- details_discrete } if(do_continuous && !do_discrete) { results_out <- results_continuous + details_out <- details_continuous } ## Handle output diff --git a/R/multi.ace_fun.R b/R/multi.ace_fun.R index a35029e9..3228da14 100755 --- a/R/multi.ace_fun.R +++ b/R/multi.ace_fun.R @@ -300,6 +300,22 @@ bind.characters <- function(continuous, discrete, order) { ordering <- matrix(c(1:ncol(bound), c(order$continuous, order$discrete)), ncol = 2, byrow = FALSE, dimnames = list(c(cont_names, disc_names), c("out", "in"))) return(bound[, names(sort(ordering[, 2, drop = TRUE]))]) } +## Bind the continuous and discrete details and reorder them +bind.details <- function(continuous, discrete, order) { + ## Reorder the details out per characters + if(length(length(discrete[[1]])) > 1) { + discrete_details <- list() + for(one_char in 1:length(discrete[[1]])) { + discrete_details[[one_char]] <- lapply(discrete, `[[`, 1) + } + if(!is.null(discrete[[1]])) { + names(discrete_details) <- names(discrete[[1]]) + } + } + ## Bind the two lists + return(c(continuous, discrete_details)[c(order$continuous, order$discrete)]) +} + diff --git a/tests/testthat/test-multi.ace.R b/tests/testthat/test-multi.ace.R index 89b1f9b1..508789ca 100755 --- a/tests/testthat/test-multi.ace.R +++ b/tests/testthat/test-multi.ace.R @@ -34,7 +34,7 @@ test_that("multi.ace works", { verbose = FALSE, parallel = FALSE, output = "list")) - expect_equal(error[[1]], "matrix must be of class matrix or list.") + expect_equal(error[[1]], "matrix must be of class matrix or list or data.frame.") error <- capture_error(multi.ace(data = matrix_complex, tree = "tree_test", @@ -58,7 +58,7 @@ test_that("multi.ace works", { verbose = FALSE, parallel = FALSE, output = "list")) - expect_equal(error[[1]], "models must be of class character or matrix or list.") + expect_equal(error[[1]], "models must be of class character or list or matrix.") error <- capture_error(multi.ace(data = matrix_complex, tree = tree_test, @@ -142,7 +142,7 @@ test_that("multi.ace works", { verbose = FALSE, parallel = FALSE, output = "something")) - expect_equal(error[[1]], "output option must be one of the following: matrix, list, combined, combined.list, combined.matrix.") + expect_equal(error[[1]], "output option must be one of the following: matrix, list, combined, combined.list, combined.matrix, dispRity.") error <- capture_error(multi.ace(data = matrix_complex, From 16cf0b715afea10164ceb8fa8b36c51671b481fe Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Tue, 23 Jan 2024 13:07:54 +0000 Subject: [PATCH 09/15] changes stops to , call. = FALSE --- R/MCMCglmm.utilities.R | 5 ++--- R/char.diff.R | 10 +++++----- R/chrono.subsets.R | 2 +- R/chrono.subsets_fun.R | 2 +- R/dispRity.R | 2 +- R/dispRity.metric.R | 2 +- R/dispRity.utilities_fun.R | 14 +++++++------- R/dispRity_fun.R | 2 +- R/dtt.dispRity_fun.R | 18 +++++++++--------- R/make.metric_fun.R | 2 +- R/match.tip.edge.R | 2 +- R/multi.ace.R | 32 +++++++++++++++++--------------- R/pair.plot_fun.R | 2 +- R/pgls.dispRity.R | 8 ++++---- R/plot.dispRity_fun.R | 2 +- R/reduce.space.R | 4 ++-- R/remove.zero.brlen.R | 2 +- R/slice.tree_fun.R | 2 +- R/slide.nodes.R | 10 +++++----- R/test.metric.R | 2 +- 20 files changed, 63 insertions(+), 62 deletions(-) diff --git a/R/MCMCglmm.utilities.R b/R/MCMCglmm.utilities.R index 24eee3b4..ea1cb6b6 100755 --- a/R/MCMCglmm.utilities.R +++ b/R/MCMCglmm.utilities.R @@ -193,10 +193,9 @@ MCMCglmm.covars <- function(MCMCglmm, n, sample){ } else { check.class(sample, c("numeric", "integer")) - ## Check for incorect samples + ## Check for incorrect samples if(length(incorect_sample <- which(sample > length(MCMCglmm.sample(MCMCglmm)))) > 0) { - #dispRity_export in: MAKE dispRity STOP STYLE - stop("Some samples are not available in the MCMCglmm object.")#dispRity_export out: + stop("Some samples are not available in the MCMCglmm object.", call. = FALSE) } } } else { diff --git a/R/char.diff.R b/R/char.diff.R index 50eec2cf..65ca3960 100755 --- a/R/char.diff.R +++ b/R/char.diff.R @@ -124,7 +124,7 @@ char.diff <- function(matrix, method = "hamming", translate = TRUE, special.toke if(matrix_class == "list") { ## Check length if(length(matrix) != 2) { - stop(paste0("When matrix argument is a list, it must contain only two elements.\nYou can convert ", as.expression(match_call$matrix), " to a matrix using:\n", as.expression(match_call$matrix), " <- do.call(rbind, ", as.expression(match_call$matrix), ")")) + stop(paste0("When matrix argument is a list, it must contain only two elements.\nYou can convert ", as.expression(match_call$matrix), " to a matrix using:\n", as.expression(match_call$matrix), " <- do.call(rbind, ", as.expression(match_call$matrix), ")"), call. = FALSE) } ## Convert into a matrix @@ -141,7 +141,7 @@ char.diff <- function(matrix, method = "hamming", translate = TRUE, special.toke ## Checking for the reserved character reserved <- grep("\\@", matrix) if(length(reserved) > 0) { - stop("The matrix cannot contain the character '@' since it is reserved for the dispRity::char.diff function.") + stop("The matrix cannot contain the character '@' since it is reserved for the dispRity::char.diff function.", call. = FALSE) } ## Method is hamming by default @@ -174,12 +174,12 @@ char.diff <- function(matrix, method = "hamming", translate = TRUE, special.toke ## Checking for the reserved character reserved <- c("\\@", "@") %in% special.tokens if(any(reserved)) { - stop("special.tokens cannot contain the character '@' since it is reserved for the dispRity::char.diff function.") + stop("special.tokens cannot contain the character '@' since it is reserved for the dispRity::char.diff function.", call. = FALSE) } ## Checking whether the special.tokens are unique if(length(unique(special.tokens)) != length(special.tokens)) { - stop("special.tokens cannot contain duplicated tokens.") + stop("special.tokens cannot contain duplicated tokens.", call. = FALSE) } ## If any special token is NA, convert them as "N.A" temporarily @@ -221,7 +221,7 @@ char.diff <- function(matrix, method = "hamming", translate = TRUE, special.toke check.class(correction, "function") test_correction <- make.metric(correction, silent = TRUE)$type if(!is.null(test_correction) && test_correction == "error") { - stop("Incorrect correction function.") + stop("Incorrect correction function.", call. = FALSE) } } diff --git a/R/chrono.subsets.R b/R/chrono.subsets.R index 59d48283..2c74d2ff 100755 --- a/R/chrono.subsets.R +++ b/R/chrono.subsets.R @@ -445,7 +445,7 @@ chrono.subsets <- function(data, tree = NULL, method, time, model, inc.nodes = F check.class(bind.data, "logical") } else { if(bind.data) { - stop(paste0("Impossible to bind the data to the trees since the number of matrices (", length(data), ") is not equal to the number of trees (", length(tree), ").")) + stop(paste0("Impossible to bind the data to the trees since the number of matrices (", length(data), ") is not equal to the number of trees (", length(tree), ")."), call. = FALSE) } } diff --git a/R/chrono.subsets_fun.R b/R/chrono.subsets_fun.R index 9743509d..70cfb7d3 100755 --- a/R/chrono.subsets_fun.R +++ b/R/chrono.subsets_fun.R @@ -6,7 +6,7 @@ get.percent.age <- function(tree, percent = 0.01) { percent <- percent + 0.01 tree_slice <- slice.tree.sharp(tree, tree$root.time - (percent * tree$root.time)) if(percent >= 100) { - stop("Impossible to find a starting point to slice the tree. This can happen if the tree has no branch length or has a \"ladder\" structure. You can try to fix that by setting specific slicing times.") + stop("Impossible to find a starting point to slice the tree. This can happen if the tree has no branch length or has a \"ladder\" structure. You can try to fix that by setting specific slicing times.", call. = FALSE) break } } diff --git a/R/dispRity.R b/R/dispRity.R index 1cb2f02e..333968f1 100755 --- a/R/dispRity.R +++ b/R/dispRity.R @@ -342,7 +342,7 @@ dispRity <- function(data, metric, dimensions = NULL, ..., between.groups = FALS ## Serial is a list, check if it contains the right information (pairs of things that exist) pairs <- unique(unlist(lapply(between.groups, length))) if(length(pairs) > 1 || pairs != 2 || max(unlist(between.groups)) > length(data$subsets)) { - stop("The provided list of groups (between.groups) must be a list of pairs of subsets in the data.") + stop("The provided list of groups (between.groups) must be a list of pairs of subsets in the data.", call. = FALSE) } list_of_pairs <- between.groups is_between.groups <- TRUE diff --git a/R/dispRity.metric.R b/R/dispRity.metric.R index 6ec18d41..82460a28 100755 --- a/R/dispRity.metric.R +++ b/R/dispRity.metric.R @@ -1116,7 +1116,7 @@ projections.tree <- function(matrix, tree, type = c("root","ancestor"), referenc } ## Sanitizing (to avoid obscure error message!) if(any(is_null <- unlist(lapply(from_to, is.null)))) { - stop(paste0("The following type argument is not recognised in projections.tree: ", paste0(type[is_null], collapse = ", "))) + stop(paste0("The following type argument is not recognised in projections.tree: ", paste0(type[is_null], collapse = ", ")), call. = FALSE) } if(all(invariables)) { diff --git a/R/dispRity.utilities_fun.R b/R/dispRity.utilities_fun.R index bead0626..33d959ff 100755 --- a/R/dispRity.utilities_fun.R +++ b/R/dispRity.utilities_fun.R @@ -50,7 +50,7 @@ check.subsets <- function(subsets, data) { if(is(subsets, "numeric") || is(subsets, "integer")) { if(any(na_subsets <- is.na(match(subsets, 1:length(data$disparity))))) { ## Subsets not found - stop(paste0(ifelse(length(which(na_subsets)) > 1, "Subsets ", "Subset "), paste0(subsets[which(na_subsets)], collapse = ", "), " not found.")) + stop(paste0(ifelse(length(which(na_subsets)) > 1, "Subsets ", "Subset "), paste0(subsets[which(na_subsets)], collapse = ", "), " not found."), call. = FALSE) } } else { if(is(subsets, "character")) { @@ -61,10 +61,10 @@ check.subsets <- function(subsets, data) { ## Check if the searched ones exist if(any(na_subsets <- is.na(match(subset_search, subset_available)))) { ## Subsets not found - stop(paste0(ifelse(length(which(na_subsets)) > 1, "Subsets ", "Subset "), paste0(subsets[which(na_subsets)], collapse = ", "), " not found.")) + stop(paste0(ifelse(length(which(na_subsets)) > 1, "Subsets ", "Subset "), paste0(subsets[which(na_subsets)], collapse = ", "), " not found."), call. = FALSE) } } else { - stop("subsets argument must be of class \"numeric\" or \"character\".") + stop("subsets argument must be of class \"numeric\" or \"character\".", call. = FALSE) } } @@ -76,12 +76,12 @@ check.subsets <- function(subsets, data) { } if(length(subsets) > length(data$subsets)) { - stop("Not enough subsets in the original data.") + stop("Not enough subsets in the original data.", call. = FALSE) } else { if(is(subsets, "numeric") || is(subsets, "integer")) { if(any(na_subsets <- is.na(match(subsets, 1:length(data$subsets))))) { ## Subsets not found - stop(paste0(ifelse(length(which(na_subsets)) > 1, "Subsets ", "Subset "), paste0(subsets[which(na_subsets)], collapse = ", "), " not found.")) + stop(paste0(ifelse(length(which(na_subsets)) > 1, "Subsets ", "Subset "), paste0(subsets[which(na_subsets)], collapse = ", "), " not found."), call. = FALSE) } } else { if(is(subsets, "character")) { @@ -89,11 +89,11 @@ check.subsets <- function(subsets, data) { subsets <- subsets[which(is.na(match(subsets, names(data$subsets))))] orthograph <- ifelse(length(subsets) == 1, "Subset ", "Subsets ") - stop(paste0(orthograph, paste0(subsets, collapse = ", "), " not found.")) + stop(paste0(orthograph, paste0(subsets, collapse = ", "), " not found."), call. = FALSE) } } else { - stop("subsets argument must be of class \"numeric\" or \"character\".") + stop("subsets argument must be of class \"numeric\" or \"character\".", call. = FALSE) } } } diff --git a/R/dispRity_fun.R b/R/dispRity_fun.R index 96f0b174..a2c18499 100755 --- a/R/dispRity_fun.R +++ b/R/dispRity_fun.R @@ -268,7 +268,7 @@ decompose.VCV <- function(one_subsets_bootstrap, fun, data, use_array, use_tree #fun(data$covar[[one_subsets_bootstrap[1]]][[1]], data$covar[[one_subsets_bootstrap[2]]][[2]]) } } else { - stop("Impossible to use tree metric in dispRity with covar (yet!).") + stop("Impossible to use tree metric in dispRity with covar (yet!).", call. = FALSE) } } diff --git a/R/dtt.dispRity_fun.R b/R/dtt.dispRity_fun.R index 12a07935..163e2d38 100755 --- a/R/dtt.dispRity_fun.R +++ b/R/dtt.dispRity_fun.R @@ -66,7 +66,7 @@ geiger.dtt.dispRity <- function(phy, data, metric, relative){ ## By array if(length(dim(td$data)) != 3){ - stop("Error in data: must be a matrix or a array of matrix (length(dim(data)) must be equal to 2 or 3).") + stop("Error in data: must be a matrix or a array of matrix (length(dim(data)) must be equal to 2 or 3).", call. = FALSE) } ## Looping through the array @@ -161,7 +161,7 @@ geiger.sim.char <- function(phy, par, nsim = 1, model = c("BM", "speciational", nbranches<-nrow(phy$edge) nspecies<-Ntip(phy) - if(length(root)>1) stop("'root' should be a single value") + if(length(root)>1) stop("'root' should be a single value", call. = FALSE) if(model%in%c("BM", "speciational")) { @@ -190,7 +190,7 @@ geiger.sim.char <- function(phy, par, nsim = 1, model = c("BM", "speciational", } for(j in 1:nchar) { m=model.matrix[[j]] - if(!root%in%c(1:nrow(m))) stop(paste("'root' must be a character state from 1 to ", nrow(m), sep="")) + if(!root%in%c(1:nrow(m))) stop(paste("'root' must be a character state from 1 to ", nrow(m), sep=""), call. = FALSE) p=lapply(el, function(l) matexpo(m*l)) for(k in 1:nsim) { @@ -221,16 +221,16 @@ geiger.make.modelmatrix <- function(m, model=c("BM", "speciational", "discrete") for(j in 1:length(m)){ #.check.Qmatrix m=unique(dim(m[[j]])) - if(length(m)>1) stop("'Q' must be a square matrix") + if(length(m)>1) stop("'Q' must be a square matrix", call. = FALSE) didx=1 + 0L:(m - 1L) * (m + 1) - if(!all(abs(rowSums(m[[j]]))<0.000001)) stop("rows of 'Q' must sum to zero") - if(!all(m[[j]][didx]<=0)) stop("diagonal elements of 'Q' should be negative") - if(!all(m[[j]][-didx]>=0)) stop("off-diagonal elements of 'Q' should be positive") + if(!all(abs(rowSums(m[[j]]))<0.000001)) stop("rows of 'Q' must sum to zero", call. = FALSE) + if(!all(m[[j]][didx]<=0)) stop("diagonal elements of 'Q' should be negative", call. = FALSE) + if(!all(m[[j]][-didx]>=0)) stop("off-diagonal elements of 'Q' should be positive", call. = FALSE) } } } else { - if(is.numeric(m)) m=as.matrix(m) else stop("Supply 'm' as a matrix of rates") - if(any(diag(m)<0)) stop("'m' appears to have negative variance component(s)") + if(is.numeric(m)) m=as.matrix(m) else stop("Supply 'm' as a matrix of rates", call. = FALSE) + if(any(diag(m)<0)) stop("'m' appears to have negative variance component(s)", call. = FALSE) } return(m) } diff --git a/R/make.metric_fun.R b/R/make.metric_fun.R index a18c791b..5ec9fb86 100755 --- a/R/make.metric_fun.R +++ b/R/make.metric_fun.R @@ -11,6 +11,6 @@ check.metric <- function(metric) { return("class.metric") } } else { - stop("Invalid metric.") + stop("Invalid metric.", call. = FALSE) } } \ No newline at end of file diff --git a/R/match.tip.edge.R b/R/match.tip.edge.R index edfe60b6..30ee3538 100755 --- a/R/match.tip.edge.R +++ b/R/match.tip.edge.R @@ -58,7 +58,7 @@ match.tip.edge <- function(vector, phylo, replace.na, use.parsimony = TRUE) { ## TODO: check number of nodes as well if(length(vector) != Ntip(phylo)[1]) { if(length(vector) != Ntip(phylo)[1]+Nnode(phylo)[1]) { - stop(paste0("The input vector must of the same length as the number of tips (", Ntip(phylo)[1], ") or tips and nodes (", Ntip(phylo)[1]+Nnode(phylo)[1] ,") in phylo.")) + stop(paste0("The input vector must of the same length as the number of tips (", Ntip(phylo)[1], ") or tips and nodes (", Ntip(phylo)[1]+Nnode(phylo)[1] ,") in phylo."), call. = FALSE) } } if(length(vector) == Ntip(phylo)[1]+Nnode(phylo)[1]) { diff --git a/R/multi.ace.R b/R/multi.ace.R index 8e31c090..5433b80d 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -179,7 +179,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Check the tree and data cleaned_data <- clean.data(matrix, tree) if(!is.na(cleaned_data$dropped_tips) || !is.na(cleaned_data$dropped_rows)) { - stop(paste0("Some names in the data or the tree(s) are not matching.\nYou can use dispRity::clean.data(", as.expression(match_call$data), ", ", as.expression(match_call$tree), ") to find out more.")) + stop(paste0("Some names in the data or the tree(s) are not matching.\nYou can use dispRity::clean.data(", as.expression(match_call$data), ", ", as.expression(match_call$tree), ") to find out more."), call. = FALSE) } ## Find the node labels (and eventually add them to the trees) @@ -306,7 +306,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec } ## Check dispRity if(output == "dispRity" && do_discrete) { - stop("Only ancestral state estimations for continuous characters can be converted into a dispRity object.\nSelect an other output method.") + stop("Only ancestral state estimations for continuous characters can be converted into a dispRity object.\nSelect an other output method.", call. = FALSE) } } @@ -338,12 +338,12 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Checking for the reserved character reserved <- c("\\@", "@") %in% special.tokens if(any(reserved)) { - stop("special.tokens cannot contain the character '@' since it is reserved for the dispRity::char.diff function.") + stop("special.tokens cannot contain the character '@' since it is reserved for the dispRity::char.diff function.", call. = FALSE) } ## Checking whether the special.tokens are unique if(length(unique(special.tokens)) != length(special.tokens)) { - stop("special.tokens cannot contain duplicated tokens.") + stop("special.tokens cannot contain duplicated tokens.", call. = FALSE) } ## If any special token is NA, convert them as "N.A" temporarily @@ -455,25 +455,25 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec models_continuous <- set.continuous.args.ace.models(models, n = n_characters_continuous) } if(do_discrete && do_continuous) { - stop("Only one model is specified but both discrete and continuous characters are detected.") + stop("Only one model is specified but both discrete and continuous characters are detected.", call. = FALSE) } } else { ## Vector of models if(length(models) != n_characters) { - stop(paste0("Incorrect number of models specified: ", length(models), " models for ", n_characters, " characters.")) + stop(paste0("Incorrect number of models specified: ", length(models), " models for ", n_characters, " characters."), call. = FALSE) } else { check.method(models, c(available_models_discrete, available_models_continuous), msg = "models applied to characters") ## Check models per character types ## Discrete if(sum(models %in% available_models_discrete) != n_characters_discrete) { - stop(paste0("Incorrect number of models specified: ", sum(models %in% available_models_discrete), " models for ", n_characters, " discrete characters.")) + stop(paste0("Incorrect number of models specified: ", sum(models %in% available_models_discrete), " models for ", n_characters, " discrete characters."), call. = FALSE) } else { ## Discrete models (valid) models_discrete <- as.list(models[models %in% available_models_discrete]) } ## Continuous if(sum(models %in% available_models_continuous) != n_characters_continuous) { - stop(paste0("Incorrect number of models specified: ", sum(models %in% available_models_continuous), " models for ", n_characters, " continuous characters.")) + stop(paste0("Incorrect number of models specified: ", sum(models %in% available_models_continuous), " models for ", n_characters, " continuous characters."), call. = FALSE) } else { ## Continuous models (valid) models_continuous <- sapply(models[models %in% available_models_continuous], set.continuous.args.ace.models, n = 1) @@ -485,7 +485,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Models is a transition matrix (discrete only) if(models_class == "matrix") { if(do_continuous) { - stop("Transition matrices can be used as models only for discrete characters.") + stop("Transition matrices can be used as models only for discrete characters.", call. = FALSE) } else { models_discrete <- replicate(n_characters_discrete, models, simplify = FALSE) } @@ -495,7 +495,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec if(models_class == "list") { if(length(models) == 1) { if(do_discrete && do_continuous) { - stop("Only one model is specified but both discrete and continuous characters are detected.") + stop("Only one model is specified but both discrete and continuous characters are detected.", call. = FALSE) } ## Set the models for discrete if(do_discrete) { @@ -532,18 +532,20 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Handle the options ## ######### - - ## options.args - warning("TODO: change castor options to just options (parsed to castor or ape)") - if(missing(options.args)) { ## No options options.args <- NULL } else { ## must be list with names check.class(options.args, "list") + options_error <- "options.args must be an unambiguous named list of options for castor::asr_mk_model() or ape::ace()." + ## Check the available names + options_avail <- c(names(formals(castor::asr_mk_model)), names(formals(ape::ace))) + + + if(is.null(names(options.args))) { - stop("options.args must be a named list of options for castor::asr_mk_model().", call. = FALSE) + stop(options_error, call. = FALSE) } } diff --git a/R/pair.plot_fun.R b/R/pair.plot_fun.R index 3aeaedf1..edd13b55 100755 --- a/R/pair.plot_fun.R +++ b/R/pair.plot_fun.R @@ -21,7 +21,7 @@ find.num.elements <- function(x) { stop(paste("The number of elements to plot is not a function of an entire number.", "To get the right number of combinations, you must satisfy:", " ncol(combn(seq(1:n_elements), 2)) == nrow(data)", - "where 'n_elements' is the number of elements compared pairwise.", sep = "\n")) + "where 'n_elements' is the number of elements compared pairwise.", sep = "\n"), call. = FALSE) } return(n) } \ No newline at end of file diff --git a/R/pgls.dispRity.R b/R/pgls.dispRity.R index 4a3d8211..a9f1accd 100644 --- a/R/pgls.dispRity.R +++ b/R/pgls.dispRity.R @@ -50,7 +50,7 @@ pgls.dispRity <- function(data, tree, formula, model = "BM", ..., optim = list() data <- add.tree(data, tree = tree, replace = TRUE) } else { if(is.null(get.tree(data))) { - stop("No tree was found in the provided data and none was provided through the tree argument.") + stop("No tree was found in the provided data and none was provided through the tree argument.", call. = FALSE) } } @@ -110,14 +110,14 @@ get.formula <- function(disparity) { group <- lapply(disparity$subsets, function(x) return(c(x$elements))) ## Check overlap if(any(table(unlist(group)) != 1)) { - stop("Some groups have overlapping elements.") + stop("Some groups have overlapping elements.", call. = FALSE) } ## Return the correct formula if(disparity$call$subsets[[1]] == "customised") { return(disparity ~ group) } else { ## Warning for time auto-correlation - stop("It is currently not possible to apply an phylogenetic linear model on dispRity data with time series.") + stop("It is currently not possible to apply an phylogenetic linear model on dispRity data with time series.", call. = FALSE) # warning("Data contains time series: the default formula used is disparity ~ time but it does not take time autocorrelation into account.", call. = FALSE) # colnames(group_table) <- "time" # return(list(formula = disparity ~ time, group = NULL, time = group_table)) @@ -165,7 +165,7 @@ get.pgls.data <- function(data) { multiple_trees <- (length(trees) > 1) multiple_datas <- (length(data_list) > 1) if(multiple_datas && multiple_trees) { - stop(paste0("Data must either same number of matrices (", length(data_list), ") and trees (", length(trees) , ") or just one tree or matrix combined with respectively with multiple matrices or trees.")) + stop(paste0("Data must either same number of matrices (", length(data_list), ") and trees (", length(trees) , ") or just one tree or matrix combined with respectively with multiple matrices or trees."), call. = FALSE) } ## Combine the data if(multiple_datas) { diff --git a/R/plot.dispRity_fun.R b/R/plot.dispRity_fun.R index 576c57f0..64f66037 100755 --- a/R/plot.dispRity_fun.R +++ b/R/plot.dispRity_fun.R @@ -1574,7 +1574,7 @@ do.plot.projection <- function(data, specific.args, cent.tend, ...) { ## Get the central tendencies if(!all(specific.args$correlation.plot %in% names(data)) && length(specific.args$correlation.plot) != 2) { - stop(paste0("correlation.plot argument must contain 2 elements from data (data contains: ", paste(names(data), collapse = ", "), ").")) + stop(paste0("correlation.plot argument must contain 2 elements from data (data contains: ", paste(names(data), collapse = ", "), ")."), call. = FALSE) } ## Remove the phylogeny part (if exists) diff --git a/R/reduce.space.R b/R/reduce.space.R index f95c8ec3..90e1ac80 100755 --- a/R/reduce.space.R +++ b/R/reduce.space.R @@ -107,11 +107,11 @@ reduce.space <- function(space, type, remove, parameters, tuning, verbose = FALS if(remove < 100) { remove <- remove/100 } else { - stop("remove must be a probability or a percentage.") + stop("remove must be a probability or a percentage.", call. = FALSE) } } } else { - stop("remove must be a probability or a percentage.") + stop("remove must be a probability or a percentage.", call. = FALSE) } ## Straight returns if remove = 0 or 1 diff --git a/R/remove.zero.brlen.R b/R/remove.zero.brlen.R index 735b2f12..ed79ce01 100755 --- a/R/remove.zero.brlen.R +++ b/R/remove.zero.brlen.R @@ -87,7 +87,7 @@ remove.zero.brlen <- function(tree, slide, verbose = FALSE) { if(any(connect_to_tip <- root_edges[,2] <= Ntip(tree))) { ## Check if that branch length is zero if(tree_bkp$edge.length[bad_edge <- which(root_edges[connect_to_tip,1] == tree_bkp$edge[,1])[connect_to_tip]] == 0) { - stop(paste0("The root of the tree is connecting to a tip with a zero branch length: neither can be slid. You can try moving the tip manually by assigning a value to the following edge:\n ", as.expression(match_call$tree), "$edge.length[",bad_edge ,"] <- your_value")) + stop(paste0("The root of the tree is connecting to a tip with a zero branch length: neither can be slid. You can try moving the tip manually by assigning a value to the following edge:\n ", as.expression(match_call$tree), "$edge.length[",bad_edge ,"] <- your_value"), call. = FALSE) } } diff --git a/R/slice.tree_fun.R b/R/slice.tree_fun.R index bd738828..cd12a2f2 100755 --- a/R/slice.tree_fun.R +++ b/R/slice.tree_fun.R @@ -144,7 +144,7 @@ slice.tree_parent.node <- function(tree, tip) { parent_node <- tree$node.label[parent_edge-Ntip(tree)] #error if not working if (length(parent_node) != 1) { - stop('No parent node found!') + stop("No parent node found!", call. = FALSE) } return(parent_node) } diff --git a/R/slide.nodes.R b/R/slide.nodes.R index b74cd7b3..44d944d8 100755 --- a/R/slide.nodes.R +++ b/R/slide.nodes.R @@ -57,7 +57,7 @@ slide.nodes <- function(nodes, tree, slide, allow.negative.root = FALSE) { ## Getting the node IDs (if character) if(node_class == "character") { if(is.null(tree$node.label)) { - stop("The tree has no node labels, provide the nodes as integers.") + stop("The tree has no node labels, provide the nodes as integers.", call. = FALSE) } nodes <- which(tree$node.label %in% nodes) + Ntip(tree) } @@ -65,14 +65,14 @@ slide.nodes <- function(nodes, tree, slide, allow.negative.root = FALSE) { check.class(tree, "phylo") ## Check whether nodes exist in the tree - if(any(nodes > (Nnode(tree)+Ntip(tree)))) stop("node(s) not found in tree.") - if(any(nodes < Nnode(tree))) stop("node(s) not found in tree.") + if(any(nodes > (Nnode(tree)+Ntip(tree)))) stop("node(s) not found in tree.", call. = FALSE) + if(any(nodes < Nnode(tree))) stop("node(s) not found in tree.", call. = FALSE) if(!allow.negative.root) { if(any(nodes == (Ntip(tree)+1))) warning(paste0("The parent of the root node (", (Ntip(tree) + 1), ") cannot be slid.")) } ## Check whether the tree has edge lengths - if(is.null(tree$edge.length)) stop("The tree has no edge lengths.") + if(is.null(tree$edge.length)) stop("The tree has no edge lengths.", call. = FALSE) ## Slide check.class(slide, c("numeric", "integer")) @@ -82,7 +82,7 @@ slide.nodes <- function(nodes, tree, slide, allow.negative.root = FALSE) { ## Catch eventual errors if(is.null(tree)) { - stop(paste0("The slide value (", slide, ") produced negative branch length(s).")) + stop(paste0("The slide value (", slide, ") produced negative branch length(s)."), call. = FALSE) } return(tree) } \ No newline at end of file diff --git a/R/test.metric.R b/R/test.metric.R index b9e0c7fd..6742fccd 100755 --- a/R/test.metric.R +++ b/R/test.metric.R @@ -144,7 +144,7 @@ test.metric <- function(data, metric, ..., shifts, shift.options, model, replica ## Check the arguments arguments <- names(formals(model)) if(length(arguments) > 1 || arguments != "data") { - stop("model function argument can only take \"data\" as an argument.") + stop("model function argument can only take \"data\" as an argument.", call. = FALSE) } } From 0638cd666482e8f3f974cee01f3ca50a9d34b9db Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Tue, 23 Jan 2024 13:36:51 +0000 Subject: [PATCH 10/15] need to add parallel back --- R/multi.ace.R | 40 ++++++++++++++++++++++----------- tests/testthat/test-multi.ace.R | 19 +++++----------- 2 files changed, 33 insertions(+), 26 deletions(-) diff --git a/R/multi.ace.R b/R/multi.ace.R index 5433b80d..90973595 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -390,8 +390,14 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec has_invariants <- TRUE ## Stop if they are only invariant characters - if(length(invariants) == n_characters_discrete) { - warning(match_call$data, " contains only invariant discrete characters.") + if(do_continuous) { + if(length(invariants) == n_characters_discrete) { + warning(match_call$data, " contains only invariant discrete characters.") + } + } else { + if(length(invariants) == n_characters_discrete) { + stop.call(call = match_call$data, " contains only invariant characters.") + } } ## Remove the characters @@ -534,19 +540,30 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ######### if(missing(options.args)) { ## No options - options.args <- NULL + options.ace <- options.castor <- options.args <- NULL } else { ## must be list with names check.class(options.args, "list") options_error <- "options.args must be an unambiguous named list of options for castor::asr_mk_model() or ape::ace()." ## Check the available names options_avail <- c(names(formals(castor::asr_mk_model)), names(formals(ape::ace))) - - - - if(is.null(names(options.args))) { + if(is.null(names(options.args)) || !all(names(options.args) %in% options_avail)) { stop(options_error, call. = FALSE) } + ## Sort the options + options.ace <- options.castor <- NULL + if(do_continuous) { + options.ace <- options.args[names(options.args) %in% names(formals(ape::ace))] + if(length(options.ace) == 0) { + options.ace <- NULL + } + } + if(do_discrete) { + options.castor <- options.args[names(options.args) %in% names(formals(castor::asr_mk_model))] + if(length(options.castor) == 0) { + options.castor <- NULL + } + } } ## Check the estimation details @@ -577,7 +594,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ## Setting the continuous characters call if(do_continuous) { ## Create the character arguments - character_continuous_args <- mapply(function(character, ace.args) return(c(x = character, ace.args)), characters_continuous, models_continuous, SIMPLIFY = FALSE) + character_continuous_args <- mapply(function(character, ace.args, options = NULL) return(c(x = character, ace.args, options)), characters_continuous, models_continuous, MoreArgs = list(options = options.ace), SIMPLIFY = FALSE) ## Create the character and tree arguments tree_character_continuous_args <- list() @@ -596,7 +613,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec } ## Set up the arguments for one tree - character_discrete_args <- mapply(make.args, characters_tables, characters_states, models_discrete, MoreArgs = list(estimation.details = details_out), SIMPLIFY = FALSE) + character_discrete_args <- mapply(make.args, characters_tables, characters_states, models_discrete, MoreArgs = list(estimation.details = details_out, castor.options = options.castor), SIMPLIFY = FALSE) ## Create the character and tree arguments tree_character_discrete_args <- list() @@ -661,7 +678,7 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec # } } else { ## Make the functions verbose - if(verbose) cat("Running ancestral states estimations:\n") + if(verbose) cat("Running ancestral states estimations:") ## Run the continuous characters if(do_continuous) { @@ -777,9 +794,6 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec dispRity = make.dispRity(data = lapply(results_out, add.tips, matrix = matrix), tree = tree) ) - ## Handle details - warning("TODO: handle details") - ## Results out if(is.null(estimation.details)) { if(length(tree) == 1) { diff --git a/tests/testthat/test-multi.ace.R b/tests/testthat/test-multi.ace.R index 508789ca..84fcf2d9 100755 --- a/tests/testthat/test-multi.ace.R +++ b/tests/testthat/test-multi.ace.R @@ -23,7 +23,6 @@ test_that("multi.ace works", { parallel = FALSE, output = "list") - error <- capture_error(multi.ace(data = "matrix_complex", tree = tree_test, models = "ER", @@ -144,7 +143,6 @@ test_that("multi.ace works", { output = "something")) expect_equal(error[[1]], "output option must be one of the following: matrix, list, combined, combined.list, combined.matrix, dispRity.") - error <- capture_error(multi.ace(data = matrix_complex, tree = tree_test, models = "ER", @@ -223,7 +221,7 @@ test_that("multi.ace works", { verbose = FALSE, parallel = FALSE, output = "list")) - expect_equal(error[[1]], "models should be list of characters or/and matrices of length 10.") + expect_equal(error[[1]], "models list must be the same length as the number of characters (10).") ## Castor options works well error <- capture_error(results <- multi.ace(data = list_matrix, @@ -231,8 +229,8 @@ test_that("multi.ace works", { verbose = FALSE, parallel = FALSE, output = "list", - castor.options = list(2))) - expect_equal(error[[1]], "castor.options must be a named list of options for castor::asr_mk_model().") + options.args = list(2))) + expect_equal(error[[1]], "options.args must be an unambiguous named list of options for castor::asr_mk_model() or ape::ace().") ## Threshold works well results <- multi.ace(data = list_matrix, @@ -315,7 +313,6 @@ test_that("multi.ace works", { expect_is(ancestral_states[[1]], "matrix") expect_equal(dim(ancestral_states[[1]]), c(11, 10)) - ## Parallel works # expect_is(multi.ace(matrix_test, tree_test, parallel = TRUE), "list") # test_verbose <- capture.output(test <- multi.ace(matrix_test, tree_test, parallel = 2, verbose = TRUE)) @@ -385,10 +382,7 @@ test_that("multi.ace works", { estimation.details = c("loglikelihood", "transition_matrix"))) expect_equal(test, c("Preparing the data:.....Done." , - "Running ancestral states estimations:" , - ".......... Done.", - "Running ancestral states estimations:" , - ".......... Done.")) + "Running ancestral states estimations:....................Done.")) # set.seed(3) # test <- capture.output(results <- multi.ace(data = matrix_complex, @@ -422,7 +416,7 @@ test_that("multi.ace works with continuous and mix", { data <- space.maker(elements = 15, dimensions = 5, distribution = rnorm, elements.name = tree$tip.label) ## Run the multi.ace on the continuous data - test <- multi.ace(data = data, tree = tree, output = "combined.matrix", verbose = TRUE) + expect_warning(test <- multi.ace(data = data, tree = tree, output = "combined.matrix", verbose = TRUE)) ## Works well for continuous expect_is(test, "matrix") @@ -436,11 +430,10 @@ test_that("multi.ace works with continuous and mix", { data <- cbind(data, "new_char2" = as.character(sample(1:2, 15, replace = TRUE))) ## Works well for mixed characters - test <- multi.ace(data = data, tree = tree, output = "combined.matrix") + expect_warning(test <- multi.ace(data = data, tree = tree, output = "combined.matrix")) expect_is(test, "data.frame") expect_equal(dim(test), c(15+14, 7)) expect_equal(sort(rownames(test)), sort(c(tree$tip.label, tree$node.label))) - classes <- character() for(i in 1:ncol(test)) { classes[i] <- class(test[, i]) From 426d465dfc3860a06752eee7e4d8f7f1a9b80a4f Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Tue, 23 Jan 2024 16:11:14 +0000 Subject: [PATCH 11/15] fancier plots --- NEWS.md | 19 ++++++------------ R/plot.dispRity_fun.R | 31 +++++++++++++++++++++++++---- tests/testthat/test-plot.dispRity.R | 15 ++++++++++++++ 3 files changed, 48 insertions(+), 17 deletions(-) diff --git a/NEWS.md b/NEWS.md index 336e7f05..a93ca8b2 100755 --- a/NEWS.md +++ b/NEWS.md @@ -1,23 +1,16 @@ -dispRity v1.8.3 (2024-01-18) +dispRity v1.8.3 (2024-01-23) ========================= ### NEW FEATURES - * Redesigned `multi.ace` to be more modular and handle both continuous and/or discrete characters. - - inc. changed default `models` argument from `"ER"` to `"ER"` or `"BM"` depending on the character types. - - inc. change `"castor.options"` to the more generic `"options.args"` - - [ ] `models` arg now takes character names + complex list (for `ape::ace`) - - [x] `models` selection is now done after sorting the characters (for default: discrete = ER, continuous = BM/REML) - - [ ] add explanations of new `models` argument in the doc - - [ ] change `castor.options` to `options` in general (with a named list: `list(castor = ..., ape = ...)`) - - [ ] new internal function for running `ape::ace` - - [ ] combine results more smartly - - [ ] `estimation.details` can now take arguments from both `castor` and `ape` - - [ ] new output `dispRity` format + * Redesigned `multi.ace` to be more modular and handle both continuous and/or discrete characters. Changes include a **change in argument name** from `castor.options` to the generic `options.args` (the options can be provided the same way as before though); and a **change in default arguments** for `models` which can now be left missing (previously was `"ER"`) and applies `"ER"` and `"BM"` for respectively discrete and continuous characters by default. + - [ ] rework `parallel` options + - [ ] update manual. - [ ] also update manual ### MINOR IMPROVEMENTS - * `custom.subsets` can now take a logical vector for the `group` argument + * `custom.subsets` can now take a logical vector for the `group` argument. + * `plot` functions doing scatter plot now centers them without changing the scale of both axes. ### BUG FIXES diff --git a/R/plot.dispRity_fun.R b/R/plot.dispRity_fun.R index 64f66037..5e61d862 100755 --- a/R/plot.dispRity_fun.R +++ b/R/plot.dispRity_fun.R @@ -750,9 +750,13 @@ do.plot.preview <- function(data, specific.args, ...) { plot_args <- get.dots(dots, plot_args, "ylab", paste0("Dimension ", specific.args$dimensions[2], " (", loading[specific.args$dimensions[2]], "%)")) ## Setting plot limits - plot_lim <- range(unlist(lapply(data$matrix[specific.args$matrix], function(matrix, dim) c(matrix[, dim]), dim = specific.args$dimensions))) - plot_args <- get.dots(dots, plot_args, "xlim", plot_lim) - plot_args <- get.dots(dots, plot_args, "ylim", plot_lim) + xrange <- range(unlist(lapply(data$matrix[specific.args$matrix], function(matrix, dim) c(matrix[, dim]), dim = specific.args$dimensions[1]))) + yrange <- range(unlist(lapply(data$matrix[specific.args$matrix], function(matrix, dim) c(matrix[, dim]), dim = specific.args$dimensions[2]))) + ## Get the centered scale range + plot_lims <- get.center.scale.range(xrange, yrange) + + plot_args <- get.dots(dots, plot_args, "xlim", plot_lims$xlim) + plot_args <- get.dots(dots, plot_args, "ylim", plot_lims$ylim) ## Get the number of colour groups n_groups <- length(data$subsets) @@ -1596,4 +1600,23 @@ do.plot.projection <- function(data, specific.args, cent.tend, ...) { ## Plot the results plot(plot_data, ...) } -} \ No newline at end of file +} + + +## Get centered scaled range for prettier plots +get.center.scale.range <- function(xrange, yrange) { + ## Get the ranges + x_diff <- diff(xrange) + y_diff <- diff(yrange) + + ## Largest range stays unchanged + if(x_diff >= y_diff) { + xlim <- xrange + ylim <- c(mean(yrange)-(x_diff/2), mean(yrange)+(x_diff/2)) + } + if(x_diff < y_diff) { + ylim <- yrange + xlim <- c(mean(xrange)-(y_diff/2), mean(xrange)+(y_diff/2)) + } + return(list(xlim = xlim, ylim = ylim)) +} diff --git a/tests/testthat/test-plot.dispRity.R b/tests/testthat/test-plot.dispRity.R index 977d0dbf..5cb75850 100755 --- a/tests/testthat/test-plot.dispRity.R +++ b/tests/testthat/test-plot.dispRity.R @@ -415,4 +415,19 @@ test_that("preview works with fuzzy matrices and trees", { expect_null(plot(data)) expect_null(plot(data, specific.args = list(tree = TRUE))) expect_null(plot(data, specific.args = list(matrix = 1, tree = 1))) +}) + +test_that("get.center.scale.range gives the correct scales" { + set.seed(1) + ## X bigger + xrange <- range(rnorm(10)) + yrange <- range(runif(10)) + test <- get.center.scale.range(xrange, yrange) + expect_gt(diff(xrange), diff(yrange)) + expect_equal(diff(test$xlim), diff(test$ylim)) + + yrange <- range(runif(10, max = 100)) + test <- get.center.scale.range(xrange, yrange) + expect_lt(diff(xrange), diff(yrange)) + expect_equal(diff(test$xlim), diff(test$ylim)) }) \ No newline at end of file From 30616cc47d6aee10f98bd95d683e67e413ffa2f5 Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Tue, 23 Jan 2024 16:11:43 +0000 Subject: [PATCH 12/15] multi ace working on example --- R/multi.ace.R | 2 ++ R/multi.ace_fun.R | 3 ++- man/multi.ace.Rd | 2 +- tests/testthat/test-multi.ace.R | 2 +- 4 files changed, 6 insertions(+), 3 deletions(-) diff --git a/R/multi.ace.R b/R/multi.ace.R index 90973595..156d8752 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -583,6 +583,8 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec } ## Check the requested details check.method(estimation.details, return_args, msg = "estimation.details") + } else { + return_args_discrete <- return_args_continuous <- NULL } ######### diff --git a/R/multi.ace_fun.R b/R/multi.ace_fun.R index 3228da14..b1ff5959 100755 --- a/R/multi.ace_fun.R +++ b/R/multi.ace_fun.R @@ -311,9 +311,10 @@ bind.details <- function(continuous, discrete, order) { if(!is.null(discrete[[1]])) { names(discrete_details) <- names(discrete[[1]]) } + discrete <- discrete_details } ## Bind the two lists - return(c(continuous, discrete_details)[c(order$continuous, order$discrete)]) + return(c(continuous, discrete)[c(order$continuous, order$discrete)]) } diff --git a/man/multi.ace.Rd b/man/multi.ace.Rd index c58e56a0..ca106c32 100755 --- a/man/multi.ace.Rd +++ b/man/multi.ace.Rd @@ -42,7 +42,7 @@ multi.ace( \item{options.args}{optional, a named list of options to be passed to function called by \code{castor::asr_mk_model}.} -\item{estimation.details}{optional, whether to also return the details for each estimation as returned by \code{castor::asr_mk_model}. This argument can be left \code{NULL} (default) or be any combination of the elements returned by \code{castor::asr_mk_model} (e.g. \code{c("loglikelihood", "transition_matrix")}).} +\item{estimation.details}{optional, whether to also return the details for each estimation as returned by \code{castor::asr_mk_model} or \code{ape::ace}. This argument can be left \code{NULL} (default) or be any combination of the elements returned by \code{castor::asr_mk_model} or \code{ape::ace} (e.g. \code{c("loglikelihood", "transition_matrix", "CI95")}).} } \value{ Returns a \code{"matrix"} or \code{"list"} of ancestral states. By default, the function returns the ancestral states in the same format as the input \code{matrix}. This can be changed using the option \code{output = "matrix"} or \code{"list"} to force the class of the output. diff --git a/tests/testthat/test-multi.ace.R b/tests/testthat/test-multi.ace.R index 84fcf2d9..fb52b40c 100755 --- a/tests/testthat/test-multi.ace.R +++ b/tests/testthat/test-multi.ace.R @@ -416,7 +416,7 @@ test_that("multi.ace works with continuous and mix", { data <- space.maker(elements = 15, dimensions = 5, distribution = rnorm, elements.name = tree$tip.label) ## Run the multi.ace on the continuous data - expect_warning(test <- multi.ace(data = data, tree = tree, output = "combined.matrix", verbose = TRUE)) + expect_warning(test <- multi.ace(data = data, tree = tree, output = "combined.matrix", verbose = FALSE)) ## Works well for continuous expect_is(test, "matrix") From 25f343494dd273e412c518d28b655022dbfa03a6 Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Wed, 24 Jan 2024 11:18:33 +0000 Subject: [PATCH 13/15] Update test-plot.dispRity.R --- tests/testthat/test-plot.dispRity.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/testthat/test-plot.dispRity.R b/tests/testthat/test-plot.dispRity.R index 5cb75850..567e15b5 100755 --- a/tests/testthat/test-plot.dispRity.R +++ b/tests/testthat/test-plot.dispRity.R @@ -417,7 +417,7 @@ test_that("preview works with fuzzy matrices and trees", { expect_null(plot(data, specific.args = list(matrix = 1, tree = 1))) }) -test_that("get.center.scale.range gives the correct scales" { +test_that("get.center.scale.range gives the correct scales", { set.seed(1) ## X bigger xrange <- range(rnorm(10)) From 20521b4fdfccf916bb9f0a53f875465c0c004017 Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Wed, 24 Jan 2024 11:18:51 +0000 Subject: [PATCH 14/15] parallel working and tested --- R/multi.ace.R | 169 ++++++++++++++++++-------------- tests/testthat/test-multi.ace.R | 19 +++- 2 files changed, 111 insertions(+), 77 deletions(-) diff --git a/R/multi.ace.R b/R/multi.ace.R index 156d8752..59f021b1 100755 --- a/R/multi.ace.R +++ b/R/multi.ace.R @@ -4,19 +4,14 @@ #' #' @param data A \code{matrix}, \code{data.frame} or \code{list} with the characters for each taxa. #' @param tree A \code{phylo} or \code{mutiPhylo} object (if the \code{tree} argument contains node labels, they will be used to name the output). -#' @param models A \code{vector} of models to be passed to \code{castor::asr_mk_model}. - - -#TODO: add info about continuous characters - - +#' @param models A \code{character} vector, unambiguous named \code{list} or \code{matrix} to be passed as model arguments to \code{castor::asr_mk_model} or \code{ape::ace} (see details). #' @param threshold either \code{logical} for applying a relative threshold (\code{TRUE} - default) or no threshold (\code{FALSE}) or a \code{numeric} value of the threshold (e.g. 0.95). See details. #' @param special.tokens optional, a named \code{vector} of special tokens to be passed to \code{\link[base]{grep}} (make sure to protect the character with \code{"\\\\"}). By default \code{special.tokens <- c(missing = "\\\\?", inapplicable = "\\\\-", polymorphism = "\\\\&", uncertainty = "\\\\/")}. Note that \code{NA} values are not compared and that the symbol "@" is reserved and cannot be used. #' @param special.behaviours optional, a \code{list} of one or more functions for a special behaviour for \code{special.tokens}. See details. #' @param brlen.multiplier optional, a vector of branch length modifiers (e.g. to convert time branch length in changes branch length) or a list of vectors (the same length as \code{tree}). #' @param verbose \code{logical}, whether to be verbose (\code{TRUE}) or not (\code{FALSE} - default). -#' @param parallel \code{logical}, whether to use parallel algorithm (\code{TRUE}) or not (\code{FALSE} - default). -#' @param output optional, see Value section below. +#' @param parallel Either a \code{logical}, whether to use parallel algorithm (\code{TRUE}) or not (\code{FALSE} - default); or directly an \code{integer} indicating the number of cores to use (note that if \code{parallel = 1}, one core will be used but the parallel integration will still be called). +#' @param output optional, see \code{Value} section below. #' @param options.args optional, a named list of options to be passed to function called by \code{castor::asr_mk_model}. #' @param estimation.details optional, whether to also return the details for each estimation as returned by \code{castor::asr_mk_model} or \code{ape::ace}. This argument can be left \code{NULL} (default) or be any combination of the elements returned by \code{castor::asr_mk_model} or \code{ape::ace} (e.g. \code{c("loglikelihood", "transition_matrix", "CI95")}). #' @@ -66,14 +61,8 @@ #' @return #' Returns a \code{"matrix"} or \code{"list"} of ancestral states. By default, the function returns the ancestral states in the same format as the input \code{matrix}. This can be changed using the option \code{output = "matrix"} or \code{"list"} to force the class of the output. #' To output the combined ancestral states and input, you can use \code{"combined"} (using the input format) or \code{"combined.matrix"} or \code{"combined.list"}. +#' If using continuous characters only, you can use the output option \code{"dispRity"} to directly output a usable \code{dispRity} object with all trees and all the data (estimated and input). #' \emph{NOTE} that if the input data had multiple character types (continuous and discrete) and that \code{"matrix"} or \code{"combined.matrix"} output is requested, the function returns a \code{"data.frame"}. - - -#TODO: add dispRity format - - - -# To output the light version to be passed to \code{dispRity} functions (a list of two elements: 1) the input \code{matrix} and 2) a list of ancestral states matrices) you can use \code{output = "dispRity"}. #' #' @examples #' set.seed(42) @@ -309,6 +298,10 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec stop("Only ancestral state estimations for continuous characters can be converted into a dispRity object.\nSelect an other output method.", call. = FALSE) } } + ## Set data.frame output to matrix + if(output == "data.frame") { + output <- "matrix" + } ## Handle the tokens # special.tokens <- character(); special.behaviours <- list() ; warning("DEBUG: multi.ace") @@ -603,6 +596,16 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec for(one_tree in 1:length(tree)) { tree_character_continuous_args[[one_tree]] <- lapply(character_continuous_args, function(character, tree) {character$phy <- tree; return(character)}, tree[[one_tree]]) } + ## Set verbose fun + if(verbose) { + fun_continuous <- function(...) { + cat(".") + return(ape::ace(...)) + } + } else { + fun_continuous <- ape::ace + } + } ## Setting the discrete characters call if(do_discrete) { @@ -631,70 +634,90 @@ multi.ace <- function(data, tree, models, threshold = TRUE, special.tokens, spec ######### if(do_parallel) { - stop("TODO: do parallel multi.ace") - ## Do parallel! - # if(do_parallel) { - # ## Remove verbose - # if(verbose) { - # cat(paste0("Running the estimation for ", length(tree), " tree", ifelse(length(tree) > 1, "s ", " "), "using ", cores, " core", ifelse(cores == 1, "...", "s..."))) - # was_verbose <- TRUE - # verbose <- FALSE - # } else { - # was_verbose <- FALSE - # } - - # ## Set up the cluster - # cluster <- parallel::makeCluster(cores) - - # ## Get the current environment - # current_env <- environment() - # ## Get the export lists - # export_arguments_list <- c("tree_args_list", - # "special.tokens", - # "invariants", - # "threshold.type", - # "threshold", - # "verbose", - # "characters_states", - # "invariant_characters_states") - # export_functions_list <- c("one.tree.ace", - # "castor.ace", - # "tree.data.update", - # "add.state.names", - # "translate.likelihood") - - # ## Export from this environment - # parallel::clusterExport(cluster, c(export_arguments_list, export_functions_list), envir = current_env) - - # ## Call the cluster - # results_out <- parLapply(cl = cluster, tree_args_list, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) - - # ## Stop the cluster - # parallel::stopCluster(cluster) - - # ## Reactivate the verbose - # if(was_verbose) { - # cat("Done.") - # verbose <- TRUE - # } - # } - } else { - ## Make the functions verbose - if(verbose) cat("Running ancestral states estimations:") + ## Remove verbose + if(verbose) { + cat(paste0("Running the estimation for ", length(tree), " tree", ifelse(length(tree) > 1, "s ", " "), "using ", cores, " core", ifelse(cores == 1, "...", "s..."))) + was_verbose <- TRUE + verbose <- FALSE + } else { + was_verbose <- FALSE + } - ## Run the continuous characters + ## Set up the cluster + cluster <- parallel::makeCluster(cores) + + ## Get the current environment + current_env <- environment() + + export_arguments_list <- export_functions_list <- character() + + if(do_discrete) { + ## Get the export lists + export_arguments_list <- c("tree_character_discrete_args", + "special.tokens", + "invariants", + "threshold.type", + "threshold", + "verbose", + "characters_states", + "invariant_characters_states") + export_functions_list <- c("one.tree.ace", + "castor.ace", + "tree.data.update", + "add.state.names", + "translate.likelihood") + } if(do_continuous) { + export_arguments_list <- c(export_arguments_list, "tree_character_continuous_args") + export_functions_list <- c(export_functions_list, "fun_continuous") + } + + ## Export from this environment + parallel::clusterExport(cluster, c(export_arguments_list, export_functions_list), envir = current_env) - ## Set verbose fun - if(verbose) { - fun_continuous <- function(...) { - cat(".") - return(ape::ace(...)) + ## Call the cluster + if(do_discrete) { + discrete_estimates <- parLapply(cl = cluster, tree_character_discrete_args, one.tree.ace, special.tokens, invariants, characters_states, threshold.type, threshold, invariant_characters_states, verbose) + } + if(do_continuous) { + continuous_estimates <- parLapply(cl = cluster, tree_character_continuous_args, lapply, function(x) do.call(fun_continuous, x)) + ## Remove the ugly call + continuous_estimates <- lapply(continuous_estimates, lapply, function(x) {x$call <- "ape::ace"; return(x)}) + ## Add node labels + # TODO: to be removed if ape update + add.nodes <- function(obj, phy) { + ## Adding node labels to $ace and $CI95 (if available) + options(warn = -1) + names <- as.integer(names(obj$ace)) + options(warn = 0) + if(!is.null(phy$node.label) && !is.na(names[1])) { + ordered_node_labels <- phy$node.label[c(as.integer(names(obj$ace))- Ntip(phy))] + names(obj$ace) <- ordered_node_labels + if(!is.null(obj$CI95)) { + rownames(obj$CI95) <- ordered_node_labels + } } - } else { - fun_continuous <- ape::ace + return(obj) } + for(one_tree in 1:length(tree)) { + continuous_estimates[[one_tree]] <- lapply(continuous_estimates[[one_tree]], add.nodes, phy = tree[[one_tree]]) + } + } + + ## Stop the cluster + parallel::stopCluster(cluster) + ## Reactivate the verbose + if(was_verbose) { + cat("Done.") + verbose <- TRUE + } + } else { + ## Make the functions verbose + if(verbose) cat("Running ancestral states estimations:") + + ## Run the continuous characters + if(do_continuous) { ## Run all the ace continuous_estimates <- lapply(tree_character_continuous_args, lapply, function(x) do.call(fun_continuous, x)) ## Remove the ugly call diff --git a/tests/testthat/test-multi.ace.R b/tests/testthat/test-multi.ace.R index fb52b40c..60ae3210 100755 --- a/tests/testthat/test-multi.ace.R +++ b/tests/testthat/test-multi.ace.R @@ -314,10 +314,10 @@ test_that("multi.ace works", { expect_equal(dim(ancestral_states[[1]]), c(11, 10)) ## Parallel works - # expect_is(multi.ace(matrix_test, tree_test, parallel = TRUE), "list") - # test_verbose <- capture.output(test <- multi.ace(matrix_test, tree_test, parallel = 2, verbose = TRUE)) - # expect_is(test, "list") - # expect_equal(test_verbose, c("Preparing the data:.....Done.", "Running the estimation for 2 trees using 2 cores...Done.")) + expect_is(multi.ace(matrix_test, tree_test, parallel = 1), "list") + test_verbose <- capture.output(test <- multi.ace(matrix_test, tree_test, parallel = 2, verbose = TRUE)) + expect_is(test, "list") + expect_equal(test_verbose, c("Preparing the data:.....Done.", "Running the estimation for 2 trees using 2 cores...Done.")) ## Examples work set.seed(42) @@ -439,4 +439,15 @@ test_that("multi.ace works with continuous and mix", { classes[i] <- class(test[, i]) } expect_equal(unique(classes), c("numeric", "character")) + + ## Works for parallel + test <- multi.ace(data = data, tree = tree, parallel = 1) + expect_is(test, "data.frame") + expect_equal(dim(test), c(14, 7)) + expect_equal(sort(rownames(test)), sort(c(tree$node.label))) + classes <- character() + for(i in 1:ncol(test)) { + classes[i] <- class(test[, i]) + } + expect_equal(unique(classes), c("numeric", "character")) }) \ No newline at end of file From e4a1d432dd0cd65ab54c602dce2c94af871f9853 Mon Sep 17 00:00:00 2001 From: Thomas Guillerme Date: Wed, 24 Jan 2024 11:20:36 +0000 Subject: [PATCH 15/15] updated multi.ace for 1.8.3 --- NEWS.md | 5 +- inst/gitbook/05_other-functionalities.Rmd | 67 ++++++++++++++++++++++- man/multi.ace.Rd | 7 ++- 3 files changed, 69 insertions(+), 10 deletions(-) diff --git a/NEWS.md b/NEWS.md index a93ca8b2..9c04d74a 100755 --- a/NEWS.md +++ b/NEWS.md @@ -1,11 +1,8 @@ -dispRity v1.8.3 (2024-01-23) +dispRity v1.8.3 (2024-01-24) ========================= ### NEW FEATURES * Redesigned `multi.ace` to be more modular and handle both continuous and/or discrete characters. Changes include a **change in argument name** from `castor.options` to the generic `options.args` (the options can be provided the same way as before though); and a **change in default arguments** for `models` which can now be left missing (previously was `"ER"`) and applies `"ER"` and `"BM"` for respectively discrete and continuous characters by default. - - [ ] rework `parallel` options - - [ ] update manual. - - [ ] also update manual ### MINOR IMPROVEMENTS diff --git a/inst/gitbook/05_other-functionalities.Rmd b/inst/gitbook/05_other-functionalities.Rmd index 01d9d148..3a14d163 100755 --- a/inst/gitbook/05_other-functionalities.Rmd +++ b/inst/gitbook/05_other-functionalities.Rmd @@ -516,7 +516,7 @@ tree.age(tree, order = "present") ## `multi.ace` -This function allows to run the `ape::ace` function (ancestral characters estimations) on multiple trees. +This function allows to run ancestral characters estimations on multiple trees. In it's most basic structure (e.g. using all default arguments) this function is using a mix of `ape::ace` and `castor::asr_mk_model` depending on the data and the situation and is generally faster than both functions when applied to a list of trees. However, this function provides also some more complex and modular functionalities, especially appropriate when using discrete morphological character data. @@ -636,7 +636,7 @@ For example, you can use different models for each character via the `models` ar ### Feeding the results to `char.diff` to get distance matrices -Finally, after running your ancestral states estimations, it is not uncommon to then use these resulting data to calculate the distances between taxa and then ordinate the results to measure disparity. +After running your ancestral states estimations, it is not uncommon to then use these resulting data to calculate the distances between taxa and then ordinate the results to measure disparity. You can do that using the `char.diff` function [described above](#char.diff) but instead of measuring the distances between characters (columns) you can measure the distances between species (rows). You might notice that this function uses the same modular token and behaviour descriptions. That makes sense because they're using the same core C functions implemented in dispRity that greatly speed up distance calculations. @@ -658,4 +658,65 @@ We can then feed these matrices directly to `char.diff`, say for calculating the distances <- lapply(ancestral_states, char.diff, method = "mord", by.col = FALSE) ``` -And we now have a list of distances matrices with ancestral states estimated! \ No newline at end of file +And we now have a list of distances matrices with ancestral states estimated! + +### Running ancestral states estimations for continuous characters + +You can also run `multi.ace` on continuous characters. +The function detects any continuous characters as being of class `"numeric"` and runs them using the `ape::ace` function. + +```{r} +set.seed(1) +## Creating three coalescent trees +my_trees <- replicate(3, rcoal(15), simplify = FALSE) +## Adding node labels +my_trees <- lapply(my_trees, makeNodeLabel) +## Making into a multiPhylo object +class(my_trees) <- "multiPhylo" + +## Creating a matrix of continuous characters +data <- space.maker(elements = 15, dimensions = 5, distribution = rnorm, + elements.name = my_trees[[1]]$tip.label) +``` + +With such data and trees you can easily run the `multi.ace` estimations. +By default, the estimations use the default arguments from `ape::ace`, knowingly a Brownian Motion (`model = "BM"`) with the REML method (`method = "REML"`; this method "first estimates the ancestral value at the root (aka, the phylogenetic mean), then the variance of the Brownian motion process is estimated by optimizing the residual log-likelihood" - from `?ape::ace`). + +```{r} +## Running multi.ace on continuous data +my_ancestral_states <- multi.ace(data, my_trees) +## We end up with three matrices of node states estimates +str(my_ancestral_states) +``` + +This results in three matrices with ancestral states for the nodes. +When using continuous characters, however, you can output the results directly as a `dispRity` object that allows visualisation and other normal dispRity pipeline: + +```{r} +## Running multi.ace on continuous data +my_ancestral_states <- multi.ace(data, my_trees, output = "dispRity") +## We end up with three matrices of node states estimates +plot(my_ancestral_states) +``` + +You can also mix continuous and discrete characters together. +By default the `multi.ace` detects which character is of which type and applies the correct estimations based on that. +However you can always specify models or other details character per characters. + + +```{r} +## Adding two discrete characters +data <- as.data.frame(data) +data <- cbind(data, "new_char" = as.character(sample(1:2, 15, replace = TRUE))) +data <- cbind(data, "new_char2" = as.character(sample(1:2, 15, replace = TRUE))) + +## Setting up different models for each characters +## BM for all 5 continuous characters +## and ER and ARD for the two discrete ones +my_models <- c(rep("BM", 5), "ER", "ARD") + +## Running the estimation with the specified models +my_ancestral_states <- multi.ace(data, my_trees, models = my_models) +``` + +Of course all the options discussed in the first part above also can apply here! \ No newline at end of file diff --git a/man/multi.ace.Rd b/man/multi.ace.Rd index ca106c32..38c85db0 100755 --- a/man/multi.ace.Rd +++ b/man/multi.ace.Rd @@ -24,7 +24,7 @@ multi.ace( \item{tree}{A \code{phylo} or \code{mutiPhylo} object (if the \code{tree} argument contains node labels, they will be used to name the output).} -\item{models}{A \code{vector} of models to be passed to \code{castor::asr_mk_model}.} +\item{models}{A \code{character} vector, unambiguous named \code{list} or \code{matrix} to be passed as model arguments to \code{castor::asr_mk_model} or \code{ape::ace} (see details).} \item{threshold}{either \code{logical} for applying a relative threshold (\code{TRUE} - default) or no threshold (\code{FALSE}) or a \code{numeric} value of the threshold (e.g. 0.95). See details.} @@ -36,9 +36,9 @@ multi.ace( \item{verbose}{\code{logical}, whether to be verbose (\code{TRUE}) or not (\code{FALSE} - default).} -\item{parallel}{\code{logical}, whether to use parallel algorithm (\code{TRUE}) or not (\code{FALSE} - default).} +\item{parallel}{Either a \code{logical}, whether to use parallel algorithm (\code{TRUE}) or not (\code{FALSE} - default); or directly an \code{integer} indicating the number of cores to use (note that if \code{parallel = 1}, one core will be used but the parallel integration will still be called).} -\item{output}{optional, see Value section below.} +\item{output}{optional, see \code{Value} section below.} \item{options.args}{optional, a named list of options to be passed to function called by \code{castor::asr_mk_model}.} @@ -47,6 +47,7 @@ multi.ace( \value{ Returns a \code{"matrix"} or \code{"list"} of ancestral states. By default, the function returns the ancestral states in the same format as the input \code{matrix}. This can be changed using the option \code{output = "matrix"} or \code{"list"} to force the class of the output. To output the combined ancestral states and input, you can use \code{"combined"} (using the input format) or \code{"combined.matrix"} or \code{"combined.list"}. +If using continuous characters only, you can use the output option \code{"dispRity"} to directly output a usable \code{dispRity} object with all trees and all the data (estimated and input). \emph{NOTE} that if the input data had multiple character types (continuous and discrete) and that \code{"matrix"} or \code{"combined.matrix"} output is requested, the function returns a \code{"data.frame"}. } \description{