Moderating the Mediation Bootstrap for Causal Inference

Moderating the Mediation Bootstrap for Causal Inference

Introduction

This paper analyses various bootstrap inference methods for indirect effects in the simple mediation model. Meditation is concerned with how the effect of a causal variable $X$ on a consequent variable $Y$ is possibly transmitted through an intervening variable $M$. The analysis of mediating processes has a long history that can be traced back to the path analysis introduced by Wright1920 and originally formulated as a statistical hypothesis by Woodworth1928. During the 1950s, mediation analysis as we know it today was developed in the social sciences, with main contributions in psychology, see for instance Rozeboom1956. The seminal paper of baron1986 laid out statistical requirements for detecting a true mediation relationship and made a distinction between mediating and moderating variables. The current approach favored by Hayes and others, see for instance Preacher2007, has shifted the focus, but both approaches employ a set of two regression models, with coefficients that are used to measure the direct and indirect (mediation) effect. The simple mediation model is given by the following bivariate recursive system (i.e. triangular system with diagonal disturbance covariance matrix):

where $y$, $x$ and $m$ denote $n\times 1$ observable vectors, while $u$ and $ v$ are $n\times 1$ non-observed error vectors. The equations could include an intercept, but we assume without loss of generality that all variables, i.e. $y$, $x$ and $m$, are expressed in deviation from their means.\footnote{ Variables in deviation from their means can be obtained as residuals after regression on a constant. In fact one may also use other variables, e.g. observable confounders, in such a preliminary regression. This only affects the degrees of freedom in the $t$-distributions below.} The\ (indirect) mediation effect is the product $\gamma =\theta _{x}\beta _{m}$, which can be estimated by $\hat{\theta}_{x}\hat{\beta}_{m}$. We are interested in constructing confidence intervals for the mediation effect.

A classic method for the construction of a confidence interval of $\gamma $ is based on the asymptotic standard normal approximation of the studentized quantity:

where $SE(\hat{\gamma})$ denotes the standard error of $\hat{\gamma}$, typically obtained by the delta method. An alternative method is to exploit the quantiles of the product of two standard normal distributions directly; see Craig1936 and Aroian1947. Craig1936 showed that this distribution of the product is symmetric with a kurtosis of 6 when the distributions are independent and both have a zero mean, i.e. $\alpha _{x}=\beta _{m}=0$. Using numerical integration, Meeker81 tabulated quantiles of the product of two normally distributed variables.

Asymptotic results may be of limited worth in small samples, especially when the distribution of $\hat{\gamma}$ is highly non-normal. The bootstrap addresses these small-sample issues like non-normality and is currently the preferred method to construct confidence intervals for $\gamma $ in practice; see for instance MacKinnon2007, MacKinnon04, Shrout2002, Montoya2017 and for an earlier contribution bollen1990. The prevalent bootstrap approach is a paired bootstrap, where resamples are drawn from the tuples $\left( y_{i},x_{i},m_{i}\right) $ with replacement. Various influential papers recommend the bias-corrected (BC)\ method, e.g. MacKinnon04 (more than 7,260 citations per May 2, 2022) state \textquotedblleft As a result, the single best method overall was the bias-corrected bootstrap ...", while also noting that \textquotedblleft The bias-corrected bootstrap did have Type I error rates that were above the robustness interval for some parameter combinations ...". Such liberal inference disqualifies the BC method as a valid inference method not only from a theoretical point of view, but also in practice with deviations that can be substantial. When testing is based on these confidence intervals, as is common practice, it is not size correct and overrejection can be substantial as the simulation results show. The validity of bootstrap methods might also be fundamentally problematic given the absence of pivotal statistics and strong dependence on nuisance parameters.

There are many other bootstrap methods, however, that might be valid, including residual, single versus double, or parametric versus non-parametric bootstrap. In addition, various different versions of constructing confidence intervals exist, including the five main approaches discussed in Section (ref). We follow up on MacKinnon04 who conclude that there are several ways to improve the confidence limits worth investigating. We address the question which methods have correct (minimal) coverage rates and which method in the plethora of bootstrap approaches is preferred. The double bootstrap in particular seems an obvious candidate to correct for coverage errors associated with single bootstrap procedures.\ It is computationally more demanding, however, and coverage of such confidence intervals has not been investigated in the mediation setting. So, we analyze the double bootstrap in detail and find that it over-corrects for specific parameter constellations and parts of the sample space. We provide an explanation for this finding.

In order to analyze the bootstrap methods we use known asymptotic results, but also derive a number of new exact (finite-sample) distributional results. Under strict normality of the errors, we derive the joint distribution of $\hat{\theta}_{x}$ and $\hat{\beta}_{m}$ and show that, conditionally on $x,$ the two estimators are independent and unbiased. This implies that $\hat{\theta}_{x}\hat{\beta}_{m}$ is mean unbiased for $\gamma $ and suggests that the commonly used bias correction is redundant and only introduces unnecessary and harmful randomness when constructing confidence intervals, although the BC\ method is based on the median instead of the mean bias. The distribution of\ $\hat{\gamma}$ is shown to be a Mellin convolution of a student-$t$ and normal distribution with a skewness that only disappears if $\theta _{x}$ and $\beta _{m}$ are zero. This distribution has fatter tails than the distribution of the product of two normal random variables due to the student-$t$ distribution.

