A Nearly Similar Powerful Test for Mediation

This paper derives a new powerful test for mediation that is easy to use. Testing for mediation effects is empirically extremely important in various scientific disciplines. A key paper in psychology, Baron and Kenny (1986) has more than 100,000 citations¹¹1 Cited by 106,782 on 13 October 2021, 90,147 on 15 January 2020, and 79,205 on 22 October 2018. and is used in many other fields. Mediation testing is important in accounting, e.g. Coletti et al. (2005), marketing, e.g. MacKenzie et al. (1986), sociology, e.g. Alwin and Hauser (1975) who used the expression indirect effect, a term commonly used in economics also.For a recent overview of mediation in economics see Huber (2020) and e.g. Heckman and Pinto (2015a, b) on treatment effects and production technology, and Imbens (2020) who extensively reviews connections between directed acyclic graphs (DAGs), potential outcomes, causal inference, instrumental variables, and mediation. Frewen et al. (2013) is exemplary for the increasing network literature, including DAGs, using mediation. The minimal selection here is hardly representative for the vast body of literature on mediation analysis. It only illustrates the breadth of its empirical relevance. Tests for mediation can have extremely low power, especially when the effect is small, or estimated with large variance. The primary purpose of this paper is to provide a new and more powerful test that is easy to use.

The aim of mediation testing is to discover if an independent variable ($X$) causes a dependent variable ($Y$) via an intervening, or mediating variable ($M$). The mediating variable is exogenous in the common experimental setting in psychology and other fields, but is also considered exogenous in other settings where assignments are random or constitute a natural experiment. The basic model is simply:

where all variables are taken in deviation from their means, or more generally after partialing out other exogenous effects. The disturbances $u$ and $v$ are assumed to be independent because of an experimental set up and more generally because no influence of $Y$ on $M$ is assumed in this type of model. This independence is a crucial identification condition, since the parameter $\theta_{2}$ cannot be estimated consistently if $M$ is endogenous. We make a convenient distributional assumption: $\left(u_{i},v_{i}\right)^{\prime}\sim II\emph{N}\left(0,diag\left(\sigma_{11},\sigma_{22}\right)\right)$, $i=1,\cdots,n,$ with $n$ the number of observations. This facilitates a likelihood analysis, but is not necessary for the asymptotic normality of the $t$-statistics that will be used.

MacKinnon et al. (2002) give a literature review and compare 14 different methods for testing the effects of a mediation variable. These methods are based on standardized measures of the product of two coefficients $\theta_{1}\theta_{2}$ or based on the difference of two related coefficients ($\tau^{\ast}-\tau$) in equations (1) and (3):

If there is a mediation effect, then $X$ influences $M,$ such that $\theta_{1}\neq 0,$ and $M$ influences $Y,$ such that $\theta_{2}\neq 0.$ If there is no mediation by $M,$ then the effect of $X$ on $Y$ is not altered by the inclusion of $M$ such that $\tau^{\ast}-\tau=0.$

Model (3) is a restricted version of (1) with $\theta_{2}=0,$ and it is straightforward therefore to show that the OLS estimates for the three models satisfy $\hat{\tau}^{\ast}=\hat{\tau}+\hat{\theta}_{1}\hat{\theta}_{2}$ and the relation $\tau^{\ast}-\tau=\theta_{1}\theta_{2}$ also holds in model interpretation terms; see Appendix A.

Figure 1 illustrates mediation, no mediation, and the related concept of moderation in terms of directed acyclic graphs. $M$ is a moderator if it changes the relation between $X$ and $Y$. In its basic form this is modeled by adding the interaction term between $X$ and $M$ to the model. When confounders $W$ are observed and included we have:

$W$ can also be included in (1) - (3), but partialed out such that $Y$, $X,$ and $M$ are residuals after regressions on $W.$ Moderation and confoundedness can be tested by the ordinary $t$ or $F$ tests on $\gamma$ and $\delta,$ but testing for mediation is less straightforward.

