*Developed by Naomi Schalken and Rens van de Schoot*

**-----------------------------------------------------------------------**

### Contents

**Exercise 1 - WAMBS checklist**

Example data 2

**The WAMBS Checklis**

When to worry, and how to Avoid the Misuse of Bayesian Statistics

**To be checked before estimating the Model**

**To be checked after estimation but before inspecting Model**

Step 1: Do you understand the priors?

Step 2: run the model and check for convergence Step 3: doubling the number of iterations and convergence diagnostics

Step 4: Does the posterior distribution histogram have enough precision?

Step 5: Do the chains exhibit a strong degree of autocorrelation

Step 6: Does the posterior distribution make substantive sense?

Step 7: Do different specification of the multivariate variance priors influence the results?

Step 8: Is there a notable effect of the prior when compared with non-informative priors?

Step 9: Are the results stable from a sensitivity analysis??

Step 10: Is the Bayesian way of interpreting and reporting model results used?

**-----------------------------------------------------------------------**

### Exercise 1 - WAMBS checklist

In this exercise we will be following the first eight steps of the When-to-Worry-and-How-to-Avoid-the-Misuse-of-Bayesian-Statistics - checklist, or the WAMBS-checklist (Depaoli and van de Schoot 2015), to analyze data of PhD-students.

#### Example Data

The model we will be using for the first exercise is based on a study about predicting PhD-delays (Van de Schoot, Yerkes, Mouw and Sonneveld 2013). Among many other questions, the researchers asked the Ph.D. recipients how long it took them to finish their Ph.D. thesis (n=304). It appeared that Ph.D. recipients took an average of 59.8 months (five years and four months) to complete their Ph.D. trajectory. For the current exercise we will answer the question why some Ph.D. recipients took longer than other Ph.D. recipients by using age as a predictor (M = 30.7, SD = 4.48, min-max = 26-69). The relation between completion time and age is expected to be non-linear. It is expected that at a certain point in your life (i.e., mid-thirties), family life takes up more of your time than when you are in your twenties or when you are older. However, we expect that if you are in your mid-thirties and you are doing a Ph.D. you also take this extra time into account. The researchers asked Ph.D. candidates about their planned graduation day according to the original contract and their actual graduation day. The average gap between the two data appeared to be 9.6 months (SD = 14.4, min-max = -31 to 69). We expect that the lag between planned and actual time spent on the trajectory is less prone to non-linearity compared to actual project time.

So, in our model the GAP is the dependent variable and AGE and AGE2 are the predictors. The data can be found in the file `data_example.txt`

.

#### Step 1: Do you understand the priors?

Before actually looking at the data we first need to think about the prior distributions and hyperparameters for our model. For the current model, we need to think about 4 priors:

- the intercept
- the two regression parameters (β1 for the relation with AGE and β2 for the relation with AGE2)
- and the residual variance (∈)

Then, we need to determine which distribution to use for the priors. Let’s use for the

- intercept a normal prior with N(μ0, σ20), where μ0 is the prior mean of the distribution and σ20 is the variance parameter
- β1 a normal prior with N(μ1, σ21)
- β2 a normal prior with N(μ2, σ22)
- ∈ an Inverse Gamma distribution with IG(κ0,θ0), where κ0 is the shape parameter of the distribution and θ0 the rate parameter

Next, we need to specify actual values for the hyperparameters of the prior distributions. Let’s say we found a previous study and based on this study the following hyperparameters can be specified:

- intercept ~ N(10,5)
- β1 ~ N(.8,5)
- β2 ~ N(0,10)
- ∈ ~ IG(0.5, 0.5)

#### Mplus

Run the `Model_1.inp`

file and look at the prior distribution plots of all parameters. Resist looking at any other results for the moment.

Fill in the following **Table 1** that you find under “Input/Output Files”. Also, you can directly download a document with all tables by clicking on the tables below!

**Table 1**

