Fossilized birth-death diffusion process

This model follows from the episodic fossilized birth-death model, where we relax the assumption of constant rates for instead allowing the per-lineage instantaneous speciation rates $\lambda(t)$ and extinction rates $\mu(t)$ to follow separate Geometric Brownian motions (GBM), such that

\[d(\text{ln}(\lambda(t)) = \alpha_{\lambda} dt + \sigma_{\lambda} d W(t), \\ d(\text{ln}(\mu(t)) = \alpha_{\mu} dt + \sigma_{\mu} d W(t),\]

where $W(t)$ is the Wiener process (i.e., standard Brownian motion), $\alpha_{\lambda}$ and $\alpha_{\mu}$ are the drift in speciation and extinction, and $\sigma_{\lambda}$ and $\sigma_{\mu}$ are the diffusion for speciation and extinction, respectively.

Simulations

One can specify a time or a number of lineages to simulate under the FBDD. For example, for $50$ species

sim_gbmfbd(50, λ0 = 1.0, μ0 = 0.9, αλ = 0.0, αμ = 0.0, σλ = 0.1, σμ = 0.1),

Or one can simulate for a given amount of time, say, $10$ time units

sim_gbmfbd(10.0, 1.0, 0.9, 0.0, 0.0, 0.1, 0.1),

Inference

One can perform inference using:

r, tv = insane_gbmfbd(tree,
                      nburn    = 1_000,
                      niter    = 50_000,
                      nthin    = 50, 
                      ofile    = "<directory>",
                      λ0_prior = (0.05, 148.41),
                      μ0_prior = (0.05, 148.41),
                      αλ_prior = (0.0, 1.0),
                      αμ_prior = (0.0, 1.0),
                      σλ_prior = (0.05, 0.05),
                      σμ_prior = (3.0, 0.1),
                      tρ       = Dict("" => 1.0))

where we have log-normal priors on the initial (root) values for $\lambda_0$ and $\mu_0$, normal priors on the drift for speciation and extinction, and Inverse Gamma prior for the speciation and extinction diffusion rates. Usually, an informative prior such as σμ_prior = (3.0, 0.1) or more is often needed.

If additional occurrences exist, they can be added as specified in Adding external fossil occurrences

Full documentation

Tapestree.INSANE.sim_gbmfbdFunction
sim_gbmfbd(n       ::Int64;
           λ0      ::Float64         = 1.0,
           μ0      ::Float64         = 0.2,
           αλ      ::Float64         = 0.0,
           αμ      ::Float64         = 0.0,
           σλ      ::Float64         = 0.1,
           σμ      ::Float64         = 0.1,
           ψ       ::Vector{Float64} = [0.1],
           ψts     ::Vector{Float64} = Float64[],
           init    ::Symbol          = :stem,
           δt      ::Float64         = 1e-3,
           nstar   ::Int64           = 2*n,
           p       ::Float64         = 5.0,
           warnings::Bool            = true,
           maxt    ::Float64         = δt*1e6)

Simulate iTfbd according to geometric Brownian motions for birth and death rates. Note that it will not necessarily evolve along all of the epochs.

source
sim_gbmfbd(t   ::Float64;
           λ0  ::Float64         = 1.0,
           μ0  ::Float64         = 0.2,
           αλ  ::Float64         = 0.0,
           αμ  ::Float64         = 0.0,
           σλ  ::Float64         = 0.1,
           σμ  ::Float64         = 0.1,
           ψ   ::Vector{Float64} = [0.1],
           ψts ::Vector{Float64} = Float64[],
           δt  ::Float64         = 1e-3,
           nlim::Int64           = 10_000,
           init::Symbol          = :crown)

Simulate iTfbd according to geometric Brownian motions for birth and death rates.

source
Tapestree.INSANE.insane_gbmfbdFunction
insane_gbmfbd(tree    ::sTf_label;
              λ0_prior::NTuple{2,Float64}     = (0.1, 148.41),
              μ0_prior::NTuple{2,Float64}     = (0.1, 148.41),
              αλ_prior::NTuple{2,Float64}     = (0.0, 1.0),
              αμ_prior::NTuple{2,Float64}     = (0.0, 1.0),
              σλ_prior::NTuple{2,Float64}     = (0.05, 0.05),
              σμ_prior::NTuple{2,Float64}     = (3.0, 0.1),
              ψ_prior ::NTuple{2,Float64}     = (1.0, 1.0),
              ψ_epoch ::Vector{Float64}       = Float64[],
              f_epoch ::Vector{Int64}         = Int64[0],
              niter   ::Int64                 = 1_000,
              nthin   ::Int64                 = 10,
              nburn   ::Int64                 = 200,
              nflush  ::Int64                 = nthin,
              ofile   ::String                = string(homedir(), "/fbdd"),
              ϵi      ::Float64               = 0.2,
              λi      ::Float64               = NaN,
              μi      ::Float64               = NaN,
              ψi      ::Float64               = NaN,
              αλi     ::Float64               = 0.0,
              αμi     ::Float64               = 0.0,
              σλi     ::Float64               = 1e-3,
              σμi     ::Float64               = 1e-3,
              pupdp   ::NTuple{7,Float64}     = (1e-3, 1e-3, 1e-3, 1e-4, 1e-4, 0.1, 0.2),
              δt      ::Float64               = 1e-3,
              survival::Bool                  = true,
              mxthf   ::Float64               = 0.1,
              prints  ::Int64                 = 5,
              stnλ    ::Float64               = 0.5,
              stnμ    ::Float64               = 0.5,
              tρ      ::Dict{String, Float64} = Dict("" => 1.0))

Run insane for fossilized birth-death diffusion fbdd.

source