cla_secsse_ml: Maximum likehood estimation for (SecSSE)

Description

Maximum likehood estimation under Several examined and concealed States-dependent Speciation and Extinction (SecSSE) with cladogenetic option

Usage

cla_secsse_ml(
  phy,
  traits,
  num_concealed_states,
  idparslist,
  idparsopt,
  initparsopt,
  idparsfix,
  parsfix,
  cond = "proper_cond",
  root_state_weight = "proper_weights",
  sampling_fraction,
  tol = c(1e-04, 1e-05, 1e-07),
  maxiter = 1000 * round((1.25)^length(idparsopt)),
  optimmethod = "subplex",
  num_cycles = 1,
  loglik_penalty = 0,
  is_complete_tree = FALSE,
  verbose = (optimmethod == "simplex"),
  num_threads = 1,
  atol = 1e-08,
  rtol = 1e-07,
  method = "odeint::bulirsch_stoer"
)

Value

Parameter estimated and maximum likelihood

Arguments

phy: phylogenetic tree of class phylo, rooted and with branch lengths.
traits: vector with trait states for each tip in the phylogeny. The order of the states must be the same as the tree tips. For help, see vignette("starting_secsse", package = "secsse").
num_concealed_states: number of concealed states, generally equivalent to the number of examined states in the dataset.
idparslist: overview of parameters and their values.
idparsopt: a numeric vector with the ID of parameters to be estimated.
initparsopt: a numeric vector with the initial guess of the parameters to be estimated.
idparsfix: a numeric vector with the ID of the fixed parameters.
parsfix: a numeric vector with the value of the fixed parameters.
cond: condition on the existence of a node root: "maddison_cond", "proper_cond" (default). For details, see vignette.
root_state_weight: the method to weigh the states: "maddison_weights", "proper_weights" (default) or "equal_weights". It can also be specified for the root state: the vector c(1, 0, 0) indicates state 1 was the root state.
sampling_fraction: vector that states the sampling proportion per trait state. It must have as many elements as there are trait states.
tol: A numeric vector with the maximum tolerance of the optimization algorithm. Default is c(1e-04, 1e-05, 1e-05).
maxiter: max number of iterations. Default is 1000 * round((1.25) ^ length(idparsopt)).
optimmethod: A string with method used for optimization. Default is "subplex". Alternative is "simplex" and it shouldn't be used in normal conditions (only for debugging). Both are called from DDD::optimizer(), simplex is implemented natively in DDD, while subplex is ultimately called from subplex::subplex().
num_cycles: Number of cycles of the optimization. When set to Inf, the optimization will be repeated until the result is, within the tolerance, equal to the starting values, with a maximum of 10 cycles.
loglik_penalty: the size of the penalty for all parameters; default is 0 (no penalty).
is_complete_tree: logical specifying whether or not a tree with all its extinct species is provided. If set to TRUE, it also assumes that all all extinct lineages are present on the tree. Defaults to FALSE.
verbose: sets verbose output; default is TRUE when optimmethod is "simplex". If optimmethod is set to "simplex", then even if set to FALSE, optimizer output will be shown.
num_threads: number of threads to be used. Default is one thread.
atol: A numeric specifying the absolute tolerance of integration.
rtol: A numeric specifying the relative tolerance of integration.
method: integration method used, available are: "odeint::runge_kutta_cash_karp54", "odeint::runge_kutta_fehlberg78", "odeint::runge_kutta_dopri5", "odeint::bulirsch_stoer" and "odeint::runge_kutta4". Default method is: "odeint::bulirsch_stoer".

Examples

Run this code

# Example of how to set the arguments for a ML search.
library(secsse)
library(DDD)
set.seed(13)
# Check the vignette for a better working exercise.
# lambdas for 0A and 1A and 2A are the same but need to be estimated
# (CTD model, see Syst Biol paper)
# mus are fixed to zero,
# the transition rates are constrained to be equal and fixed 0.01
phylotree <- ape::rcoal(31, tip.label = 1:31)
#get some traits
traits <-  sample(c(0,1,2), ape::Ntip(phylotree), replace = TRUE)
num_concealed_states <- 3
idparslist <- cla_id_paramPos(traits,num_concealed_states)
idparslist$lambdas[1,] <- c(1,1,1,2,2,2,3,3,3)
idparslist[[2]][] <- 4
masterBlock <- matrix(5,ncol = 3,nrow = 3,byrow = TRUE)
diag(masterBlock) <- NA
diff.conceal <- FALSE
idparslist[[3]] <- q_doubletrans(traits,masterBlock,diff.conceal)
startingpoint <- bd_ML(brts = ape::branching.times(phylotree))
intGuessLamba <- startingpoint$lambda0
intGuessMu <- startingpoint$mu0
idparsopt <- c(1,2,3)
initparsopt <- c(rep(intGuessLamba,3))
idparsfix <- c(0,4,5)
parsfix <- c(0,0,0.01)
tol <- c(1e-04, 1e-05, 1e-07)
maxiter <- 1000 * round((1.25) ^ length(idparsopt))
optimmethod <- 'subplex'
cond <- 'proper_cond'
root_state_weight <- 'proper_weights'
sampling_fraction <- c(1,1,1)
model <- cla_secsse_ml(
 phylotree,
 traits,
 num_concealed_states,
 idparslist,
 idparsopt,
 initparsopt,
 idparsfix,
 parsfix,
 cond,
 root_state_weight,
 sampling_fraction,
 tol,
 maxiter,
 optimmethod,
 num_cycles = 1,
 num_threads = 1,
 verbose = FALSE)
# [1] -90.97626

Run the code above in your browser using DataLab