### ----------------------------------------------------

### Contents

**Exercise 1 - WAMBS checklist**

Example data

**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 = -3 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.dat`

.

#### 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, τ0), where μ0 is the prior mean of the distribution and τ0 is the precision parameter
- β1 a normal prior with N(μ0, τ0)
- β2 a normal prior with N(μ0, τ0)
- e an Inverse Gamma distribution with IG(κ0,θ0), where θ0 is the shape parameter of the distribution and κ0 the rate parameter

Note: In Bugs/Jags they use tau, τ (= precision), instead of sigma^2, σ2 (= variance) or sigma σ (= SD). If you have information about the variance (σ2) then you take the inverse by taking the quotient: 1/σ2. For example, if you want a prior variance of 100 in Jags, you do 1/variance = 1/100 = 0.01.

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)

To make some nice plots of the prior distribution we created a function in R. First copy paste the following syntax into R (without making any changes):

`prior.plot <- function(type="normal", mean=0, variance=100, shape=0, scale=1, sec.min=-50, sec.max=50, step=.1, label=label) {`

x <- seq(sec.min, sec.max, by = step)if (type == "normal") {

prior.d <- dnorm(x,mean = mean, sd = sqrt(variance))

}if (type == "gamma") {

prior.d <- dgamma(x,shape = shape, scale = scale)

}

`plot(x, prior.d, type="l", xlab = label)`

abline(v = mean, col="red")

abline(v = (mean+2*sqrt(variance)), col = "blue")

abline(v = (mean-2*sqrt(variance)), col = "blue")

}

Input:

- Type = either “normal” or “gamma”
- Mean = the mean of the distribution
- Variance = the variance of the distribution
- Shape = for gamma: the shape of the distribution
- Scale = for gamma: the scale of the distribution
- sec.min/.max and step = make the interval shown on the x-axis smaller or wider and take smaller or larger steps
- label = label of the x-axis

For example, to obtain a plot for the intercept with a normal prior with a mean of 10 and a variance of 100, use the syntax:

`prior.plot(type = "normal", mean = 10, variance = 100, label="intercept")`

To plot all four priors use this code, but fill in values for X:

`par(mfrow=c(2,2))`

prior.plot(type = "normal", mean = X, variance = X, label="intercept")

prior.plot(type = "normal", mean = X, variance = X, label="beta 1")

prior.plot(type = "normal", mean = X, variance = X, label="beta 2")

prior.plot(type = "gamma", shape = X, scale = X, sec.min=0, sec.max=0.5, step=.001, label="res.var")

par(mfrow=c(1,1))

Now you can fill in the first table of the WAMBS-checklist (see the next page):

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

Load the data and make sure you set the working directory to the folder where your data is stored.

To set your working directory to the right folder use:

`setwd("C:/your/folder/here")`

For example:

`setwd("E:/exercise-day-4")`

To load the data

`data <- read.table("data_example.dat")`

Variables in the model:

- y1 = project time
- y2 = lag between planned PhD duration and actual PhD duration
- x1 = age of participant (centered)
- x2 = age squared

To make the code we will use later a bit easier, we give names to the variables in the dataset:

`colnames(data)<- c("y1", "y2", "x1", "x2")`

Specify the regression model you want to analyze and run the regression with the blavaan() function. The priors for the intercept, the regression coefficients and the variance are specified in the regression model below. Do not forget to change the labels in red (make sure to use the precision for the prior variance - to go from variance to precision do 1/variance = precision).

`model.regression <- '#regression`

y2 ~ prior("dnorm(X, X)")*x1 + prior("dnorm(X, X)")*x2#variance

y2 ~~ prior("dgamma(X, X)")*y2#intercept

y2 ~ prior("dnorm(X, X)")*1'

`fit.bayes <- blavaan(model.regression, data=data, n.chains = 2, burnin = 100, sample = 100, adapt=100)`

In the blavaan() function, you first specify the model object and then the dataset with the variables of interest. Then the number of chains is given. After that you specify the number of burn-in iterations, the size of the sample to draw from the posterior distribution and the number of adaptive iterations.

To obtain summary estimates of the parameters use

`summary(fit.bayes, fit.measures=TRUE, ci = TRUE)`

To obtain diagnostic plots to check convergence for all parameters use

`plot(fit.bayes, pars= 1:4, plot.type="trace")`

Now, you can fill in columns 1 and 2 of Table 2 (the information for the other columns will be added later in the exercise):

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

To obtain the G-R diagnostic and plots use

`mcmc.list <- blavInspect(fit.bayes, what="mcmc")`

gelman.diag(mcmc.list)

gelman.plot(mcmc.list)

Look at the Upper CI/Upper limit this should be close to 1. Note: The G-R diagnostic is also automatically given in the summary of blavaan of Step 2 under the column PSRF. To obtain the Geweke diagnostic use

`geweke.diag(mcmc.list)`

This syntax will provide a z-score, where a z-score higher than the absolute value of 1.96 is associated with a p-value of < .05 (two-tailed). These values should therefore be lower than the absolute value of 1.96.

If Geweke indicates that the first and last part of a sample from a Markov chain are not drawn from the same distribution, it may be useful to discard the first few iterations to see if the rest of the chain has “converged”. This plot shows what happens to Geweke’s z-score when successively larger numbers of iterations are discarded from the beginning of the chain after burnin.

`par(mfrow=c(4,2))`

geweke.plot(mcmc.list, auto.layout=FALSE)

**Did the model reach convergence?**

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

… and run the model with the following code:

`fit.bayes2 <- blavaan(model.regression, data=data, n.chains = 2, burnin = X, sample = X, adapt=X)`

summary(fit.bayes2, fit.measures=TRUE, ci = TRUE)

Although you could compute the bias by hand, I created an easy function to quickly compute bias of parameter estimates between two models. First copy past the following syntax into R (without making any changes):

`bias.blavaan <- function(model1, model2) {`

bias <- as.vector(1)

for (i in 1:nrow(model1)) {

bias[i] <- (model2[i,1]-model1[i,1])/model1[i,1]*100

}

bias <- cbind(bias, as.vector(model1[,1]), as.vector(model2[,1]))

rownames(bias) <- c("beta1", "beta2","variance", "intercept")

colnames(bias) <- c("% bias", "est model 1", "est model 2")

print(bias)

}

