Mediation analysis aims to assess if, and how, a certain exposure influences an outcome of interest through intermediate variables. This problem has recently gained a surge of attention due to the tremendous need for such analyses in scientific fields. Testing for the mediation effect is greatly challenged by the fact that the underlying null hypothesis (i.e. the absence of mediation effects) is composite. Most existing mediation tests are overly conservative and thus underpowered.
To overcome this significant methodological hurdle, we develop an adaptive bootstrap testing framework that can accommodate different types of composite null hypotheses in the mediation pathway analysis. Applied to the product of coefficients (PoC) test and the joint significance (JS)  test, our adaptive testing procedures provide type 1 error control under the composite null, resulting in much improved statistical power compared to existing tests. Both theoretical properties and numerical examples of the proposed methodology are discussed.

§ INTRODUCTION

Mediation analysis plays a crucial role in investigating
the underlying mechanism or pathway between  an exposure and an outcome through an intermediate variable called a  mediator  [MacKinnon, 2008, VanderWeele, 2015].
It decomposes the “total effect” of an exposure on an outcome into an indirect effect that is through a given mediator and a direct effect, not through the mediator.
The former holds the key to uncovering the exposure-outcome mechanism and is often known as the mediation effect.
The mediation effect was initially studied under structural equation models (SEMs) in social sciences [Sobel, 1982, Baron and Kenny, 1986] and has  been given formal causal definitions  [Robins and Greenland, 1992, Pearl, 2001, Imai et al., 2010] within the counterfactual framework  [Imbens and Rubin, 2015].
the presence or absence of the mediation effect
can facilitate a deeper understanding of the underlying causal pathway from the exposure to the outcome and can give essential  insights into
intervention consequences, e.g.,   manipulating the mediator to change the exposure-outcome mechanism.
As a result, it is of interest to apply  mediation analysis
in many scientific fields,
such as psychology [MacKinnon and Fairchild, 2009, Valeri and VanderWeele, 2013],
genomics [Zhao et al., 2014, Huang and Pan, 2016, Huang, 2018, Guo et al., 2022], and
epidemiology [Barfield et al., 2017, Fulcher et al., 2019], among others.

To analyze the mediation effect, one classical setting  models the relationship between the exposure, the potential mediator, and the outcome as a directed acyclic graph; see Figure <ref>.

[shape=circle,draw, thick] (S) at (-2.7,0) $S$;
[shape=circle,draw, thick] (M) at (0,0) $M$;
[shape=circle,draw, thick] (Y) at (2.7,0) $Y$;
[shape=circle,draw, thick] (X) at (0, -1.6) $X$;

[->, thick, red] (S.east) – node[midway,above] $\alpha_S$ (M.west);
[->, thick, red] (M.east) – node[midway,above] $\beta_M$ (Y.west);
[->, thick] (X.north) –  (M.south);
[->, thick] (X) –  (Y);
[->, thick] (X) –  (S);
[->, thick] (S) .. controls +(up:15.5mm) and + (up:14.5mm)  .. node[midway,above] $\tau_S$  (Y);

Directed Acyclic Graph for Mediation Analysis. The exposure is $\S$; the mediator is $\M$; the outcome is $\Y$; the potential confounders are $X$.

let $\alphaS$ parametrize the causal effect of the exposure on the mediator, and $\betaM$ parametrize the causal effect of the mediator on the outcome conditioning on the exposure.
Then in the classical linear SEM   [Sobel, 1982, Baron and Kenny, 1986],
the causal mediation effect is proportional to  $\alphaS \betaM$
under suitable identification assumptions  [Imai et al., 2010].
More generally,
this product expression $\alphaS \betaM$ may also appear in the  causal mediation effect  under many other models,
such as generalized linear models and survival analysis models [VanderWeele and Vansteelandt, 2010, VanderWeele, 2011, Huang and Cai, 2016].
the important scientific question of whether or not the mediation effect is absent can be formulated as the hypothesis testing problem $H_0: \alphaS \betaM = 0$ against $H_A: \alphaS\neq 0$ and $\betaM \neq 0$ [MacKinnon, 2008].
Note that $H_0: \alphaS\betaM=0$ holds if and only if $\alphaS=0$ or $\betaM = 0$,
corresponding to two parameter sets $\mathcal{P}_{\alpha}=\{(\alphaS,\betaM) : \alphaS = 0\}$
and $\mathcal{P}_{\beta}=\{(\alphaS,\betaM) : \betaM = 0\}$, respectively.
It follows that the parameter set of $H_0:\alphaS\betaM=0$ is the union of two sets  $\mathcal{P}_{\alpha}$ and $\mathcal{P}_{\beta}$.
We visualize  $\mathcal{P}_{\alpha}$, $\mathcal{P}_{\beta}$, and their union  $\mathcal{P}_{\alpha}\cup \mathcal{P}_{\beta}$
in Figures <ref>–<ref>, respectively.

minor tick num=0,
axis line style = line width=1pt, draw=gray!45,
grid style=line width=.1pt, draw=gray!45,
axis lines = middle,
major tick style=draw=none,
xlabel=$\alphaS$,ylabel=$\betaM$,label style =
at=(ticklabel cs:1.18), font=]
[-stealth, green!60!black, line width=1.7pt] coordinates (0,-1) (0,1);
[-stealth, line width=0.5pt, gray] coordinates  (0.24,0.44) (0.02,0.44);
[-, line width=.5pt, draw=gray!45] coordinates  (-1,-1) (-1,1);
[-, line width=.5pt, draw=gray!45] coordinates  (-1,-1) (1,-1);
[-, line width=.5pt, draw=gray!45] coordinates  (1,-1) (1,1);
[-, line width=.5pt, draw=gray!45] coordinates  (1,1) (-1,1);

[green!60!black, align=center, thick] at (axis cs: 0.61,0.46) $\alphaS=0$;

minor tick num=0,
axis line style =  line width=1pt, draw=gray!45,
grid style=line width=.1pt, draw=gray!45,
axis lines = middle,
minor tick length=0.1pt,
major tick style=draw=none,
xlabel=$\alphaS$,ylabel=$\betaM$,label style =
at=(ticklabel cs:1.18), font=]
[-stealth, blue!85!black, line width=1.6pt] coordinates
(-1,0) (1,0);
[-stealth, line width=0.5pt, gray] coordinates  (0.58,0.34) (0.58,0.05);
[-, line width=.5pt, draw=gray!45] coordinates  (-1,-1) (-1,1);
[-, line width=.5pt, draw=gray!45] coordinates  (-1,-1) (1,-1);
[-, line width=.5pt, draw=gray!45] coordinates  (1,-1) (1,1);
[-, line width=.5pt, draw=gray!45] coordinates  (1,1) (-1,1);

[blue!85!black, align=center, thick] at (axis cs:  0.58,0.46) $\betaM=0$;

minor tick num=0,
line width=1pt, draw=gray!45,
grid style=line width=.1pt, draw=gray!45,
axis lines = middle,
major tick style=draw=none,
xlabel=$\alphaS$,ylabel=$\betaM$,label style =
at=(ticklabel cs:1.18), font=]
[-, line width=.5pt, draw=gray!45] coordinates  (-1,-1) (-1,1);
[-, line width=.5pt, draw=gray!45] coordinates  (-1,-1) (1,-1);
[-, line width=.5pt, draw=gray!45] coordinates  (1,-1) (1,1);
[-, line width=.5pt, draw=gray!45] coordinates  (1,1) (-1,1);

[-stealth, red, line width=1.5pt] coordinates (-1,0) (1,0);
[-stealth, line width=0.5pt, gray] coordinates  (0.6,0.45) (0.6,0.05);
[red, align=center, thick] at (axis cs:  0.52,0.7) $\alphaS\betaM=0$;
[-stealth, line width=0.5pt, gray] coordinates (0.38,0.44) (0.02,0.3);
[-stealth, red, line width=1.5pt] coordinates (0,-1) (0,1);

minor tick num=0,
axis line style =  line width=1pt, draw=gray!45,
grid style=line width=.1pt, draw=gray!45,
axis lines = middle,
major tick style=draw=none,
xlabel=$\alphaS$,ylabel=$\betaM$,label style =
at=(ticklabel cs:1.18), font=]
[-, line width=.5pt, draw=gray!45] coordinates  (-1,-1) (-1,1);
[-, line width=.5pt, draw=gray!45] coordinates  (-1,-1) (1,-1);
[-, line width=.5pt, draw=gray!45] coordinates  (1,-1) (1,1);
[-, line width=.5pt, draw=gray!45] coordinates  (1,1) (-1,1);

[-stealth, line width=0.5pt, gray] coordinates  (0.45,0.37) (0.11,0.11);
[orange, align=center, thick] at (axis cs: 0.5,0.79) $\alphaS=0$;
[orange, align=center, thick] at (axis cs: 0.5,0.53) $\betaM=0$;

[regular polygon,regular polygon sides=3, orange, draw, fill=orange, inner sep=1.7pt] at (axis cs: 0,0) ;

Visualization of parameter spaces of $(\alphaS, \betaM)$ under different constraints

To test $H_0:\alphaS\betaM=0$,
a broad class of methods
is  based on the  product of coefficients (PoC)
$\hatalphaSn \hatbetaMn$, where $\hatalphaSn$ and $\hatbetaMn$ denote the sample estimates
of parameters $\alphaS$ and $\betaM$, respectively.
One popular  PoC method
is Sobel's test  [Sobel, 1982], which is a Wald-type test and approximates the variance of $\hatalphaSn \hatbetaMn$ by the first-order Delta method.
In addition,  the joint significance (JS) test [Fritz and MacKinnon, 2007], also known as the MaxP test, is another widely used test that rejects the $H_0$ of no mediation effect
if both $\hatalphaSn$ and $\hatbetaMn$ pass a certain cutoff of statistical significance. [Liu et al., 2021] pointed out that
the MaxP test is a kind of likelihood ratio test under  normality assumptions.

Although there are
various procedures available for testing mediation effects,
properly controlling the type 1 error remains a challenge due to the  intrinsic structure of the null parameter space.
In particular,  $H_0: \alphaS\betaM=0$ is  composed of three different parameter cases:
(i) $\alphaS = 0$ and $ \betaM \neq 0$;
(ii) $\alphaS \neq 0$ and $\betaM = 0$;
(iii) $\alphaS = 0$ and $\betaM = 0$.
Case (iii), illustrated in Figure <ref>,
is a singleton given by the intersection set $\mathcal{P}_{\alpha}\cap \mathcal{P}_{\beta}$.
Under case (iii), both parameters $\alphaS$ and $\betaM$ are fixed at $0$,
whereas cases (i) and (ii) have one fixed parameter and the other parameter to be estimated.
This intrinsic difference leads to distinct asymptotic behaviors of test statistics.
Since the underlying truth is typically unknown in practice, it is difficult to obtain  proper $p$-values under the composite null hypothesis.

Particularly, in the popular Sobel's test and the MaxP test, the asymptotic distributions of the test statistics under cases (i) and (ii) are known to be different from those under case (iii). These tests have been shown to be overly conservative in case (iii), because statistical inference is carried out according to the asymptotic distributions determined in cases (i) and (ii) [MacKinnon et al., 2002, Fritz and MacKinnon, 2007].
This issue has gained a surge of attention
in recent genome-wide  epidemiological studies,
where for the majority of omics markers,
it holds that $\alphaS =\betaM = 0$, and

the classical tests are generally underpowered [Barfield et al., 2017].
Some recent work    [Liu et al., 2021, Dai et al., 2020, Du et al., 2022]  utilized the relative proportions of the cases (i)–(iii) in the population, but they rely on accurate estimation of the true proportions.
[Huang, 2019, Huang, 2019]  adjusted the  composition of $H_0:\alphaS\betaM=0$ through the variances of test statistics but required that the non-zero
coefficients are weak and sparse, which can be violated when the sample size is large.
Another line of related research [Sampson et al., 2018, Djordjilović et al., 2019, Djordjilović et al., 2020, Derkach et al., 2020] used a screening step
to control the family-wise error rate or the false discovery rate for a large group of hypotheses,
but they did not directly provide proper  $p$-values for each of the  composite null hypotheses.
[Van Garderen and Van Giersbergen, 2022] proposed to  construct a critical region for testing that can nearly    control the type 1 error at one prespecified significance level.
[Miles and Chambaz, 2021]  construct a rejection region that can achieve type 1 error control at significance level $\omega$ with $\omega^{-1}$ a positive integer.
Despite these developments, the fundamental issue of correctly characterizing the distributions of test statistics  to obtain
well-calibrated $p$-values under a composite null hypothesis remains an important challenging  problem in the current literature of mediation analyses.

In this paper, we develop a new hypothesis testing methodology to  address the challenge of obtaining uniformly distributed $p$-values under the composite null hypothesis
of no mediation effect.
Particularly, we propose an adaptive bootstrap procedure that can flexibly accommodate different types of null hypotheses.
In the current literature,
the nonparametric bootstrap is   directly applied to the PoC test statistic  $\hatalphaSn\hatbetaMn$,
which has been, unfortunately, found numerically to be
overly conservative when $\alphaS=\betaM=0$  [Fritz and MacKinnon, 2007, Barfield et al., 2017].
This paper unveils analytically the reason for the failure of the conventional  nonparametric bootstrap method, which stems from non-regular limiting behaviors
of the PoC test statistic at the neighborhood of  $(\alphaS, \betaM)=(0,0)$.
To overcome the non-regularity near $(\alphaS, \betaM)=(0,0)$,
we derive an explicit representation of the asymptotic distribution of the PoC test statistic through a local model,
and perform a consistent bootstrap estimation by incorporating suitable thresholds.
In addition, for the JS test,
we show that the conventional nonparametric bootstrap also fails to control type I error properly, which can be fixed by an adaptive bootstrap test similar to the procedure of the PoC test.
For both the PoC test and the JS test,
the proposed methods can  circumvent the nonstandard  limiting behaviors of the test statistics
and therefore uniformly adapt to different types of null cases  of no mediation effect.

The structure of this paper is as follows. In Section  <ref>,
we briefly review the basic problem setting and several popular testing methods in the literature.
In Section <ref>, we introduce the adaptive bootstrap method that can be applied to the representative PoC and JS tests under classical linear SEMs.
In Section <ref>, we conduct extensive simulation studies to compare the finite sample performances of the proposed tests with  popular counterparts.
In Section <ref>, we apply our adaptive bootstrap tests to investigate the mediation pathways of metabolites on the association of  environmental exposures with a health outcome.
In Section <ref>, we develop extensions of the  adaptive bootstrap, including joint testing of multivariate mediators and testing mediation effects under non-linear models.
We conclude the paper and discuss interesting extensions in Section <ref>.

Notation. For two sequences of real numbers $\{a_n\}$ and $\{b_n\}$,
we write $a_n=o(b_n)$ if $\lim_{n\to \infty} a_n/b_n=0$. We let $\xrightarrow{d}$ denote convergence in distribution.
We let $\overset{d^*}{\leadsto} $ denote   bootstrap consistency relative to the Kolmogorov-Smirnov distance;
see an introduction of this consistency notion in Section 23 of [van der Vaart, 2000]. To ensure clarity, we also provide the definitions of all the convergence modes in Section <ref> of the Supplementary  Material.

§ HYPOTHESIS TESTS OF NO MEDIATION EFFECT

Mediation Analysis Model.   To examine the mediation effect of the exposure $\S$ on the outcome $\Y$ through the intermediate variable $\M$, the causal inference literature  utilizes the counterfactual framework  [VanderWeele, 2015].
In particular, let $M(s)$ denote the potential value of the mediator under the exposure $\S=s$, and let $\Y(s, m)$ denote the potential outcome that would have been observed if $\S$ and $\M$ had been set to $s$ and $m$, respectively.
Throughout the paper, we adopt the  Stable Unit Treatment Value Assumption [Rubin, 1980], so that $\M=\M (\S)$ and $Y=Y (\S,\M(S))$.
Then the mediation effect or the natural indirect effect of $\S=s$ versus $s^*$ [Imai et al., 2010] is defined as

\begin{align}
\mathrm{E}\big\{ \Y(s, \M(s)) - \Y(s, \M(s^*))  \big\}. \label{eq:MEpotential}
\end{align}

For  ease of illustration, we consider the popular linear Structural Equation Model (SEM)  [MacKinnon, 2008, VanderWeele, 2015]:

\begin{eqnarray}
\M &=& \alphaS \S +  \X^\mytrans \alphaX +\eM, \label{eq:fullmod1} \\  \Y   &=&  \betaM \M + \X^\mytrans  \betaX + \DE \S + \eY,   \notag
\end{eqnarray}

where $\X$ denotes a vector of confounders with the first element being 1 for the intercept,
and $\eY$ and $\eM$ are independent error terms with mean zero and finite variances $\sigmaeY^2$ and $\sigmaeM^2$, respectively.
We assume that there are $n$ independent and identically distributed (i.i.d.) observations, $\{(\S_i, \X_i, \M_i, \Y_i): i=1,\ldots, n\}$,  sampled from Model (<ref>).
The independence of $\eY$ and $\eM$  holds under the following no unmeasured confounding assumptions.
In particular,
let $A\perp  B \mid C$ denote the independence of $A$ and $B$ conditional on $C$,
and we assume that for all levels of $s, s^*$ and $m$,
(i) $\Y(s,m) \perp \S \mid \{\X =\boldsymbol{x}\}$, no confounder for the relation of $\Y$ and $\S$;  (ii) $\Y(s,m) \perp \M \mid \{\S=s, \X=\boldsymbol{x}\}$, no confounder for the relation of $\Y$ and $\M$ conditioning on $\S=s$; (iii) $\M(s)\perp \S \mid \{\X=\boldsymbol{x}\}$, no confounder for the relation of $\M$ and $\S$;  (iv) $\Y(s,m) \perp \M(s^*) \mid \{\X = \boldsymbol{x}\}$,
no confounder for the $\M$-$\Y$ relation that is affected by $\S$
[VanderWeele and Vansteelandt, 2009].
Under these assumptions,
the model can be visualized by the directed acyclic graph in Figure <ref>, and
the mediation effect $\eqref{eq:MEpotential}$ equals $\alphaS\betaM (s-s^*)$.

the scientific goal of detecting the presence of a mediation effect can be formulated as the following hypothesis testing problem:

\begin{align*}
H_0: \alphaS  \betaM = 0 \quad \text{\ versus\ } \quad H_A: \alphaS \betaM  \neq 0.
\end{align*}

This null hypothesis is composite and
can be decomposed into three disjoint cases:

\begin{equation}
H_0: \begin{cases}
H_{0,1}:  & \alphaS = 0, ~\betaM \neq 0; \\
H_{0,2}: & \alphaS \neq 0,~ \betaM =0; \\
H_{0,3}: & \alphaS = \betaM =0,
\end{cases} \label{eq:compositenull}
\end{equation}

