2.1 M-regression

To alleviate the effects of influential observations in the least-squares fitting, we adopt the M-regression to estimate the regression parameters, which can be regarded as a generalization of the maximum likelihood estimation as follows:

(4) $$ \begin{align} (\hat{\beta}_{1}, \hat{c})^{T}&=\arg\min_{\beta_{1}, c}\sum_{i=1}^{n}\rho(Y_{i}-\beta_{1}-cX_{i}), \end{align} $$

(5) $$ \begin{align} (\hat{\beta}_{2}, \hat{a})^{T}&=\arg\min_{\beta_{2}, a}\sum_{i=1}^{n}\rho(M_{i}-\beta_{2}-aX_{i}), \end{align} $$

(6) $$ \begin{align} (\hat{\beta}_{3}, \hat{c^{\prime}}, \hat{b})^{T}&=\arg\min_{\beta_{3}, c^{\prime},b} \sum_{i=1}^{n}\rho(Y_{i}-\beta_{3}-c^{\prime}X_{i}-bM_{i}), \end{align} $$

where $\rho (\cdot )$ is the loss function with three properties: (i) non-negativity such that $\rho (\epsilon )\geq 0$ with $\rho (0)=0$ , (ii) symmetricity such that $\rho (\epsilon )=\rho (-\epsilon )$ , and (iii) monotonicity such that $\rho (\epsilon )\geqslant \rho (\epsilon ^{\prime })$ for any $|\epsilon |\geqslant |\epsilon ^{\prime }|$ .

Let $\psi (\epsilon ) = (d/d\epsilon )\rho (\epsilon )$ be the first derivative of the loss function, referred to as the influence curve. Let also $X=(X_{1}, \ldots , X_{n})^{T}$ , $M=(M_{1}, \ldots , M_{n})^{T}$ , $Y=(Y_{1}, \ldots , Y_{n})^{T}$ , $I=(1, \ldots , 1)^{T}$ , $\tilde {\boldsymbol {X}}=(I, X)$ , and $\check {\boldsymbol {X}}=(I, X, M)$ . For large samples, we further assume that $\boldsymbol {U}$ is the limiting matrix of $(n^{-1}\tilde {\boldsymbol {X}}^{T}\tilde {\boldsymbol {X}})^{-1}$ , and $\boldsymbol {V}$ is the limiting matrix of $(n^{-1}\check {\boldsymbol {X}}^{T}\check {\boldsymbol {X}})^{-1}$ . Then by Huber and Ronchetti (Reference Huber and Ronchetti2009), we have the following asymptotic normality for the M-estimators of the regression parameters.

Lemma 1. For the mediation model linked with (1)–(3), under the regularity conditions given on pages 163–164 of Huber and Ronchetti (Reference Huber and Ronchetti2009), the M-estimators in (4)–(6) are all normally distributed:

$$ \begin{align*}\sqrt{n}(\hat{c}-c)&\sim N\left(0, \frac{\mathrm{E}_{\epsilon_{1}}[\psi^{2}]}{(\mathrm{E}_{\epsilon_{1}}[\psi^{\prime}])^{2}} \boldsymbol{U}_{[2, 2]}\right), & \sqrt{n}(\hat{a}-a)&\sim N\left(0, \frac{\mathrm{E}_{\epsilon_{2}}[\psi^{2}]}{(\mathrm{E}_{\epsilon_{2}}[\psi^{\prime}])^{2}} \boldsymbol{U}_{[2, 2]}\right), \\ \sqrt{n}(\hat{c^{\prime}}-c^{\prime})&\sim N\left(0,\frac{\mathrm{E}_{\epsilon_{3}}[\psi^{2}]} {(\mathrm{E}_{\epsilon_{3}}[\psi^{\prime}])^{2}}\boldsymbol{V}_{[2, 2]}\right), & \sqrt{n}(\hat{b}-b)&\sim N\left(0, \frac{\mathrm{E}_{\epsilon_{3}}[\psi^{2}]} {(\mathrm{E}_{\epsilon_{3}}[\psi^{\prime}])^{2}}\boldsymbol{V}_{[3, 3]}\right). \end{align*} $$

