<!--
%\VignetteEngine{knitr::knitr}
%\VignetteIndexEntry{A quick-start guide to running models in JAGS}
-->

```{r, echo = FALSE, include=FALSE}
library(runjags)
if(require(rjags)){
	load.module('dic')
}
library(parallel)
```

# A quick-start guide to using the runjags package

This document is intended as a short introduction to fitting MCMC models using JAGS and runjags, focussing on learning by example.  The intention is to keep this document as brief as possible, but a slightly longer version with more examples is available from [here](https://runjags.sourceforge.net/quickjags.html).  Additional information and more details are available in the [user guide](https://runjags.sourceforge.net/userguide.pdf).  For detailed documentation of individual functions, see the runjags manual for that function.

## Preparation

The runjags packages requires a working installation of [JAGS](https://mcmc-jags.sourceforge.io/).  To test the installation, it is recommended to use the testjags() function:

```{r}
testjags()
```

This should verify that your JAGS installation has been found.  The rjags and modeest packages are not strictly required, but are reccomended packages and should be installed where possible.

```{r, echo = FALSE, include=FALSE}
runjags.options(inits.warning=FALSE, rng.warning=FALSE, blockignore.warning=FALSE, blockcombine.warning=FALSE, silent.jags=TRUE, silent.runjags=TRUE)
```

## Running models using run.jags

A recommended introductory text for individuals not familiar with the BUGS language is:

Lunn D, Jackson C, Best N, Thomas A, Spiegelhalter D (2012). The BUGS book: A practical introduction to Bayesian analysis. CRC press.

To illustrate how to use run.jags we will look at the Salmonella example from chapter 6.5.2 of this book, with thanks to Lunn, Jackson, Best, Thomas, and Spiegelhalter 2012, for kind permission to reproduce this code.  The model below was taken from the WinBUGS example, modified to include a mandatory pair of curly brackets to demarcate the Data and Inits blocks, and adding an optional second Inits block with different starting values for 2 chains:

```{r}
filestring <- "
Poisson model...

model {
for (i in 1:6) {
for (j in 1:3) {
y[i,j] ~ dpois(mu[i])
}
log(mu[i]) <- alpha + beta*log(x[i] + 10) + gamma*x[i]
}
for (i in 1:6) {
y.pred[i] ~ dpois(mu[i])
}
alpha ~ dnorm(0, 0.0001)
beta ~ dnorm(0, 0.0001)
gamma ~ dnorm(0, 0.0001)
}

Data{
list(y = structure(.Data = c(15,21,29,16,18,21,16,26,33,
27,41,60,33,38,41,20,27,42),
.Dim = c(6, 3)),
x = c(0, 10, 33, 100, 333, 1000))
}

Inits{
list(alpha = 0, beta = 1, gamma = -0.1)
}

Inits{
list(alpha = 10, beta = 0, gamma = 0.1)
}
"
```

The following call to run.jags reads, compiles, runs and returns the model information along with MCMC samples and summary statistics:

```{r, results='hide'}
results <- run.jags(filestring, monitor=c('alpha','beta','gamma'))
```

Note that this model is run from an embedded character vector here for convinience, but run.jags will also accept a file path as the main argument in which case the information will be read from the specified file.  There is a print method associated with runjags class objects:

```{r}
results
```

This method displays relevant summary statistics - the effective sample size (SSeff) and convergence diagnostic (psrf) is shown for each stochastic variable.  It is also advisable to plot the results using the default plot method to assess convergence graphically:

```{r, results='hide', include=FALSE}
plot(results, layout=c(3,4), file='plots.pdf', height=10, width=8)
```

```{r, results='hide'}
plot(results, layout=c(3,4))
```

![](plots.pdf)

In this case the chains have converged, but there is a large amount of auto-correlation and therefore a small effective sample size, which can be reduced by loading the GLM module in JAGS:

```{r, results='hide'}
resultswithglm <- run.jags(filestring, monitor=c('alpha','beta','gamma'), modules='glm')
```
```{r}
resultswithglm
```
This second model has a smaller AC.10 (autocorrelation with lag 10) and therefore larger SSeff (effective sample size), so there will be much less Monte Carlo error associated with these estimates comapred to the first model.

## Using data and initial values from R

To facilitate running models using data that is already in R, runjags will look for a number of comment tags in the model file to indicate the variables that should be extracted from the R working environment.  For example the following model uses N, Y and X as data, and coef, int and precision are given initial values, all of which are taken from R.  In addition, the model indicates to monitor coef, int and precision so that no monitor vector is required in the call to run.jags:

```{r, results='hide'}
# Simulate the data:
set.seed(1)
X <- 1:100
Y <- rnorm(length(X), 2*X + 10, 1)
N <- length(X)

model <- "
	model {
        for(i in 1 : N){ #data# N
        Y[i] ~ dnorm(true.y[i], precision) #data# Y
        true.y[i] <- (coef * X[i]) + int #data# X
	}
  	coef ~ dunif(-1000,1000)
  	int ~ dunif(-1000,1000)
  	precision ~ dexp(1)
  	#inits# coef, int, precision
  	#monitor# coef, int, precision
}"
# A function to return initial values for each chain:
coef <- function(chain) return(switch(chain, "1"= -10, "2"= 10))
int <- function(chain) return(switch(chain, "1"= -10, "2"= 10))
precision <- function(chain) return(switch(chain, "1"= 0.01, "2"= 100))

# Run the simulation:
results <- run.jags(model, n.chains = 2)
```

The output of this function can be extended using:

```{r, results='hide'}
results <- extend.jags(results, sample=5000)
```

This function call takes an additional 5000 samples and combines them with the original simulation.  It is also possible to automatically extend a model until it has converged (as assessed by Gelman and Rubin's convergence diagnostic) using the autorun.jags and autoextend.jags functions, although manual inspection of trace plots is always recommended to ensure chains have really converged.  Note also that each of the functions run.jags, extend.jags, autorun.jags and autoextend.jags allows a method argument to be specified, which allows separate chains to be run on parallel processors where available.  Full details are provided are in the user guide, but the most used options are 'rjags', 'interrputible', 'parallel' and 'background'.  

## Generating a GLMM template

One of the most difficult things about getting started with JAGS or BUGS is the requirement to write the model yourself.  Looking at examples is a good way to learn, but it is even easier when the example is tailored towards your specific model.  The template.jags function creates a model template file based on an lme4-style formula interface.  We can look at an example from the help file for lm:

```{r, results='hide'}
ctl <- c(4.17,5.58,5.18,6.11,4.50,4.61,5.17,4.53,5.33,5.14)
trt <- c(4.81,4.17,4.41,3.59,5.87,3.83,6.03,4.89,4.32,4.69)
group <- gl(2, 10, 20, labels = c("Ctl","Trt"))
weight <- c(ctl, trt)
D9 <- data.frame(weight, group)

# The JAGS equivalent:
model <- template.jags(weight ~ group, D9, n.chains=2, family='gaussian')
```

The model is written to disk by template.jags, and includes all of the data and initial values necessary to run the model.  Examining this model file is a good way to learn how the BUGS syntax works, as well as some of the options allowed by the runjags package.  To run the models:

```{r, results='hide'}
JAGS.D9 <- run.jags(model)
lm.D9 <- lm(weight ~ group, data=D9)
```
And to compare results:

```{r}
JAGS.D9
summary(lm.D9)
```

Note that lm reports sigma and JAGS the precision - to make them more comparable we could use a mutate function:

```{r, results='hide'}
JAGS.D9 <- run.jags(model, mutate=list(prec2sd, 'precision'))
```

And focus on the results from the summary that we are interested in:

```{r}
summary(JAGS.D9, vars=c('regression_precision.sd', 'intercept', 'group_effect'))
```

We can also compare the estimated residuals using for example:

```{r}
summary(residuals(lm.D9) - residuals(JAGS.D9, output='mean'))
```


## Model evaluation

The Deviance Information Criterion (DIC) is typically used to compare models fit using JAGS and BUGS, and is displayed as part of the print method shwon above.  This can also be extracted from a runjags model using the extract method, along with other information such as the samplers used:

```{r}
extract(JAGS.D9, what='samplers')
```

Another possible method to assess the fit of Bayesian models is using leave-one-out cross-validation, or a drop-k study.  A utility function is included within runjags to facilitate this - for more details see the help file for the drop.k function or the [extended version of this vignette](https://runjags.sourceforge.net/quickjags.html).


## Options

There are a large number of global options that can be controlled using the runjags.options function.  For an explanation of these see the help file for the relevant function:

```{r, eval=FALSE}
?runjags.options
```


## Further information

There is a [sourceforge repository](https://runjags.sourceforge.net/) associated with the runjags package, which contains additional versions of the package and internal JAGS module.  This page also contains the [extended version of this vignette](https://runjags.sourceforge.net/quickjags.html) which contains additional material and examples not found in the [version hosted on CRAN](https://CRAN.R-project.org/package=runjags).  

This package is in constant development, so if you find any issues or the code does not behave as expected please either email the pacckage developer or visit the [bug-tracker page](https://github.com/ku-awdc/runjags/issues) to report the problem.  General comments and suggestions on how you find using this package are also very welcome.


## Citation

Devleopment of this R package has consumed a great deal of time and energy, and is a somewhat peripheral activity in my current academic position.  In order to justify continued devleopment of the software I need to be able to demonstrate that it is being used within the general research community.  So, if you use this software as part of a publication, please remember to cite the software as follows:

```{r}
citation('runjags')
```

It is also important to cite JAGS and R!  



## This vignette was built with:

```{r}
sessionInfo()
```

