Sensitivity Plots for Confounder Bias in the Single Mediator Model

Abstract

Background:

Causal inference continues to be a critical aspect of evaluation research. Recent research in causal inference for statistical mediation has focused on addressing the sequential ignorability assumption; specifically, that there is no confounding between the mediator and the outcome variable.

Objectives:

This article compares and contrasts three different methods for assessing sensitivity to confounding and describes the graphical depiction of these methods.

Design:

Two types of data were used to fully examine the plots for sensitivity analysis. The first type was generated data from a single mediator model with a confounder influencing both the mediator and the outcome variable. The second was data from an actual intervention study. With both types of data, situations are examined where confounding has a large effect and a small effect.

Subjects:

The nonsimulated data were from a large intervention study to decrease intentions to use steroids among high school football players. We demonstrate one situation where confounding is likely and another situation where confounding is unlikely.

Conclusions:

We discuss how these methods could be implemented in future mediation studies as well as the limitations and future directions for these methods.

Keywords

mediation indirect effects causal inference confounder bias sensitivity analysis

In 1986, Baron and Kenny published their seminal article on distinguishing and testing for moderators and mediators (Baron and Kenny 1986). Although most of the methods in that article were not new (cf. James and Brett 1984; Judd and Kenny 1981; Hyman 1955), it provided a clear framework for testing mediation and moderation within the social sciences. Baron and Kenny (1986, 1176) defined a mediator as any variable which “accounts for the relation between the predictor and the criterion.” Theory-driven interventions highlight several causal mechanisms that are hypothesized to change via the intervention, which then change the outcome variable of interest. The utility in testing for mediating variables is that the magnitude of the effect of the mediating variables can be used to improve interventions to target the most influential mediating variables and remove or reduce strategies to change other mediating variables that did not significantly affect the outcome. This theoretical approach to intervention development should improve the effectiveness and efficiency of subsequent interventions (MacKinnon and Dwyer 1993).

Randomized experiments are often considered a necessary component when making any causal inference. However, randomized experiments focusing on group differences in outcomes generally do not reveal the causal mediating mechanisms of the phenomenon being studied. The goal of mediation analysis is to investigate the underlying mechanisms of an intervention. Evaluation researchers use randomized experiments to determine if an evaluation program affects the hypothesized mediator and if the mediator affects the targeted outcome. However, work on causal inference in mediation models shows that randomizing the treatment alone is not adequate to reach accurate causal conclusions about mediating processes (Holland 1988; Robins and Greenland 1992). Therefore, it is critical to consider causal estimation issues when testing mediating mechanisms. The purposes of this article are to describe causal estimation issues in the current use of mediation analysis and to illustrate several sensitivity analyses for investigating confounder bias applied to mediation analysis of an anabolic steroid prevention program.

Linear Regression Approach to Mediation

The conventional mediation analysis estimates the parameters of the following three equations:

Y = i_{1} + c X + e_{1}