The outline of the paper is as follows. In Section (ref), we will derive finite-sample properties of $\hat{\gamma}$ assuming errors are normally distributed. Although the attractiveness of the non-parametric bootstrap lies in the fact that no distributional assumptions are required, the stylized Gaussian setup will act as a benchmark: if the bootstrap does not work in this setup, then it will neither work in the non-Gaussian setup. Section (ref) describes the residual and paired bootstrap and the various confidence intervals. Section (ref) contains a more elaborate exposition of the double bootstrap for confidence intervals since it has not yet\ been investigated or applied in the mediation setting. The Monte Carlo results are shown in Section (ref), where the finite-sample results derived in Section (ref) are useful for explaining some of the observed results. Concluding remarks are given in Section (ref).

Relevant Exact Finite-Sample Results

The bootstrap can be interpreted as a simulation method in which population parameters are (implicitly or explicitly) replaced by sample analogs. The residual bootstrap resamples regression residuals and could be carried out easily in the simple mediation model given in ((ref))-((ref)) recursively generating $m$ and $y$, while keeping $x$ fixed; see equations ((ref)) and ((ref)). On the other hand, the paired bootstrap directly resamples from the observed data; see equation ((ref)). The former approach respects the endogeneity of $y$ and $m$, and explicitly treats $x$ as exogenous, whereas the latter approach only does so implicitly. Although we find that the paired bootstrap works similar as the residual bootstrap, the properties are more directly investigated using a parametric bootstrap approach for the residual bootstrap. In order to characterize the parametric bootstrap, we derive the joint finite-sample distribution of $\hat{\theta}_{x}$ and $\hat{ \beta}_{m},$ and the distribution of their product under the normality assumption: $u_{i}\sim N(0,\sigma _{u}^{2})$ independent of $v_{i}\sim N(0,\sigma _{v}^{2})$ for $i=1,...,n$.

The obvious estimator of the indirect effect is the product of $\hat{\alpha} _{x}$ and $\hat{\beta}_{m}$ from the two regressions involving different covariates: $\hat{\theta}_{x}$ conditional on $x,$ and $\hat{\beta}_{m}$ conditional on both $x$ and $m$. Given the independence of the two error terms, the system in ((ref))-((ref)) is called recursive. Therefore, all model parameters can be estimated consistently and efficiently using ordinary least squares (OLS); see inter alia Rothenberg1964. Conditional on $x,$ equation ((ref)) satisfies all classical Gaussian linear regression assumptions and the maximum likelihood estimator for the conditional mean parameters equals the OLS estimator. The same holds for equation ((ref)) conditional on both $x$ and $m$. We therefore have the standard exact distributional results that, conditional on $x$:

and, conditional on both $x$ and $m$, with $X=[x:m]$ an $n\times 2$ matrix:

We show next in Proposition (ref) that the estimators $\hat{\theta}_{x}$ and $\hat{\beta}_{m}$ are independent, which leads to analytical and numerical simplifications, and that the estimator $ \hat{\beta}_{m}$ itself has a scaled $t_{(n-2)}$-distribution, rather than, as usual, its $t$-ratio. Since the student-$t$ distribution has fatter tails than the normal distribution, especially when the degrees of freedom are small, inference based on the product of two normal distributions can be misleading.

The probability density function (pdf) of $\hat{\gamma}=\hat{\alpha}_{m}\hat{ \beta}_{m}$ can in principle be derived using:

see e.g. Moodetal1974Intro, but this does not lead to a closed-form expression. Equation ((ref)) nevertheless useful provides a convenient way to numerically determine the density and its associated probabilities. We use equation ((ref)) in Section (ref) for the parametric bootstrap distribution with $ (\theta _{x},\beta _{m},\sigma _{v}^{2},\sigma _{u}^{2},x^{\prime }x)$ evaluated at $(\hat{\theta}_{x},\hat{\beta}_{m},s_{v}^{2},s_{u}^{2},x^{ \prime }x)$.

The next proposition\ gives the first three moments of the estimator $\hat{ \gamma}=\hat{\theta}_{x}\hat{\beta}_{m}$. In particular, equation ((ref)) shows that $\hat{\gamma}$ is mean unbiased. This does not imply that it is also median unbiased, however, due to the skewness given in equation ((ref)). The distribution of $\hat{\gamma}$ is skewed when $\theta _{x}\neq 0$ and $\beta _{m}\neq 0$, although the distributions of the individual estimators $\hat{\theta}_{x}$ and $\hat{\beta }_{m}$ are both symmetric. If $\theta _{x}\beta _{m}>0$, the distribution of $\hat{\gamma}$ positively skewed, while it is negatively skewed if $\alpha _{x}\beta _{m}<0$. Note that the BC method uses the median to bias correct the confidence interval.

The expression for the variance is new to the literature. It can be decomposed in terms of different orders. Assuming $x^{\prime }x=O(n)$, or in probability if $x_{i}$ is i.i.d., the first two terms in ((ref)) are of size $O(n^{-1})$ when $\theta _{x}\neq 0$ and $\beta _{m}\neq 0$. The last term in ((ref)) is of smaller magnitude $O(n^{-2})$, and is always larger than zero, even when $ \theta _{x}=\beta _{m}=0$.

The next proposition considers the distributions of the appropriate centered $t$-statistic for $\theta _{x}$ under $H_{0}:\theta _{x}=\theta _{x}^{0}$ and $\beta _{m}$ under $H_{0}:\beta _{x}=\beta _{x}^{0}$, where the superscript $0$ indicates the true value. Since the system is recursive, the $t$-statistics for $\theta _{x}$ and $\beta _{m}$ have a $t_{n-1}$ and $ t_{n-2}$-distribution conditional on the regressors in the model. Although the variance of $\hat{\beta}_{m}$ in the second model depends on the residuals $\hat{u}$ of the first model, the $t$-distributions are still independent from each other as shown in the Proposition (ref) .

