likelihood_MSBD_unresolved: Likelihood calculation for unresolved trees
In ML.MSBD: Maximum Likelihood Inference on Multi-State Trees

Description Usage Arguments Details Value Examples

View source: R/likelihood_MSBD_unresolved.R

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.

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)

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

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 = "(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)