Y = i_{2} + c^{'} X + b M + h X M + e_{2}

M = i_{3} + a X + e_{3},

where Y is the dependent variable, M is the mediator, and X is the independent variable (MacKinnon 2008; see Figure 1). Note that Y and M are treated as continuous normally distributed variables and linear relations were assumed between variables. The coefficient c represents the total effect of X on Y; c' is the coefficient relating X to Y adjusted for the effects of M (i.e., the direct effect of X on Y); b is the coefficient relating the mediator to the dependent variable adjusted for the effects of the independent variable; a is the coefficient relating the independent variable to the mediating variable; e ₁, e ₂, and e ₃ represent unexplained or error variances, which are assumed to be independent across equations; i ₁, i ₂, and i ₃ are the intercepts; and the coefficient h in Equation 2 reflects the interaction effect of X and M on Y.

Figure 1.

Single mediator model.

The two commonly used approaches to test for mediation are Baron and Kenny’s (1986) test of mediation (also called a causal steps test; MacKinnon et al. 2002) and the product of coefficients method. In the Baron and Kenny method, a series of statistical tests are conducted to decide if the data are consistent with mediation. The second method is to calculate the point estimate of the mediated effect as the product of coefficients, ab; and then obtain asymmetric confidence intervals (MacKinnon 2008). Then, it is concluded whether the relation between X and Y is mediated through M, indicating a causal chain if the confidence interval does not contain zero.

Potential Outcomes Framework

The potential outcomes approach to treatment effects (Holland 1986, 1988; Rubin 1974) makes a distinction between an individual’s observed and counterfactual outcomes. Let variable X be an evaluation program with level x (x = 1 for the treatment, x = 0 for the control) and variable Y the outcome variable. This framework considers all possible conditions that an individual could serve including both treatment and control conditions, even though for observed data an individual can only serve in one of the two groups. If an individual is assigned to the treatment group, the potential outcome, Y(1), is in fact the observed outcome Y. The counterfactual outcome for that individual is the value he or she would obtain on the outcome variable if he or she had been assigned to the control group, Y(0). The corresponding individual causal effect is then equal to Y(1) − Y(0). Because it is not possible for the same person to appear in both groups at the same time (referred as the “fundamental problem of causal inference” by Holland 1986), the individual causal effect is not estimable but it is possible to compute a causal effect averaged across individuals in each group. The average causal effect computed as E[Y(1) − Y(0)] represents the difference in average Y between the groups and is a causal estimator under assumptions, primarily, that individuals have been randomized to the two conditions.

Now suppose a potential mediating variable M with level m mediates the relation between X and Y. Including M leads to the formulation of the following effects: controlled direct effect (CDE), natural direct effect (NDE), and natural indirect effect (NIE; Pearl 2001; Robins and Greenland 1992). Let Y(x, m) denote the potential outcome for an individual under the intervention level x and mediator level m. X takes the value x = 0 for the control group and x = 1 for the treatment group. If the actual value of M is m for an individual, then the counterfactual value of M for that individual is denoted as m′.

The average CDE of treatment X on outcome Y is then the direct effect of treatment on the outcome at a fixed level of the mediator at m:

C D E = Y (1, m) - Y (0, m) .

The average NDE of X on Y is different from the average CDE in that M is set to the level M(x), the level that would have naturally occurred under one of the conditions of X. For instance, in the case of M(0), the NDE is the effect of treatment X on outcome Y when X did not influence the mediator M (or the participants were assigned the mediator level under the control condition).

N D E = Y (1, M (0)) - Y (0, M (0)) .

Following the same approach, the average NIE is the effect of the treatment on outcome when changing the level of M. In other words, the NIE is the effect of treatment on outcome when the level of M is changed when X is set to a certain value (0 in this case).

N I E = Y (0, M (1)) - Y (0, M (0)) .

The sum of the NDE and NIE equals the total effect. However, if there is an interaction between X and M, there are two NDEs; one at the mediator level for the control group (pure NDE) and another at the mediator level for the treatment group (total NDE). Pearl (2001) demonstrates that the decomposition of total effect into the NDE and NIE holds even in models with interactions (i.e., the pure NDE plus the total NIE equals the total effect). Although the NIE outlined previously allows for interactions between X and M, the Baron and Kenny (1986) approach typically assumes no interaction between the two variables.

Assumptions for Identification of the Direct and Indirect Effects

The following two assumptions are required for CDE to be identified (Pearl 2001; VanderWeele 2010):

i. No unmeasured confounder for the relation between X and Y; and

ii. No unmeasured confounder for the relation between M and Y.

The following two assumptions are required in addition to the two assumptions mentioned above for NDE and NIE to be identified:

iii. No unmeasured confounder for the relation between X and M; and

iv. No measured or unmeasured M to Y confounder affected by treatment.

Assumptions (i) and (iii) refer to the ignorability of treatment assignment (i.e., treatment assignment is independent of potential outcomes for the mediator and outcome), given observed pretreatment confounders. This assumption is usually satisfied with randomization of X.

Assumption (ii) refers to the ignorability of the mediator and the potential outcomes for Y, given the observed treatment and pretreatment confounders. This assumption is hard to meet because randomization of M is usually not plausible for many studies (i.e., the mediator status is not randomly assigned but rather self-selected by participants). Even conditioning on observed confounders for the relation between M and Y, there can still be unobserved confounders. This assumption is strong and is usually ignored in studies even though a causal interpretation of a mediated effect is not possible without its existence (MacKinnon, Kisbu-Sakarya, and Gottschall 2013; MacKinnon 2008, chapter 13). Assumption (iv) is made by the methods described in this article and methods to investigate this assumption are only starting to be developed (Coffman and Zhong 2012).

To summarize, the linear regression approach, as well as several other potential outcome methods including the NDE and NIE, assumes sequential ignorability which consists of the ignorability of the treatment assignment and the ignorability of the mediator. In other words, with successful random assignment of X, c and a represent causal effects (with some assumptions) but b and c' do not have a causal interpretation without further assumptions. For instance, given that the individuals are randomly assigned to the evaluation program, there should be no confounders of the X to Y, and X to M relations. However, since participants are not randomly assigned to values of M, the relation between M and Y cannot be considered causal since there can be one or more confounders that can account for the effect of M on Y.

Sensitivity Analyses to Assess Confounder Bias

Given the causal inference challenges in mediation analysis, researchers often want to assess how robust the results are to the effect of unmeasured variables. Sensitivity analysis is a method for judging how sensitive a model is to confounding. One can even argue that a statistical analysis is not complete without sensitivity analysis (Imai, Keele, and Tingley 2010). Issues related to confounding have been a source of concern for researchers across disciplines with some of the more acrimonious arguments in the biomedical literature. For example, Cornfield et al. (1959/2009) argued that although the literature at the time indicated that there was a relation between smoking and cancer, researchers should be wary in making any causal claims. The authors also posit that the research evidence at the time suggested that smokers were nine times more likely to develop cancer than nonsmokers and any confounding variables would have to be nine times more prevalent in the smoking than the nonsmoking population. The idea of identifying alternative variables related to the outcome with a magnitude that nullifies the hypothesized or observed relation is central to many sensitivity analyses.

Methods for assessing sensitivity to confounding in regression have been available for some time (cf. Frank 2000; Greenland 1996; Lin, Psaty, and Kronmal 1998; Mauro 1990; Rosenbaum 1986). Some critics have argued that these methods are cumbersome and difficult to interpret, resulting in a lack of use over the years (Frank 2000). Increased interest in mediation has made sensitivity analysis a salient issue because of the assumption of sequential ignorability. In order to make sensitivity analysis for mediation more accessible, we describe three methods for assessing the sensitivity of results to confounder bias and corresponding plots. A goal of all methods discussed here is to obtain information on when a confounding variable effect would be large enough to reduce the mediated effect to zero; however, researchers can choose any criterion that is pertinent to their research interests. The three methods were chosen based on their ease of use and their ability to display the sensitivity information graphically.

The first approach described by VanderWeele (2010) is based on the relation of the confounder to Y and the difference in proportion of persons with the confounder prevalence between treatment groups at the same level of the mediator. The second method presents confounder bias as correlated error terms between the error in the mediator model and the error in the outcome model (Imai, Keele, and Yamamoto 2010). The Imai, Keele, and Yamamoto (2010) and VanderWeele methods use counterfactual definitions of mediated effects as described by Robins and Greenland (1992) and Pearl (2001, 2012). A third method is based on the correlations of a potential confounder with study variables and is adapted from Mauro (1990).

Binary confounder method

VanderWeele (2010) tests the sensitivity of the mediated effect to the violation of the assumption that there is no unmeasured confounding affecting the relation between the mediator and the outcome. The formulas for the bias due to the confounder are based on the expected potential outcome differences. The confounder bias plot displays the value of NIE as a function of two parameters: γ and δ. The γ coefficient corresponds to the effect of an unobserved binary confounder on Y for individuals with the same value of M. The δ coefficient is more complicated but describes the problem with typical mediation analysis. The δ coefficient is the difference in prevalence of the confounder variable, U, for persons with the same value of M in the treatment and control groups. The bias in NIE due to an unmeasured confounder variable (U), conditional on measured covariates (C) is formulated as (see VanderWeele 2010 for computational details):

B i a s N I E (x) = - δ γ .

The extent to which a confounder affects the indirect effect estimate can be obtained by subtracting the value γδ from the potentially biased estimate of the indirect effect. It should be noted that the computation of the bias as shown above makes the following assumptions: first, there is no correlation between U and X when conditional on C. Second, there is no interaction of U with the variables X and C. Third, there is no relation between U and C when conditioned on X and M.

VanderWeele’s formulas to compute the confounder bias in mediation analysis are general so that they can be applied in case of various estimation methods such as marginal structural models, structural mean models, and clustered data settings. Additionally, this method only requires the researcher to input the mediated effect and not the raw data.

Correlated residuals method

Imai, Keele, and Yamamoto (2010) proposed a sensitivity analysis method applicable when using coefficients from the regression approach to estimate the mediated effect. Their method measures the sensitivity of the mediated effect to pretreatment unmeasured confounders. The method assesses how robust the results are to the departure from the sequential ignorability assumption.

If the sequential ignorability assumption holds, then residuals in Equations 2 and 3 (i.e., e ₂ and e ₃) should be uncorrelated. Thus, the sensitivity parameter in this method is the correlation (ρ) between e ₂ and e ₃. When the assumption of sequential ignorability is met, ρ is equal to zero; ρ being different from zero indicates departure from sequential ignorability. The results of the sensitivity analysis show how much the indirect effect changes as a function of ρ. In other words, the NIE, called δ(x) in Imai, Keele, and Yamamoto’s (2010) notation, is plotted as a function of the correlation ρ between e_i ₂ and e_i ₃. The following equation from Imai, Keele, and Yamamoto (2010) states that the average causal indirect effect given as ρ is identified and given by:

\overset{ˉ}{δ} (0) = \overset{ˉ}{δ} (1) = \frac{(a σ_{1})}{σ_{3}} {\tilde{ρ} - ρ \sqrt{(1 - {\tilde{ρ}}^{2}) / (1 - ρ^{2})}},

where σ _j ² $\equiv$ var (e_ij ) for j = 2,3 and $\tilde{ρ} \equiv C o r r (e_{i 1}, e_{i 3})$

When the sequential ignorability assumption may not hold, the best possible bounds for the sensitivity parameter ρ and the mediated effect are (−1, 1) and (−∞, ∞) respectively. The R software package mediation (Imai et al. 2010) can accommodate various models such as linear regression models, generalized linear models, ordered response models, generalized additive models, and quantile regression models. The program first runs the causal mediation analysis, then conducts the sensitivity analysis, and provides the plots to interpret the sensitivity analysis results. The plots illustrate the estimated average mediated effect as ρ increases, though there is no cutoff value for the bias that would indicate a problematic level of susceptibility of the results to confounding. Additionally, Version 7 of Mplus (Muthén and Muthén 2012) can produce the plots from Imai, Keele, Tingley, and Yamamoto (2010).

Left out variables error (L.O.V.E.)

Another way of assessing potential bias in mediation analysis is the L.O.V.E. method (Mauro 1990). This method classifies the error in a regression model that is due to omitted variables as a “misspecification error” and uses correlational techniques to determine the bias. The technique uses the correlations among the predictor (or the mediator), dependent, and omitted variables to detect the potential bias.

If the omitted variable (U) is added to the single mediator model with X as the treatment, M as the mediator, and Y as the outcome variable, the equation with standardized variables becomes

Y = b^{*} M + {c^{'}}^{*} X + g^{*} U + e,

where g is the effect of the omitted variable on Y. The regression coefficient b* is equal to

b^{*} = \frac{r_{y m} (1 - r_{x u}^{2}) + r_{y x} (r_{m u} r_{x u} - r_{m x}) + r_{y u} (r_{m x} r_{x u} - r_{m u})}{1 + 2 r_{m x} r_{x u} r_{m u} - r_{m u}^{2} - r_{x u}^{2} - r_{m x}^{2})} .

The estimates of the mediator M (or the treatment X) effects on the dependent variable are biased if there is a correlation between the omitted variable U and other predictors in the model. However, we usually do not know the values of these correlations. But the limits of these correlations are known using the mathematically possible values. Consider correlations between three variables X_i , X_j , and X_k . The correlation between X_i and X_j can only be in the following range:

r_{i k} r_{j k} \pm \sqrt{(1 - r_{i k}^{2}) (1 - r_{j k}^{2})} .

Using the mathematically possible range of the correlations among the omitted variable and other variables in the model, one can obtain a range of values for b*. When that range is known, the question of interest is how large does the difference between b* in the model without the omitted variable and b* in the model including the omitted variable need to be in order to create substantial bias in the estimate of the mediated effect. For each of the methods outlined here, a substantially different estimate of the coefficient is one that makes a truly nonsignificant effect appear statistically significant and vice versa. To numerically determine substantial difference, one can use the increment in the proportion of variance explained by the mediator and test it against the residual variance in an F-ratio (k = the number of predictors in the model including the omitted variable; N = sample size; R ² = proportion of variance in the dependent variable Y explained by the model):,

F (1, N - k - 1) = \frac{p r o p . (v a r)}{(1 - R^{2}) / (N - k - 1)}, w h e r e

p r o p . (v a r) = \frac{B_{p}^{2}}{(1 - r_{x u}^{2}) / (1 + 2 r_{m x} r_{x u} r_{m u} - r_{m u}^{2} r_{x u}^{2} r_{m x}^{2})} .

This test produces a quadratic equation with two roots that are the limits of the range of the correlation between the omitted variable and the dependent variable that makes the effect of the mediator on the dependent variable nonsignificant. Using this method, researchers can conduct an omitted variable analysis that produces an output for the possible values of the correlations between the omitted variable and the mediator (r_mu ), treatment (r_xu ), and outcome variables (r_yu ) that may lead to a substantial bias in the estimator for the mediator (or treatment) effect arising from not including the omitted variable.

Plots of confounder bias

Three methods of investigating confounder bias with plots were used. A way to investigate the sensitivity of model results to confounding variables is to assess how large a confounder effect on the M to Y relation must be in order to invalidate conclusions about mediation. In other words, the plots portray the minimum strength of the relation of a potential confounder U and both M and Y necessary in order to make the true b path, and consequently the true mediated effect ab, zero.

The three methods described above use different sensitivity parameters to describe how a confounder would affect mediation results. Note that these methods require untestable assumptions of the use of these confounder bias assessment techniques such as that the unmeasured confounder does not interact with measured confounders (VanderWeele 2010).

SAS and R programs were written to make the plots from VanderWeele (2010), which only require the value of the indirect effect as input. SAS and R programs were written to create the L.O.V.E. plot using the three correlations among X, M, and Y as input. The R package mediation, written by Imai, Keele, Tingley, and Yamamoto (2010), was used to make the plot as a function of the correlation between error terms. The interpretation of the VanderWeele and the L.O.V.E. plots is similar to that of the Imai et al. plot of ρ; that is, given the input values for the respective methods, the values necessary to reduce the mediated effect to zero were generated and plotted. The ease of implementing each method varies. The Imai, Keele, Tingley, and Yamamoto method requires the raw data, the L.O.V.E method requires the correlations among X, M, and Y, and the VanderWeele method only requires the mediated effect (ab). In order to illustrate how the graphs are used and interpreted, real data were used to demonstrate when confounding may and may not be a potential concern. R (R Core Development Team 2012) was used to create the graphs in the following examples. Appendix contains the code and directions necessary for researchers to use their data to create similar graphs. Table 1 gives a brief overview of each of the methods in terms of necessary software, input parameters, and under what conditions each method would be useful.

Table 1.

Description of Application Methods for Three Sensitivity Analysis Methods.

Sensitivity Analysis	Software	Necessary Input	Conditions When Can be Used
Correlated residuals (Imai et al. 2010)	R	Raw data specifying outcome variable, mediator, treatment, and covariates	Continuous outcome variable; continuous mediator, binary outcome, binary mediator, multiple mediators, interaction between X and M
VanderWeele (2010)	R SAS	Mediated effect (e.g., the product of a × b)	Any model that produces the mediated effect ab
L.O.V.E (Mauro 1990)	R SAS	Correlations among X, M, and Y including r_xm , r_xy , r_my	Single mediator models with continuous dependent and mediator variables

In order to show how the plots look under circumstances where the population parameters are known, two sets of data were generated to show situations where confounder bias had a large effect and a small effect on a single mediator model. For each situation, a sample of 10,000 data points were generated randomly from a normal distribution with the following standardized regression parameters: for the large effect of confounder bias a = .59, b = 0, c′ = 0, f = .59, and g = .59; for the small effect of confounder bias a = .59, b = 0, c′ = 0, f = .14, and g = .14, where f and g represent the effect on the confounder on the mediator and the outcome, respectively. The values for small (.14) and large (.59) are based on Cohen’s (1988) conventions for small and large effects for multiple regression corresponding to 2% and 26% of the variance explained in the outcome variable. The b path was set to zero to demonstrate a situation when there is no mediated effect through the hypothesized variables, resulting in the confounded accounting for the mediated effect. All data were generated in R (R Core Development Team 2012). Figure 2 shows what each of the three plots looks like for the two examples.

Figure 2.

Sensitivity plots for generated data.

The plots in the left column of Figure 2 represent the situation where there is a small effect of confounding and the plots in the right column represent the situation where there is a large effect of confounding. In the Imai, Keele, Tingley, and Yamamoto (2010) plots, the monotonic curve represents the different values of the average causal mediated effect at different values of ρ. The dotted line represents the value of the average causal mediated effect when ρ is zero. The further the dotted line is from the x-axis and the further the monotonic curve is from the y-axis when it crosses the x-axis, the less likely confounder bias is a problem in the model.

For the L.O.V.E plots, the x-axis represents the correlations between an unmeasured confounder and the outcome variable (Y) and the y-axis represents the correlations between an unmeasured confounder and the mediator variable (M). The curved line represents what values of r _um and r _uy (denoted RUM and RUY, respectively, in the graph) are necessary to make the mediated effect become zero. The further the line is from the axes, the less likely that a confounder is influencing the model.

For the VanderWeele (2010) plots, the x-axis represents values of delta (δ) and the y-axis represents values of gamma (γ). The curved line represents the values of δ and γ that would cause the mediated effect to become zero. The further this line is from the axes, the less likely that a confounder is influencing the model.

Examples

Below is a description and interpretation for each of the sensitivity plots using data from the Athletes Training and Learning to Avoid Steroids (ATLAS) intervention (Goldberg et al. 1996). ATLAS was an intervention designed to prevent steroid use among high school athletes. Between pretest and posttest, the experimental schools received a 7-week, 14-session prevention program, while the control schools received only a pamphlet about steroid use. For clarity of the presentation, only complete cases from baseline to first follow-up and baseline to 1 year follow-up were analyzed in this presentation (Figure 3).

Figure 3.

Sensitivity plots for the Athletes Training and Learning to Avoid Steroids data.

The participants were 1,506 adolescent football players from 31 high schools in Oregon and Washington. The average age at baseline was 15 years and 10 months and 78.3% were White, 5.3% were African American, 3.3% were Asian, 2.7% were Hispanic, 0.8% were Native American, and 9.6% had a mixed ethnic heritage. Of the initial 1,506 participants, 81.4% (n = 1,226) were present at posttest. From pretest to 1-year follow-up, 57.7% (n = 869) of the sample was retained.

Thirty-four schools were randomized to conditions after matching on win/loss records and socioeconomic status. Three schools in the experimental group dropped out prior to the start of the study. The hierarchical structure was not included in the statistical analysis for this article. The players were measured in late August 1994, when the football season started (pretest), in November 1994 near the end of the football season (posttest), and again in late August 1995 (1-year follow-up). Graduating seniors were assessed just prior to graduation in spring 1995 as they would not be available for assessment in August. Their data were combined with the August 1995 data to comprise the 1-year follow-up assessment (N Program = 94, N Control = 146).

The measures were developed based on previous research on androgenic and anabolic steroids with adolescents. More on the ATLAS project can be found in Goldberg et al. (1996) and more on ATLAS mediation analysis can be found in Mackinnon et al. (2001). The mediating variable for all analyses was team as an information source (TIS) measured at posttest. The dependent variable for the first example was strength training self-efficacy (STSE) at posttest and the dependent variable for the second example was nutrition behavior (NB) at 1-year follow-up.

Example 1

L.O.V.E. plot

The range of correlations between the variable TIS at posttest (M) and other observed covariates in the data set is from −.13 to .52, and the correlations between STSE at posttest (Y) and observed covariates in the data set range from −.11 to .53 (this information is not displayed on the graph; it was computed from the observed data). The line in the L.O.V.E plot displays when the value of ab is close to zero. If we look at a single point on the line of ab = 0 on the graph, for example (r _um = .5, r _uy = .6), this means that the correlation between TIS and some potentially unmeasured confounder would have to be .5 and the correlation between STSE and some unmeasured confounder would have to be .6 for the true mediated effect ab to equal zero, indicating confounding. Although TIS and STSE can attain a correlation as high as .5 with observed covariates and thus potentially with some unmeasured variable as well, it is highly unlikely that substantive researchers who are experts on the topics being studied would overlook such a possible influential variables. On the other hand, if the correlations of TIS source with an unmeasured variable, and STSE with the same unmeasured variable necessary for true ab to equal zero were smaller, it would be easier to believe that researchers could have overlooked such a variable and thus the possibility of confounding of the mediated effect would be much more plausible.

VanderWeele plot

For the plot using VandeWeele’s (2010) method, delta (δ) represents the difference in prevalence of the binary confounder U in the two levels of binary X (exposed/treated and unexposed/untreated) for a given level of M, and gamma (γ) represents the difference in expected outcome, STSE, when comparing the presence of U (U = 1) and the absence of U (U = 0) conditional on a given value of the mediator, TIS. The smaller the values of δ and γ for which ab = 0, the more susceptible the mediated effect is to confounding by an unmeasured binary variable. For example, this plot shows that the difference in prevalence between the groups of exposed and unexposed participants would have to be 42% and the difference in the expected value of STSE between U = 1 and U = 0 would have to be 0.4 in order for the observed value of the mediated effect to be due to confounding.

Correlated residuals method

The results of Imai, Keele, and Yamamoto (2010) R package to assess sensitivity of the mediated effect states that ρ would have to equal .35 for the true mediated effect to be zero. Although there is no known scale for evaluating what qualifies as a small, medium, and large correlation between residuals, the above value of ρ was the highest among all ρ values for all mediated effects in the ATLAS data. Thus, by comparison, the mediated effect of exposed/unexposed group membership (X) on STSE (Y) through TIS (M) requires the highest correlation between residuals in order to invalidate the existence of a mediated effect, which means that this model is one where the observed mediated effect is less likely to be due to an unmeasured confounder.

Example 2

L.O.V.E. plot

As before, the range of correlations between the variable TIS at posttest (M) and other observed covariates in the data set is from −.13 to .52. The correlations between STSE at 1-year follow-up (Y) and observed covariates in the data set range from −.14 to .56 (this information is not displayed on the graph; it was computed from the observed data). Following the same logic for interpreting L.O.V.E. plots as before, one can pick a point on the line ab = 0, for example (r _um = .27, r _uy = .4), and argue that a potential confounder producing correlations well within the range of existing correlations in the observed data (.01–.56) with TIS (M) and NB (Y) was unmeasured and could account for the observed mediated effect.

Binary confounder plot

In this example, the δ and γ values are lower than the first example, meaning that the mediated effect is more susceptible to confounding than that in the previous example. The difference in prevalence between the groups of exposed and unexposed participants would have to be 20% (half of the difference in the previous example) and the difference in the expected value of NB (Y) between U = 1 and U = 0 would have to be .3 in order for the observed value of the mediated effect to be due to confounding. There are no set values for what constitutes realistic and concerning values of δ and γ; however, substantive researchers might have an intuition necessary to make such judgments in their respective areas. Failing that, researchers could review the relevant literature for potential confounding variables.

Correlated residuals plot

For the relation between exposed/unexposed group membership (X) and NB at 1-year follow-up (Y) through TIS at posttest (M), ρ is .1. Comparing the same model where the dependent variable was STSE and ρ was .35, it can be seen that the former observed mediated effect is at a higher risk of being due to the presence of an unmeasured confounder than the latter. Again, although there is no known metric that specifies a value of ρ that is considered small, by comparison one can say that the mediated effect from group to NB is more susceptible to potential confounders than the mediated effect from exposed/unexposed group membership to STSE. Even without a comparator, because the values of ρ bounded from −1 to 1, values further from zero indicate confounding.

Discussion

Comparative interpretations of the parameters across the three methods are difficult because of the different meanings of each of the parameters. However, graphing the parameters provides a method of interpretation that, although subjective, allows for the assessment of confounding across the different approaches. For all three approaches, the closer the curve is to the x-axis, the more likely confounding occurred in the mediation process.

Theory is an integral component of interpreting both these plots and the corresponding parameters when the plots suggest that confounding is likely to have occurred. Specifically, the researcher must determine likely confounders of the results as well as the magnitude of the confounding effect. In Example 2, a likely confounder could have been acquiescence bias. Respondents may have been responding in such a way that caused them to provide more socially desirable responses for both M and Y.

These methods can easily be applied post hoc to mediation analyses. Ideally, the researcher would use all of the methods to determine the likelihood of confounding and report these findings in the results section. By providing this information, the authors will either provide evidence for the validity of their findings or information about likely confounders that need to be addressed in future research. In either case, this approach to sensitivity analysis of mediation models serves as a much-needed assessment for determining if the sequential ignorability assumption holds.

Although the most obvious method for dealing with confounders is to include all relevant variables (including confounders) in the mediation model, there are several other approaches to dealing with confounding to help provide more accurate tests. These methods include using propensity scores (Coffman 2011; Jo et al. 2011), instrumental variables (Angrist, Imbens, and Rubin 1996; Jo 2008), principal stratification (Frangakis and Rubin 2002; Rubin 2004; Gallop et al. 2009; Jo 2008), structural nested mean models (Robins 1994; Ten Have et al. 2007), and the “phantom” variable approach (Finkel 1995) used in structural equation models. A thorough explanation of these methods is beyond the scope of this article, but these methods can be used in conjunction with the sensitivity analyses described previously to verify and support the assumption of sequential ignorability.

Experimental design is another way to investigate causal mediation. Imai, Tingley, and Yamamoto (2013) have outlined several experimental design methods for evaluating causality in mediation models and suggested alternative encouragement designs that address the problem of manipulating the mediator in mediation analysis. These experimental methods are also limited as noted by comments on the article. These limitations include several strong and potentially untenable assumptions made in order to determine the causal effect. Another criticism of these approaches is they involve randomizing individuals to the mediator, which is often not possible, particularly in evaluation research. Sensitivity analyses described in this article are important for the interpretation of the results from these experimental studies.

Limitations

Although these sensitivity analysis methods provide a critical first step in determining the plausibility of the assumption of sequential ignorability, there are some limitations. First, there is not yet a clear metric that can be used to compare the different parameters in each of the methods. The ρ value is not clearly comparable to δγ nor the correlations between a confounder and M and Y. Also, VanderWeele’s method assumes confounders to be dichotomous and the δγ parameters are interpreted as such. Determining how these parameters relate to each other will help to provide a more objective measure of when confounding exists and will help make comparisons across methods more interpretable. Second, the examples provided here used only single mediator models with continuous outcomes and no covariates. The VanderWeele (2010) and Imai, Keele, and Yamamoto (2010) methods can incorporate both covariates and dichotomous outcomes, while the L.O.V.E. plots would require additional programming to conduct these analyses such as using correlations partialed for additional variables.

Last, the examples provided here show consensus among the three plots for each of the two examples, but this may not always be the case. Because each of the methods use different parameters to assess confounding, there is likely to be circumstances where one method suggests that a result is sensitive to confounding and another method may not. Future research should focus on understanding when and why this may occur.

The debate regarding confounder bias can be acrimonious because there is no one best method for assessing sensitivity or controlling for confounding. These issues become more contentious as the mediation models become more complex with the inclusion of binary mediator and outcome variables, multiple mediators, and interactions between X and M (see Imai, Keele, and Yamamoto 2010; VanderWeele 2010). There is some progress in this area. For example, Imai and Yamamoto (2013) have recently demonstrated that if multiple parallel mediators can be shown to be independent after conditioning on the treatment and covariates, the average causal mediated effect can still be estimated via a method described by Imai, Keele, and Yamamoto (2010). Over time, as our understanding of these methods in applied settings increases, their use will also increase, ultimately leading to a consensus regarding how to address confounder bias.

Footnotes

Appendix

sensitivity.bias = function(data.sens, delta.uv, gamma.uv, x, m, y, name.m, name.x){

#The following function can be used to create plots assessing sensitivity to confounder bias.

#Each of the values that need to be inputted are listed below. Values need to be inputted in this order for the plots to work correctly.

#Each plot is saved to the working directory as a jpeg file. The file names are “Love plot.jpg”, “VanderWeele plot.jpg”,

#and “Imai plot.jpg”

#data.sens = data to be used for the Imai plots

#delta.uv = This is the upper value of delta for the data set and VanderWeele plot

#gamma.uv = This is the upper value of gamma for the example data set and VanderWeele plot

#”name.m” = Variable name for the mediator in the dataset listed in data.sens. Variable name must occur between quotation marks.

#”name.x” = Variable name for the treatment in the dataset listed in data.sens. Variable name must occur between quotation marks.

model.m=lm(m∼x, na.action=na.omit, data=data.sens)

sum.model.m=summary(model.m)

a=sum.model.m$coefficients[2,1]

model.y=lm(y∼x+m, na.action=na.omit, data=data.sens)

sum.model.y=summary(model.y)

b=sum.model.y$coefficients[3,1]

med=a*b

RYX=cor(y, x)

RMX=cor(m, x)

RYM=cor(y, m)

#**********Imai et al. Plot************

attach(data.sens)

name.x=as.character(name.x)

name.m=as.character(name.m)

install.packages(“mediation”)

library(mediation)

Imai=mediate(model.m, model.y, sims=2000, treat=name.x, mediator=name.m)

Imai.sens=medsens(Imai, rho.by=.01, eps=.01)

jpeg(“Imai plot.jpg”)

plot(Imai.sens, ylim=c(-1,1))

dev.off()

detach(data.sens)

#*********VanderWeele Plot**************

# This program plots VanderWeele's (2010) sensitivity analysis,

#VanderWeele, T. J. (2010). Bias formulas for sensitivity analysis for direct

#and indirect effects. Epidemiology, 21, 540-551.

comb.g=gamma.uv*1000+1

comb.d=delta.uv*1000+1

comb.t=comb.g*comb.d

t=0

gamma.vw=rep(0, comb.t)

delta.vw=rep(0, comb.t)

ab.vw=rep(0, comb.t)

gamma=seq(0, gamma.uv, by=.001)

cases.gamma=comb.g

delta=seq(0, delta.uv, by=.001)

cases.delta=comb.d

for(i in 1: cases.gamma){

for(j in 1: cases.delta){

t=t+1

gamma.vw[t]=gamma[i]

delta.vw[t]=delta[j]

ab.vw[t]=med

}

abnew=ab.vw-(gamma.vw*delta.vw)# The value abnew represents the mediated effect adjusted for potential values of bias.

rabnew=round(abnew, digits=3)# The user can choose the precision of the estimate by changing the value of digits=3.

VanderWeele=data.frame(ab.vw, gamma.vw, delta.vw, abnew, rabnew)

new=VanderWeele[which(VanderWeele$rabnew==0),]

install.packages(“Hmisc”)

library(Hmisc)

jpeg(“VanderWeele plot.jpg”)

#The following lines produce a plot of combinations of values of delta and gamma for which the true mediated effect is zero.

plot(gamma.vw∼delta.vw, data=new, col=2, pch=1, xlab=“delta”, ylab=“gamma”, ylim=c(0, gamma.uv), xlim=c(0, delta.uv), lab=c(10,10,10), tcl=-.5)

minor.tick(nx=10, ny=10, tick.ratio=.4)

dev.off()

#************LOVE PLOTS********************

# The following code plots the combinations of values of correlations of the confounder U with M and Y for which the true mediated effect is zero.

# It was adapted from Mauro (1990).

#The code below creates the values for the correlations between u/y, u/m, and u/x and puts them in a data frame; the user need not change any values.

t=0

ruy.love=rep(0,10011001)

rum.love=rep(0,10011001)

rux.love=rep(0,10011001)

RUM=seq(0,1, by=.001)

cases.rum=1001

RUY=seq(0,1, by=.0001)

cases.ruy=10001

for(i in 1: cases.rum){

for(j in 1: cases.ruy){

t=t+1

if (med < 0){

gamma.vw[t]=-gamma[i]

} else {

gamma.vw[t]=gamma[i]

}

delta.vw[t]=delta[j]

ab.vw[t]=med # This is the observed value of the mediated effect; it can be changed by the user.

}

}love.cor=data.frame(rum.love, ruy.love, rux.love)

RUM=love.cor$rum.love

RUY=love.cor$ruy.love

RUX=love.cor$rux.love

# Values of a, b, cpr that account for the influence of a potential confounder U.

CPR=(RYX*(1-RUM^2)+RYM*(RUX*RUM-RMX)+RUY*(RUX*RUM- RUX))/(1+2*(RMX*RUM*RUX)-RUX^2-RUM^2-RMX^2)

B=(RYM*(1-RUX^2)+RYX*(RUM*RUX-RMX)+RUY*(RMX*RUX- RUM))/(1+2*(RMX*RUX*RUM)-RUM^2-RUX^2-RMX^2)

A=(RMX-RUM*RUX)/(1-RUX^2)

#Observed standardized values of a, b, cpr.

CPRBIASED=(RYX-RYM*RMX)/(1-RMX^2)

BBIASED=(RYM-RYX*RMX)/(1-RMX^2)

ABIASED=RMX

BIASCPR=CPRBIASED-CPR

BIASB=BBIASED-B

BIASA=ABIASED-A

#The following code computes standardized ab that accounts for potential bias, standardized biased (observed) ab, and the potential bias of the mediated effect.

TRUEAB=A*B

BIASEDAB=ABIASED*BBIASED

BIASAB=BIASEDAB-TRUEAB

RTRUEAB= round(TRUEAB, digits =2)# This line of code rounds the estimate of ab that accounts for bias; the user can change the value of digits=2.

love.cor=data.frame(RUM, RUY, RTRUEAB, BIASAB)

#This line of code gets rid of the unnecessary variables which speeds up the processing.

rm(t, ruy.love, rum.love, rux.love, RUM, RUY, RYX, RMX, RYM, CPR, A, B, CPRBIASED, BBIASED, ABIASED, BIASCPR, BIASA, BIASB, TRUEAB, BIASEDAB, BIASAB, RTRUEAB, RUX)

#The code below creates a data frame with combinations of correlations where AB is equal to 0.

for.zero=love.cor[which(love.cor$RTRUEAB==0),]

rm(love.cor, cases.rum, cases.ruy)

#Below is the code to create the plot.

# The name of the axes can be changed by changing the name in the quotation marks for xlab and ylab;

#the former should refer to the correlation of the potential confounder and the dependent variable,

#and the latter should refer to the correlation of the potential confounder with the mediator.

jpeg(“LOVE plot.jpg”)

plot(RUM∼RUY, data=for.zero, xlab=“RUY”, ylab=“RUM”, col=2, pch=1, xlim=c(0,1), ylim=c(0,1), lab=c(10,10,10), tcl=-.5)

dev.off()

}

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This research was supported in part by the National Institute on Drug Abuse, Grant No. R01DA009757 and by the National Institute of Mental Health Grant No T32 MH18387. Collection of the ATLAS data was supported by National Institute on Drug Abuse, Grant No. R01DA07356. Part of this research was presented at the 2012 Harvard University Frontiers of Causal Inference Conference.

References

Angrist

J. D.

Imbens

G. W.

Rubin

D. B.

. 1996. “Identification of Causal Effects Using Instrumental Variables (with discussion).” Journal of the American Statistical Association 91:444–72.

Baron

R. M.

Kenny

D. A.

. 1986. “The Moderator-mediator Variable Distinction in Social Psychological Research: Conceptual, Strategic, and Statistical Considerations.” Journal of Personality and Social Psychology 51:1173–82.

Coffman

D. L.

2011. “Estimating Causal Effects in Mediation Analysis using Propensity Scores.” Structural Equation Modeling: A Multidisciplinary Journal 18:357–69.

Coffman

D. L.

Zhong

. 2012. “Assessing Mediation Using Marginal Structural Models in the Presence of Confounding and Moderation.” Psychological Methods 17:642–64.

Cohen

1988. Statistical Power Analyses for the Behavioral Sciences. 2nd ed. Mahwah, NJ: Erlbaum.

Cornfield

Haenszel

Hammond

E. C.

Lilienfeld

A. M.

Shimkin

M. B.

Wynder

E. L.

. 2009. “Smoking and Lung Cancer: Recent Evidence and a Discussion of Some Questions.” International Journal of Epidemiology 38:1175–91 (Reprinted from Journal of the National Cancer Institute , pp. 173–203, 1959).

Finkel

S. E.

1995. Causal Analysis with Panel Data, vol. 105. Thousands, CA: Sage, Incorporated.

Frangakis

C. E.

Rubin

D. B.

. 2002. “Principal Stratification in Causal Inference.” Biometrics 58:21–29.

Frank

K. A.

2000. “Impact of a Confounding Variable on a Regression Coefficient.” Sociological Methods & Research 29:147–94.

10.

Gallop

Small

D. S.

Lin

J. Y.

Elliot

M. R.

Joffe

Ten Have

T. R.

. 2009. “Mediation Analysis with Principal Stratification.” Statistics in Medicine 28:1108–30.

11.

Goldberg

Elliot

D. L.

Clarke

G. N.

MacKinnon

D. P.

Zoref

Moe

Wolf

S. L.

. 1996. “The Adolescents Training and Learning to Avoid Steroids (ATLAS) Prevention Program: Background and Results of a Model Intervention.” Archives of pediatrics & adolescent medicine 150:713.

12.

Greenland

1996. “Basic Methods for Sensitivity Analysis of Biases.” International Journal of Epidemiology 25:1107–16.

13.

Holland

P. W.

1986. “Statistics and Causal Inference.” Journal of the American Statistical Association 81:945–60.

14.

Holland

P. W.

1988. “Causal Inference, Path Analysis, and Recursive Structural Equation Models.” Sociological Methodology 18:449–84.

15.

Hyman

H. H.

1955. Survey Design and Analysis: Principles, Cases, and Procedures. Glencoe, IL: Free Press.

16.

Imai

Keele

Tingley

. 2010. “A General Approach to Causal Mediation Analysis.” Psychological Methods 15:309–34.

17.

Imai

Keele

Tingley

Yamamoto

. 2010. “Advances in Social Science Research Using R. In Causal mediation analysis using R, edited by Vinod

H. D.

, 129–54. New York, NY: Springer.

18.

Imai

Keele

Yamamoto

. 2010. “Identification, Inference, and Sensitivity Analysis for Causal Mediation Effects.” Statistical Science 25:51–71.

19.

Imai

Tingley

Yamamoto

. 2013. “Experimental Designs for Identifying Causal Mechanisms.” Journal of the Royal Statistical Society: Series A (Statistics in Society) 176:5–51.

20.

Imai

Yamamoto

. 2013. “Identification and Sensitivity Analysis for Multiple Causal Mechanisms: Revisiting Evidence from Framing Experiments.” Political Analysis 21:141–71.

21.

James

L. R.

Brett

J. M.

. 1984. “Mediators, Moderators, and Tests for Mediation.” Journal of Applied Psychology 69:307–21.

22.

2008. “Causal Inference in Randomized Experiments with Mediational Processes.” Psychological Methods 13:314–36.

23.

Stuart

MacKinnon

D. P.

Vinokur

. 2011. “The Use of Propensity Scores in Mediation.” Multivariate Behavioral Research 46:425–52.

24.

Judd

C. M.

Kenny

D. A.

. 1981. “Process Analysis Estimating Mediation in Treatment Evaluations.” Evaluation Review 5:602–19.

25.

Lin

D. Y.

Psaty

B. M.

Kronmal

R. A.

. 1998. “Assessing the Sensitivity of Regression Results to Unmeasured Confounders in Observational Studies.” Biometrics 54:948–63.

26.

MacKinnon

D. P.

2008. Introduction to Statistical Mediation Analysis. Mahwah, NJ: Erlbaum.

27.

MacKinnon

D. P.

Dwyer

J. H.

. 1993. “Estimating Mediated Effects in Prevention Studies.” Evaluation Review 17:144–58.

28.

MacKinnon

D. P.

Goldberg

Clarke

G. N.

Elliot

D. L.

Cheong

Lapin

Moe

E. L.

Krull

J. L.

. 2001. “Mediating Mechanisms in a Program to Reduce Intentions to Use Anabolic Steroids and Improve Exercise Self-efficacy and Dietary Behavior.” Prevention Science 2:15–28.

29.

MacKinnon

D. P.

Kisbu-Sakarya

Gottschall

. 2013. “Developments in Mediation Analysis.” In Oxford Handbook of Quantitative Methods, edited by Little

T. D.

, 338–60. New York: Oxford University Press.

30.

MacKinnon

D. P.

Lockwood

C. M.

Hoffman

J. M.

West

S. G.

Sheets

. 2002. “A Comparison of Methods to Test Mediation and Other Intervening Variable Effects.” Psychological Methods 7:83–104.

31.

Mauro

1990. “Understanding L.O.V.E. (left out variables error): A Method for Estimating the Effects of Omitted Variables.” Psychological Bulletin 108:314–29.

32.

Muthén

L. K.

Muthén

B. O.

. 2012. Mplus (Version 7.0) [Computer Software]. Los Angeles, CA: Muthén & Muthén.

33.

Pearl

2001. Direct and Indirect Effects. Proceedings of the 17th Conference on Uncertainty in Artificial Intelligence, 411–20. San Francisco, CA: Morgan Kaufmann.

34.

Pearl

2012. “The Causal Mediation Formula: A Guide to the Assessment of Pathways and Mechanisms.” Prevention Science 13:426–36.

35.

R Core Development Team. 2012. R (Version 2.15.1) [Computer software]. Vienna, Austria: R Foundation for Statistical Computing. http://www.R-project.org.

36.

Robins

J. M.

1994. “Correcting for Non-compliance in Randomized Trials Using Structural Nested Mean Models.” Communications in Statistics - Theory and Methods 23:2379–412.

37.

Robins

J. M.

Greenland

. 1992. “Identifiability and Exchangeability for Direct and Indirect Effects.” Epidemiology 3:143–55.

38.

Rosenbaum

P. R.

1986. “Dropping Out of High School in the United States: An Observational Study.” Journal of Educational and Behavioral Statistics 11:207–24.

39.

Rubin

D. B.

1974. “Estimating Causal Effects of Treatments in Randomized and Nonrandomized Studies.” Journal of Educational Psychology 66:688.

40.

Rubin

D. B.

2004. “Direct and Indirect Causal Effects Via Potential Outcomes.” Scandinavian Journal of Statistics 31:161–70.

41.

Ten Have

T. R.

Joffe

M. M.

Lynch

K. G.

Brown

G. K.

Maisto

S. A.

Beck

A. T.

. 2007. “Causal Mediation Analysis with Rank Preserving Models.” Biometrics 63:926–34.

42.

VanderWeele

T. J.

2010. “Bias Formulas for Sensitivity Analysis for Direct and Indirect Effects.” Epidemiology 21:540–51.