Finally, based on the M-estimators in (4)–(6), we can define two new estimators of the indirect effect: one is the difference estimator $\hat {c}-\hat {c^{\prime }}$ and the other is the product estimator $\hat {a}\hat {b}$ .

2.2 Solution to M-regression

For a general loss $\rho (\cdot )$ , noting that the M-estimator may not have an explicit expression, a numerical solution is often required. To present our algorithm, we will focus only on (4) since the same algorithm can be extended to solve (5) and (6) as well. Differentiating the objective function $\sum _{i=1}^{n}\rho (Y_{i}-\beta _{1}-cX_{i})$ with respect to $\beta _{1}, c$ and setting the partial derivatives to be zero, it yields a system of two estimating equations as

$$ \begin{align*} \sum_{i=1}^{n}\psi(Y_{i}-\beta_{1}-cX_{i})&=0, \\ \sum_{i=1}^{n}\psi(Y_{i}-\beta_{1}-cX_{i})X_{i}&=0. \end{align*} $$

Further by introducing the weight function $w(e)=\psi (e)/e$ , the estimating equations can be rewritten as

$$ \begin{align*} \sum_{i=1}^{n}w_{i}\times(Y_{i}-\beta_{1}-cX_{i})&=0, \\ \sum_{i=1}^{n}w_{i}\times(Y_{i}-\beta_{1}-cX_{i})X_{i}&=0, \end{align*} $$

where $w_{i}=w(Y_{i}-\beta _{1}-cX_{i})$ . Solving these two equations is equivalent to minimizing

$$ \begin{align*} \sum_{i=1}^{n}w_{i}\times(Y_{i}-\beta_{1}-cX_{i})^2, \end{align*} $$

which is a weighted LS problem. Moreover, an iteratively reweighted least-squares (IRLS) algorithm can be appropriate to obtain the numerical solution of the regression coefficients, because the weights depend on the regression coefficients, and the regression coefficients in turn depend on the weights (Holland & Welsch, Reference Holland and Welsch1977). To also handle the multiple-minima problem, in case it has, we choose several different points in the parameter space as the initial estimates, in such a way to get a higher confidence to obtain the true global minimum (Green, Reference Green1984). More specifically, the IRLS algorithm for our problem is as follows.

2.3 Error conditions for model consistency

When is the product of parameters $ab$ equal to the difference in parameters $c-c^{\prime }$ in population? This is an important question in mediation analysis since it uncovers the relationship between the indirect, direct and total effects (Wang et al., Reference Wang, Yu, Zhou, Tong and Liu2023; Wang & Yu, Reference Wang and Yu2023; Yuan & MacKinnon, Reference Yuan and MacKinnon2014).

Note that the three regression equations, (1)-(3), are interrelated in the simple mediation model. By substituting (2) into (3), we have

(7) $$ \begin{align} Y_{i}&=\beta_{3}+c^{\prime}X_{i}+b(\beta_{2}+aX_{i}+\epsilon_{2, i})+\epsilon_{3, i} \notag \\ &=(\beta_{3}+b\beta_{2})+(c^{\prime}+ab)X_{i}+\epsilon_{i}, \end{align} $$

where $\epsilon _{i}=b\epsilon _{2, i}+\epsilon _{3, i}$ . Assume that $\epsilon _{2, i}$ and $\epsilon _{3, i}$ are independent and symmetrically distributed with median $0$ , then $\epsilon _{i}$ is also symmetric with $\mathrm {Med}[\epsilon _i]=0$ (see Proposition 1 in Wang and Yu (Reference Wang and Yu2023)). In addition, let $\epsilon _{1,i}$ also be symmetrically distributed with $\mathrm {Med}[\epsilon _{1,i}]=0$ . Then by (1) and (7),