The best well known and commonly used test for mediation by Sobel (1982) is a Wald-type test of the form $\ \hat{\theta}_{1}\hat{\theta}_{2}/SE(\hat{\theta}_{1}\hat{\theta}_{2})$, with $SE(\hat{\theta}_{1}\hat{\theta}_{2})$ an estimate of the standard error of the product $\hat{\theta}_{1}\hat{\theta}_{2}.$ It is available in standard statistical packages such as SAS, R, Stata, or SPSS. It has good properties when either $\theta_{1}$ or $\theta_{2}$ is large and the standard errors of $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ are small, but if the two $t$-tests for testing $\theta_{1}=0$ and $\theta_{2}=0$ tend to be small, properties deteriorate. For parameter values under the null, the Null Rejection Probability (NRP) can be very close to zero and, under the alternative, power can fall far below the size (highest NRP) of $5\%$ that we use throughout. All tests considered in MacKinnon et al. (2002) suffer from these problems. The distributions of the test statistics considered in the literature depend on the value of the parameters under the null. As a consequence none of these tests is similar, meaning that rejection probabilities are not constant on the boundary of the null hypothesis. In fact rejection probabilities under alternatives close to the origin, i.e. power, can be much lower than size and these tests are biased.

Much effort in the literature has gone into improving well-known test statistics, such as the Wald statistic, without a satisfactory solution. The bootstrap is invalid (see Van Garderen and Van Giersbergen (2021)) and cannot salvage these statistics. The key step in our approach is to move away from test statistics and consider the critical region in the sample space directly and optimize the flexible boundary we define.

We make two main contributions. Our main theoretical contribution in Theorem 2 shows that a similar mediation test exists if and only if the level of the test is the reciprocal of an integer, i.e. $1/\alpha\in\mathbb{N},$ or $\alpha=0$. Hence, for practical levels such as $1\%,$ $5\%$, or $10\%$ an exact similar test exists. The proof is constructive and the test is unique within the class considered. Unfortunately, the critical region is objectionable in that it includes an area near the origin where both $t$-statistics are arbitrarily close to zero. Such values do not provide overwhelming evidence against the null and Perlman and Wu (1999) coined the term “Emperor’s New Tests” for similar tests with undesirable properties like these. Insistence on similarity can render LR tests $\alpha$-inadmissible, cf. Lehmann and Romano (2005, Section 6.7), but Perlman and Wu (1999) give examples where similar tests have extremely undesirable properties, yet inadmissible LR tests still provide reasonable answers. In the mediation setting the LR test is also inadmissible, but does not provide a satisfactory answer in this case. It is much better than the Wald test as we will see, but nevertheless suffers from extremely poor power properties for parameter values close to the origin.

So a better test is called for and our second, more important contribution therefore is practical. We construct a new simple test for mediation that is uniformly more powerful than the LR test without the undesirable properties sketched by Perlman and Wu (1999), and that is nearly similar and practically unbiased.

It is extremely simple to use:

1. 1.

   Order the absolute values of the common $t$-statistics from the basic OLS regressions (1) and (2) for testing $\theta_{1}=0$ or $\theta_{2}=0:$ $\left|t\right|_{(1)}=\min\left\{\left|t_{1}\right|,\left|t_{2}\right|\right\}$, $\left|t\right|_{(2)}=\max\left\{\left|t_{1}\right|,\left|t_{2}\right|\right\}$

2. 2.

   Reject if $\left|t\right|_{(1)}>g(\left|t\right|_{(2)})$ using Table 1 (or employing the code in Appendix E).

That the test can be based on elementary $t$-statistics is important for ease of application, but is theoretically justified below by sufficiency and invariance arguments. The testing problem is invariant to permutations of parameters and statistics and sign changes. We show that the ordered absolute $t$-statistic, $(\left|t\right|_{(1)},\left|t\right|_{(2)}),$ is a maximal invariant. The critical region is a subset of the relevant sample space which is an octant in $\mathbb{R}^{2}$. This can also be justified asymptotically under very weak assumptions and different estimation methods.

