Effect size measure for mediation analysis with a multicategorical predictor

Many currently available effect size measures for mediation have limitations when the predictor is nominal with three or more categories. The mediation effect size measure υ was adopted for this situation. A simulation study was conducted to investigate the performance of its estimators. We manipulated several factors in data generation (number of groups, sample size per group, and effect sizes of paths) and effect size estimation [different R-squared (R2) shrinkage estimators]. Results showed that the Olkin–Pratt extended adjusted R2 estimator had the least bias and the smallest MSE in estimating υ across conditions. We also applied different estimators of υ in a real data example. Recommendations and guidelines were provided about the use of this estimator.


Introduction
Mediation analysis has enabled behavioral researchers to better understand the mechanistic relationships between variables. Many researchers are specifically interested in the role of mediators (M) that account for the relation between a predictor variable (X; independent variable) and an outcome variable (Y; dependent variable). A variety of effect size measures have been developed for mediation analysis. However, most of these effect size measures have limitations, including non-monotonicity and spurious inflation. Lachowicz et al. (2018) developed a new effect size measure, upsilon (υ), which has overcome these two limitations. The current study extends its work by applying this effect size metric to mediation models with a multicategorical predictor.
A simple mediation model shown in Figure 1A represents a mediation process in which a predictor X indirectly influences an outcome Y through a mediator M. This causal sequence suggests that X exerts a direct effect on M (a path), which, in turn, causally affects outcome Y (b path). X can have a direct effect on and Y (c ′ path), in which X directly influences Y. (1) This model can be estimated by a set of regression models or by structural equation modeling when the effects are linear and M and Y are treated as continuous. Equations 1, 2 are required to estimate the effects of the a path and b path, respectively. The indirect effect of X on Y is the product of a and b. The direct effect is c ′ (Eq 2). The total effect is c in Figure 1B and Eq 3, which equals the sum of X's direct and indirect effects on Y (Eq 4).
In this simple mediation model, predictor X can be dichotomous, or it can be treated as continuous. However, mediation analysis with a multicategorical predictor is common (Kalyanaraman and Sundar, 2006). When X is multicategorical, it can be expressed by applying coding strategies in regression analysis (Hayes and Preacher, 2014). When there are k groups comprising a multicategorical predictor X, (k − 1) coded variables are computed to represent each group. Different coding strategies are available, and the choice of strategy depends on the research question. Hayes and Preacher (2014) suggested using a dummy or contrast coding for a multicategorical predictor X in mediation analysis. In dummy coding, (k − 1) dummy variables (D i , i = 1, . . . , k − 1) are constructed, where D i is set to 1 to represent the cases in group i; otherwise, it is set to 0. The kth group is the reference group and is coded 0 in all D i s (see Table 1 for a three-group example). The effect of D i is the difference between the ith group and the reference group. In contrast coding, D i is constructed such that its effect represents the difference between the ith group and the average of the (i + 1)th group to the kth group (see Table 1 for a three-group example). A mediation model with a multicategorical predictor X is expressed in Figure 2A and Equations 5, 6.
Hayes termed a i b the relative indirect effect and c ′ i the relative direct effect. Equation 7 and Figure 2B show the relative total effect of X on Y, c i , which represents the sum of the corresponding relative direct and indirect effects (Equation 8). The effects are "relative" because they are based on comparisons between the group i and the other group (s), depending on the coding strategy of D i . For statistical inference of the multicategorical mediation analysis, we can use the non-parametric bootstrapping method, Monte Carlo, or product of moments methods to calculate the confidence intervals of relative indirect effects (Hayes and Preacher, 2014). Creedon et al. (2016) developed an omnibus test of the indirect effect of X on Y through M. They suggest using R-squared (R 2 ) to capture the overall effect of X on M from Equation 5 and construct the confidence interval for the product of this R 2 and b to provide an omnibus test of the indirect effect. They tested the performance of this method using Smith, Wherry-1, Wherry-2, Olkin-Pratt, Pratt, and Claudy R 2 shrinkage estimators for R 2 (Creedon et al., 2016), which are supposed to produce less biased estimates for R 2 , and the non-parametric bootstrapping method for .
/fpsyg. .  N is sample size. p is number of predictors. R 2 is the sample squared multiple correlation. r is the sample correlation.

