likelihood_MSBD: Likelihood calculation for randomly sampled trees
In ML.MSBD: Maximum Likelihood Inference on Multi-State Trees

Description Usage Arguments Details Value Examples

View source: R/likelihood_MSBD.R

Calculates the negative log likelihood of a multi-states model given a tree. This function is designed to work with constant extant and/or extinct sampling.

likelihood_MSBD(tree, shifts, gamma, lambdas, mus, lambda_rates = NULL,
  stepsize = NULL, uniform_weights = TRUE, p_lambda = 0, p_mu = 0,
  rho = 1, sigma = 0, rho_sampling = TRUE, add_time = 0,
  unresolved = FALSE)

`tree`	Phylogenetic tree (in ape format) to calculate the likelihood on.
`shifts`	Matrix describing the positions (edges and times) of shifts. See 'Details'.
`gamma`	Rate of state change.
`lambdas`	Birth rates of all states.
`mus`	Death rates of all states.
`lambda_rates`	Rates of decay of birth rate for all states. To use exponential decay, `stepsize` should also be provided.
`stepsize`	Size of the step to use for time discretization with exponential decay, default NULL. To use exponential decay, `lambda_rates` should also be provided.
`uniform_weights`	Whether all states are weighted uniformly in shifts, default TRUE. If FALSE, the weights of states are calculated from the distributions `p_lambda` and `p_mu`. See 'Details'.
`p_lambda`	Prior probability distribution on lambdas, used if `uniform_weights = FALSE`.
`p_mu`	Prior probability distribution on mus, used if `uniform_weights = FALSE`.
`rho`	Sampling proportion on extant tips, default 1.
`sigma`	Sampling probability on extinct tips (tips are sampled upon extinction), default 0.
`rho_sampling`	Whether the most recent tips should be considered extant tips, sampled with sampling proportion `rho`. If FALSE, all tips will be considered extinct tips, sampled with sampling probability `sigma`. Should be TRUE for most macroevolution datasets and FALSE for most epidemiology datasets.
`add_time`	The time between the most recent tip and the end of the process (>=0). This is an internal variable used in calculations for unresolved trees.
`unresolved`	Whether this tree is the backbone of an unresolved tree. This is an internal variable used in calculations for unresolved trees.

It is to be noted that all times are counted backwards, with the most recent tip positioned at 0.

The 'shifts' matrix is composed of 3 columns and a number of rows. Each row describes a shift: the first column is the index of the edge on which the shift happens, the second column is the time of the shift, and the third column is the index of the new state. For example the row vector (3,0.5,2) specifies a shift on edge number 3, at time 0.5, towards the state that has parameters lambdas[2], lambda_rates[2] and mus[2].

The weights w are used for calculating the transition rates q from each state i to j: q(i,j)=γ*w(i,j). If uniform_weights = TRUE, w(i,j)=1/(N-1) for all i,j, where N is the total number of states. If uniform_weights = FALSE, w(i,j)=pλ(λj)pμ(μj)/sum(pλ(λk)pμ(μk)) for all k!=i where the distributions pλ and pμ are provided by the inputs p_lambda and p_mu.

The value of the negative log likelihood of the model given the tree.

# Input a phylogeny
tree <- ape::read.tree(text = "(((t4:0.7293960718,(t1:0.450904974,t3:0.09259337652)
        :0.04068535892):0.4769176776,t8:0.1541864066):0.7282000314,((t7:0.07264320855,
        (((t5:0.8231869878,t6:0.3492440532):0.2380232813,t10:0.2367582193):0.5329497182,
        t9:0.1016243151):0.5929288475):0.3003101915,t2:0.8320755605):0.2918686506);")

# Calculate the log likelihood under a constant birth-death model (i.e, no shifts) 
# with full extant & extinct sampling
likelihood_MSBD(tree, shifts = c(), gamma = 0, lambdas = 10, mus = 1, sigma = 1)
# Calculate the log likelihood under a multi-states model with 2 states 
# and full extant & extinct sampling
likelihood_MSBD(tree, shifts = matrix(c(2,1.8,2), nrow = 1), 
                gamma = 0.05, lambdas = c(10, 6), mus = c(1, 0.5), sigma = 1)
# Calculate the log likelihood under a multi-states model with 2 states and exponential decay 
# with full extant & extinct sampling
likelihood_MSBD(tree, shifts = matrix(c(2,1.8,2), nrow = 1), 
                gamma = 0.05, lambdas = c(10, 6), mus = c(1, 0.5), 
                sigma = 1, stepsize = 0.01, lambda_rates = c(0.1, 0.1))