and the alternative hypothesis is $H_A: \alphaS \neq 0 \text{ and } \betaM \neq 0.$

Composite null problems similar to  (<ref>) can occur in  settings other than Model (<ref>); the latter is considered to demonstrate the essential analytic details useful for possible generalizations.

Similar issues have also been observed in many other scenarios, including partially linear models [Hines et al., 2021],
survival analysis [VanderWeele, 2011], and high-dimensional models [Zhou et al., 2020].
The analytic details of the methodology development in this paper can pave the path for useful generalizations to other important statistical models and applications.

To test  the composite null (<ref>), various methods have been proposed, and a comprehensive survey can be found in  [MacKinnon et al., 2002].
There are two representative classes of tests:
(I) the product of coefficients (PoC) test, which corresponds to a Wald-type test,
and (II) the joint significance (JS) test, which  is the likelihood ratio test under  normality of the error terms [Liu et al., 2021].
The first class of methods examine the PoC: $\hatalphaSn \hatbetaMn$, where  $\hatalphaSn$ and $\hatbetaMn$  denote
consistent estimates of $\alphaS$ and $\betaM$, respectively.
One common practice is to apply a normal approximation
to $\hatalphaSn \hatbetaMn$ divided by  its standard error, where
[Sobel, 1982] derives the standard error formula by   the  first-order Delta method.
In addition to the large-sample approximation,
the bootstrap has also been used to evaluate the distribution of  $\hatalphaSn\hatbetaMn$  [MacKinnon et al., 2004, Fritz and MacKinnon, 2007].
(II)  The JS test, also known as the MaxP test, rejects $H_0: \alphaS \betaM = 0$ if  $\max\{p_{\alphaS}, p_{\betaM}\} < \omega$,
where $\omega$ is a prespecified significance level,
and $p_{\alphaS}$  and $p_{\betaM}$ denote the $p$-values for $H_0: \alphaS=0$ (the link $\S \rightarrow \M$)  and $H_0: \betaM=0$ (the link $\M \rightarrow \Y$), respectively.
Despite their popularity, these methods have been found numerically to be  overly conservative  under $H_{0,3}$ in (<ref>)   [MacKinnon et al., 2002, Barfield et al., 2017].
See a further discussion on the non-regular asymptotic behaviors of statistics underlying the conservatism in Section <ref>.

Similar issues have also been broadly recognized for Wald tests in various   statistical problems including three-way contingency table analysis and  factor analysis
[Glonek, 1993, Dufour et al., 2013, Drton and Xiao, 2016].
However, characterizing non-regular asymptotic behaviors under the singular null hypothesis $H_{0,3}$  is still insufficient to address intrinsic technical challenges in testing (<ref>).
In particular, the composite null  (<ref>) includes not only the singular case $H_{0,3}$ but also the other two non-singular cases $H_{0,1}$ and $H_{0,2}$.
Since a test statistic follows different distributions under different null cases, and the underlying true null case is unknown,
it is difficult to obtain uniformly distributed $p$-values through one simple asymptotic  distribution under (<ref>).
To address this technical difficulty,
we adopt, justify, and evaluate an adaptive bootstrap procedure.
For both  Wald-type PoC test and non-Wald JS test, we will show that the proposed procedure can naturally adapt to the three types of null hypotheses in (<ref>) and yield uniformly distributed $p$-values across all different null cases.

§ ADAPTIVE BOOTSTRAP TESTS

In this section, we propose a general  Adaptive Bootstrap (AB)  procedure for testing the composite null hypothesis (<ref>).
For illustration, we apply the  adaptive bootstrap to the representative PoC test and show it can address the  non-regularity issue.
We emphasize that this general strategy can be applied in a wide range of scenarios. We also derive adaptive bootstrap procedure  in  other examples, including the Joint Significance test (Section <ref> in the Supplementary Material), joint testing of multivariate mediators (Section <ref>) and testing mediation effect under the generalized linear models  (Sections <ref> and <ref>.).

To conduct hypothesis testing or estimate confidence intervals for statistics whose limiting distributions deviate  from the normal,
a simple and powerful approach is to apply the bootstrap resampling technique.
However, the classical bootstrap is not a panacea,  and on some occasions it can fail to work properly, including unfortunately the non-regular scenarios considered in this paper.
In particular, it has been observed through simulation studies that the classical  bootstrap technique is overly conservative under $H_{0,3}: \alphaS=\betaM=0$  [MacKinnon et al., 2002, Barfield et al., 2017].
We next unveil the key insights  underlying the failure of the classical bootstrap, which motivates our use of the adaptive bootstrap.

Non-Regularity of the PoC Test.   When  $(\alphaS, \betaM)\neq (0,0)$, one of the first-order gradients $\frac{\partial \alphaS \betaM}{\partial \alphaS} = \betaM$ and $\frac{\partial \alphaS \betaM}{\partial \betaM} = \alphaS$ is  nonzero.
Thus the Delta method  can be applied to support
the use of Sobel's test (based on asymptotic normality) and classical bootstrap [Barfield et al., 2017].
However, under $H_{0,3}: (\alphaS, \betaM)=(0,0)$,
the gradients $\frac{\partial \alphaS \betaM}{\partial \alphaS} =\frac{\partial \alphaS \betaM}{\partial \betaM} = 0$,
and validity of Sobel's test and the classical bootstrap cannot be  obtained as above.

We next illustrate the non-regular limiting behavior of the PoC $\hatalphaSn \hatbetaMn$ under $H_{0,3}$. For ease of exposition, consider a  special case of (<ref>):
$\M = \alphaS \S +\eM,$ and $\Y = \betaM \M +  \eY.$
Let $(\hatalphaSn,\hatbetaMn)$ denote the ordinary least squares estimators of $(\alphaS,\betaM)$, and let $(\hatalphaSn^*, \hatbetaMn^*)$ the corresponding nonparametric bootstrap estimators.
Here and throughout this paper,
we use the superscript “$*$" to indicate estimators obtained from the nonparametric bootstrap, namely “bootstrap in pairs" in the regression setting.
By classical asymptotic theory [van der Vaart, 2000], under mild conditions,

\begin{align}\label{eq:simplecltsep}
\sqrt{n}(\hatalphaSn - \alphaS, \, \hatbetaMn - \betaM)^{\mytrans} \xrightarrow{d} (Z_{\S,0}, \, Z_{\M,0})^{\mytrans},
\end{align}

$(Z_{S,0},Z_{M,0})^{\top}$ denotes a mean-zero normal random vector with a covariance matrix given by that of the random vector
$(\eM \S /V_{\S,0}, \eY \M/V_{\M,0})^{\mytrans}$,   $V_{\S,0}=\mathrm{E}(\S^2)$,  $V_{\M,0}=\mathrm{E}(\M^2)$.

\begin{align}\label{eq:simplecltsepboot}
\sqrt{n}(\hatalphaSn^* - \hatalphaSn, \hatbetaMn^* - \hatbetaMn)^{\mytrans}
\overset{d^*}{\leadsto}
(Z_{\S,0}', \, Z_{\M,0}')^{\mytrans},
\end{align}

where $(Z_{\S,0}', Z_{\M,0}')$ is an independent copy of  $(Z_{\S,0}, Z_{\M,0})$ in (<ref>) under the same distribution.
By (<ref>),
$n( \hatalphaS \hatbetaM - \alphaS \betaM )\xrightarrow{d} Z_{\S,0}Z_{\M,0}$ under $H_{0,3}$,
with the convergence rate $n$ different from the standard parametric $\sqrt{n}$ rate.

By (<ref>)–(<ref>),
we have
$n(\hatalphaSn^*\hatbetaMn^* - \hatalphaSn \hatbetaMn) =  n\{(\hatalphaSn^* - \hatalphaSn)\hatbetaMn + (\hatbetaMn^* - \hatbetaMn)\hatalphaSn +(\hatalphaSn^* - \hatalphaS )(\hatbetaMn^* -  \hatbetaMn)\} \overset{d^*}{\leadsto}  {Z_{\S,0}'Z_{\M,0}} + {Z_{\S,0}Z_{\M,0}'} +{Z_{\S,0}'Z_{\M,0}'}.$
We can see that the limit of $n(\hatalphaSn^*\hatbetaMn^* - \hatalphaSn \hatbetaMn)$ is different from that of
$n( \hatalphaSn \hatbetaMn - \alphaS \betaM )$,
implying inconsistency of the classical nonparametric  bootstrap.

Adaptive Bootstrap of the PoC Test.  To address the challenge of correctly evaluating the distribution of the  PoC statistic,
we utilize the local asymptotic analysis framework.
Intuitively, the goal is to evaluate if a small change in the target parameters leads to little change on the limit of the statistics.
To this end,
given targeted parameters $(\alphaS,\betaM)$,
we define their locally perturbed counterparts as $\alphaSn = \alphaS + n^{-1/2}\balpha$, and $\betaMn = \betaM + n^{-1/2}\bbeta$, respectively, where
$(\balpha, \bbeta)$ denote the local parameters of perturbations from  our targeted coefficients $(\alphaS, \betaM). $

We then frame the problem under the local linear SEM as follows:

\begin{eqnarray}
\M &=& \alphaSn \S +  \X^\mytrans \alphaX +\eM, \hspace{1.5em}  \Y=\betaMn \M + \X^\mytrans  \betaX + \DE \S + \eY,\label{eq:fullmodlocal2}
\end{eqnarray}

where $\eM$ and $\eY$ are independent error terms with mean zero and finite variances.

Fixing the target parameters $(\alphaS,\betaM)$, according to [van der Vaart, 2000]  the formulation given in (<ref>) may also be viewed as a local statistical experiment with  local parameters $(\balpha, \bbeta)$ under which we are interested in examining the limit of test statistics.
Note that with the local parameters $(\balpha,\bbeta)=(0,0)$, (<ref>) reduces to the original model  (<ref>) with  the parameters $(\alphaS,\betaM)$.
Our inference goal remains the same: that is, to test the underlying true coefficients $(\alpha_{\S}, \beta_{\M})$.  Technically, we consider a $\sqrt{n}$-vicinity of local neighboring values $(\alpha_{\S,n}, \beta_{\M, n})$ only for the theoretical investigation of  local asymptotic behaviors.
Such an idea has also been used for studying
other non-regularity issues  <cit.>.

To examine  the limit  of $\hatalphaSn \hatbetaMn - \alphaSn \betaMn$ under (<ref>), we assume the following general regularity  Condition <ref>.

(C1.1) $\mathrm{E}(\eM| \X, \S )=0$ and $\mathrm{E}(\eY| \X, \S, \M )=0$.
(C1.2) $\mathrm{E}(\boldsymbol{D}\boldsymbol{D}^{\mytrans})$
is a positive definite matrix with bounded eigenvalues,
where $\boldsymbol{D} = (\X^{\mytrans}, \M, \S)^{\mytrans}$.
(C1.3) The second moments of $(\eM, \eY, \Sperp, \Mperpp, \eM\Sperp, \eY\Mperpp)$ are finite,
where  $\Sperp = \S - \X^{\mytrans}Q_{1,\S}$  with  $Q_{1,\S}= \{\mathrm{E}(\X \X^{\mytrans})\}^{-1}\times \mathrm{E}(\X \S)$, and $\Mperpp =\M - \tilde{\X}^{\mytrans} Q_{2,\M}$   with  $\tilde{\X}=(\X^{\mytrans},\S)^{\mytrans}$ and $Q_{2,\M} = \{\mathrm{E}(\tilde{\X}^{\mytrans}\tilde{\X})\}^{-1}\times \mathrm{E}(\tilde{\X}\M)$.

Similarly to our above discussions under the simplified model,  Theorem <ref> establishes the limits of
$\sqrt{n}\times (\hatalphaSn \hatbetaMn - \alphaS \betaM)$  and $n\times (\hatalphaSn \hatbetaMn - \alphaS \betaM)$
when $(\alphaS, \betaM) \neq (0,0)$ and $(\alphaS, \betaM) = (0,0)$, respectively.

Assume Condition <ref>.
Under the local model (<ref>),

(i) when $(\alphaS, \betaM) \neq (0,0)$,
\sqrt{n}\times (\hatalphaSn \hatbetaMn - \alphaSn \betaMn) \xrightarrow{d} \alphaS {Z_{\M}} + \betaM {Z_{\S}};

(ii) when $(\alphaS, \betaM) = (0,0)$,
n\times (\hatalphaSn \hatbetaMn - \alphaSn \betaMn) \xrightarrow{d} \balpha {Z_{\M}} + \bbeta {Z_{\S}} + {Z_{\M}Z_{\S}},

where  $(Z_S, Z_M)^{\mytrans}$ is a mean-zero normal random vector with a covariance matrix given by that of the random vector $(\eM \Sperp/V_S, \eY \Mperpp/V_M)^{\mytrans}$ with $V_S=\mathrm{E}(\Sperp^2)$, and $V_M=\mathrm{E}(\Mperpp^2)$.

Theorem <ref> suggests the
of $\sqrt{n}(\hatalphaSn\hatbetaMn - \alphaSn \betaMn )$  is not uniform with respect to $(\alphaS , \betaM)$, and the non-uniformity occurs around $(\alphaS , \betaM)=(0,0)$.
On the other hand,
in the neighborhood of  $(\alphaS , \betaM)=(0,0)$, the limit of $n(\hatalphaSn\hatbetaMn - \alphaSn \betaMn )$ is continuous as a function of $(\balpha, \bbeta)^{\mytrans} \in \realnumber^2$ into the space of distribution functions. Therefore, using this local limit in the bootstrap, we expect better finite-sample accuracy, compared to the classical nonparametric bootstrap that does not take into account the local asymptotic behavior.

Moreover, to discern the null cases,
we will consider a decomposition of the statistic.
The idea is to isolate the possibility that $(\alphaS, \betaM ) \neq (0,0)$ by comparing the absolute value of the standardized statistics  $T_{\alpha,n}=\sqrt{n}\hatalphaSn/\hatsigmaalphan$ and $T_{\beta, n}=\sqrt{n}\hatbetaMn /\hatsigmabetan$ to some thresholds,
where $\hatsigmaalphan$  and $\hatsigmabetan$ denote the
sample  standard deviations  of $\hatalphaSn$ and $\hatbetaMn$, respectively.
Specifically, we decompose

\begin{align}
\hatalphaSn \hatbetaMn - \alphaSn \betaMn=&~ (\hatalphaSn \hatbetaMn - \alphaSn \betaMn )\times  (1-\myindicator_{\alphaS, \lambda_n}\myindicator_{\betaM, \lambda_n} )\label{eq:proddecomp1}  \\
& +(\hatalphaSn \hatbetaMn - \alphaSn \betaMn)\times  \myindicator_{\alphaS, \lambda_n}\myindicator_{\betaM, \lambda_n}\notag
\end{align}

with the indicators
$\myindicator_{\alphaS, \lambda_n}=\myindicator \{ |T_{\alpha,n}| \leq \lambda_{n}, \alphaS = 0\}$ and $\myindicator_{\betaM, \lambda_n}=\myindicator\{|T_{\beta,n}| \leq \lambda_{n}, \betaM = 0\}$,
where $\myindicator \{E\}$ represents the indicator function of an event $E$, and $\lambda_n$ is a certain threshold to be specified.
When $(\alphaS, \betaM) \neq (0,0)$, the classical  bootstrap is consistent for the first term in (<ref>).
For the second term in (<ref>),
we next develop a bootstrap statistic motivated by  Theorem <ref> (ii).

To construct the bootstrap statistic,
we introduce some notation following the convention in the empirical process literature [van der Vaart, 2000].
Particularly, throughout the paper,
$P_n$ denotes the population probability measure of $(\S, \X, \M, \Y)$,
$\mathbb{P}_n$ denotes the empirical measure with respect to the i.i.d. observations $\{(\S_i, \X_i, \M_i, \Y_i): i=1,\ldots, n\}$, and $\mathbb{P}_n^*$ denotes the nonparametric bootstrap version of $\mathbb{P}_n$.
For any measurable function $f(\S, \X, \M, \Y)$, we define the empirical process $\mathbb{G}_n f = \sqrt{n}(\mathbb{P}_n -P_n)f=\sqrt{n}[n^{-1}\sum_{i=1}^n f(\S_i, \X_i, \M_i, \Y_i) -\mathrm{E}\{f(\S, \X, \M, \Y)\} ]$, and its nonparametric bootstrap version is $\mathbb{G}_n^*=\sqrt{n}(\mathbb{P}_n^*- \mathbb{P}_n)$.
With the above notation, we define the sample versions of
$Q_{1,\S}, Q_{2,\M}, \Sperp,$ and $\Mperpp$ in Condition <ref> as
$\hat{Q}_{1,\S}= \{\mathbb{P}_n(\X \X^{\mytrans})\}^{-1} \mathbb{P}_n(\X \S)$,
$\hat{Q}_{2,\M} = \{\mathbb{P}_n( \tilde{\X}  \tilde{\X}^{\mytrans})\}^{-1} \mathbb{P}_n( \tilde{\X} \M)$,
$\hatSperp = \S - \X^{\mytrans}\hat{Q}_{1,\S}$, and $\hatMperpp = \M - \tilde{\X}^{\mytrans}\hat{Q}_{2,\M},$
where we use “$\hat{\quad}$” to denote the sample counterparts in this paper.
Similarly, we define their nonparametric bootstrap counterparts $(Q_{1,\S}^*, Q_{2,\M}^*, \Sperpboot, \Mperppboot )$
by replacing $\mathbb{P}_n$ with $\mathbb{P}_n^*$ in the above definitions.

When $(\alphaS, \betaM)=(0,0)$, motivated by Theorem <ref> (ii),
we construct a bootstrap statistic $\mathbb{R}_{n}^*(\balpha, \bbeta)$ as a bootstrap counterpart of of $\balpha Z_M + \bbeta Z_S+ Z_MZ_S$.
In particular, $\mathbb{R}_{n}^*(\balpha, \bbeta) = {\balpha \mathbb{Z}_{\M,n}^*}  +{\bbeta\mathbb{Z}_{\S,n}^*} + {\mathbb{Z}_{\S,n}^*\mathbb{Z}_{\M,n}^* },$
$\mathbb{Z}_{\S,n}^*= \mathbb{G}_n^*(\hatenM \Sperp^*)/\mathbb{V}_{\S, n}^*$, $\mathbb{Z}_{\M,n}^*=\mathbb{G}_n^*(\hatenY \Mperpp^*)/\mathbb{V}_{\M, n}^*,$ $\mathbb{V}_{\S, n}^* = \mathbb{P}_n^*\{(\Sperp^*)^2\},$ $\mathbb{V}_{\M, n}^* = \mathbb{P}_n^*\{(\Mperpp^*)^2\}$,
$\hatenM$ and $\hatenY$ denote the sample residuals obtained from the ordinary least squares regressions of the two models in (<ref>).
When $(\alphaS, \betaM)\neq (0,0)$, we still consider
the classical nonparametric bootstrap estimator $\hatalphaSn^* \hatbetaMn^*$.
To develop an adaptive bootstrap test, we utilize the decomposition (<ref>) and
propose to replace the indicators $\myindicator_{\alphaS, \lambda_n}$ and $\myindicator_{\betaM, \lambda_n}$ in (<ref>) by