Bootstrap Inference

We consider two main bootstrap approaches that are used in the regression model: (i) the paired bootstrap proposed by Efron79 and (ii) the residual bootstrap first analyzed in Bickel/Freedman81. The paired bootstrap is the one generally applied in papers on mediation for bootstrap inference. The paired-bootstrap results are then interpreted conditional on $ x$, but given the causal structure assumed in the mediation setup as $ x\rightarrow m\rightarrow y$, it might be more intuitive to also consider the residual bootstrap; see Hall92 for a more detailed discussion about the different assumptions underlying the paired and residual bootstrap in a regression context.

If $w_{i}=(y\,_{i},x_{i},m_{i})$ denotes the vector containing the $i$-th observation, then the general idea for constructing bootstrap confidence intervals can be summarized as follows:

There are several ways to construct the bootstrap sample $w_{1}^{\ast },...,w_{n}^{\ast }$ in step 1. The paired bootstrap simply resamples from the original $w_{1},...,w_{n}$ with probabilities:

A residual bootstrap generates $w_{1}^{\ast },...,w_{n}^{\ast }$ by resampling bootstrap errors from the residuals, or from a fitted (parametric) distribution, and subsequently constructing $w_{i}^{\ast }$ according to the estimated model. So in the mediation model, bootstrap errors $u^{\ast }$ and $v^{\ast }$ are drawn and the bootstrap observables $ w_{i}^{\ast }=(y_{i}^{\ast },m_{i}^{\ast },x_{i})$ constructed using the estimated parameter values as:

In the non-parametric bootstrap, $u_{i}^{\ast }$ and $v_{i}^{\ast }$ are sampled with replacement from the rescaled OLS\ residuals $\sqrt{n/(n-3)} \hat{u}_{i}$ and $\sqrt{n/(n-2)}\hat{v}_{i}$ respectively for $i=1,...,n$, with the rescaling as originally suggested by Efron82. In the parametric bootstrap, one might draw $u_{i}^{\ast }\sim N(0,\hat{\sigma} _{u}^{2})$ independent of $v_{i}^{\ast }\sim N(0,\hat{\sigma}_{v}^{2})$, using the estimated variances, instead of resampling residuals. The residual bootstrap allows for clear-cut conditioning by keeping $x$ fixed\ in ((ref)) and ((ref)).

The following main methods for constructing confidence intervals have been \ presented in the bootstrap literature: (i) basic (ii) percentile (iii) bias-corrected (BC) percentile (iv) bias-corrected and accelerated (BC$_{a}$) and (v) percentile-$t$ methods.

(i) The basic method for constructing a two-sided equal-tailed $ (1-\alpha )$ confidence interval is based on the idea that the distribution of $\hat{\gamma}-\gamma $ can be approximated by $\hat{\gamma}^{\ast }-\hat{ \gamma}$ leading to the following interval:

where $q_{\alpha }^{\ast }$ denotes the $\alpha $-quantile of $\hat{\gamma} ^{\ast }-\hat{\gamma}$, i.e. $\mathbb{P}^{\ast }[\hat{\gamma}^{\ast }-\hat{ \gamma}\leq q_{\alpha }^{\ast }]=\alpha $; see davison1997bootstrap. Due to the fact that $q_{\alpha }^{\ast }\equiv ( \hat{\gamma}^{\ast }-\hat{\gamma})_{\alpha }=\hat{\gamma}_{\alpha }^{\ast }- \hat{\gamma}$ with $\hat{\gamma}_{\alpha }^{\ast }$ the $\alpha $-quantile of the bootstrap distribution, we can write the basic confidence interval alternatively as:

Comparing the quantiles with those of the percentile and percentile-$t$ explains why this interval is also known as the hybrid interval; see e.g. Shao/Tu96. Note that the lower confidence limit is based on the upper tail of the bootstrap distribution, while the upper confidence limit uses the lower tail. Hence, asymmetry in the basic confidence interval is opposite to the asymmetry of the percentile interval.\ This is an attractive feature because if $\hat{\gamma}^{\ast }- \hat{\gamma}$ is positively skewed, this suggests that larger values of $ \hat{\gamma}$ could more easily be generated by smaller values of $\gamma $ than the other way round.

(ii)\ The two-sided percentile confidence interval is given by:

where $\hat{\gamma}_{\alpha }^{\ast }$ denotes the $\alpha $-quantile of the bootstrap distribution, i.e. $\mathbb{P}^{\ast }[\hat{\gamma}^{\ast }\leq \hat{\gamma}_{\alpha }^{\ast }]=\alpha $; see e.g. Efron81.

(iii) Efron81 also introduces the BC percentile method as an improvement to correct for estimation bias. It uses the proportion of bootstrap replications less than the original estimate $\hat{\gamma}$:

where $\Phi ^{-1}(\cdot )$ denotes the inverse function of the standard normal distribution function, $\mathbbm{1\{\cdot \}}$ denotes the indicator function, and $\hat{\gamma}_{b}^{\ast }$ is the $b$-th bootstrap realization. So, $\hat{z}_{0}$ measures the bootstrap approximation of the median bias of $\hat{\gamma}^{\ast }$ in normal units. If exactly half of the $\hat{\gamma}_{b}^{\ast }$ is less than $\hat{\gamma}$, then $ \hat{z}_{0}=0$. Although Proposition (ref) has established that the estimator $\hat{\gamma}$ is mean unbiased, suggesting that $\hat{z}_{0}$ is close to zero, the correction as defined by Efron81 uses the median, which differs from the mean because of the skewness.

