inst/doc/MCMCvis.R

In MCMCvis: Tools to Visualize, Manipulate, and Summarize MCMC Output

## ---- eval = FALSE------------------------------------------------------------
#  library(rjags)
#  
#  # create JAGS model
#  mf <- "
#  model {
#  for (i in 1:10)
#  {
#    y[i] ~ dnorm(mu, 0.01);
#  }
#  mu ~ dnorm(0, 0.01)
#  }
#  "
#  
#  data <- list(y = rnorm(10))
#  
#  jm <- rjags::jags.model(
#    textConnection(mf),
#    data = data,
#    n.chains = 3)
#  
#  jags_out <- rjags::coda.samples(
#    jm,
#    variable.names = 'mu',
#    n.iter = 500)

## -----------------------------------------------------------------------------
library(MCMCvis)

## ---- eval = FALSE------------------------------------------------------------
#  MCMCsummary(jags_out, round = 2)

## ---- eval = FALSE------------------------------------------------------------
#  ##     mean   sd  2.5%   50% 97.5% Rhat n.eff
#  ## mu -0.28 2.97 -6.13 -0.14  5.22    1  1397

## ---- eval = FALSE------------------------------------------------------------
#  library(nimble)
#  
#  # create NIMBLE model
#  mf2 <- nimbleCode({
#    for (i in 1:10) {
#      y[i] ~ dnorm(mu, 0.01);
#    }
#    mu ~ dnorm(0, 0.01)
#  })
#  
#  nimble_out <- nimbleMCMC(
#    code = mf2,
#    constants = list(N = 10),
#    data = data,
#    inits = list(mu = 0),
#    nchains = 4,
#    niter = 500)

## ---- eval = FALSE------------------------------------------------------------
#  MCMCsummary(nimble_out, round = 2)

## ---- eval = FALSE------------------------------------------------------------
#  ##     mean   sd  2.5%   50% 97.5% Rhat n.eff
#  ## mu -0.03 2.99 -5.76 -0.1  5.88    1  2000

## ---- eval = FALSE------------------------------------------------------------
#  library(rstan)
#  
#  # create Stan model
#  sm <- "
#  data {
#    real y[10];
#  }
#  parameters {
#    real mu;
#  }
#  model {
#    for (i in 1:10) {
#      y[i] ~ normal(mu, 10);
#    }
#    mu ~ normal(0, 10);
#  }
#  "
#  
#  stan_out <- stan(
#    model_code = sm,
#    data = data,
#    iter = 500)

## ---- eval = FALSE------------------------------------------------------------
#  MCMCsummary(stan_out, round = 2)

## ---- eval = FALSE------------------------------------------------------------
#  ##       mean   sd  2.5%   50% 97.5% Rhat n.eff
#  ## mu   -0.51 2.82 -6.06 -0.36  5.07 1.01   414
#  ## lp__ -0.45 0.61 -2.27 -0.23 -0.01 1.00   508

## ---- message = FALSE---------------------------------------------------------
MCMCsummary(MCMC_data)

## ---- message = FALSE---------------------------------------------------------
MCMCsummary(MCMC_data, round = 2)

## -----------------------------------------------------------------------------
MCMCsummary(MCMC_data, 
            params = 'alpha', 
            round = 2)

## -----------------------------------------------------------------------------
MCMCsummary(MCMC_data, 
            params = 'alpha[1]', 
            ISB = FALSE, 
            exact = TRUE, 
            round = 2)

## -----------------------------------------------------------------------------
MCMCsummary(MCMC_data, 
            params = 'alpha\\[1\\]', 
            ISB = FALSE, 
            exact = FALSE, 
            round = 2)

## -----------------------------------------------------------------------------
MCMCsummary(MCMC_data, 
            params = 'all', 
            excl = 'alpha', 
            ISB = FALSE, 
            exact = FALSE, 
            round = 2)

## -----------------------------------------------------------------------------
MCMCsummary(MCMC_data, 
            params = 'alpha', 
            Rhat = TRUE, 
            n.eff = TRUE, 
            round = 2)

## -----------------------------------------------------------------------------
MCMCsummary(MCMC_data, 
            params = 'alpha', 
            Rhat = TRUE, 
            n.eff = TRUE, 
            probs = c(0.1, 0.5, 0.9), 
            round = 2)

## -----------------------------------------------------------------------------
MCMCsummary(MCMC_data, 
            params = 'alpha', 
            Rhat = TRUE, 
            n.eff = TRUE, 
            HPD = TRUE, 
            hpd_prob = 0.8,
            round = 2)

## -----------------------------------------------------------------------------
MCMCsummary(MCMC_data, 
            params = 'alpha', 
            Rhat = TRUE, 
            n.eff = TRUE, 
            round = 2, 
            func = function(x) ecdf(x)(-10), func_name = "ecdf-10",
            pg0 = TRUE)

## -----------------------------------------------------------------------------
MCMCpstr(MCMC_data, 
         params = 'alpha', 
         func = mean,
         type = 'summary')

## -----------------------------------------------------------------------------
MCMCpstr(MCMC_data, func = function(x) quantile(x, probs = c(0.01, 0.99)))

