Skip to contents

Maximum likehood estimation for (SecSSE)

Source: R/secsse_ml.R
cla_secsse_ml.Rd

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 = "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
)

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) 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. 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 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::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).

Value

Parameters estimated and maximum likelihood

Examples

# 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))
#> You are optimizing lambda0 mu0 
#> You are fixing lambda1 mu1 
#> Optimizing the likelihood - this may take a while. 
#> The loglikelihood for the initial parameter values is -67.5669039297515.
#> 
#> Maximum likelihood parameter estimates: lambda0: 15.604290, mu0: 15.874464, lambda1: 0.000000, mu1: 0.000000: 
#> Maximum loglikelihood: 25.282018
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-03, 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 <- 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)
#> Note: you set some transitions as impossible to happen.
#> 1 15.6042898684248 15.6042898684248 15.6042898684248 -121.487688555362 initial 
#> 2 16.3845043618461 16.3845043618461 14.2366516329269 -115.87210515215 expand 
#> 3 16.3845043618461 16.3845043618461 14.2366516329269 -115.87210515215 reflect 
#> 4 15.8564618714564 16.1164115482592 13.4435417361893 -112.708850473379 expand 
#> 5 17.8616456533575 16.661129036425 11.7804033067207 -106.427345504873 expand 
#> 6 16.661129036425 18.8748937365321 10.7024010667367 -102.649991462868 expand 
#> 7 17.5464711985784 18.8748937365321 8.83574291518734 -96.8407395467205 expand 
#> 8 21.2672928057985 23.7499826152845 6.8803975853617 -92.1784042042333 expand 
#> 9 19.2372760644554 34.9835268796627 5.34799248335495 -90.1181535926296 expand 
#> 10 22.6933092609064 33.8535718517169 4.81579185777643 -89.8916337567571 reflect 
#> 11 22.6933092609064 33.8535718517169 4.81579185777643 -89.8916337567571 reflect 
#> 12 22.6933092609064 33.8535718517169 4.81579185777643 -89.8916337567571 contract outside 
#> 13 22.6933092609064 33.8535718517169 4.81579185777643 -89.8916337567571 reflect 
#> 14 22.6933092609064 33.8535718517169 4.81579185777643 -89.8916337567571 contract inside 
#> 15 22.6933092609064 33.8535718517169 4.81579185777643 -89.8916337567571 reflect 
#> 16 22.6933092609064 33.8535718517169 4.81579185777643 -89.8916337567571 contract inside 
#> 17 22.6933092609064 33.8535718517169 4.81579185777643 -89.8916337567571 contract inside 
#> 18 22.6933092609064 33.8535718517169 4.81579185777643 -89.8916337567571 contract inside 
#> 19 21.5954238570094 40.0611420950045 4.64166201307183 -89.8907285434026 contract inside 
#> 20 21.5954238570094 40.0611420950045 4.64166201307183 -89.8907285434026 contract inside 
#> 21 22.09252160483 38.7838963866625 4.68738843123826 -89.8892657403911 contract inside 
#> 22 22.09252160483 38.7838963866625 4.68738843123826 -89.8892657403911 reflect 
#> 23 22.4146047038665 36.0332588409853 4.73313430169762 -89.8867774721457 contract inside 
#> 24 22.4146047038665 36.0332588409853 4.73313430169762 -89.8867774721457 reflect 
#> 25 22.4146047038665 36.0332588409853 4.73313430169762 -89.8867774721457 contract inside 
#> 26 22.4146047038665 36.0332588409853 4.73313430169762 -89.8867774721457 reflect 
#> 27 23.1619267394853 33.7939708537398 4.67617743982144 -89.8846394388549 expand 
#> 28 22.9967027764633 32.9570145766782 4.76019916146665 -89.8845176846956 reflect 
#> 29 25.9077673007118 28.8858477767871 4.66368227158788 -89.8827694880962 expand 
#> 30 27.7298236564981 25.5641274625481 4.63398478686598 -89.8786893066793 expand 
#> 31 31.395418958883 22.2163034259215 4.70407868110388 -89.8743829342342 expand 
#> 32 31.395418958883 22.2163034259215 4.70407868110388 -89.8743829342342 contract inside 
#> 33 33.1020565748236 20.0301019807209 4.72372769998256 -89.8702086639009 expand 
#> 34 33.1020565748236 20.0301019807209 4.72372769998256 -89.8702086639009 reflect 
#> 35 33.1020565748236 20.0301019807209 4.72372769998256 -89.8702086639009 contract outside 
#> 36 33.1020565748236 20.0301019807209 4.72372769998256 -89.8702086639009 contract outside 
#> 37 33.1020565748236 20.0301019807209 4.72372769998256 -89.8702086639009 reflect 
#> 38 33.5964443113748 19.2007230053455 4.73425605791087 -89.8687094544291 reflect 
#> 39 29.7262450758438 21.852133724945 4.71093041222478 -89.8681251857471 reflect 
#> 40 29.7262450758438 21.852133724945 4.71093041222478 -89.8681251857471 contract inside 
#> 41 30.575466286766 20.4365061188632 4.69143428701633 -89.8641760452911 expand 
#> 42 30.575466286766 20.4365061188632 4.69143428701633 -89.8641760452911 reflect 
#> 43 30.575466286766 20.4365061188632 4.69143428701633 -89.8641760452911 reflect 
#> 44 26.902167800435 21.4498889292796 4.68802298364535 -89.8561108566437 expand 
#> 45 26.902167800435 21.4498889292796 4.68802298364535 -89.8561108566437 reflect 
#> 46 26.902167800435 21.4498889292796 4.68802298364535 -89.8561108566437 reflect 
#> 47 23.2912492321529 22.1288807686561 4.6175356044978 -89.8453679007587 expand 
#> 48 25.1791704795815 18.8060270739465 4.67267401030954 -89.8381783694394 expand 
#> 49 19.1938996609292 23.7069393574904 4.65026318942184 -89.8333475178622 expand 
#> 50 16.4805975240266 21.1666443055712 4.56592684998421 -89.8164897620991 expand 
#> 51 14.9769117628143 19.1592610465612 4.65280890366182 -89.8023543431036 expand 
#> 52 14.9769117628143 19.1592610465612 4.65280890366182 -89.8023543431036 reflect 
#> 53 14.9769117628143 19.1592610465612 4.65280890366182 -89.8023543431036 reflect 
#> 54 16.063236451579 16.7816791253283 4.60155492072445 -89.7960968434354 reflect 
#> 55 16.063236451579 16.7816791253283 4.60155492072445 -89.7960968434354 reflect 
#> 56 15.9381640876111 16.4024847909521 4.67219593044267 -89.7956270297493 contract outside 
#> 57 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 reflect 
#> 58 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 contract inside 
#> 59 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 reflect 
#> 60 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 contract inside 
#> 61 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 contract inside 
#> 62 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 reflect 
#> 63 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 reflect 
#> 64 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 reflect 
#> 65 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 reflect 
#> 66 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 contract inside 
#> 67 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 contract inside 
#> 68 13.8854168337677 14.3888583685557 4.61844421632923 -89.7902278102615 contract inside 
#> 69 14.3019139360004 14.479947109031 4.632966660664 -89.7901344817794 contract outside 
#> 70 14.3019139360004 14.479947109031 4.632966660664 -89.7901344817794 reflect 
#> 71 14.4426137546406 14.3165211694009 4.62620711474184 -89.7900934687737 contract inside 
#> 72 14.4426137546406 14.3165211694009 4.62620711474184 -89.7900934687737 reflect 
#> 73 14.4426137546406 14.3165211694009 4.62620711474184 -89.7900934687737 contract inside 
#> 74 14.4426137546406 14.3165211694009 4.62620711474184 -89.7900934687737 contract inside 
#> 75 14.4426137546406 14.3165211694009 4.62620711474184 -89.7900934687737 reflect 
#> 76 14.2521569922074 14.2700219864961 4.62904905591931 -89.7900838525261 contract inside 
#> Optimization has terminated successfully. 
# [1] -90.9763