multi.ace: Ancestral states estimations with multiple trees

Description

Fast ancestral states estimations run on multiple trees using the Mk model from castor::asr_mk_model.

Usage

multi.ace(
  data,
  tree,
  models,
  ml.collapse = list(type = "majority"),
  special.tokens,
  special.behaviours,
  brlen.multiplier,
  verbose = FALSE,
  parallel = FALSE,
  output,
  options.args,
  estimation.details = NULL
)

Value

Returns a "matrix" or "list" of ancestral states. By default, the function returns the ancestral states in the same format as the input matrix. This can be changed using the option output = "matrix" or "list" to force the class of the output. To output the combined ancestral states and input, you can use "combined" (using the input format) or "combined.matrix" or "combined.list". If using continuous characters only, you can use the output option "dispRity" to directly output a usable dispRity object with all trees and all the data (estimated and input). NOTE that if the input data had multiple character types (continuous and discrete) and that "matrix" or "combined.matrix" output is requested, the function returns a "data.frame".

If using the option output = "multi.ace", the function returns a "multi.ace" object that can be recycled in the multi.ace function. This can be useful for testing different sampling strategies of threshold methods without having to recalculate all the ancestral trait estimations.

Arguments

data: A matrix, data.frame or list with the characters for each taxa. Or a multi.ace list (see details).
tree: A phylo or multiPhylo object (if the tree argument contains node labels, they will be used to name the output).
models: A character vector, unambiguous named list or matrix to be passed as model arguments to castor::asr_mk_model or ape::ace (see details).
ml.collapse: A named list for specifying how to collapse maximum likelihood ancestral state estimates. Must contain $type which can be "majority" (default; returns the discrete state(s) with highest likelihood, using uncertainty token if tied, and a point estimate (mean) for continuous traits), "relative" (returns discrete states with scaled likelihood >= max likelihood - 1/n_states), "threshold" (returns discrete states >= a specified threshold), or "sample" (samples discrete states by their scaled likelihoods and continuous traits using specified sampling functions). See details for additional elements.
special.tokens: optional, a named vector of special tokens to be passed to grep (make sure to protect the character with "\\"). By default special.tokens <- c(missing = "\\?", inapplicable = "\\-", polymorphism = "\\&", uncertainty = "\\/"). Note that NA values are not compared and that the symbol "@" is reserved and cannot be used.
special.behaviours: optional, a list of one or more functions for a special behaviour for special.tokens. See details.
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 tree).
verbose: logical, whether to be verbose (TRUE) or not (FALSE - default).
parallel: Either a logical, whether to use parallel algorithm (TRUE) or not (FALSE - default); or directly an integer indicating the number of cores to use (note that if parallel = 1, one core will be used but the parallel integration will still be called).
output: optional, see Value section below.
options.args: optional, a named list of options to be passed to function called by castor::asr_mk_model.
estimation.details: optional, whether to also return the details for each estimation as returned by castor::asr_mk_model or ape::ace. This argument can be left NULL (default) or be any combination of the elements returned by castor::asr_mk_model or ape::ace (e.g. c("loglikelihood", "transition_matrix", "CI95")).

Author

Thomas Guillerme, Caleb Scutt

Details

Depending on the type of characters models argument can be either:

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 models = "ER" applies the Equal Rates model to all characters (assuming they are all discrete characters).
a vector of model names to apply to different type of characters (see below for the list). For example models = c("ER", "ER", "BM") applies the Equal Rates model to the two first characters (discrete) and the "BM" model to the third character (continuous).
a transition "matrix" to be applied to all characters (if discrete). For example models = matrix(0.2, 2, 2).
an single named list of arguments to be applied to all characters by passing it to ape::ace (if continuous). For example models = list(method = "GLS", corStruct = corBrownian(1, my_tree)).
an un-ambiguous list of arguments to be passed to either castor::asr_mk_model (discrete characters) or ape::ace (continuous characters). For example 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 castor::asr_mk_model are:

"ER" for all equal rates
"SYM" for symmetric rates
"ARD" all rates are different
"SUEDE" equal stepwise transitions (e.g. for meristic/counting characters)
"SRD" different stepwise transitions

See directly castor::asr_mk_model for more models.

The available built-in models and methods for continuous characters in ape::ace are:

"BM" model: for a default Brownian Motion with the "REML" method
"REML" method: for a default Brownian Motion with the "REML" method (same as above)
"ML" method: for a default Brownian Motion with the "ML" method
"pic" method: for a default Brownian Motion with the "pic" (least squared) method

special.behaviours allows to generate a special rule for the special.tokens. The functions should can take the arguments character, all_states with character being the character that contains the special token and 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.

missing = function(x,y) y
inapplicable = function(x,y) y
polymorphism = function(x,y) strsplit(x, split = "\\&")[[1]]
uncertainty = function(x,y) strsplit(x, split = "\\/")[[1]]