$$ \begin{align*} \mathrm{Med}[Y_{i}|X_{i}]&=\beta_{1}+cX_{i}+\mathrm{Med}[\epsilon_{1, i}|X_{i}], \\ \mathrm{Med}[Y_{i}|X_{i}]&=(\beta_{3}+b\beta_{2})+(c^{\prime}+ab)X_{i}+\mathrm{Med}[\epsilon_{i}|X_{i}]. \end{align*} $$

Noting also that the random errors are independent of the corresponding regressors as assumed in Section 2.1, we have $\mathrm {Med}[\epsilon _{1, i}|X_i]=\mathrm {Med}[\epsilon _{1, i}]=0$ and $\mathrm {Med}[\epsilon _{i}|X_{i}]=\mathrm {Med}[\epsilon _{i}]=0$ , and moreover,

$$ \begin{align*} \beta_{1}+cX_{i}\equiv (\beta_{3}+b\beta_{2})+(c^{\prime}+ab)X_{i}, \quad i=1,\ldots,n, \end{align*} $$

which further yields that $\beta _{1}=\beta _{3}+b\beta _{2}$ and $c=c^{\prime }+ab$ . Finally, by comparing (1) and (7), we also have $\epsilon _{i}=\epsilon _{1, i}$ . For convenience, we summarize the above result in Theorem 1.

2.4 Inference based on confidence interval

There are two estimators for the indirect effect: $\hat {c}-\hat {c^{\prime }}$ and $\hat {a}\hat {b}$ . Unlike the equivalence of the two LS estimators (MacKinnon et al., Reference MacKinnon, Warsi and Dwyer1995; Wang et al., Reference Wang, Yu, Zhou, Tong and Liu2023), the two M-estimators of the indirect effect for a general loss are not the same in general, that is, $\hat {a}\hat {b}\neq \hat {c}-\hat {c^{\prime }}$ . Simulation studies show that the product estimator is often more efficient than the difference estimator (see Appendix A0). Interestingly, the same conclusion can also be seen when the LAD loss is applied (Wang & Yu, Reference Wang and Yu2023). In view of this, we thus consider the null hypothesis $H_{0}: ab=0$ . To test whether $ab=0$ , there are two common methods in the literature including the parameter method (Sobel, Reference Sobel1982) and the nonparametric resampling method (MacKinnon et al., Reference MacKinnon, Lockwood and Williams2004; Preacher & Selig, Reference Preacher and Selig2012).

To move forward, our first method is to perform a robust Sobel test. Given the robust estimates $\hat {a}$ and $\hat {b}$ , we define the robust test statistic as

$$ \begin{align*} Z=\frac{\hat{a}\hat{b}}{\widehat{\mathrm{SE}}_{Sobel}}, \end{align*} $$

where $\widehat {\mathrm {SE}}_{Sobel}=\sqrt {\hat {a}^{2}\times \widehat {\mathrm {SE}}_{b}^{2}+\hat {b}^{2}\times \widehat {\mathrm {SE}}_{a}^{2}}$ , and $\widehat {\mathrm {SE}}_{a}$ and $\widehat {\mathrm {SE}}_{b}$ are the standard errors (SEs) of $\hat {a}$ and $\hat {b}$ , respectively. Following Theorem 1, the two SEs can be estimated by

