likelihood_MSBD {ML.MSBD}

Likelihood calculation for randomly sampled trees

Description

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.

Usage

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)

Arguments

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.

Details

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}=\gamma*w_{i,j}. If uniform_weights = TRUE, w_{i,j} = \frac{1}{N-1} for all i,j, where N is the total number of states. If uniform_weights = FALSE, w_{i,j} = \frac{p_\lambda(\lambda_j)p_\mu(\mu_j)}{sum_{k \ne i}p_\lambda(\lambda_k)p_\mu(\mu_k)} where the distributions p_\lambda and p_\mu are provided by the inputs p_lambda and p_mu.

Value

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

Examples

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