\begin{align}
\myindicator_{\alphaS, \lambda_{n}}^*  = \myindicator \{ |T_{\alpha,n}^*|\leq \lambda_{n}, ~|T_{\alpha,n}|\leq \lambda_{n} \} \quad \text{and} \quad  	\myindicator_{\betaM, \lambda_{n}}^* = \myindicator \{ |T_{\beta,n}^*|\leq \lambda_{n},~ |T_{\beta,n}|\leq \lambda_{n} \}, \label{eq:indicatoralphabetasep}
\end{align}

where $T_{\alpha,n}^*=\sqrt{n}\hatalphaSn^*/\hatsigmaalphan^*$ and $T_{\beta, n}^*=\sqrt{n}\hatbetaMn^* /\hatsigmabetan^*$
denote the classical nonparametric bootstrap versions of  $T_{\alpha, n}$ and $T_{\beta, n}$, respectively.
Following the decomposition in (<ref>), we define a statistic

\begin{align*}
U^*=(\hatalphaSn^* \hatbetaMn^*  - \hatalphaSn \hatbetaMn ) \times (1- \myindicator_{\alphaS, \lambda_n}^* \myindicator_{\betaM, \lambda_n}^*  ) + n^{-1} \mathbb{R}_{n}^*( \balpha, \bbeta ) \times \myindicator_{\alphaS, \lambda_n}^*  \myindicator_{\betaM, \lambda_n}^*,
\end{align*}

termed as Adaptive Bootstrap (AB) test statistic in this paper.
Theorem <ref> below establishes the bootstrap consistency of $U^*$.

Assume the conditions of Theorem <ref> are satisfied. When  $\lambda_n=o(\sqrt{n})$ and $\lambda_n\to \infty$ as $n\to \infty$,

$$c_nU^* \overset{d^*}{\leadsto} c_n(\hatalphaSn \hatbetaMn - \alphaS \betaM),$$

where $c_n$ is a non-random scaling factor satisfying

\begin{align}\label{eq:cndef}
c_n=	\begin{cases}
\sqrt{n}, & \text{ when }(\alphaS, \betaM) \neq (0,0) \\
n, &  \text{ when } (\alphaS, \betaM) = (0,0)
\end{cases}.
\end{align}

Theorem <ref> suggests that
under the original model (<ref>), i.e., $(\balpha, \bbeta)=(0,0)$,
the AB statistic $U^*$ is a consistent bootstrap estimator for $\hatalphaSn \hatbetaMn - \alphaS \betaM$ with a proper scaling.
for any fixed targeted parameters $(\alphaS, \betaM)$, in their local neighborhoods,
i.e.,  $(\balpha, \bbeta)\neq (0,0)$, the bootstrap consistency still holds as a smooth function of  $(\balpha, \bbeta)$.
Intuitively, this suggests that  a small change in the target parameters does not affect the consistency property, and $U^*$ is “regular” under the local model.

In practice,  without knowing which case is the true null we rely on $U^*$ as the bootstrap statistic for $\hatalphaSn \hatbetaMn - \alphaSn \betaMn$ generally.
This strategy is viable because with a given finite  sample size $n$,
using $\sqrt{n}U^*$ for bootstrapping  $\sqrt{n}(\hatalphaSn \hatbetaMn - \alphaSn \betaMn)$ is  equivalent to using $nU^*$ for bootstrapping  $n(\hatalphaSn \hatbetaMn - \alphaSn \betaMn)$.
Therefore, as desired, $U^*$ will approximate well the distribution of $\hatalphaSn \hatbetaMn - \alphaSn \betaMn$ regardless of the underlying null case.

As a comparison, we also discuss the naive non-parametric bootstrap when $(\alphaS, \betaM) = (0,0)$. Specifically, we obtain the following expression (in Remark <ref> of the Supplementary Material),

\begin{align} \label{eq:classboot00expand}
n(\hatalphaSn^* \hatbetaMn^* - \hatalphaSn \hatbetaMn)%\notag\\
=&~\mathbb{R}_{n}^*(\balpha, \bbeta)+  \mathbb{Z}_{\S,n}\mathbb{Z}_{\M,n}^* + \mathbb{Z}_{\M,n}\mathbb{Z}_{\S,n}^*,
\end{align}

$\mathbb{Z}_{\S,n}= \mathbb{G}_n(\eM \hatSperp)/\mathbb{V}_{\S, n}$, $\mathbb{Z}_{\M,n}=\mathbb{G}_n(\eY \hatMperpp)/\mathbb{V}_{\M, n},$ $\mathbb{V}_{\S, n} = \mathbb{P}_n(\hatSperp^2),$ and $\mathbb{V}_{\M, n} = \mathbb{P}_n(\hatMperpp^2)$.
In addition to the term $\mathbb{R}_{n}^*(\balpha, \bbeta)$, (<ref>) has two extra random terms $\mathbb{Z}_{\S,n}\mathbb{Z}_{\M,n}^* + \mathbb{Z}_{\M,n}\mathbb{Z}_{\S,n}^*$, which suggests that using (<ref>) in the bootstrap would not be consistent.

The issue of the classical bootstrap being inconsistent at  $(\alphaS, \betaM) = (0,0)$ is circumvented by the proposed local bootstrap statistic $\mathbb{R}_n^*(\balpha, \bbeta )$.

Adaptive Bootstrap Test Procedure.
We introduce a consistent bootstrap test procedure for $\hatalphaSn \hatbetaMn$
based on
Theorem <ref>.
Given a nominal level $\omega$, let $q_{\omega/2}$ and $q_{1-\omega/2}$ denote the lower and upper $\omega/2$ quantiles, respectively, of the bootstrap estimates $U^*$. If $\hatalphaSn \hatbetaMn$ falls outside  the interval $(q_{\omega/2}, q_{1-\omega/2} )$, we reject the composite null (<ref>), and conclude that the mediation effect is statistically significant at the level $\omega$.
We clarify that the goal is to test the underlying true coefficients $(\alpha_{\S}, \beta_{\M})$. The reason to consider their $\sqrt{n}$-local coefficients $(\alpha_{\S,n}, \beta_{\M, n})$ is merely for theoretical investigation of local asymptotic behaviors. Therefore, to test (<ref>) under the original model (<ref>), it suffices to calculate $U^*$ with $\balpha=\bbeta=0$.
We point out that the rejection region in the adaptive procedure may also be constructed  through the asymptotic distribution as an alternative to the bootstrap; nevertheless,
the proposed bootstrap procedure is more flexible and does not rely on a particular form of the limiting distributions,  and therefore, it can be easily extended under various mediation models; see more examples in Section <ref>.

Choice of the Tuning Parameters.
Under the conditions of Theorem <ref>, which specify   $\lambda_n=o(\sqrt{n})$ and $\lambda_n\to \infty$ as $n\to \infty,$ we have
$\lim_{n\to \infty} \Pr(|T_{\alpha, n}|>\lambda_n, |T_{\beta, n}|>\lambda_n \mid  \alphaS = \betaM = 0 ) = 0$,
suggesting that $\myindicator_{\alphaS, \lambda_n}\myindicator_{\betaM, \lambda_n}$ can provide a consistent test for $\alphaS=\betaM=0$.
If $\lambda_n$ remains bounded as $n\to \infty$,
$U^*$ asymptotically reduces to $\hatalphaSn^*\hatbetaMn^*-\hatalphaSn\hatbetaMn$, i.e., the classical nonparametric bootstrap procedure, which  may be conservative.
In the simulation experiments, we set $\lambda_n=\lambda\sqrt{n}/\log n$ and find that a fixed constant $\lambda$, e.g.,  $\lambda=2$ can give a good performance.
In general settings,
we can choose the tuning parameter through the double bootstrap [Chen, 2016]; see Section <ref> of the Supplementary Material for more implementation details.

Our proposed adaptive procedures examine the non-regular asymptotic behaviors of test statistics
through local models.
In effect, the idea of local models may be traced back to  econometrics [Andrews, 2001] and was utilized in other statistical problems, such as classification and post-selection inference
[Laber and Murphy, 2011, McKeague and Qian, 2015, McKeague and Qian, 2019, Wang et al., 2018].
we emphasize that there are unique statistical challenges of
testing mediation effects.
First, in terms of the parameter space,
the null hypothesis of no mediation effect is essentially a union of individual hypotheses.
This results in a  non-standard shape of the null parameter space, on which both regular and non-regular asymptotic behaviors can occur, as illustrated in  Figure  <ref>.
Second, in terms of the behavior of the estimator,
we unveil a fundamental zero-gradient phenomenon.
This is caused by the special form of the product statistic and cannot be directly addressed by  the existing adaptive procedures mentioned above.
Third, in terms of the models,
the mediation analysis involves a system of  structural equations.
Ignoring the model structure in the implementation could lead to slow computation;  see Section <ref> of the Supplementary Material for more  details on computation.
Due to these unique challenges, new developments  in methodology, theory, and computation are necessary.

In addition to the Wald-type PoC test,
we also address the non-regularity issue of the non-Wald joint significance/maxP test through our proposed adaptive bootstrap.
It is noteworthy that non-regular behaviors of the JS and PoC tests under the singleton  $H_{0,3}$ are distinct, as the two statistics take different forms.
Particularly, PoC statistic has the zero-gradient issue discussed above, whereas JS statistic has a certain inconsistent convergence issue.
Despite that difference, we can similarly develop an adaptive bootstrap for  the JS test and obtain uniformly distributed $p$-values.
Refer to the detail in Section  <ref> of  the Supplementary Material.
This suggests that our proposed adaptive bootstrap is not restricted to the Wald-type test and may be further generalized to other tests with similar circumstances.

It is worth noting that the proposed strategy can be generalized to deal with multiple   mediators under suitable identifiability conditions.
In the following,
we delve into three scenarios of practical importance.

(i) We consider the group-level joint Mediation Effect (ME) via a set of mediators $\boldsymbol{M}=(M_1,\ldots, M_J)$ shown by the red path in Figure <ref> below.
This type of joint ME has been considered in the literature by [Huang and Pan, 2016] and [Hao and Song, 2022], among others.
We  generalize the AB method to test the joint mediation  effect in Section <ref>.

(ii) We consider multiple mediators that are causally uncorrelated [Jérolon et al., 2020] or governed by  the parallel path model [Hayes, 2017].
In this case,
the indirect effect of one single mediator
can be identified under the known identifiability assumptions outlined in [Imai and Yamamoto, 2013]. In particular, under
the multivariate  linear SEM (<ref>) with no
causal interplay between mediators,
the null hypothesis of
no individual indirect effect via one mediator, say, $M_1$, could be formed as
$H_0:\alpha_{S,1}\beta_{M,1}=0$, illustrated in Figure <ref> below.
To apply the AB test to $\alpha_{S,1}\beta_{M,1}$,
we note that (<ref>) can be equivalently rewritten as
$ and
Y=\beta_{M,1}M_1+ \boldsymbol{M}_{(-1)}\boldsymbol{\beta}_{(-1)} + \X \boldsymbol{\beta}_{\X} +\tau_S S+\epsilon_Y,$
where  $\boldsymbol{\beta}_{(-1)}=(\beta_{M,2},\ldots, \beta_{M,J})^{\top}$ and $ \boldsymbol{M}_{(-1)}=(M_2,\ldots, M_J)^{\top}$.
This form resembles (<ref>), and
the AB method in Section <ref> can be employed to test $\alpha_{S,1}\beta_{M,1}=0$ by adjusting $(\boldsymbol{M}_{(-1)},\X)$ in the outcome model.
We provide details including the identification assumptions in Section <ref> of the Supplementary Material.

(iii) When the mediators are causally correlated, evaluating individual indirect effects along different posited paths requires correct specification of the mediators' causal structure  [VanderWeele et al., 2014].
To relax such stringent assumptions,
researchers have proposed alternative methods, one of which is to examine the interventional indirect effects specific to each distinct mediator [Loh et al., 2021].
the interventional indirect effect via a target mediator $M_1$ is supposed to capture all of the exposure-outcome effects that are mediated by $M_1$ as well as any other mediators causally preceding $M_1$;
see a diagram in Figure <ref> below.
Under a typical class of linear and additive mean models, estimators of interventional indirect effects take the same product form   of coefficients as that in the above Case (ii).
Thus the proposed AB method in Section <ref> can be similarly applied with little effort.
We provide relevant details including the definition and identification assumptions of the interventional indirect effects in Section <ref> of the Supplementary Material.

[shape=circle, thick] (S) at (-2.2,0) $S$;

(M_all) at (0,0) $
\boldsymbol{M}=\left\lbrack\begin{aligned}
\vdots\ \, \\
\end{aligned}
\right\rbrack$;
[shape=circle, thick] (Y) at (2.2,0) $Y$;

[shape=circle, thick] (X) at (0,-2.1) $\boldsymbol{X}$;

[-stealth, thick, red] (S.east) – node[midway,above] $\boldsymbol{\alpha}_S$ (M_all.west);
[-stealth, thick, red] (M_all.east) – node[midway,above] $\boldsymbol{\beta}_M$ (Y.west);
[-stealth, thick,black!20] (X.west) .. controls + (left:8mm) and + (down:20mm) ..  node[midway,above]  (S);
[-stealth, thick,black!20] (X) – (M_all);
[-stealth, thick,black!20] (X.east) .. controls + (right:8mm) and + (down:20mm) .. node[midway,above]  (Y);

[-stealth, thick] (S) .. controls + (up:20mm) and + (up:20mm)  ..  (Y);

Joint ME via a group of mediators $(M_1,\ldots, M_J)$

[shape=circle, thick] (S) at (-2.2,0) $S$;

(M_all) at (0,0) $M_1$;

[shape=circle, thick] (Y) at (2.2,0) $Y$;

[thick] (X) at (0,-2.1) $\boldsymbol{X}$;

[thick] (Mother) at (-1.2,-1.2) $\boldsymbol{M}_{(-1)}$;

[-stealth, thick, red] (S.east) – node[midway,above] $\alpha_{S,1}$ (M_all.west);

[-stealth, thick, red] (M_all.east) – node[midway,above] $\beta_{M,1}$ (Y.west);
[-stealth, black!20, thick] (X.west) .. controls + (left:8mm) and + (down:20mm) .. node[midway,above]  (S);
[-stealth, thick, black!20] (X.north) – (M_all.south);

[dashed, thick, black!40] (Mother.north) – (M_all.215);

[-stealth, thick, black!20] (X.175) – (Mother.south);

[-stealth, thick, black!20] (X.east) .. controls + (right:8mm) and + (down:20mm)  ..  node[midway,above]  (Y);

[-stealth, thick] (S) .. controls + (up:20mm) and + (up:20mm)  ..  (Y);

Individual ME via one mediator $M_1$

[shape=circle, thick] (S) at (-2.2,0) $S$;

(M_all) at (0,0) $M_1$;

[shape=circle, thick] (Y) at (2.2,0) $Y$;

[thick] (X) at (0,-2.1) $\boldsymbol{X}$;

[thick] (Mother) at (-1.2,-1.4) $\boldsymbol{M}_{1,\text{pre}}$;

[-stealth, thick, red] (S.east) – node[midway,above]  (M_all.west);

[-stealth, thick, red] (M_all.east) – node[midway,above]  (Y.west);

[-stealth, black!20, thick] (X.west) .. controls + (left:8mm) and + (down:20mm) .. node[midway,above]  (S);
[-stealth, thick, black!20] (X.north) – (M_all.south);

[-stealth, thick, red] (S.east) – node[midway,above]  (Mother.102);

[-stealth, thick, red] (Mother.78) – (M_all.215);

[-stealth, thick, black!20] (X.175) – (Mother.south);

[-stealth, thick, black!20] (X.east) .. controls + (right:8mm) and + (down:20mm)  ..  node[midway,above]  (Y);

[-stealth, thick] (S) .. controls + (up:20mm) and + (up:20mm)  ..  (Y);

Interventional ME via one mediator $M_1$

Path diagram of the mediation model with multiple mediators: Dashed lines represent possible non-causal correlations or independence, and solid arrowed lines represent possible causal relationships.
In Panel (c), $\boldsymbol{M}_{1,\operatorname{pre}}$ represents mediators that are causally preceding $M_1$.

§ NUMERICAL EXPERIMENTS

In this section, we conduct simulation experiments to evaluate the finite-sample performance of the proposed adaptive bootstrap PoC and JS tests.
Particularly, we generate data through the following model:

\begin{eqnarray}
\M &=&  \alphaS \S + \alpha_I+ \alpha_{X,1} X_1 + \alpha_{X,2} X_2  +\eM, \label{eq:simmodel}  \\
\Y & = & \betaM \M +  \beta_I + \beta_{X,1} X_1 + \beta_{X,2} X_2  + \DE \S + \eY.  \notag
\end{eqnarray}

In the model (<ref>), the exposure variable $\S$ is simulated from the Bernoulli distribution with the success probability 0.5; the covariate $X_1$ is continuous and simulated from a standard normal distribution $\mathcal{N}(0,1)$; the covariate $X_2$ is discrete and  simulated from  the Bernoulli distribution with the success probability 0.5; two error terms $\eM$ and $\eY$ are simulated independently from $\mathcal{N}(0,\sigma_{\eM}^2)$ and $ \mathcal{N}(0,\sigma_{\eY}^2)$, respectively.
We set the  parameters  $(\alpha_I, \alpha_{X,1}, \alpha_{X,2}) = (1, 1, 1)$, $(\beta_I, \beta_{X,1}, \beta_{X,2}) = (1, 1, 1)$, $\DE=1$, and $\sigma_{\eY}=\sigma_{\eM}=0.5$.
Moreover,  we consider sample sizes $n\in \{200, 500\}$, and set the bootstrap sample size at 500.

