Likelihood for SecSSE model, using Rcpp Loglikelihood calculation for the cla_SecSSE model given a set of parameters and data using Rcpp
Source:R/secsse_loglik.R
cla_secsse_loglik.RdLikelihood for SecSSE model, using Rcpp Loglikelihood calculation for the cla_SecSSE model given a set of parameters and data using Rcpp
Usage
cla_secsse_loglik(
parameter,
phy,
traits,
num_concealed_states,
cond = "proper_cond",
root_state_weight = "proper_weights",
sampling_fraction,
setting_calculation = NULL,
see_ancestral_states = FALSE,
loglik_penalty = 0,
is_complete_tree = FALSE,
take_into_account_root_edge = FALSE,
num_threads = 1,
method = "odeint::bulirsch_stoer",
atol = 1e-08,
rtol = 1e-07,
display_warning = TRUE,
use_normalization = TRUE
)Arguments
- parameter
list where first vector represents lambdas, the second mus and the third transition rates.
- phy
phylogenetic tree of class
phylo, rooted and with branch lengths. Alternatively, multiple phylogenetic trees can be provided as themultiPhyloclass.- 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 amultiPhyloset 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.
- 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 vectorc(1, 0, 0)indicates state 1 was the root state. When using amultiPhyloobject, 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
multiPhyloobject, sampling fraction should be list where each entry in the list corresponds to the sampling proportion for each tree.- setting_calculation
argument used internally to speed up calculation. It should be left blank (default :
setting_calculation = NULL).- see_ancestral_states
Boolean for whether the ancestral states should be shown? Defaults to
FALSE.- 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 toFALSE.- 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)
- num_threads
number of threads to be used. Default is one thread.
- 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".- atol
A numeric specifying the absolute tolerance of integration.
- rtol
A numeric specifying the relative tolerance of integration.
- display_warning
display a warning if necessary
- use_normalization
normalize the density vector during integration, more accurate but slower (default = TRUE)
Examples
rm(list=ls(all=TRUE))
library(secsse)
set.seed(13)
phylotree <- ape::rcoal(12, tip.label = 1:12)
traits <- sample(c(0,1,2),ape::Ntip(phylotree),replace=TRUE)
num_concealed_states <- 3
sampling_fraction <- c(1,1,1)
phy <- phylotree
# the idparlist for a ETD model (dual state inheritance model of evolution)
# would be set like this:
idparlist <- cla_id_paramPos(traits,num_concealed_states)
lambd_and_modeSpe <- idparlist$lambdas
lambd_and_modeSpe[1,] <- c(1,1,1,2,2,2,3,3,3)
idparlist[[1]] <- lambd_and_modeSpe
idparlist[[2]][] <- 0
masterBlock <- matrix(4,ncol=3,nrow=3,byrow=TRUE)
diag(masterBlock) <- NA
idparlist [[3]] <- q_doubletrans(traits,masterBlock,diff.conceal = FALSE)
# Now, internally, clasecsse sorts the lambda matrices, so they look like:
prepare_full_lambdas(traits,num_concealed_states,idparlist[[1]])
#> [[1]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 1 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0
#> [6,] 0 0 0 0 0 0 0 0 0
#> [7,] 0 0 0 0 0 0 0 0 0
#> [8,] 0 0 0 0 0 0 0 0 0
#> [9,] 0 0 0 0 0 0 0 0 0
#>
#> [[2]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0 0 0 0 0 0 0 0 0
#> [2,] 0 1 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0
#> [6,] 0 0 0 0 0 0 0 0 0
#> [7,] 0 0 0 0 0 0 0 0 0
#> [8,] 0 0 0 0 0 0 0 0 0
#> [9,] 0 0 0 0 0 0 0 0 0
#>
#> [[3]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 1 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0
#> [6,] 0 0 0 0 0 0 0 0 0
#> [7,] 0 0 0 0 0 0 0 0 0
#> [8,] 0 0 0 0 0 0 0 0 0
#> [9,] 0 0 0 0 0 0 0 0 0
#>
#> [[4]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 2 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0
#> [6,] 0 0 0 0 0 0 0 0 0
#> [7,] 0 0 0 0 0 0 0 0 0
#> [8,] 0 0 0 0 0 0 0 0 0
#> [9,] 0 0 0 0 0 0 0 0 0
#>
#> [[5]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 2 0 0 0 0
#> [6,] 0 0 0 0 0 0 0 0 0
#> [7,] 0 0 0 0 0 0 0 0 0
#> [8,] 0 0 0 0 0 0 0 0 0
#> [9,] 0 0 0 0 0 0 0 0 0
#>
#> [[6]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0
#> [6,] 0 0 0 0 0 2 0 0 0
#> [7,] 0 0 0 0 0 0 0 0 0
#> [8,] 0 0 0 0 0 0 0 0 0
#> [9,] 0 0 0 0 0 0 0 0 0
#>
#> [[7]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0
#> [6,] 0 0 0 0 0 0 0 0 0
#> [7,] 0 0 0 0 0 0 3 0 0
#> [8,] 0 0 0 0 0 0 0 0 0
#> [9,] 0 0 0 0 0 0 0 0 0
#>
#> [[8]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0
#> [6,] 0 0 0 0 0 0 0 0 0
#> [7,] 0 0 0 0 0 0 0 0 0
#> [8,] 0 0 0 0 0 0 0 3 0
#> [9,] 0 0 0 0 0 0 0 0 0
#>
#> [[9]]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
#> [1,] 0 0 0 0 0 0 0 0 0
#> [2,] 0 0 0 0 0 0 0 0 0
#> [3,] 0 0 0 0 0 0 0 0 0
#> [4,] 0 0 0 0 0 0 0 0 0
#> [5,] 0 0 0 0 0 0 0 0 0
#> [6,] 0 0 0 0 0 0 0 0 0
#> [7,] 0 0 0 0 0 0 0 0 0
#> [8,] 0 0 0 0 0 0 0 0 0
#> [9,] 0 0 0 0 0 0 0 0 3
#>
# which is a list with 9 matrices, corresponding to the 9 states
# (0A,1A,2A,0B,etc)
# if we want to calculate a single likelihood:
parameter <- idparlist
lambda_and_modeSpe <- parameter$lambdas
lambda_and_modeSpe[1,] <- c(0.2,0.2,0.2,0.4,0.4,0.4,0.01,0.01,0.01)
parameter[[1]] <- prepare_full_lambdas(traits,num_concealed_states,
lambda_and_modeSpe)
parameter[[2]] <- rep(0,9)
masterBlock <- matrix(0.07, ncol=3, nrow=3, byrow=TRUE)
diag(masterBlock) <- NA
parameter [[3]] <- q_doubletrans(traits,masterBlock,diff.conceal = FALSE)
cla_secsse_loglik(parameter, phy, traits, num_concealed_states,
cond = 'maddison_cond',
root_state_weight = 'maddison_weights', sampling_fraction,
setting_calculation = NULL,
see_ancestral_states = FALSE,
loglik_penalty = 0)
#> Warning: Deduced names and order of used states to be: 1, 2, 0
#> if this is incorrect, consider passing states as matching numeric
#> ordering, e.g. 1 for the first state, 2 for the second etc.
#> [1] -42.18407
# LL = -42.18407