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::bulirsch_stoer",
  use_normalization = TRUE
)

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::bulirsch_stoer".

use_normalization

normalize the density vector during integration, more accurate but slower (default = TRUE)

Value

Parameter 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.487683644503 initial 
#> 2 16.3845043618461 16.3845043618461 14.2366516329269 -115.872102944723 expand 
#> 3 16.3845043618461 16.3845043618461 14.2366516329269 -115.872102944723 reflect 
#> 4 16.1164115482592 15.8564618714564 13.4435417361893 -112.708848926855 expand 
#> 5 16.661129036425 17.8616456533575 11.7804033067207 -106.427344244374 expand 
#> 6 18.8748937365321 16.661129036425 10.7024010667367 -102.649990642589 expand 
#> 7 18.8748937365321 17.5464711985784 8.83574291518734 -96.8407390014834 expand 
#> 8 23.7499826152845 21.2672928057985 6.8803975853617 -92.1784039041511 expand 
#> 9 34.9835268796627 19.2372760644554 5.34799248335495 -90.118153445057 expand 
#> 10 33.8535718517169 22.6933092609064 4.81579185777643 -89.8916336320649 reflect 
#> 11 33.8535718517169 22.6933092609064 4.81579185777643 -89.8916336320649 reflect 
#> 12 33.8535718517169 22.6933092609064 4.81579185777643 -89.8916336320649 contract outside 
#> 13 33.8535718517169 22.6933092609064 4.81579185777643 -89.8916336320649 reflect 
#> 14 33.8535718517169 22.6933092609064 4.81579185777643 -89.8916336320649 contract inside 
#> 15 33.8535718517169 22.6933092609064 4.81579185777643 -89.8916336320649 reflect 
#> 16 33.8535718517169 22.6933092609064 4.81579185777643 -89.8916336320649 contract inside 
#> 17 33.8535718517169 22.6933092609064 4.81579185777643 -89.8916336320649 contract inside 
#> 18 33.8535718517169 22.6933092609064 4.81579185777643 -89.8916336320649 contract inside 
#> 19 40.0611420950045 21.5954238570094 4.64166201307183 -89.8907284337715 contract inside 
#> 20 40.0611420950045 21.5954238570094 4.64166201307183 -89.8907284337715 contract inside 
#> 21 38.7838963866625 22.09252160483 4.68738843123826 -89.8892656096984 contract inside 
#> 22 38.7838963866625 22.09252160483 4.68738843123826 -89.8892656096984 reflect 
#> 23 36.0332588409853 22.4146047038665 4.73313430169762 -89.8867773517614 contract inside 
#> 24 36.0332588409853 22.4146047038665 4.73313430169762 -89.8867773517614 reflect 
#> 25 36.0332588409853 22.4146047038665 4.73313430169762 -89.8867773517614 contract inside 
#> 26 36.0332588409853 22.4146047038665 4.73313430169762 -89.8867773517614 reflect 
#> 27 33.7939708537398 23.1619267394853 4.67617743982144 -89.8846392966006 expand 
#> 28 32.9570145766782 22.9967027764633 4.76019916146665 -89.8845175271118 reflect 
#> 29 28.8858477767871 25.9077673007118 4.66368227158788 -89.882769361436 expand 
#> 30 25.5641274625481 27.7298236564981 4.63398478686598 -89.8786891959933 expand 
#> 31 22.2163034259215 31.395418958883 4.70407868110388 -89.8743828184079 expand 
#> 32 22.2163034259215 31.395418958883 4.70407868110388 -89.8743828184079 contract inside 
#> 33 20.0301019807209 33.1020565748236 4.72372769998256 -89.8702085435553 expand 
#> 34 20.0301019807209 33.1020565748236 4.72372769998256 -89.8702085435553 reflect 
#> 35 20.0301019807209 33.1020565748236 4.72372769998256 -89.8702085435553 contract outside 
#> 36 20.0301019807209 33.1020565748236 4.72372769998256 -89.8702085435553 contract outside 
#> 37 20.0301019807209 33.1020565748236 4.72372769998256 -89.8702085435553 reflect 
#> 38 19.2007230053455 33.5964443113748 4.73425605791087 -89.8687093432768 reflect 
#> 39 21.852133724945 29.7262450758438 4.71093041222478 -89.8681250737671 reflect 
#> 40 21.852133724945 29.7262450758438 4.71093041222478 -89.8681250737671 contract inside 
#> 41 20.4365061188632 30.575466286766 4.69143428701633 -89.8641759306454 expand 
#> 42 20.4365061188632 30.575466286766 4.69143428701633 -89.8641759306454 reflect 
#> 43 20.4365061188632 30.575466286766 4.69143428701633 -89.8641759306454 reflect 
#> 44 21.4498889292796 26.902167800435 4.68802298364535 -89.8561107507925 expand 
#> 45 21.4498889292796 26.902167800435 4.68802298364535 -89.8561107507925 reflect 
#> 46 21.4498889292796 26.902167800435 4.68802298364535 -89.8561107507925 reflect 
#> 47 22.1288807686561 23.2912492321529 4.6175356044978 -89.8453677814197 expand 
#> 48 18.8060270739465 25.1791704795815 4.67267401030954 -89.8381782639351 expand 
#> 49 23.7069393574904 19.1938996609292 4.65026318942184 -89.8333474179573 expand 
#> 50 21.1666443055712 16.4805975240266 4.56592684998421 -89.8164896536086 expand 
#> 51 19.1592610465612 14.9769117628143 4.65280890366182 -89.8023542183295 expand 
#> 52 19.1592610465612 14.9769117628143 4.65280890366182 -89.8023542183295 reflect 
#> 53 19.1592610465612 14.9769117628143 4.65280890366182 -89.8023542183295 reflect 
#> 54 16.7816791253283 16.063236451579 4.60155492072445 -89.7960967620457 reflect 
#> 55 16.7816791253283 16.063236451579 4.60155492072445 -89.7960967620457 reflect 
#> 56 16.4024847909521 15.9381640876111 4.67219593044267 -89.7956269434704 contract outside 
#> 57 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 reflect 
#> 58 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 contract inside 
#> 59 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 reflect 
#> 60 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 contract inside 
#> 61 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 contract inside 
#> 62 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 reflect 
#> 63 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 reflect 
#> 64 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 reflect 
#> 65 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 reflect 
#> 66 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 contract inside 
#> 67 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 contract inside 
#> 68 14.3888583685557 13.8854168337677 4.61844421632923 -89.7902277035283 contract inside 
#> 69 14.479947109031 14.3019139360004 4.632966660664 -89.790134364218 contract outside 
#> 70 14.479947109031 14.3019139360004 4.632966660664 -89.790134364218 reflect 
#> 71 14.3165211694009 14.4426137546406 4.62620711474184 -89.7900933522345 contract inside 
#> 72 14.3165211694009 14.4426137546406 4.62620711474184 -89.7900933522345 reflect 
#> 73 14.3165211694009 14.4426137546406 4.62620711474184 -89.7900933522345 contract inside 
#> 74 14.3165211694009 14.4426137546406 4.62620711474184 -89.7900933522345 contract inside 
#> 75 14.3165211694009 14.4426137546406 4.62620711474184 -89.7900933522345 reflect 
#> 76 14.2700219864961 14.2521569922074 4.62904905591931 -89.790083734614 contract inside 
#> Optimization has terminated successfully. 
# [1] -90.9763