likelihood_MSBD_unresolved {ML.MSBD}

Likelihood calculation for unresolved trees

Description

Calculates the negative log likelihood of a multi-states model given a tree. This function is designed to work with unresolved trees, where tips represent collapsed clades. This sampling scheme is not recommended for epidemiology datasets. The MRCA times of collapsed clades and the number of collapsed lineages need to be provided for all tips. If neither is provided the function will default to random sampling. Extinct tips can be present outside of the unresolved parts, but not below the time(s) set for tcut.

Usage

likelihood_MSBD_unresolved(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, lineage_counts = c(), tcut = NULL)

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.

lineage_counts

Number of lineages collapsed on each tip. Should be set to 1 for extinct tips.

tcut

Times of clade collapsing for each tip (i.e time of the MRCA of all collapsed lineages). Can be a single number or a vector of length the number of tips.

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 = "(t3:0.9703302342,((t4:0.1999577823,(t2:0.1287530271,
        (t7:0.08853561159,(t8:0.07930237712,t9:0.07930237712):0.009233234474):0.04021741549):
        0.07120475526):0.4269919425,(((t10:0.0191876225,t5:0.0191876225):0.04849906822,
        t6:0.06768669072):0.1672340445,t1:0.2349207353):0.3920289896):0.3433805094);")

# Calculate the log likelihood under a constant birth-death model (i.e, no shifts) 
# with unresolved tips
likelihood_MSBD_unresolved(tree, shifts = c(), gamma = 0, lambdas = 10, mus = 1, 
                           lineage_counts = c(2,5,1,3,1,1,1,1,2,6), tcut = 0.05)
# Calculate the log likelihood under a multi-states model with 2 states and unresolved tips
likelihood_MSBD_unresolved(tree, shifts = matrix(c(2,0.7,2), nrow = 1), 
                           gamma = 0.05, lambdas = c(10, 5), mus = c(1, 1), 
                           lineage_counts = c(2,5,1,3,1,1,1,1,2,6), tcut = 0.05)