1. Markov chain diagnostics & reproducibility
No sample is perfect. Markov chains are no different. In this video, we'll explore some simple Markov chain diagnostics and discuss best practices for Markov chain reproducibility.
2. Markov chain output
Consider the first 1000 iterations of our `m` Markov chain. The trace plot (left) illustrates the *longitudinal* behavior of the chain. The density plot (right) illustrates the *distribution* of the values visited along this chain and *approximates* the posterior model of `m`.
3. Questions to consider
The fact that this is a mere approximation begs the following questions:
"What does a "good" Markov chain look like?";
"How accurate is the Markov chain approximation of the posterior?";
and "For how many iterations should we run the Markov chain?". RJAGS provides a set of diagnostic tools that provide insight into these questions.
4. Diagnostic: trace plots
First, re-consider our Markov chain trace plot.
5. Diagnostic: trace plots
Note the stability in the long-run behavior of this chain. The Markov chain trend captured by the blue model smooth appears roughly constant. This is good!
Contrast this with another chain in which the long-run trend has not yet stabilized after 1000 iterations. Perhaps it will continue traversing downward. Perhaps it will swing back up. Without knowing, we should increase the chain length or number of iterations until the chain stabilizes.
6. Diagnostic: multiple chains
Running multiple *parallel* chains provides another diagnostic. By default, when we compile a model using `jags.model()`, the subsequent simulation returns *one* Markov chain.
7. Diagnostic: multiple chains
With the inherent randomness in simulation, a second chain would take a different path.
8. Diagnostic: multiple chains
Compiling multiple, say 4 chains, illustrates a range of alternative Markov chain results. The 4 chains here are indeed different, but enjoy similar features and trends. Again, this similarity reflects stability in our Markov chain simulation.
9. Diagnostic: multiple chains
In contrast, *these* 4 chains take discrepant paths, thus would provide discrepant approximations of the posterior. In the face of such instability, it would be a mistake to stop our simulation here. Rather, we should increase the length of our chain.
10. Diagnostic: standard error
Finally, summary statistics of the simulation output provide additional and crucial Markov chain diagnostics.
Focusing on table 1 for now, notice that the mean `m` chain value was 29.10. This is an *estimate* of the posterior mean of `m`.
The corresponding *naive standard error*, 0.2836, measures the potential *error* in this estimate. This error is calculated by dividing the standard deviation of the `m` chain values by the square root of the number of iterations (or sample size). This is a *naive* approach to calculating error in a setting with dependent samples, but it suits our current purposes fine.
11. Diagnostic: standard error
Let's examine this concept in pictures. On a density plot of the `m` chain values, the red line represents the chain mean, our estimate of the posterior mean.
Zooming in, notice the error bars of +/- 2 standard errors around the mean. If you deem these error bars to be too wide, no problem! You can simply increase the number of Markov chain iterations. Smaller errors will follow.
12. Markov chain work flow
To summarize the Markov chain work flow:
after defining, compiling, and simulating a Bayesian model,
you should perform diagnostic checks of the simulation. This should include the examination of trace plots, multiple chain output, and standard errors.
Last but not least, you must *finalize* your simulation.
13. Finalizing the Markov chain: Reproducibility
Reproducibility and preservation of your work are crucial here. For reproducible simulation output, you must set the seed of the RJAGS random number generator. This works differently than in base `R`. In RJAGS, the seed is specified when you compile the model. To this end, there are many options - here we used the `inits` argument to set the random number generator to Wichmann-Hill and the seed to 1989.
14. Let's practice!
Time to practice your new diagnostic and reproducibility skills.