Skip to contents

Maximum likehood estimation for (SecSSE)

Source: R/secsse_ml.R
secsse_ml.Rd

Maximum likehood estimation under Several examined and concealed 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
)

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

Parameter estimated and maximum likelihood

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))
#> 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,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)
#> Note: you set some transitions as impossible to happen.
#> 1 15.6042898684248 15.6042898684248 15.6042898684248 3.12085797368496 -55.0829367206248 initial 
#> 2 16.1826565594201 16.1826565594201 14.2366516329269 3.23679263471176 -53.5955799417801 expand 
#> 3 16.1826565594201 16.1826565594201 14.2366516329269 3.23679263471176 -53.5955799417801 reflect 
#> 4 15.6744474204993 15.0418248706591 14.0813764954283 3.43976602145215 -51.0733848195086 expand 
#> 5 15.6744474204993 15.0418248706591 14.0813764954283 3.43976602145215 -51.0733848195086 reflect 
#> 6 15.7897266624575 16.3431211350014 12.1251758810308 3.51957331290552 -48.3309421257716 expand 
#> 7 15.7897266624575 16.3431211350014 12.1251758810308 3.51957331290552 -48.3309421257716 reflect 
#> 8 15.4239724497177 14.7480876746975 11.1658804598539 4.18949505076817 -43.7404845144724 expand 
#> 9 15.4239724497177 14.7480876746975 11.1658804598539 4.18949505076817 -43.7404845144724 reflect 
#> 10 16.5037295388755 15.9290391482382 8.60698610560038 4.62942837051011 -40.0935241249407 expand 
#> 11 16.5037295388755 15.9290391482382 8.60698610560038 4.62942837051011 -40.0935241249407 reflect 
#> 12 15.8415238453205 13.3794084219105 7.17017313753921 6.95418165007756 -34.1209268706555 expand 
#> 13 15.8415238453205 13.3794084219105 7.17017313753921 6.95418165007756 -34.1209268706555 reflect 
#> 14 15.8415238453205 13.3794084219105 7.17017313753921 6.95418165007756 -34.1209268706555 reflect 
#> 15 17.5043421500954 13.3823278970245 4.87199464662075 18.4611852376756 -32.8526378038631 expand 
#> 16 15.9467686168646 12.5755678121138 4.78859464019957 27.3897761206681 -30.4983662829612 reflect 
#> 17 15.9467686168646 12.5755678121138 4.78859464019957 27.3897761206681 -30.4983662829612 reflect 
#> 18 15.9467686168646 12.5755678121138 4.78859464019957 27.3897761206681 -30.4983662829612 reflect 
#> 19 15.9467686168646 12.5755678121138 4.78859464019957 27.3897761206681 -30.4983662829612 contract inside 
#> 20 15.147576617686 10.5613404064817 4.69711763222928 94.2610341671819 -28.3896435001547 expand 
#> 21 15.147576617686 10.5613404064817 4.69711763222928 94.2610341671819 -28.3896435001547 contract inside 
#> 22 15.147576617686 10.5613404064817 4.69711763222928 94.2610341671819 -28.3896435001547 reflect 
#> 23 13.2944179344031 11.785197674832 6.05959674478364 12.5257381509452 -27.9360254115791 expand 
#> 24 12.427159941311 10.025417328563 5.11292486227123 201.453133415635 -25.6515083796731 expand 
#> 25 10.9905939475982 8.34578335972368 6.17006883864307 90.4771031319024 -23.1011255930613 expand 
#> 26 9.75994290755567 8.57564220737894 6.49398223369463 26.0606430096442 -21.6551083931462 expand 
#> 27 9.75994290755567 8.57564220737894 6.49398223369463 26.0606430096442 -21.6551083931462 reflect 
#> 28 9.75994290755567 8.57564220737894 6.49398223369463 26.0606430096442 -21.6551083931462 contract inside 
#> 29 9.75994290755567 8.57564220737894 6.49398223369463 26.0606430096442 -21.6551083931462 reflect 
#> 30 8.10099790792233 7.09013800751324 8.67472282764109 27.1921182249967 -20.7696634619746 reflect 
#> 31 8.10099790792233 7.09013800751324 8.67472282764109 27.1921182249967 -20.7696634619746 contract inside 
#> 32 8.10099790792233 7.09013800751324 8.67472282764109 27.1921182249967 -20.7696634619746 reflect 
#> 33 8.10099790792233 7.09013800751324 8.67472282764109 27.1921182249967 -20.7696634619746 reflect 
#> 34 8.10099790792233 7.09013800751324 8.67472282764109 27.1921182249967 -20.7696634619746 reflect 
#> 35 7.6296917700292 6.57122597333637 7.92660128042042 32.8090393064141 -19.4764177788075 expand 
#> 36 7.6296917700292 6.57122597333637 7.92660128042042 32.8090393064141 -19.4764177788075 contract inside 
#> 37 7.6296917700292 6.57122597333637 7.92660128042042 32.8090393064141 -19.4764177788075 reflect 
#> 38 7.6296917700292 6.57122597333637 7.92660128042042 32.8090393064141 -19.4764177788075 reflect 
#> 39 7.33360372296142 5.97316058830953 7.68427636024989 118.816682619588 -19.4198789511935 reflect 
#> 40 7.33360372296142 5.97316058830953 7.68427636024989 118.816682619588 -19.4198789511935 reflect 
#> 41 6.82280472299841 5.92679948615982 7.53991257390515 158.727081937072 -19.0941238059158 reflect 
#> 42 6.56993279361279 5.43595841112662 7.97018235660804 38.4164204484102 -18.1404642986635 expand 
#> 43 6.56993279361279 5.43595841112662 7.97018235660804 38.4164204484102 -18.1404642986635 reflect 
#> 44 6.56993279361279 5.43595841112662 7.97018235660804 38.4164204484102 -18.1404642986635 contract inside 
#> 45 6.56993279361279 5.43595841112662 7.97018235660804 38.4164204484102 -18.1404642986635 reflect 
#> 46 6.56993279361279 5.43595841112662 7.97018235660804 38.4164204484102 -18.1404642986635 reflect 
#> 47 6.56993279361279 5.43595841112662 7.97018235660804 38.4164204484102 -18.1404642986635 reflect 
#> 48 5.28548785597528 4.55759418497378 9.04362818071381 21.1941663790621 -17.6156315313864 expand 
#> 49 7.02879026644281 5.68729611088657 5.91795875462147 16.7956725521271 -17.4832388535357 expand 
#> 50 6.0027044935016 5.07393969003788 5.8849948419789 32.6143423991358 -16.7040171530715 expand 
#> 51 6.0027044935016 5.07393969003788 5.8849948419789 32.6143423991358 -16.7040171530715 reflect 
#> 52 6.0027044935016 5.07393969003788 5.8849948419789 32.6143423991358 -16.7040171530715 reflect 
#> 53 6.59731881865314 5.16880316840877 4.59242490608013 16.5576292279009 -16.6550838615707 reflect 
#> 54 4.89250652444853 4.03050834801968 5.62433227241967 20.6451732738895 -16.2362569999038 reflect 
#> 55 4.89250652444853 4.03050834801968 5.62433227241967 20.6451732738895 -16.2362569999038 reflect 
#> 56 4.89250652444853 4.03050834801968 5.62433227241967 20.6451732738895 -16.2362569999038 reflect 
#> 57 4.89250652444853 4.03050834801968 5.62433227241967 20.6451732738895 -16.2362569999038 reflect 
#> 58 4.89250652444853 4.03050834801968 5.62433227241967 20.6451732738895 -16.2362569999038 reflect 
#> 59 4.89250652444853 4.03050834801968 5.62433227241967 20.6451732738895 -16.2362569999038 contract inside 
#> 60 4.89250652444853 4.03050834801968 5.62433227241967 20.6451732738895 -16.2362569999038 contract inside 
#> 61 4.89250652444853 4.03050834801968 5.62433227241967 20.6451732738895 -16.2362569999038 reflect 
#> 62 5.36292098261263 4.29755626460154 4.9432116196164 20.9554132411357 -16.2305970410823 contract outside 
#> 63 5.43119190496216 4.42130450226141 4.63166814110711 24.5958372312959 -16.2221011236137 contract inside 
#> 64 4.7322264841979 3.85798700599787 6.07104455700087 29.2728976901882 -16.2133957895249 reflect 
#> 65 4.91809644345243 3.99201132494822 5.06667446879964 26.1345881287268 -16.1965873713107 contract inside 
#> 66 4.91809644345243 3.99201132494822 5.06667446879964 26.1345881287268 -16.1965873713107 contract inside 
#> 67 4.91809644345243 3.99201132494822 5.06667446879964 26.1345881287268 -16.1965873713107 reflect 
#> 68 5.11626629597035 4.17542026220493 5.03471485009209 25.7817285958635 -16.1950928224326 contract inside 
#> 69 5.06776822078151 4.09944794960308 5.25844675572359 23.3631227319687 -16.1940427432698 contract outside 
#> 70 4.87393343865206 3.96942443425438 5.59519716752651 26.6069763274549 -16.1884001191538 contract inside 
#> 71 4.87393343865206 3.96942443425438 5.59519716752651 26.6069763274549 -16.1884001191538 contract inside 
#> 72 4.87393343865206 3.96942443425438 5.59519716752651 26.6069763274549 -16.1884001191538 contract inside 
#> 73 4.87393343865206 3.96942443425438 5.59519716752651 26.6069763274549 -16.1884001191538 reflect 
#> 74 4.77262513746732 3.89796625572618 5.60040120501363 26.6663710202517 -16.1854157446027 reflect 
#> 75 4.77262513746732 3.89796625572618 5.60040120501363 26.6663710202517 -16.1854157446027 contract inside 
#> 76 4.77262513746732 3.89796625572618 5.60040120501363 26.6663710202517 -16.1854157446027 contract inside 
#> 77 4.91209832956955 3.99632097953225 5.3466965313345 25.5730292012906 -16.1853514137255 contract inside 
#> 78 4.8684042341106 3.96366018361655 5.39995406717573 24.8674024723561 -16.1846893516531 contract outside 
#> 79 4.86939969565057 3.96095533029654 5.48048289360783 25.0778517814404 -16.1846697130798 contract inside 
#> 80 4.86939969565057 3.96095533029654 5.48048289360783 25.0778517814404 -16.1846697130798 reflect 
#> 81 4.81480199115496 3.92563417491044 5.51825896546419 26.0353970092068 -16.1843446039365 contract inside 
#> 82 4.81480199115496 3.92563417491044 5.51825896546419 26.0353970092068 -16.1843446039365 reflect 
#> 83 4.80442051858571 3.91305955878167 5.51631017722255 25.8035376268681 -16.1841907517101 contract inside 
#> 84 4.83960328595387 3.94124851189487 5.4661965593553 25.2304664626269 -16.1841177735429 contract inside 
#> 85 4.83960328595387 3.94124851189487 5.4661965593553 25.2304664626269 -16.1841177735429 contract inside 
#> 86 4.79009685844859 3.90314059407521 5.56113479702454 25.560280848137 -16.1840525040837 contract inside 
#> 87 4.79009685844859 3.90314059407521 5.56113479702454 25.560280848137 -16.1840525040837 contract inside 
#> 88 4.79009685844859 3.90314059407521 5.56113479702454 25.560280848137 -16.1840525040837 contract inside 
#> 89 4.79282808766046 3.90698915212831 5.52269364406131 25.7357051253013 -16.184037588417 reflect 
#> 90 4.8211760660312 3.92740183197262 5.49704836790131 25.4497979455624 -16.1839944791021 contract inside 
#> 91 4.8211760660312 3.92740183197262 5.49704836790131 25.4497979455624 -16.1839944791021 reflect 
#> 92 4.78565141183831 3.90096240033502 5.54238582920928 25.4517694533707 -16.1839799786709 reflect 
#> 93 4.78565141183831 3.90096240033502 5.54238582920928 25.4517694533707 -16.1839799786709 reflect 
#> 94 4.78565141183831 3.90096240033502 5.54238582920928 25.4517694533707 -16.1839799786709 contract inside 
#> 95 4.81265691520433 3.92321817841846 5.49000194047387 25.5498069297059 -16.1839479536819 reflect 
#> 96 4.81265691520433 3.92321817841846 5.49000194047387 25.5498069297059 -16.1839479536819 reflect 
#> 97 4.81265691520433 3.92321817841846 5.49000194047387 25.5498069297059 -16.1839479536819 reflect 
#> 98 4.78987270355293 3.90682029657606 5.55085462395537 25.4411289935749 -16.183893523778 expand 
#> 99 4.80571214245693 3.91934757404329 5.52042309257121 25.6285622976534 -16.1838769884004 reflect 
#> Optimization has terminated successfully. 
# model$ML
# [1] -16.47099