Shrinkage estimator Formula
confidence intervals (R 2 × b). While there is no mathematical proof that the product of R 2 and b quantifies the overall mediation effect properly, this method has shown satisfactory results in maintaining the type I error rate at a set level (Creedon et al., 2016). In addition to statistical inferences, reporting on effect sizes is highly encouraged or mandated by many peer-reviewed journals and organizations, including the American American Psychological Association (2010, p. 33) and the International Committee of Medical Journal Editors via the Consolidated Standard of Reporting Trials (CONSORT; Moher et al., 2010). Reporting on the effect size of an indirect effect is challenging because mediation analysis is a non-standard regression-based analysis, such that standardized mean differences, correlation coefficients, and proportions of variance explained are effect size metrics that are not sufficient to capture the entirety of an indirect effect (Lachowicz et al., 2018). Preacher (2011) have listed the desiderata for a good effect size measure: (a) An effect size should have an interpretable scale, (b) its confidence interval can be calculated, and (c) its sample estimation of the population parameter should be unbiased, consistent, and efficient. Wen (2015) have since added a desideratum that (d) the function of the effect size measure should be a monotonic representation, either in raw or absolute form, of the quantified effect.
In this study, we focus on a mediation effect size measure υ developed by Lachowicz et al. (2018), which is a modification of MacKinnon's (2008) R 2 4.5 effect size measure. MacKinnon (2008) recommends three proportions of variance-explained measures termed R 2 4.5 , R 2 4.6 , and R 2 4.7 . The subscripts (4.5), (4.6), and (4.7) are based on the original equation numbers in MacKinnon (2008), and these notations have continued to be referenced in subsequent literature (e.g., Lachowicz et al., 2018;Preacher, 2011). In the simple mediation model in Figure 1A, these effect sizes are calculated as follows (Equations 9-11): R 2 4.6 = (r 2 MX ) r 2 YM.X where r 2 YM is the squared correlation of Y and M. r 2 YX is the squared correlation of Y and X. R 2 Y,MX is the squared multiple correlations of Y jointly explained by M and X. r 2 MX is the squared correlation between M and X. r 2 YM.X is the squared partial correlation of Y and M controlling for X. R 2 4.5 is the explained variance in Y jointly by M and X. R 2 4.6 is the proportion of variance in Y accounted for solely by M, weighted by the proportion of explained variance in M by X. R 2 4.6 is difficult to interpret on an R 2 metric since it is the product of two proportions of variance from different variables (Preacher, 2011). R 2 4.7 is R 2 4.6 divided by R 2 Y,MX , which is interpreted as the proportion of explained variance in Y jointly explained by M and X.
Among the three effect size measures mentioned previously, Lachowicz et al. (2018) modified R 2 4.5 and proposed a new measure υ to address the problem of spurious inflation. To illustrate the issue of spurious inflation, one needs to consider the elements of the simple mediation model in Equations 1-3. It is assumed that Y is dependent on M, and Y and M are mutually dependent on X, given that all other assumptions hold (temporal precedence, covariation of the cause and effect, etc.,). The zero-order correlation between M and Y has two components: (a) the conditional correlation between M and Y independent of X, and (b) the correlation due to the mutual dependence of Y and M on X. The first component is often referred to as true correlation, and the latter component is considered spuriously inflated. A true correlation between M and Y is needed to create a mediation effect size without spurious inflation (Lachowicz et al., 2018). Lachowicz et al. decomposed the correlation between Y and M (r YM ): Frontiers in Psychology frontiersin.org . /fpsyg. . where β a is the standardized a path. β c ′ is the standardized c ′ path. β b is the standardized b path, which captures the true correlation between M and Y in the mediation model. Lachowicz et al. pointed out that β a β c ′ is the component that is spuriously inflated, indicating this is a key limitation of R 2 4.5 . If there is no indirect effect, either β b is zero (i.e., r YM = β a β c ′ ), or the inflated term, β a β c ′ , is zero (i.e., β a = 0; r YM = β b ). If there is no direct effect, the spuriously inflated term is zero (i.e., r YM = β b ). As a result, r YM cannot distinguish whether an indirect effect is present. Lachowicz et al. (2018) developed the effect size υ by removing the spurious inflation from R 2 4.5 ; υ measures the variance in Y explained jointly by M and X, correcting for the spurious inflation between M and Y on X (Equation 13). The term (r YM − β a β c ′ ) 2 is the squared true correlation between M and Y. R 2 Y,MX − r 2 YX is the difference between the total variance in Y explained by M and X (R 2 Y,MX ) and the total variance in Y explained solely by X (r 2 YX ). Equation 14 is equal to Equation 13, which replaces (r YM − β a β c ′ ) 2 with the squared β 2 b from Equation 12. Equation 15 is also equivalent to Equations 13 and 14 (Lachowicz et al., 2018).
According to Lachowicz et al., υ has numerous desirable properties as an effect size measure of the indirect effect: (a) It is standardized (scale-free), (b) it is independent of sample size, (c) its function in the absolute value of the indirect effect is monotonic, and (d) its confidence interval can be constructed. Using Equation 15, Lachowicz et al. proposed sample estimators of υ for the simple mediation model ( Figure 1A). The first estimator isυ = β 2 aβ 2 b where• is the sample estimator. Their simulation study found that it was upwardly biased. Based on the relationship that is the expectation function, B is the unstandardized regression coefficient, and σ 2 B is the variance ofB, they propose the second estimator: where σ 2 a and σ 2 b are the variances ofâ andb, respectively.σ 2 X and σ 2 Y are the sample variances of X and Y. The confidence interval of υ can be calculated via non-parametric bootstrapping. Their simulation study found that this estimator had an acceptable bias with only four experimental conditions resulting in percent relative biases > 5% at N = 100. The confidence interval of υ can be constructed using this estimator via non-parametric bootstrapping; their simulation showed that the bootstrapped confidence interval .