$$ \begin{align*} \widehat{\mathrm{SE}}_{a}&=\left(\frac{n^{-1}\sum_{i=1}^{n}\psi^{2}(M_{i}-\hat{\beta}_{2}-\hat{a}X_{i}) [(\tilde{\boldsymbol{X}}^{T}\tilde{\boldsymbol{X}})^{-1}]_{[2, 2]}} {[n^{-1}\sum_{i=1}^{n}\psi^{\prime}(M_{i}-\hat{\beta}_{2}-\hat{a}X_{i})]^{2}}\right)^{1/2}, \\ \widehat{\mathrm{SE}}_{b}&=\left(\frac{n^{-1}\sum_{i=1}^{n}\psi^{2}(Y_{i}-\hat{\beta}_{3}-\hat{c^{\prime}}X_{i}-\hat{b}M_{i}) [(\check{\boldsymbol{X}}^{T}\check{\boldsymbol{X}})^{-1}]_{[3, 3]}} {[n^{-1}\sum_{i=1}^{n}\psi^{\prime}(Y_{i}-\hat{\beta}_{3}-\hat{c^{\prime}}X_{i}-\hat{b}M_{i})]^{2}}\right)^{1/2}. \end{align*} $$

Moreover, the normal-based $(1-\alpha )\%$ CI of $ab$ can be constructed as

$$ \begin{align*} [\hat{a}\hat{b}-z_{1-\alpha/2}\widehat{\mathrm{SE}}_{Sobel}, \: \hat{a}\hat{b}+z_{1-\alpha/2}\widehat{\mathrm{SE}}_{Sobel}], \end{align*} $$

where $\alpha $ is the significance level, and $z_{1-\alpha /2}$ represents the $(1-\alpha /2)$ quantile of the standard normal distribution. Note however that, when a and b are small, the sampling distribution of $\hat {a}\hat {b}$ may not be normal (MacKinnon et al., Reference MacKinnon, Lockwood and Williams2004; Wang et al., Reference Wang, Yu, Zhou, Tong and Liu2023). Thus to obtain an accurate CI, critical values of the distribution of $\hat {a}\hat {b}$ can be obtained by Monte Carlo simulation study (Meeker et al., Reference Meeker, Cornwell and Aroian1981; Meeker & Escobar, Reference Meeker and Escobar1994). In fact, one can easily obtain these critical values via inputting $\hat {a}$ , $\hat {b}$ , $\widehat {\mathrm {SE}}_{a}$ and $\widehat {\mathrm {SE}}_{b}$ into an R procedure medci() which was introduced by Tofighi and MacKinnon (Reference Tofighi and MacKinnon2011).

Our second method to construct CI is the bootstrap method based on resampling. The bootstrap method is nonparametric and robust in the sense that it does not need to estimate the SEs. First, we repeatedly resample the original dataset with replacement (Efron & Tibshirani, Reference Efron and Tibshirani1993); second, we estimate the indirect effect for each bootstrap sample using our proposed Huber method; third, we construct the CI by the percentile bootstrap (PRCT) as $[q_{\alpha /2}, q_{1-\alpha /2}]$ , where $q_{\alpha /2}$ is the $\alpha /2$ quantile of the empirical distribution of the indirect effect. To adjust and remove the potential estimation bias, the bias-corrected and accelerated bootstrap (BCa) is an important variation (Efron, Reference Efron1987; Efron & Tibshirani, Reference Efron and Tibshirani1993). In general, the BCa method can yield a more accurate CI than the PRCT method when the true parameter value is not the median of the distribution of the bootstrap estimates (MacKinnon et al., Reference MacKinnon, Lockwood and Williams2004).

3.1 Huber Loss

The Huber loss, as defined in Huber (Reference Huber1964), is given as

