secsse_ml: Maximum likehood estimation for (SecSSE)

In secsse: Several Examined and Concealed States-Dependent Speciation and Extinction

Maximum likehood estimation for (SecSSE)

Description

Maximum likehood estimation under Several examined and concealed trait States dependent Speciation and Extinction (SecSSE)

Usage

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 = "simplex",
  num_cycles = 1,
  loglik_penalty = 0,
  is_complete_tree = FALSE,
  take_into_account_root_edge = FALSE,
  verbose = FALSE,
  num_threads = 1,
  atol = 1e-08,
  rtol = 1e-07,
  method = "odeint::runge_kutta_cash_karp54",
  use_normalization = TRUE,
  return_root_state = FALSE,
  see_ancestral_states = FALSE
)

Arguments

phy	phylogenetic tree of class phylo, rooted and with branch lengths. Alternatively, multiple phylogenetic trees can be provided as the multiPhylo class.
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"). When providing a multiPhylo set of multiple phylognies, traits should be a list where each entry in the list corresponds to the matching phylogeny on that position.
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), ⁠"equal_weights"'. or ⁠"stationary_weights"⁠It can also be specified for the root state: the vector⁠c(1, 0, 0)⁠ indicates state 1 was the root state. When using a⁠multiPhylo' object, root_state_weight should be list where each entry in the list corresponds to the root_state_weight for each tree.
sampling_fraction	vector that states the sampling proportion per trait state. It must have as many elements as there are trait states. When using a multiPhylo object, sampling fraction should be list where each entry in the list corresponds to the sampling proportion for each tree.
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 "simplex". Alternative is "subplex". 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.
take_into_account_root_edge	if TRUE, the LL integration is continued along the root edge. This also affects conditioning (as now, conditioning no longer needs to assume a speciation event at the start of the tree)
verbose	sets verbose output; default is TRUE.
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	ODE integration method. Choose from: "odeint::runge_kutta_cash_karp54", "odeint::runge_kutta_fehlberg78", "odeint::runge_kutta_dopri5", "odeint::bulirsch_stoer" and "odeint::runge_kutta4". Default method is: "odeint::runge_kutta_cash_karp54".
use_normalization	normalize the density vector during integration, more accurate but slower (default = TRUE)
return_root_state	if TRUE, returns the state of the system at the root, this can be useful to use as the starting point of a simulation. When used in ML, after finishing the ML optimization, the found optimum is evaluated one more time to retrieve the root state (to avoid having to store the root state every ML evaluation).
see_ancestral_states	Boolean for whether the ancestral states for each of the internal nodes should be output. Defaults to FALSE.

Value

A list with the following elements $MLpars: the maximum likelihood parameter estimates $ML: the maximum likelihood of the data (phylogeny + tip states) given the parameters (speciation, extinction, transition rates). $conv: whether the optimization converged or not If see_ancestral_states = TRUE, then there will be two additional elements: $ancestral_states: a matrix with the probabilities of each state at the internal nodes $states: a matrix with the probabilities E, D (normalized) and S that are used in the calculations. The ancestral_states matrix is a submatrix of this matrix. This matrix is mostly used for package developers. If return_root_state = TRUE, then there will be one additional element: $root_state: vector with probabilities of each state at the root. This vector is the same as the top row of $ancestral_states We have used the shorthand description of "probabilities of each state", but technically, the probabilities are the normalized probabilities D of the data given each state at the internal nodes.

Examples

# Example of how to set the arguments for a ML search.
library(secsse)
library(DDD)
set.seed(13)
# lambdas for 0A and 1A and 2A are the same but need to be estimated
# mus are fixed to
# the transition rates are constrained to be equal and fixed 0.01
phylotree <- ape::rcoal(31, tip.label = 1:31)
traits <-  sample(c(0,1,2), ape::Ntip(phylotree),replace=TRUE)#get some traits
num_concealed_states<-3
idparslist <- id_paramPos(traits, num_concealed_states)
idparslist[[1]][c(1,4,7)] <- 1
idparslist[[1]][c(2,5,8)] <- 2
idparslist[[1]][c(3,6,9)] <- 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 <- DDD::bd_ML(brts = ape::branching.times(phylotree))
intGuessLamba <- startingpoint$lambda0
intGuessMu <- startingpoint$mu0
idparsopt <- c(1,2,3,5)
initparsopt <- c(rep(intGuessLamba,3),rep((intGuessLamba/5),1))
idparsfix <- c(0,4)
parsfix <- c(0,0)
tol <- c(1e-02, 1e-03, 1e-03)
maxiter <- 1000 * round((1.25)^length(idparsopt))
optimmethod <- 'simplex'
cond <- 'proper_cond'
root_state_weight <- 'proper_weights'
sampling_fraction <- c(1,1,1)
model<-secsse_ml(
phylotree,
traits,
num_concealed_states,
idparslist,
idparsopt,
initparsopt,
idparsfix,
parsfix,
cond,
root_state_weight,
sampling_fraction,
tol,
maxiter,
optimmethod,
num_cycles = 1,
verbose = FALSE)
# model$ML
# [1] -16.47099

secsse documentation built on May 15, 2026, 5:06 p.m.