Current study
Because of the desirable properties of υ, we applied it in a simple mediation model with a multicategorical predictor ( Figure 2A). We conducted a simulation study to investigate the performance of a sample estimator of υ. Using Equation 15, the υ for each relative indirect effect a i b equals β 2 a i β 2 b , where β a i is the standardized a i path and β a i = a i σ D i /σ M . In addition, we believe that researchers would be interested in calculating υ for the overall indirect effect of X on Y through M. Equation 15 and its sample estimators are difficult to apply in this situation because of multiple a i paths; thus, we chose Equation 14 and proposed an estimator of υ:υ . Table 2 summarizes the formulas of the shrinkage estimators to adjust R 2 and r 2 (i.e., R 2 in single-predictor regression). Based on the results from previous simulation studies on the performance of the shrinkage estimators (Raju et al., 1999;Yin, 2001;Walker, 2007;Wang and Thompson, 2007;Shieh, 2008;Creedon et al., 2016), Pratt, Ezekiel, Smith, Wherry-2, Walker, and Olkin-Pratt extended formulas performed well in at least one study. Nevertheless, none of these studies directly tested the mediation effect size measures. It is difficult to determine which formula is the best option for the sample adjustments of R 2 Y,MX and r 2 YX when estimating υ. We conducted a simulation study to examine the following shrinkage estimators: Claudy, Ezekiel, Olkin-Pratt, Olkin-Pratt extended, Pratt, Smith, Walker, and Wherry.