(iv) The BC$_{a}$ interval proposed by Efron87 not only corrects for bias, but also for skewness by the so-called acceleration constant $a$. There are various ways to estimate the acceleration constant $ a $, but a commonly used estimate based on jackknife values is given by:

where $\hat{\gamma}_{(-i)}$ denotes the $i$-th jackknife value based on the sample information excluding the $i$-th observation $ w_{i}=(y_{i},x_{i},m_{i})$ and $\bar{\gamma}_{(\cdot )}=1/n\sum_{i=1}^{n} \hat{\gamma}_{(i)}$ the average of the $n$ jackknife values $\hat{\gamma} _{(1)},...,\hat{\gamma}_{(n)}$. This is also the standard implementation used by, inter alia, the R-packages Lavaan and Boot. A two-sided $(1-\alpha )$ BC$_{a}$ confidence interval can now be defined as the interval:

where the quantiles are based on the probabilities:

If the estimated skewness is 0, then $\hat{a}=0$ and the BC$_{a}$ interval reduces to the BC interval.

(v) Finally, a percentile-$t$ confidence interval is defined by:

where $\tau _{\alpha }^{\ast }$ denote the $\alpha $-quantile of the studentized root $(\hat{\gamma}^{\ast }-\hat{\gamma})/SE(\hat{\gamma}^{\ast })$; see Efron82.

Only the BC$_{a}$ and the percentile-$t$ methods are second-order accurate and are said to achieve asymptotic refinement, see for instance Hall92. However, the accuracy of this latter method in practice depends on the accuracy of the standard error $SE(\hat{\gamma})$. Note that a well-behaved standard error for $\hat{\gamma}$ is problematic; see for instance simulation evidence in mackinnon2002 for a variety of choices. The usual sobel1982 formula:

is used when reporting results for the percentile-$t$ method. Results are also reported based on the Jacknife standard error:

which is related to the denominator of the estimated acceleration constant $ a $ shown in ((ref)). Note that the asymmetry of the percentile-$t$ method is in the same direction as the basic/hybrid method.

Double Bootstrap Methods

Simulation results reported in the mediation literature, for instance MacKinnon04 and its follow-up study Fritzetal2012, show that bootstrap confidence intervals are liberal, i.e. the probability coverage is larger than the nominal $1-\alpha $ coverage, when the indirect effect $ \gamma $ is small. This does not lead to invalid inference, but to confidence intervals that are too wide and therefore to very low probabilities of rejecting the null of no mediation. In particular, when testing the null hypothesis of no mediation, $H_{0}:\gamma =0$, by checking whether the value zero is included by the confidence interval, a liberal interval leads to a low rejection probability of the null. In fact rejection probabilities are very much lower than the significance level, which is a serious problem given that establishing a mediation effect is usually the primary purpose of this type of analysis.

