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: 15.6534. Trying with the next bulk.
#> No convergence yet (steps count: 11000). Gelman-Rubin's R: 35.3876. Trying with the next bulk.
#> No convergence yet (steps count: 16000). Gelman-Rubin's R: 44.6823. Trying with the next bulk.
#> No convergence yet (steps count: 21000). Gelman-Rubin's R: 45.5207. Trying with the next bulk.
#> No convergence yet (steps count: 26000). Gelman-Rubin's R: 42.1064. Trying with the next bulk.
#> No convergence yet (steps count: 31000). Gelman-Rubin's R: 43.9147. Trying with the next bulk.
#> No convergence yet (steps count: 36000). Gelman-Rubin's R: 43.6296. Trying with the next bulk.
#> No convergence yet (steps count: 41000). Gelman-Rubin's R: 39.8554. Trying with the next bulk.
#> No convergence yet (steps count: 46000). Gelman-Rubin's R: 39.8168. Trying with the next bulk.
#> No convergence yet (steps count: 51000). Gelman-Rubin's R: 40.8008. Trying with the next bulk.
#> No convergence yet (steps count: 56000). Gelman-Rubin's R: 40.9831. Trying with the next bulk.
#> No convergence yet (steps count: 61000). Gelman-Rubin's R: 39.6546. Trying with the next bulk.
#> No convergence yet (steps count: 66000). Gelman-Rubin's R: 39.9588. Trying with the next bulk.
#> No convergence yet (steps count: 71000). Gelman-Rubin's R: 41.4681. Trying with the next bulk.
#> No convergence yet (steps count: 76000). Gelman-Rubin's R: 42.2324. Trying with the next bulk.
#> No convergence yet (steps count: 81000). Gelman-Rubin's R: 42.7171. Trying with the next bulk.
#> No convergence yet (steps count: 86000). Gelman-Rubin's R: 42.7701. Trying with the next bulk.
#> No convergence yet (steps count: 91000). Gelman-Rubin's R: 43.7239. Trying with the next bulk.
#> No convergence yet (steps count: 96000). Gelman-Rubin's R: 44.0830. Trying with the next bulk.
#> No convergence yet (steps count: 101000). Gelman-Rubin's R: 43.7994. Trying with the next bulk.
#> No convergence yet (steps count: 106000). Gelman-Rubin's R: 43.2222. Trying with the next bulk.
#> No convergence yet (steps count: 111000). Gelman-Rubin's R: 44.0337. Trying with the next bulk.
#> No convergence yet (steps count: 116000). Gelman-Rubin's R: 44.1949. Trying with the next bulk.
#> No convergence yet (steps count: 121000). Gelman-Rubin's R: 44.8269. Trying with the next bulk.
#> No convergence yet (steps count: 126000). Gelman-Rubin's R: 44.4841. Trying with the next bulk.
#> No convergence yet (steps count: 131000). Gelman-Rubin's R: 44.3707. Trying with the next bulk.
#> No convergence yet (steps count: 136000). Gelman-Rubin's R: 43.7705. Trying with the next bulk.
#> No convergence yet (steps count: 141000). Gelman-Rubin's R: 44.1808. Trying with the next bulk.
#> No convergence yet (steps count: 146000). Gelman-Rubin's R: 43.4634. Trying with the next bulk.
#> No convergence yet (steps count: 151000). Gelman-Rubin's R: 43.3565. Trying with the next bulk.
#> No convergence yet (steps count: 156000). Gelman-Rubin's R: 43.6057. Trying with the next bulk.
#> No convergence yet (steps count: 161000). Gelman-Rubin's R: 43.0680. Trying with the next bulk.
#> No convergence yet (steps count: 166000). Gelman-Rubin's R: 42.8754. Trying with the next bulk.
#> No convergence yet (steps count: 171000). Gelman-Rubin's R: 43.3233. Trying with the next bulk.
#> No convergence yet (steps count: 176000). Gelman-Rubin's R: 43.4753. Trying with the next bulk.
#> No convergence yet (steps count: 181000). Gelman-Rubin's R: 43.2369. Trying with the next bulk.
#> No convergence yet (steps count: 186000). Gelman-Rubin's R: 43.0555. Trying with the next bulk.
#> No convergence yet (steps count: 191000). Gelman-Rubin's R: 43.2616. Trying with the next bulk.
#> No convergence yet (steps count: 196000). Gelman-Rubin's R: 43.2817. Trying with the next bulk.
#> No convergence yet (steps count: 2e+05). Gelman-Rubin's R: 43.4492. Trying with the next bulk.
#> No convergence reached after 2e+05 steps (995 final count of samples).
# }