The new test is constructed by varying the boundary of the critical region, defined as a function $g$, such that the test is almost similar and minimizing the distance to the power envelope surface. It is based on a new general, so called varying-$g$ method that can be applied to other testing problems with nuisance parameters more generally to obtain near similar tests. It does not require a choice of mixture distribution, nor the construction of least favorable distributions, cf. Andrews and Ploberger (1994), Andrews et al. (2006, 2008), Elliott et al. (2015), Guggenberger et al. (2019). It can be given a random critical value interpretation as in Moreira and Mourão (2016) since any critical region in a higher-dimensional space has a boundary for one statistic in terms of the remaining statistics. The critical region that we construct is fixed however, not random, avoids simulation, and our approach appears to lend itself better to multivariate extensions, as will be shown for dimension three.

We use and develop numerical methods that avoid simulations and use numerical integration instead. With the required computing completed for the mediation problem, practitioners can simply use Table 1 or the computer code provided in the appendix. In fact no further reading on the motivation and derivation of the test is required for its implementation.

Section 5 shows the ease of implementation using an interesting application by Alan et al. (2018) on educational attainment. The test confirms that the negative effect on girls’ attainment of 1-year exposure to teachers with traditional attitudes, is mediated through students’ own gender role beliefs. Neither the procedure in Alan et al. (2018) nor the LR test reject the no mediation hypothesis in this case, but our new test does.

A further empirical illustration on union sentiment among southern nonunion textile workers is provided in Section 7. Different mediation channels are tested involving two mediating variables and requires an extension of our methods. We therefore consider general hypotheses of the form $H_{0}:\theta_{1}\cdots\theta_{K}=0$. We give the relevant distributions of maximal invariants that can be used to derive the critical regions that are nearly similar and do so explicitly for three dimensions.

For practitioners the major advantage of our test is that there is a better chance of formally showing that there is a mediation effect. Our test has better power, especially when the two channeling effects are small or less accurately estimated. Given the enormous interest in testing for mediation and the fact that our test can have close to $5\%$ more power than standard level $5\%$ tests, many unpublished examples will exist where it can now be concluded that there is a statistically significant mediation effect.

The joint density of $\left(Y,M\right)$ given $X$ in equations (1) and (2) can be written as:

with $\lambda_{1}=\left(\tau,\theta_{2},\sigma_{11}\right)^{\prime}$, $\lambda_{2}=\left(\theta_{1},\sigma_{22}\right)^{\prime}$, and $\lambda=\left(\lambda_{1}^{\prime},\lambda_{2}^{\prime}\right)^{\prime}$. The parameters $\lambda_{1}$, and $\lambda_{2}$ vary freely as a result of the triangular structure of the model. The mediation variable is the endogenous variable in (2), but is exogenous for $\theta_{2}$ in (1) since $Y$ is not causal for $M$. For a sample of $n$ independent observations the loglikelihood equals the sum of two normal loglikelihoods corresponding to (1) and (2)²²2This can easily be extended to include more regressors/covariates. Instrumental variables can also be used, but note that $X$ and $M$ appear in both equations and in the standard setup $u$ and $v$ are independent because of the experimental interpretation of $M$. All that is required is the asymptotic normality of the estimators and $t$-statistics.:

