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.

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_is (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_ib the relative indirect effect and ci′ 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²) to capture the overall effect of X on M from Equation 5 and construct the confidence interval for the product of this R² 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² shrinkage estimators for R² (Creedon et al., 2016), which are supposed to produce less biased estimates for R², and the non-parametric bootstrapping method for confidence intervals (R² × b). While there is no mathematical proof that the product of R² 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) R4.52 effect size measure. MacKinnon (2008) recommends three proportions of variance-explained measures termed R4.52, R4.62, and R4.72. 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):

where rYM2 is the squared correlation of Y and M. rYX2 is the squared correlation of Y and X. RY,MX2 is the squared multiple correlations of Y jointly explained by M and X. rMX2 is the squared correlation between M and X. rYM.X2 is the squared partial correlation of Y and M controlling for X. R4.52 is the explained variance in Y jointly by M and X. R4.62 is the proportion of variance in Y accounted for solely by M, weighted by the proportion of explained variance in M by X. R4.62 is difficult to interpret on an R² metric since it is the product of two proportions of variance from different variables (Preacher, 2011). R4.72 is R4.62 divided by RY,MX2, 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 R4.52 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):

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 R4.52. If there is no indirect effect, either β_b is zero (i.e., rYM=β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 R4.52; υ 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 (rYM−βaβc′)2 is the squared true correlation between M and Y. (RY,MX2-rYX2) is the difference between the total variance in Y explained by M and X (RY,MX2) and the total variance in Y explained solely by X (rYX2). Equation 14 is equal to Equation 13, which replaces (rYM−βaβc′)2 with the squared βb2 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 υ^=β^a2β^b2 where •^ is the sample estimator. Their simulation study found that it was upwardly biased. Based on the relationship that E(B2^)=B2+σB2, where E(•) is the expectation function, B is the unstandardized regression coefficient, and σB2 is the variance of B^, they propose the second estimator:

where σa2 and σb2 are the variances of â and b^, respectively. σ^X2and σ^Y2 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 had an acceptable coverage rate (i.e., between 92.5 and 97.5%) at N = 250.

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_ib equals βai2 βb2, where β_ai is the standardized a_i path and β_ai = a_i(σ_Di/σ_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 υ:

Following Equation 16, we chose (b^2-σb2)(σ^M2σ^Y2) in the hope of getting a less biased estimate of βb2. R~Y,MX2 and r~YX2 are the shrinkage estimators of RY,MX2 and rYX2, respectively. Following Creedon et al. (2016), we hoped that the shrinkage estimators would provide less biased estimates of RY,MX2 and rYX2. Therefore, all the components in Equation 14 (i.e., population υ) were adjusted for small-sample biases in Equation 17.

Table 2 summarizes the formulas of the shrinkage estimators to adjust R² and r² (i.e., R² 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 RY,MX2 and rYX2 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.

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:

In all conditions, the residual of the mediator M, e_M, was a normal variable with variance determined by RMX2, 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 RY,MX2 and rYX2 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, RY,MX2, and rYX2 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 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 ci′ 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² shrinkage method across a number of groups, group sizes, and sizes of a_i, b, and ci′ paths. ANOVAs and post hoc pairwise comparison results of each υ~ estimator are provided in the online supplements. All the R² 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² 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.

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₁ = 26), Hispanic (group 2; n₂ = 37), Afro-Hispanic (group 3; n₃ = 12), and others (group 4: n₄ = 15). Dummy coding was used and non-Hispanic white was the reference group. The mediator, perceived discrimination (M₁ = 5.04, SD₁ = 1.34; M₂ = 5.27, SD₂ = 0.93; M₃ = 5.25, SD₃ = 0.97; M₄ = 5.40, SD₄ = 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₁ = 77.44, SD₁ = 19.26; M₂ = 61.88, SD₂ = 35.57; M₃ = 61.90, SD₃ = 32.44; M₄ = 52.35, SD₄ = 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² shrinkage methods. Similar to simulation results, all R² 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].

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² shrinkage methods produced slightly less biased υ estimates. The R² shrinkage methods performed similarly on adjusting R~Y,MX2 and r~YX2. R² 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:

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.

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

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.

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.