This final post about the STAN platform will specifically focus on data visualizations that can come from STAN models. In particular, we will explore these visualizations by hand with the popular shinystan package.

As we already know, the STAN platform typically uses particular Markov Chain Monte Carlo (MCMC) algorithms: the Hamiltonian Monte Carlo (HMC) or the No-U-Turn sampler (NUTS). It is also typically more efficient than other popular samplers, such as the Gibbs sampler or the Metropolis-Hasting algorithm. This is due to updates of the posterior converging more quickly to the stationary posterior distribution.

We have already introduced STAN in the three previous posts:
introduction with rstan, with rstanarm and with bayesplot. Now, the goal is to use all learned knowledge from before to build a shiny application through the shinystan framework.

Several vignettes are available for the shinystan package: For examples, see here or here for the deployment to shinapps server. The solutions to this set of exercises can be found here.

Exercise 1
In order to have a first visualization of the structure of a shiny application from STAN, run the demo app of the shinystan package. This demo app comes from the “Meta analysis” chapter of the STAN manual. You can also see it in this rstan vignette in the example section.

Exercise 2
We want to simulate data in order to build our own shiny application with the shinystan package.
We will take the same kind of examples in the previous set of exercises (see exercises with the “Bayesian inference” and “MCMC” tags) with the two-parameters Gumbel distribution coming from the Extreme Value Theory. This kind of posterior distribution is often hard to evaluate, hence why we need MCMC samplers.
Simulate a sample size of 100 from a Gumbel distribution (relying on the evd package) with location and scale parameters set to 10 and 3, respectively.
Do not forget to set the global options that allows you to automatically save a bare version of a compiled Stan program. Saving a bare version will ensure it does not need to be recompiled, along with executing multiple chains in parallel by taking all the available cores of your machine (use the parallel package.)

Exercise 3
Write the STAN model in your R script that estimates the generated sample by following the pertaining blocks of a STAN code:

– The data block that will gather the data used in the STAN model: n (the number of data points; they cannot be negative) and y (the observations of interest.)
– The parameters block: mu and sigma (cannot be negative) that represent the two parameters of a Gumbel distribution that we will have to estimate.
– The model block: defines the prior distribution for the two parameters: weakly informative centered on their sample estimate, for example, the empirical mean and the empirical SD of the sample with a variance of 25. Then, create the vector of interest from these two parameters.
– The generated quantities block: defines the posterior predictive values (y_pred) from the model.

Do not forget to constrain the upper and lower bounds in the data and parameter block declarations.

Note that it is recommended to write the STAN code in a separate file (“.stan” extension) and then call it directly through the stan() function instead of writing it in a character string, as done here. That would be easier to debug, read, etc.

Exercise 4
Define the named data list that will be used inside the STAN data block written above.

Exercise 5
Run the STAN model using the stan() function and use the following input parameters:
– The STAN code defined in Exercise 3
– The data list defined in Exercise 4
– 4 different chains
– 1000 iterations per chain
– A warm-up phase of 200. This phase defines the number of iterations that are used by the sampler for the adaptation phase before sampling begins; it is different than the burn-in phase we have seen in the Metropolis or Gibbs samplers.
Use all the available cores of your machine and make it reproducible (take the seed 1234.)

Exercise 6
Now, launch the shiny application that will allow you to have full access to the MCMC visualizations and diagnostics. Specify that you want to visualize the application in your own Rstudio viewer pane.

Exercise 7
Now, do the same as Exercise 6, but by viewing the application on your localhost directly in your default browser.

Exercise 8
Now, explore the shiny application and do the following:
a. Do a Trace-plot for the location parameter with all 4 of the chains
b. Do a Trace-plot for the location parameter, but visualize only the 3rd chain
c. Do a Trace-plot for the scale parameter with the 4th chain only. Moreover, take the log10 transformation of this parameter.
d. Find the standard deviation summary of the sampler parameters. Include the warm-up phase and display 3 decimals.
e. Check the effective sample sizes, the monte-carlo standard error and posterior SD, and the Gelman-Rubin Rhat. Set the warning levels to 50%, 15% and 1.01, respectively.
f. Compute the Partial Auto-correlation for every chain and for every parameter. (HINT: check the “Show/Hide” Options.) Then save the plot as a pdf.
After that, close the app.

Exercise 9
Now that you have diagnosed the convergence of the generated Markov chains, we can go through the estimation of the model. Go to the “estimate” tab in the main panel:
a. Compute the 95% and 70% posterior intervals for the parameter sigma. Include an indicator of the Rhat (Gelman Rubin diagnostic) in the plot.
b. Display the summary table of the parameter estimates. Highlight the scale parameter and order the table by effective sample size. Then, download the table.
c. Finally, go to the “explore” tab and check the bi-variate plot of the posterior draws of the two parameters.

Exercise 10 (BONUS)
Since this feature is still experimental, some bugs may still appear.
Check the fit of the model with the the posterior predictive distribution. Rely on the generated quantities of your STAN model (Exercise 3.)
You need to use the “PPcheck” tab and relaunch the app on your browser. It will check if you have both the predictive posterior samples and the sample data in a good format.