Authors: John Kidd (1*Department of Biostatistics, Gillings School of Global Public Health, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina, U.S.A. .), Annie Green Howard (1*Department of Biostatistics, Gillings School of Global Public Health, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina, U.S.A. .; 2Carolina Population Center, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina, U.S.A. .), Heather M. Highland (3Department of Epidemiology, Gillings School of Global Public Health, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina, U.S.A. .), Penny Gordon-Larsen (2Carolina Population Center, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina, U.S.A. .; 4Department of Nutrition, Gillings School of Global Public Health, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina, U.S.A. .), Michael Patrick Bancks (5Department of Epidemiology and Prevention, Wake Forest School of Medicine, Winston-Salem, North Carolina, U.S.A..), Mercedes Carnethon (6Department of Preventive Medicine, Northwestern University, Chicago, Illinois, U.S.A..), Dan-Yu Lin (1*Department of Biostatistics, Gillings School of Global Public Health, University of North Carolina at Chapel Hill, Chapel Hill, North Carolina, U.S.A. .)
Categories: Article, Mediation pathway, Confidence intervals, Joint significance test, Mediation analysis, Missing data
Source: Statistical methods & applications
Authors: John Kidd, Annie Green Howard, Heather M. Highland, Penny Gordon-Larsen, Michael Patrick Bancks, Mercedes Carnethon, Dan-Yu Lin
Mediation analysis seeks to determine whether an independent variable affects a response directly or whether it does so indirectly, by way of a mediator or mediators. Scenarios that assume a single mediation are often overly simplistic, and analyses that include multiple mediators are becoming more common, particularly with the incorporation of high-dimensional data. Surprisingly, however, little attention has been given to multiple mediator and interaction effects. In this article, we propose new methods for testing the null hypothesis of no indirect effect with multiple mediators and interaction effects. We allow the estimators of the path effects to be possibly correlated; we also consider the practice of using confidence intervals to determine whether a mediation effect is zero. We compare the performance of our proposed method with existing methods through extensive simulation studies. Finally, we provide an application to data from the Coronary Artery Risk Development in Young Adults (CARDIA) study.
Mediation analysis seeks to determine, when assuming causal relationships, if an independent variable affects a response variable directly or whether it does so indirectly, through secondary variables called mediators (Alwin and Hauser, 1975). Mediation analysis is now used in a variety of situations, including scenarios with multiple mediators and with missing data (Lin et al, 2020; VanderWeele and Vansteelandt, 2013). Mediation analysis can also be used to discover how multiple mediators work together to create indirect effects, examining whether they act in ordered or unordered pathways (Taylor et al, 2008).
Detecting the indirect effect can be difficult (Biesanz et al, 2010; Huang and Pan, 2016; Vanderweele and Vansteelandt, 2009; Zhong et al, 2019). We consider the scenario where where the independent variable is independent of the response given all measured covariates, known as the no unmeasured confounding assumption (Cho and Huang, 2019; Imai et al, 2010; Vanderweele, 2011). In this scenario, the product of coefficients method posits that the indirect effect is the product of the effects of the path from the independent variable to the potential mediators, and then to the response (Alwin and Hauser, 1975; MacKinnon et al, 2007). When the product of the pathway effects is equal to zero, the null hypothesis of no indirect effect holds. This occurs for any scenario where at least one effect along the path is equal to zero. The asymptotic distribution of the estimator of the indirect effect depends on the number of effects that are truly zero (Kisbu-Sakarya et al, 2014; Wang, 2018), which is unknown.
Various methods for testing the null hypothesis that the indirect effect is zero have been proposed for the single-mediator scenario. The Sobel test utilizes the delta method to derive an approximate normal distribution of the product (Sobel, 1982). A joint significance test, called maxP, conducts tests of significance on the independent pathways comprising the indirect effect, using the maximum p-value of the two tests as the p-value for the overall hypothesis test (Cohen and Cohen, 1983; MacKinnon et al, 2002). Two recently developed methods, the PS-test and ASQ-test, use the cumulative probabilities of independent test statistics from the path effects and incorporate rejection regions to increase power over maxP (Kidd and Lin, 2023). In addition to hypothesis test-based methods, some methods obtain a confidence interval and determine whether the interval includes zero (Li et al, 2006; MacKinnon et al, 2004).
In this paper, we extend various continuous single-mediator methods for determining whether indirect effects are zero to scenarios with interaction effects and multiple mediators. We also propose a transformation approach for cases where nonzero correlations between the estimators of the path effects (which can occur in scenarios with missing data) violate the independence assumption of the maxP, PS-test, and ASQ-test. We demonstrate that the proposed methods have controlled type I error rates and compare their power to that of confidence interval approaches through extensive simulation studies. Finally, we apply our methods to data from the Coronary Artery Risk Development in Young Adults (CARDIA) study (Friedman et al, 1988).
Let Y denote the outcome and let G denote the exposure of interest. Let S denote the mediating variable of interest in the single-mediator setting (see Figure 1a), and S1 and S2 denote the two mediating variables of interest in unordered or ordered-mediators settings (Figures 1b and 1c).
In the single-mediator scenario, we relate S to G and additional variables through the linear regression model (1)S=αGG+αWTW+ε, where W is a vector that includes the unit component and any additional variables (e.g., ethnicity, age, sex), αG and αW are regression parameters, and ε is zero-mean normal with variance τ2. With α=αG,αWTT and X=(G,WT)T, equation (1) can then be expressed as (2)S=αTX+ε.
The outcome Y is related to S by the linear regression model (3)Y=βGG+βWTW+βSS+βSGSG+ϵ, where all β terms are regression parameters, and ϵ is zero-mean normal with variance σ2. With β=βG,βWT,βS,βSGT and Z=G,WT,S,SGT, equation (3) can then be expressed as (4)Y=βTZ+ϵ.
For an unordered-mediators scenario, we consider S=S1,S2T and relate G to S by the multivariate linear regression model (5)S=αGG+αWTW+ε, where vector αG=α1G,α2GT and matrix αW=α1W,α2W are regression parameters, and ε is bivariate normal with mean zero and covariance matrix τ12ρτ1τ2ρτ1τ2τ22. Define α to be the matrix with columns α1 and α2, where αj for j=1,2 is equal to αjG,αjWTT. Equation (5) can then be expressed as (6)S=αTX+ε.
The outcome Y is then related to the two unordered mediators by the linear regression model (7)Y=βGG+βWTW+βSTS+βSGTSG+ϵ, where βS=βS1,βS2T and βSG=βS1G,βS2GT are regression parameters, and ϵ is zero-mean normal with variance σ2. With β=βG,βWT,βST,βSGTT and Z=(G,WT,ST,STGT, equation (7) can be expressed as equation (4).
For two ordered mediators, we must consider three models. The first mediator, S1, is related to G by (8)S1=αGG+αWTW+ε1, where αG and αWT are as in the single-mediator scenario, and ε1 is zero-mean normal with variance τ12. With X1=G,WTT, we can express equation (8) as (9)S1=αTX1+ε1.
We then relate the second mediator, S2, to S1 using the linear regression model (10)S2=γGG+γWTW+γS1S1+γS1GS1G+ε2, where γG,γW,γS1, and γS1G are regression parameters and ε2 is zero-mean normal with variance τ22. Let γ=γG,γWT,γS1,γS1GT and X2=G,WT,S1,S1GT. Equation (10) can then be expressed as (11)S2=γTX2+ε2.
The outcome Y is then related to S2 by the linear regression model (12)Y=βGG+βWTW+βS1S1+βS1GS1G+βS2S2+βS2GS2G+ϵ, where the β terms are regression parameters and ϵ is as in the previous mediation scenarios. For all scenarios, W may contain different variables in the various equations. With β=βG,βWT,βS1,βS1G,βS2,βS2GT and Z=G,WT,S1,S1G,S2,S2GT, equation (12) can be expressed as equation (4).
In the various mediation scenarios, the indirect effect is denoted as δ, and is the product of the various path effects. For the single-mediator scenario, δ=αGβS+βSG. For the unordered-mediators scenario, the indirect effect through mediator j is denoted as δj, and δj=αjG(βSj+βSjG) for j=1,2. In the two ordered-mediators scenario, δ=αGγS1+γS1GβS2+βS2G, with simple extensions to three or more (Taylor et al, 2008). When interactions are included, the value of the indirect effect is specific to the level of the value of G.
Testing the null hypothesis H0:δ=0 (or δj=0) is complex, because the indirect effect will be equal to zero when a single term of the product is equal to zero. Thus, the null distribution includes all possibilities where one or more pathway terms are equal to zero.
Utilizing the delta method, the Sobel test is able to incorporate both interaction terms and multiple mediators. However, the Sobel test often has overly conservative type I error in the single mediator case. This problem is more pronounced for ordered-mediators scenarios, in which additional pathway terms may be equal to zero, making the asymptotic distribution derived by the delta method incorrect. Normal approximations of δ, although simple and commonly used, often have suboptimal performance.
Here, we extend joint significance-based methods, such as maxP, the PS-test, and the ASQ-test (Kidd and Lin, 2023). Where interaction effects are present, the tests of significance or cumulative probabilities must account for both the main and interaction effect (i.e., βS and βSG). Once the various significance tests have been performed and the cumulative probabilities found, the respective rejection criterias can be considered.
For mediation scenarios with multiple mediators, the extensions of the joint significance-based methods are different for unordered versus ordered-mediators scenarios. For unordered mediators, the indirect effect through a single mediator is found as in the single-mediator scenario. For an ordered-mediators scenario, three p-values or cumulative probabilities must be considered.
The extension for maxP simply uses the maximum of the p-values from the independent tests. For the ASQ-test, which compares the cumulative probability of two test statistics in the single-mediator scenario to a rejection region comprised of squares that align at diagonal corners (Kidd and Lin, 2023), the rejection region must be modified in order to maintain the correct type I error rate. Specifically, we use cubes with sides measuring α/2 beginning at the eight corners of the area defined by the cumulative probabilities and aligning at diagonal verticies. The cubes ascend (or descend) towards (0.5, 0.5, 0.5), or where each test statistic is equal to the null value (see Figure 2a). The extension is limited by the same compatibility concern as the single-mediator scenario and only certain significance levels can be tested, resulting in a p-value threshold.
As in the ASQ-test, the PS-test uses a rejection region within the space defined by the cumulative probabilities (Kidd and Lin, 2023). Whereas the single-mediator scenario finds the smallest value of α such that the rejection region comprised of a corner square and central-bands rejects the null hypothesis, the ordered-mediators scenario requires an additional dimension. The rejection region becomes cubes with sides of length α/2 in the eight corners, and central bands extending from the corner cubes towards the center of the space (see Figure 2b). The bands incorporate instances where the cumulative probabilities (or one minus any number of the cumulative probabilities) are all within α/4. As in the single-mediator scenario, the smallest value of α is found such that the test rejects the null hypothesis, and that value is the p-value for the test.
The corner cubes of the ASQ-test and PS-test represent where the null hypothesis would be rejected by maxP. The additional regions of the ASQ-test and PS-test ensure that they will be at least as powerful as maxP. These additional regions cause the null hypothesis to be rejected for cases where the cumulative probabilities of the path effects are very similar, but the signal is not strong enough to be found by maxP, such as when the sample size is small.
The ASQ-test, PS-test, and maxP all require that the terms of the product for δ be independent. Under the assumption of no unmeasured confounding (Imai et al, 2010), the error terms from the various linear regression models above are independent. This assumption is often accompanied by the assumption that the regression parameters from the different models are independent as well, but this may not be true, particularly when data is missing and alternative estimation approaches are used (Lin et al, 2020; Little and Rubin, 2002). This produces nonzero covariances, and joint significance-based tests are therefore inappropriate.
To address this limitation, we propose a transformation approach. Let θˆ be a vector containing the terms used to find the indirect effect; θˆ may contain separate terms for the interaction terms, or may contain a combination of terms (i.e., in the single-mediator scenario, θˆ may have separate elements for βS and βSG, or may have a single element for βS and βSG). In addition, let Σ denote the covariance matrix of θˆ, and let C denote the lower triangular Cholesky decomposition of Σ. We then find (13)θˆ*=C-1θˆ.
The components of θˆ* are independent with unit variances and can be used in the maxP, PS-test, and ASQ-test.
In addition to joint significance methods such as maxP, ASQ-test, and PS-test, significance in mediation analysis can also be determined using confidence intervals, with the indirect effect declared to be nonzero if the null value of δ does not lie within the confidence limits (Li et al, 2006; MacKinnon et al, 2004). While the Sobel test has an associated confidence interval, the test’s limitations make it less effective than methods such as Monte Carlo (MC) and bootstrap confidence intervals. MC confidence intervals utilize the known distributions of the estimators to generate many estimates of the path effects, which are then used to generate estimates of δˆ; the confidence limits are the appropriate quantiles (Preacher and Selig, 2012; Tofighi and MacKinnon, 2016). Because we generate the path estimates, no assumption about the distribution of δ is necessary. MC confidence intervals outperform other normal approximations, but the estimates of the effects require more complex data generation, and require knowledge of the multivariate distribution of the effects. In contrast, bootstrap confidence intervals find many estimates of δ and require no distributional assumptions (MacKinnon et al, 2004), but they are computationally expensive.
However, confidence intervals are not designed for hypothesis testing and do not fully address the complexity of the null distribution of the indirect effect. Though confidence intervals give additional information, hypothesis tests are more effective for determining whether an indirect effect is zero. We compare the performance of confidence intervals to hypothesis tests in the following section.
To compare the performance of the various methods described here, we conducted extensive simulation studies. For all mediation scenarios, Z is generated with five Z1=0,1, and2 , with probabilities p2,2p(1-p), and (1-p)2, respectively; Z2=1;Z3 is standard normal; Z4 is Bernoulli with a probability of success equal to 0.5; and Z5 is standard uniform. This setup represents a genetic scenario, with Z1 representing G as a SNP genotype under the Hardy-Weinberg equilibrium with a minor allele frequency p,Z2 representing the intercept term, Z3 representing a principal component for ancestry, Z4 representing sex, and Z5 representing the normalized age. We use Z3 to simulate population stratification, letting p=e0.5Z3/1+e0.5Z3, then set X,X1, and X2 to Z in the various mediation scenarios. We generate mediators and outcomes from equations (2–11) with intercept and error variance terms set to 1. Missing values are generated under a missing-completely-at-random scenario. Values of the sample size, path effects, missing percentages, and correlation (in the unordered mediation scenario) are varied. Additional simulations with Z1 representing a continuous independent variable yielded similar results and are omitted. Results are based on 20,000 replicates.
For the single-mediator scenario, we consider the power and type I error rates shown in Figure 3. Whereas the Sobel test has very conservative type I error rates, the PS-test, ASQ-test, and maxP have type I error rates much closer to the nominal significance level of α=0.05. The PS-test and ASQ-test have the highest power for all scenarios and are nearly indistinguishable. The PS-test, ASQ-test, and maxP are all more powerful than the confidence interval approaches, though the difference decreases as αG or βS+βSG increase. Whereas the power of all methods decreases as missingness increases (Figure 3f), the PS-test and ASQ-test retain more power than the other tested methods. The Sobel test is the most conservative method.
For two unordered mediators, we present the results for only one of the two mediation paths, as the type I error and power are nearly identical (Figure 4). All shown scenarios have the missingness of the two mediators set to 25%. Similar to the single-mediator scenario, the PS-test and ASQ-test are the most powerful with controlled type I error rates. For unordered mediators, the differences in power between the PS-test, ASQ-test, maxP, bootstrap confidence intervals, and MC confidence intervals decrease as the sample size increases (Figure 4d). The power of the methods can be limited by a single parameter value (Figure 4f). Again, the Sobel test is the most conservative, and it requires larger path effects to reach the same power as the other methods.
For two ordered mediators, all results have the missingness of the two mediators set to 25% (Figure 5). There is very little difference between the PS-test, ASQ-test, and maxP, and all three methods are more powerful than the confidence interval approaches and the Sobel test. Additionally, we see a similar increase in power in the plots for both γS1+γS1G and βS2+βS2G (Figures 5c and d).
In summary, none of the methods have inflated type I error; the PS-test and ASQ-test have the highest power, and the power of confidence intervals is generally lower than all the hypothesis testing methods except for the Sobel test, which is overly conservative for all scenarios. We therefore recommend the use of the PS-test and ASQ-test for testing whether the indirect effect is zero.
The CARDIA study is a prospective cohort study originally designed to examine cardiovascular disease risk factors in Black and white young adults. The study began in 1985–1986 and enrolled 5,115 Black and white men and women aged 18 to 30 from communities in Birmingham, Alabama; Chicago, Illinois; Minneapolis, Minnesota; and Oakland, California, followed over time. Over 35 years later, CARDIA continues to provide insights into various measurements of health and cardiovascular disease (Friedman et al, 1988).
Follow-up examinations have been completed at various time points, including 20, 25, and 30 years after enrollment. Various health-related measurements are collected at the examinations including physical activity (Sternfeld et al, 1999); lipid measurements, such as high and low density lipoprotein (HDL and LDL), triglycerides, and fasting glucose; and body mass index (BMI). We created and included other variables such as age, sex, race, community of enrollment, smoking status, caloric and alcohol intake, education level, and an a priori diet score (Friedman et al, 1988). Some participants were excluded from our those with diabetes by year 30, those on medication, or those with missing covariate or response measurements. As the ASQ-test does not return an exact p-value and performs similarly to the PS-test, we exclude it from the presented results.
We consider the potential indirect effect of physical activity on fasting glucose through BMI, HDL, LDL, the ratio of triglycerides over HDL, and the ratio of triglycerides over LDL. It is thought that physical activity affects future measurements of the mediators, which may influence future measurements of fasting glucose. We consider the effect of physical activity at year 20 on possible mediators at year 25, controlling for additional variables at year age, sex, race, community of enrollment, smoking status, caloric and alcohol intake, education level, and diet. We further consider the effect of this mediator at year 25 on fasting glucose at year 30, controlling for the same set of additional variables as measured at year 25. As early data collection and study drop-off may be associated with racial and educational groups, we repeat the analysis for all participants and within additional divisions (white, Black, post-high school education, no post-high school education, and a cross between these racial and educational levels). We consider each mediator in a single-mediator scenario, and each possible unordered and ordered pair in the corresponding mediation scenario. We find that the PS-test and maxP often yield similar results, with the PS-test never being more conservative than maxP (as seen in the Simulation Section) therefore, only results with differences between maxP and the PS-test are shared.
One example in the single-mediator scenario is when BMI is used as the mediator for the Black no post-high school education group. In this analysis,the PS-test yields smaller p-values than all other methods. With a sample size of167 and an estimated indirect effect of 0.002, the p-value of the maxP is 0.0910 and the p-value of the PS-test is 0.0445. The 95% bootstrap confidence interval is (−0.0009, 0.0070), and the 95% MC confidence interval is (−0.0015, 0.0074).
For the multiple mediators scenarios, we present only the unordered mediation results shown in Table 1; this is because scenarios where the PS-test yields smaller p-values than maxP are less common in ordered mediation scenarios (see Figure 5). We see that the maxP and PS-test are often similar, but in some scenarios, the PS-test yields smaller p-values than maxP, most noticeably in the white no post-high school education subgroup with BMI and HDL as the mediators. In the unordered scenario with BMI and the ratio of triglyceride over LDL, the PS-test finds a p-value of 0.0178 for the mediation through BMI, where maxP finds a p-value of 0.1463. In the scenario with BMI and the ratio of triglyceride over LDL, the ratio of triglyceride over LDL has a p-value of 0.0009 with the PS-test and a p-value of 0.0043 with maxP. As discussed above and seen in the simulation results shown in Figure 3, the power of the PS-test over other methods is improved in smaller sample sizes, such as the education level subgroup analyses.
In Table 1, we see additional scenarios where hypothesis tests provide stronger evidence than confidence intervals for determining whether an indirect effect is zero. For the white, no post-high school education subgroup with the mediators of BMI and HDL, the 95% confidence intervals for BMI include zero, but the PS-test returns a p-value less than 0.05. Confidence intervals may miss findings that could be found with the hypothesis testing methods.
Detecting an indirect effect in mediation analysis is a complex problem requiring careful approaches. We have considered methods for testing indirect effects in scenarios involving multiple mediators and possible interactions, and have proposed a method for handling correlations between the estimators of the path effects, such as those introduced by missing data approaches. We have demonstrated the increased power of our methods over other hypothesis testing methods and confidence intervals, and have applied the proposed and existing methods to a potential mediation scenario in the CARDIA study.
We have compared the performance of hypothesis testing methods to that of confidence interval methods in determining the significance of an indirect effect. Although we have shown that confidence intervals are less powerful for some mediation settings, they do provide additional information not yielded by significance tests alone, and they should not be discounted or ignored in mediation analysis. This is particularly relevant in the case of three or more sequential mediators, where joint significance methods have very low power. However, hypothesis testing procedures are preferable in our considered scenarios when the main objective is to determine whether an indirect effect is zero.
With the increased use of mediation analysis in large genetic studies, there has been increased research into improved methods for analyzing high-dimensional data (Dai et al, 2020; Huang, 2018). The PS-test and ASQ-test (Kidd and Lin, 2023) can be easily applied to some of these methods, such as those developed by Zhong et al.(Zhong et al, 2019), by converting the calculated p-values into cumulative probabilities and then utilizing the PS-test or ASQ-test. We are currently researching more direct applications of the methods proposed here to high-dimensional genetic data.