Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
84,460 characters · 13 sections · 55 citation commands
A Nearly Similar Powerful Test for Mediation
\centerline{Version October, 2021}
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, Baron1986 has more than 100,000 citations\footnote{ \ 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. Coletti2005, marketing, e.g. Mackenzie1986, sociology, e.g. Alwin1975 who used the expression indirect effect, a term commonly used in economics also.For a recent overview of mediation in economics see huber2020mediation and e.g. Heckman2015a, Heckman2015b on treatment effects and production technology, and imbens2020potential who extensively reviews connections between directed acyclic graphs (DAGs), potential outcomes, causal inference, instrumental variables, and mediation. frewen2013perceived 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.
Mackinnon2002 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 ((ref)) and ((ref)):
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 ((ref)) is a restricted version of ((ref)) 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 (ref).
Figure (ref) 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 ((ref)) - ((ref)), 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 Sobel1982 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 Mackinnon2002 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 VanGarderenVanGiersbergen2021) 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 (ref) 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 Perlman1999 coined the term \textquotedblleft Emperor's New Tests\textquotedblright\ for similar tests with undesirable properties like these. Insistence on similarity can render LR tests $\alpha $-inadmissible, cf. Lehmann2005, but Perlman1999 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 Perlman1999, and that is nearly similar and practically unbiased.
It is extremely simple to use:
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\vert t\right\vert _{(1)},\left\vert t\right\vert _{(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. Andrews1994, Andrews2006,Andrews2008, Elliott2015, Guggenberger2019 . It can be given a random critical value interpretation as in moreira2016critical 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 (ref) 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 (ref) shows the ease of implementation using an interesting application by Alan2018 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 Alan2018 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 (ref). 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 ((ref)) and ((ref)) 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 ((ref)), but is exogenous for $\theta _{2}$ in ((ref)) 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 ((ref)) and ((ref))\footnote{ This 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. HillierVanGarderenVanGiersbergen2021 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 Dufour2017 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 Lehmann2005.
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 Vaughan1972:
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\vert T\right\vert $ 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 (ref) 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 (ref) would then lead to a density in terms of noncentral Chi-squared distributions.
Standard mediation test statistics used in practice have distributions that depend on the parameter values under the null. The rejection probabilities are therefore not constant and the tests are biased with power dropping below the size of the test, especially in a neighborhood of the origin. We illustrate the issue for the classic Wald and LR tests.
The null hypothesis $\theta _{1}\theta _{2}=0$ defines a manifold that is almost everywhere continuously differentiable, with the exception of the origin which is a so-called \textquotedblleft double point" where the two restrictions $\theta _{1}=0$ and $\theta _{2}=0$, each defining a one-dimensional line, coincide. The widely used Sobel1982 test equals the square root of the Wald test. Glonek1993 derives the asymptotic distribution for the Wald test statistic:
As a consequence the asymptotic $5\%$ critical value for $W$ when both $ \theta _{1}=0$ and $\theta _{2}=0$ is $\frac{1}{4}\chi _{1}^{2}(0.95),$ but jumps to $\chi _{1}^{2}(0.95)$ i.e. the usual Chi-squared critical value for one restriction, for any other value. The discrete jump in the asymptotic distribution from the origin to any other fixed parameter is remarkable and shows explicitly that the distribution depends heavily on the parameter values under the null. This discontinuity in the asymptotic distribution\ and dependence on the parameter also invalidates bootstrap procedures and they are oversized. For an NRP of $5\%$ at the origin the critical value should be $0.96$ but this would lead to over-rejection for other values under the null and the test would be oversized (size $>30\%$) and invalid. One could consider drifting sequences of parameter values to investigate the behavior of the Wald statistic near the origin, but that would not solve the problem. The problematic behavior of the Wald test under the null with singularities is well documented by DrtonXiao2016 and drton2009likelihood. No satisfactory solution has been found in the preceding decades to salvage the Wald statistic, see e.g. Dufour2017 . This prompted our investigation and to propose an alternative solution.
The LR test was shown by vanGiersbergen2014\ to equal:
and rejects when both $H_{0}^{\theta _{1}}:\theta _{1}=0$ and $H_{0}^{\theta _{2}}:\theta _{2}=0$ are rejected by basic $t$-tests. In Mackinnon2002 this is referred to as the test for joint significance, but not identified as the LR test. The rejection probability for critical value $cv$ is:
by independence of $T_{1}$ and $T_{2}$. These rejection probabilities are monotonically increasing in the absolute values of $\theta _{1}$ and $\theta _{2}$. Correct size is therefore obtained by choosing the critical value of the test by letting $\theta _{1}\rightarrow \infty $ when $\theta _{2}=0,$ or $\theta _{2}\rightarrow \infty $ if $\theta _{1}=0,$ to guarantee that the rejection probability under the null is always smaller than or equal to the nominal size. The asymptotic $5\%$ critical value is therefore the usual $1.96$. The NRP will depend on the values of $\theta _{1}$ and $\theta _{2}$ and vary between the following two extremes:
where $z_{0.025}$ is the upper $2.5\%$ percentile of the standard normal distribution. For an NRP of $5\%$ at the origin $\left( \theta _{1},\theta _{2}\right) =\left( 0,0\right) ,$ the critical value should equal $ cv_{LR00}=1.217$. This leads to massive over-rejection if only one parameter is zero and the other much larger. The test with this critical value is oversized $($size $>20\%)$ and invalid.
The third classic test, the Lagrange Multiplier (LM) or score test, is even more problematic because its definition depends on parameter values under the null and there are three different versions depending on which $\theta $ , or both $\theta $'s are zero.
All these classic tests are functions of two $t$-statistics. Their distributions, as well as their NRPs, clearly depend on the parameter values under the null and the tests are not similar. A test is called similar on the boundary of $H_{0}$ if the probability of rejecting the null is constant for all parameter values on the boundary of $H_{0}$ and $H_{1}$. For mediation, this boundary equals $H_{0}$ itself and consists of the horizontal and vertical axes of the $\left( \theta _{1},\theta _{2}\right) $ space. None of the classic tests is similar and in a neighborhood of the origin the NRPs are close to zero. As a result the power in a neighborhood of the origin is also close to zero and far below the size of the test and the tests are biased since there are parameter values with probability of rejection under the alternative lower than under the null.
The behavior and construction of the classic test statistics is problematic. Given that no satisfactory adjustments of classic test statistics have been found, despite considerable efforts over recent decades, a different approach is required.
In order to derive an alternative test procedure we shift the focus from the test statistic to the critical region (CR). A critical region defines a test statistic of course, but choosing a class of tests, such as Wald, LR, or LM tests, restricts the shape of the critical region. For the same reason the tests focusing on improving the standard error of $\hat{ \theta}_{1}\hat{\theta}_{2}$ or $\left( \hat{\tau}^{\ast }-\hat{\tau}\right) $ analyzed in Mackinnon2002 restrict possible shapes of the critial region.
We construct a new test procedure by constructing the critical region directly by determining its shape in the two-dimensional sample space of the $t$-statistics used in the construction of the tests. We consider critical regions that are bounded by a function $g$ and reject when $\left\vert T\right\vert _{\left( 1\right) }>g(\left\vert T\right\vert _{\left( 2\right) }).$ We impose some weak regularity conditions. In particular we assume that $g$ is\ a c\`{a}dl\`{a}g function from $ \mathbb{R} _{0}^{+}$ (including 0) to $ \mathbb{R} _{0}^{+}.$ This will allow $g$ to have jumps, but limits the number of jumps to countably many. We also insist that $g$ is weakly increasing. This assures that a rejection (acceptance) for a realization $(\left\vert t\right\vert _{\left( 1\right) },\left\vert t\right\vert _{\left( 2\right) }) $ is not reversed when either $\left\vert t\right\vert _{\left( 1\right) } $ or $\left\vert t\right\vert _{\left( 2\right) }$ is increased (decreased). This reasoning also motivates a CR that is\ (topologically) simply connected. Denote the set of weakly increasing c\`{a}dl\`{a}g functions by $\mathbb{D}\left( \mathbb{R} _{0}^{+}, \mathbb{R} _{0}^{+}\right) $, as in the\ common Skorokhod space notation, but here with weak monotonicity.
Justification for only using $t$-statistics is by sufficiency and invariance. First, the MLE \ \ $\hat{\lambda}=\left( \hat{\tau},\hat{\theta} _{2},\hat{\sigma}_{11},\hat{\theta}_{1},\hat{\sigma}_{22}\right) ^{\prime }$ \ is a complete minimal sufficient statistic. The model constitutes a full exponential model since the dimensions of the minimal sufficient statistic and the parameter space are equal; see vanGarderen1997. Second, $ T_{1}$ and $T_{2}$ have distributions under the null that are independent of the nuisance parameters $\tau ,$ $\sigma _{11},$ and $\sigma _{22}.$ HillierVanGarderenVanGiersbergen2021 shows that, also in finite samples, $T=\left( T_{1},T_{2}\right) ^{\prime }$ is a maximal invariant under an appropriate group of transformations generalizing the scale invariance of the $t$-statistics. Theorem (ref) shows that as a consequence of the permutation and reflection invariance only $ 1/8^{th}$ of the two-dimensional sample space of $T$ needs to be considered and we define the critical region in the first octant (east to north-east). The other seven octants in $\mathbb{R}^{2}$ follow by symmetry. The test defined by $CR_{g}$ is indeed invariant to permutations, reflections, and scale transformations. The domain of $g\left( \cdot \right) $ can therefore be restricted to the non-negative real line and bounded by the $45^{\circ }$ line: $g\left( x\right) \leq x.$
We can put the definition of a similar test in terms of the boundary function $g\left( \cdot \right) $, noting that $H_{0}$ is itself the boundary of $H_{0}$ and $H_{1}$:
The boundary functions defined by the Wald (Sobel) and LR are not similar. Figures (ref) and (ref) show the boundary functions that define the critical region in terms of $\left( T_{1},T_{2}\right) $ for the Wald and LR test. We show the boundaries for two critical values: one such that for large $\left\vert \theta _{1}\right\vert $ or $\left\vert \theta _{2}\right\vert $ the NRP is $5\%$ asymptotically. This value is the usual $3.84$ for the Wald test and $1.96$ for the LR\ test. The second, smaller critical value is such that the NRP is $5\%$ when $\theta _{1}=\theta _{2}=0.$ This value is $0.96$ for the Wald test and $1.217$ for the LR test. The rejection probabilities are shown as a function of the noncentrality parameter $\mu _{2}=\theta _{2}/\sigma _{\theta _{2}}$ for given $\mu _{1}=0,$ such that $H_{0}$ holds. For the LR test with critical value $1.96$ the NRP goes down to $0.0025=0.05^{2}$ when $\theta _{1}=\theta _{2}=0.$ In the second case, with the smaller critical value $1.217$, the NRP is $5\%$ by construction when $\theta _{1}=\theta _{2}=0,$ but this is not a valid test since for other values the NRPs are much higher than the nominal size of $5\%$. The same situation will occur when constructing point optimal invariant tests. The Wald test is considerably worse than the LR test with lower NRP over a wider range of $\mu _{1}$, see figures (ref) and (ref).
The trinity of classic tests is clearly nonsimilar. The question is if we can do much better. Does there exist a similar test or is this a problem that is intrinsically unsolvable? The next theorem, and the main theoretical contribution, answers this question.
For common significance levels, including $1\%$, $5\%$, and $10\%$, the theorem proves that exact similar tests exist. The proof is given in Appendix (ref) and exploits the symmetries of the problem and the completeness of the normal distribution. The proof is constructive showing that the $g$ function must be a step function and is unique within the class of weakly increasing c\`{a}dl\`{a}g functions. For $\alpha =0.25$ the exact similar boundary is shown in Figure (ref).
The CR shows a clear parallel with Figure 2 of berger1989. The step function starts horizontal and defines a CR with objectionable features. In particular, the CR will include all $T$ such that $0<\left\vert T\right\vert _{\left( 2\right) }<\Phi ^{-1}\left( \frac{1}{2}+\frac{\alpha }{2}\right) $. For $\alpha =0.05$ this corresponds to rejection when both $t$-statistics are smaller than $0.0627$ in absolute value (but non-zero) \ (for $\alpha =0.10$ and $0.01$ smaller than $0.1257$ and $0.0125$ respectively). Such $t$ -statistics close to zero cannot be characterized as strong evidence against the null, in favor of both $\theta _{1}$ and $\theta _{2}$ being non-zero. Nevertheless, the test is size correct and renders the LR test $\alpha $ -inadmissible. So it is an \textquotedblleft Emperor's New Tests \textquotedblright\ in the terminology of Perlman1999 and does not provide a satisfactory solution to the problem of finding a better test. A second unattractive feature of this exact similar boundary is that for increasing values of the test statistics parallel and close to the diagonal, the decision alternates between rejection and acceptance, despite the fact that evidence against the null appears monotonically increasing.
Given the uniqueness of the similar test, we relax the strict similarity requirement and consider a class of near similar tests with NRPs that differ from $\alpha $ by no more than $\epsilon $. Within this class of so called $ \epsilon $-similar tests given in Definition (ref) below, we determine a test that avoids the objectionable properties of the exact similar test, but achieves good power properties. Within this class of near similar tests we determine the power envelope and determine the test that minimizes the distance between its power surface and this power envelope. This new test is therefore optimal in this sense. It is easy to implement using Table (ref) or using the R-code in Appendix (ref).
The first step in the composition of this optimal test is a new general method for constructing near similar tests.
A general method for constructing near similar tests involves three generic steps:
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.
Step 1 in the varying $g$-method is to determine the relevant sample space for the testing problem. The first dimensional reduction is to the MLE which is minimal sufficient and complete. The second reduction to $\left( T_{1},T_{2}\right) $ follows from location-scale type invariance shown in HillierVanGarderenVanGiersbergen2021. Permutation and reflection symmetries further reduce the relevant sample space to one octant according to Theorem (ref). The maximal invariant is the ordered absolute $t$-statistic $(\left\vert T\right\vert _{\left( 1\right) },\left\vert T\right\vert _{\left( 2\right) })=\left( \min \left( \left\vert T_{1}\right\vert ,\left\vert T_{2}\right\vert \right) ,\max \left( \left\vert T_{1}\right\vert ,\left\vert T_{2}\right\vert \right) \right) $ with density given in Lemma \ (ref) . For general $\mu $, the rejection probability (RP) for the $g$-method is given by:
This two-dimensional integral can be reduced to a one-dimensional integral by expressing the inner integral in terms of the CDF of the standard normal distribution $\Phi \left( \cdot \right) $:
The RP thus simplifies to a one-dimensional integral by integrating formula ( (ref)) for $b=g(t_{2})$ over $t_{2}\in \lbrack 0,\infty )$, which greatly improves efficiency and accuracy.
All NRPs in the paper were determined by numerical integration over $ t_{2}\in \lbrack 0,\mu _{2}+9]$ \ under the null with $\mu _{1}=0$. The numerical integration was performed in Julia, see Bezanson17julia, using the function quadgk that is based on adaptive Gauss-Kronrod quadrature. For the final $g$-function, the estimated upper bound on the absolute error for the calculated NRP is approximately $ 10^{-8}$.
The $g$-boundary is generally determined by an algorithm and Appendix (ref) shows the basic implementation of the varying $g$ -method using linear splines with $J+2$ knots, with the first and last knots fixed. In spite of its simplicity, it leads to big improvements even for small values of $J$. Figure (ref) illustrates the construction of the $g$-function for a fixed number of grid points $J=6$ and the resulting $CR_{g}$ in the sample space of $( \left\vert T\right\vert _{\left( 1\right) },\left\vert T\right\vert _{\left( 2\right) }) $, the East to North-East octant of $ \mathbb{R} ^{2}$.
Figure (ref) shows the NRPs of the test in comparison to the LR\ and Wald (Sobel) tests. There is a remarkable gain in the lowest NRP, and therefore local power, from $0.25\%$ to $4.999\%$. Already for $J=2$ there is a large improvement and after $J=8$ gains are small, and beyond $ J=16$ there was hardly any improvement.
Insistence on similarity can have negative consequences for the power in general, but not here. The power envelope with or without (near) similarity restriction are very close. The power surface of our optimal test is also close. Even the basic test with $J=16$ has very good power in comparison with the Sobel and LR tests and is uniformly better for all values of the noncentrality parameter $\mu $. In a neighborhood of the origin with $\mu =0, $ it is essentially $5\%$ points higher.
The RP $\pi _{g}\left( \mu _{1},\mu _{2}\right) $ defined in Equation ((ref)) is the NRP if $\mu _{1}$ and/or $\mu _{2}$ equal $0$ (the null hypothesis). When both are non-zero,\ $H_{0}$ is false and $\pi _{g}\left( \mu _{1},\mu _{2}\right) $ is the power of the test defined by $CR_{g}$. Figure (ref) illustrates the power in the $ 45^{\circ },$ $\mu _{1}=\mu _{2},$ direction. Power in other directions is also superior to the Wald (Sobel) and LR tests.
There is a straightforward explanation for the additional power. The Wald and LR\ test both reject far less than $5\%$ near the origin. The critical region can therefore be extended and the power increased without failing the size condition. In the origin the NRPs are close to $0\%$ for the LR\ and Wald (Sobel) tests. By extending the critical region we can therefore gain almost $5\%$ power without violating the size condition. The LR test has some attractive features including that rejection for a particular $\left( t_{1},t_{2}\right) $ implies rejection for larger values of $t_{1}$ and/or $ t_{2}$ also. This is intuitive since the evidence against the null is increasing. A disadvantage is however that one never rejects when either $ t_{1}$ or $t_{2}$ is smaller than $1.96$ and this causes extreme conservativeness that can be resolved by adding area to the critical region. It should also be noted that the relevant distribution of the order statistic $(\left\vert T\right\vert _{\left( 1\right) },\left\vert T\right\vert _{\left( 2\right) })$ is quite different from the distribution of the product of two standard normals for $\left( T_{1},T_{2}\right) $ or their absolute values.
Comparison to the Sobel and LR tests is limited because they both perform very poorly for small values of $\mu $. The absolute quality, or even near optimality, of the new $g$-test can only be assessed by comparing the power surface of the test to the power envelope, or a tight upper bound thereof, for a class of tests that satisfy appropriate invariance-, size and almost similarity restrictions. Since the exact similar invariant test is objectionable when it exists, we introduce a class $\Gamma _{\alpha ,\epsilon }$ of near similar tests with NRPs that deviate less than $ \epsilon $ from the $\alpha $ level and an operational (super)class $\Gamma _{\alpha ,\epsilon }^{\mathbb{M}_{0}}\supseteq \Gamma _{\alpha ,\epsilon }$ as follows.
For $\epsilon =0$ the boundary functions in $\Gamma _{0}$ would be similar. We consider $\alpha =0.05$ in all our numerical examples, but for general $ 1/\alpha \notin \mathbb{N} $ the set would be empty according to Theorem (ref). For $\epsilon =\alpha $, on the other hand, $\Gamma _{\alpha ,\alpha }$ contains all $g$-based tests that satisfy the size condition. Our interest is in $\epsilon $ close to $0$, when $\Gamma _{\alpha ,\epsilon }$ contains boundaries that are almost similar. The minimum value of $\epsilon $ for which $\Gamma _{\alpha ,\epsilon }$ is not empty is 0 for the mediation problem when $1/\alpha \in \mathbb{N} $, but in general depends on the testing problem considered and larger than $ 0$ if no similar test exists.
The class $\Gamma _{\alpha ,\epsilon }^{\mathbb{M}_{0}}$ can be thought of as a discretization of $\Gamma _{\alpha ,\epsilon }$ in the sense that a grid of points under the null is considered. It imposes less restrictions and enforces near similarity conditions on a finite number of points only. As a consequence it may contain boundaries that do not satisfy the size condition for points that are not in $\mathbb{M}_{0}.$ Obviously $\Gamma _{\alpha ,\epsilon }$\ $\subseteq \Gamma _{\alpha ,\epsilon }^{\mathbb{M} _{0}}$ since the size and NRP conditions also hold for the points in $ \mathbb{M}_{0}$.
Within the class $\Gamma _{\alpha ,\epsilon }$ there is no unique solution. As a consequence one has to choose a boundary function from $\Gamma _{\alpha ,\epsilon },$ or in practice from $\Gamma _{\alpha ,\epsilon }^{\mathbb{M} _{0}},$ to obtain an operational test. For the construction of the power envelope we can select the test in $\Gamma _{\alpha ,\epsilon }^{\mathbb{M} _{0}}$ that maximizes the power against a particular point $\left( \mu _{1},\mu _{2}\right) $ in the alternative. This test is a Point Optimal Invariant Near Similar (POINS) test. The critical region of this test varies with $\left( \mu _{1},\mu _{2}\right) $ and no uniformly most powerful test exists within the class $\Gamma _{\alpha ,\epsilon }$. It can be used however, to construct an upper bound for the power envelope.
For notational simplicity we have suppressed the dependence on $\epsilon $ and $\mathbb{M}_{0}.$ Since $\Gamma _{\alpha ,\epsilon }$\ $\subseteq \Gamma _{\alpha ,\epsilon }^{\mathbb{M}_{0}}$ and elements of $\Gamma _{\alpha ,\epsilon }^{\mathbb{M}_{0}}$ do not necessarily satisfy the size condition for all parameter values it follows that $\bar{\pi}\left( \mu _{1},\mu _{2}\right) \geq \pi \left( \mu _{1},\mu _{2}\right) ,$ because fewer conditions are imposed.\ Choosing a finer grid for $\mathbb{M}_{0}$ will force $\bar{\pi}\left( \mu _{1},\mu _{2}\right) \ $closer to $\pi \left( \mu _{1},\mu _{2}\right) $, at least in the additional points in $\mathbb{M}_{0}$ where the size condition is now required to hold. Also note that the \textquotedblleft point\textquotedblright\ optimal $g$ that maximizes power for the point $\left( \mu _{1},\mu _{2}\right) ,$ may have undesirable features such as including parts of the axes in the critical region, even though such observations are perfectly in line with the null hypothesis.
We determine $\bar{\pi}\left( \mu _{1},\mu _{2}\right) $ numerically for $ \alpha =0.05$ by maximizing the power directly by selecting critical region points in the sample space that maximize the probability of rejection when the true density has parameter $\left( \mu _{1},\mu _{2}\right) $, under the side conditions that the NRP$\in \left[ 0.05-\epsilon ,0.05\right] $ for all parameters $\left( 0,\mu _{0}\right) \in \mathbb{M}_{0}$. The sample space is decomposed into $285,150$ squares\ and a modern optimization routine is used to determine which squares should be included in the critical or acceptance region in order to maximize the power while at the same time satisfying the approximate similarity condition. This is repeated for a grid of $\left( \mu _{1},\mu _{2}\right) $ points. So for each point on the grid the POINS\ critical region is determined and the power recorded. Appendix (ref) gives details of the algorithm and the optimization routine that can deal with a large number of variables and side conditions.
By dropping the near similarity restriction ($0.05-\epsilon \leq NRP)$ in the same algorithm, we can construct a power envelope for nonsimilar tests. The maximal difference from the (higher) nonsimilar power surface is $2\%$ points when power is around $40\%,$ showing that the power loss due to the similarity requirement is small. The calculated $\bar{\pi}\left( \mu _{1},\mu _{2}\right) $ surface enables us to construct a correctly sized optimal test derived in the next section.
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 (ref). This function is given in Appendix (ref) and R-code is also provided there. For ease of implementation we give values of $g\left( t\right) $ in Table (ref). Figure (ref) 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 (ref) 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 (ref), 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 (ref). 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 (ref) 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 (ref) 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 Alan2018 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: \textquotedblleft 1-year exposure\textquotedblright\ (at most one year), \textquotedblleft 2-3 year exposure\textquotedblright\ (more than one year and at most three years) and \textquotedblleft 4-year exposure\textquotedblright\ (at most four years).\ Alan2018 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 (ref) 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 Alan2018, although they use the approach by imai2013experimental. 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 Alan2018 use. Using the new test, however, we have $\left\vert t\right\vert _{(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$.
In a further empirical illustration below, mediation may be via channels that involve two mediating variables. This requires a multivariate extension of the new test. The necessary invariance and distribution theory was given in Section (ref) but a further aspect is the coherency between solutions in different dimensions. Consider the null hypothesis $ H_{0}:\theta _{1}\theta _{2}\cdots \theta _{K}=0$ in $K$ dimensions. If it were known that $\theta _{K}\neq 0,$ then the null hypothesis reduces to $ H_{0}:\theta _{1}\theta _{2}\cdots \theta _{K-1}=0.$ This implies that the critical region for the $K-1$ corresponding $t$-statistics should reduce to the solution found for $K-1$ dimensions when $\left\vert \mu _{K}\right\vert $ is large. For large values of $\left\vert T_{K}\right\vert $ (very small $ p $-values) it is essentially known that $\mu _{K}$ and $\theta _{K}$ are non-zero. The probability of rejection will effectively depend only on the $ K-1$ other $t$-values. In two dimensions this means that as $t\rightarrow \infty $ the boundary function $g\left( t\right) \rightarrow 1.96$, which is the one-dimensional solution for testing $H_{0}:\theta _{1}=0$ when $\alpha =5\%.$ In three dimensions it means that the solution must reduce to the $g$-test derived in Section (ref). We make this requirement explicit in the definition below.
We used a multivariate spline generalization to implement the varying $g$ -method based on barycentric coordinates in three dimensions. It resulted in a maximum of $0.2\%$ points difference from $5\%$. But imposing dimensional coherency was complicated and the problem suffers from the curse of dimensionality: the dimension of the integral increases with $K$ and the number of knots necessary to define $g$ impedes optimization. For practical purposes we therefore propose a simple solution that exploits the dimensional coherency inductively and weighs the LR test in dimension $K$ with the solution obtained in dimension $K-1$. First note that one could satisfy the coherency condition in three dimensions by simply rejecting when $\left\vert T\right\vert _{\left( 1\right) }>g(\left\vert T\right\vert _{\left( 2\right) }),$ irrespective of $\left\vert T\right\vert _{\left( 3\right) }.$ This results in an invalid test however, because it is oversized with a maximum NRP of $7.2\%$ when $K=3$. The LR test on the other hand is conservative, in particular near the origin. A practical solution therefore is to use a weighted average between the liberal and conservative boundary. We use weights that depend on the largest absolute $t$-statistic. In particular for $K=3$
with weight function $w(t_{3})$ a linear spline with $0\leq w(t_{3})\leq 1,$ and $\lim_{t_{3}\rightarrow \infty }w(t_{3})=1$. The test rejects if $ \left\vert T\right\vert _{\left( 1\right) }>g(\left\vert T\right\vert _{\left( 2\right) },\left\vert T\right\vert _{\left( 3\right) }).$ Minimizing deviation ofthe NRP from the significance level and imposing the size restriction, results in a spline with knots:
leading to a maximum of $0.13\%$ points difference in NRPs from $5\%$ and never exceeding $5\%$. The solution is shown in Figure (ref).
For a numerical illustration, we consider the recursive model of union sentiment among southern nonunion textile workers as used by Bollen1990. The model:
is a simplified version of McDonald1984 and discussed in some detail by Bollen1989. It analyses the direct and indirect effects of tenure and age on union sentiment via deference and/or labor activism. Tenure $x_{1}$ is measured in log of years working in a particular textile mill and age $x_{2}$ is measured in years. The variables sentiment towards unions $y$, deference/submissiveness to managers $m_{1}$, and support for labor activism $m_{2}$, are measures based on 7, 4, and 9 survey questions respectively. The disturbances ($u$, $v_{1}$ and $v_{2}$) are assumed to be uncorrelated across equations and individuals. When they are normally distributed, ML estimation of the system reduces to OLS applied to each equation separately due to the recursive structure.
We use a selection of 100 observations out of the original 173 and focus on three alternative theories of the indirect effects from age to union sentiment: two competing parallel effects that the age effect is mediated by increased deference in which case $i_{1}=\alpha _{32}\beta _{13}$ quantifies the indirect effect. The alternative mediation channel is that activism mediates such that $i_{2}=\alpha _{22}\beta _{12}$ is the indirect effect. The third channel is a serial effect that age affects deference, which in turn affects activism, which in turn affects union sentiment such that $ i_{3}=\alpha _{32}\beta _{12}\beta _{23}$ measures the indirect effect. Figure (ref) illustrates the three mediation channels. The OLS estimates of the coefficients of the structural equations and their $ t$-statistics are shown in Table (ref).
The point estimates of the indirect effects and their $t$-statistics based on the delta method are shown in Table (ref). For the $g$-test we need the absolute order statistics, evaluate $g,$ and compare. For $H_{0}^{i_{1}}:\alpha _{32}\beta _{13}=0$\ we observe $|t(\hat{ \beta}_{13})|=1.838>1.770=g(1.902)=g(|t(\hat{\alpha}_{32})|),$ and hence reject. For $H_{0}^{i_{2}}:\alpha _{22}\beta _{12}=0$ \ we have $|t(\hat{ \alpha}_{22})|=2.709>1.960=g(7.120)=g(|t\left( \hat{\beta}_{12}\right) |)$ and also reject. Testing the last null hypothesis $H_{0}^{i_{3}}:\alpha _{32}\beta _{12}\beta _{23}=0$ requires the three-dimensional solution given in Figure (ref). We have $|t(\hat{\alpha} _{32})|=1.902<1.96=g_{2}(3.582,7.120)=g_{2}(|t\left( \hat{\beta} _{23}|\right) ,|t\left( \hat{\beta}_{12}\right) |)$ and do not reject.
The Sobel test with critical value $1.96$ concludes that $i_{2}$ is significant but does not find enough evidence for the $i_{1}$ mediation channel. The new $g$-test in contrast, concludes that $i_{1}$ is also significant. Both $t$-values in this case are smaller than $1.96$, so the LR test would not reject either. The two $t$-values are of comparable magnitude and the $g$-test finds a significant mediation effect. For implementation of the $g$-test only the relevant $t$-statistics are required. The absolute values are ordered and the smallest value compared with the value of the $g$ -function evaluated at the largest absolute $t$-value. This can be looked up in Table (ref) (possibly using linear interpolation) or one can use the R code provided in Appendix (ref). \footnote{ The bootstrap is a popular alternative for testing mediation. Because of the asymmetry of the distribution involved this is carried out through alternative confidence intervals of the indirect effect. See e.g. MacKinnon2004 and Preacher2008. The bootstrap is not valid. Simulations we carried out showed that bootstrap tests for mediation based on generally preferred BCa confidence intervals can have sizes of 8% when $ n=100$ and higher for $n$ smaller.}
For $i_{3}$ both tests draw the same conclusion. The three $t$-values involved are not of comparable magnitude and the $t$-statistics for $\beta _{12}$ and $\beta _{23}$ are so large that rejecting the null $H_{0}:\alpha _{32}\beta _{12}\beta _{23}=0$ essentially depends on whether $\alpha _{32}$ is zero.\ The corresponding absolute $t$-value of $1.90$ is too small to warrant such conclusion.
This paper proposes a new near similar more powerful mediation test that is simple to used based on two ordinary $t$-statistics.\ The mediation problem is empirically extremely important in many different fields. Theoretically we solve an interesting statistical problem which has generated results dating back to Craig1936 and still continues today with contributions on poor performance of the Wald statistic, construction of similar tests, and hypotheses with singularities. A main theoretical contribution has been the derivation of an exact similar test that is unique within the general class considered (with CR that is topologically simply connected and has weak monotone c\`{a}dl\'{a}g boundary). This exact test has unattractive statistical properties, however, leading us to consider a class of nearly similar tests and choose an attractive test within it. By relaxing the strict similarity condition we are able to construct a near similar test that has superior power and avoids the disagreeable properties of the exact similar test and would please even statistically erudite emperors in Perlman1999. The new test can also be justified asymptotically under much weaker conditions and other estimation methods.
The new test is constructed using a new general method we propose for constructing tests that are approximately similar. This varying-$g$ method considers a flexible critical region boundary and minimizes the difference from the level $\alpha $ of the rejection probabilities at a number of points on the boundary of the null hypothesis. Conceptually and practically this was very simple and straightforward to implement. It does not require, as in other approaches, a choice of mixture distribution, nor the construction of least favorable distributions. Numerically it is also attractive in terms of convergence properties and avoids the need for simulations. The appropriate distribution theory for the mediation case and higher dimensional extensions is explicitly given. All our calculations are done using numerical integration using the distribution of the maximal invariant.
The new method is applicable to many other testing problems with nuisance parameters. It is remarkable that this simple method works so well and can deliver substantial improvements. Even the simplest linear interpolation implementation for the mediation hypothesis increases power by almost $5\%$ points when mediation effects are small.
We have calculated a power envelope upper bound for the mediation testing problem that is very tight. Using this result, we were able to construct a test that is optimal within the class of near similar tests. It minimizes the total difference between its power surface and the power envelope bound. It results in a point wise difference less than $0.0062$ for all alternative parameter points considered. This implies that the test is practically optimal even if a more general class of possible test construction is considered. A power envelope for nonsimilar tests showed that power loss due to the similarity requirement is minimal since the maximum power loss is less than $2\%$ points (when maximum power is around $40\%$).
The optimal $g$-test satisfies the size condition. The critical region is strictly larger than the LR and Wald critical regions and is therefore strictly and uniformly more powerful. The classic tests are therefore not admissible and their bootstrapped variants are not valid. For large values of the standardized coefficients the power difference becomes negligible, but when mediation effects are small or have relatively large standard errors, the power can be close to $5\%$ points higher than these classic $ 5\% $-level tests. This has important consequences for empirical work. It enables researcher to prove mediation effects earlier in circumstances that one could not show mediation before due to extreme conservativeness of standard tests near the origin.