In Giappo/sls: Likelihood model for diversification processes with single lineage rate shifts

Legenda:

This is the pipeline to follow for the sls project. The aim of the project is to correctly estimate the likelihood for a process with a rate shift occurring on a single lineage at a given time.

Load the library:

library(sls, quietly = TRUE)

Step 1: Simulations

First we simulate the process

lambdas <- c(0.3, 0.6)
mus <- c(0.2, 0.1)
cond <- 3
seed <- 42
sim_pars <- c(lambdas[1], mus[1], lambdas[2], mus[2])
l_2 <- sls::sim_get_standard_l_2(crown_age = 5, shift_time = 2)
set.seed(seed)
sim <- sls::sls_sim(
  lambdas = lambdas,
  mus = mus,
  cond = cond,
  l_2 = l_2
)
sim

We can also plot both the main clade and the subclade

par(mfrow = c(1, 2))
ape::plot.phylo(
  DDD::L2phylo(sim$l_tables[[1]], dropextinct = FALSE)
)
ape::plot.phylo(
  DDD::L2phylo(sim$l_tables[[2]], dropextinct = FALSE)
)

Step 2: calculate the likelihood

We can calculate the (log) likelihood function for this main and sub clade

pars_m <- c(0.2, 0.1)
pars_s <- c(0.5, 0.4)
brts_m <- sim$brts[[1]]
brts_s <- sim$brts[[2]]
pars <- c(pars_m, pars_s)
brts <- list(brts_m, brts_s)

loglik <- sls::loglik_sls_p(
  pars = pars,
  brts = brts,
  cond = cond
)
loglik

We have many ways to calculate the same likelihood. For example we can integrate the Q-equation

loglik <- sls::loglik_sls_q(
  pars = pars,
  brts = brts,
  cond = cond
)
loglik

The idea of the project is to provide a correct likelihood as opposed to the ones currently available to solve these kind of problems, namely calculating the likelihood for a tree where a rate shift occurs on a specific lineage at some point in time. You can verify that the likelihood functions provided by older methods do actually yield different values for the same tree.

# Nee et al.
loglik <- sls::loglik_sls_q_nodiv(
  pars = pars,
  brts = brts,
  cond = cond
)
loglik

# Bisse
loglik <- sls::loglik_bisse_shift(
  pars = pars,
  brts = brts,
  cond = cond
)
loglik

# DDD
loglik <- sls::loglik_ddd(
  pars = pars,
  brts = brts,
  cond = cond
)
loglik

Step 3: maximizing the likelihood

The final goal is to illustrate how the new likelihood can improve the parameter inference using maximum likelihood

# This might take some time to run
if (1 == 2) {
  mle <- sls::sls_ml(
    loglik_function = loglik_sls_p,
    brts = brts,
    cond = cond,
    verbose = FALSE
  )
  mle
}

Step 4: simplifying the process

The entire routine can be easily executed using the function sls_main. This function can accommodate for the use of multiple likelihood functions. For example here we can directly compare the old Nee likelihood with the new sls one.

# This might take some time to run
if (1 == 2) {
  mle <- sls::sls_main(
    sim_pars = sim_pars,
    l_2 = l_2,
    seed = seed,
    loglik_functions = c(loglik_sls_p, loglik_sls_p_nodiv),
    cond = cond,
    verbose = FALSE
  )
  mle
}

Giappo/sls documentation built on Feb. 1, 2021, 9:55 a.m.