Methods
A simulation study was conducted to evaluate the performance of sample estimator υ in Equation 17 under finite samples. The simulation was based on the mediation model in Figure 2A. Five factors were manipulated: (1) A number of groups in X: k = 3, 4, 5. The conditions followed those in Creedon et al. (2016).
. /fpsyg. . (2) Sample size per group, n = 10, 20, 25, 50, 100. The range of n covered small-to-large sample sizes across groups. (3) Effect size of a i paths were manipulated as the mean difference between adjacent groups on M. The mean difference was set in terms of Cohen's d = 0, 0.2, 0.5, 0.8, representing null, small, medium, and large effects, respectively (Cohen, 1992). In all conditions, the residual of the mediator M, e M , was a normal variable with variance determined by R 2 MX , and the residual of the outcome Y, e Y , was a standard normal variable. The simulation used a full factorial design with a total of 720 conditions (3 × 5 × 4 × 4 × 3). For each condition, 10,000 replications were created. For each condition and replication, nine R 2 Y,MX and r 2 YX estimators were used to estimate υ: unadjusted, Claudy, Ezekiel, Olkin-Pratt, Olkin-Pratt extended, Pratt, Smith, Walker, and Wherry. The unadjusted sample estimates of υ were calculated using Equation 14, in which β b , R 2 Y,MX , and r 2 YX were not adjusted. The rest of the estimators were calculated using Equation 17. The 95% confidence interval of υ was constructed using non-parametric bootstrapping with 1,000 bootstrap samples. The simulation study was conducted in R (Version 3.5.3; Windows system), and the packages boot 1.3-20 (Canty, 2017), dummies 1.5.6 (Brown, 2012), effsize 0.7.4 (Torchiano, 2018), and lm.beta 1.5-1 (Behrendt, 2014) were utilized.
To evaluate the performance of sample estimators of υ, the following outcomes were used: bias, standardized bias, mean squared error (MSE), and coverage rate. For any parameter θ , bias is the difference between the expectation of the sample estimatesθ and the parameter (Equation 18).
In each condition, the population value of υ was calculated using Equation 14. The online Supplementary material present the population value of υ in each simulation condition.
Standardized bias is the bias divided by the standard deviation of the sample estimatesθ (Equation 19). Standardized bias within ±0.4 can be regarded as acceptable (Collins et al., 2001).
Mean squared error is the expected squared difference between the sample estimates and the parameter. It is equal to the sum of the variance of sample estimates and squared bias (Equation 20). An unbiased sample estimator would produce an MSE equal to the variance of the estimator.
Coverage rate is the proportion of the sample in which the population parameter is contained within the 95% confidence interval across replications in a condition. An acceptable coverage rate is between 92.5 and 97.5% (Bradley, 1978).
We expected that the unadjustedυ would produce upwardly biased estimates, particularly for small sample size (N) and small effect size conditions. The bias of unadjustedυ would decrease with increasing N and effect size. For the shrinkage adjustedυ, it was expected that the estimates would have acceptable bias and .
/fpsyg. . that the changes in N would not change the bias. We expected that the MSE of the unadjustedυ would decrease as N increased since larger sample sizes decrease both sampling error and bias. For the same reason, we expected that the MSE of the shrinkage adjustedυ would also decrease as N increased. Finally, for both unadjustedυ and shrinkage adjustedυ, we expected that the coverage rate would approach the acceptable range as N increased. We conducted one-way analyses of variance (ANOVAs) to study the effects of each manipulated factor (number of groups, sample size per group, effect size of a i paths, size of b path, and effect size of c ′ i paths) on the bias, standardized bias, MSE, and coverage rate. For each significant factor, we conducted post-hoc pairwise comparisons with Tukey's HSD tests and produced boxplots by different conditions within the factor. Table 3 shows the bias, standardized bias, MSE, and coverage rate ofυ using different R 2 shrinkage method across a number of groups, group sizes, and sizes of a i , b, and c ′ i paths. ANOVAs and post hoc pairwise comparison results of eachυ estimator are provided in the online supplements. All the R 2 shrinkage methods performed similarly. The Olkin-Pratt extended method had the least bias and the lowest MSE. The Smith method yielded the least standardized bias. However, the unadjusted method and all the R 2 shrinkage methods produced very large standardized biases and were beyond the acceptable range (>0.4). None of the estimators produced satisfactory coverage rates (>0.9) across conditions. The unadjustedυ had the highest coverage rate among all the sample estimators. This finding was consistent with Lachowicz et al. (2018). Based on the results, we decided to focus on the Olkin-Pratt extended shrinkage method (υ OPE ) hereafter. The online supplements provide the results of ANOVAs and post hoc pairwise comparisons of each manipulated factor on the bias, standardized bias, MSE, and coverage rate ofυ OPE .