As a consequence the Maximum Likelihood Estimators (MLEs) for $\theta_{1}$ and $\theta_{2}$ are the basic OLS estimators for the two equations separately. Furthermore, both observed and expected Fisher information matrices will be block diagonal in terms of $\lambda_{1}$ and $\lambda_{2}$ as well as in $\left(\tau,\theta_{2}\right)^{\prime}$, $\sigma_{11}$, $\theta_{1},$ and $\sigma_{22}.$ As a result the standard $t$-statistics $T_{1}$ and $T_{2}$ for $\theta_{1}$ and $\theta_{2}$ respectively are asymptotically independent and normally distributed with means $\mu_{1}\equiv\theta_{1}^{0}/\sigma_{\theta_{1}}$, $\mu_{2}\equiv\theta_{2}^{0}/\sigma_{\theta_{2}}$ where $\theta_{1}^{0}$, $\theta_{2}^{0}$ denote the true parameter values and $\sigma_{\theta_{1}}$, $\sigma_{\theta_{2}}$ the standard deviations of the OLS estimators: $\ \left(T-\mu\right)\overset{d}{\rightarrow}\emph{N}\left(0,I_{2}\right)$, with $T=\left(T_{1},T_{2}\right)^{\prime}.$

Restricting attention to $T$ can be justified by statistical sufficiency, since the MLE is minimal sufficient and complete, and by invariance since the values of $\sigma_{11},\sigma_{22}$ and $\tau$ do not affect whether $H_{0}:$ $\theta_{1}\theta_{2}=0$ is true or not. Hillier et al. (2021) show that $T$ is a maximal invariant under a relevant group of transformations. The testing problem has two more obvious symmetries. The problem is not affected by sign changes or permutations. This also holds in higher dimensions, e.g. when mediation is through a chain of effects as in our empirical illustration. If $X\rightarrow M^{\left(0\right)}\rightarrow\cdots\rightarrow M^{\left(K-1\right)}\rightarrow Y,$  then $K$ parameters are required to be non-zero for this channel to operate. In $K$ dimensions the null hypothesis that at least one parameter is zero and the alternative is all $K$ parameters non-zero:

There are many other hypotheses including collapsibility in contingency tables, testing indirect effects, channels in DAGs, see Dufour et al. (2017) for a range of examples.

In the multivariate setting we also assume that the estimator $\hat{\theta}$ is normally distributed with known covariance matrix $\Omega=diag\left(\sigma_{1}^{2},...,\sigma_{K}^{2}\right)$. Hence, if we let $T=\Omega^{-1/2}\hat{\theta}$ and $\mu=\Omega^{-1/2}\theta^{0}$ such that $T_{k}=\hat{\theta}_{k}/\sigma_{\theta_{k}}$ are the $t$-ratios and $\mu_{k}=\theta_{k}/\sigma_{\theta_{k}}$ are the non-centrality parameters, we assume:

The testing problem is invariant to reordering the parameters (permutations) and sign changes (reflections) of the $K$ parameters $\left\{\theta_{i}\right\}_{i=1}^{K}$. The group of permutations, $\mathbf{G}_{1}$ say, has $K!$ elements and the group of sign changes, $\mathbf{G}_{2}$ say, has $2^{K}$ elements (two possible signs for each element). The groups $\mathbf{G}_{1}$ and $\mathbf{G}_{2}$ have only the identity element in common, but are otherwise non-overlapping. The full group $\mathbf{G=G}_{1}\mathbf{\times G}_{2}$ generated by $\mathbf{G}_{1}$ and $\mathbf{G}_{2}$ therefore has $K!2^{K}$ elements. Exploiting the invariance and symmetry properties of the problem reduces the domain of integration by a factor $K!2^{K}$. This is important because all optimizations require probabilities calculated by numerical integration. The density after a sign change in $T_{k}$ is obtained by a corresponding sign change in $\mu_{k}$ and for a permutation of $T$ also $\mu$ permutes accordingly. Hence for any element $g$$\in\mathbf{G}$ we have $\mathbf{g}\cdot T\sim$ $\emph{N}\left(\mathbf{g}\cdot\mu,I_{K}\right)$ or $P_{\mathbf{g}\mu}\left[\mathbf{g}T\in A\right]=P_{\mu}\left[T\in A\right]$ so the distribution is invariant; see Lehmann and Romano (2005).

The Wald and LR tests are functions of the maximal invariant and so will our new test. The density is required for probability calculations and optimizations for the new test. It is easily derived for arbitrary dimension using Equation (6) of Vaughan and Venables (1972):

