Lecture 6

class: center, middle, inverse, title-slide

# Lecture 6
## Introduction of Bayesian analysis using JAGS
### WILD6900 (Spring 2021)

---

## Readings

> Kéry & Schaub 38-40; 58-64

---
class: inverse, middle, center

# From custom MCMC to BUGS

---
## The BUGS language

#### **B**ayesian **A**nalysis **U**sing **G**ibbs **S**ampling

--
#### Language/program invented in the 1990's by epidemiologists in Cambridge

--
#### Later became WinBUGS

- First customizable MCMC software

--
- Revolutionized the use of Bayesian methods in applied statistics

---
## The BUGS language

#### Since the development of WinBUGS, other Bayesian software programs have been developed:

- OpenBugs

- [JAGS](https://sourceforge.net/projects/mcmc-jags/files/)

- [Nimble](https://r-nimble.org/)

- [Stan](https://mc-stan.org/)

--
#### For the remainder of the course, we will fit models using JAGS

- **J**ust **A**nother **G**ibbs **S**ampler

- Uses BUGS language (easy to learn, lots of online documentation)

- Often out-performs WinBUGS

- Available for all operating systems

---
## The BUGS language

#### Last week, we learned how to:

--
- Determine the *full conditional distributions* for a linear regression model

--
- Write a custom MCMC sampler to produce samples from the joint posterior distribution

--
### Given a statistical model and user-specified prior distributions, JAGS does these steps for you!

--
- Possible to fit arbitrarily complex models `$^1$`

--
- "Frees the modeler in you"

???

`$^1$` remember, just because you *can* fit it doesn't mean it will return robust inferences. As we have already discussed, complex models often suffer from lack of identifiability and poor sampling of the joint posterior distribution. But, if you can write out a valid statistical model, JAGS can create a sampler and draw samples from the joint posterior distribution

---
## Running JAGS from R

### JAGS is a stand alone software program

- Can be run from GUI or command line

--
### JAGS can also be run from `R` using the `jagsUI` package (among others)

--
- Write model in `R` script and save as `.jags` file

--
- Provide `jagsUI` with data, initial values, and MCMC settings

--
- model run in JAGS

--
- model output brought back in to `R` for diagnostics/analysis/visualization

---
## The BUGS language

### Very similar to `R` (but more limited)

- Limited ability to vectorize operations

--
### If you can write your model on the whiteboard, you can write it in JAGS

- Stochastic relationships represented by `~`

- Deterministic relationships represented by `<-`

---
## The BUGS language

### Linear regression model

`$$\LARGE y_i = \alpha + \beta \times x_i + \epsilon_i$$`

`$$\LARGE \epsilon_i \sim Normal(0, \tau)$$`

---
## The BUGS language

### Linear regression model

`$$\LARGE y_i = \underbrace{\alpha + \beta \times x_i}_{Deterministic} + \underbrace{\epsilon_i}_{Stochastic}$$`

???

Remember that `$\tau = \frac{1}{\sigma^2}$`

---
## The BUGS language

### Linear regression model

`$$\underbrace{\LARGE \mu_i = \alpha + \beta \times x_i}_{Deterministic}$$`

`$$\underbrace{\LARGE y_i \sim Normal(\mu_i, \tau)}_{Stochastic}$$`

Remember that these equations define the *likelihood* of our data given values of `$\alpha$`, `$\beta$`, and `$\tau$`

---
## The BUGS language

### Linear regression model

To specify a fully Bayesian model, we also need to define the priors:

`$$\LARGE [\alpha] \sim Normal(\alpha|0, 0.001)$$`

`$$\LARGE [\beta] \sim Normal(\alpha|0, 0.001)$$`

`$$\LARGE [\tau] \sim Gamma(\tau|0.01, 0.01)$$`

---
## The BUGS language

### Linear regression model

```
model{
 ## Priors
 alpha ~ dnorm(0, 0.001)
 beta ~ dnorm(0, 0.001)
 tau ~ dgamma(.001,.001) # Precision
 sigma <- 1/sqrt(tau) # Calculate sd from precision
 
 ## Likelihood
 for(i in 1:N){
 mu[i] <- alpha + beta * x[i]
 y[i] ~ dnorm(mu[i], tau)
 }
} #end of model
```

---
## Writing model files

```r
sink(file="jags/linear_regression.jags")
cat("
 model{
 ## Priors
 alpha ~ dnorm(0, 0.001)
 beta ~ dnorm(0, 0.001)
 tau ~ dgamma(.001,.001) # Precision
 sigma <- 1/sqrt(tau) # Calculate sd from precision
 ## Likelihood
 for(i in 1:N){
 mu[i] <- alpha + beta * x[i]
 y[i] ~ dnorm(mu[i], tau)
 }
 } #end of model
 ", fill=TRUE)
sink()
```

---
## Preparing the data

```r
## Read simulated data frame
dat <- readRDS("data/sim_seed_counts.rds")

## Store data for JAGS as list
jags_data <- list(y = dat$y, x = dat$visits.c, N = nrow(dat))

## Create function that returns random initial values
jags_inits <- function(){list(alpha = runif(1, 200, 300), 
 beta = runif(1, 25, 75), 
 tau = runif(1))}

## Tell JAGS which parameters we want to monitor
params <- c("alpha", "beta", "tau", "sigma")
```

---
## Run the model

```r
## Run the model
jags_fit <- jagsUI::jags(data = jags_data, inits = jags_inits, 
 parameters.to.save = params, 
 model.file = "jags/linear_regression.jags", 
 n.chains = 3, n.iter = 10000, 
 n.burnin = 2500, n.thin = 1)
```

---
## Diagnostics

```r
print(jags_fit)
```

```
## JAGS output for model '/Library/Frameworks/R.framework/Versions/4.0/Resources/library/WILD6900/jags/linear_regression.jags', generated by jagsUI.
## Estimates based on 3 chains of 10000 iterations,
## adaptation = 100 iterations (sufficient),
## burn-in = 2500 iterations and thin rate = 1,
## yielding 22500 total samples from the joint posterior. 
## MCMC ran for 0.032 minutes at time 2021-01-13 09:41:26.
## 
## mean sd 2.5% 50% 97.5% overlap0 f Rhat n.eff
## alpha 249.462 0.600 248.296 249.462 250.645 FALSE 1 1 9616
## beta 49.656 0.590 48.507 49.652 50.796 FALSE 1 1 22500
## tau 0.016 0.002 0.013 0.016 0.020 FALSE 1 1 22500
## sigma 7.884 0.428 7.098 7.866 8.783 FALSE 1 1 22500
## deviance 1218.809 2.454 1215.986 1218.177 1225.131 FALSE 1 1 22500
## 
## Successful convergence based on Rhat values (all < 1.1). 
## Rhat is the potential scale reduction factor (at convergence, Rhat=1). 
## For each parameter, n.eff is a crude measure of effective sample size. 
## 
## overlap0 checks if 0 falls in the parameter's 95% credible interval.
## f is the proportion of the posterior with the same sign as the mean;
## i.e., our confidence that the parameter is positive or negative.
## 
## DIC info: (pD = var(deviance)/2) 
## pD = 3 and DIC = 1222 
## DIC is an estimate of expected predictive error (lower is better).
```

---
## Diagnostics

```r
jagsUI::traceplot(jags_fit)
```