A Markov Chain Monte Carlo posterior simulation should be visually inspected to assess convergence.
# load CNPBayes suppressMessages(library(CNPBayes)) # load packages for manipulating and visualizing data suppressMessages(library(dplyr)) suppressMessages(library(tidyr)) suppressMessages(library(ggplot2))
The number of starting values, burnin MCMC simulatioms, and MCMC simulations
after burnin are controlled by a
McmcParams object. Here, we specify a small
number so that the example runs quickly. As discussed in overview, four types of models are possible--SB, MB, SBP, and MBP. Below we fit only the SBP model with 2 components.
set.seed(1) N <- 200 n <- 81 lrr <- c(rnorm(100, -0.5, sd=0.1), rnorm(100, 0, sd=0.1)) mp <- McmcParams(iter=600, burnin=10, nStarts=4) mlist <- gibbs(model="SBP", mp=mp, dat=lrr, k_range=c(2, 2))
Had more than 1 model been fit, the models in
mlist would be sorted by decreasing value of the marginal likelihood. Note, the marginal likelihood is only meaningful if the models have converged. By default, the function
gibbs requires that the effective sample size (number of independent MCMC draws) for all parameter chains is at least 500 and the multivariate Gelman Rubin diagnostic to be less than 1.2. Otherwise, starting values for 4 (
nStarts) new models will be independently initialized and the burnin and number of
thin iterations will be increased by a factor of 2. This process is repeated until the maximum number of burnin iterations has been reached (see
max_burnin) or the effective size and Gelman Rubin criteria are both satisfied. For each model returned, the
nStarts chains will be combined. In addition to the above automated diagnostics, we suggest plotting the chains as illustrated below.
model <- mlist[] figs <- ggChains(model) figs[] figs[]
The marginal likelihood of the models can be retrieved using the
m.lik <- round(sapply(mlist, marginal_lik), 1) m.lik
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.