The noncentral Chi-distribution $\chi\left(t,\mu\right)$ with one degree of freedom equals the folded normal distribution and if $T\sim N\left(\mu,1\right),$ then the density of $\left|T\right|$ can be written as $f_{|T|}(t,\mu)=\,\sqrt{2/\pi}\exp\{-\tfrac{1}{2}(t^{2}+\mu^{2})\}\cosh(\mu\ t),$ for $t\geq 0$. Substitution in Lemma 1 and simplifying gives the following result that is the basis for the numerical calculations that follow:

The ordered squared $t$-statistic could also be used as maximal invariant. Lemma 1 would then lead to a density in terms of noncentral Chi-squared distributions.

A general method for constructing near similar tests involves three generic steps:

1. 1.

   Define a flexible boundary $g$ for the critical region in the relevant sample space.

2. 2.

   Define a criterion function $Q\left(g\right)$ that penalizes the deviation of the NRP from the level $\alpha$ for a grid of parameter values under the null (and possibly restrictions on $g$ and other aspects deemed relevant).

3. 3.

   Systematically vary and determine $g$ such that it minimizes the criterion function and is therefore as close to similarity as possible in the metric defined by $Q$.

The relevant sample space is determined by the particular testing problem at hand and may have been reduced by sufficiency, invariance, or other principles, to dimension $k$, say. The boundary $g$ of the critical and acceptance region is then of dimension $\left(k-1\right)$, but may consist of disjoint parts if the critical and/or acceptance region are not simply connected in a topological sense. There are various possibilities to define $g$ flexibly, but we will use splines.

The criterion function may include aspects other than similarity, for instance smoothness and monotonicity of $g$, convexity of the critical or acceptance regions, or even rejection probabilities under alternatives. Consequently, Step 3 will generally be a constraint optimization problem. The systematic variation of $g$ is intended to be in line with the optimization routine used to minimize $Q$ as in, e.g. a Newton-Raphson-type procedure.

For mediation testing with $\alpha=0.05$ an exact test exists, but the varying $g$-method may not find it because (i) $g$ is not flexible enough (ii) $Q\left(g\right)$ includes criteria other than NRPs (iii) numerical difficulties determining the jumps (iv) and finally restrictions to purposely exclude objectionable CRs.

An explicit implementation of the varying $g$-method to the mediation problem is given next.

Having determined an upper bound to the power envelope, we can determine a $g$-boundary function with a power surface as close as possible to this upper bound. This optimal test is found using the algorithm given in Appendix D. This function is given in Appendix E and R-code is also provided there. For ease of implementation we give values of $g\left(t\right)$ in Table 1. Figure 7 shows the optimal $g$-boundary test for the mediation problem.

The optimal $CR_{g}$ includes a narrow region close to the $45^{\circ}$ line where both $t$-statistics are of the same magnitude. This is expedient for two reasons. First, because mediation requires both $\theta_{1}$ and $\theta_{2}$ to be non-zero. The best possibility of detecting this is along the $45^{\circ}$ line as illustrated by the power surface in Figure 9 showing highest power on the diagonal. The optimal $CR_{g}$ does exclude both $t$-statistics smaller than $0.1$, unlike the unappealing region of the exact test. Second, the near similarity condition requires additional critical region area in the left corner of the octant because NRPs are particularly low for small parameter values. The increased power is naturally linked to the increase in Type I error, but correct size of a test by definition merely requires that this is not larger than $5\%$. Nevertheless, size (NRP)/power trade-off exists as well as other compromises that can be assessed using critical region analysis. For instance, it may seem less intuitive that rejection is not monotonic in $t_{1}$ and $t_{2}$ since an increase in both $t_{1}$ and $t_{2}$ represents increased evidence against the null. The LR and Wald tests are monotonic in this sense, but lead to a reduction in power to nearly zero for small parameter values. No observed value $t$ of $T$ will ever lie on the horizontal or vertical axis. Any observed $t$ is therefore more likely given an alternative parameter value than a value under the null. It is therefore desirable to add area to the LR critical region even if this results in a non-convex critical region or acceptance region. One could cogitate about the very narrow region close to the diagonal and whether the acceptance should not continue along the $45^{\circ}$ line further than $0.1,$ until e.g. $1.217$ as in Figure 2(b), but the new $g$-boundary is the optimal solution to a well-defined problem.

