Question: Write down the null and alternative hypotheses that represent this question. Which hypothesis do you deem more likely?

Question: Have all your data been loaded in correctly? That is, do all data points substantively make sense? If you are unsure, go back to the .csv-file to inspect the raw data.

MODEL RESULTS

                                                    Two-Tailed
                    Estimate       S.E.  Est./S.E.    P-Value

 Means
    DIFF               9.967      0.790     12.622      0.000
    AGE               31.676      0.375     84.433      0.000
    AGE2            1050.217     35.916     29.241      0.000

Step 1: Do you understand the priors?

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, there are four priors:

• the intercept
• the two regression parameters (\(\beta_1\) for the relation with AGE and \(\beta_2\) for the relation with AGE2)
• the residual variance (\(\in\))

We first need to determine which distribution to use for the priors. Let’s use for the

• intercept a normal prior with \(\mathcal{N}(\mu_0, \sigma^{2}_{0})\), where \(\mu_0\) is the prior mean of the distribution and \(\sigma^{2}_{0}\) is the variance parameter
• \(\beta_1\) a normal prior with \(\mathcal{N}(\mu_1, \sigma^{2}_{1})\)
• \(\beta_2\) a normal prior with \(\mathcal{N}(\mu_2, \sigma^{2}_{2})\)
• \(\in\) an Inverse Gamma distribution with \(IG(\kappa_0,\theta_0)\), where \(\kappa_0\) is the shape parameter of the distribution and \(\theta_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 \(\sim \mathcal{N}(-35, 20)\)
• \(\beta_1 \sim \mathcal{N}(.8, 5)\)
• \(\beta_2 \sim \mathcal{N}(0, 10)\)
• \(\in \sim IG(.5, .5)\) This is an uninformative prior for the residual variance, which has been found to perform well in simulation studies.

It is a good idea to plot these distribution to see how they loo and waht expected delay would be given these priors. With these priors the regression formula would be: \(delay=-35+ .8*age + 0*age^2\). In the Blavaan version of the WAMBS checklist we explain how you can easily do that in R.

To run a multiple regression with Mplus, you first specify the model (after the MODEL line using the ON command for regression coefficients, and the [] command for the intercept/means). As default Mplus does not run a Bayesian analysis, so you would have to change the ESTIMATOR to BAYES under ANALYSIS in the input file and then look at the output under MODEL RESULTS. In Mplus, the priors are set in the under the MODEL PRIORS command using the ~ sign, an N for a normal distribution, and an IG for the inverse gamma distribution. To see what default priors Mplus uses we can have a look at the Appendix A of this manual.

The priors are presented in code as follows:

TITLE: Bayesian analysis;

DATA: FILE IS phd-delays_nonames.csv;

VARIABLE: NAMES ARE diff child sex Age Age2; ! All the variables in the dataset

USEVARIABLES ARE diff Age Age2; ! The variables we use in this analysis

ANALYSIS:
ESTIMATOR IS bayes; ! Specify that we want to use a Bayesian analysis
Bseed = 23082018; !specify a seed value for reproducing the results
CHAINS = 3; ! set the number of chains we want to use

MODEL: 
[diff] (intercept);       ! specify that we want an intercept
! this model would also work without this line, but this way it is possible to easily set a prior
diff ON Age (Beta_Age);   ! Regression coefficient 1. 
diff ON Age2(Beta_Age2);  ! Regression coefficient 2 
! You need to name these regression coefficients to later set priors
diff (e);                 ! Error variance

MODEL PRIORS:
  Beta_Age ~ N(.8, 5); ! These are the different priors
  Beta_Age2 ~ N(0, 10);
  intercept ~  N(-35, 20);
  e ~ IG(0.5, 0.5);


OUTPUT: 
  cinterval(hpd); !to request CIs based on higher posterior density intervals
  tech8; !to request additional information about priors and convergence
  
PLOT: 
    TYPE IS PLOT2; !to request for plots

To check the prior distribution in Mplus you have to first run the analysis adding PLOT: TYPE IS PLOT2; to the input file. Then, in the top menu of Mplus, go to Plot > View Plots > Bayesian prior parameter distribution > view > select the parameter you would like to see the prior of > kernel density > OK Inspect all priors. (Be careful the actual density you see, is the density of the posterior, the prior is summarized in the different coloured lines. This seems to be a bit of a bug in Mplus)

Step 2: Run the model and check for convergence

To run a multiple regression in Mplus with a specific number of iteratation, you first specify the model as before but add FBITERATIONS = xxx; to the ANALYSIS part of the input file.

2. Does the trace-plot exhibit convergence?

By default Mplus discards half of the iteration as burnin period. After running this model we can check the trace plots by going to Plot > View Plots > Bayesian posterior parameter trace plots > view > select the parameter you would like to see the prior of > OK

First, run the model with only 100 iterations. Inspect the traceplots, what do you conclude?

TITLE: Bayesian analysis;

DATA: FILE IS phd-delays_nonames.csv;

VARIABLE: NAMES ARE diff child sex Age Age2; ! All the variables in the dataset

USEVARIABLES ARE diff Age Age2; ! The variables we use in this analysis

ANALYSIS:
ESTIMATOR IS bayes; ! Specify that we want to use a Bayesian analysis
Bseed = 23082018; !specify a seed value for reproducing the results
CHAINS = 3; ! set the number of chains we want to use
FBITERATIONS = 100 ! this is how we set the number of iterators

MODEL: 
[diff] (intercept);       ! specify that we want an intercept
! this model would also work without this line, but this way it is possible to easily set a prior
diff ON Age (Beta_Age);   ! Regression coefficient 1. 
diff ON Age2(Beta_Age2);  ! Regression coefficient 2 
! You need to name these regression coefficients to later set priors
diff (e);                 ! Error variance

MODEL PRIORS:
  Beta_Age ~ N(.8, 5); ! These are the different priors
  Beta_Age2 ~ N(0, 10);
  intercept ~  N(-35, 20);
  e ~ IG(0.5, 0.5);


OUTPUT: 
  cinterval(hpd); !to request CIs based on higher posterior density intervals
  tech8; !to request additional information about priors and convergence
  
PLOT: 
    TYPE IS PLOT2; !to request for plots

Let’s specify a fixed number of iterations by adding FBITERATIONS = 2000; to the ANALYSIS part. Inspect the traceplots again.

We can double check check if the chains converged by having a look at the convergence diagnostics, the Gelman and Rubin diagnostic. To have a look at the Gelman-Rubin Diagnostic (PSRF) we can check this table under the TECH8outout of Mplus. The Gelman-Rubin Diagnostic shows the PSRF values (using the within and between chain variability), which should be close to 1 for ALL iterations after burn-in. Is this the case?.

   TECHNICAL 8 OUTPUT FOR BAYES ESTIMATION

     CHAIN    BSEED
     1        23082018
     2        147894
     3        830180

                     POTENTIAL       PARAMETER WITH
     ITERATION    SCALE REDUCTION      HIGHEST PSR
     100              1.000               1
     200              1.000               1
     300              1.000               1
     400              1.000               1
     500              1.001               2
     600              1.002               2
     700              1.002               2
     800              1.002               2
     900              1.002               2
     1000             1.001               2
     1100             1.000               2
     1200             1.001               1
     1300             1.001               1
     1400             1.001               1
     1500             1.001               1
     1600             1.000               2
     1700             1.000               1
     1800             1.000               1
     1900             1.000               1
     2000             1.000               1

3. Does convergence remain after doubling the number of iterations?

As is recommended in the WAMBS checklist, we double the amount of iterations to check for local convergence. To this by changing this part of the input.

FBITERATIONS = 4000 ! this is how we set the number of iterators

You should again have a look at the above-mentioned convergence statistics, but we can also compute the relative bias to inspect if doubling the number of iterations influences the posterior parameter estimates (\(bias= 100*\frac{(model \; with \; double \; iteration \; – \; initial \; converged \; model )}{initial \; converged \; model}\)). In order to preserve clarity we just calculate the bias of the two regression coefficients.

You should combine the relative bias in combination with substantive knowledge about the metric of the parameter of interest to determine when levels of relative deviation are negligible or problematic. For example, with a regression coefficient of 0.001, a 5% 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 < |5|%, then do not worry;
• if relative deviation > |5|%, then rerun with 4x nr of iterations.

Question: calculate the relative bias. Are you satisfied with number of iterations, or do you want to re-run the model with even more iterations?

MODEL RESULTS 

                    Estimate       

 DIFF       ON
    AGE                2.143       
    AGE2              -0.021       

 Intercepts
    DIFF             -36.290       

 Residual Variances
    DIFF             196.606    

MODEL RESULTS

                    Estimate      

 DIFF       ON
    AGE                2.141      
    AGE2              -0.021       

 Intercepts
    DIFF             -36.194  
    
Residual Variances
    DIFF             196.343

4. Does the posterior distribution histogram have enough information?

By having a look at the histogram we can check if they contain enough information. To plot the histograms go to Plot > View Plots > Bayesian posterior parameter distributions > view > curve type: histogram > select the parameter you would like to see the histogram of > OK. Also, compare these plots with the ones obatained from your first model.

Question: What can you conclude about distribution histograms?

5. Do the chains exhibit a strong degree of autocorrelation?

To obtain information about autocorrelations you can use (Plot > View Plots > Bayesian autocorrelation plots > view > select the parameter you would like to see the autocorrelation of > OK)

Question: What can you conclude about these autocorrelation plots?

6. Do the posterior distributions make substantive sense?

We plot the posterior distributions and see if they are unimodel (one peak), if they are clearly centered around one value, if they give a realistic estimate and if they make substantive sense compared to our prior beliefs (priors). Here we plot the posteriors of the regression coefficients. To plot the densities go to Plot > View Plots > Bayesian posterior parameter distributions > view > curve type: Kernel densirty > select the parameter you would like to see the density of > OK

Question: What is your conclusion; do the posterior distributions make sense?

7. Do different specification of the variance priors influence the results?

To understand how the prior on the residual variance impacts the posterior results, we compare the previous model with a model where different hyperparameters for the Inverse Gamma distribution are specified. To see what default priors Mplus uses we can have a look at the Appendix A of this manual: “The default for variance parameters, i.e., variance covariance blocks of size 1 is \(IG(-1,0)\)”. in theTECH8output this is summarized as:

  PRIORS FOR ALL PARAMETERS            PRIOR MEAN      PRIOR VARIANCE     PRIOR STD. DEV.

     Parameter 4~IG(-1.000,0.000)          infinity            infinity            infinity

So far we have used the –\(\in \sim IG(.5, .5)\) prior, but we can also use the –\(\in \sim IG(.01, .01)\) prior and see if doing so makes a difference. To quantify this difference we again calculate a relative bias. We change the hyperparameters in the regression model that was specified in step 2 using the ~ command. After rerunning the analysis we can then calculate a bias to see the influence of different priors

Question: Are the results robust for different specifications of the prior on the residual variance (compute the relative bias)?

MODEL PRIORS:
  Beta_Age ~ N(.8, 5);
  Beta_Age2 ~ N(0, 10);
  intercept ~  N(-35, 20);
  e ~ IG(.01, .01);

MODEL RESULTS

                                
                    Estimate    
 DIFF       ON
    AGE                2.143      
    AGE2              -0.021       

 Intercepts
    DIFF             -36.287      

 Residual Variances
    DIFF             197.186      

MODEL RESULTS

                                Posterior  One-Tailed         95% C.I.
                    Estimate       S.D.      P-Value   Lower 2.5%  Upper 2.5%  Significance

 DIFF       ON
    AGE                2.143       0.208      0.000       1.740       2.561      *
    AGE2              -0.021       0.003      0.000      -0.026      -0.015      *

 Intercepts
    DIFF             -36.290       4.187      0.000     -45.020     -28.414      *

 Residual Variances
    DIFF             196.606      15.468      0.000     168.534     227.870      *

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

To see what default priors Mplus uses we can have a look at the Appendix A of this manual. “The default prior for intercepts, regression slopes, and loading parameters is \(\beta \sim \mathcal{N}(0, \infty)\).” This means in practice that we have a normal prior with a super wide variance, and allow probability mass across the entire parameter space. The default priors can also be found in the TECH8 output of Mplus

   PRIORS FOR ALL PARAMETERS            PRIOR MEAN      PRIOR VARIANCE     PRIOR STD. DEV.

     Parameter 1~N(0.000,infinity)           0.0000            infinity            infinity
     Parameter 2~N(0.000,infinity)           0.0000            infinity            infinity
     Parameter 3~N(0.000,infinity)           0.0000            infinity            infinity
     Parameter 4~IG(-1.000,0.000)          infinity            infinity            infinity

Question: What is your conclusion about the influence of the priors on the posterior results?

To anwser this, rerun the model without the prior specifications in the input file and compute the relative bias.

MODEL RESULTS

                               
                    Estimate      

 DIFF       ON
    AGE                2.672       
    AGE2              -0.026       

 Intercepts
    DIFF             -47.409      

 Residual Variances
    DIFF             198.461    

9. Are the results stable from a sensitivity analysis?

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

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 unstable 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.”

For more information on this topic, please also refer to this paper.

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

For a summary on how to interpret and report models, please refer to https://www.rensvandeschoot.com/bayesian-analyses-where-to-start-and-what-to-report/

MODEL RESULTS

                                Posterior  One-Tailed         95% C.I.
                    Estimate       S.D.      P-Value   Lower 2.5%  Upper 2.5%  Significance

 DIFF       ON
    AGE                2.143       0.208      0.000       1.740       2.561      *
    AGE2              -0.021       0.003      0.000      -0.026      -0.015      *

 Intercepts
    DIFF             -36.290       4.187      0.000     -45.020     -28.414      *

 Residual Variances
    DIFF             196.606      15.468      0.000     168.534     227.870      *

In the current model we see that:

• The estimate for the intercept is -36.29 [-45.02; -28.41]
• The estimate for the effect of \(age\) is 2.14 [1.74 ; 2.56]
• The estimate for the effect of \(age^2\) is -0.021 [-0.026; -0.015]

We can see that none of 95% Posterior HPD Intervals for these effects include zero, which means we are can be quite certain that all of the effects are different from 0. Futhermore, yhe third column of the MODEL RESULTS table gives the one-tailed p-value based on the posterior distribution. For a positive estimate, the p-value is the proportion of the posterior distribution that is below zero. For a negative estimate, the p-value is the proportion of the posterior distribution that is above zero. So the a p- value of 0.000 means a probability is around >.999 that the corresponding parameter is not 0.

References

Depaoli, S., & Van de Schoot, R. (2017). Improving transparency and replication in Bayesian statistics: The WAMBS-Checklist. Psychological Methods, 22(2), 240.

Link, W. A., & Eaton, M. J. (2012). On thinning of chains in MCMC. Methods in ecology and evolution, 3(1), 112-115.

van Erp, S., Mulder, J., & Oberski, D. L. (2017). Prior sensitivity analysis in default Bayesian structural equation modeling.

Van de Schoot, R., & Depaoli, S. (2014). Bayesian analyses: Where to start and what to report. European Health Psychologist, 16(2), 75-84.