## -----------------------------------------------------------------------------
ex <- MCMCpstr(MCMC_data, type = 'chains')
dim(ex$alpha)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCtrace(MCMC_data, 
          params = c('beta[1]', 'beta[2]', 'beta[3]'), 
          ISB = FALSE, 
          exact = TRUE,
          pdf = FALSE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCtrace(MCMC_data, 
          params = 'beta', 
          type = 'density', 
          ind = TRUE, 
          pdf = FALSE)

## ---- eval = FALSE, fig.align = 'center'--------------------------------------
#  MCMCtrace(MCMC_data,
#            pdf = TRUE,
#            open_pdf = FALSE,
#            filename = 'MYpdf',
#            wd = 'DIRECTORY_HERE')

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCtrace(MCMC_data, 
          params = c('beta[1]', 'beta[2]', 'beta[3]'), 
          ISB = FALSE, 
          exact = TRUE, 
          iter = 100, 
          ind = TRUE, 
          pdf = FALSE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
# note that the same prior used for all parameters
# the following prior is equivalent to dnorm(0, 0.001) in JAGS
PR <- rnorm(15000, 0, 32)

MCMCtrace(MCMC_data, 
          params = c('beta[1]', 'beta[2]', 'beta[3]'),
          ISB = FALSE,
          exact = TRUE,
          priors = PR,
          pdf = FALSE,
          Rhat = TRUE,
          n.eff = TRUE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCtrace(MCMC_data, 
          params = c('beta[1]', 'beta[2]', 'beta[3]'),
          ISB = FALSE,
          exact = TRUE,
          priors = PR,
          pdf = FALSE,
          post_zm = FALSE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
PPO <- MCMCtrace(MCMC_data,
                 params = c('beta[1]', 'beta[2]', 'beta[3]'),
                 ISB = FALSE,
                 exact = TRUE,
                 priors = PR,
                 plot = FALSE,
                 PPO_out = TRUE)
PPO

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCtrace(MCMC_data,
          params = c('beta[1]', 'beta[2]', 'beta[3]'),
          ISB = FALSE,
          exact = TRUE,
          priors = PR,
          pdf = FALSE,
          Rhat = TRUE,
          n.eff = TRUE,
          xlab_tr = 'This is the x for trace',
          ylab_tr = 'This is the y for trace',
          main_den = 'Custom density title',
          lwd_den = 3,
          lty_pr = 2,
          col_pr = 'green',
          sz_txt = 1.3,
          sz_ax = 2,
          sz_ax_txt = 1.2,
          sz_tick_txt = 1.2,
          sz_main_txt = 1.3)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
# generating values for each parameter used to simulate data
GV <- c(-10, -5.5, -15)

MCMCtrace(MCMC_data,
          params = c('beta[1]', 'beta[2]', 'beta[3]'),
          ISB = FALSE,
          exact = TRUE,
          gvals = GV,
          pdf = FALSE)

## -----------------------------------------------------------------------------
ex <- MCMCchains(MCMC_data, params = 'beta')

## -----------------------------------------------------------------------------
apply(ex, 2, mean)

## -----------------------------------------------------------------------------
ex2 <- MCMCchains(MCMC_data, 
                  params = 'beta',  
                  mcmc.list = TRUE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCplot(MCMC_data, 
         params = 'beta', 
         ci = c(50, 90))

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCplot(MCMC_data, 
         params = 'beta', 
         ci = c(50, 80), 
         HPD = TRUE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCplot(MCMC_data, 
         params = 'beta', 
         ref_ovl = TRUE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCplot(MCMC_data, 
         params = 'beta', 
         rank = TRUE, 
         xlab = 'PARAMETER ESTIMATE', 
         guide_lines = TRUE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCplot(MCMC_data, 
         params = 'beta', 
         rank = TRUE, 
         horiz = FALSE, 
         ylab = 'PARAMETER ESTIMATE')

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCplot(object = MCMC_data, 
         object2 = MCMC_data2, 
         params = 'beta', 
         offset = 0.1,
         ref_ovl = TRUE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCplot(MCMC_data, 
         params = 'beta', 
         xlim = c(-60, 40),
         xlab = 'My x-axis label',
         main = 'MCMCvis plot',
         labels = c('First param', 'Second param', 'Third param', 
                    'Fourth param', 'Fifth param', 'Sixth param'),
          col = c('red', 'blue', 'green', 'purple', 'orange', 'black'),
          sz_labels = 1.5,
          sz_med = 2,
          sz_thick = 7,
          sz_thin = 3,
          sz_ax = 4,
          sz_main_txt = 2)

## -----------------------------------------------------------------------------
# note that the same prior used for all parameters
# the following prior is equivalent to dnorm(0, 0.001) in JAGS
PR <- rnorm(15000, 0, 32)

PPO <- MCMCtrace(MCMC_data,
                 params = c('beta[1]', 'beta[2]', 'beta[3]'),
                 ISB = FALSE,
                 exact = TRUE,
                 priors = PR,
                 plot = FALSE,
                 PPO_out = TRUE)

## ---- fig.width = 7, fig.height = 5, fig.align = 'center'---------------------
MCMCplot(MCMC_data, 
         params = c('beta[1]', 'beta[2]', 'beta[3]'),
          ISB = FALSE,
          exact = TRUE,
          xlim = c(-60, 35))
        
# each parameter is a y-unit of 1
for (i in 1:NROW(PPO)) {
  text(x = 10, y = NROW(PPO) - i + 1, paste0('PPO: ', PPO[i,2], '%'), pos = 4, col = 'red')
}

## ---- eval = FALSE------------------------------------------------------------
#  MCMCdiag(model_fit,
#           round = 3,
#           file_name = 'model-summary-YYYY-MM-DD',
#           dir = 'Results',
#           mkdir = 'model-YYYY-MM-DD',
#           add_field = '1.0',
#           add_field_names = 'Data version',
#           save_obj = TRUE,
#           obj_name = 'model-fit-YYYY-MM-DD',
#           add_obj = list(DATA, sessionInfo()),
#           add_obj_names = c('model-data-YYYY-MM-DD', 'session-info-YYYY-MM-DD'),
#           cp_file = c('model.stan', 'fit-model.R'),
#           cp_file_names = c('model-YYYY-MM-DD.stan, fit-model-YYYY-MM-DD.R'))
#