In simulation studies, we compare eight testing methods:
the adaptive bootstrap for the PoC test (PoC-AB),
the classical nonparametric bootstrap for the PoC test (PoC-B),
Sobel's test (PoC-Sobel),
the  adaptive bootstrap  for the JS test (JS-AB),
the classical nonparametric bootstrap for the JS test (JS-B),
the MaxP test (JS-MaxP),
the nonparametric bootstrap method in the causal mediation analysis R package  [Tingley et al., 2014] (CMA),
and the method in
[Huang, 2019] (MT-Comp).
It is noteworthy that [Huang, 2019]'s MT-Comp made specific model assumptions,
which are not fully compatible with our simulation  settings, and we include this method just for the purpose of comparison.
Some other methods <cit.> relied on estimating the relative proportions of the three cases, which is  not directly  applicable here  and thus not included.

§.§ Null Hypotheses: Type 1 Error Rates

Setting 1: Under a fixed type of null.   In the first setting, we
simulate data under a fixed null hypothesis over 2000 Monte Carlo replications  to  estimate the distribution of $p$-values.
Particularly, we consider
three types of null hypotheses below:

\begin{align}
H_{0,1}: (\alphaS, \betaM)=(0,0.5),\quad \ H_{0,2}: (\alphaS, \betaM)=(0.5,0), \quad \ H_{0,3}: (\alphaS, \betaM)=(0,0). \label{eq:simcompnull}
\end{align}

We draw the Q-Q plots with $n=200$ in  Figure <ref>.
QQ-plots under $n=500$ are similar and presented in Figure  <ref> of the Supplementary Material.

In Figure <ref>, three subfigures in the first row present the results of the PoC tests under three fixed  nulls $H_{0,1}, H_{0,2}$, and $H_{0,3}$, respectively, and  three subfigures in  the second row present the corresponding results of the JS tests, respectively.

Q-Q plots of $p$-values under the fixed null with $n = 200$.
$H_{0,1}: (\alphaS, \betaM)=(0, 0.5)$

$H_{0,2}: (\alphaS, \betaM)=(0.5, 0)$

$H_{0,3}: (\alphaS, \betaM)=(0,0)$

Figure <ref>  shows that for the PoC type of tests,
under $H_{0,1}$ or $H_{0,2}$,
the PoC-AB, the  PoC-B, and the PoC-Sobel can correctly approximate the distribution of the PoC test statistic.
However, under $H_{0,3}$,
the PoC-B and the PoC-Sobel  become conservative, while
the proposed PoC-AB still approximates the distribution of the PoC statistic well.
Similarly, for the JS type of tests,
under $H_{0,1}$ or $H_{0,2}$,
the JS-AB, the JS-B, and the JS-MaxP all work well.
In contrast, under $H_{0,3}$,
the JS-B inflates,
and the JS-MaxP becomes conservative,
while the JS-AB still exhibits a good performance.
In addition, Figures <ref> and <ref> also display the results of both [Huang, 2019]'s MT-Comp and the causal mediation analysis R package CMA [Tingley et al., 2014] for comparison.
We observe that the MT-Comp  properly controls the type 1 error under $H_{0,3}$, but fails to do so   under $H_{0,2}$ and $H_{0,3}$ with inflated type I errors.  This may be because the models considered in [Huang, 2019] are not compatible with our simulation settings.
On the other hand, the causal mediation R package [Tingley et al., 2014] produces uniformly distributed $p$-values under $H_{0,1}$ and $H_{0,2}$,
but is conservative under $H_{0,3}$.  This means that the R package CMA test is underpowered.

Setting 2: Under a random type of null.   In the second setting, we simulate data
over 2000 Monte Carlo replications, where  in each replication, the null hypothesis is not fixed but   randomly  selected from   $H_{0,1}$–$H_{0,3}$ in  (<ref>).
Specifically, for $(H_{0,1},H_{0,2}, H_{0,3})$,
we consider
three selection probabilities (I) $(1/3, 1/3, 1/3)$, (II) $(0.2, 0.2, 0.6)$, and (III) $(0.05, 0.05, 0.9)$,  respectively.
We provide QQ-plots of $p$-values with $n=200$ in Figure <ref>, and QQ-plots under  $n=500$ are similar and  provided in Figure <ref> of the Supplementary Material.
In Figure <ref>, three subfigures in the first row present the results of the PoC tests with three null selection  probabilities (I)–(III), respectively,
and three subfigures in the second row present the corresponding results of the JS test, respectively.

Figures <ref> shows that
the adaptive bootstrap procedures for the PoC and JS tests perform well under different settings.
The  PoC-B test, PoC-Sobel's test, the JS-MaxP test, and the R package CMA [Tingley et al., 2014] are conservative, and they become more conservative as the probability of choosing $H_{0,3}$ increases.
We mention that in many biological studies such as genomics, $H_{0,3}$  predominates the null cases,
hence these  tests that are conservative under $H_{0,3}$ may not be preferred.
Moreover, the JS-B test and the MT-Comp method can have inflated type 1 errors.  The performance of JS-B becomes worse as the proportion of $H_{0,3}$ rises, while  the MT-Comp method deteriorates as the proportions of $H_{0,1}$ and $H_{0,2}$ increase.

Q-Q plots of $p$-values under the mixture of  nulls: $n = 200$.
(I) $(1/3, 1/3, 1/3)$

(II) $(0.2, 0.2, 0.6)$

(III) $(0.05, 0.05, 0.9)$

§.§ Alternative Hypotheses: Statistical  Power

In this subsection, we evaluate the statistical  power  of the proposed AB tests under alternative hypotheses.
Particularly, we simulate data under two settings:
(I) fix $\alphaS=\betaM$ for the convenience of pictorial presentation, which   takes various values  beginning from zero;
(II) fix the size of the mediation effect $\alphaS \betaM$,  and vary the ratio $\alphaS/ \betaM$.
In the setting (I), we consider $n \in \{200, 500\}$, and then plot the empirical  rejection rates, based on 500 Monte Carlo replications, versus the signal size of $\alphaS$, which is equal to $ \betaM$ in the setting (I).
In the setting (II), we fix  $\alphaS \betaM=0.04$ when $n=200$, and $\alphaS \betaM=0.015$ when $n=500$.
Then we plot the empirical  rejection rates versus the ratio $\alphaS/\betaM$.
The results in the two settings (I) and (II) are shown in Figures  <ref> and  <ref>, respectively.

