runjags (version 0.9.0-4)

autorun.jags: Run a User Specified Bayesian MCMC Model in JAGS with Automatically Calculated Run Length and Convergence Diagnostics

Description

Runs a user specified JAGS (similar to WinBUGS) model from within R, returning a list of the MCMC chain(s) along with convergence diagnostics, autocorrelation diagnostics and monitored variable summaries. Chain convergence over the first run of the simulation is assessed using the Gelman and Rubin's convergence diagnostic. If necessary, the simulation is extended to improve chain convergence (up to a user-specified maximum time limit), before the required sample size of the Markov chain is calculated using the Raftery and Lewis's diagnostic. The simulation is extended to the required sample size dependant on autocorrelation. This function is provided as an educational tool only, and is not a replacement for manually assessing convergence and monte carlo error for real-world applications. JAGS is called using the lower level function run.jags.

Usage

autorun.jags(model=stop("No model supplied"), 
   monitor = stop("No monitored variables supplied"), 
   data=NA, n.chains=2, inits = replicate(n.chains, NA), 
   startburnin = 5000, startsample = 10000,
   psrf.target = 1.05, crash.retry = 1,
   thin.sample = FALSE, jags = findjags(), 
   silent.jags = FALSE, interactive=TRUE, max.time=Inf, 
   adaptive=list(type="burnin", length=200))

Arguments

model
a character string of the model in the JAGS language. No default.
monitor
a character vector of the names of variables to monitor. No default.
data
either a named list or a character string in the R dump format containing the data. If left as NA, the model will be run without external data.
n.chains
the number of chains to use with the simulation. More chains will improve the sensitivity of the convergence diagnostic, but will cause the simulation to run more slowly. The minimum (and default) number of chains is 2.
inits
a character vector with length equal to either n.chains or length 1. Each element of the vector must be a character string in the R dump format representing the initial values for that chain. If inits is of length 1, all chains will be run using the sam
startburnin
the number of initial updates to discard before sampling. Only used on the initial run before checking convergence. Default 5000 iterations.
startsample
the number of samples on which to assess convergence. More samples will give a better chance of allowing the chain to converge, but will take longer to achieve. Also controls the length of the pilot chain used to assess the required sampling length. Th
psrf.target
the value of the point estimate for the potential scale reduction factor of the Gelman Rubin statistic below which the chains are deemed to have converged (must be greater than 1). Default 1.05.
crash.retry
the number of times to re-attempt a simulation if the model returns an error. Default 1 retry (simulation will be aborted after the second crash).
thin.sample
option to thin the final MCMC chain(s) before calculating summary statistics and returning the chains. Thinning heavily autocorrelated, very long chains allows summary statistics to be calculated more quickly. If TRUE, the chain is thinned to as close t
jags
the system call or path for activating JAGS. Default calls findjags() to attempt to locate JAGS on your system.
silent.jags
should the JAGS output be suppressed? (logical) If TRUE, no indication of the progress of individual models is supplied. Default FALSE.
interactive
option to allow the simulation to be interactive, in which case the user is asked if the simulation should be extended when run length and convergence calculations are performed and the extended simulation will take more than 1 minute. The function will
max.time
the maximum time the function is allowed to run for, as a character string including units or as an integer in which case units are taken as seconds. Ignored if interactive==TRUE. If the function thinks that the next simulation extension will result in
adaptive
a list of advanced options controlling the length of the adaptive mode of each simulation. Extended simulations do not require an adaptive phase, but JAGS prints a warning if one is not performed. Reduce the length of the adpative phase for very time co

Value

  • A list including the following elements:
  • mcmcan MCMC object representing the model output for all chain(s). Renamed pilot.mcmc if the simulation is aborted before convergence or before the required sampling length is achieved
  • end.statethe end state of the last simulation extension performed. Can be used as initial values to extend the simulation further if required
  • req.samplesthe minimum sample size required for the Markov chain as calculated using the Raftery and Lewis's diagnostic
  • req.burninthe minimum burnin period for the chains to stabilise as calculated using the Raftery and Lewis's diagnostic
  • samples.to.convthe number of sampled iterations discarded due to poor convergence. The total number of iterations performed for the simulation is equal to req.samples + req.burnin + samples.to.conv (unless the simulation was aborted before reaching the required sampling length)
  • summarythe summary statistics for the monitored variables. Renamed pilot.summary if pilot.only==TRUE or if the simulation is aborted before the required sampling length is achieved
  • psrfthe Gelman Rubin statistic for the monitored variables (output of gelman.diag())
  • autocorrthe autocorrelation diagnostic for the monitored variables (output of autocorr.diag())

See Also

run.jags autorun.jagsfile raftery.diag gelman.diag autocorr.diag

Examples

Run this code
# run a model to calculate the intercept and slope of the expression y = m x + c, assuming normal observation errors for y:

# Simulate the data
x <- 1:100
y <- rnorm(length(x), 2*x + 10, 1)

# Model in the JAGS format
model <- "model {
for(i in 1 : N){
Y[i] ~ dnorm(true.y[i], precision);
true.y[i] <- (m * X[i]) + c;
}
m ~ dunif(-1000,1000);
c ~ dunif(-1000,1000);
precision ~ dexp(1);
}"

# Convert the data to a named list
data <- list(X=x, Y=y, N=length(x))

# Run the model
results <- run.jags(model=model, monitor=c("m", "c", "precision"), data=data)

# Analyse the results
results$summary

Run the code above in your browser using DataCamp Workspace