Model estimation using Markov Chain Monte Carlo

Source: R/aphylo_mcmc.R, R/formulas.R

The function is a wrapper of fmcmc::MCMC().

APHYLO_DEFAULT_MCMC_CONTROL

aphylo_mcmc(
  model,
  params,
  priors = uprior(),
  control = list(),
  check_informative = getOption("aphylo_informative", FALSE),
  reduced_pseq = getOption("aphylo_reduce_pseq", TRUE)
)

APHYLO_PARAM_DEFAULT

Format

An object of class list of length 6.

An object of class numeric of length 9.

Arguments

model

A model as specified in aphylo-model.

params

A vector of length 7 with initial parameters. In particular psi[1], psi[2], mu[1], mu[2], eta[1], eta[2] and Pi.

priors

A function to be used as prior for the model (see bprior).

control

A list with parameters for the optimization method (see details).

check_informative

Logical scalar. When TRUE the algorithm stops with an error when the annotations are uninformative (either 0s or 1s).

reduced_pseq

Logical. When TRUE it will use a reduced peeling sequence in which it drops unannotated leafs. If the model includes eta this is set to FALSE.

Value

An object of class aphylo_estimates.

Details

APHYLO_DEFAULT_MCMC_CONTROL lists the default values for the MCMC estimation:

  • nsteps: 1e4L

  • burnin: 5e3L

  • thin : 10L

  • nchains : 2L

  • multicore : FALSE

  • conv_checker : fmcmc::convergence_auto(5e3)

For more information about the MCMC estimation process, see fmcmc::MCMC().

Methods base::print(), base::summary(), stats::coef, stats::window(), stats::vcov(), stats::logLik(), predict(), and the various ways to query features of the trees via Ntip() are available post estimation.

The vector APHYLO_PARAM_DEFAULT lists the starting values for the parameters in the model. The current defaults are:

  • psi0: 0.10

  • psi1: 0.05

  • mu_d0: 0.90

  • mu_d1: 0.50

  • mu_s0: 0.10

  • mu_s1: 0.05

  • eta0: 1.00

  • eta1: 1.00

  • Pi: 0.50

See also

Other parameter estimation: aphylo_mle()

Examples

# Using the MCMC ------------------------------------------------------------

# \donttest{

set.seed(1233)
# Simulating a tree
tree <- sim_tree(200)

# Simulating functions
atree <- raphylo(
  tree = tree,
  psi  = c(.01, .03),
  mu_d = c(.05, .02),
  Pi   = .5
)

# Running the MCMC
set.seed(1231)

ans_mcmc <- aphylo_mcmc(
  atree ~ mu_d + psi + eta + Pi,
  control = list(nsteps = 2e5, burnin=1000, thin=200)
)
#> Warning: While using multiple chains, a single initial point has been passed via `initial`: c(0.1, 0.05, 0.9, 0.5, 1, 1, 0.5). The values will be recycled. Ideally you would want to start each chain from different locations.
#> No convergence yet (steps count: 6000). Gelman-Rubin's R: 16.3895. Trying with the next bulk.
#> No convergence yet (steps count: 11000). Gelman-Rubin's R: 46.8403. Trying with the next bulk.
#> No convergence yet (steps count: 16000). Gelman-Rubin's R: 41.3444. Trying with the next bulk.
#> No convergence yet (steps count: 21000). Gelman-Rubin's R: 42.6481. Trying with the next bulk.
#> No convergence yet (steps count: 26000). Gelman-Rubin's R: 41.3417. Trying with the next bulk.
#> No convergence yet (steps count: 31000). Gelman-Rubin's R: 39.4372. Trying with the next bulk.
#> No convergence yet (steps count: 36000). Gelman-Rubin's R: 40.7279. Trying with the next bulk.
#> No convergence yet (steps count: 41000). Gelman-Rubin's R: 39.6160. Trying with the next bulk.
#> No convergence yet (steps count: 46000). Gelman-Rubin's R: 40.1712. Trying with the next bulk.
#> No convergence yet (steps count: 51000). Gelman-Rubin's R: 41.0175. Trying with the next bulk.
#> No convergence yet (steps count: 56000). Gelman-Rubin's R: 41.4682. Trying with the next bulk.
#> No convergence yet (steps count: 61000). Gelman-Rubin's R: 41.7795. Trying with the next bulk.
#> No convergence yet (steps count: 66000). Gelman-Rubin's R: 42.4355. Trying with the next bulk.
#> No convergence yet (steps count: 71000). Gelman-Rubin's R: 44.0403. Trying with the next bulk.
#> No convergence yet (steps count: 76000). Gelman-Rubin's R: 45.7104. Trying with the next bulk.
#> No convergence yet (steps count: 81000). Gelman-Rubin's R: 46.1605. Trying with the next bulk.
#> No convergence yet (steps count: 86000). Gelman-Rubin's R: 46.7461. Trying with the next bulk.
#> No convergence yet (steps count: 91000). Gelman-Rubin's R: 47.1607. Trying with the next bulk.
#> No convergence yet (steps count: 96000). Gelman-Rubin's R: 46.7380. Trying with the next bulk.
#> No convergence yet (steps count: 101000). Gelman-Rubin's R: 45.3756. Trying with the next bulk.
#> No convergence yet (steps count: 106000). Gelman-Rubin's R: 45.2792. Trying with the next bulk.
#> No convergence yet (steps count: 111000). Gelman-Rubin's R: 45.7542. Trying with the next bulk.
#> No convergence yet (steps count: 116000). Gelman-Rubin's R: 45.9003. Trying with the next bulk.
#> No convergence yet (steps count: 121000). Gelman-Rubin's R: 45.5821. Trying with the next bulk.
#> No convergence yet (steps count: 126000). Gelman-Rubin's R: 44.3891. Trying with the next bulk.
#> No convergence yet (steps count: 131000). Gelman-Rubin's R: 44.3884. Trying with the next bulk.
#> No convergence yet (steps count: 136000). Gelman-Rubin's R: 43.4947. Trying with the next bulk.
#> No convergence yet (steps count: 141000). Gelman-Rubin's R: 43.2876. Trying with the next bulk.
#> No convergence yet (steps count: 146000). Gelman-Rubin's R: 43.1664. Trying with the next bulk.
#> No convergence yet (steps count: 151000). Gelman-Rubin's R: 43.0121. Trying with the next bulk.
#> No convergence yet (steps count: 156000). Gelman-Rubin's R: 42.7654. Trying with the next bulk.
#> No convergence yet (steps count: 161000). Gelman-Rubin's R: 42.3587. Trying with the next bulk.
#> No convergence yet (steps count: 166000). Gelman-Rubin's R: 42.1333. Trying with the next bulk.
#> No convergence yet (steps count: 171000). Gelman-Rubin's R: 42.3992. Trying with the next bulk.
#> No convergence yet (steps count: 176000). Gelman-Rubin's R: 42.8027. Trying with the next bulk.
#> No convergence yet (steps count: 181000). Gelman-Rubin's R: 42.5090. Trying with the next bulk.
#> No convergence yet (steps count: 186000). Gelman-Rubin's R: 42.5404. Trying with the next bulk.
#> No convergence yet (steps count: 191000). Gelman-Rubin's R: 42.7285. Trying with the next bulk.
#> No convergence yet (steps count: 196000). Gelman-Rubin's R: 43.2007. Trying with the next bulk.
#> No convergence yet (steps count: 2e+05). Gelman-Rubin's R: 42.9600. Trying with the next bulk.
#> No convergence reached after 2e+05 steps (995 final count of samples).
# }