To compare the two models first retrieve the parameter estimates from the models you want to compare and then put these in a matrix for both models.

`summary1 <- blavInspect(fit.bayes, what="postmean")`

summary2 <- blavInspect(fit.bayes2, what="postmean")

model1 <- as.matrix(summary1)

model2 <- as.matrix(summary2)

Then use:

`bias.blavaan(model1, model2)`

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

**What can be concluded about the bias between the two models?**

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

Histograms can be obtained in two different ways. The easiest one is as follows:

`par(mfrow(2,2))`

plot(fit.bayes2, pars=1:4, plot.type="histogram")

Note: The par(mfrow(2,2)) code is added to make sure that the four histograms fit into the screen. Run these two lines all at once.

Another way to obtain the histograms, is with the following code:

`mcmc.list.2 <- blavInspect(fit.bayes2, what="mcmc")`

model.mcmc.2 <- as.matrix(mcmc.list.2)

par(mfrow=c(2,2))

hist(model.mcmc.2[,1], breaks = 50, main = "Intercept")

hist(model.mcmc.2[,2], breaks = 50, main = "Age")

hist(model.mcmc.2[,3], breaks = 50, main = "Age squared")

hist(model.mcmc.2[,4], breaks = 50, main = "Residual variance")

Copy-paste the histograms (for one of the two options) 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 there are also two possible codes (make sure to refer to fit.bayes2). The first one with the blavaan plot:

`par(mfrow(2,2))`

plot(fit.bayes2, pars=1:4, plot.type="autocorr", col='blue')

Or the second one with the autocorr.plot() function (based on the mcmc.list which is created in step 4 with the second option):

`par(mfrow(2,2))`

autocorr.plot(mcmc.list.2)

Copy-paste the plots (one of the two options) 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 use one of the following codes:

`par(mfrow(2,2))`

plot(fit.bayes2, pars=1:4, plot.type="density")

Or (make sure to refer to mcmc.list.2):

`par(mfrow(2,2))`

plot(mcmc.list.2, trace=FALSE, density=TRUE)

Copy-paste the plots (one of the two options) to Table 2. **What is your conclusion; do the posterior distributions make sense?**

#### Step 7: Do different specifications 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 Gamma distribution are specified.

The specification of the hyperparameters .5, .5 was inspired by this paper but appears to be a very informative prior. A non-informative specification would be to use .001, .001

Specify number of iterations and number of burn-in iterations:

`fit.bayes.3 <- blavaan(model.regression.2, data=data, n.chains = 2, burnin = X, sample = X, adapt=X)`

Change the hyperparameters in the regression model that was specified in Step 2:

`model.regression.2 <- '#regression`

y2 ~ prior("dnorm(0.8, 0.2)")*x1 + prior("dnorm(0, 0.1)")*x2

#variances

y2 ~~ prior("dgamma(X , X)")*y2

#intercepts

y2 ~ prior("dnorm(10, 0.2)")* 1'

Save the output and obtain the summary statistics:

`fit.bayes.3 <- blavaan(model.regression.2, data=data, n.chains = 2, burnin = 250, sample = 250, adapt=250)`

summary(fit.bayes.3, fit.measures=TRUE, ci = TRUE)

Save the parameter values of model 3 in the right format:

`summary3 <- blavInspect(fit.bayes.3, what="postmean")`

model3 <- as.matrix(summary3)

Calculate the bias between model 2 and model 3:

`bias.blavaan(model2, model3)`

**Are there any substantial differences?**

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

Specify number of iterations and number of burn-in iterations:

`fit.bayes.4 <- blavaan(model.regression.3, data=data, n.chains = 2, burnin = X, sample = X, adapt=X)`

Change the priors for the intercept and regression coefficients in the regression model that was specified in the previous step:

`model.regression.3 <- '#regression`

y2 ~ prior("dnorm(X , X)")*x1 + prior("dnorm(X, X)")*x2#variances

y2 ~~ prior("dgamma(.001, .001)")*y2

`#intercepts`

y2 ~ prior("dnorm(X, X)")* 1'

Did this non-informative model reach convergence (think about all the previous steps)? Save the output and obtain the summary statistics:

`fit.bayes.4 <- blavaan(model.regression.3, data=data, n.chains = 2, burnin = 250, sample = 250, adapt=250)`

summary(fit.bayes.4, fit.measures=TRUE, ci = TRUE)

Save the parameter values of model 4 in the right format:

`summary4 <- blavInspect(fit.bayes.4, what="postmean")`

model4 <- as.matrix(summary4)

Calculate the bias between model 2 and model 3:

`bias.blavaan(model2, model4)`

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

**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*

*Geweke, J. Evaluating the accuracy of sampling-based approaches to calculating posterior moments. In *Bayesian Statistics 4* **(ed JM Bernado, JO Berger, AP Dawid and AFM Smith). Clarendon Press, Oxford, UK.*

*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**

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

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

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

**First Bayesian Inference**