When confidence intervals are liberal, a second-level bootstrap can possibly be used to estimate the overcoverage and correct for it. Such a procedure is called a double bootstrap and, despite its great potential, has hardly been investigated in the mediation setting and principal reason to investigate it here. The main idea is to adjust the quantiles used in the confidence intervals. We use the one-sided percentile method to illustrate the approach. Let $\mathcal{I}_{1}(\alpha ;\mathcal{X},\mathcal{X}^{\ast })=(-\infty ,\hat{\gamma}_{1-\alpha }^{\ast })$ denote the original percentile interval based on sample information $\mathcal{X}$, and resample information $\mathcal{X}^{\ast }$, as a function of the nominal coverage $ 1-\alpha $. The true coverage probability, denoted $\pi (\alpha )=\mathbb{P} [\gamma \in \mathcal{I}_{1}(\alpha ;\mathcal{X},\mathcal{X}^{\ast })]$ could differ significantly from $1-\alpha $. Let $\delta _{\alpha }$ denote the `correct nominal' coverage such that, when used in the procedure, has $\pi (\delta _{a})=1-\alpha $. In general, $\delta _{\alpha }$ is unknown since $ \pi (\alpha )$ is unknown, but $\pi (\alpha )$ can be estimated by:

In practice, this is estimated by simulation, based on $B$ observed first-level bootstrap samples $\mathcal{X}_{1}^{\ast },...,\mathcal{X} _{B}^{\ast }:$

However, the distribution of $\mathcal{X}_{b}^{\ast \ast }$ given $\mathcal{X }_{b}^{\ast }$ is generally unknown, but can be estimated in turn by $C,$ second-level bootstrap samples $\mathcal{X}_{b1}^{\ast \ast },...,\mathcal{X} _{bC}^{\ast \ast }$ from $\mathcal{X}_{b}^{\ast }$. If $\hat{\delta}_{\alpha }$ solves $\hat{\pi}_{B}(\hat{\delta}_{\alpha })=1-\alpha $, then the double bootstrap confidence interval for $\gamma $ is $\mathcal{I}_{2}(\hat{\delta} _{\alpha };\mathcal{X},\mathcal{X}^{\ast })$.

Following davison1997bootstrap, define:

The simulation-based approximation to $\hat{\delta}_{\alpha }$ is given by:

where $\tilde{u}_{[1]}^{\ast }\leq ...\leq \tilde{u}_{[B]}^{\ast }$. Finally, the simulation-based approximation of the double bootstrap percentile interval is given by:

Similarly, for the basic confidence interval in ((ref)), we only have to modify $\tilde{u}_{b}^{\ast }$ to:

so that the double bootstrap basic interval is given by:

Two-sided intervals can be obtained by the set difference of two one-sided intervals, i.e. $\mathcal{J}_{1}(\alpha _{1},\alpha _{2})=\mathcal{I} _{1}(1-\alpha _{1}/2)\backslash \mathcal{I}_{1}(\alpha _{2}/2)$, leading to the following two-sided double bootstrap intervals:

where $\tilde{\delta}_{\alpha /2}$ and $\tilde{\delta}_{1-\alpha /2}$ are two appropriate order statistics. It is also possible to consider a $ (1-\alpha )$ two-sided percentile confidence interval via a two-sided version of $\mathcal{I}_{1}(\alpha )$ directly, i.e. $\mathcal{J}_{1}(\alpha )=(\hat{\gamma}_{\alpha /2}^{\ast },\hat{\gamma}_{1-\alpha /2}^{\ast })$; see for instance Lee1996, although they use $\alpha $ to denote the coverage instead of $(1-\alpha )$.

Figure (ref) shows the simulated coverage based on the double bootstrap of the two-sided percentile for one particular sample: $(\alpha _{x},\beta _{m})=(0,0)$, $(\hat{\theta}_{x},\hat{\beta}_{m})=(0.051,-0.070)$ and $n=50$. The double bootstrap in this case suggests using a 77% confidence level $(\tilde{\delta}_{\alpha }=0.77)$ to obtain a 95% confidence interval. Given that $(\hat{\theta}_{x},\hat{\beta}_{m})$ is close to $(0,0)$, it is reasonable to shorten the confidence interval.

The total number of bootstrap replications equals $B\cdot C$, where $B$ and $ C$ are the number of first- and second-level bootstrap simulations. To ease the computational burden, we exploit the fact that $\tilde{\delta}_{\alpha /2}$ and $\tilde{\delta}_{1-\alpha /2}$ are based on quantiles in both tails of $\tilde{u}^{\ast }$; see also Nankervis2005. Note that $\tilde{u} _{b}^{\ast }$ can be interpreted as a $p$-value. If the bootstrap distribution is centered around $\hat{\gamma}_{b}^{\ast }$, we expect $ \tilde{u}_{b}^{\ast }$ to be large/small when $\hat{\gamma}_{b}^{\ast }$ is far to the left/right of $\hat{\gamma}$. Hence, after sorting $\hat{\gamma} _{b}^{\ast }$, we only carry out the double bootstrap for the $\tfrac{1}{2}M$ smallest and $\tfrac{1}{2}M$ largest values of $\hat{\gamma}_{b}^{\ast }$. Since only $M,$ instead of $B$ values of $\tilde{u}_{b}^{\ast }$ are determined, $\tilde{\delta}_{\alpha /2}$ and $\tilde{\delta}_{1-\alpha /2}$ are based on the $(B/M\cdot \alpha /2)$ and $(1-B/M\cdot \alpha /2)$ quantiles of $\tilde{u}^{\ast }$. In this way, the number of bootstrap replications is reduced from $B\cdot C$ to $M\cdot C$.

Simulation Setup and Results

We have simulated the various bootstrap methods extensively using the \raisebox{-.2\height} programming language, see bezanson17julia. The chosen sample sizes $ n\in \{25,50,100,500\}$ are broadly relevant in various subject areas. The $ x $-vector is drawn from a standard normal distribution, rescaled to have a sample variance of 1 and kept fixed in all simulations since inference is conditional on $x$. The errors $v_{i}$ and $u_{i}$ are independently drawn from the standard normal distribution. For the parameters $\theta _{x}$ and $ \beta _{m}$, we follow MacKinnon04 to indicate the strength of the effects: $0.0$ (none), $0.14$ (small), $0.39$ (medium) and $0.59$ (large) and $\beta _{m}\geq \theta _{x}$. The number of Monte Carlo simulation, $REP$ , is set to 10,000. Since the BC$_{a}$ and double bootstrap intervals adjust the levels of the quantiles, the number of bootstrap replications is taken higher than the usual $1,000$; see e.g. Booth1998. So for each sample, the bootstrap distribution is based on $B=1,999$ first-level bootstrap samples, while the $M=1,000$ second-level bootstrap $p$-values are based on $C=1,000$ second-level bootstrap samples; see Figure (ref) for an illustration. Hence, the bootstrap $p$-values $\tilde{u} _{b}^{\ast }\in \{0.0\%,0.1\%,...,99.9\%,100.0\%\}$ are multiples of $0.1\%$ and $\tilde{\delta}_{\alpha /2}(B+1)$ and $\tilde{\delta}_{1-\alpha /2}(B+1)$ are integers. In this way, no interpolation is needed when constructing double-bootstrap confidence intervals; see Halletal2000.

We report the percentage that confidence intervals are to the left and to the right of $\gamma $ and therefore do not contain the true value. We refer to them as non-coverage rejection frequencies (ncRFs). For $95\%$ equal-tailed confidence intervals, these percentage points should equal $ 2.5\%$, but only approximately, due to simulation error. Given that each ncRF is based on 10,000 trials, it is not significantly different from $ 2.5\% $ (at the 95% confidence level) if its value is contained in the interval:

The ncRFs for $n=100$ are shown in Table (ref), while the results for the other sample sizes can be found in tables (ref) - (ref) in the Appendix. An asterisk (*) after a ncRF indicates that it is significantly different from $ 2.5\%$.

We begin with the results of the percentile method. In line with earlier findings in the literature, the ncRFs of this method when $\theta _{x}=0$ and $\beta _{m}$ small are extremely low. Figure (ref) illustrates this fact by showing the 10,000 estimated values of $(\hat{\alpha }_{x},\hat{\beta}_{m})$ in the simulation for $n=100$ and $\theta _{x}=\beta _{m}=0$ as dots. They are colored red if the residual bootstrap interval based on this realization does not contain the true value $\gamma =0$. There are 18 red realizations out of 10,000. This corresponds to a ncRF of $0.18\%$ and far lower than the nominal $5.0\%$, which has serious consequences for the power of the test based on this confidence interval.

In order to explain this poor behavior, we plot, in\ the same figure, green lines as the boundary of an area having $95\%$ probability of $(\hat{\alpha} _{x},\hat{\beta}_{m})$ lying inside. This is based on the exact distribution of $\hat{\gamma}=\hat{\theta}_{x}\hat{\beta}_{m},$ shown in Figure (ref) as dashed green lines, and it is determined by numerical integration using equation\ ((ref)). Quantiles for $\hat{ \gamma}$ are $\pm 0.02236$ such that:

The green boundary lines in Figure (ref) are the restriction $|\hat{\theta}_{x}\hat{\beta}_{m}|=0.02236$ and essentially 5% of the realizations lie outside it. For the vast majority of this 5% of $( \hat{\theta}_{x},\hat{\beta}_{m})$ realizations with $|\hat{\theta}_{x}\hat{ \beta}_{m}|>0.02236$, the bootstrap generates percentile confidence intervals for $\gamma $ that includes 0. The red colored dots are the values for which the percentile intervals exclude the true value $\gamma =0$. It is clear that this is nowhere near 5%: only $0.18\%$ exclude $\gamma =0$ and coverage is $99.82\%$ instead of the nominal $95\%$. The reason is that the bootstrap distributions for $\hat{\gamma}\neq 0$ are very asymmetric and changes substantially with the estimated parameters, both in location and in their shape. When a sample is drawn and $(\hat{\theta}_{x},\hat{\beta}_{m})$ calculated, the bootstrap distribution approximates the pdf of $\hat{\alpha} _{x}\hat{\beta}_{m}$ with parameter values $(\hat{\theta}_{x},\hat{\beta} _{m},s_{v}^{2},s_{u}^{2})$, which differs substantially from their true values $\left( 0,0,1,1\right) $. This dependence on the parameters can be extreme, as seen in the asymptotic distribution of Sobel's test statistic: if $\gamma =0$ and $\left( \theta _{x},\beta _{m}\right) =(0,0)$ then the asymptotic distribution of the Wald statistic for testing $\gamma =0$ is $ \frac{1}{4}\chi _{1}^{2}$, but $\chi _{1}^{2}$ if $\gamma =0$ and $\left( \theta _{x},\beta _{m}\right) \neq (0,0)$; see glonek1993.

Figure (ref) shows three distributions: the distribution of $ \hat{\gamma}$ for the true parameter values $(0,0,1,1)$ as dashed green lines, the parametric bootstrap distribution of $\hat{\gamma}$ based on equation ((ref)) with parameter values $(\hat{\theta}_{x}, \hat{\beta}_{m},s_{v}^{2},s_{u}^{2})$ in purple, and very close to it, the non-parametric (residual) bootstrap distribution for that particular sample as a light blue histogram. These last two are for one particular realization $(\hat{\theta}_{x},\hat{\beta} _{m},s_{v}^{2},s_{u}^{2})=(0.2216,0.2477,0.9668,1.0913)$ which is the purple star in Figure (ref). When $\left( \theta _{x},\beta _{m}\right) =\left( 0,0\right) ,$ the true distribution of $\hat{\gamma}$ is symmetric, but one will always estimate $(\hat{\theta}_{x},\hat{\beta} _{m})\neq \left( 0,0\right) .$ This will lead to an asymmetric bootstrap distribution with skewness as in formula ((ref)) with $\left( \theta _{x},\beta _{m}\right) =(\hat{\theta}_{x},\hat{\beta}_{m})$.

The red lines in Figure (ref) are determined such that the appropriate limit of the parametric bootstrap confidence interval equals $ \gamma $ (=0). For $\hat{\gamma}>0$, these lines represent the values for $( \hat{\theta}_{x},\hat{\beta}_{m})$ such that the 2.5%-quantile of $\hat{ \gamma}^{\ast }$ is $\gamma =0$ based on numerical integration of equation ( (ref)) with $(\theta _{x},\beta _{m},\sigma _{v}^{2},\sigma _{u}^{2})=(\hat{\theta}_{x},\hat{\beta}_{m},1,1)$, i.e. $(\hat{\theta}_{x}, \hat{\beta}_{m})$ such that

We expect all dots outside the red lines, away from the origin, to not include the true value $\gamma =0$, since the distribution used to determine these boundaries could be interpreted as the parametric bootstrap distribution with knowledge of the nuisance parameters $(\sigma _{v}^{2},\sigma _{u}^{2})$. If the true parametric bootstrap distribution does not significantly depend on the values of the estimated nuisance parameters, we expect the bootstrap distribution based on the true nuisance values to be an accurate approximation. This is indeed the case, since most dots outside the red lines are colored red. The difference between the green and the red lines is that the green line is based on the\ quantiles of $\hat{ \gamma}$ for the single point $(\theta _{x},\beta _{m})=(0,0)$ and the red lines are based on quantiles of the parametric bootstrap distribution of $ \hat{\gamma}^{\ast }$ for all possible values of $(\hat{\theta}_{x},\hat{ \beta}_{m})$. The root of the coverage problem of the percentile method is the extreme dependence of the shape, rather than location, of the bootstrap distribution of $\hat{\gamma}^{\ast }$ on $(\hat{\theta}_{x},\hat{\beta} _{m}) $.

We continue with the results for the percentile interval for $\gamma \neq 0$ : the ncRFs are asymmetric, such that the ncRF is higher on the left of $ \gamma $ than on the right of $\gamma $. This asymmetry becomes less as the sample size increases. Even for $n=100$ ($500)$, most (half) of the ncRFs are significantly different from $2.5\%$. The sum of ncRFs, referred to as total ncRFs, should be around $5\%$, i.e. inside $(4.573\%,5.423\%)$ based on $10,000$ simulations (with 95% confidence). We see that only for $n=500$ that ncRFs are not significantly above $5\%$. Comparing the residual with the paired bootstrap, we observe that the ncRFs for the residual bootstrap are somewhat closer to $2.5\%$ than the paired bootstrap. Hence, exploiting the correct structure as done by the residual bootstrap is noticeable, but the improvement is marginal.

Next, we discuss the results for the basic interval. The results for these intervals are worse than reported for the percentile intervals: in general the ncRFs are lower/higher for small/large values of $\gamma $ and they are also more asymmetric. For instance, when $(\theta _{x},\beta _{m})=(.0,.59)$ and $n=100$, the ncRFs of the basic interval are 1.0 and 1.3 compared to 2.3 and 3.0 of the percentile interval.

The BC and BC$_{a}$ intervals, advocated by inter alia MacKinnon04, do not seem to perform much better than the percentile intervals, but we confirm their finding that they are liberal, i.e. coverage rates below the required 95%. When these intervals are used in testing, the overrejection is clear for $\theta _{x}$ small/medium and $\beta _{m}$ medium: for $ (\theta _{x},\beta _{m})=(0.14,0.14)$ and $n=100$, a ncRF larger than $8\%$ is found for the paired bootstrap and $7.6\%$ for the residual bootstrap. There is hardly any difference between the BC and BC$_{a}$ intervals, due to the estimated acceleration constants in a small interval around 0.

Although the percentile-$t$ intervals theoretically improve an order of magnitude upon the accuracy of the percentile intervals, the ncRFs in the simulation vary substantially with $\gamma $: e.g. for $n=100,$ when $\gamma $ is small, the ncRFs are significantly smaller than $2.5\%$ (but conservative intervals do not violate the stated 95%), for large values of $ \gamma $ ncRFs are close to $2.5\%$, but for the intermediate values $ (\theta _{x},\beta _{m})=(0.14,0.14)$ the total ncRFs are over $18\%$ for the paired bootstrap and $17.7\%$ for the residual bootstrap. These coverage rates worse than $82.3\%$ (instead of the required $95\%$) disqualify the percentile-$t$ method. No substantial difference is observed between the percentile-$t$ based on Sobel's or Jackknife standard errors. Apparently neither one is able to appropriately standardize $(\hat{\gamma}^{\ast }-\hat{ \gamma})$ and turn it into a proper pivotal root.

Finally, the double bootstrap results show that this method, in spite of the promising results in other applications reported in the literature, see for instance Shi1992, Letson1998, McKnight2000, Chronopoulos2015, and Montoya2017, is not able to make the required adjustments. In fact, for medium values of $\gamma $, the second-level bootstrap seems to aggravate the high left ncRFs for the percentile method. The effect on the basic method appears to be even larger.

To investigate this unexpected behavior, Figure (ref) shows the double-bootstrap correction as function of $(\hat{\theta}_{x},\hat{\beta}_{m})$ for $n=100$. The double bootstrap percentile interval can be written as $(\hat{\gamma}_{\alpha _{1}}^{\ast },\hat{\gamma}_{\alpha _{2}}^{\ast })$, where $\alpha _{1}$ and $ \alpha _{2}$ denote the $2.5\%$ and $97.5\%$ percentiles of the $p$-values $ u^{\ast }$ based on the double bootstrap. For a grid of $(\hat{\theta}_{x}, \hat{\beta}_{m})$-values, the parametric bootstrap approximation, assuming $ (\sigma _{v}^{2},\sigma _{u}^{2})=(1,1),$ is used to determine the values of $\alpha _{1}$ and $\alpha _{2}$. When $(\hat{\theta}_{x},\hat{\beta}_{m})$ is close to the origin, there is a substantial double bootstrap correction: $ \alpha _{1}$ is close to $15\%$ and $\alpha _{2}$ close to $85\%$. So, when $ (\hat{\theta}_{x},\hat{\beta}_{m})$ is close to the origin, the 95% double-bootstrap percentile interval uses $(\hat{\gamma}_{0.15}^{\ast },\hat{ \gamma}_{0.85}^{\ast })$, which is much smaller than the single bootstrap percentile interval $(\hat{\gamma}_{0.025}^{\ast },\hat{\gamma} _{0.975}^{\ast })$: the smaller the interval, the higher the probability that it excludes the true value $\gamma $ leading to a severe increase in the ncRF.

This behavior can be seen in Figure (ref) where values of $(\hat{\theta}_{x},\hat{\beta}_{m})$ that lead to bootstrap confidence intervals that exclude the true value $\gamma $ are shown as red colored dots. The upper scatter plots are for $(\theta _{x},\beta _{m})=(0.14,0.14)$, while the lower ones are for $(\theta _{x},\beta _{m})=(0.14,0.39)$ and $n=100$ in each case. In all scatter plots, the true value $(\theta _{x},\beta _{m})$ is represented by the blue diamond, while the blue lines represent all values $(\hat{\theta}_{x},\hat{\beta}_{m})$ such that their product equals this same true value $\gamma =\alpha _{x}\beta _{m}$. The intervals in the left/right scatter plots are based on single/double bootstrap. In all plots, we observe two clusters of red $(\hat{ \theta}_{x},\hat{\beta}_{m})$-realizations away from the upper right blue line that lead to non-coverage as expected. However, for $(\theta _{x},\beta _{m})=(0.14,0.14)$ there are (white) dots even further away that do not lead to non-coverage. The reason is that $(-\theta _{x},-\beta _{m})$ results in the same $\gamma $ value and and hence the lower left blue line. For the double bootstrap we observe far more red non-coverage points near the origin. The explanation is provided by Figure (ref) which shows the large (over)correction close to the origin and leading to increased probability of excluding the true value as just described.

For comparison we also show the results for $(\theta _{x},\beta _{m})=(0.14,0.39)$ when there are hardly any $(\hat{\theta}_{x},\hat{\beta} _{m})$-values close to the origin and therefore the double-bootstrap correction is negligible. So in conclusion, the double bootstrap is counter productive: it either overcorrects when $(\hat{\theta}_{x},\hat{\beta}_{m})$ is close to the origin, or hardly corrects when $(\hat{\theta}_{x},\hat{\beta }_{m})$ is further afield.

Conclusions

Reliable inference on mediation effects is crucial in empirical research and non-trivial from a statistical perspective. The mediation effect $\gamma $ can be estimated as the product of two estimators $\hat{\theta}_{x}$ and $ \hat{\beta}_{m}$, with a distribution that is highly dependent on the values of the two parameters $\theta _{x}$ and $\beta _{m}$, in particular when one of them is 0. Currently the preferred method of inference is based on resampling joint observations, i.e. the paired bootstrap. It is well known that this method has various problems. Confidence intervals are very conservative when both $\theta _{x}$ and $\beta _{m}$ are small, i.e. their coverage probability is much larger than the nominal $95\%$ coverage. This implies that the intervals are (much) too large. In terms of hypothesis testing, this leads to low power, implying that it is harder to establish statistical significance of mediation effects. It should be noted that, even in the idealized parametric setting, mediation tests such as the Likelihood Ratio and Wald (Sobel) test suffer from extremely low power when the mediation effects are small.

There are many alternatives to the paired bootstrap that we include and investigate in this paper. Some have been investigated previously, but we refine some of the results by detailing different constellations of $\alpha _{x}$ and $\beta _{m}$. In certain cases this leads to opposite conclusions from aggregate analysis where various values combined. In particular when $ \gamma =0,$ e.g. when $\theta _{x}=0,$ the value of $\beta _{m}$ has a very large impact on the various distributions (see section 2), simulation results, and even which method is preferred. For instance, although on the whole, bias correcting the bootstrap may appear a good idea, for medium values of $\beta _{x}$ this actually renders confidence intervals invalid since the coverage is smaller than the stated level. Even smaller than would be acceptable under the Bradley1978 liberal robustness criterion that is sometimes employed.

In this paper, the main question addressed is whether the double bootstrap is able to solve the problem of conservative coverage of the single bootstrap. This iterated bootstrap seems a logical solution, but had not yet been investigated. We show that it overcorrects when parameter estimates are small, which can lead to undercoverage, and for large estimates, hardly corrects at all. Hence it provides no solution. This holds for both the paired- and the residual bootstrap. In a single bootstrap setting the residual bootstrap performs slightly better than the paired bootstrap. The explanation we offer is that it exploits the structure of the model, but that makes it susceptible to misspecification, whereas the paired bootstrap is robust against e.g. heteroskedasticity; see Shao/Tu96.

To analyze and explain different simulation results, the finite-sample distribution of $\hat{\gamma}=\theta _{x}\beta _{m}$ is derived assuming normality of the errors. This distribution is used as a benchmark and is very useful in explaining various findings in the simulation results. The result that $\hat{\beta}_{m}$ conditional on $x$ has a student-$t$ distribution in the simple mediation model is new to the literature.

The research turns out to be a cautionary tale about the appropriate choice of bootstrap to use. The simulations results suggest that only the percentile method based on the single bootstrap is able to control the coverage probability. The bias correction methods seem to introduce unnecessary randomness leading to conservative coverage probability for moderate non-zero value of $\gamma $. Comparing the residual to the paired bootstrap, there is not much difference between the two. The double-bootstrap correction seems to be large, but unfortunately in the wrong direction. Based on graphical methods, it appears that this correction is substantial when $\hat{\theta}_{x}$ and $\hat{\beta}_{m}$ are small. However, small values of $\hat{\theta}_{x}$ and $\hat{\beta}_{m}$ only lead to non-coverage for larger values of $\theta _{x}$ and $\beta _{m}$ that do not require a correction. Stated otherwise, the double-bootstrap correction is large when it should be small and vice versa.

Overall, not every bootstrap provides the panacea that current practice seems to suggest. Moreover, none of the bootstrap methods solves the well-known problem of extreme conservative coverage that leads to extremely low rejection probabilities when testing for mediation when effects are small or inaccurately estimated. A different non-bootstrap solution for this problem is given by vanGarderen2022.