Figures <ref> and  <ref> show that for the three PoC tests,  the PoC-AB has higher power than that of the classical nonparametric bootstrap, and both are more powerful than the Sobel's test.
Similarly, for the JS tests,
the JS-AB has higher power than that of the  classical bootstrap, and both have higher power than the MaxP test. In addition, the JS-B test   has slightly inflated type 1  errors when $(\alphaS,  \betaM)=(0,0)$, which is consistent with the results in Figure <ref>.
Among the three classical methods
(Sobel's test, the MaxP test, and the PoC-B), the MaxP test seems to achieve the best balance between the type 1 error and the statistical power, while Sobel's test has the lowest power.
These findings are consistent with those reported in the current literature [MacKinnon et al., 2002, Barfield et al., 2017].
[Huang, 2019]'s  MT-Comp test
has shown seriously inflated type 1 errors
in Figure <ref>, and therefore is not a fair competitor in our considered settings  despite its high power. Overall, it is clear that the proposed PoC-AB and JS-AB tests are superior over these existing methods, with the most robust control of type I error and highest power.

Empirical rejection rates (power) versus the signal strength of $\alphaS = \betaM$.

Empirical rejection rates (power) versus the ratio $\alphaS/\betaM$.

§ EXTENSIONS

The adaptive bootstrap in Section <ref> offers a general strategy that can be extended in a wide range of scenarios beyond the model (<ref>).
We next examine three examples, including  testing the joint mediation effect of multivariate mediators in Section <ref>,
testing the mediation effect in terms of odds ratio for a binary outcome in Section <ref>,
and  testing the mediation effect in terms of risk difference when the outcome is continuous, and the mediator follows a generalized linear model in Section <ref>.
In each scenario, we present details in the order of (1) Model, (2) Non-regularity issue, (3) Asymptotic theory and adaptive bootstrap, and (4) Numerical results.

§.§ Testing Joint Mediation Effect of  Multivariate Mediators

When the number of mediators is large,
it can also be of interest to conduct  group-based mediation analyses for a set of mediators
[VanderWeele and Vansteelandt, 2014, Daniel et al., 2015, Huang and Pan, 2016, Sohn and Li, 2019, Hao and Song, 2022];
also see a review  in [Blum et al., 2020].
In this section, we show that the proposed AB method can be generalized to  test joint mediation effects.

(1) Model.
As an extension of (<ref>),
we consider the multivariate linear SEM [VanderWeele and Vansteelandt, 2014, Huang and Pan, 2016, Hao and Song, 2022],

\begin{align}
\M_j =\alpha_{\S,j} \S +  \X^\mytrans \alpha_{\X,j} +\error_{\M,j},	\hspace{2em}\Y  = \sum_{j=1}^{\nM} \beta_{\M,j} \M_j + \X^\mytrans  \betaX + \DE \S + \eY, \label{eq:fullmodmult1}
\end{align}

where $\X$ denotes a vector of confounders with the first element being 1 for the intercept,
$\eY$ and $\boldsymbol{\epsilon}_M:=(\epsilon_{M,1},\ldots, \epsilon_{M,J})^{\top}$ are independent error terms with mean zero, $\mathrm{var}(\eY)=\sigmaeY^2$, and $\mathrm{cov}(\boldsymbol{\epsilon}_M)=\boldsymbol{\Sigma}_M$.
Assume identification conditions similar to those in Section <ref> (see Condition <ref> in the Supplementary Material). The joint mediation effect through the group of mediators $\boldsymbol{M}$ is $\mathrm{E}\{\Y(s, \boldsymbol{\M}(s)) -\Y(s, \boldsymbol{\M}(s^*))\}=(s-s^*)\boldsymbol{\alpha}_{\S}^{\mytrans}\boldsymbol{\beta}_{\M}$  [Huang and Pan, 2016], where $\boldsymbol{\alpha}_{S}=(\alpha_{S,1},\ldots, \alpha_{S,J})^{\top}$ and  $\boldsymbol{\beta}_{M}=(\beta_{M,1},\ldots, \beta_{M,J})^{\top}$.

(2) Non-Regularity Issue.
We are interested in $H_0:$  joint mediation effect  $=0$, which is equivalent to $H_0: \boldsymbol{\alpha}_{\S}^{\mytrans}\boldsymbol{\beta}_{\M}=0$.
Similarly to Section <ref>,
when $(\boldsymbol{\alpha}_{\S}, \boldsymbol{\beta}_{\M}) \neq \boldsymbol{0}$, i.e., there exists at least one coefficients $\alpha_{S,j}\neq 0$ or $\beta_{M,j}\neq 0$,
we have
${\partial  (\boldsymbol{\alpha}_{\S}^{\mytrans}\boldsymbol{\beta}_{\M})}/{\partial \alpha_{S,j}}=\beta_{M,j}\neq 0$ or ${\partial  (\boldsymbol{\alpha}_{\S}^{\mytrans}\boldsymbol{\beta}_{\M})}/{\partial \beta_{M,j}} =\alpha_{S,j}\neq 0$.
However, when $(\boldsymbol{\alpha}_{\S}, \boldsymbol{\beta}_{\M}) = \boldsymbol{0}$, i.e.,
$\alpha_{S,j}=\beta_{M,j}=0$ for all $j\in \{1,\ldots, J\}$,
${\partial (\boldsymbol{\alpha}_{\S}^{\mytrans}\boldsymbol{\beta}_{\M})}/{\partial \alpha_{S,j}} = {\partial  (\boldsymbol{\alpha}_{\S}^{\mytrans}\boldsymbol{\beta}_{\M})}/{\partial \beta_{M,j}}=0$ for all $j\in\{1,\ldots, J\}$.
We expect that  a non-regularity issue similar to that in Section <ref> would occur when $(\boldsymbol{\alpha}_{\S}, \boldsymbol{\beta}_{\M}) = \boldsymbol{0}$.
This issue is also illustrated by numerical experiments in Section <ref> of the Supplementary Material.

(3) Asymptotic Theory and Adaptive Bootstrap.
To better understand the non-regularity issue, we similarly consider a local linear SEM $\M_j =\alpha_{\S,j,n} \S +  \X^\mytrans \alpha_{\X,j} +\error_{\M,j},$ and $\Y  = \sum_{j=1}^{\nM} \beta_{\M,j,n} \M_j + \X^\mytrans  \betaX + \DE \S + \eY$, where  ${\alpha}_{\S, j, n}={\alpha}_{\S,j} +n^{-1/2}{b}_{\alpha,j}$ and ${\beta}_{\M, j, n}={\beta}_{\M,j} + n^{-1/2}{b}_{\beta,j}$.

Under Conditions <ref> and  <ref> (the latter is a regularity condition on the design matrix similar to Condition <ref>), and the local model,

when $(\boldsymbol{\alpha}_{\S}, \boldsymbol{\beta}_{\M}) \neq \boldsymbol{0}$,
\sqrt{n}\times (\hat{\boldsymbol{\alpha}}_{\S,n}^{\mytrans}\hat{\boldsymbol{\beta}}_{\M,n}-{\boldsymbol{\alpha}}_{\S,n}^{\mytrans}{\boldsymbol{\beta}}_{\M,n})\xrightarrow{d} {\boldsymbol{\alpha}}_{\S}^{\mytrans} \vec{Z}_{\M}+\boldsymbol{\beta}_{\M}^{\mytrans}\vec{Z}_{\S};$

(ii) when $(\boldsymbol{\alpha}_{\S}, \boldsymbol{\beta}_{\M}) = \boldsymbol{0}$,
n\times (\hat{\boldsymbol{\alpha}}_{\S,n}^{\mytrans}\hat{\boldsymbol{\beta}}_{\M,n}-{\boldsymbol{\alpha}}_{\S,n}^{\mytrans}{\boldsymbol{\beta}}_{\M,n})\xrightarrow{d}  \boldsymbol{b}_{0,\alpha}^{\mytrans}\vec{Z}_{\M}+\boldsymbol{b}_{0,\beta}^{\mytrans}\vec{Z}_{\S} + \vec{Z}_{\S}^{\mytrans}\vec{Z}_{\M} ,$

where $(\vec{Z}_{\S}, \vec{Z}_{\M})$ are defined to be multivariate counterparts of $({Z}_{\S},{Z}_{\M})$ in Theorem <ref>, and the detailed definitions are given  in Section <ref> of the Supplementary Material.

To present the theory of bootstrap consistency, we define the  multivariate counterparts of
${\mathbb{R}}_{n}^*(\balpha, \bbeta)$ in Section <ref> as
$\vec{\mathbb{R}}_{n}^*(\boldsymbol{b}_{\alpha}, \boldsymbol{b}_{\beta})$.
The detailed forms are given in Section <ref> of the Supplementary Material.
Similarly to $U^*$ in Section <ref>,
we define the AB statistic under the multivariate setting as

\vec{U}^*= (\hat{\boldsymbol{\alpha}}_{\S,n}^{*\mytrans} \hat{\boldsymbol{\beta}}_{\M, n}^*  - \hat{\boldsymbol{\alpha}}_{\S,n}^{\mytrans} \hat{\boldsymbol{\beta}}_{\M, n}) \times (1- 	\vec{\myindicator}_{\lambda_n}^* ) + n^{-1} \vec{\mathbb{R}}_{n}^*( \boldsymbol{b}_{\alpha}, \boldsymbol{b}_{\beta} ) \times 	\vec{\myindicator}_{\lambda_n}^*,

\vec{\myindicator}_{\lambda_n}^*=\myindicator\{\, \max\{|T_{\alpha,j,n}|,\, |T_{\alpha,j,n}^*|,\, |T_{\beta,j,n}|,\, |T_{\beta,j,n}^*|:\, 1\leqslant j\leqslant \nM \}\leqslant \lambda_n \, \},
where $T_{\alpha,j,n}=\sqrt{n}\hat{\alpha}_{\S, j}/\hat{\sigma}_{\alphaS,j,n}$ and $T_{\beta, j, n}=\sqrt{n}\hat{\beta}_{\M,j,n} /\hat{\sigma}_{\betaM,j,n}$ denote the sample T-statistics of the two coefficients $\alpha_{S,j}$ and $\beta_{M,j}$, respectively, and
$T_{\alpha,j, n}^*=\sqrt{n}\hat{\alpha}_{\S, j}^*/\hat{\sigma}_{\alphaS,j,n}^*$ and $T_{\beta, j, n}^*=\sqrt{n}\hat{\beta}_{\M,j,n}^* /\hat{\sigma}_{\betaM,j,n}^*$
denote the bootstrap counterparts of the two sample  T-statistics.
We establish bootstrap consistency for the joint AB statistic $\vec{U}^*$ below.

Under the conditions of 	Theorem <ref>, when the tuning parameter $\lambda_n$ satisfies $\lambda_n=o(\sqrt{n})$ and $\lambda_n\to \infty$ as $n\to \infty,$

$$c_n\vec{U}^* \overset{d^*}{\leadsto} c_n(\hat{\boldsymbol{\alpha}}_{\S,n}^{\mytrans} \hat{\boldsymbol{\beta}}_{\M, n} - {\boldsymbol{\alpha}}_{\S,n}^{\mytrans} {\boldsymbol{\beta}}_{\M, n}),$$

where $c_n$ is specified as in (<ref>).

Based on Theorem <ref>, we can develop an AB test  similar to that in Section <ref>.

(4) Numerical Results.
To evaluate the performance of the joint AB test,
we conduct numerical experiments, detailed in Section <ref> of the Supplementary Material.
We compare the AB test with the classical bootstrap and two tests in [Huang and Pan, 2016]: the Product Test based on Normal Product distribution (PT-NP) and  the
Product Test based on Normality (PT-N).
We observe results similar to those in Section <ref>.
Specifically, under $H_0: \boldsymbol{\alpha}_S^{\top}\boldsymbol{\beta}_M=0$, when $(\boldsymbol{\alpha}_S, \boldsymbol{\beta}_M)\neq \mathbf{0}$, both the proposed AB test and the compared methods yield uniformly distributed $p$-values.
However, when  $(\boldsymbol{\alpha}_S, \boldsymbol{\beta}_M)= \mathbf{0}$, the compared methods become overly conservative, whereas the AB test still produces uniformly distributed $p$-values.
Under $H_A$, the AB test can achieve higher empirical power than the compared methods.
Besides simulations, we also provide an exemplary data analysis in Section  <ref> of the Supplementary Material.

§.§ Non-Linear Scenario I: Binary  Outcome and General Mediator

(1) Model.
Suppose the outcome is binary, and consider the model

\begin{align} \label{eq:extendmodel1}
P(\Y=1 \mid \S, \M, \X ) =  &~\mathrm{logit}^{-1}( \betaM \M + \X^\mytrans  \betaX + \DE \S), \\
\mathrm{E}(\M \mid \S, \X )=&~ h^{-1}(\alphaS \S +  \X^\mytrans \alphaX ), \notag
\end{align}

where $h^{-1}(\cdot)$ is the inverse of a canonical link function in generalized linear models.
Under Model (<ref>), since the outcome is binary, it is conventional to define the mediation effect as the odds ratio  [VanderWeele and Vansteelandt, 2010].
Specifically, under the identification assumption given in Section <ref>,  the conditional natural indirect effect (mediation effect) can be identified as

\begin{align*}
\mathrm{OR}^{\mathrm{NIE}}_{s\mid s^*}(s,\boldsymbol{x})=\frac{P\left\{Y(s, M(s))=1 \mid  \X=\boldsymbol{x} \right\} /\left\{1-P\left(Y(s, M(s))=1 \mid  \X=\boldsymbol{x}\right)\right\}}{P\left\{Y(s, M(s^*))=1 \mid  \X=\boldsymbol{x}\right\} /\left\{1-P\left\{Y(s, M(s^*))=1 \mid \X =\boldsymbol{x} \right\}\right\}},
\end{align*}

where $M(s)$ denotes the potential value of the mediator under the exposure $\S=s$,  and $Y(s,m)$ denotes the potential outcome that would have been observed if $\S$ and $\M$ had been set to $s$ and $m$, respectively.
Under $H_0$ of no mediation effect,

\begin{align} \label{eq:ornull}
H_0: \mathrm{OR}^{\mathrm{NIE}}_{s\mid s^*}(s,\boldsymbol{x})=1 \ \Leftrightarrow \ &~\log \mathrm{OR}^{\mathrm{NIE}}_{s\mid s^*}(s,\boldsymbol{x}) =0\\
\ \Leftrightarrow \ &~P\left\{Y(s, M(s))=1 \mid  \X=\boldsymbol{x} \right\}-P\left\{Y(s, M(s^*))=1 \mid  \X =\boldsymbol{x}\right\}=0, \notag
\end{align}

where the second equivalence follows from the strict increasing monotonicity of the function $x/(1-x)$ when $0<x<1$.

We consider natural indirect/mediation effects conditioning on covariates $\X=\boldsymbol{x}$ following [VanderWeele and Vansteelandt, 2010].
Alternatively, [Imai et al., 2010] proposed to examine the average NIE that marginalizes the distribution of $\X$.
Examining the conditional NIE is mainly for technical convenience.
The conditional NIE $=0$ for all $\boldsymbol{x}$ can give a sufficient condition for the average NIE $=0$.
Conclusions of conditional NIE may be obtained for average NIE similarly.
Please see Remark <ref> in the Supplementary Material for more details.

(2) Non-Regularity Issue.
The null hypothesis of no mediation effect (<ref>) looks different from $H_0:\alphaS\betaM=0$ under the linear SEMs in Section <ref>. Nevertheless, we can show that the non-regularity issue similar to that in Section <ref> would still arise.  This is formally stated as Proposition <ref> below.

Under the model (<ref>),    Condition <ref> (a general regularity condition on the link function $h^{-1}(\cdot)$ and the distribution of $M$), and identification conditions in Section <ref>,

* $H_0$ (<ref>) holds for $s\neq s^*$ if and only if $\alphaS=0$ or $\betaM=0$.

* For simplicity of notation, let $\mathrm{NIE}$ be a shorthand for $\log \mathrm{OR}^{\mathrm{NIE}}_{s\mid s^*}(s,\boldsymbol{x})$. We have

(i)  $\frac{\partial   \mathrm{NIE}}{\partial \alphaS} \Big|_{\betaM = 0}= \frac{\partial  \mathrm{NIE}}{\partial \betaM} \Big|_{\alphaS = 0}=0$,   (ii) $\frac{\partial  \mathrm{NIE}}{\partial \alphaS} \Big|_{\alphaS=0, \betaM\neq 0}\neq 0$,   (iii) $\frac{\partial  \mathrm{NIE}}{\partial \betaM} \Big|_{\alphaS\neq 0, \betaM = 0}\neq 0$.

It is interesting to see that even though the conditional  mediation effect $\mathrm{OR}^{\mathrm{NIE}}_{s\mid s^*}(s,\boldsymbol{x})$ does not take a product form,
a non-regularity issue caused by zero gradient can still arise under $H_0$ in (<ref>), which is similar to the PoC statistic in Section <ref>.
Proposition <ref> implies that
when $\alphaS=\betaM=0$,  the first-order Delta method cannot be directly applied to the inference of $\mathrm{NIE}$, which is different from the scenarios when $\alphaS\neq 0$ or $\betaM\neq 0$.
Therefore, we expect that the ordinary estimator of NIE can behave differently under different types of null hypotheses, and a non-regularity issue can occur.
This phenomenon is indeed demonstrated by numerical experiments in Section <ref> of the Supplementary Material.

(3) Asymptotic Theory and Adaptive Bootstrap.
For ease of presentation, we next derive  asymptotic theory under a special case of (<ref>), where the mediator is binary and follows a logistic regression model.
We point out that the analysis in this section can  be readily extended to cases where the mediator $M$ follows a linear model or other canonical generalized linear models.
let $\M$ and $\Y$ be Bernoulli random variables with mean values in (<ref>), and $h^{-1}(x)=\mathrm{logit}^{-1}(x)$.
In this case, $\log
\mathrm{OR}^{\text{NIE}}_{s\mid s^*}(s,\boldsymbol{x})=l(P_s)-l(P_{s^*}),
where  $P_s:=P\left\{Y(s, M(s))=1 \mid  \X=\boldsymbol{x} \right\} $,  $P_{s^*}:=P\left\{Y(s, M(s^*))=1 \mid  \X=\boldsymbol{x} \right\},
and $l(x)=\log\frac{x}{1-x}$.
Similarly to Section <ref>,
we are interested in understanding how the local limiting behaviors of $\alphaS$ and $\betaM$ coefficients change. To this end, we consider a general  local logistic model:

\begin{align} \label{eq:extendmodel1local}
\mathrm{E}(\M \mid \S, \X )=g(\alphaSn \S +  \X^\mytrans \alphaX),\hspace{1.5em}
\mathrm{E}(\Y \mid \S, \M, \X )  = g( \betaMn \M + \X^\mytrans  \betaX + \DE \S), %\notag
\end{align}

where $\alphaSn=\alphaS+\balpha/\sqrt{n}$, $\betaMn=\betaM+\bbeta/\sqrt{n}$, and $g(x)=\mathrm{logit}^{-1}(x)=e^x/(1+e^x)$.
Under the local model  (<ref>), we have for $\iota \in \{s,s^*\}$,

\begin{align}
P_{\iota}= &~ g\big( \iota \times \alphaSn   +\boldsymbol{x}^{\mytrans}\alphaX\big)\times d_{\beta,n} + P_*, \label{eq:pspsstardef}
\end{align}

where $ d_{\beta,n} =g(\betaMn  + \boldsymbol{x}^{\mytrans} \boldsymbol{\beta}_{\X} + \tau_S s ) -g( \boldsymbol{x}^{\mytrans} \boldsymbol{\beta}_{\X} + \tau_S s )$ and $P_*=g(\boldsymbol{x}^{\mytrans} \boldsymbol{\beta}_{\X} + \tau_S s )$. (Please see   the proof of  Theorem  <ref> for the derivations.)
For simplicity of notation, let NIE be a shorthand of $\log \mathrm{OR}^{\text{NIE}}_{s\mid s^*}(s, \boldsymbol{x})$, and  by (<ref>), $H_0\Leftrightarrow$ NIE $=0$.
Let $\widehat{\text{NIE}}=l(\hat{P}_s)-l(\hat{P}_{s^*})$ denote an estimator of NIE,
where $\hat{P}_s$ and $\hat{P}_{s^*}$ are defined similarly to (<ref>) with $({\alpha}_{S,n}, {\boldsymbol{\alpha}}_{\boldsymbol{X}}, {\beta}_{M,n}, {\boldsymbol{\beta}}_{\boldsymbol{X}}, {\tau}_S)$ replaced by their corresponding sample regression coefficient estimators $(\hat{\alpha}_{S}, \hat{\boldsymbol{\alpha}}_{\boldsymbol{X}}, \hat{\beta}_{M}, \hat{\boldsymbol{\beta}}_{\boldsymbol{X}}, \hat{\tau}_S)$.

Assume $P_s$ and $ P_{s^*}\in (0,1)$ and Condition <ref> in the Supplementary Material (a regularity condition on the design matrix similar to Condition <ref>).
Under the local model (<ref>) and $H_0:$ $\alphaS\betaM=0$,

(i) when  $(\alphaS,\betaM)\neq \mathbf{0}$, $
\sqrt{n}(\widehat{\mathrm{NIE}}-\mathrm{NIE})\xrightarrow{d}	 ( d_{\alpha} Z_{\beta} +   d_{\beta}  Z_{\alpha})\gamma_{0};

(ii) when  $(\alphaS,\betaM)=\mathbf{0}$,
n(\widehat{\mathrm{NIE}}-\mathrm{NIE})\xrightarrow{d}	( d_{\balpha} Z_{\beta} +   d_{\bbeta}Z_{\alpha}+ Z_{\alpha}Z_{\beta})\gamma_{0},

where $d_{\alpha}=g(\alphaS s + \boldsymbol{x}^{\mytrans}\alphaX)-g(\alphaSn s^* + \boldsymbol{x}^{\mytrans}\alphaX)$,  $ d_{\beta} =g(\betaM + \boldsymbol{x}^{\mytrans} \boldsymbol{\beta}_{\X} + \tau_S s ) -g(\boldsymbol{x}^{\mytrans} \boldsymbol{\beta}_{\X} + \tau_S s )$,
$d_{\balpha}=g'(\boldsymbol{x}^{\mytrans}\alphaX)(s-s^*) \balpha$, $
d_{\bbeta}=g'(\boldsymbol{x}^{\mytrans}\betaX + \tau_S s)(s-s^*) \bbeta$,
$(Z_{\alpha},Z_{\beta})$ represent bivariate mean-zero  normal distributions specified in Lemma <ref>, and $\gamma_0 = \{P_*(1-P_*)\}^{-1}$
is a non-zero constant with $P_*$ given in (<ref>).

We next study consistency of bootstrap estimators.
Let $\widehat{\mathrm{NIE}}^*$ denote the classical nonparametric bootstrap estimator  of $\mathrm{NIE}$.
Specifically, $\widehat{\mathrm{NIE}}^*=l( \hat{P}_s^*)-l(\hat{P}_{s^*}^*)$,
where $ \hat{P}_s^*$, and $\hat{P}_{s^*}^*$ are defined similarly to (<ref>) with $({\alpha}_{S,n}, {\boldsymbol{\alpha}}_{\boldsymbol{X}}, {\beta}_{M,n}, {\boldsymbol{\beta}}_{\boldsymbol{X}}, {\tau}_S)$ replaced by  their classical nonparametric bootstrap estimators
(\hat{\alpha}_S^*, \hat{\boldsymbol{\alpha}}_{\boldsymbol{X}}^*, \hat{\beta}_M^*, \hat{\boldsymbol{\beta}}_{\boldsymbol{X}}^*, \hat{\tau}_S^*).
Motivated by Theorem <ref>, we define the AB statistic

$$U_{e,1}^*= (\widehat{\mathrm{NIE}}^* -\widehat{\mathrm{NIE}} ) \times (1- \myindicator_{\alphaS, \lambda_n}^* \myindicator_{\betaM, \lambda_n}^*  )+n^{-1} \big( d_{\balpha}  \mathbb{Z}_{\beta}^* +   d_{\bbeta} \mathbb{Z}_{\alpha}^*+ \mathbb{Z}_{\alpha}\mathbb{Z}_{\beta}^*\big)\hat{\gamma}_{0}^* \times \myindicator_{\alphaS, \lambda_n}^*  \myindicator_{\betaM, \lambda_n}^*,$$

where $\myindicator_{\alphaS, \lambda_{n}}^*$ and $\myindicator_{\betaM, \lambda_{n}}^*$ are defined similarly to (<ref>).
The following theorem proves  consistency of the AB statistic $U_{e,1}^*$, based on which we can develop an AB test  similar to that in Section <ref>.

Under the conditions of Theorem <ref>, when the tuning parameter $\lambda_n$ satisfies $\lambda_n=o(\sqrt{n})$ and $\lambda_n\to \infty$ as $n\to \infty,$
$c_n U_{e,1}^* \overset{d^*}{\leadsto}c_n( \widehat{\mathrm{NIE}}- \mathrm{NIE}),$
where $c_n$ is specified as in (<ref>).

(4) Numerical Results.
We conduct simulation studies to compare the AB and the classical non-parametric bootstrap under the model (<ref>).
The detailed results are provided in Section <ref> of the Supplementary Material.

Our findings align closely with those presented
in Section <ref>.
Specifically, under $H_0:{\alpha}_S{\beta}_M=0$, when $({\alpha}_S, {\beta}_M)\neq \mathbf{0}$,
both the proposed AB test and the classical non-parametric bootstrap  yield uniformly distributed $p$-values.
However, when  $({\alpha}_S,{\beta}_M)= \mathbf{0}$, the classical bootstrap becomes overly conservative, whereas the AB test still yields uniformly distributed $p$-values.
Under $H_A$, the AB test can achieve higher empirical power than the classical bootstrap.

§.§ Non-Linear Scenario II: Linear Outcome and General Mediator

(1) Model.
Suppose the outcome follows a  linear model, and  consider

\begin{align} \label{eq:extendmodel2}
\mathrm{E}(\M \mid \S, \X )= h^{-1}(\alphaS \S +  \X^\mytrans \alphaX),  \hspace{1.5em} \mathrm{E}(\Y \mid \S, \M, \X )  =  \betaM \M + \X^\mytrans  \betaX + \DE \S,%\notag
\end{align}

where $h^{-1}(\cdot)$ can be the inverse of a canonical link function.
Similarly to the non-linear Scenario I, we examine the conditional natural indirect effect/mediation effect defined as the risk difference:

\begin{align}
\mathrm{NIE}_{s\mid s^*}(s,\boldsymbol{x}) := &~\mathrm{E}\big\{\Y(s, \M(s))- \Y(s, \M(s^*)) \mid \X =\boldsymbol{x} \big\} \notag\\
=&~\betaM\left\{ h^{-1}\big(\alphaS s+ \boldsymbol{x}^{\mytrans}\alphaX  \big) -h^{-1}\big(\alphaS s^*+ \boldsymbol{x}^{\mytrans}\alphaX \big) \right\}. \label{eq:h0binarymediator}
\end{align}

(2) Non-Regularity Issue.
We are interested in testing  $H_0: \mathrm{NIE}_{s\mid s^*}(s,\boldsymbol{x})=0$, which
looks different from $H_0:\alphaS\betaM=0$ in Section <ref>. Nevertheless, we can show that  the non-regularity issue similar to that in Section <ref> would arise.  This is formally stated as Proposition <ref> below.

Under the model (<ref>), assume $h^{-1}(\cdot)$ is strictly monotone, and the identification conditions in Section <ref> hold.
Let $\mathrm{NIE}$ be a shorthand for $\mathrm{NIE}_{s\mid s^*}(s,\boldsymbol{x})$ in (<ref>).

* $H_0:$ $\mathrm{NIE}=\eqref{eq:h0binarymediator}=0$ holds if and only if  $\alphaS=0$ or $\betaM=0$.

(i)  $\frac{\partial   \mathrm{NIE}}{\partial \alphaS} \Big|_{\betaM = 0}= \frac{\partial  \mathrm{NIE}}{\partial \betaM} \Big|_{\alphaS = 0}=0$.  (ii) $\frac{\partial  \mathrm{NIE}}{\partial \alphaS} \Big|_{\alphaS=0, \betaM\neq 0}\neq 0$.  (iii) $\frac{\partial  \mathrm{NIE}}{\partial \betaM} \Big|_{\alphaS\neq 0, \betaM = 0}\neq 0$.

Similarly to Proposition <ref>, Proposition <ref> implies that
a non-regularity issue caused by  zero gradient would arise  under $H_0$.
Specifically, the ordinary estimator of NIE can behave differently  when $\alphaS=\betaM=0$, and when one of $\alphaS$ and $\betaM \neq 0$.
This is similar to the PoC  statistic in Section <ref> and the odds ratio in Section <ref>.

(3) Asymptotic Theory and Adaptive Bootstrap.
For ease of presentation, we next derive  asymptotic theory under a specific instance of (<ref>).
Specifically,  the mediator $M$ is a Bernoulli random variable with its conditional mean given in (<ref>) and $h^{-1}(x)=\mathrm{logit}^{-1}(x)=e^x/(1+e^x)$,
and the outcome $Y$ follows the linear model in (<ref>).
The analysis in this section can be readily extended when the mediator $M$ follows other canonical generalized linear models.
As we are interested in how the local limiting behavior of $\alphaS$ and $\betaM$ coefficients change,
we consider the following general local  model

\begin{align} \label{eq:extendmodel2local}
\mathrm{E}(\M \mid \S, \X )=\mathrm{logit}^{-1}(\alphaSn \S +  \X^\mytrans \alphaX),\hspace{2em} Y  = 	\betaMn \M + \X^\mytrans  \betaX + \DE \S +\eY.
\end{align}

where $\alphaSn=\alphaS+\balpha/\sqrt{n}$, and $\betaMn=\betaM+\bbeta/\sqrt{n}$.

Assume Condition <ref> in the Supplementary Material (a regularity condition on the design matrix similar to Condition <ref>).
Under model (<ref>),

(i) when   $(\alphaS, \betaM)\neq \mathbf{0}$, $\sqrt{n}(\widehat{\mathrm{NIE}}-\mathrm{NIE}) \xrightarrow{d} d_{\alpha} Z_{\beta}+ \betaM Z_{\alpha} $;

(ii) when $(\alphaS, \betaM)=\mathbf{0}$, $n(\widehat{\mathrm{NIE}}-\mathrm{NIE})\xrightarrow{d}  d_{\balpha}  Z_{\beta} + \bbeta Z_{\alpha} +Z_{\alpha}Z_{\beta}$,

where $d_{\alphaS}=g(\alphaS s + \boldsymbol{x}^{\mytrans}\alphaX)-g(\alphaSn s^* + \boldsymbol{x}^{\mytrans}\alphaX)$,
$d_{\balpha}=g'(\boldsymbol{x}^{\mytrans}\alphaX)(s-s^*) \balpha$,
$Z_{\alpha}$ represents a normal distribution specified in Lemma <ref>,
and $Z_{\beta}$ is redefined to be a mean-zero normal distribution with a covariance same as  the random vector $V_M^{-1}\epsilon_YM_{\perp'}$, where $V_M$ and $M_{\perp'}$ are defined in Theorem <ref>.

We next establish bootstrap consistency theory.
Similarly to Section <ref>,
let $\widehat{\mathrm{NIE}}^*$ denote the nonparametric bootstrap estimator  of $\mathrm{NIE}$.
In particular, we redefine $\widehat{\mathrm{NIE}}^*=\hat{\beta}_M^*\{g(\hat{\alpha}_S^*s + \boldsymbol{x}^{\top}\hat{\boldsymbol{\alpha}}_{\boldsymbol{X}}^* ) -g(\hat{\alpha}_S^*s^* + \boldsymbol{x}^{\top}\hat{\boldsymbol{\alpha}}_{\boldsymbol{X}}^* )\}$,
$(\hat{\alpha}_S^*, \hat{\boldsymbol{\alpha}}_{\boldsymbol{X}}^*, \hat{\beta}_M^*)$ denotes the classical nonparametric bootstrap estimators of   $({\alpha}_S, {\boldsymbol{\alpha}}_{\boldsymbol{X}}, {\beta}_M)$.
Motivated by Theorem <ref>,  we define the AB statistic

U_{e,2}^*=(\widehat{\mathrm{NIE}}^* -\widehat{\mathrm{NIE}} )  (1- \myindicator_{\alphaS, \lambda_n}^* \myindicator_{\betaM, \lambda_n}^*  ) +n^{-1}(d_{\balpha}  \mathbb{Z}_{\beta}^* + \bbeta  \mathbb{Z}_{\alpha}^* +\mathbb{Z}_{\alpha}^*\mathbb{Z}_{\beta}^*)  \myindicator_{\alphaS, \lambda_n}^*  \myindicator_{\betaM, \lambda_n}^*,

where $\myindicator_{\alphaS, \lambda_{n}}^*$ and $\myindicator_{\betaM, \lambda_{n}}^*$ are defined similarly to (<ref>).
The following theorem establishes consistency of the AB statistic $U_{e,2}^*$.

Under conditions of Theorem <ref>, when the tuning parameter $\lambda_n$ satisfies $\lambda_n=o(\sqrt{n})$ and $\lambda_n\to \infty$ as $n\to \infty,$
$c_nU_{e,2}^*\overset{d^*}{\leadsto} c_{n}( \widehat{\mathrm{NIE}}- \mathrm{NIE}),$
where $c_n$ is specified as in (<ref>).

(4) Numerical Results.
We conduct simulation studies to compare the AB and the classical non-parametric bootstrap under the model (<ref>).
The detailed results are provided in Section <ref> of the Supplementary Material.
The obtained results are very similar to those in Section <ref> and  Section <ref> Part 4),
and therefore, we refrain from repeating the details here.

