A default Bayesian hypothesis test for mediation

Abstract

In order to quantify the relationship between multiple variables, researchers often carry out a mediation analysis. In such an analysis, a mediator (e.g., knowledge of a healthy diet) transmits the effect from an independent variable (e.g., classroom instruction on a healthy diet) to a dependent variable (e.g., consumption of fruits and vegetables). Almost all mediation analyses in psychology use frequentist estimation and hypothesis-testing techniques. A recent exception is Yuan and MacKinnon (Psychological Methods, 14, 301–322, 2009), who outlined a Bayesian parameter estimation procedure for mediation analysis. Here we complete the Bayesian alternative to frequentist mediation analysis by specifying a default Bayesian hypothesis test based on the Jeffreys–Zellner–Siow approach. We further extend this default Bayesian test by allowing a comparison to directional or one-sided alternatives, using Markov chain Monte Carlo techniques implemented in JAGS. All Bayesian tests are implemented in the R package BayesMed (Nuijten, Wetzels, Matzke, Dolan, & Wagenmakers, 2014).

Appendixes

Appendix 1. JAGS code

JAGS code for correlation

JAGS code for partial correlation

Appendix 2. Testing the correctness of our JAGS implementation

To assess the correctness of our JAGS implementation, we compared the analytical results for the two-sided Bayes factor against the Savage-Dickey density ratio results based on the MCMC samples from JAGS. The distribution that fit the posterior samples best⁴ is the nonstandardized t-distribution with the following density:

$$ p\left(x\left|\nu, \mu, \sigma \right.\right)=\frac{\varGamma \left(\frac{\nu +1}{2}\right)}{\varGamma \left(\frac{\nu }{2}\right)\sqrt{\left(\pi \nu \sigma \right)}}{\left(1+\frac{1}{\nu }{\displaystyle {\left(\frac{x-\mu }{\sigma}\right)}^2}\right)}^{-\frac{\nu +1}{2}}, $$

(28)

with ν degrees of freedom, location parameter μ, and scale parameter σ. With the samples of the parameter of interest, we can estimate ν, μ, and σ and, thus, the exact shape of the distribution and the exact height of the distribution at the point of interest.

We checked the fit of this distribution and the performance of the SD method in a small simulation study. We considered the following sample sizes: N = 20, 40, 80, or 160. We simulated correlational data by drawing N values for X from a standard normal distribution, and conditional on X, we simulated values for Y according to the following equation:

$$ {Y}_i={\beta}_0+\tau {X}_i+\varepsilon, $$

(29)

where the subscript i denotes subject i and τ represents the relation between X and Y. For each of the four sample sizes, we generated 100 data sets, in each of which τ was drawn from a standard uniform distribution.

Next, we tested the correlation in each data set with both the analytical Bayesian correlation test and the SD method with the nonstandardized t-distribution and compared the results. The results are shown in Fig. 3. The figure shows that the proposed SD method performs well: The Bayes factors of the analytical test and the SD method are similar for all sample sizes and correlations.

Fig. 3

Natural logarithm of the Bayes factors for correlation obtained with analytical calculations (x axis) or obtained with the SD method based on a nonstandardized t-distribution (y axis) for different sample sizes (N). The graphs show fewer points as the samples grow larger, because in these situations, there are more extreme Bayes factors that fall outside the axis limits. We restricted the graphs, since it is most important that the lower Bayes factors lie on the diagonal; it is not important whether a Bayes factor is 2,000 or 3,000, since it is overwhelming evidence in any case