Permutation-based methods for mediation analysis in studies with small sample sizes

Proposed approach 1: permutation test of the IERM

The permutation test of the IERM makes use of the Freedman and Lane approach for permuting the residuals (Freedman & Lane, 1983) from the reduced model. Using their general framework, the IERM approach fits full models for Eqs. (1) and (2) in order to obtain an estimate for the indirect effect (Eq. 3). Two reduced models are fit excluding the parameters used to calculate the indirect effect and residuals from each are calculated and permuted. This approach is intended to break the associations between M and Y, and X and M, without disturbing associations with C, thus creating an estimate of the null distributions for α₁ and γ₃, and of ${\mathrm{\alpha }}_1\ast {\mathrm{\gamma }}_3$.

The proposed IERM approach implements the following steps:

1. Fit the full models: $Y={\hat{\mathrm{\gamma }}}_0+{\hat{\mathrm{\gamma }}}_{1}X+{\hat{\mathrm{\gamma }}}_{2}C+{\hat{\mathrm{\gamma }}}_{3}M+e_Y$ and $M={\hat{\mathrm{\alpha }}}_0+{\hat{\mathrm{\alpha }}}_{1}X+{\hat{\mathrm{\alpha }}}_{2}C+e_M$

2. Estimate the original estimate for the indirect effect, ${\left({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3\right)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$

3. Fit reduced models: $Y={\hat{\mathrm{\gamma }}}_{0\left(r\right)}+{\hat{\mathrm{\gamma }}}_{1\left(r\right)}X+{\hat{\mathrm{\gamma }}}_{2\left(r\right)}C+e_{Y\left(r\right)}$ and $M={\hat{\mathrm{\alpha }}}_{0\left(r\right)}+{\hat{\mathrm{\alpha }}}_{2\left(r\right)}C+e_{M\left(r\right)}$

4. Using the reduced models from Step 3, estimate $\hat{Y}$, $e_{Y\left(r\right)}$, $\hat{M}$, and $e_{M\left(r\right)}$

5. Permute residuals from the reduced models, now labeled ${e_Y}^{\ast }$ and ${e_M}^{\ast }$, n_perm times and for each permutation, calculate $Y^{\ast }=\hat{Y}+{e_Y}^{\ast }$ and $M^{\ast }=\hat{M}+{e_M}^{\ast }$.

6. For each permutation, fit the regression models from step 1, replacing Y and M with Y* and M* respectively. Calculate ${{\hat{\mathrm{\alpha }}}_1}^{\ast }\ast {{\hat{\mathrm{\gamma }}}_3}^{\ast }$.

7. The absolute value of $({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3{)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$ is compared to the distribution of absolute values of ${{\hat{\mathrm{\alpha }}}_1}^{\ast }\ast {{\hat{\mathrm{\gamma }}}_3}^{\ast }$ (for a two-tailed test). The p-value is calculated as the proportion of values in the distribution that have equal or greater magnitudes than $({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3{)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$, that is, whose absolute values are greater than or equal to the absolute value of $({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3{)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$.

Proposed approach 2: permutation supremum test under reduced models

The PSRM also makes use of the Freedman and Lane permutation approach (Freedman & Lane, 1983), and incorporates a supremum testing approach for mediation by Wagner et al. (2017) to evaluate the composite null hypothesis ${\mathrm{\alpha }}_1\ast {\mathrm{\gamma }}_3=0$, ${\mathrm{\alpha }}_1=0$ and ${\mathrm{\gamma }}_3=0$. The PSRM differs from the IERM above in that the individual coefficients along the proposed pathway, ${\hat{\mathrm{\alpha }}}_1$ and ${\hat{\mathrm{\gamma }}}_3$, are tested along with the indirect effect. In order to conclude that mediation is present, ${\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3$, ${\hat{\mathrm{\alpha }}}_1$ and ${\hat{\mathrm{\gamma }}}_3$ must be significantly different from zero.

The proposed PSRM differs from the IERM in Steps 2, 6 and 7, and implements the following steps:

1. Fit the full models: $Y={\hat{\mathrm{\gamma }}}_0+{\hat{\mathrm{\gamma }}}_{1}X+{\hat{\mathrm{\gamma }}}_{2}C+{\hat{\mathrm{\gamma }}}_{3}M+e_Y$ and $M={\hat{\mathrm{\alpha }}}_0+{\hat{\mathrm{\alpha }}}_{1}X+{\hat{\mathrm{\alpha }}}_{2}C+e_M$

2. Estimate the original estimates ${{\hat{\mathrm{\alpha }}}_1}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$, ${{\hat{\mathrm{\gamma }}}_3}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$ and $({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3{)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$.

3. Fit reduced models: $Y={\hat{\mathrm{\gamma }}}_{0\left(r\right)}+{\hat{\mathrm{\gamma }}}_{1\left(r\right)}X+{\hat{\mathrm{\gamma }}}_{2\left(r\right)}C+e_{Y\left(r\right)}$ and $M={\hat{\mathrm{\alpha }}}_{0\left(r\right)}+{\hat{\mathrm{\alpha }}}_{2\left(r\right)}C+e_{M\left(r\right)}$

4. Using the reduced models from Step 3, estimate $\hat{Y}$, $e_{Y\left(r\right)}$, $\hat{M}$, and $e_{M\left(r\right)}$

5. Permute residuals from the reduced models, now labeled ${e_Y}^{\ast }$ and ${e_M}^{\ast }$, n_perm times and for each permutation, calculate $Y^{\ast }=\hat{Y}+{e_Y}^{\ast }$ and $M^{\ast }=\hat{M}+{e_M}^{\ast }$.

6. For each permutation, fit the regression models from step 1, replacing Y and M with Y* and M* respectively. Estimate ${{\hat{\mathrm{\alpha }}}_1}^{\ast }\ast {{\hat{\mathrm{\gamma }}}_3}^{\ast }$, ${{\hat{\mathrm{\alpha }}}_1}^{\ast }$ and ${{\hat{\mathrm{\gamma }}}_3}^{\ast }$.

7. The absolute value of $({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3{)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$ is compared to the distribution of absolute values of ${{\hat{\mathrm{\alpha }}}_1}^{\ast }\ast {{\hat{\mathrm{\gamma }}}_3}^{\ast }$ (for a two-tailed test). The p-value is calculated as the proportion of values in the distribution that have equal or greater magnitudes than $({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3{)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$. The absolute value of ${{\hat{\mathrm{\alpha }}}_1}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$ is compared to the distribution of absolute values of ${{\hat{\mathrm{\alpha }}}_1}^{\ast }$ (for a two-tailed test). The p-value is calculated as the proportion of values in the distribution that have equal or greater magnitudes than ${{\hat{\mathrm{\alpha }}}_1}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$. Similarly, a p-value is obtained for ${{\hat{\mathrm{\gamma }}}_3}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$. The null hypothesis is rejected only if $({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3{)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$, ${{\hat{\mathrm{\alpha }}}_1}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$ and ${{\hat{\mathrm{\gamma }}}_3}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$ are significantly different from zero, and the significance level is the supremum of the significance levels of the individual tests.

Approach for comparison: the permutation test of the IEFM

The permutation test for the IEFM makes use of ter Braak’s method (Manly, 1997; ter Braak, 1992), fitting full models for both regression models and estimating the indirect effect. Rather than estimating both full and reduced models as in the PSRM and IEFM, only the full models are fitted in this approach. Residuals are permuted and used to create a sampling distribution of the test statistic and estimate confidence limits for the indirect effect (Koopman et al., 2015; Taylor & MacKinnon, 2012).

The IEFM method implements the following steps:

1. Fit the full models: $Y={\hat{\mathrm{\gamma }}}_0+{\hat{\mathrm{\gamma }}}_{1}X+{\hat{\mathrm{\gamma }}}_{2}C+{\hat{\mathrm{\gamma }}}_{3}M+e_Y$ and $M={\hat{\mathrm{\alpha }}}_0+{\hat{\mathrm{\alpha }}}_{1}X+{\hat{\mathrm{\alpha }}}_{2}C+e_M$

2. Estimate the original estimate for the indirect effect $({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3{)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}$

3. Permute residuals from the full models, now labeled ${e_Y}^{\ast }$ and ${e_M}^{\ast }$, n_perm times and for each permutation, calculate $Y^{\ast }=\hat{Y}+{e_Y}^{\ast }$ and $M^{\ast }=\hat{M}+{e_M}^{\ast }$.

4. For each permutation, fit the regression models from step 1, replacing Y and M with Y* and M* respectively. Calculate ${{\hat{\mathrm{\alpha }}}_1}^{\ast }\ast {{\hat{\mathrm{\gamma }}}_3}^{\ast }$.

5. Confidence bounds are estimated using the $\left({\displaystyle \frac{\mathrm{\omega }}{2}}\right)\ast 100$ and $\left(1-{\displaystyle \frac{\mathrm{\omega }}{2}}\right)\ast 100$ percentiles of the distribution of ${{\hat{\mathrm{\alpha }}}_1}^{\ast }\ast {{\hat{\mathrm{\gamma }}}_3}^{\ast }$, where ω corresponds to the desired alpha level.

6. The null hypothesis $({\hat{\mathrm{\alpha }}}_1\ast {\hat{\mathrm{\gamma }}}_3{)}_{\mathrm{o}\mathrm{r}\mathrm{i}\mathrm{g}}=0$, is rejected if 0 is not contained within the confidence bounds.

Simulation studies

Simulation studies were used to empirically evaluate the behavior of the two proposed approaches (IERM and PSRM) and the existing IEFM method (Koopman et al., 2015; Taylor & MacKinnon, 2012). To ensure scenarios of weak and strong confounding, normally distributed data were generated for X, M, Y and C from a range of correlation structures by multiplying the Cholesky decomposed matrix of the correlation structure with a matrix of independent and random normally distributed values. A subset of all possible correlations of 0, 0.15, 0.3 and 0.6 between the four variables with positive definite correlation structures was evaluated (see Tables S1 and S2 for complete list of scenarios and the corresponding regression coefficients). Correlation strengths were selected to simulate comparable regression coefficients to prior studies on permutation methods for mediation (Taylor & MacKinnon, 2012); the simulated regression coefficients are reported in the Online Supplemental 1. For example, one condition evaluated was a scenario where C was a weak confounder between X and Y, with a correlation structure of $\begin{array}{c}\begin{array}{c}X\\ C\end{array}\\ \begin{array}{c}M\\ Y\end{array}\end{array}\left[\begin{array}{cccc}1& 0.15& 0.15& 0.15\\ \dots & 1& 0& 0.15\\ \dots & \dots & 1& 0.15\\ \dots & \dots & \dots & 1\end{array}\right]$

Three covariate or confounding scenarios were considered in this study: (1) C as a covariate, associated with outcome Y but not with the exposure X or mediator M (Fig. 2A), (2) C as a confounder of the exposure–outcome (X–Y) relationship (Fig. 2B), and (3) C as a confounder of the mediator–outcome (M–Y) relationship (Fig. 2C). For each scenario, C was assessed as either a “weak” or “strong” covariate or confounder, with correlation of 0.15 or 0.6 respectively, resulting in a total of six covariate scenarios evaluated.

Under the null of no indirect effect, type I error rates were evaluated at a significance level of α = 0.05. The null hypothesis of no indirect effect may be simulated from the following three scenarios: (1) both ${\mathrm{\alpha }}_1=0$ and ${\mathrm{\gamma }}_3=0$, (2) ${\mathrm{\alpha }}_1=0$ and ${\mathrm{\gamma }}_3\ne 0$, or (3) ${\mathrm{\alpha }}_1\ne 0$ and ${\mathrm{\gamma }}_3=0$. All three null possibilities were evaluated for C as a weak or strong covariate or confounder. All other correlations were held at 0.15 for the “weak” scenario, and 0.6 for the “strong” scenario, resulting in a total of 27 conditions (13 conditions for scenario a, eight conditions for scenario b, and six conditions for scenario c) per sample size being evaluated under the null.

Under the alternative hypothesis, both ${\mathrm{\alpha }}_1\ne 0$ and ${\mathrm{\gamma }}_3\ne 0$. Correlation between the exposure X and the mediator M, ρ_XM, or correlation between the mediator M and outcome Y, ρ_MY, varied from 0.15, 0.3 and 0.6. For each of the six covariate scenarios, ρ_XM and ρ_MY were varied as 0.15, 0.3 or 0.6, while all other correlations were held at 0.15 for the “weak” scenario, and 0.6 for the “strong” scenario, resulting in a total of 51 conditions (18 conditions for scenario a, 18 conditions for scenario b, and 15 conditions for scenario c) per sample size being evaluated under the alternative. Similar to Koopman et al. (2015), we considered sample sizes of 30 and 100. For each simulation condition, 1,000 replicates were run with 10,000 permutations per replicate. Simulation code is available at https://github.com/kroehlm/Permutation_Mediation_Test.

Type I error

Type I error rates with C as a confounder between the exposure–outcome (X–Y) relationship are shown in Table 1. Results are shown in Table 1 for the null scenario (1) in which no relationship exists between the exposure–mediator (X–M) and mediator-outcome (M–Y) (i.e., ${\mathrm{\alpha }}_1=0$, ${\mathrm{\gamma }}_3=0$) and the null scenario (2) where no relationship exists between the exposure–mediator (X–M) (i.e., ${\mathrm{\alpha }}_1=0$ and ${\mathrm{\gamma }}_3\ne 0$). The null scenario (3) where no relationship exists between mediator-outcome (M–Y) (i.e., ${\mathrm{\alpha }}_1\ne 0$ and ${\mathrm{\gamma }}_3=0$) are available in the Online Supplemental 1. Each row of Table 1 represents a single simulation condition, and for each condition, the table includes the correlation conditions among all variables, the corresponding α₁ and γ₃ coefficient values, and the empirical type I error rates for the three different permutation applications.

Type I rates were generally slightly higher when C was a strong confounder compared to when it was a weak confounder. The results in the first row of Table 1 represent the scenario where C is a weak confounder between the X and Y relationship, the type I error rate for the test of the indirect effect under the reduced models, the IERM approach, was 0.044. The proposed PSRM approach and the existing IEFM approach had more conservative type I error rates of 0.005 and 0.006, respectively. Rates for all three methods increased for conditions with a strong confounder.

The IERM approach had type I rates at the nominal level when both α₁ and γ₃ were zero. When either α₁ or γ₃ was nonzero, the type I error rates exceeded the nominal levels, with rates increasing with increasing size of the nonzero coefficient, and reaching levels ${\mathrm{\gamma }}_3\ne 0$ ${\mathrm{\alpha }}_1=0$ ${\mathrm{\alpha }}_1=0$, ${\mathrm{\gamma }}_3=0$ ${\mathrm{\alpha }}_1=0$ ${\mathrm{\gamma }}_3\ne 0$ ${\mathrm{\alpha }}_1=0$ ${\mathrm{\alpha }}_1=0$ ${\mathrm{\gamma }}_3=0$ ${\mathrm{\gamma }}_3\ne 0$ ${\mathrm{\alpha }}_1=0$ ${\mathrm{\gamma }}_3=0$ ${\mathrm{\alpha }}_1=0$ ${\mathrm{\gamma }}_3\ne 0$ ${\mathrm{\alpha }}_1=0$ ${\mathrm{\gamma }}_3=0$as high as 0.653 for a sample size of 30 and 0.783 for a sample size of 100. Both the proposed PSRM approach and the existing IEFM approach had type I error rates well below the nominal level when both α₁ and γ₃ were zero. When either α₁ or γ₃ was nonzero, the type I error rates increased with increasing size of the nonzero coefficient, and was generally around the nominal level when the nonzero coefficient was large. However, when the sample size was small (i.e., 30), nominal levels were exceeded a handful of times for the IEFM approach and only once for the proposed PSRM approach. For a larger sample size, type I error rates remained below any expected deviations from nominal levels with one exception; the error rates were slightly inflated for the scenario where C was a weak confounder between X and Y, and γ₃ was large for the IEFM test (condition 12 in Table 1) but not for the proposed PSRM approach. Across the 54 null conditions evaluated, type I error rates were exceeded in 3.7% (2) of the conditions for the PSRM, and in 14.8% (8) of the conditions for the IEFM. With a smaller sample size, the IEFM approach occasionally had slightly higher rates than those for the proposed PSRM approach. However, for the larger size of n = 100, the error rates between the two tests were very similar, with neither consistently being larger than the other.

Power

Results for power when C was a confounder between the X and Y relationship are shown in Table 2 (similar results for the other covariate scenarios are in Tables S5 and S6). Due to the highly inflated error rates for the proposed IERM approach, we excluded results of power from this approach. As with the type I error rates, there were no major differences of power based on the covariate scenario within a permutation method. As with Table 1, each row represents a single simulation condition, and for each condition, the table includes the correlations among all variables, the corresponding α₁ and γ₃, and the power for the two permutation approaches. The results in the first row represent the scenario where C is a weak confounder between the X and Y relationship, under the alternative condition that both α₁ and γ₃ are nonzero. As a weak confounder between X and Y, the correlations between X–C and C–Y were simulated to be 0.15, with zero correlation between C and M. The correlation between X and Y was also simulated to be 0.15. The correlations between X–M and M–Y were both simulated to be 0.15, with corresponding coefficients of ${\mathrm{\alpha }}_1=0.1535$ and ${\mathrm{\gamma }}_3=0.1335$. For this condition, the power for the proposed PSRM approach was 0.014, and power for IEFM approach was 0.019.

As expected, power was lowest for small values of a coefficient, and increased with increasing values of the coefficients. For the small sample size (i.e., n = 30), power was typically larger for the IEFM approach, especially for larger coefficient values. In the 51 alternative conditions evaluated for this sample size, power for the PSRM was only equal to the IEFM in one condition, and was, on average, 2.2% lower than power for the IEFM. Once the sample size was increased to 100, however, there were no differences in power between the two tests; see Fig. 3. Power for the PSRM was equal or better than the IEFM in 25 of the 51 conditions for this sample size, with a mean difference between the two methods of 0.31%. Figure 4 displays the change in power of the proposed PSRM approach for both weak and strong confounders, based on their coefficient values, for n = 100 where C was a confounder between the X and Y relationship.

Comparison to bootstrap based approaches for mediation analysis

Prior evaluations of the type I error rates for the IEFM approach were not inflated as indicated by simulation studies (Koopman et al., 2015; Taylor & MacKinnon, 2012). Upon completion of our studies and finding the IEFM and PSRM did exceed nominal rates under a small set of conditions, we conducted a post-hoc study to compare these methods with bootstrap methods under the four most extreme null conditions we observed, as well as six alternative conditions. The purpose of these studies was to provide a side-by-side comparison of the methods reported by Koopman et al. (2015) and the permutation methods evaluated here under identical conditions with weak and strong confounding. Simulations were carried out as described above, and the indirect effect was evaluated by the PSRM, IEFM, and the three bootstrap methods evaluated in Koopman et al. (2015): the percentile bootstrap (PB), the bias corrected bootstrap (BCB), and the bias corrected accelerated bootstrap (BCAB). With respect to type I error, results from these studies (Table 3) indicate the PB approach performs best among the bootstrap methods, while the BCB and BCAB approaches exceeded nominal levels under all four scenarios. Overall, the PSRM approach performed better than the IEFM and PB approaches, maintaining error rates at the nominal level in all but one condition. Comparisons of power among the five methods are presented in Table 4. For all conditions evaluated, power was highest for the BCB and BCAB, the two methods which exhibited inflated type I errors under all four null conditions. Power was slightly higher for the IEFM when compared with the PSRM and PB when n = 30, and comparable among the three methods when n = 100.

Example: how framing of media stories influences attitudes regarding immigration policy

We now illustrate the application of the PSRM and IEFM with an empirical example based on the framing data of Brader, Valentino & Suhay (2008). In this study, participants were randomly assigned to different media stories about immigration with either positive or negative framing and asked about their attitudes and political behavior with respect to immigration policy. They hypothesized that anxiety mediates the relationship between framing and whether a participant would agree to send a letter about immigration policy to his or her member of Congress. Their analysis, controlling for education, age, income, and sex, suggested anxiety does act as a significant mediator between the two variables. Here, we will extend the primary analysis by considering the role of framing on attitudes toward increased immigration, a four-point item with larger vales indicating more negative attitudes. Despite randomization, there was some unevenness of distribution across education groups with respect to treatment exposure, so we stratify analyses by education level. Similar to Brader, Valentino & Suhay (2008) we adjusted for age, income and sex.

Results from our extended analysis suggest that anxiety may mediate the relationship between framing and attitudes toward immigration in some, but not all, education groups. For subjects with a high school education, both the p-value from the PSRM and 95% CI from the IEFM were close to their respective decision points, but did not exceed thresholds to reject the null hypothesis (Table 5). For all education levels, inference between the two permutation methods did not differ.