§ DATA ANALYSIS

We illustrate an application of our proposed method to the analysis of data from a cohort study “Early Life Exposures in Mexico to ENvironmental Toxicants" (ELEMENT) [Perng et al., 2019].
One of the central interests in this scientific study concerns the mediation effects of metabolites, in particular, the family of lipids,
on the association between environmental exposure and children growth and development. In the literature of environmental health sciences, exposure to endocrine-disrupting chemicals (EDCs) such as phthalates have been found to be detrimental to children's health outcomes. Such findings of  direct associations need to be further assessed for possible mediation effects through metabolites, because environmental toxicants such as phthalates can alter metabolic profiles at the molecular level.

Our illustration focuses on the outcome of body mass index (BMI) and exposure to one phthalate, MEOHP (a chemical in food production and storage).  BMI is a widely used biomarker in pediatric research to measure childhood obesity.
The dataset contains 382   adolescents aged 10-18 years old living in Mexico   City.
Our mediation analysis involves a set of 149 lipids that are hypothesized to have potential mediation effects on children's growth and development.  Our goal is to identify the mediation pathways of exposure to MEOHP $\rightarrow$ lipids $\rightarrow$ BMI.
Two key potential confounders, gender and age,  are included throughout the analyses.

It is worth noting that adjusting for gender and age may not be sufficient for proper confounding adjustments. To conduct a more plausible causal analysis and interpretation, a further investigation is deemed necessary to rigorously assess the underlying causal assumptions such as a sensitivity analysis for the sequential ignorability assumption.
In our analyses,
we compare the results of six tests:
JS-AB, JS-MaxP, PoC-AB, PoC-B, PoC-Sobel, and CMA, which have been compared in our simulation studies  in Section <ref>.
In particular,
all the bootstrap methods (including JS-AB, PoC-AB, PoC-B, CMA)
are conducted based on $10^4$ bootstrap resamples.
Here we no longer include the JS-B test
and the MT-Comp method, as they are known to have  inflated type 1 errors according to our simulation studies in Section <ref>.

The joint model generates $p$-values of 15 top lipids selected by a preliminary  screening step.

Lipid     JS-AB      JS-MaxP      PoC-AB       PoC-B     PoC-Sobel      CMA
X16.0.LYSO.PC_3     0.2897     0.5402     0.4195     0.5433     0.6295     0.5680

FA.7.0-OH_1     0.3755     0.6095     0.1156     0.6146     0.6179     0.5820

FA.7.0-OH_2     0.0654     0.2515     0.1148     0.3921     0.4057     0.4440

LPC.18.3_1     0.4730     0.6815     0.5068     0.6706     0.7047     0.6620

GLY     0.1824     0.4290     0.1405     0.5822     0.5611     0.5980

GLY-H2O     0.5448     0.7465     0.2202     0.7753     0.7537     0.8020

LAURIC.ACID     0.0042     0.0623     0.0027     0.0512     0.1403     0.0560

X16.0.LYSO.PC_2     0.9836     0.9966     0.6871     0.8015     0.9966     0.8600

FA.16.0-OH_1     0.4795     0.6900     0.2792     0.6966     0.6992     0.6840

LPC.16.1_3     0.1738     0.4156     0.1907     0.5553     0.4922     0.5740

LPC.18.2_2     0.1732     0.4105     0.1793     0.5009     0.4842     0.4720

URSODIOL_2     0.6889     0.8278     0.7346     0.9274     0.8388     0.9460

DECENOYL     0.0882     0.2916     0.0589     0.3891     0.3269     0.3940

FA.12.0-DiC     0.5502     0.7416     0.5294     0.8026     0.7545     0.7900

FA.18.0-DiC     0.0095     0.0963     0.0158     0.1462     0.1899     0.1740

Abbreviations:  dicarboxylic (DiC); fatty acid (FA); hydroxyl (OH); lysophosphatidylcholine (LPC);  glycochenodeoxycholate (GLY); decenoylcarnitine (DECENOYL).

As the sample size is limited compared to the large number of mediators, we first apply
a screening analysis to identify a subset of lipids as potential candidates.
We then jointly model the chosen lipids in the second step of our analysis.
To mitigate the potential issues arising from double dipping the data, we adopt a random data splitting approach by dividing the dataset into two distinct parts, each dedicated to one of the two respective analytic tasks.

In the first screening step, we examine the effect along the path MEOHP $\to$ lipid $\to$ BMI for one lipid at a time, and the corresponding  $p$-values are obtained with the six tests, respectively.

For each test, we select
a proportion ($q\%$) of lipids with the smallest $p$-values.
The second step
examines the path MEOHP $\to$ selected lipids $\to$ BMI, with the selected lipids being modeled jointly.

To test the mediation effect through a target lipid $M$ within the selected set,
we adjust for non-target mediators within the outcome model, following the discussions on Page  page:discussionmultimed; please see
more details in  Section <ref> of the Supplementary Material.

Subsequently, we select lipids based on their $p$-values obtained in the second step, after adjusting for multiple comparisons with controlled false discovery rate (FDR)  [Benjamini and Hochberg, 1995].
In our analysis, we explore a range of $q$ values $\{5, 10, 15, 20, 25\}$ and observe very similar results, indicating the robustness of our approach to the choice of the screening threshold in the first step.
We next present the results obtained with $q = 10$ (i.e., 15 selected lipids based on their $p$-values), while results for other $q$ values are detailed in Section <ref>  of the Supplementary Material.

As an illustrative example,
we first present the results from a single random split in Table <ref>.
Table <ref> provides the corresponding $p$-values for the lipids selected by at least one test in the second step of the analysis.

In this instance, the non-AB tests fail to detect any lipids. In contrast, the PoC-AB test identifies lauric acid (L.A) and FA.7.0-OH_1 (FA.7) while controlling the FDR at 0.10, and the JS-AB test selects both L.A and FA.7 when the FDR is controlled at 0.05 and 0.10, respectively.

To gauge the variability of results across random splits, we repeat the data-splitting analysis 400 times.

As shown in Figure <ref>, L.A and  FA.7 are the two most frequently selected mediators in our analysis.

Furthermore, the  AB tests exhibit a notably higher chance of selecting L.A compared to the non-AB tests.
This aligns with our observations from simulations in Section <ref>,  suggesting that the AB tests can attain higher power than their non-AB counterparts.
Lauric acid is a saturated fatty acid and is found in many vegetable fats and in coconut and palm kernel oils [Dayrit, 2015].
The results suggest
that the exposure to MEOHP may influence the  process of breaking down fat tissue in the human body, leading to obesity and other adverse health outcomes.

Lipids selected in the second step.

Lipids     JS-AB      JS-MaxP      PoC-AB       PoC-B     PoC-Sobel      CMA
L.A     0.0017 ($\times *$)     0.0399     0.0043 ($*$)     0.0406     0.1254     0.0426

FA.7     0.0008 ($\times *$)      0.0146     0.0090 $(*)$     0.0236     0.0937     0.0208

(Abbreviations: LAURIC.ACID (L.A); FA.7.0-OH_1 (FA.7). $p$-values with $(\times)$ and $(*)$ indicate that the lipid specified by the row is selected by the method specified by the column under 0.05 and 0.10 FDR levels, respectively.)

Since the first screening step considers one mediator at a time,  we also conduct sensitivity analyses to evaluate the effects of the unadjusted mediators  similarly to [Liu et al., 2021].
We use the procedure  proposed by [Imai et al., 2010], which  utilized the idea that  the error term in the M-S model and that in the Y-M model are likely to be correlated if the sequential ignorability assumption is violated and vice versa.
The detailed results are provided in Section <ref> in the Supplementary Material.
As a brief summary, the sensitivity analysis suggests that our  first  screening analysis could be robust to unadjusted  mediators.

Times of mediators being selected in Step 2 by the six tests with  FDR$=0.10$ over 400 random splits of the data.

§ DISCUSSION

This paper proposes a new adaptive framework  for testing composite null hypotheses in mediation pathway analysis.
The method incorporates a consistent pre-test threshold into the bootstrap procedure, which helps  circumvent the non-regularity issue arising from the composite null hypotheses.
If at least one of the two coefficients is significant, the procedure would reduce to the classical nonparametric bootstrap; otherwise, it approximates the local asymptotic behavior of the statistics.
Our proposed strategy accommodates different types of null hypotheses under various models.
Particularly, we have established similar results for both the individual and joint mediation effects under classical linear structural equation models,
and we have generalized the conclusions under generalized linear models.
Through comprehensive simulation studies,
we have demonstrated that the adaptive tests can properly and robustly control the type 1 error under different types of null hypotheses and improve the statistical power.

The proposed methodology offers an exemplary analytic toolbox that can be broadly extended to handle other problems of similar types involving composite null hypotheses.
There are several interesting future research directions that are worth exploration.
First, the non-regularity issue can similarly arise in other scenarios, such as  survival analysis  [VanderWeele, 2011, Huang and Pan, 2016],  different data types [Sohn and Li, 2019], partially linear models  [Hines et al., 2021], and models with exposure-mediator interactions; see more discussions in Section <ref> of the Supplementary Material.

These complicated models
require special care in
the causal interpretation of mediation effects and in the implementation of the bootstrap procedure, warranting further investigation.

when the dimension of mediators and covariates becomes high,
it is of interest to extend the adaptive bootstrap under high-dimensional mediation models for both individual and joint mediation effects
[Zhou et al., 2020].
Similarly to our discussions on adjusting multivariate mediators at the end of Section <ref>,
we might apply the adaptive bootstrap after properly adjusting high-dimensional covariates.
In the data analysis, we have applied the marginal screening to reduce the dimension of mediators, which might potentially overlook the complicated causal dependence among mediators.
When mediators have potential  causal dependence,
[Shi and Li, 2022] proposed to
first estimate a directed acyclic graph of mediators and develop a testing procedure that can control the type 1 error to be less than or equal to the nominal level.
It would be of interest to extend our proposed AB under such settings to mitigate potential conservatism.
Third, the proposed AB strategy can also be utilized to examine the replicability across independent studies [Bogomolov and Heller, 2018], which is fundamental to scientific discovery.
Specifically, let $\beta_i, i=1,\ldots, K$, denote the true signals from $K$ independent studies, respectively.
Testing whether the signals in these $K$ studies are all significant  corresponds to   $H_0: \prod_{i=1}^K \beta_i=0$ versus $H_A: \prod_{i=1}^K \beta_i\neq 0$.
Moreover, for two studies with true signals $\beta_1$ and $\beta_2$,
to investigate whether the effects of both studies are  significant in the same direction,
one can formulate
the hypothesis testing problem as $H_0: \beta_1\beta_2 \leq 0$ versus $H_A: \beta_1\beta_2>0$.
For these testing problems,  the null hypotheses are composite.

To  properly control the type 1 error,
the adaptive strategy proposed in this paper may serve as a valuable building block, while additional effort is needed to analyze those different cases carefully.
Last, in our data analysis, all measurements are obtained cross-sectionally at one given clinical visit within a time window of approximately three months.
To further study potential long-term influences of toxicant exposures, it may be of interest to
investigate how the mediation effects might vary over time.
Such time-varying mediation effects may be naturally analyzed in the scenario of longitudinal studies that collect time-varying measurements.
This is a very challenging research field with only minimal investigation in the current literature [Bind et al., 2016].
Extending the proposed AB method to analyze time-varying mediation effects would be a compelling future direction.

§ ACKNOWLEDGMENTS

We are grateful to the joint editors, Dr. Daniela Witten and Dr. Aurore Delaigle, an associate editor, and three anonymous referees for their helpful comments and suggestions.
This work is partially supported by  NSF DMS-1811734, DMS-2113564,  SES-1846747, SES-2150601,  NIH R01ES024732,
R01ES033656, and Wisconsin Alumni Research Foundation.

Conflict of interest: None declared.

§ SUPPLEMENTARY MATERIAL

We provide more results in the online Supplementary Material, including all the proofs, hyperparameter tuning, additional results on simulations and data analysis, and details of extensions (Joint Significance test, multiple-mediator settings, and non-linear models).
For reproducibility of our results, we provide the R package and code on the GitHub repository:  [He et al., 2023].

§ DATA AVAILABILITY

Due to privacy restrictions, we are unable to directly share the raw data publicly but they may be obtained offline according to a formal data request procedure outlined in the University of Michigan Data Use Agreement protocol. To satisfy the need of reproducibility, instead, we have introduced a pseudo-dataset with  added noise  on the GitHub repository: [He et al., 2023].

[Andrews, 2001]
Andrews, D. W. (2001) Testing when a parameter is on the boundary of the
maintained hypothesis.
Econometrica, 69, 683–734.

[Barfield et al., 2017]
Barfield, R., Shen, J., Just, A. C., Vokonas, P. S., Schwartz, J., Baccarelli,
A. A., VanderWeele, T. J. and Lin, X. (2017) Testing for the indirect effect
under the null for genome-wide mediation analyses.
Genetic Epidemiology, 41, 824–833.

[Baron and Kenny, 1986]
Baron, R. M. and Kenny, D. A. (1986) The moderator–mediator variable
distinction in social psychological research: Conceptual, strategic, and
statistical considerations.
Journal of Personality and Social Psychology, 51,

[Benjamini and Hochberg, 1995]
Benjamini, Y. and Hochberg, Y. (1995) Controlling the false discovery rate: a
practical and powerful approach to multiple testing.
Journal of the Royal Statistical Society: Series B
(Methodological), 57, 289–300.

[Bickel and Freedman, 1981]
Bickel, P. J. and Freedman, D. A. (1981) Some asymptotic theory for the
The annals of statistics, 9, 1196–1217.

[Bind et al., 2016]
Bind, M.-A., Vanderweele, T., Coull, B. and Schwartz, J. (2016) Causal
mediation analysis for longitudinal data with exogenous exposure.
Biostatistics, 17, 122–134.

[Blum et al., 2020]
Blum, M. G., Valeri, L., François, O., Cadiou, S., Siroux, V., Lepeule,
J. and Slama, R. (2020) Challenges raised by mediation analysis in a
high-dimension setting.
Environmental Health Perspectives, 128, 055001.

[Bogomolov and Heller, 2018]
Bogomolov, M. and Heller, R. (2018) Assessing replicability of findings across
two studies of multiple features.
Biometrika, 105, 505–516.

[Chen, 2016]
Chen, S. X. (2016) Peter Hall's contributions to the bootstrap.
The Annals of Statistics, 44, 1821–1836.

[Dai et al., 2020]
Dai, J. Y., Stanford, J. L. and LeBlanc, M. (2020) A multiple-testing procedure
for high-dimensional mediation hypotheses.
Journal of the American Statistical Association, 1–16.

[Daniel et al., 2015]
Daniel, R., De Stavola, B., Cousens, S. and Vansteelandt, S. (2015) Causal
mediation analysis with multiple mediators.
Biometrics, 71, 1–14.

[Dayrit, 2015]
Dayrit, F. M. (2015) The properties of lauric acid and their significance in
coconut oil.
Journal of the American Oil Chemists' Society, 92,

[Derkach et al., 2020]
Derkach, A., Moore, S. C., Boca, S. M. and Sampson, J. N. (2020) Group testing
in mediation analysis.
Statistics in Medicine, 39, 2423–2436.

[Djordjilović et al., 2020]
Djordjilović, V., Hemerik, J. and Thoresen, M. (2020) On optimal two-stage
testing of multiple mediators.
arXiv preprint arXiv:2007.02844.

[Djordjilović et al., 2019]
Djordjilović, V., Page, C. M., Gran, J. M., Nøst, T. H., Sandanger,
T. M., Veierød, M. B. and Thoresen, M. (2019) Global test for
high-dimensional mediation: Testing groups of potential mediators.
Statistics in Medicine, 38, 3346–3360.

[Drton and Xiao, 2016]
Drton, M. and Xiao, H. (2016) Wald tests of singular hypotheses.
Bernoulli, 22, 38–59.

[Du et al., 2022]
Du, J., Zhou, X., Hao, W., Liu, Y., Jennifer, S. and Mukherjee, B. (2022)
Methods for large-scale single mediator hypothesis testing: Possible choices
and comparisons.
arXiv preprint arXiv:2203.13293.