Functions in the list must be named following the special token of concern (e.g. missing), have only x, y as inputs and a single output a single value (that gets coerced to integer automatically). For example, the special behaviour for the special token "?" can be coded as: special.behaviours = list(missing = function(x, y) return(NA) to make ignore the character for taxa containing "?".

When using the parallel option (either through using parallel = TRUE by using the number of available cores minus on or manually setting the number of cores - e.g. parallel = 5), the castor::asr_mk_model function will use the designated number of cores (using the option Nthreads = <requested_number_of_cores>). Additionally, if the input tree is a "multiPhylo" object, the trees will be run in parallel for each number of cores, thus decreasing computation time accordingly (e.g. if 3 cores are requested and tree contains 12 "phylo" objects, 4 different "phylo" objects will be run in parallel on the 3 cores making the calculation around 3 times faster).

Additional listed elements for ml.collapse = list():

list(type = "majority", tie.breaker = TRUE) randomly selects one state when highest likelihoods are tied;
list(type = "threshold" , threshold = 0.95) can take any numeric value (0-1);
list(type = "sample", sample = n, sample.fun = list(fun, param)) where n is an integer > 1 that denotes the number of matrices to sample continuous and/or discrete states across. sample.fun specifies a function and parameters for the sampling of continuous traits. The default is sample.fun = list(fun = runif, param = list(min = min, max = max)) for applying a random uniform sampling (runif) with the parameters (the minimum and the maximum are applied using respectively the min and max functions on the estimated data). For applying different samplings to different traits, you can use a list of arguments in the sample format as sample.fun (e.g. sample.fun = list(trait_uniform = list(fun = runif, param = list(min = min, max = max)), trait_normal = list(fun = rnorm, param = list(mean = mean, sd = function(x)return(diff(range(x))/4))) - here the standard deviation is calculated as a quarter of the 95% CI range). Note that if type = "sample", sample = 1, the mean point estimate will be returned if continuous and majority rule applied if discrete.

It is also possible to just run the estimations without sampling but for sampling later (i.e. just run the estimations; save them as a multi.ace object and then rerun them into the function). You can do that by using the option output = "multi.ace". Using this option ignores the following options: ml.collapse, and estimation.details.

Examples

Run this code

set.seed(34)
## A simple example:
## A random tree with 10 tips
tree <- rcoal(10)
## Setting up the parameters
my_rates = c(rgamma, rate = 10, shape = 5)

## A random Mk matrix (10*50)
matrix_simple <- sim.morpho(tree, characters = 50, model = "ER", rates = my_rates,
                            invariant = FALSE)

## Run a basic ancestral states estimations
ancestral_states <- multi.ace(matrix_simple, tree)
ancestral_states[1:5, 1:5]

## A random continuous trait matrix
matrix_continuous <- space.maker(elements = 10, dimensions = 2,
                                 distribution = rnorm, 
                                 elements.names = tree$tip.label)
## A basic estimations for discrete data
ancestral_traits <- multi.ace(matrix_continuous, tree)

## A more complex example for discrete data
## Create a multiple list of 5 trees
multiple_trees <- rmtree(5, 10)

## Modify the matrix to contain missing and special data
matrix_complex <- matrix_simple
matrix_complex[sample(1:length(matrix_complex), 50)] <- "-"
matrix_complex[sample(1:length(matrix_complex), 50)] <- "0%2"
matrix_complex[sample(1:length(matrix_complex), 50)] <- "?"
matrix_complex[1:5,1:5]

## Set a list of extra special tokens
my_spec_tokens <- c("weirdtoken" = "%")

## Set some special behaviours for the "weirdtoken" and for "-" and "?"
my_spec_behaviours <- list()
## Inapplicable tokens "-" are ignored
my_spec_behaviours$inapplicable <- function(x,y) return(NA)
## Missing tokens "?" are considered as all states
my_spec_behaviours$missing      <- function(x,y) return(y)
## Weird tokens are considered as state 0 and 3
my_spec_behaviours$weirdtoken   <- function(x,y) return(c(1,2))

## Create a random branch length modifier to apply to each tree
branch_lengths <- rnorm(18)^2

## Setting a list of model ("ER" for the 25 first characters and then "SYM")
my_models <- c(rep("ER", 25), rep("SYM", 25))

## Run the ancestral states on all the tree with multiple options
ancestral_states <- multi.ace(matrix_complex, multiple_trees,
                              verbose = TRUE,
                              models = my_models,
                              ml.collapse = list(type = "majority", tie.breaker = TRUE),
                              special.tokens = my_spec_tokens,
                              special.behaviours = my_spec_behaviours,
                              brlen.multiplier = branch_lengths,
                              output = "combined.matrix")

## The results for the the two first characters for the first tree
ancestral_states[[1]][, 1:2]

if (FALSE) {
## The same example but running in parallel
ancestral_states <- multi.ace(matrix_complex, multiple_trees,
                              verbose = TRUE,
                              models = my_models,
                              ml.collapse = list(type = "majority", tie.breaker = TRUE),
                              special.tokens = my_spec_tokens,
                              special.behaviours = my_spec_behaviours,
                              brlen.multiplier = branch_lengths,
                              output = "combined.matrix",
                              parallel = TRUE)
}

## Saving your estimations as a "multi.ace" object
ancestor_estimates <- multi.ace(matrix_simple, multiple_trees, output = "multi.ace")

## Estimating ancestral states by sampling distributions
## Discrete characters:
sampled_ancestors <- multi.ace(ancestor_estimates,
     ml.collapse = list(type = "sample",
                        sample = 25))

## Continuous characters:
## Sampling 100 matrices from the ancestral estimations
sampled_traits <- multi.ace(matrix_continuous, tree,
     ml.collapse = list(type = "sample",
                        sample = 100))
## Setting up two specific sampling functions
sample_fun <- list(
     unif_sampl = list(fun = runif,
                       param = list(min = min, max = max)),
     norm_sampl = list(fun = rnorm,
                       param = list(mean = mean,
                                    sd = function(x)return(diff(range(x))/4)))
)
## Sampling 100 matrices with different distribution functions
sampled_traits <- multi.ace(matrix_continuous, tree,
     ml.collapse = list(type = "sample",
                        sample = 100,
                        sample.fun = sample_fun))

Run the code above in your browser using DataLab