Bias
The size of the b path had the strongest effect on the bias of υ OPE , F(3, 716) = 448.22, p < 0.0001; η 2 = 0.65. Figure 3A shows the boxplots of bias ofυ OPE by different b path conditions. The bias ofυ OPE at b = 0.59 was significantly higher than all the other b path conditions. The bias at b = 0.39 was also significantly higher than the b = 0.15 and b = 0 conditions. A higher effect of b path was associated with a more positive biasυ OPE .
The effect size of a i paths had the second strongest effect on the bias ofυ OPE , F(3, 716) = 39.13, p < 0.0001, η 2 = 0.14. Figure 3B shows the boxplots of bias ofυ OPE by different effect sizes of a i paths. All pairwise comparisons were significant (ps < 0.05) except (d = 0.2 vs. d = 0). The bias ofυ OPE at d = 0.8 was significantly lower than d = 0.5, 0.2, and 0. Bias ofυ OPE was the lowest when d = 0.5. The effect size of c ′ i paths were also affected by the bias of υ OPE with F(2, 717) = 8.26, p < 0.001 η 2 = 0.02. Figure 3C shows that the larger effect sizes of c ′ i paths had a larger bias ofυ OPE . The number of groups had a small effect on the bias ofυ OPE , F(2, 717) = 3.20, p < 0.05, η 2 = 0.009 ( Figure 3D). The bias ofυ OPE at groups of k = 5 was significantly smaller than that at k = 3, and the bias was the lowest at k = 5. Group size had no significant effect on the bias ofυ OPE , F(4, 715) = 0.57, p = 0.69, η 2 = 0.003.

Standardized bias
The size of the b path had an effect on the standardized bias of υ OPE , F(3, 716) = 25.73, p < 0.0001, η 2 = 0.10.υ OPE had the most positive standardized bias when b = 0.59 with ( Figure 4A). When b = 0.39,υ OPE had more than 50% of cases with positive standardized bias > 0.4. When b = 0.15 and 0, the median standardized bias was within the ±0.4 criterion, yet some cases were beyond this range. The effect size of a i paths had an effect on the standardized bias of υ OPE , F(3, 716) = 19.98, p < 0.0001, η 2 = 0.08. Figure 4B shows that the standardized bias ofυ OPE decreased when the effect size of a i paths increased. When d = 0 and 0.2, the standardized bias of υ OPE was > 0.4. When d = 0.5, the median standardized bias was < 0.4. When d = 0.8, the standardized bias ofυ OPE was the lowest and the median standardized bias was close to 0. The effect size of c ′ i paths did not have a significant effect on the standardized bias of υ OPE , F(2, 717) = 0.25, p = 0.78, η 2 = 0.001.
The number of groups had no significant effect on the standardized bias ofυ OPE , F(2, 717) = 2.84, p = 0.06, η 2 = 0.008. Group size had a significant effect on the standardized bias ofυ OPE , F(4, 715) = 7.66, p < 0.0001, η 2 = 0.04. Figure 4C shows that the standardized bias decreased when group size increased. The standardized bias ofυ OPE was the most positive when n = 10 (median > 0.4), and its standardized bias was significantly more positive than those in other conditions (ps < 0.05). The median standardized bias was >0.4 when group size = 25 and 100. When group size = 100 and 200, the median standardized bias was <0.4.