[Dufour et al., 2013]
Dufour, J.-M., Renault, E. and Zinde-Walsh, V. (2013) Wald tests when
restrictions are locally singular.
arXiv preprint arXiv:1312.0569.

[Frisch and Waugh, 1933]
Frisch, R. and Waugh, F. V. (1933) Partial time regressions as compared with
individual trends.
Econometrica, 387–401.

[Fritz and MacKinnon, 2007]
Fritz, M. S. and MacKinnon, D. P. (2007) Required sample size to detect the
mediated effect.
Psychological Science, 18, 233–239.

[Fulcher et al., 2019]
Fulcher, I. R., Shi, X. and Tchetgen, E. J. T. (2019) Estimation of natural
indirect effects robust to unmeasured confounding and mediator measurement
Epidemiology, 30, 825–834.

[Glonek, 1993]
Glonek, G. (1993) On the behaviour of wald statistics for the disjunction of
two regular hypotheses.
Journal of the Royal Statistical Society: Series B
(Methodological), 55, 749–755.

[Guo et al., 2022]
Guo, X., Li, R., Liu, J. and Zeng, M. (2022) High-dimensional mediation
analysis for selecting DNA methylation loci mediating childhood trauma and
cortisol stress reactivity.
Journal of the American Statistical Association, 1–32.

[Hao and Song, 2022]
Hao, W. and Song, P. X.-K. (2022) A simultaneous likelihood test for joint
mediation effects of multiple mediators.
Statistica Sinica.

[Hayes, 2017]
Hayes, A. F. (2017) Introduction to mediation, moderation, and
conditional process analysis: A regression-based approach.
Guilford publications.

[He et al., 2023]
He, Y., Song, P. X.-K. and Xu, G. (2023) ABtest.

[Hines et al., 2021]
Hines, O., Vansteelandt, S. and Diaz-Ordaz, K. (2021) Robust inference for
mediated effects in partially linear models.
Psychometrika, 1–24.

[Huang, 2018]
Huang, Y.-T. (2018) Joint significance tests for mediation effects of
socioeconomic adversity on adiposity via epigenetics.
The Annals of Applied Statistics, 12, 1535–1557.

[Huang, 2019]
— (2019a) Genome-wide analyses of sparse mediation effects under
composite null hypotheses.
The Annals of Applied Statistics, 13, 60–84.

[Huang, 2019]
— (2019b) Variance component tests of multivariate mediation
effects under composite null hypotheses.
Biometrics, 75, 1191–1204.

[Huang and Cai, 2016]
Huang, Y.-T. and Cai, T. (2016) Mediation analysis for survival data using
semiparametric probit models.
Biometrics, 72, 563–574.

[Huang and Pan, 2016]
Huang, Y.-T. and Pan, W.-C. (2016) Hypothesis test of mediation effect in
causal mediation model with high-dimensional continuous mediators.
Biometrics, 72, 402–413.

[Imai et al., 2010]
Imai, K., Keele, L. and Tingley, D. (2010a) A general approach to
causal mediation analysis.
Psychological Methods, 15, 309–334.

[Imai et al., 2010]
Imai, K., Keele, L. and Yamamoto, T. (2010b) Identification,
inference and sensitivity analysis for causal mediation effects.
Statistical Science, 25, 51–71.

[Imai and Yamamoto, 2013]
Imai, K. and Yamamoto, T. (2013) Identification and sensitivity analysis for
multiple causal mechanisms: Revisiting evidence from framing experiments.
Political Analysis, 21, 141–171.

[Imbens and Rubin, 2015]
Imbens, G. W. and Rubin, D. B. (2015) Causal inference in statistics,
social, and biomedical sciences.
Cambridge University Press.

[Jérolon et al., 2020]
Jérolon, A., Baglietto, L., Birmelé, E., Alarcon, F. and Perduca, V.
(2020) Causal mediation analysis in presence of multiple mediators uncausally
The International Journal of Biostatistics, 17,

[Judd and Kenny, 1981]
Judd, C. M. and Kenny, D. A. (1981) Estimating the effects of social
Cambridge University Press.

[Laber and Murphy, 2011]
Laber, E. B. and Murphy, S. A. (2011) Adaptive confidence intervals for the
test error in classification.
Journal of the American Statistical Association,
106, 904–913.

[Liu et al., 2021]
Liu, Z., Shen, J., Barfield, R., Schwartz, J., Baccarelli, A. A. and Lin, X.
(2021) Large-scale hypothesis testing for causal mediation effects with
applications in genome-wide epigenetic studies.
Journal of the American Statistical Association, 1–15.

[Loh et al., 2021]
Loh, W. W., Moerkerke, B., Loeys, T. and Vansteelandt, S. (2021) Disentangling
indirect effects through multiple mediators without assuming any causal
structure among the mediators.
Psychological Methods.

[MacKinnon, 2008]
MacKinnon, D. (2008) Introduction to Statistical Mediation Analysis.
Multivariate Applications Book Series. Taylor & Francis.

[MacKinnon and Fairchild, 2009]
MacKinnon, D. P. and Fairchild, A. J. (2009) Current directions in mediation
Current Directions in Psychological Science, 18,

[MacKinnon et al., 2002]
MacKinnon, D. P., Lockwood, C. M., Hoffman, J. M., West, S. G. and Sheets, V.
(2002) A comparison of methods to test mediation and other intervening
variable effects.
Psychological Methods, 7, 83–104.

[MacKinnon et al., 2004]
MacKinnon, D. P., Lockwood, C. M. and Williams, J. (2004) Confidence limits for
the indirect effect: Distribution of the product and resampling methods.
Multivariate Behavioral Research, 39, 99–128.

[McKeague and Qian, 2015]
McKeague, I. W. and Qian, M. (2015) An adaptive resampling test for detecting
the presence of significant predictors.
Journal of the American Statistical Association,
110, 1422–1433.

[McKeague and Qian, 2019]
— (2019) Marginal screening of 2$\times$2 tables in large-scale case-control
Biometrics, 75, 163–171.

[Miles and Chambaz, 2021]
Miles, C. H. and Chambaz, A. (2021) Optimal tests of the composite null
hypothesis arising in mediation analysis.
arXiv preprint arXiv:2107.07575.

[Pearl, 2001]
Pearl, J. (2001) Direct and indirect effects.
In Proceedings of the Seventeenth Conference on Uncertainty
in Artificial Intelligence, UAI’01, 411–420. San Francisco, CA, USA:
Morgan Kaufmann Publishers Inc.

[Perng et al., 2019]
Perng, W. et al. (2019) The early life exposure in mexico to environmental
toxicants (ELEMENT) project.
British Medical Journal Open, 9, e030427.

[Robins and Greenland, 1992]
Robins, J. M. and Greenland, S. (1992) Identifiability and exchangeability for
direct and indirect effects.
Epidemiology, 3, 143–155.

[Rubin, 1980]
Rubin, D. B. (1980) Randomization analysis of experimental data: The fisher
randomization test comment.
Journal of the American statistical association,
75, 591–593.

[Sampson et al., 2018]
Sampson, J. N., Boca, S. M., Moore, S. C. and Heller, R. (2018) FWER and FDR
control when testing multiple mediators.
Bioinformatics, 34, 2418–2424.

[Shi and Li, 2022]
Shi, C. and Li, L. (2022) Testing mediation effects using logic of boolean
Journal of the American Statistical Association,
117, 2014–2027.

[Sobel, 1982]
Sobel, M. E. (1982) Asymptotic confidence intervals for indirect effects in
structural equation models.
Sociological Methodology, 13, 290–312.

[Sohn and Li, 2019]
Sohn, M. B. and Li, H. (2019) Compositional mediation analysis for microbiome
The Annals of Applied Statistics, 13, 661–681.

[Tingley et al., 2014]
Tingley, D., Yamamoto, T., Hirose, K., Keele, L. and Imai, K. (2014) Mediation:
R package for causal mediation analysis.
Journal of Statistical Software, 59, 1–38.

[van der Vaart, 2000]
van der Vaart, A. W. (2000) Asymptotic statistics, vol. 3 of
Cambridge Series in Statistical and Probabilistic Mathematics.
Cambridge University Press.

[Valeri and VanderWeele, 2013]
Valeri, L. and VanderWeele, T. J. (2013) Mediation analysis allowing for
exposure–mediator interactions and causal interpretation: theoretical
assumptions and implementation with SAS and SPSS macros.
Psychological Methods, 18, 137–155.

[Van Garderen and Van Giersbergen, 2022]
Van Garderen, K. and Van Giersbergen, N. (2022) A nearly similar powerful test
for mediation.
Tech. rep.

[VanderWeele, 2011]
VanderWeele, T. J. (2011) Causal mediation analysis with survival data.
Epidemiology, 22, 582–585.

[VanderWeele, 2015]
— (2015) Explanation in causal inference: methods for mediation and
Oxford University Press.

[VanderWeele and Vansteelandt, 2009]
VanderWeele, T. J. and Vansteelandt, S. (2009) Conceptual issues concerning
mediation, interventions and composition.
Statistics and its Interface, 2, 457–468.

[VanderWeele and Vansteelandt, 2010]
— (2010) Odds ratios for mediation analysis for a dichotomous outcome.
American Journal of Epidemiology, 172, 1339–1348.

[VanderWeele and Vansteelandt, 2014]
— (2014) Mediation analysis with multiple mediators.
Epidemiologic Methods, 2, 95–115.

[VanderWeele et al., 2014]
VanderWeele, T. J., Vansteelandt, S. and Robins, J. M. (2014) Effect
decomposition in the presence of an exposure-induced mediator-outcome
Epidemiology (Cambridge, Mass.), 25, 300.

[Vansteelandt and Daniel, 2017]
Vansteelandt, S. and Daniel, R. M. (2017) Interventional effects for mediation
analysis with multiple mediators.
Epidemiology (Cambridge, Mass.), 28, 258.

[Wang et al., 2018]
Wang, H. J., McKeague, I. W. and Qian, M. (2018) Testing for marginal linear
effects in quantile regression.
Journal of the Royal Statistical Society: Series B
(Statistical Methodology), 80, 433–452.

[Zhao et al., 2014]
Zhao, S. D., Cai, T. T. and Li, H. (2014) More powerful genetic association
testing via a new statistical framework for integrative genomics.
Biometrics, 70, 881–890.

[Zhou et al., 2020]
Zhou, R. R., Wang, L. and Zhao, S. D. (2020) Estimation and inference for the
indirect effect in high-dimensional linear mediation models.
Biometrika, 107, 573–589.

Supplementary Material for
“Adaptive Bootstrap Tests for Composite Null Hypotheses in the Mediation Pathway Analysis”

In this Supplementary Material,
Section <ref> provides rigorous definitions of notation used in the main text and proofs.
Section <ref> presents details of the Adaptive Bootstrap for the Joint Significance Test, mentioned on Page page:adjst of the main text.

Section <ref> provides proofs of  Theorems <ref>–<ref> and Theorems
Section <ref> provides detailed theoretical and numerical results on multiple mediators  supplementary to Figure <ref> and Section <ref> in the main text.
Section <ref>  provides detailed theoretical and numerical results on non-linear models supplementary to Sections <ref> and <ref>.
Section <ref> discuss implementation details of a data-driven procedure for choosing the tuning parameter.
Section <ref> provides additional numerical experiments and data analysis results.
Section <ref> provides a clarification of the partially linear model mentioned in Section <ref> of the main text.

§ DEFINITIONS AND NOTATION

Notation on Convergence
For two sequences of real numbers $\{a_n\}$ and $\{b_n\}$,
we write $a_n=o(b_n)$ if $\lim_{n\to \infty} a_n/b_n=0$.
We say $a_n\to \infty$ as $n\to \infty$ if the value of $a_n$ can be arbitrarily large by taking $n$ to be sufficiently large.
Given a sequence of random variables $\{X_n\}$ and a random variable $X$,
we let $X_n \xrightarrow{\mathrm{P}} X$ represent  the convergence in probability, i.e.,
for any $\epsilon>0$, $\lim_{n\to \infty} \mathrm{P}(|X_n-X|>\epsilon)=0$.
We let $X_n \xrightarrow{a.s.} X$
represent the almost sure convergence, i.e., $ \mathrm{P}(\, \lim_{n\to \infty} X_n = X )=1$.
We let $X_n \xrightarrow{d} X$
represent the convergence in distribution, i.e.,
$\lim_{n\to \infty} F_n(x) = F(x)$ for every $x$ at which $F(x)$ is continuous,
where $F_n(x)$ and $F(x)$ represent the cumulative distribution functions of $X_n$ and $X$, respectively.

Bootstrap Consistency
We next introduce
the definition of bootstrap consistency relative to the Kolmogorov-Smirnov distance; also see,  Section 23.2 of [van der Vaart, 2000].
Let $\hat{\theta}_n$ be an estimator of some parameter $\theta$ attached to the distribution $P_n$ of the observations.
Let $\mathbb{P}_n$ be an estimate of the underlying distribution $P_n$ of the observations.
Let $\hat{\theta}_n^*$  denote the bootstrap estimator for  $\hat{\theta}_n$ and are obtained according to $\mathbb{P}_n$ in the same way $\hat{\theta}_n$ computed from the true observations with distribution $P_n$.
We write
$\sqrt{n}(\hat{\theta}_n^*-\hat{\theta}_n) \overset{d^*}{\leadsto} \sqrt{n}(\hat{\theta}_n - \theta)$ if

\begin{align}\label{eq:bootconsistdef}
\sup_{x\in \mathbb{R}} \left| P\left( \sqrt{n}(\hat{\theta}_n - \theta) \leqslant x \mid P_n \right) - P\left( \sqrt{n}(\hat{\theta}_n^*-\hat{\theta}_n)\leqslant x  \mid \mathbb{P}_n\right)  \right| \xrightarrow{\mathrm{P}} 0.
\end{align}

In the following proofs,
when the target $\sqrt{n}(\hat{\theta}_n - \theta)$ converges in distribution to a continuous distribution function $F$,
(<ref>) is also equivalent to that for every $x\in \mathbb{R}$,

\begin{align*}
\text{if}\ \    P\left( \sqrt{n}(\hat{\theta}_n - \theta) \leqslant x \mid P_n \right) \to F(x),\quad\  \text{then}\quad P\left( \sqrt{n}(\hat{\theta}_n^*-\hat{\theta}_n)\leqslant x \mid \mathbb{P}_n \right) \xrightarrow{\mathrm{P}} F(x).
\end{align*}

This type of consistency implies the asymptotic consistency of confidence intervals.
Moreover, we let $\hat{\theta}_n^* \overset{\mathrm{P}^*}{\leadsto} \theta$ denote the convergence in conditional probability, i.e., for any $\epsilon>0$,

\begin{align*}
P(|\hat{\theta}_n^* - {\theta}|>\epsilon \mid \mathbb{P}_n ) \xrightarrow{a.s.} 0.
\end{align*}

§ ADAPTIVE BOOTSTRAP  FOR THE JOINT SIGNIFICANCE TEST

In addition to the PoC test, another popular class of methods is the joint significance test [MacKinnon et al., 2002],
which is useful for combining a series of tests for a set of links in a causal chain [Judd and Kenny, 1981, Baron and Kenny, 1986] such as a directed acyclic graph.
Specifically, the JS test, also known as the MaxP test, rejects $H_0: \alphaS \betaM = 0$ if   $p_{\text{MaxP}}\equiv \max\{p_{\alphaS}, p_{\betaM}\} < \omega$,
where $\omega$ is a prespecified significance level,
and $p_{\alphaS}$  and $p_{\betaM}$ denote the $p$-values for $H_0: \alphaS=0$ (the link $\S \rightarrow \M$)  and $H_0: \betaM=0$ (the link $\M \rightarrow \Y$), respectively.
Despite the ease of operation,  the MaxP test has also been found to be overly conservative under  $H_{0,3}$    [Barfield et al., 2017].
To see this analytically,
note that when $p_{\alphaS}$ and $p_{\betaM}$ are asymptotically independent,
under $H_{0,3}$, $\Pr(p_{\text{MaxP}}<\omega) \to  \Pr(p_{\alphaS}<\omega) \Pr(p_{\betaM}<\omega) \to \omega^2 <\omega$,
which suggests that  the MaxP test is always conservative under $H_{0,3}$ even if the sample size goes to infinity.

In this subsection, we focus on the adaptive bootstrap for the JS test.

As  discussed in Section <ref>, the popular MaxP test that rejects $H_0: \alphaS \betaM =0$ using the rule $p_{\mathrm{MaxP}}<\omega$
can be conservative.

To correctly evaluate the distribution of  $p_{\mathrm{MaxP}}$ under the composite null,
we next develop the corresponding adaptive bootstrap test procedure for the JS test.

For ease of deriving bootstrap consistency,
instead of directly examining  $p_{\mathrm{MaxP}}$,

we investigate the corresponding statistic based on the t-statistics, which usually have asymptotic normality.

Specifically, we introduce the statistic $\sqrt{n}\hatmaxpn = H(T_{\alpha, n}, T_{\beta, n})$,
where $T_{\alpha,n}=\sqrt{n}\hatalphaSn/\hatsigmaalphan $ and $ T_{\beta, n}=\sqrt{n}\hatbetaMn /\hatsigmabetan$ are the standardized statistics of the two coefficients, respectively, and

\begin{align*}
H(t_1,t_2)=(t_1, t_2) h(t_1, t_2)\   \text{ with }\ 	h(t_1, t_2)=\biggr( \myindicator\biggr\{\underset{k=1,2}{\operatorname{min}}\,|t_k|=|t_1|\biggr\}, ~\myindicator\biggr\{\underset{k=1,2}{\operatorname{min}}\, |t_k|=|t_2| \biggr\} \biggr)^{\mytrans}.
\end{align*}

With this construction, the value of $\sqrt{n}\hatmaxpn$ equals the one in $\{T_{\alpha, n}, T_{\beta, n}\}$ that has a smaller absolute value, and $|\sqrt{n}\hatmaxpn| = \min\{|T_{\alpha, n}|, |T_{\beta, n}| \}$.
When  $T_{\alpha, n}$ and $T_{\beta, n}$ are asymptotically normal,

$\sqrt{n}\hatmaxpn$ equals the t-statistic

whose asymptotic $p$-value equals $p_{\mathrm{MaxP}}$.

This equivalence motivates us to further derive
a valid and non-conservative $p$-value
for $\sqrt{n}\hatmaxpn$ in the JS test.
We correspondingly define the centering parameter for $\hatmaxpn$ as $\maxpnull=H(\alphaS/\sigmaalpha, \betaM/ \sigmabeta )$.
Particularly, $\maxpnull$ follows the same form as $\hatmaxpn$ but replacing $(\hatalphaSn, \hatbetaMn, \hatsigmaalphan, \hatsigmabetan )$ in $\hatmaxpn$ with their population versions $(\alphaS, \betaM, \sigmaalpha, \sigmabeta )$,

