Assessing Omitted Confounder Bias in Multilevel Mediation Models

ABSTRACT To draw valid inference about an indirect effect in a mediation model, there must be no omitted confounders. No omitted confounders means that there are no common causes of hypothesized causal relationships. When the no-omitted-confounder assumption is violated, inference about indirect effects can be severely biased and the results potentially misleading. Despite the increasing attention to address confounder bias in single-level mediation, this topic has received little attention in the growing area of multilevel mediation analysis. A formidable challenge is that the no-omitted-confounder assumption is untestable. To address this challenge, we first analytically examined the biasing effects of potential violations of this critical assumption in a two-level mediation model with random intercepts and slopes, in which all the variables are measured at Level 1. Our analytic results show that omitting a Level 1 confounder can yield misleading results about key quantities of interest, such as Level 1 and Level 2 indirect effects. Second, we proposed a sensitivity analysis technique to assess the extent to which potential violation of the no-omitted-confounder assumption might invalidate or alter the conclusions about the indirect effects observed. We illustrated the methods using an empirical study and provided computer code so that researchers can implement the methods discussed.

Mediation analysis has become popular in identifying and testing causal mechanisms underlying psychological processes (Judd & Kenny, 1981;MacKinnon, 2008). For example, one can test whether the effect of student academic achievement (ACH; antecedent variable, X) on student career aspirations for the future (FUT; outcome variable, Y) is mediated (completely or partially) by student academic self-concept (ASC; mediator, M). Figure 1a depicts linear relationships between the variables in this example, where a is the effect of X on M; b is the effect of M on Y controlling for X; c is the (direct) effect of X on Y controlling for M; and c is the total effect of X on Y in Figure 1b. For this model, the indirect effect is equivalently represented as the product of coefficients, ab, or the difference in coefficients, c − c (Pearl, 2012).
Because many questions in psychology that involve mediation processes are in the context of clustered data (e.g., clients nested within therapists, employees nested within supervisors, observations nested within person), multilevel mediation is an extremely important tool that is only recently becoming widely used because of the advances in methodological explications and software to implement these complicated models (e.g., Lanaj, Johnson, & Barnes, 2014;Sturgeon, Zautra, & Arewasikporn, 2014;Tofighi & Thoemmes, 2014;Wang et al., CONTACT Davood Tofighi dtofighi@psych.gatech.edu School of Psychology J. S. Coon Bldg., Georgia Institute of Technology,  Cherry Street, Atlanta, Georgia -. Color versions of one or more of the figures in the article can be found online at www.tandfonline.com/hmbr. Supplemental data for this article can be accessed at tandfonline.com/hmbr. 2013). As an example, consider the study by Nagengast and Marsh (2012), which examined the indirect effect of student ACH on student FUT through student ASC (self-perception of student academic skills). In this scenario, students (Level 1 units) are clustered (nested) within schools (Level 2 units). Because of the clustering, student data within schools tend to be correlated. Because of the potential lack of independence among student observations, traditional (single-level) techniques such as ordinary least squares (OLS) regression would produce invalid inference about the coefficients (e.g., underestimated standard errors; Raudenbush & Bryk, 2002;Snijders & Bosker, 2012). Instead, researchers recommend multilevel mediation analysis because it (a) correctly estimates the standard errors of the coefficients, yielding more accurate inference of the indirect effects (Kenny, Korchmaros, & Bolger, 2003;Krull & MacKinnon, 1999, 2001, (b) can estimate indirect effects separately at the student level (Level 1) and school level (Level 2; Preacher, 2011;Tofighi & Thoemmes, 2014), (c) can test whether the Level 2 indirect effect differs from the Level 1 indirect effect (Marsh, 1987;Marsh, Kuyper, Morin, Parker, & Seaton, 2014;Marsh, Trautwein, Lüdtke, Köller, & Baumert, 2005;Pituch & Stapleton, 2012), and (d) can test whether a Level 2 variable mod- erates the Level 1 indirect effects (Bauer, Preacher, & Gil, 2006). Despite the increasing popularity of multilevel mediation analysis, there remain key unresolved issues in the specification and interpretation of important quantities. A salient issue that has received little attention is a disagreement about the proper way to estimate and interpret an indirect effect. Unlike the single-level mediation model, the equality between ab or c − c does not hold in multilevel mediation analysis under certain conditions (Tofighi, West, & MacKinnon, 2013). Some authors argue that c − c is the proper estimate (Bauer et al., 2006;Kenny et al., 2003;Muthén & Muthén, 1998-2012Preacher, 2011); other authors recommend ab as the correct estimate (Tofighi & Thoemmes, 2014;Zhang, Zyphur, & Preacher, 2009).
We argue that such a discrepancy in defining an indirect effect is because of the lack of development in the area of specification assumptions necessary to define and estimate a causal unbiased estimate of an indirect effect in a multilevel mediation model. Whereas the specification assumptions have garnered attention in singlelevel mediation analysis (e.g., Imai, Keele, & Tingley, 2010;Judd & Kenny, 1981;Pearl, 2001;Valeri & Vander-Weele, 2013;VanderWeele, 2010), extending the specification assumptions to the multilevel mediation models has received little attention (Tofighi et al., 2013). One important part of the specification assumptions is the noomitted-confounder assumption (Valeri & VanderWeele, 2013;VanderWeele, 2010). The no-omitted-confounder assumption states that there must be no omitted common causes (i.e., confounders) of the observed variables in a mediation model. The no-omitted-confounder assumption is critical to estimate an unbiased estimate of the indirect effect. Because the no-omitted-confounder assumption is untestable, the indirect effect is potentially biased (Holland, 1988). That is, we cannot rule out the biasing effect of the potential confounders on the indirect effects. It is unfortunate that ramifications of potential violation of the no-omitted-confounder assumption have received little attention in the multilevel mediation literature.
In this article, we explicate the specification assumptions necessary to define and estimate an unbiased causal indirect effect in a 1 → 1 → 1 multilevel mediation model. In the notation 1 → 1 → 1, the first, second, and third number indicate that the antecedent, mediator, and outcome variable, respectively, are measured at Level 1 (Krull & MacKinnon, 2001). 1 We show that a violation of the specification assumptions, especially the no-omittedconfounder assumption, is a source of major disagreement in defining and estimating an indirect effect. That is, when the no-omitted-confounder assumption is violated, the two methods of defining and estimating indirect effect are no longer equal. In addition, estimates of coefficients would be biased, including the estimates of indirect effects, because the indirect effects are functions of the estimates of coefficients. Because the no-omittedconfounder assumption is untestable, we propose a sensitivity analysis technique to assess the extent to which any potential violation of the no-omitted-confounder assumption would change the conclusions about the indirect effects.

Unresolved issues in multilevel mediation analysis
The implications of violating the specification assumptions, including the no-omitted-confounder assumption, have received little attention in multilevel mediation analysis. An exception is the work by Tofighi et al. (2013), which examined the biasing effects of omitting a Level 2 (school-level), not a Level 1 (student-level), confounder in a 1 → 1 → 1 model. Tofighi et al. assumed that the omitted Level 2 confounder was correlated with the mediator and outcome variable, not the antecedent variable. They found that the omitted Level 2 confounder would bias the estimates of the indirect effect at Level 1 and bias the variance-covariance estimates at Level 2. Tofighi et al. concluded that the interpretation of indirect effects would become ambiguous in the presence of the omitted confounder.
Several important issues remain unresolved. First, Tofighi et al. (2013) considered the biasing effects of a Level 2 (e.g., school-level), not a Level 1 (e.g., studentlevel) omitted confounder. Second, Tofighi et al. studied a simplified case in which the Level 2 omitted confounder was assumed to influence the mediator and outcome variable but not the antecedent variable. Cases in which the omitted confounder influences all of the observed variables at both levels of analysis were not considered. Therefore, Tofighi et al. 's results are limited to randomized experimental studies in which omitted confounders only exist at Level 2.
We discuss three scenarios in which omitted confounders can arise in practice (Lash, Fox, & Fink, 2009): 1. Observational or experimental mediation studies with a single unmeasured, but known, confounder at Level 1. This scenario might occur when a theory exists about the relevant Level 1 confounder and yet the confounder is unmeasured. This is likely to happen in archival data sets, in which a researcher might not have had control of the choice of variables being measured. 2. Observational or experimental mediation studies with an unknown confounder or confounders. This scenario occurs when a researcher does not have a theory about the nature of omitted variables. For example, theoretical development and substantive theory have not progressed enough to identify all of the possible omitted confounders. 3. Observational or experimental mediation studies with multiple unmeasured confounders, where a theory identifies important unmeasured confounders and their relationships with the observed variables. These confounders are neither measured nor included in the study (e.g., an archival data set). In this article we present a framework that can examine the biasing effect of omitted confounder(s) in Scenarios 1 and 2. We do not address Scenario 3 because it is beyond the scope of this article.
Another unresolved issue is how to assess the extent to which indirect effect estimates might be biased. For this, we develop a method to assess sensitivity of indirect effect estimates to the potential violations of the no-omitted-confounder assumption. Because the noomitted-confounder assumption is not testable (Holland, 1988), epidemiologists have recommended that researchers probe the sensitivity of indirect effects to potential violation of the no-unmeasured-confounder assumption (Blakely, 2002;Hafeman, 2011). Sensitivity analysis helps researchers determine how sensitive the estimates of indirect effects are to the potential violations of the no-omitted-confounder assumption.

Goals
We first study the multilevel biasing effects of an omitted Level 1 confounder on the interpretation and estimation of indirect effects at both Level 1 and Level 2 for a two-level 1 → 1 → 1 mediation model with random intercepts and slopes. A key goal is to assess the biasing effects of an omitted Level 1 confounder that is correlated at both Levels 1 and 2 with all the variables in a multilevel model. We extend the single-level mediation analysis framework (Cox, Kisbu-Sakarya, Miocević, & MacKinnon, 2013) to analytically examine the multilevel biasing effects by decomposing Level 1 confounder(s) into two orthogonal components. We analytically derive expressions quantifying the magnitudes of confounder bias on the Level 1 and Level 2 coefficients and indirect effects as well as the Level 2 variance of the random indirect effect.
The second goal is to assess the extent to which multilevel omitted confounder bias would change conclusions about Level 1 and Level 2 indirect effects, which we do using a sensitivity analysis. The sensitivity analysis offers estimates of the amount of bias in key estimates that help bracket the likely magnitude of the true indirect effect had it been modeled with the confounders included. The sensitivity analysis addresses the following questions: 1. How large is the relationship between the omitted confounders and X, M, and Y that would yield zero (or nonsignificant) estimates of indirect effects at Level 1 or 2? 2. What are the estimates of indirect effects when omitted confounders are (not/moderately/strongly) correlated with X, M, and Y? We will show how the sensitivity analysis can answer these important research questions in the context of the aforementioned empirical example.
In the single-level single-mediator model, there are several techniques to address the potential violation of the no-omitted-confounder assumption (see MacKinnon & Pirlott, 2015; for a detailed discussion, see Cox et al., 2013). We extend the sensitivity analysis method developed for the single-level single-mediator model (Cox et al., 2013) to the 1 → 1 → 1 model to assess potential bias for the two types of the omitted confounder(s) discussed in Scenarios 1 and 2. Cox et al. 's sensitivity analysis is an extension of Mauro's (1990) technique that used the correlation of an omitted variable with the observed variables in OLS regression to assess changes in the conclusions about regression coefficients. The advantage of our method is that it extends the technique based on the correlations of omitted confounder(s) with the observed variables to both observational and experimental studies (less X) in multilevel mediation studies.

Background: Multilevel mediation model
Before proceeding further, we present background and necessary equations to specify a two-level, 1 → 1 → 1 model with random intercepts and slopes. The results presented in our study are general, but for concreteness we use the aforementioned education example throughout the manuscript. Consider again the empirical example by Nagengast and Marsh (2012). In this 1 → 1 → 1 model, the random slopes capture the heterogeneity of Level 1 (Within) effects across schools; the random intercepts model the Level 2 (Between) relationships. Other two-level mediation models such as 2 → 1 → 1 and 2 → 2 → 1 have also been proposed (Krull & MacKinnon, 2001). We focus on the 1 → 1 → 1 model with random intercepts and slopes because it is the most detailed and complex model. It has the greatest number of both fixed (at least six coefficients, three for each level of analysis) and random effects (ten Level 2 covariances, five Level 2 variances, and two Level 1 variances) compared to other similar models.
For our educational example, we derive analytical results using centering within cluster 2 (CWC2) with latent cluster means for the following reasons, some of which are outlined by Marsh et al. (2009); CWC2 centers variables within cluster (school) and adds the cluster (school) means as covariates into the model at Level 2. First, our research question guided the choice of centering strategy (Enders & Tofighi, 2007). In our example, we are interested in whether a student ASC is a positive function of student ACH and a negative function of school-average ACH score. Students with higher ACH scores are expected to have higher ASC. However, academically selective schools with high-ability students might negatively influence student ASC. Such a differential effect in the school-level (Between) and studentlevel (Within) ACH→ASC is a contextual effect known as the big-fish-little-pond effect (BFLPE; Marsh, 1987). In addition, we chose CWC2 to clearly show the effect of an omitted confounder bias on Between, Within, and cross-level effects in the model. Using CWC2 simplified the analytical results because it allowed us to decompose Between and Within effects for the observed variables as well as the confounder(s). We chose using latent cluster means instead of observed cluster means to be consistent with results in Marsh et al. (2009). Marsh et al. recommended using latent mean centering for this problem. In addition, using latent cluster means is likely to produce less attenuated estimate of Between effects (Lüdtke et al., 2008). 2  We provide a more extensive treatment of centering of the predictors in  →  →  models in the supplemental materials.
Equations for the 1 → 1 → 1 model The first set of equations decompose X ij , M ij , and Y ij into the orthogonal Level 2 and Level 1 latent variables using CWC2 with latent cluster means, where i denotes a student and j denotes a school: In Equation (1), η Xj is the Level 2 latent school (cluster) mean on ACH; this is a random intercept that is school specific. The within-school (Within) latent component, η Xij , measures the deviation of each student's ACH score from his or her school's latent mean. This score shows the standing of each student relative to his or her school's level. For the mediator and outcome variables in (2) and (3), respectively, the Between components, η Mj and η Yj , represent latent school means on ASC and FUT, respectively. For M ij and Y ij , the Within components η Mij and η Yij represent student i's deviation score from his or her school's latent mean, respectively. For example, η Mij represents the standing of a particular student's ASC relative to his or her latent school mean on ASC.
A common practice in specifying a multilevel mediation model is to write equations for Level 1 (e.g., student) and Level 2 (e.g., school). Separately, the Level 1 equations for the population relationships at the student level are The Level 2 equations for the relationships at the school level are One can also specify the following Level 1 and 2 equations to estimate the population total effect of ACH on FUT: Equations (4), (5), and (11) describe population relationships between observed variables for each Level 1 unit i (e.g., student). For example, Equation (4) shows that student ACH predicts ASC. Coefficient a j is the latent effect of ACH on ASC for Level 2 unit j (e.g., school). The subscript j indicates that this effect can vary across schools. Because a j is unobserved and varies across schools, it is also called a random effect (coefficient); b j is the random effect of ASC on FUT, while controlling for ACH; c j is the random effect of ACH on FUT, while controlling for ASC; c j is the total random effect of ACH on FUT. The εs denote the Level 1 residuals.
Because we used CWC2, Level 2 Equations (6), (7), and (12) describe the Level 2 (Between) population coefficients. The Between coefficients are denoted by the subscript "B": a B is the Between effect of ACH on ASC; b B is the Between effect of ASC on FUT, while controlling for ACH; c B is the Between effect of ACH on FUT, while controlling for ASC; and c B is the Between total effect of ACH on FUT.
Finally, the Level 2 Equations (8)-(10) and (13) describe the relationships between random coefficients and population-average coefficients (i.e., the expected values of the random coefficients) for each Level 2 unit j. The population-average coefficients are the Level 1 (Within) coefficients denoted by the subscript "W". For example, Equation (8) shows that a j equals the population-average Within coefficient, a W , which is the expected value of all a j s across all the schools, plus a deviation from the population-average coefficient for school j, denoted by u aj . Henceforth, we call the population-average Within coefficients simply Within coefficients. The Within coefficients, b W , c W , and c W , have similar interpretation as the population average for their respective random coefficients. The u terms, which denote deviations from the Within coefficients, are also referred to as Level 2 residuals.

Distributional assumptions
In a 1 → 1 → 1 model, distributional assumptions about the residuals are critical in obtaining unbiased estimates of indirect effects (Bauer et al., 2006;Tofighi et al., 2013). Here, we provide a general description of the distributional assumptions.
The vector of Level 1 (Within) residuals is denoted by ε = (ε Mi j , ε Yi j ) T , where T denotes vector transpose. The vector of residuals is assumed to have a bivariate normal distribution with the means of zero and a 2×2 Within the variance-covariance matrix denoted by W : where the terms σ 2 ε Mi j and σ 2 ε Yi j are variances and σ ε Mi j ,ε Yi j is the covariance between Level 1 residuals across the equations associated with M and Y.
The vector of Level 2 residual terms, u = (u M j , u a j , u Y j , u b j , u c j ) T , contains five residual terms associated with M and Y in Equations (6)-(10), two for the random intercepts and three for random slopes. The vector is assumed to have a multivariate normal distribution with the means of zero and a 5×5 Between variance-covariance matrix, B , as follows: The off-diagonal elements are the covariances between pairs of the residual terms (or, equivalently, random effects). Considering B , we delineate subtle, but important, differences between two types of covariances between residuals in a 1 → 1 → 1 model (Bauer et al., 2006). The 1 → 1 → 1 model is a bivariate model that consists of two Level 1 equations with two dependent variables (i.e., M and Y). The covariances in the 1 → 1 → 1 model contain both within-equation and betweenequation covariances. Interpretation of within-equation covariances in a 1 → 1 → 1 model is similar to that in a conventional univariate multilevel model (Raudenbush & Bryk, 2002). For example, u a j and u Mj both appear in the Level 2 equations for the mediator M. When an appropriate centering strategy is applied, Snijders and Bosker (2012, Chapter 5) recommend estimating and interpreting the within-equation covariances. For example, withinequation covariance σ u a j ,u M j is the covariance between the random intercept and slope predicting the same dependent variable, M. A negative value for this covariance could indicate that the higher the school-average ACH, the smaller the relationship between student ACH and ASC. Between-equation covariances, however, are specific to a multivariate model such as the 1 → 1 → 1 model. In a conventional multilevel model, there are no betweenequation covariances. Unlike within-equation covariance, a nonzero between-equation covariance might not have a substantive interpretation (Tofighi et al., 2013).

Level 1 and Level 2 indirect effects
In multilevel analysis, indirect effects can be defined at each level. First, we discuss the random indirect effect that is the product of two random coefficients: a j b j . Using CWC2, we define the random indirect effect as a schoolspecific, Within indirect effect for school j that varies randomly around the population-average Within indirect effect, a W b W , or, simply, Within indirect effect.
A point of contention is whether the expected value of the random indirect effect equals the (Within) indirect effect. Several authors broadly define indirect effect as the expected value of the random indirect effect (Bauer et al., 2006;Kenny et al., 2003;Preacher, 2015;Preacher, Zyphur, & Zhang, 2010). When CWC2 is applied, these authors define the Within indirect effect as follows: where σ a j ,b j is the covariance between the random effects a j and b j , or equivalently, between the respective Level 2 residuals (Tofighi et al., 2013). However, Tofighi et al. showed that when there is an omitted confounder at Level 2, the expression does not estimate an unbiased estimate of the Within indirect effect. Next, we provide a set of specification assumptions necessary to compute an unbiased estimate of an indirect effect.

Specification assumptions
We extend the specification assumptions from the singlelevel mediation literature in different methodological traditions (e.g., path analysis and counter-factual framework) to the multilevel mediation analysis (James, 1980;James & Brett, 1984;Judd & Kenny, 1981;McDonald, 1997;Pearl, 2001;Valeri & VanderWeele, 2013;Van-derWeele, 2010). Investigating the correct specification assumptions in mediation analysis is more challenging for a 1 → 1 → 1 model than for a single-level model. The challenges result from the complex nature of the multilevel data. For example, in a 1 → 1 → 1 model (a) the relationships might exist at two levels of analysis; (b) the interaction effects might occur at either level or across the levels; or (c) each variable might be measured at Level 1 or 2, resulting in models with distinct interpretations of indirect effects (e.g., e.g., 2 → 1 → 1 model; Krull & MacKinnon, 1999, 2001Pituch & Stapleton, 2012). As discussed, scaling predictors plays an important role in estimation and interpretation of the model parameters (Enders & Tofighi, 2007;Lüdtke et al., 2008;Preacher, 2011;Zhang et al., 2009). Finally, distributional assumptions about the residuals are more complicated in a multilevel model. For the 1 → 1 → 1 model, as shown previously, one needs to make distributional assumptions about seven residuals instead of two residuals in the single-mediator model in Figure 1.

A correctly specified model
For a correctly specified 1 → 1 → 1 model, the following set of specification assumptions is necessary to have an unbiased estimate of indirect effects: 1. Correct functional forms: The functional relationships as well as the causal order of the relationships between X, M, and Y are correctly specified. 2. Validity and reliability: The observed variables X, M, and Y are perfectly valid and perfectly reliable measures of the respective constructs. 3. No method bias: Method bias might occur when two or more variables in a study are measured using the same method of measurement (Campbell & Fiske, 1959;Podsakoff, MacKenzie, Lee, & Podsakoff, 2003;Richardson, Simmering, & Sturman, 2009). For example, if the observed variables in a mediation study are collected using a self-report questionnaire, a researcher might suspect that "at least some of the observed covariation between them may be due to the fact that they share the same method of measurement" (Podsakoff, MacKenzie, & Podsakoff, 2012, p. 540). The no method bias assumption states that there exists no shared method variance biasing the relationships between X, M, and Y at either level. 4. No omitted confounder: There is no omitted confounder of the X, M, and Y relationships at either level of analysis. More formally, this assumption states that there must be no omitted confounders of the relationships between X, M, and Y that are posited in the 1 → 1 → 1 model in Equations (1)-(10). The no-omitted-confounder assumption has also been termed the "no omitted variables" assumption (MacKinnon & Pirlott, 2015), "orthogonality of residuals" assumption (McDonald, 1997), and "sequential ignorability" assumption (Imai et al., 2010).

No interaction: No interaction effect exists
between the observed variables as well as between the observed and omitted variables at each level and across levels. This assumption is also referred to as the "no-interaction, " "linearity, " or "constanteffect" assumption (Pearl, 2012). Finally, as we discussed previously, the distributional assumptions about multilevel models must also be met for proper statistical inference (e.g., confidence interval [CI] and p values; Raudenbush & Bryk, 2002).
Assumption 1 emphasizes the longitudinal nature of a mediation model in that X is measured before M and M is measured before Y (Davis, 1985;MacKinnon, 2008). Cross-sectional mediation studies in which this condition is not met can lead to invalid results (Maxwell & Cole, 2007;Maxwell, Cole, & Mitchell, 2011). When the above assumptions hold, Tofighi et al. (2013) showed that This result indicates that both methods of calculating Level 1 and Level 2 indirect effects are equivalent. It also implies that both the product-of-coefficient and difference-in-coefficients methods will produce an unbiased estimate of indirect effects at the respective levels of analysis. Assumptions 2, 3, and 5, which are not discussed by Tofighi et al. (2013), are also critical in obtaining equivalent unbiased estimates of the indirect effects using either method. When these assumptions hold, an unbiased, causal estimate of the Within indirect effect is In addition, covariance between a j and b j will be zero: σ a j ,b j = 0. That is, when these assumptions hold, we can define and estimate a causal estimate of the Within indirect effect. However, as will be shown later, if Assumptions 3 and 5 are violated, the relationships in Equations (17) and (18) do not hold. More important, Equation (16) does not provide a causal, unbiased estimate of Within indirect effect.

A misspecified model
To investigate the implications of violation of the noomitted-confounder assumption, we consider a "misspecified" 1 → 1 → 1 model in which the omitted confounder(s) of X, M and Y relationships may exist at both levels. We consider the types of omitted confounder(s) in Scenarios 1 and 2: single unmeasured confounder and unknown confounder(s). In terms of the specification assumptions, we assume that Assumptions 4 and 5 are violated, whereas Assumptions 1-3 hold. For this misspecified model, estimating biasing effects of omitted confounder(s) poses the following challenges.

Challenge 
One key challenge that arises when a potential confounder is correlated with the observed variables at both levels of analysis is that the resulting compounded, multilevel biasing effects can be comprised of Level 1, Level 2, and cross-level effects. That is, a confounder may not only bias the relationships at Level 1 and Level 2, but may also serve as a potential moderator of the Level 1 relationships. This result implies that the no-interaction assumption can also be violated. This violation occurs because the Level 1 omitted variable that varies at both Levels 1 and 2 can also moderate the Within coefficients that substantially vary across Level 2 units (e.g., across schools). In this case, the Within indirect effect would be moderated by omitted confounders at Level 2 (e.g., school characteristics). Thus, we can have an omitted cross-level interaction effect. In addition, omitted confounders at the individual level (e.g., student characteristics) can also moderate the Level 1 indirect effect. In this case, we can have an omitted Level 1 interaction effect.
We investigate the important case of the omitted confounder effect on cross-level interaction. 3 A cross-level effect is of substantive importance in various areas of psychological research and related fields (Raudenbush & Bryk, 2002;Snijders & Bosker, 2012). Omitted crosslevel interaction effects are of substantive interest in school settings; for example, educational researchers are often interested in estimating environmental factors (e.g., school-district characteristics) affecting student performance and what school-level characteristics would moderate student-level indirect effects.

Challenge 
A second challenge is finding a mathematical framework that computes and estimates the biasing effects of omitted confounder(s). The mathematical framework needs to consider that potential biasing effects of an omitted confounder are likely to be a weighted composite of unobserved Level 1 and Level 2 correlations (Kreft, de Leeuw, & Aiken, 1995). To estimate the weighted composite effect of the omitted variable, values for three unobserved quantities are needed: (a) Level 2 correlations of an omitted confounder with the observed variables, (b) Level 1 correlations of an omitted confounder with the observed variables, and (c) the weights for the composite effect. Mauro (1990) described techniques to assess the potential biasing effect of an omitted variable in OLS regression. Cox et al. (2013) extended these techniques to assess omitted confounder bias in a single-level single-mediator model. We extend this framework for 1 → 1 → 1 models. We describe three critical stages in the extended framework: First, we introduce the concept of a latent proxy variable to model confounder bias in Scenarios 1 and 2. Second, we discuss the concept of "augmented" model. Third, we describe the necessity of reexpressing relationships between the latent proxy variable and observed variables to derive analytic results.

Latent proxy variable
To address Challenge 1, we introduce a single latent variable, Z ij , that serves as a proxy for omitted confounder(s) that may exist at both Levels 1 and 2. First, we assume that the latent proxy variable, Z ij , is potentially correlated with all of the observed variables (i.e., the noomitted-confounder assumption is violated); the sign and  Deriving analytic results for both omitted cross-level and Level  interaction effects is beyond the scope of this manuscript.
the magnitude of the correlations may vary. Second, we assume that the latent proxy variable linearly influences the observed variables in the model. That means that potential nonlinear relationship (e.g., quadratic relationship) between the omitted and observed variables are not modeled.
The use of Z ij as a proxy for the potential omitted confounder(s) is justified because a single latent variable is sufficient to account for spurious correlations (covariances) between (residuals of) three observed variables (Brewer, Campbell, & Crano, 1970). By spurious we mean extraneous correlation not accounted for by the posited mediation model. When the no-omitted-confounder assumption is violated, between-equation covariances in W in (14) and B in (15) will be nonzero; 4 otherwise, the between-equation covariances will be zero if specification assumptions 1-5 hold. Introducing Z ij accounts for the spurious correlations between observed variables, thereby rendering the between-equation covariances zero. In addition, because the potential omitted confounders are likely to explain heterogeneity of indirect effects across Level 2 units, we investigate whether Z ij moderates the Level 1 indirect effect, violating the nointeraction assumption. It should be noted that Z ij serves as a proxy for both Level 1 and Level 2 confounders.
The latent proxy variable is suitable to address confounder-bias in Scenarios 1 and 2. In Scenario 1, for a single known unmeasured confounder, a researcher is more likely to find plausible values for the correlations of the unmeasured confounder with the observed variables from the literature or experts. A researcher can use this latent proxy variable along with the range of the correlation values to evaluate the potential biasing effect on the conclusions about the indirect effect.
For the unknown confounders in Scenario 2, the latent proxy variable assumes that the effect of all unknown confounders can be represented by a single latent variable. Making this assumption is reasonable because researchers have made similar types of assumptions about unknown causes when specifying disturbances in single-level structural equation models (SSEM). Disturbance is "the set of unspecified causes of the effect variable. Analogous to an error or residual in a prediction equation ... The disturbance is treated as a latent variable" (Kenny, 2011).

Augmented model
Challenge 2 was to find a mathematical framework that models correlational structure at Levels 1 and 2 as well as the composite weights for multilevel effects of Z ij on the model parameters. We begin by treating the correlational  Except, in a rare situation that other sources of spurious between-equation correlations exist (e.g., common method effect) such that the sum of the spurious correlations would become zero.
structure as known. We propose the following mathematical framework. First, we augment the misspecified 1 → 1 → 1 model by adding the latent proxy variable Z ij . We term this new model the "augmented" model. To address the challenge of estimating composite weights, we use the concept of partitioning a Level 1 variable into Level 1 and Level 2 components from the multilevel centering literature (Enders & Tofighi, 2007). That is, we decompose Z ij into two orthogonal components and then study the biasing effects of each component on the level-specific and cross-level effects: where η Zij and η Zj are the Level 1 and Level 2 components of the latent proxy variable, respectively. This approach permits the biasing effects of the unmeasured confounder(s) to be examined at different levels of analysis. The final augmented 1 → 1 → 1 model with the orthogonal components of Z ij is shown in Figure 2. We specify the Level 1 and 2 equations for the augmented model. The Level 1 equations for the augmented model with the random intercepts and slopes are The Level 2 equations are We now are in a position to study the multilevel biasing effects of omitted confounder(s) using the augmented model.

Reexpressing equations for the latent proxy variable
An important part of the analytic derivation was to reexpress the relations between the latent proxy variable and the antecedent and mediator variables, such that the latent proxy variable is the dependent variable and the observed variables are predictors. More specifically, we reexpressed the Within and Between linear relationships of X ij and M ij Figure . The augmented model with random intercepts and slopes. ACH = academic achievement; ASC = academic selfconcept; FUT = career aspirations for the future; Z ij , a latent proxy variable for omitted confounder (e.g., academic interest), is decomposed into the Within (η Zij ) and Between (η Zj ) latent variables. At the Within part of the graph, filled circles lying in the middle of the arrows represent random slopes. Filled circles at the end of arrows represent random intercepts, whose relationships are further depicted in the Between part of the graph. Note that within-equation correlation (e.g., intercept-slope correlation, σ u a j ,u M j ) for M and Y are not shown. The between-equation covariances (e.g., σ u a j ,u bj ) are zero.
to Z ij as follows: where ss and qs are regression coefficients, and Zi j and u Z j are the residual terms. We detail the interpretation of these coefficients in the next section. We use the coefficients in the aforementioned equations in the analytic results. Before proceeding, it is important to compare the relationships of the coefficients q ZM and s ZM in Equations (29) and (30) to the coefficients l MZ and d MZ in Equations (21) and (25). These four equations model the linear correlation between X ij , M ij , and Z ij at Level 1 or 2, and thus are related. For example, in Equation (30), the coefficient s ZM is proportional to the partial correlation between η Zj and η Mj , while controlling for η Xj . In Equation (25), the coefficient d MZ is also proportional to the partial correlation between η Zj and η Mj , while controlling for η Xj . As a result, the coefficients s ZM and d MZ are proportional to one another. Similarly, q ZM and l MZ are proportional to the partial correlation between η Zij and η Mij , while controlling for η Xij . The next step to derive analytic results is to omit the latent proxy variable from the augmented model and then compute the bias in the resulting misspecified model. The coefficients in the resulting misspecified model are denoted using an asterisk and are termed "biased" (unadjusted). The coefficients in the augmented model (e.g., a B ) are termed "bias corrected" (adjusted). The full analytic derivation can be found in the Appendix.

Analytic results
To make the analytic results more concrete, we present a summary of the analytical results and interpret them using our educational example. Because the values of ACH are not randomized, as is often the case with independent variables in psychology and related fields, the estimates of indirect effects at Levels 1 and 2 can be biased because of the unobserved confounders. For our example, there is evidence that student academic interest is correlated with both ACH and ASC (Marsh et al., 2005). For simplicity, we assume that student academic interest is the only unmeasured confounder as in Scenario 1, which is potentially correlated with the observed variables at Levels 1 and 2. Thus, Z ij denotes academic interest; η Zj denotes Between (school-average, or simply school) academic interest; and η Zij denotes Within (student) academic interest. School academic interest, η Zj , is indicative of the general academic interest at the school level. Student academic interest, η Zij , indicates each student's academic interest relative to the school latent mean academic interest.

Biased between coefficients
Biased (unadjusted) Between coefficients a * B , b * B , and c * B are as follows: For the coefficient a * b , which represents the effect of school ACH on school ASC, the biasing effect of the omitted variable equals d MZ d ZX . The size of d MZ depends on the magnitude of the partial correlation between school academic interest (η Zj ) and school ASC, while controlling for school ACH (see Equation [25]); the size of d ZX depends on the correlation between school academic interest and school ACH (see Equation [23]). As the absolute values of the correlation coefficients increase, the amount of bias will also increase. If the correlation coefficients have the same sign, the bias (E[â * B ] − a b ) is positive; otherwise, the bias is negative. In our hypothetical example, one can expect the correlation coefficients between academic interest and ACH and ASC to be positive at school level. As school ACH increases, school academic interest level also increases, which, in turn, would elevate the level of ASC at the school level.
For b * B , the biasing effect of the omitted variable equals d YZ s ZM . In Equation (30), s ZM is proportional to the correlation between academic interest and ASC, while controlling for ACH at school level. In Equation (24), d YZ represents the Between effect of academic interest on FUT, while controlling for ACH and ASC. Again, as the magnitudes of the correlation coefficients increase, the amount of bias also increases. The direction of the bias depends on the sign of the correlation coefficients. Finally, the amount of bias in the coefficient c * B depends on the product of two quantities: d YZ (mentioned already) and s ZX . In Equation (30), s ZX is proportional to the correlation between academic interest and ACH, while controlling for ASC at school level.

Biased within coefficients
Biased (unadjusted) Within coefficients a * W , b * W , and c * W are as follows: where μ η X j and μ η M j are the expected values of η Xj and η Zj , respectively. As can be seen in (34), there exist two sources to potentially bias the estimate of coefficient a * W when academic interest, Z ij , varies at both school and student level. The first source of bias is the second term on the right-hand side of Equation (34), d aZ (d 0Z + d ZX μ η X j ), whose magnitude depends in part on (a) d aZ (see Equation [26]), which is a cross-level effect of school academic interest on the Within effect of student ACH on ASC, and (b) d 0Z + d ZX μ η X j , which is the (conditional) mean of school academic interest, predicted by school ACH. The conditional mean value can be interpreted as the expected ("typical") score of academic interest for a school with an average (typical) ACH score (i.e., μ η X j ). The second source of bias is quantified by the term l MZ l ZX , which represents the effect of the Within part of the omitted confounder, η Zij , on a W . The coefficient l MZ is proportional to the partial correlation between student ASC and student academic interest, while controlling for student ACH (see Equation [21]). The coefficient l ZX is proportional to the correlation between student academic interest and student ACH (see Equation [20]). Similarly, two potential sources of bias for b * W , as shown in Equation (35), are as follows. Note that the second term, d bZ (s 0 + s ZX μ η X j + s ZM μ η M j ), consists of d bZ and (s 0 + s ZX μ η X j + s ZM μ η M j ); d bZ (see Equation [27]) is a cross-level effect of school academic interest on the effect of student ASC on student FUT controlling for student ACH and academic interest; (s 0 + s ZX μ η X j + s ZM μ η M j ) is the (conditional) mean of school academic interest predicted by school ACH and ASC in Equation (30). The third term, l YZ q ZM , is a result of omitting the within part of student academic interest. The coefficient l YZ is proportional to the partial correlation between student FUT and academic interest, while controlling for student ACH and ASC (see Equation [22]); the coefficient q ZM is proportional to the partial correlation between student academic interest and ASC, while controlling for ACH (see Equation [29]).
Finally, as shown in Equation (36), two sources of bias influence c * W when academic interest is omitted from the model. The first source is d c Z (s 0 + s ZX μ η X j + s ZM μ η M j ); d c Z is the cross-level effect of school academic interest on the effect of student ACH on FUT, while controlling for student ASC and academic interest; (s 0 + s ZX μ η X j + s ZM μ η M j ) is the (conditional) mean of school academic interest predicted by school ACH and ASC (see Equation [30]). The second source of bias is the term l YZ q ZX , which quantifies the effect of omitting student academic interest on c W . The coefficient l YZ is the partial correlation between student FUT and academic interest controlling for student ACH and ASC (see Equation [22]); the coefficient q ZX is the partial correlation between student academic interest and ACH, while controlling for student ASC (see Equation [29]) .

Biased indirect effects
Of special importance in estimating a multilevel mediation model is to obtain estimates of the Within and Between indirect effects. The biased Within indirect effect equals the expected value of the biased random indirect effect (see the Appendix): In Equation (37), μ η Z j and σ 2 η Z j are the mean and variance of η Zj , respectively, of Z. Equation (37) shows that the expected value of a * j b * j is biased. The amount of bias E[a * j b * j ] − a W b W consists of the three following terms. First, the term a W l YZ q ZX + b W l ZX l MZ + l MZ l ZX l YZ q ZM quantifies the product of biases of a W and b W . It can be shown that the quantity a W b W + (a W l YZ q ZX + b W l ZX l MZ + l MZ l ZX l YZ q ZM ) equals the product of a * W and b * W , which are the expected values of a * j and b * j , respectively. The second term, k 1 μ η X j , is the expected value of the linear moderated effect of the between part of the omitted confounder (η Zj ) on the Within indirect effect. The term linear refers to the fact that the cross-level bias term, k 1 μ η X j , is a function of the first-order power of expected value of η Zj . Finally, the term d aZ d bZ (μ 2 η X j + σ 2 η X j ) shows the average quadratic moderated effect of the omitted confounder on the Within indirect effect. The quadratic moderated effect is a function of the second-order power of the expected value of η Xj . This quadratic effect is a result of the omitted variable that affects both the mediator and outcome variable.
Another result of omitting a Level 1 variable is that random coefficients a * j and b * j will be correlated in the misspecified model. The covariance between a * j and b * j is as follows: The values d aZ and d bZ represent the moderating effects of the Between part of the omitted confounder on the mediator and outcome variable, respectively. As the magnitude of the moderating effects becomes larger, so too does the covariance between the random coefficients a * j and b * j . In addition, we can rewrite the expected value of the biased random indirect effect in the misspecified model. Based on Equations (37) and (38), the expected value of the biased random indirect effect is Equation (39) shows that the expected value of the biased random indirect effect equals the product of the expected values of a * j and b * j plus the covariance between the biased random coefficients. This covariance quantifies the amount of between-cluster bias induced by omitting η Zj , the between-cluster component of the omitted variable. This covariance is a spurious covariance, which biases the mean of the random indirect effect as a result of omitting η Zj . More important, we can estimate this particular bias from the covariance between a * j and b * j . That is, a multilevel data structure provides enough information to estimate σ a * j ,b * j , which is part of the bias in (39). As a result, we can obtain a less biased estimate of the Within indirect effect by subtracting σ a *

Finally, the biased Between indirect effect is
where a * B and b * B are the biased Between coefficients in Equations (31) and (32), respectively.

Sensitivity analysis
The analytical results presented above are general in that they were not derived according to a specific distributional assumption about the latent proxy variable. To identify the bias calculation formula and make the numerical results tractable, we first assume that both Level 1 and Level 2 latent proxy variables have been scaled to have a mean of zero and standard deviation of one. This assumption also simplifies deriving formulas for biascorrected (adjusted) coefficients. The resulting simplified formulas are shown in the Appendix. Second, we assume that a plausible range of values for the correlation coefficients between Z and the variables X, M, and Y, r Z = (r ZX , r ZM , r ZY ) T are available at both Level 1 and Level 2.
Obtaining a plausible range of values is not trivial. Although, the alternative of considering the correlations is itself problematic because the assumption is that there is no omitted confounder. For Scenario 1, where the the single omitted confounder is known, but unmeasured, we believe the best strategy for obtaining plausible values is to use values reported in the relevant literature if they are available. If such values are not available, a second best strategy is to use the substantive knowledge of experts in the area. A Bayesian approach has been developed for eliciting plausible values of parameters from experts (see, e.g., Gill, 2015, Chapter 5). If no prior research is available, or the omitted confounders are unknown as in Scenario 2, a third strategy is to use general suggestions in the research area for "small, " "medium, " and "large" correlation values, r = ±.1, ± .3, and ± .5, respectively (Cohen, 1978). The permutations of the plausible values of the correlation coefficients would result in several r Z s. For example, r ZX , r ZM , and r ZY could be set to ±.1, ± .3, or ± .5, which would yield eight permutations of correlation values in r Z .
Given the range of the plausible values of correlations between the latent proxy variable and the observed values, we wrote computer code in R software (see the supplemental materials) that yields a range of bias-corrected estimates of Level 1 and Level 2 indirect effects. The range of the plausible values for the indirect effects is conditional on the hypothetical, but plausible, correlation values between the latent proxy variable and the observed variables. We use the bias-corrected estimates to answer the questions posed earlier. "At what values, if any, of r Z Note. SD = standard deviation; ACH = academic achievement; ASC = academic self-concept; FUT = career aspirations for the future. Total number of schools was N = ;  students were sampled from each school resulting in a total sample size of n =  students.
would the conclusions about the estimates of the Level 1 and Level 2 indirect effects change in a meaningful way?" Another potential question might be "What is the effect of small to moderate values of the correlation of the omitted confounder with the observed values on the estimates of indirect effects?" The point of a sensitivity analysis in this context is to assess whether the conclusions drawn might be highly sensitive or robust to various combinations of the range of plausible values for the correlations.

Empirical example
We use a hypothetical empirical example to conduct a sensitivity analysis and to show the practical implications of the analytical results. Consider the 1 → 1 → 1 model for the hypothetical example used throughout the manuscript. Data for the empirical example were simulated using the results from Nagengast and Marsh (2012). The Between and Within sample correlations as well as the standard deviation (SD) are shown in Table 1. In general, it is expected that the school-level and student-level relationships between ACH and ASC differ (Marsh et al., 2014). As previously mentioned, this contrast is termed BFLPE (Marsh, 1987). A 1 → 1 → 1 model with random intercepts and slopes according to Equations (1)-(10) was estimated using Mplus (Muthén & Muthén, 1998-2012 5 . The results are shown in Table 2. First, we assume that the model is correctly specified in that Assumptions 1-5, especially the no-interaction and no-omitted-confounder assumptions, are met. The RMediation package (Tofighi & MacKinnon, 2011) was used to calculate the confidence intervals (CIs) using the analytical solution to the distribution-of-theproduct method as this method takes into account the potential skewness and high kurtosis in the sampling distribution of the indirect effect (MacKinnon, Lockwood, Hoffman, West, & Sheets, 2002). In this case, the Level 2  Data file as well as Mplus input and output files are available in the supplemental materials. Table . Results for empirical example when the  →  →  is assumed to be correctly specified (N = , n = , ).
Note. SE = standard error; CI = confidence interval; X = academic achievement (ACH); M = academic self-concept (ASC); Y = career aspirations for the future (FUT); N = number of schools; n = total number of students; subscripts "W, " "B, " and "j" denote "Within, " "Between, " and "Student j, " respectively. However, because the values of ACH and ASC were not randomized, other confounders could potentially bias the estimates of the indirect effects at Levels 1 and 2. For example, there is some evidence that student academic interest is correlated with both ACH and ASC (Marsh et al., 2005). We then conducted a sensitivity analysis to assess the robustness of the estimated indirect effects to the potential violation of the no-omittedconfounder assumption. For simplicity in this illustration, we assume that student academic interest is the single unmeasured, but known, confounder, as in Scenario 1. We also assume that this confounder could potentially influence the observed variables at both levels, thus biasing the relationships at each level and across levels. As a result, we can use a single latent proxy variable, Z, to model the effect of the confounder on the indirect effect estimates.
In assessing confounder bias, we first look at the potential cross-level bias. Given that σ a j ,b j = 0.004, p = .8, we assumed for simplicity in this illustration that the crosslevel moderating effect of the Between part of the confounder on the Level 1 indirect effect is zero in the population. The assumption of no-cross-level bias is based on the analytical result in Equation (38). In Equation (38), we showed that nonzero σ a j ,b j represents the moderating effects of η Zj on the Level 1 indirect effect. When the  Table ).
covariance is zero, however, we can conclude that potential Level 2 omitted confounder does not moderate the Level 1 indirect effect.
Next, we calculate the bias-corrected (adjusted) estimates of Level 1 and Level 2 indirect effects using the analytic results. As previously discussed, the bias correction is a function of the correlation between the latent variable (a proxy for confounders) and X, M, and Y at Levels 1 and 2. As the correlation values take on different plausible values, we calculate a range of bias-corrected estimates for the indirect effects. For illustrative purposes, we choose the following plausible range of values for r ZX = 0, .1, and .3. One may also choose additional values, but for this illustration, these three values revealed enough information as will be discussed; r ZM and r ZY ranged between −.5 and .5. We did not display the negative values of r ZX , which would produce the same range of values for the bias-corrected indirect effects.
To facilitate the interpretation of the results of the sensitivity analysis, we created sensitivity contour plots. 6 A sensitivity contour plot depicts the values of the bias-corrected (adjusted) indirect effect across the range of the plausible values of the correlation between latent proxy variable (Z) and ASC (M) and correlation between latent proxy variable and FUT (Y) at each level of correlation between the latent proxy variable and ACH (X). The ranges of the x-axis and y-axis demonstrate the range of plausible values for r ZM and r ZY , respectively. Contour lines on each plot show different values of bias-corrected (adjusted) indirect effects. Each contour line connects all of the points with the same magnitude of the biascorrected indirect effect. That is, each line contains all of the combinations of r ZM and r ZY values that produce the same bias-corrected indirect effect. Each contour line is also labeled with a numerical value of the bias-corrected indirect effect. Two adjacent contour lines display distinct values of bias-corrected indirect effects. Note that not all  We wrote R code to create sensitivity contour plots and produce numeric ranges of adjusted indirect effects in the supplemental materials.  Table ).
the values of the bias-corrected indirect effects are shown. Instead, to show the values of the bias-corrected indirect effect between adjacent contour lines, each region is filled with a color. There is a legend to the right of each plot showing a correspondence between the greyscale gradation (or spectrum of colors) and the range of biascorrected indirect effect values. In interpreting the contour plots, one may be interested in the range of the bias-corrected indirect effects in the areas located in the middle and corners of each plot. The middle area contains the biascorrected values corresponding to zero to "small" values of r ZM and r ZY . One might conclude that the contour lines in the middle area represent the bias-corrected values when the magnitude of correlation between the confounder and M and Y ranges from zero to small. The corners, on the other hand, contain the bias-corrected values for a combination of the larger correlation values (e.g., .4 to .5) for r ZM and r ZY . Contour lines in these areas might represent the bias-corrected indirect effect values when the magnitude of correlation between the confounder and M and Y ranges from "medium" to "large. " Figure 3 shows three sensitivity contour plots corresponding to three plausible values for the correlation between the latent proxy variable and ACH: r ZX = 0, .1, and .3. The value labels on the contour lines of the plot for r ZX = 0 show that the bias-corrected indirect effect ranges from approximately 0.015 to 0.035. For r ZX = .1 and .3, we can see both positive and negative values of the bias-corrected indirect effects. The R code also provides a more precise range of the bias-corrected indirect effect values for each plot; the bias-corrected values of the Between indirect effect ranged from 0.017 to 0.032 for r ZX = 0, −0.087 to 0.152 for r ZX = .1, and −0.372 to 0.447 for r ZX = .3 given the range of the plausible values considered. Recall that the estimate of the indirect effect of school ACH on school FUT was 0.019, 95% CI [−0.150, 0.188]), when we assumed the no-omitted-confounder assumption held. Comparing these (potentially biased) values of the indirect effect to the potential range of the values of the bias-corrected indirect effect produced in the sensitivity analysis, it is unlikely that the conclusion about the nonsignificant Level 2 indirect effect would have changed had we included the omitted confounder.
For the Level 1 indirect effect, the sensitivity contour plots are presented in Figure 4. The value labels on contour lines of the plots for r ZX = 0 and .1 show that the biased-corrected indirect effects are above 0. For r ZX = .3, however, the labels of the contour plot show the values ranging from 0 to .20. Recall that the contour plot itself shows only select values of the bias-corrected indirect effect. Again, we used the R code to get a more precise range of the bias-corrected indirect effect. The biascorrected estimates of the indirect effect ranged from 0.017 to 0.085 for r ZX = 0, .011 to .121 for r ZX = .1, and −0.016 to 0.226 for r ZX = .3. Recall that when we assumed that the no-omitted-confounder assumption held, the (potentially biased) estimate of the indirect effect was 0.037, 95% CI [0.014, 0.062]). Given the bias-corrected estimates from the sensitivity analysis were above 0, the significant result for the Level 1 indirect effect of ACH on FUT appears to be robust when the correlation between Z (a proxy for confounders) and X ranges from 0 to .1 (small) and correlations between Z and M and Z and Y range from 0 to .5 (large); that is, given the plausible range of values for the correlation with the confounder, the values of the Level 1 indirect effect would have been positive and statistically significant had the potential confounder been included in the model. However, when the correlation between Z and ACH is .3 (medium), the range of the values for the bias-corrected indirect effect appears to contain both positive and negative values. Upon further investigation, the ranges of r ZM and r ZY that caused zero or negative values of the biascorrected indirect effect were .43 to .50 and −.5 to .5, respectively. This indicates that if one were to assume that Z is "moderately" correlated with ACH and "strongly" correlated with ASC, then the conclusion about positive indirect effect would likely be invalid. Otherwise, the conclusion about the Level 1 indirect effect appears to be robust against the biasing effect of the potential confounder. In the present case, the likelihood that such relationships might exist is a substantive judgment based on prior research.