Mean squared error
The size of the b path had a significant and large effect on the MSE ofυ OPE , F(3, 716) = 301.85, p < 0.0001, η 2 = 0.56. MSE of υ OPE increased when b increased ( Figure 5A). The MSE at b = 0.59 was significantly higher than those in other conditions (b = 0.39, 0.15, 0; ps < 0.05), and the MSE at b = 0.39 was significantly higher than those at b = 0.15 and 0 (ps < 0.05). When b = 0 and 0.15, the MSE were close to zero. The effect sizes of a i and c  Figure 5B supports the hypothesis that the MSE ofυ OPE would decrease as group size increased in general. MSE ofυ OPE was the highest when group size = 10, and its MSE was significantly higher than those in other conditions (ps < 0.05). When group size = 25, 50, and 100, the median MSEs were close to zero.

Coverage rate
The size of the b path had a significant effect on the coverage rate ofυ OPE , F(3, 716) = 99.18, p < 0.0001, η 2 = 0.29. All pairwise comparisons were significant (ps < 0.05). Figure 6A shows that the . /fpsyg. . coverage rate decreased when b increased, and only when b = 0 didυ OPE reach satisfactory coverage rate (>0.9). The effect size of a i paths had a significant effect on the coverage rate ofυ OPE , F(3, 716) = 2.24, p < 0.0001, η 2 = 0.08. Figure 6B shows that coverage rate at d = 0 had a significantly lower coverage rate than all other conditions (d = 0.2, 0.5, 0.8; ps < 0.05). For d = 0.2, 0.5, and 0.8 conditions,υ OPE did not have satisfactory coverage rates (<0.9). The effect size of c ′ i paths did not have a significant effect on the coverage rate ofυ OPE , F(2, 717) = 0.01, p = 0.99, η 2 < 0.0001. None of the conditions for c ′ i paths reached a majority of satisfactory coverage rates.
The number of groups had no significant effect on the coverage rate ofυ OPE , F(2, 717) = 1.83, p = 0.16, η 2 = 0.005. Group size had a significant effect on the coverage rate ofυ OPE , F(4, 715) = 36.26, p < 0.0001, η 2 = 0.17. Figure 6C shows that the coverage rate decreased when the group size increased. The coverage rate at group size = 10 was the highest and was significantly higher than those in other conditions (ps < 0.05). Group size = 10 reached the satisfactory cutoff with a mean coverage rate of 0.91. The coverage was the lowest at group size = 200.
Estimating υ without finite sample adjustment to β b Results suggest that the performance ofυ OPE was influenced by the size of the b path. According to Equation 16, the finite sample adjustment ofυ has two parts: (1) Adjusting β b and (2) adjusting R 2 Y,MX and r 2 YX . R 2 Y,MX and r 2 YX were unlikely to be influenced by the size of b path. Therefore, we conducted additional analyses which estimated υ using Equation 14 without finite sample adjustment to β b and with the Olkin-Pratt extended shrinkage method for R 2 Y,MX and r 2 YX . The results are provided in the online supplements. Similar to the previous results of theυ OPE with adjusted β b (Equation 16), the size of the b path had significant effects on the bias, standardized bias, and MSE of this estimator. We conclude that this estimator did not further improveυ OPE .