#### Step 2: run the model and check for convergence

Open Mplus and then open the input file `Model_2.inp`

.

Variables in the data file:

- months = project time
- lag = lag between planned PhD duration and actual PhD duration
- age = age of participant (centered)
- age_sq = age squared

We specified the number of iterations by adding the following line to the ANALYSIS part of the input:

`FBITERATIONS = 200;`

In the input file, we have started writing down the code to set the model priors:

`MODEL PRIORS:`

a ~ N(X, X);

b ~ N(X, X);

c ~ N(X, X);

d ~ IG(X, X);

Look at the labels in the MODEL part of the syntax (a-d) and match the priors we decided on before to the Xs in the input file.

We also added the commands

`OUTPUT: TECH8;`

PLOTS: TYPE = PLOT2;

to the model syntax. The first command includes the iterations and a fit statistic for those iterations in the output file. The second command tells Mplus to create the plots we need to inspect convergence of the model visually.

Run the model and inspect the plots by clicking on Plots > View plots and choose one of the options (hint: you are interested in posterior trace plots). Fill out the first two columns of **Table 2** with copies of your plots (the final three will be filled out later).

**Table 2**

#### Step 3: doubling the number of iterations and convergence diagnostics

To obtain convergence diagnostics to check convergence for all parameters we use the Gelman and Rubin diagnostic (Gelman and Rubin 1992). This diagnostic should be close to 1.00 (at least < 1.05) to indicate convergence. You can find this diagnostic in the TECH8 output under Potential Scale Reduction (scroll down). **Did the model reach convergence?**

We want to be absolutely sure the model reached convergence and therefore we re-run the model with more than twice the number of iterations and compute the relative bias. First, re-run the model with 500 iterations:

`FBITERATIONS = 500;`

Save the new input file as `Model_3.inp`

and run the model. Fill out the table below using the formula to compute the relative bias:

**Table 3**

You should then use substantive knowledge about the metric of the parameter of interest, as well as substantive interpretations of the amount of fluctuation exhibited between chains, to determine when levels of relative deviation are negligible or problematic. For example, with a regression coefficient of 0.001, a 10% relative deviation level might not be substantively relevant. However, with an intercept parameter of 50, a 10% relative deviation level might be quite meaningful. The specific level of relative deviation should be interpreted in the substantive context of the model. Some examples of interpretations are:

- if relative deviation is < |1|%, then do not worry;
- if relative deviation > |1|%, then rerun with 4x nr of iterations.

Did the model reach convergence, or do you want to re-run the model with even more iterations? Continue until you are satisfied.

#### Step 4: Does the posterior distribution histogram have enough precision?

To check the histograms of Model_3 go to Plots > View Plots and find the posterior histograms for all parameters of interest. Copy-paste the histograms to Table 2.**What is your conclusion; do we need even more iterations?**

#### Step 5: Do the chains exhibit a strong degree of autocorrelation

To check the autocorrelation plots of Model_3 go to Plots > View Plots and find the autocorrelation plots for both chains. Copy-paste the plots to Table 2.**What is your conclusion; do we need even more iterations?**

#### Step 6: Does the posterior distribution make substantive sense?

To check the Kernel density plots of Model_3 go to Plots > View Plots and find the posterior distribution kernel density plots for all parameters of interest. Copy-paste the plots to Table 2. **What is your conclusion; do the posterior distributions make sense?**

#### Step 7: Do different specification of the multivariate variance priors influence the results?

To understand how the prior on the residual variance impacts the results, we compare the previous model with a model where different hyperparameters for the Inverse Gamma distribution are specified.

The specification of the hyperparameters .5, .5 was inspired by this paper (Van de Schoot, Broere, Perryck, Zondervan-Zwijnenburg and Van Loey 2015) but appears to be a very informative prior. A non-informative specification would be to use .001,.001 (Asparouhov and Muthén 2010).