Conclusion
With the growing popularity of multilevel mediation analysis in applied research, it is critical to probe the effects of underlying assumptions. To draw valid inference about an indirect effect in a mediation model, a set of specifications must be met. The first contribution of our manuscript is that we extended the specification assumptions from the single-level mediation literature to the multilevel mediation analysis of a 1 → 1 → 1 model. One of the specification assumptions is the noomitted-confounder assumption, which means that there are no common causes of hypothesized causal relationships in the mediation model. When the specification assumptions 1-5 hold, one can compute an unbiased estimate of each of the indirect effects. On the other hand, when the no-omitted-confounder assumption is violated, inference about the indirect effects can be severely biased; thus, the results are misleading. A formidable challenge is that the no-omitted-confounder assumption is not testable. Thus, previous research recommends assessing the extent to which potential violation of the no-omitted-confounder assumption might invalidate or alter the conclusions about the indirect effects actually observed.
To this aim, we proposed a framework to analytically examine the potential biasing effects of omitted Level 1 confounder(s) that were correlated with all of the observed variables at both Level 1 and Level 2 in a twolevel 1 → 1 → 1 mediation model with random intercepts and slopes. We discussed two scenarios about the types of the omitted confounder(s) at Level 1 that might arise in practice: (a) a single unmeasured, but known, confounder, and (b) one or multiple unknown confounders. We used a single latent proxy variable to model these two types of confounders.
Our analytic results show that omitting Level 1 confounder(s) can yield misleading results about key quantities of interest such as Level 1 and Level 2 indirect effects.
One key finding was that omitting Level 1 confounder(s) can exert compounded, biasing effects on the estimation and interpretation of the indirect effects. In addition, we showed that potential biasing effects of the omitted confounder(s) on the Level 1 estimates can be further decomposed into (a) a biasing effect due to the Within part and (b) a biasing cross-level effect due to the Between part of the omitted confounder(s). A second key result was that we presented formulas that quantified the potential biasing effects of the omitted confounder(s) in terms of the (partial) correlation between the omitted and observed variables in the model. The analytic result shows that the bias-corrected indirect effects, the estimates of the true indirect effects had the confounder(s) been included in the model, are a function of correlations between latent proxy variable and the observed variables.
Our results indicate that when the no-omittedconfounder assumption is violated, the covariance between a j and b j , σ a j ,b j becomes nonzero. As a result, the equality between two methods of calculating Level 1 indirect effect does not hold: c W − c W = a W b W . On the other hand, when the specification assumptions hold, the equality holds: c W − c W = a W b W . In addition, a nonzero covariance term σ a j ,b j might signal the presence of an omitted confounder. The reason we cannot make this a sufficient condition claim is that there might exist sources other than an omitted confounder that might cause a nonzero (spurious) covariance between a j and b j . By spurious, we mean an extraneous covariance not accounted for by the posited 1 → 1 → 1 model. In the specification assumption section, we have identified three sources that might cause nonzero spurious covariance: (a) measurement error, (b) common method effects, and (c) omitted confounders. If we assume that the three general categories cover all possible sources of spurious covariance, we can make a more definitive statement about potential sources of nonzero covariance between a j and b j . More specifically, if we can reasonably rule out the existence of the method and measurement error effects in a study, then nonzero covariance can signal the influence of an omitted confounder at Level 2. Note that only the Between part of Level 1 confounder causes nonzero covariance between a j and b j . An omitted confounder may exist solely at the Within level. Zero covariance between a j and b j cannot rule out the existence of the Within part of omitted confounder(s).
Then, for applied researchers, we developed a sensitivity analysis that assesses the robustness of the indirect effects to the violation of the no-omitted-confounder assumption. This approach integrates the investigation of the potential biasing effects of Level 1 omitted variables into the existing framework of multilevel mediation analysis including the multilevel structural equation modeling (MSEM). The initial estimate of the model parameters can be performed using any available software package capable of fitting a 1 → 1 → 1 model (e.g., Mplus). We also wrote computer code in R that conducts sensitivity analysis using the initial estimates and produces the sensitivity contour plots given the range of the plausible values for the correlation between a latent proxy variable and observed values provided by the researcher. The sensitivity contour plot illustrates potential bias-corrected (adjusted) estimates of the indirect effect across a range of plausible values for the correlations. In some cases, the plot may clearly indicate that the indirect effect is robust to plausible magnitudes of violation of the noomitted-confounder assumption. In other cases, the plot may clearly indicate that the indirect effect is not robustit will no longer be significant, even given small values of r ZX , r ZM , and r ZY correlations. Finally, in some cases, as in the present illustration, indirect effect estimates will be robust over a wide range of, but not all, potential values of confounding. In such cases, the likelihood that such relationships might exist will be a substantive judgment based on prior research.
The results of our study do not depend on a specific centering strategy for the predictors. However, we recommend centering the predictors according to the results of Enders and Tofighi (2007). Indeed, the choice of centering depends on research questions (Enders, 2013;Enders & Tofighi, 2007). Second, it is important to note that one can also use centering at the grand mean 2 (CGM2) strategy to conduct sensitivity analysis; CGM2 centers a Level 1 predictor at the grand mean while adding the cluster means on the predictor to the model at Level 2. The results of our study hold because these two centering strategies, CWC2 and CGM2, are "equivalent" in that we are (essentially) estimating the same model, only shifted, which would result in the same sample likelihood function (Kreft et al., 1995). That is, the coefficient estimates in CGM2 are a linear transformation of the coefficient estimates in CWC2 and vice versa. In CGM2, the coefficient associated with the cluster mean predictor is the contextual effect, an estimate of the differential effect of Between and Within effects (Blalock, 1984); however, in CWC2 the coefficient is an estimate of the Between effect. In both CWC2 and CGM2, the coefficient associated with the centered Level 1 predictor is an estimate of the Within effect. Because of the equivalence between the two centering strategies, the results of our study hold using either centering strategy.
Although the present work focused on the 1 → 1 → 1 model, the present analysis can be extended to other mediation models in which at least one of the variables is measured at Level 2. Examples of the models that can be fit using available software are the 1 → 1 → 2 and 2 → 1 → 1 mediation models. These additional mediation models are simpler, in terms of the number of fixed and random effects that must be considered, than the 1 → 1 → 1 model considered here. One complication may occur because one of the variables is measured at Level 2 in these models: the researcher may not be able to obtain separate estimates of the Level 1 and Level 2 effects (Pituch & Stapleton, 2012). For such models, the current analytic results may be used to examine the biasing effects of the confounder on the indirect effect. Finally, the present results can be specialized to 1 → 1 → 1 models in which randomization of X occurs at Level 1. Randomization greatly simplifies the sensitivity analysis because r ZX can be assumed to be zero.
Additional topics remain for future study. One extension of the current analytical results would be to include additional covariates that (incompletely) control for confounders in the model. The addition of covariates that are theoretically related to X, M, and Y can enhance the causal interpretation of the model. The sensitivity analysis procedure may still be applied by extending the method suggested by Mauro (1990). One limitation of the current study is that we only considered the analysis of two-level data that is commonly found in applied settings. However, the present analysis can be extended to consider multilevel mediation models with data structures involving three levels of nesting. Another limitation is that we assumed the omitted confounder(s) are linearly correlated to the observed variable in the 1 → 1 → 1 model. In addition, we did not consider a case in which there might exist multiple known, but unmeasured, confounders, as in Scenario 3. Extending the present analysis to address these cases remains a topic for future study.

Article information
The expression in Equation (A9) provides an estimate of the Within indirect effect at the grand mean of Z ij . It should be noted that μ Z i j = μ η Z j = 0. This result holds because of the orthogonal decomposition of Z ij into latent variables in (19); μ η Xi j = 0 and Z ij is a CGM score. Thus, we can use the grand mean of Z ij and η Zj interchangeably.

Simplified analytic results for sensitivity analysis
To identify the bias calculation formula and make the numerical results tractable, we first assume that both Level 1 and Level 2 latent proxy variables have been scaled to have a mean of zero and standard deviation of one. This assumption also simplifies deriving formulas for the bias-corrected (adjusted) coefficients. The simplified results for the analytical results for the Between coefficients remain unchanged: The expected values of the random coefficients are simplified as follows: In addition, the covariance between a * j and b * j is simplified as follows: Finally, we can calculate the simplified bias-corrected (adjusted) results by isolating the desired quantities on the right side of Equations (A10)-(A15). For example,