Empirical example
We present an empirical example to help facilitate the interpretation of the υ effect size measure with a multicategorical predictor. The data were acquired through the NIMH-funded (PI: Rivera Mindt; K23MH079718) Medication Adherence Study at the Icahn School of Medicine at Mount Sinai (ISMMS) in New York City. We were interested in the mediation process between race/ethnicity and antiretroviral medication adherence via perceived racial and ethnic discrimination. There were 90 participants with completed measures. Race/ethnicity was a multicategorical variable with four groups: non-Hispanic white (group 1; n 1 = 26), Hispanic (group 2; n 2 = 37), Afro-Hispanic (group 3; n 3 = 12), and others (group 4: n 4 = 15). Dummy coding was used and non-Hispanic white was the reference group. The mediator, perceived discrimination (M 1 = 5.04, SD 1 = 1.34; M 2 = 5.27, SD 2 = 0.93; M 3 = 5.25, SD 3 = 0.97; M 4 = 5.40, SD 4 = 0.91), was measured by the Perceived Ethnic Discrimination Questionnaire-Community Version Brief (PEDQ-CVB, scaled from 1 to 5; Brondolo et al., 2005). The outcome variable, medication adherence (M 1 = 77.44, SD 1 = 19.26; M 2 = 61.88, SD 2 = 35.57; M 3 = 61.90, SD 3 = 32.44; M 4 = 52.35, SD 4 = 32.57), was calculated as the percentage (0 to 100%) of doses taken on schedule as assessed by the Medication Event Monitoring System (MEMS; Group AARDEX, 2022). Table 4 presents the results for unadjustedυs and finite sample adjustedυs by correctingβ b and different R 2 shrinkage methods. Similar to simulation results, all R 2 shrinkage methods were performed equally, withυ ranging from 0.146 to 0.147, meaning the variance of medication adherence explained by race/ethnicity via perceived discrimination was equal to 0.15. All adjustedυs had similar 95% bootstrapped confidence intervals [−0.23, 1.55].

Discussion
The goal of this study was to extend the mediation effect size υ developed by Lachowicz et al. (2018) to mediation models involving a multicategorical predictor. Theoretically, υ has many desirable properties as a mediation effect size, particularly for its monotonicity, and bootstrapping can be used to construct its confidence interval. Their simulation showed that the unadjusted and finite sample adjusted sample estimators were consistent effect size measures. Furthermore, this effect size measure is standardized with an invariance of a linear transformation, so it is independent of the predictor scales, the mediator, and the outcome variable. We applied υ to a mediation model with a multicategorical predictor. In this scenario, our simulation results showed that υ did not retain some of the desiderata asserted by Lachowicz et al. (2018). Based on our results, the size of b path was the most important factor that negatively influenced the accuracy of the sample estimator of υ. In our data generation process, there were no large differences between the effect sizes of the a path and those of the b path. On further scrutiny, the performance of the sample estimator with uncorrectedβ b was similar to the sample estimator with theβ b correction. Therefore, the large effect of b path remained undiscovered. R 2 shrinkage methods produced slightly less biased υ estimates. The R 2 shrinkage methods performed similarly on adjustingR 2 Y,MX andr 2 YX . R 2 shrinkage methods might not be the source of high values of bias, standardized bias, and MSE.
Sinceυ OPE achieved higher performance among all the sample estimators of υ across all the experimental conditions, it is recommended that researchers useυ OPE for simple mediation models involving a multicategorical predictor, as illustrated earlier. However, researchers should be cautious about the scenarios mentioned below: (1) When the size of b path reaches 0.39-0.59,υ OPE could become an upwardly biased effect size measure. A large b path (e.g., b = 0.59) could also result in a large value of MSE inυ OPE .
(2)υ OPE does not perform well and will be positively biased when a i paths have small effects (i.e., Cohen's d close or equal to 0). (3) Researchers should pay attention to small group sizes.υ OPE could have high standardized bias and high MSE when n = 10.
There are several limitations to this study. First, further analyses should be conducted on the relationship between the effect of the size of b path and the bias of the sample estimator of υ. Second, since this effect size measure is appropriate for the mediation model specified in this study, future studies should further develop this effect size measure to be suitable for more complex mediation models, such as models with multicategorical mediators or outcomes, models with moderations or latent variables, and models with multilevel or longitudinal data structures. Finally, since there are more rigorous assumptions that need to be made to justify an indirect effect as a causal effect, it is possible for researchers to introduce unknown bias into the estimation of υ. Future studies should investigate the performance of the estimator when such moderate violations of assumptions occur.

Data availability statement
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics statement
The studies involving human participants were reviewed and approved by NIMH-funded study (Grant# K23MH079718; Principal Investigator: M. Rivera Mindt) at the Icahn School of Medicine at Mount Sinai (ISMMS) https://reporter.nih.gov/ search/XK3XEZYfbEuZQESg5r-YNA/project-details/7681039. The patients/participants provided their written informed consent to participate in this study.