$$ \begin{align*} \rho_{H}(e)&=\left\{\begin{array}{ll} \frac{1}{2}e^{2}, & \text{if}\quad |e|\leq k, \\ k|e|-\frac{1}{2}k^{2}, &\text{if}\quad |e|>k, \end{array}\right. \\ \psi_{H}(e)&=\left\{\begin{array}{ll} e, & \text{if}\quad |e|\leq k, \\ k\times \text{sgn}(e), & \text{if}\quad |e|>k, \end{array}\right. \end{align*} $$

where $k>0$ is the tuning parameter. A smaller value of k produces more resistance to outliers, but at the expense of lower efficiency when the error is normal. For instance, by letting $k=1.345\sigma $ with $\sigma $ being the standard deviation of the error, it will yield a $95\%$ efficiency for the normal errors, which is also resistant to outliers with a breakdown point of 5.8%. Moreover, the standard deviation $\sigma $ can be estimated robustly by the median absolute deviation (MAD) as

$$ \begin{align*} \hat{\sigma}_{M\!A\!D}=\mathrm{Med}\{|e_{i}|\}/0.6745. \end{align*} $$

For any error $\epsilon $ , we denote $\tau =\sigma _{\psi }^{2}/B_{\psi }^{2}$ as the asymptotic variance of the Huber estimator (Huber, Reference Huber1964), where $\sigma _{\psi }^{2}=\mathrm {E}[\psi ^{2}(\epsilon )]$ and $B_{\psi }=\mathrm {E}[\psi ^{\prime }(\epsilon )]$ . We then minimize the $\tau $ value to determine the optimal $\rho (\cdot )$ . For the Huber loss with a given k, we have

$$ \begin{align*} B_{\psi}(k)&=\int_{-k}^{k}\mathrm{d}F(\epsilon), \\ \sigma_{\psi}^{2}(k)&=\int_{-k}^{k}\epsilon^{2}\mathrm{d}F(\epsilon)+k^{2}(1-B_{\psi}(k)), \end{align*} $$

where $F(\cdot )$ is the cumulative distribution function of $\epsilon $ .

3.2 Optimal Tuning constant

As is known, the tuning parameter k of the Huber loss can have a great impact on the estimation efficiency. When the error is normally distributed without contamination, the best choice of k is $\infty $ . On the other hand, when the error follows a heavy-tailed distribution such as the t distribution, then k tends to be a small value close to 0.

We adopt a numerical method proposed by Wang et al. (Reference Wang, Lin, Zhu and Bai2007) to select the optimal tuning constant, which minimizes the asymptotic variance of the estimator. For the Huber loss, the optimal k minimizes the efficiency factor $\tau $ with a three-step procedure as follows. First, we compute $\tau (k)$ for a range of k values, i.e., $0\le k\leq K$ by $0.001$ , where K is a positive number, e.g., $K=4$ . Second, we select the optimal k as

$$ \begin{align*} k_{opt}=\arg\min_{0<k\leq K}\tau(k). \end{align*} $$

Lastly, we compute the minimum value $\tau (k_{opt})$ . In Appendix B, we provide an R procedure to obtain the optimal tuning constant with a known error distribution.

For ease of reference, we also list the optimal $k_{opt}$ and $\tau (k_{opt})$ in Table 1 for some error distributions, including the standard normal distribution $N(0, 1)$ , the Laplace distribution Laplace $(0, 1)$ , the mixed normal distribution $0.9N(0, 1)+0.1N(0, \sigma ^{2})$ with $\sigma =3$ or $10$ , and the t distribution with $1$ or $2$ degrees of freedom. In general, the Huber loss with the optimal tuning parameter k is more efficient than the LS and LAD losses, since the less $\tau $ is, the more efficient the loss is. Moreover, to intuitively reflect the variation trend of $\tau (k)$ as k varies, we also plot the $\tau (k)$ function for a normal mixed and $t_{1}$ distributions in Figure 2. It is evident that the value of $\tau (k)$ varies dramatically along with the k value.

Figure 2 $\tau (k)$ is plotted for $0.9N(0, 1)+0.1N(0, 3^{2})$ (left) and $t_{1}$ (right). The corresponding red lines are $\tau (1.489)=1.296$ and $\tau (0.395)=2.278$ , respectively.

Table 1 Optimal k and $\tau (k)$ for various error distributions and loss functions

3.3 Nonparametric selection of Tuning constant

Following (4) and letting $e_{i}=Y_{i}-\hat {\beta }_{1}-\hat {c}X_{i}$ be the residuals, we propose to estimate $\tau $ nonparametrically by

(8) $$ \begin{align} \hat{\tau}(k)=\frac{\hat{\sigma}_{\psi}^{2}(k)}{\hat{B}_{\psi}^{2}(k)}, \end{align} $$

where $\hat {\sigma }_{\psi }^{2}(k)=n^{-1}\sum _{i=1}^{n}\psi ^{2}(e_{i})$ and $\hat {B}_{\psi }(k)=n^{-1}\sum _{i=1}^{n}\psi ^{\prime }(e_{i})$ . More specifically for the Huber loss, we have

$$ \begin{align*} \hat{B}_{\psi}(k)&=\frac{1}{n}\sum_{i=1}^{n}\mathrm{{I}}(|e_{i}|\leq k),\\ \hat{\sigma}_{\psi}^{2}(k)&=\frac{1}{n}\sum_{i=1}^{n}\left\{e_{i}^{2} \mathrm{{I}}(|e_{i}|\leq k)+k^{2}\mathrm{{I}}(|e_{i}|>k)\right\}, \end{align*} $$

where $\mathrm {{I}}$ is the $0-1$ indicator function.

We propose a data-driven procedure that determines the optimal $\hat {k}$ by minimizing $\hat {\tau }(k)$ , which is, in fact, similar to Wang et al. (Reference Wang, Lin, Zhu and Bai2007) for a linear regression model with a scale parameter $\sigma $ . Our new procedure is summarized in Algorithm 2.

Note that in the algorithm, we have specified the maximum allowable k as $3\hat {\sigma }_{M\!A\!D}$ , which is often treated as sufficient since the probability that the errors fall within the interval $[-3\hat {\sigma }_{M\!A\!D}, 3\hat {\sigma }_{M\!A\!D}]$ is as large as 99.73% for the normal errors. To further investigate the performance of the proposed method on the selection of tuning constant k, we conduct a simulation study and report the results in Table B of the Appendix. When the sample size is large, the selected tuning constant k is very close to the theoretical one. Moreover, we note that the standard deviation of the tuning constant decreases dramatically as the sample size increases. These findings coincide with the conclusion in Wang et al. (Reference Wang, Lin, Zhu and Bai2007).

4.1 Efficiency of the LS, LAD, and Huber estimators

The MSE( $\times 10^3$ ) and standard deviation (SD $\times 10^3$ ) of the LS, LAD, and Huber estimators are presented in Table 2 for various designs. Comparing the MSE of the three estimators, we have two main findings. First, the MSE and SD of the three estimators decrease as the sample size increases. Second, the MSE and SD of the Huber estimator are always the smallest or close to the smallest. When the error follows $N(0,1)$ (or Laplace $(0,1)$ ), the LS (or LAD) estimator provides the optimal estimation. In these two cases, the Huber estimator performs very close to the performance of the optimal estimator. While for $ 0.9N(0,1)+0.1N(0,10^2) $ and $ t_2 $ , the MSE of the Huber estimator is the smallest among the three estimators. To conclude, the Huber estimator is more efficient than the LAD estimator when the error distribution is normal, and is more robust than the LS estimator when the error distribution is non-normal.

Table 2 MSE ( $\times 10^3$ ) and SD ( $\times 10^3$ labeled below MSE) for the LS, LAD, and Huber estimators

Note: Note that the bold font indicates the samllest MSE among the three estimators under one set of experimental conditions.

4.2 Type I error rate and power

We now apply the Sobel Z, PRCT and BCa methods to construct the 95% CI. Note that the medium effect sizes ( $a=b=0.39$ ) will yield a high power even when the sample size is moderate ( $n=200$ ). Thus to save space, we omit the simulation for the large effect size.

Table 3 report the type I error rates of the three estimators under various designs. When the sample size is large, i.e., $n=1000$ , we note that the type I error rates of the LS, LAD, and Huber estimators are all controlled in most cases. One exception is the LS estimator with the CIs constructed by the BCa method, which was also observed by Fritz et al. (Reference Fritz, Taylor and MacKinnon2012) with an explanation that the increased type I error rate is a function of an interaction between the nonzero effect size and the sample size. Another notable situation is that the type I error rate of the Huber loss Sobel test is slightly too high for the mixed normal and $t_2$ under the small and moderate sample sizes. Possible reasons can be, e.g., the standard error used for the Sobel test $\widehat {\mathrm {SE}}_{Sobel}$ is affected by the Optimizer’s curse (Smith & Winkler, Reference Smith and Winkler2006), and/or there is a potential gap between the optimal tuning constant and the one determined by Algorithm 2 in the small sample size. In Appendix E, we have also conducted another simulation study to assess their effect on the standard error used for the Sobel test. The results indicate that the $\widehat {\mathrm {SE}}_{Sobel}$ is indeed influenced by the optimizer’s curse, whereas its effect will diminish as the sample size increases. At the same time, the Huber estimator with the fixed $k=1.345$ performs better than the Huber estimator with the selected tuning constant (Huber-SEL) in the case of small sample size. Observing this, when the Huber-SEL estimator fails to yield satisfactory results, we suggest to take a moderate tuning constant, i.e., $k=1.345$ , as an alternative.

Table 3 Type I error rates (%) of the LS, LAD and Huber estimators for various designs

Note: Note that the bold font indicates the excessive type I error rate which exceeds 6.8% since with 1000 independent simulation runs, the type I error rate of a test with level 0.05 is expected lie in the interval $[2.3\%,6.8\%]$ with probability 0.99, using the normal approximation.

Following the same designs, we report the power of the three estimators in Table 4. For the normal errors, it is evident that the LS estimator not only controls the type I error rate but also achieves the highest power among the three estimators. Nevertheless, for the non-normal errors, the LS estimator is notably lacking in statistical power especially for the mixed normal distribution (e.g., $a=b=0.14$ , $0.9N(0,1)+0.1N(0,10^2)$ , and $n=1000$ ). In addition, despite that the LAD estimator is the most robust method with respect to the outliers, it however suffers from the efficiency loss and consequently yields a lower power (e.g., $a=b=0.14$ , $N(0,1)$ , and $n=1000$ ). In contrast, the Huber estimator makes a trade-off between the efficiency and robustness, in which its power is close to the largest and, meanwhile, it also controls the type I error rate below 5% regardless of the error distribution.

Table 4 Power (%) of the LS, LAD and Huber estimators for various designs

Note: Note that the bold font indicates the maximal empirical power among the three estimators under one set of experimental conditions.

5.1 Pathways to desistance study

Our first study is to uncover the causal mechanisms between mental health and violent offending among serious adolescent offenders (Kim et al., Reference Kim, Harris and Lee2024). In criminology, one possible mechanism is that individuals with mental health issues may be more likely to experience victimization, and this, in turn, may lead to their committing a serious crime. Our data comes from the Pathways to Desistance (PTD) study, which consists of 1354 serious juvenile offenders in two sites, including the Maricopa County in Arizona (N = 654) and Philadephia County in Pennsylvania (N = 700), over the years from 2000 to 2010 (Mulvey et al., Reference Mulvey, Schubert and Piquero2013). Focusing on the data of baseline interviews, our study contains a total of 1195 respondents after the data cleansing.

Consider the linear mediation model,

$$ \begin{align*} \textit{Expvic}_i &= \beta_{2} + a \textit{Health}_i + \boldsymbol{\delta}_{1}^T\boldsymbol{Z}_i + \varepsilon_{2,i}\\ \textit{Offend}_i &= \beta_{3} + c^\prime \textit{Health}_i + b \textit{Expvic}_i + \boldsymbol{\delta}_{2}^T\boldsymbol{Z}_i+\varepsilon_{3,i} \end{align*} $$

where Health (mental health) is the independent variable, Expvic (experienced victimization) is the mediating variable, Offend (violent effending) is the response variable. In addition, $\boldsymbol {Z}$ denotes the matrix of other controlled variables including age, gender, enthnicity, family structure, parental warmth, alcohol, marijuana, gang membership, parental hostility, and unsupervised routine activities. We summarize the type and the measure of these variables in Appendix F.

To assess the normality assumption for the errors, we compute the skewness and kurtosis of the residuals of y after regressing on x and m and the residuals of m after regressing on x, and then report them in Table 5. These values, together with the Kolmogorov–Smirnov (KS) test, clearly suggest a violation of the normality assumption. In view of this, we thus apply our new method to this dataset and also compare it with the existing methods for mediation analysis. Table 6 reports the indirect effects and the 95% CIs constructed by the Sobel Z, PRCT and BCa methods. From the results, we note that the three estimators produce similar and statistically significant indirect effects, whereas the Huber estimator yields the shortest CI.

Table 5 Skewness and kurtosis of two regression residuals and the Kolmogorov-Smirnov test for the pathways to desistance study-

Table 6 The indirect effect estimates and their 95% CIs based on the LS, LAD and Huber estimators for the pathways to desistance study

5.2 Action planning study

Our second study is to investigate the relationship between action planning and physical activity. In psychology, it is known that the action planning can promote the physical activity, yet the underlying mechanism between them is often unclear. To explore it, an illustrative study has recently been conducted to investigate the action planning promoting the physical activity mediated by the automaticity (Maltagliati et al., Reference Maltagliati, Sarrazin, Isoard-Gautheur, Pelletier, Rocchi and Cheval2023), in which a total of 135 participants over 18 years from the tertiary industry were recruited. Participants were asked to wear an accelerometer Actigraph GT3X+, which records their physical activity behaviors and the time of these activities on a notebook for a total of seven days. More specifically in their study, the action planning is the independent variable, measured by four-item Likert scales ranging from 1 (completely disagree) to 6 (full agree). And the automaticity is the mediating variable, measured by four-item of Self-Reported Habit Index ranging from 1 (strongly disagree) to 7 (strongly agree).

Consider the linear mediation model,

(9) $$ \begin{align} \textit{Auto}_i &= \beta_{2} + a \textit{Plan}_i + \delta_1 \textit{Sex}_i + \delta_2 \textit{BMI}_i + \delta_3 \textit{Ill}_i + \varepsilon_{2,i}, \end{align} $$

(10) $$ \begin{align} \textit{PA}_i &= \beta_{3} + c^{\prime}\textit{Plan}_i + \delta_4 \textit{Sex}_i + \delta_5 \textit{BMI}_i + \delta_6 \textit{Ill}_i + b \textit{Auto}_i + \varepsilon_{3,i}, \end{align} $$

where Auto, Plan, Sex, BMI, Ill, and PA represent the automaticity, action plan of exercise, gender, body mass index, illness, and physical activity of the respondent, respectively.

To assess the normality assumption for the errors, we also compute the skewness and kurtosis of the two residuals, and then report them in Table 7. These values, together with the KS test, suggest a serious violation of the normality assumption for the $y-m,x$ regression residuals. Based on this, we also apply the proposed method to the dataset and then report the result in Table 8. First of all, the three methods produce positive indirect effects from 0.6600 to 0.7594. While for the CIs, only the LAD method shows insignificant outcome in the PRCT CI. At the same time, the Huber loss also yields the shortest CI among the three methods.

Table 7 Skewness and kurtosis of two regression residuals and the Kolmogorov–Smirnov test in action planning study

Table 8 The indirect effect estimates and their 95% CIs based on the LS, LAD, and Huber losses for the action planning study