and $\maxpnull=0$ under the composite null  (<ref>).

Simulation studies in Section <ref>  show that directly applying the classical nonparametric bootstrap to $\sqrt{n}(\hatmaxpn-\theta_0)$  fails to provide proper type 1 error control.
We next analytically unveil the non-standard limiting behavior of $\sqrt{n}(\hatmaxpn-\theta_0)$  before introducing our adaptive bootstrap test.

Non-Regularity of the JS Test.  The non-regular limiting behavior of $\hatmaxpn$ is caused by
the nonuniform convergence of the indicator vector $h(T_{\alpha, n}, T_{\beta, n})$ under different types of nulls.

Under $H_{0,1}$ or $H_{0,2}$,
$h(T_{\alpha, n}, T_{\beta, n})$ converges to $h(\alphaS/\sigmaalpha, \betaM/\sigmabeta )$ consistently, and $\hatmaxpn$ is asymptotically normal.
However, under $H_{0,3}$, $h(T_{\alpha, n}, T_{\beta, n})$ does not converge to $h(\alphaS/\sigmaalpha, \betaM/\sigmabeta )$,
and  $\hatmaxpn$ does not have a normal limit.
More specifically, we characterize the limiting distribution of $\sqrt{n}(\hatmaxpn-\maxpnull)$ considering a  special case of (<ref>):
$\M = \alphaS \S +\eM,$ and $\Y = \betaM \M +  \eY$,  and assuming $\sigmaalpha=\sigmabeta=1$.
Under mild conditions,
by the strong law of large numbers,
$(\hatalphaSn, \hatbetaMn, \hatsigmaalpha,  \hatsigmabeta) \xrightarrow{a.s.} (\alphaS, \betaM, 1, 1)$.

Then  we can write

\begin{align}
\sqrt{n}(\hatmaxpn-\maxpnull)=&~\sqrt{n}(\alphaS, \betaM)\times \big\{ h(\hatalphaSn , \hatbetaMn ) -h(\alphaS, \betaM ) \big\}	 \label{eq:jstsimple} \\
&~+\sqrt{n}\big( \hatalphaSn - \alphaS, \hatbetaMn - \betaM \big)\times h(\hatalphaSn , \hatbetaMn) + o_{P_n}(1),\notag
\end{align}

where $o_{P_n}(1)$ represents a small order term converging to 0 in probability.

Under $H_{0,1}$ or $H_{0,2}$, (or more generally, when $|\alphaS|\neq |\betaM|$),

$h(t_1,t_2)$ is continuous at $(t_1,t_2)=(\alphaS, \betaM)$ by the construction of $h(t_1,t_2)$;

the continuous mapping theorem then implies that $h(\hatalphaSn , \hatbetaMn) \xrightarrow{a.s.} h(\alphaS, \betaM ).$
Under $H_{0,3}$, we have  $h(\hatalphaSn , \hatbetaMn)=h(\sqrt{n}\hatalphaSn, \sqrt{n}\hatbetaMn )$,
with $\sqrt{n}(\hatalphaSn, \hatbetaMn ) \xrightarrow{d} (Z_{\S,0}, \, Z_{\M,0})$ by (<ref>).

With these results, by (<ref>) and Slutsky's lemma,
we have $\sqrt{n}(\hatmaxpn - \maxpnull ) \xrightarrow{d} U_{2,0}$, where

\begin{align}\label{eq:u20def}
U_{2,0}= \begin{cases}
\displaystyle (Z_{\S,0},Z_{\M,0}) \times h(\alphaS, \betaM), &  \ \text{if }\, |\alphaS|\neq |\betaM|; \\[10pt]
\displaystyle (Z_{\S,0},Z_{\M,0}) \times h (Z_{\S,0},Z_{\M,0}), &  \ \text{if }\, (\alphaS, \betaM) = (0,0).
\end{cases}
\end{align}

We can see that  the distribution of $\sqrt{n}(\hatmaxpn - \maxpnull)$ does not converge uniformly with respect to $(\alphaS, \betaM),$ and the nonuniformity occurs at the neighborhood of $(\alphaS, \betaM)=(0, 0).$
This discontinuity phenomenon leads to the failure of classical nonparametric bootstrap, which was also  demonstrated by the simulation studies  in Section <ref>.

Adaptive Bootstrap of the JS Test.  Similar to Section <ref>,
to develop a consistent bootstrap procedure, we need to understand the limiting behavior of $\sqrt{n}(\hatmaxpn - \maxpnull )$
in a local neighborhood of $(\alphaS, \betaM)=(0,0)$.
To achieve this,

again we consider  the local linear SEM  (<ref>),
where we recall that  $\alphaSn = \alphaS + n^{-1/2}\balpha$ and $\betaMn = \betaM + n^{-1/2}\bbeta.$
Then we correspondingly define the centering parameter under the local linear SEM as
$\maxpnulln = H( \alphaSn/\sigmaalpha, \betaMn/\sigmabeta )$.
Note that
$\maxpnulln$ takes the same form as $\maxpnull$, except that $\alphaS$ and $\betaM$ are replaced by $\alphaSn$ and $\betaMn$, respectively.
We present the limiting distribution of $\sqrt{n}(\hat{\theta}_n -\maxpnulln)$ under Model (<ref>) in  the following theorem.

Assume Condition <ref> holds  and   $|\alphaS/\sigmaalpha| \neq |\betaM/ \sigmabeta|$ when $(\alphaS, \betaM) \neq (0,0)$.
Then, under the local linear SEM (<ref>), $\sqrt{n}(\hat{\theta}_n -\maxpnulln)	\xrightarrow{d} U_{2}$, as $n \rightarrow \infty$ with

\begin{eqnarray*}
U_{2} =   \begin{cases}
\displaystyle \biggr(\frac{Z_{\S}}{\sigma_{\alphaS}}, \frac{Z_{\M}}{ \sigma_{\betaM}} \biggr)\times h\biggr(\frac{\alphaS}{\sigmaalpha }, \frac{\betaM}{\sigmabeta} \biggr), &\ \text{if }\, (\alphaS, \betaM) \neq (0,0); \\[12pt]
\displaystyle H (K_{b,\S},K_{b,\M}) -H\biggr(\frac{\balpha}{ \sigmaalpha }, \frac{\bbeta}{\sigmabeta} \biggr),
& \ \text{if }\, (\alphaS, \betaM) = (0,0), \\[8pt]
%\displaystyle  (K_{b,\S},K_{b,\M})\times h (K_{b,\S},K_{b,\M}) -\biggr(\frac{\balpha}{ \sigmaalpha }, \frac{\bbeta}{\sigmabeta} \biggr) h \biggr(\frac{\balpha}{ \sigmaalpha }, \frac{\bbeta}{\sigmabeta} \biggr)
% & \ \text{if }\, (\alphaS, \betaM) = (0,0). \\[8pt]
%~ \displaystyle -\biggr(\frac{\balpha}{ \sigmaalpha }, \frac{\bbeta}{\sigmabeta} \biggr) h \biggr(\frac{\balpha}{ \sigmaalpha }, \frac{\bbeta}{\sigmabeta} \biggr)
\end{cases}
\end{eqnarray*}

where $(Z_{\S},  Z_{\M})$ are defined the same as in Theorem <ref>, and

\begin{align*}
K_{b,\S}= \frac{\balpha+Z_{\S}}{\sigma_{\alphaS}}, \hspace{3em} K_{b,\M}= \frac{\bbeta+Z_{\M}}{\sigma_{\betaM}}.
\end{align*}

The assumption $|\alphaS/\sigmaalpha| \neq |\betaM/ \sigmabeta|$ when $(\alphaS, \betaM) \neq (0,0)$ is satisfied under the composite null (<ref>),
and is made mainly for the technical simplicity in the proof.
When $(\alphaS, \betaM)=(0,0)$,

$H (K_{b,\S},K_{b,\M}) - H({\balpha}/{ \sigmaalpha }, {\bbeta}/{\sigmabeta} )$ in Theorem <ref> can be equivalently written as

\begin{align*}
\biggr(\frac{Z_{\S}}{\sigma_{\alphaS}}, \frac{Z_{\M}}{\sigma_{\betaM}} \biggr) h(K_{b,\S}, K_{b,\M} )+ \biggr(\frac{\balpha}{ \sigmaalpha }, \frac{\bbeta}{\sigmabeta} \biggr) \biggr\{h(K_{b,\S}, K_{b,\M} ) -h\biggr(\frac{\balpha}{ \sigmaalpha }, \frac{\bbeta}{\sigmabeta} \biggr)\biggr\}.
\end{align*}

Comparing the above expression to the form of $U_2$ when $(\alphaS, \betaM )\neq (0,0)$,
we can see $U_2$ is discontinuous with respect to the parameters $(\alphaS, \betaM)$.
On the other hand, at the $\sqrt{n}$-neighborhood of  $(\alphaS, \betaM)=(0,0)$, Theorem <ref> further shows that the limiting distribution of  $\sqrt{n}(\hatmaxpn - \maxpnulln)$ is  continuous as a function of $(\balpha, \bbeta )^{\mytrans}\in \realnumber^2$ into the space of distribution functions.
Therefore, the non-regularity at $(\alphaS, \betaM)=(0,0)$ can be  explained by the dependence
of the limiting distribution on the (nonidentifiable) local parameters $(\balpha, \bbeta)$.

Similarly to Section <ref>,
we expect that developing  a bootstrap test using the local asymptotic results in Theorem <ref> will improve the finite-sample accuracy,

whereas the classical bootstrap methods that do not take into account the local asymptotic behaviors will fail to provide consistent estimates of the distribution of $\sqrt{n}(\hatmaxpn - \maxpnulln)$.

Similarly to the PoC test in Section <ref>,
to overcome the local discontinuity issue,  we isolate the possibility that $(\alphaS, \betaM )\neq (0,0)$ by comparing the standardized statistics $|T_{\alpha,n}|$ and $|T_{\beta,n}|$ to certain threshold $\lambda_n$. Specifically, we decompose

\begin{align}
\sqrt{n}(\hatmaxpn - \maxpnulln )= \sqrt{n}(\hatmaxpn - \maxpnulln)\times(1-\myindicator_{\alphaS, \lambda_n}\myindicator_{\betaM, \lambda_n} )+\sqrt{n}(\hatmaxpn - \maxpnulln)\times \myindicator_{\alphaS, \lambda_n}\myindicator_{\betaM, \lambda_n}, \label{eq:decompform}
\end{align}

where $\myindicator_{\alphaS, \lambda_n}$ and $\myindicator_{\betaM, \lambda_n}$ are defined same as  those in (<ref>).

When $(\alphaS, \betaM) \neq (0,0)$
and $|\alphaS/\sigmaalpha| \neq |\betaM/ \sigmabeta|$,
the classical nonparametric bootstrap estimator $\sqrt{n}(\hatmaxpn^*- \hatmaxpn)$, where  $\sqrt{n}\hatmaxpn^*=H(T_{\alpha, n}^*, T_{\beta, n}^*)$, is consistent for the first term $\sqrt{n}(\hatmaxpn-\maxpnulln)$ in (<ref>).
For the second term in (<ref>),
we have $(\alphaS, \betaM) = (0,0)$  and can write
$\sqrt{n}(\hatmaxpn - \maxpnulln)=\mathbb{R}_{2,n}(\balpha, \bbeta )$,

where $ \mathbb{R}_{2,n}(\balpha, \bbeta )= H(\mathbb{K}_{b,\S}, \mathbb{K}_{b,\M})-H({\balpha}/{ \sigmaalpha }, {\bbeta}/{\sigmabeta}),$

\begin{align*}
\mathbb{K}_{b,\S}=\frac{ \balpha+\mathbb{Z}_{\S,n}}{\hatsigmaalphan }, \hspace{2em} \mathbb{K}_{b,\M}=\frac{\bbeta+\mathbb{Z}_{\M,n}}{\hatsigmabetan },
\end{align*}

and $(\mathbb{Z}_{\S,n}, \mathbb{Z}_{\M,n})$ are defined same as those in Section <ref>.
It follows that all parts of $\mathbb{R}_{2,n}(\balpha, \bbeta )$ can be viewed as smooth functions of $\mathbb{P}_n$.
Similarly to Section <ref>, it is reasonable to expect that a consistent bootstrap can
be constructed using  the nonparametric bootstrap version of $\mathbb{R}_{2,n}(\balpha, \bbeta )$.
Specifically, we define
$\mathbb{R}_{2,n}^*(\balpha, \bbeta)=H(\mathbb{K}_{b,\S}^*, \mathbb{K}_{b,\M}^*)-H({\balpha}/{\hatsigmaalphan}, {\bbeta}/{\hatsigmabetan} ),$

\begin{align*}
\mathbb{K}_{b,\S}^*=\frac{ \balpha+\mathbb{Z}_{\S,n}^*}{ \hatsigmaalphan^*},  	\hspace{2em}\mathbb{K}_{b,\M}^*=\frac{ \bbeta+\mathbb{Z}_{\M,n}^*}{ \hatsigmabetan^*}. \notag
\end{align*}

We are ready to develop our adaptive bootstrap test based on (<ref>).
Similarly to Section <ref>,
we replace the indicators $\myindicator_{\alphaS, \lambda_{n}}$ and $\myindicator_{\betaM, \lambda_{n}}$  in (<ref>) with $\myindicator_{\alphaS, \lambda_{n}}^*$ and $\myindicator_{\betaM, \lambda_{n}}^*$ in (<ref>), respectively.
Then we
define our proposed  adaptive bootstrap test statistic as

\begin{align*}
U_2^*= \sqrt{n}(\hatmaxpn^* - \hatmaxpn) \times (1- \myindicator_{\alphaS, \lambda_n}^* \myindicator_{\betaM, \lambda_n}^* ) + \mathbb{R}_{2,n}^*(\balpha, \bbeta) \times \myindicator_{\alphaS, \lambda_n}^* \myindicator_{\betaM, \lambda_n}^*.
\end{align*}

Theorem <ref> below proves that $\sqrt{n}(\hatmaxpn - \maxpnulln )$ can be consistently bootstrapped using the $U_2^*$ above.

Assume that the conditions in Theorem <ref> are satisfied,
and that the tuning parameter $\lambda_n$ satisfies $\lambda_n=o(\sqrt{n})$ and $\lambda_n\to \infty$ as $n\to \infty.$

Then under the local linear SEM (<ref>),
conditionally on the data,
the adaptive test statistic $U_2^*\overset{d^*}{\leadsto} \sqrt{n}(\hatmaxpn - \maxpnulln )$.

Theorem <ref>  establishes the bootstrap consistency of $U_2^*$ for $\sqrt{n}(\hatmaxpn - \maxpnulln )$,
which is different from Theorem <ref> that derives  $nU_1^*$ and $\sqrt{n}U_1^*$ separately for $n(\hatalphaS \hatbetaM - \alphaSn \betaMn) $ and $\sqrt{n}(\hatalphaS \hatbetaM - \alphaSn \betaMn) $ when $(\alphaS,\betaM)=(0,0)$ and $(\alphaS,\betaM)\neq (0,0)$.
This difference is attributed to the distinct   non-regular limiting behaviors in  the PoC test and the JS test.
Despite this difference, based on Theorem <ref>,
we can develop an adaptive bootstrap test procedure
for the JS test
similarly to  that in Section <ref>.
In particular,
under the composite null (<ref>),
we consider  $(\balpha, \bbeta)= (0,0)$ in $U_2^*$ to estimate  the distribution of $\sqrt{n}(\hatmaxpn^*-\hatmaxpn)$.
For a given nominal level $\omega$, we redefine $q_{\omega/2}$ and $q_{1-\omega/2}$ as the lower and upper $\omega/2$ quantiles, respectively, of $U_{2}^*$.
If $\sqrt{n}\hatmaxpn$ falls outside  the interval $(q_{\omega/2}, q_{1-\omega/2} )$, then we reject the composite null, and conclude that the mediation effect is statistically significant.

In addition,
we also choose the   value of the tuning parameter $\lambda_n$ such that $\lambda_n=o(\sqrt{n})$ and $\lambda_n\to \infty$ as $n\to \infty$ following  the discussions in Section <ref>.

Moreover, similarly to Remark <ref>,
we  emphasize that the
analysis of the JS test statistic  is distinct from the  existing literature as we consider different problem settings and testing a composite null hypothesis, which necessitates new theoretical and methodological developments.

§ PROOFS OF THEOREMS <REF>–<REF> AND THEOREMS

In this section, we develop preliminary results for the follow-up proofs in Section <ref>.
We provide proofs of Theorems <ref>, <ref>, <ref>, and <ref>
in Sections <ref>–<ref>, respectively.

§.§ Preliminary for the Proofs

In the following proofs, we use the variables defined as:

\begin{align}
&\S_{\perp} = \S - \X^{\mytrans}Q_{1,\S}, \quad \quad \quad  \M_{\perp} = \M - \X^{\mytrans} Q_{1,\M},  \label{eq:newdefvariables} \\
&\M_{\perp'} = \M - \tilde{\X}^{\mytrans} Q_{2,\M}, \quad \quad \Y_{\perp'} = \Y - \tilde{\X}^{\mytrans}Q_{2,\Y}, \notag
\end{align}

where $\tilde{\X}=(\X^{\mytrans},\S)^{\mytrans}$,

\begin{align}
&Q_{1,\M} = \{P_n(\X \X^{\mytrans})\}^{-1} P_n(\X \M), \quad \quad Q_{1,\S}= \{P_n(\X \X^{\mytrans})\}^{-1} P_n(\X \S), \label{eq:popqmomentdef} \\
&  Q_{2,\M} = \{P_n( \tilde{\X}  \tilde{\X}^{\mytrans})\}^{-1} P_n( \tilde{\X} \M), \quad \quad Q_{2,\Y}= \{P_n( \tilde{\X}  \tilde{\X}^{\mytrans})\}^{-1} P_n( \tilde{\X} \Y).  \notag
\end{align}

We mention that $\Sperp$ and $\Mperpp$ are defined same as those in Condition <ref>,
where both $P_n(\cdot)$ used in (<ref>) and  $\mathrm{E}(\cdot)$ used in Condition <ref> denote the expectation with respect to the distribution of $(\S,\M,\X,\Y)$.

Based on (<ref>), for each index $i \in \{1,\ldots, n\}$, we define $\Sperpi =\S_i-\X_i^{\mytrans}Q_{1,\S}$,
and also define $(\Mperpi, \Mperppi, \Yperppi)$ similarly.

Moreover, for the definitions in (<ref>) and (<ref>),
by replacing $P_n(\cdot)$ with $\mathbb{P}_n(\cdot)$,