cla_secsse_ml: Maximum likehood estimation for (SecSSE)

Description Usage Arguments Value Note Examples

View source: R/cla_secsse_ml.R

Description

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

Usage

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
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)),
  use_fortran = TRUE,
  methode = "ode45",
  optimmethod = "simplex",
  num_cycles = 1,
  run_parallel = FALSE,
  loglik_penalty = 0,
  is_complete_tree = FALSE,
  func = if (is_complete_tree) "cla_secsse_runmod_ct_d" else "cla_secsse_runmod",
  verbose = (optimmethod == "subplex")
)

Arguments

phy

phylogenetic tree of class phylo, ultrametric, rooted and with branch lengths.

traits

a vector with trait states for each tip in the phylogeny.

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

id of parameters to be estimated.

initparsopt

initial guess of the parameters to be estimated.

idparsfix

id of the fixed parameters.

parsfix

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 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

maximum tolerance. Default is "c(1e-04, 1e-05, 1e-05)".

maxiter

max number of iterations. Default is "1000 *round((1.25)^length(idparsopt))".

use_fortran

Should the Fortran code for numerical integration be called? Default is TRUE.

methode

method used for integration calculation. Default is "ode45".

optimmethod

method used for optimization. Default is "simplex".

num_cycles

number of cycles of the optimization (default is 1).

run_parallel

should the routine to run in parallel be called?

loglik_penalty

the size of the penalty for all parameters; default is 0 (no penalty)

is_complete_tree

whether or not a tree with all its extinct species is provided

func

function to be used in solving the ODE system. Currently only for testing purposes.

verbose

sets verbose output; default is verbose when optimmethod is 'submplex'

Value

Parameter estimated and maximum likelihood

Note

To run in parallel it is needed to load the following libraries when windows: apTreeshape, doparallel and foreach. When unix, it requires: apTreeshape, doparallel, foreach and doMC

Examples

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
# 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)
traits <-  sample(c(0,1,2), ape::Ntip(phylotree),replace=TRUE) #get some traits
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))
use_fortran <- FALSE
methode <- "ode45"
optimmethod <- "simplex"
run_parallel <- FALSE
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,
#  use_fortran,
#  methode,
#  optimmethod,
#  num_cycles = 1,
#  run_parallel)
# [1] -90.97626

secsse documentation built on July 16, 2021, 9:06 a.m.