Save a new input file `Model_4.inp`

and change the inverse gamma prior to the numbers we specified above. Run the model in Mplus.

Compute the bias between this model and `Model_3.inp`

.

**Table 4**

**Are there any substantial differences in the results?**

#### Step 8: Is there a notable effect of the prior when compared with non-informative priors?

To see what the influence of the other priors is on the model estimates, remove the prior specification in the input file for the normal distributions (intercept, beta1, beta2). This way, Mplus will use its default prior N(0,1010) again. Save the input file as `Model_5.inp`

. Run the model and, for now, only look at the prior distribution plots. Mplus provides an error message about the availability of the prior plots. What does this error message indicate?

Now look at the rest of the model output and posterior plots. Did this non-informative model reach convergence (think about all the previous steps)?

If so, compare this model to the model with informative priors `Model_3.inp`

and calculate the bias.

**Table 5**

**What is your conclusion about the influence of the priors on the posterior results?****Which results will you use to draw conclusion on?**

#### Step 9: Are the results stable from a sensitivity analysis??

From the original paper:

“If informative or weakly-informative priors are used, then we suggest running a sensitivity analysis of these priors. When subjective priors are in place, then there might be a discrepancy between results using different subjective prior settings. A sensitivity analysis for priors would entail adjusting the entire prior distribution (i.e., using a completely different prior distribution than before) or adjusting hyperparameters upward and downward and re-estimating the model with these varied priors. Several different hyperparameter specifications can be made in a sensitivity analysis, and results obtained will point toward the impact of small fluctuations in hyperparameter values. [.] The purpose of this sensitivity analysis is to assess how much of an impact the location of the mean hyperparameter for the prior has on the posterior. [.] Upon receiving results from the sensitivity analysis, assess the impact that fluctuations in the hyperparameter values have on the substantive conclusions. Results may be stable across the sensitivity analysis, or they may be highly instable based on substantive conclusions. Whatever the finding, this information is important to report in the results and discussion sections of a paper. We should also reiterate here that original priors should not be modified, despite the results obtained.”

If you still have time left, you can adjust the hyperparameters upward and downward and re-estimating the model with these varied priors to check for robustness.

#### Step 10: Is the Bayesian way of interpreting and reporting model results used?

After Points 1-9 have been carefully considered and addressed, model interpretation and the reporting of results become the next concern. For inspiration how to do this, check Appendix A of the WAMBS-paper or read the paper “Bayesian analyses: where to start and what to report” available in the dropbox folder under ‘literature’.

**References**

*Asparouhov, T., & Muthén, B. (2010). Bayesian analysis of latent variable models using Mplus. Retrieved June 17, 2014, from http://www.statmodel2.com/download/BayesAdvantages18.pdf.*

Depaoli, S., and van de Schoot, R. (2015). Improving transparency and replication in Bayesian Statistics: The WAMBS-Checklist. *Psychological Methods*.

*Gelman, A and Rubin, DB (1992) Inference from iterative simulation using multiple sequences,** **Statistical Science**,** **7**, 457-511*

*Van de Schoot, R., Broere, J.J., Perryck, K., Zondervan-Zwijnenburg, M., & Van Loey, N. (2015). **Analyzing Small Data Sets using Bayesian Estimation: The case of posttraumatic stress symptoms following mechanical ventilation in burn survivors. European Journal of Psychotraumatology.*

Van de Schoot, R., Yerkes, M., Mouw, J. & Sonneveld, H. (2013). What took them so long? Explaining PhD delays among Doctoral Candidates. PLoS One, 8(7), e68839. doi: 10.1371/journal.pone.0068839.

### Other tutorials you might be interested in

**The WAMBS-Checklist**

RJAGS: How to avoid and when to worry about the misues of Bayesian Statistics

Blavaan: How to avoid and when to worry about the misuse of Bayesian Statistics

-----------------------------------------------------------

**First Bayesian Inference**