The narrow region of the optimal $CR_{g}$ is a strict extension of the $CR_{LR},$ which itself is strictly larger than the Sobel (Wald) $CR_{W}$. Since the new test is constructed to satisfy the size condition we have the following:

The NRP as a function of the noncentrality parameter $\mu$ is shown in Figure 8. The difference from $5\%$ is less than $10^{-5}$ and so small that the scale had to be magnified greatly, to an extend that prevents comparison with the LR and Sobel tests in the same graph.

The power of the new $g$-test is very close to the power envelope (upper bound). The maximal difference is $0.00614$. This has important implications. First, the upper bound is tight as claimed earlier. Second, the upper bound of the power envelope and power surface of the $g$-test look almost identical when graphed. Figure 9 therefore shows only the power surface of the optimal $g$-test. Finally, the new $g$-test is optimal for all intents and purposes in a larger class of tests. It is optimal by construction within the class of near similar tests $\Gamma_{\alpha,\epsilon}^{\mathbb{M}_{0}}$, but given the closeness of its power surface to the (non)similar power envelope, there cannot exist any near similar test that has additional power more than $0.00614$, even if construction is based on a different method. More generally, nonsimilar tests can have $2\%$ points more power, but at the possible cost of odd rejection regions and low power for other parameter values that are not used in the construction of the test. The optimal $g$-test has good properties for all parameter values.

The power surface in Figure 9 shows only the first quadrant of the parameter space of $\left(\mu_{1},\mu_{2}\right).$ The other three quadrants follow by simple permutations and reflections of the parameters.

To illustrate the usage of the new test, we consider the mediation analysis in Alan et al. (2018) on the effect of elementary school teachers’ gender beliefs, either traditional or progressive, on student mathematical and verbal achievements. They exploit the unique institutional features of Turkish data that provides a natural experiment in the random allocation of teachers to schools. Information consists of approximately 4,000 third- and fourth-grade students and their 145 teachers. The data are available in their online appendix. Children are divided into three groups depending on the length of exposure to a participating teacher: “1-year exposure” (at most one year), “2-3 year exposure” (more than one year and at most three years) and “4-year exposure” (at most four years). Alan et al. (2018) consider three potential mediators, but we focus on students’ own gender role beliefs. In the notation of equations (1)-(3), $X$ is a dummy whether the teacher is classified as traditional or progressive. The mediating variable $M$ is the student’s own belief on gender roles. We focus on the verbal test scores as the dependent variable $Y$. Table 2 shows the estimates for the indirect effect for girls based on the full sample and the three exposure groups after controlling for school fixed effects, student characteristics, family characteristics, teacher characteristics, teacher styles and teacher effort.

The results for the full sample are similar to the values shown in Table 5 and Table 6 of Alan et al. (2018), although they use the approach by Imai et al. (2013). We consider three different exposure duration groups. The $t$-ratios of $\hat{\theta}_{1}$ are not significant at the $5\%$ level for more than 1-year exposure. For the 1-year exposure, however, it is significant, but the $t$-ratio $t_{2}=-1.941$ of $\hat{\theta}_{2}$ is not. So, the LR test would not find a significant effect and neither does the simulation based method Alan et al. (2018) use. Using the new test, however, we have $\left|t\right|_{(1)}=1.941$ and $|t|_{(2)}=2.052$, such that $g(2.05)=1.9175$ and consequently $|t|_{(1)}>g(|t|_{(2)})$ and the proposed test finds the mediation effect significant at the $5\%$ level. The R function in Appendix E can be used if greater precision is desired, e.g. $g(2.052)=1.9195$.