Sensitivity analysis for unobserved confounding in causal mediation analysis allowing for effect modification, censoring and truncation

Causal mediation analysis is used to decompose the total effect of an exposure on an outcome into an indirect effect, taking the path through an intermediate variable, and a direct effect. To estimate these effects, strong assumptions are made about unconfoundedness of the relationships between the exposure, mediator and outcome. These assumptions are difficult to verify in a given situation and therefore a mediation analysis should be complemented with a sensitivity analysis to assess the possible impact of violations. In this paper we present a method for sensitivity analysis to not only unobserved mediator-outcome confounding, which has largely been the focus of previous literature, but also unobserved confounding involving the exposure. The setting is estimation of natural direct and indirect effects based on parametric regression models. We present results for combinations of binary and continuous mediators and outcomes and extend the sensitivity analysis for mediator-outcome confounding to cases where the continuous outcome variable is censored or truncated. The proposed methods perform well also in the presence of interactions between the exposure, mediator and observed confounders, allowing for modeling flexibility as well as exploration of effect modification. The performance of the method is illustrated through simulations and an empirical example.


Introduction
To estimate causal direct and indirect effects of an exposure on an outcome a key assumption is unconfoundedness of the relationships between exposure, mediator, and outcome. Since unconfoundedness is difficult to verify in a given situation results should be accompanied by a sensitivity analysis to gauge the impact of violations on the estimated effects (Rosenbaum 2010, Chap. 14).
In mediation analysis the focus has been predominantly on sensitivity against violations of no unobserved mediator-outcome confounding. The argument used is that confounding related to the exposure could be handled by randomization or by adjusting for a ''sufficiently rich'' set of pre-exposure confounders, while confounding related to the mediator is more difficult to design or adjust away. However, in many applications the exposure cannot be randomized and it is often difficult to guarantee that a sufficiently rich set of pre-exposure confounders has been adjusted for.
Different approaches have been suggested for sensitivity analysis to unobserved mediator-outcome confounding. Among these are methods based on correcting estimates and confidence intervals (CIs) using a bias factor based on the specification of the relationships between the unobserved confounder and the mediator, outcome and/or exposure (VanderWeele 2010;Hafeman 2011;le Cessie 2016). An alternative approach using the correlation between the error terms in the parametric regression models for the mediator and outcome as the sensitivity parameter was suggested by Imai et al. (2010a) and implemented in the R (R Core Team 2017) package mediation (Tingley et al. 2014(Tingley et al. , 2019. This approach involves deriving expressions for the direct and indirect effects that take this correlation into account. It offers sensitivity analysis to unobserved mediatoroutcome confounding for continuous mediators and outcomes as well as when either the mediator or the outcome is binary, with the caveat that the binary outcome model cannot include any exposure-mediator interactions. A similar approach was suggested by Lindmark et al. (2018) for cases when both the mediator and outcome are binary and probit models are used for estimation. Instead of deriving expressions for the direct and indirect effects this approach incorporates correlations between error terms of the mediator, outcome and exposure assignment models into the estimation of the model parameters upon which the direct and indirect effects estimates are based. This approach is able to take into account not only mediator-outcome confounding but also exposuremediator and exposure-outcome confounding. It is also flexible in that a sensitivity analysis can be performed also in the presence of interactions involving the exposure, mediator and observed confounders. The latter allows richer model specification and also enables performing sensitivity analyses in situations where the investigation of effect heterogeneity in different subpopulations is of interest.
Estimation of direct and indirect effects is further complicated when there is censoring or truncation of the data. In the context of structural equation models for estimation of mediation effects Wang and Zhang (2011) showed that censoring leads to both reduced accuracy and precision, especially when it is the outcome variable that is censored. They suggested a tobit mediational model to account for censored data in the estimation of effects. However, as their approach was not in the context of causal mediation analysis no attention was given to assumptions about unconfoundedness or related sensitivity analyses. The mediation package (Tingley et al. 2014(Tingley et al. , 2019 allows estimation of causal mediation effects when the outcome is censored based on a tobit model but does not provide an accompanying sensitivity analysis method. Aside from these examples, most research into methods for mediation analysis in the presence of censoring of the outcome has taken place in the context of time-to-event outcomes (see e.g. Lange and Hansen (2011);VanderWeele (2011) and VanderWeele (2015), Chap. 4) including suggestions for sensitivity analyses to unobserved mediator-outcome confounding (Tchetgen Tchetgen 2011;VanderWeele 2013). The related but more severe issue of truncation, i.e. when the outcome is not recorded at all for certain values has to our knowledge not been examined within the mediation literature.
In this paper we extend the sensitivity analysis method to unobserved mediatoroutcome confounding and confounding involving the exposure for parametric estimation of direct and indirect effects introduced in Lindmark et al. (2018) to include cases with continuous mediators and/or outcomes. We also suggest sensitivity analysis methods for unobserved mediator-outcome confounding for the more complicated settings when the outcome is censored or truncated, building on the tobit model for censored outcomes (Tobin 1958) and its equivalent for truncated outcomes (Hausman and Wise 1977). We illustrate the performance of the method through simulations and present an empirical example. The approach is implemented in the R package sensmediation (Lindmark 2019) with the exception of the suggested methods for censored or truncated outcomes where we provide R code for the analyses performed in this paper.
The paper is structured as follows. In Sect. 2.1 direct and indirect effects are defined using the counterfactual framework for mediation (Robins and Greenland 1992;Pearl 2001) and the assumptions required for identification are presented. In Sect. 2.2 the general idea behind the sensitivity analysis method is presented and parametric estimators of direct and indirect effects with accompanying sensitivity analyses for different combinations of continuous and binary mediators and outcomes are suggested. In Sect. 2.3 corresponding results are presented for cases where the outcome is censored or truncated. The simulation scenarios are outlined in Sect. 3.1 with simulation results in Sect. 3.2 and an empirical example in Sect. 3.3. Finally, we summarize the findings and discuss limitations and further developments in Sect. 4.

Identification and assumptions
Let Z be an exposure, Y an outcome, and M a mediator of the exposure-outcome relationship (see Fig. 1).
Let M i ðzÞ denote the potential value of the mediator for individual i under exposure level z, Y i ðz; mÞ, the potential outcome for individual i under exposure level z and mediator level m and Y i ðz; M i ðz 0 ÞÞ, the composite potential outcome if the exposure Z i were set to the value z and the mediator M i were set to its value under exposure level Z i ¼ z 0 .
We define the natural direct effect contrasting two exposure levels z 1 and z 0 , as the effect on Y of changing Z from z 0 to z 1 if the mediator were allowed to vary as it would naturally if all individuals in the population were under exposure level z.
The natural indirect effect is defined as the effect on Y when, keeping the exposure fixed at z in the population, allowing the mediator to change from its potential value when z 0 to its potential value when z 1 . If we make a composition assumption, i.e. that Y i ðzÞ ¼ Y i ðz; M i ðzÞÞ (Van-derWeele and Vansteelandt 2009), the total effect TE z 1 ;z 0 ¼ E Y i ðz 1 Þ À Y i ðz 0 Þ ½ can be decomposed as either TE z 1 ;z 0 ¼ NDE z 1 ;z 0 ðz 0 Þ þ NIE z 1 ;z 0 ðz 1 Þ or TE z 1 ;z 0 ¼ NDE z 1 ;z 0 ðz 1 Þ þ NIE z 1 ;z 0 ðz 0 Þ. Using terminology introduced by Robins and Greenland (1992) the former decomposition is into the pure natural direct and total natural indirect effect and the latter decomposition into the total natural direct and pure natural indirect effects.
Often we have a binary exposure taking the values Z ¼ 1 if exposed and Z ¼ 0 if unexposed.
The most common decomposition is then To identify natural direct and indirect effects from observed data, we assume consistency, so that for an individual i with observed exposure Z i ¼ z we have that and Vansteelandt 2009). Together with the composition assumption this implies We also assume no interference, i.e. that the exposure level of one individual does not have an effect on the mediator or the outcome of another individual (De Stavola et al. 2015). Finally, we make crucial assumptions about unconfoundedness: Fig. 1 A directed acyclic graph showing the relationships between exposure Z, mediator M, and outcome Y Assumption 1 Sequential ignorability (Imai et al. 2010a) 1.
i.e., there is no unobserved confounding of the exposure-mediator and exposure-outcome relationship given the observed preexposure covariates X i . 2.
i.e., given X i and the observed exposure Z i there is no confounding of the mediator-outcome relationship. where 0\P , and all x 2 X (the support of X) and m 2 M (the support of M).
Note that Assumption 1 implies a so-called cross-world independence assumption (see e.g. VanderWeele et al. 2014), i.e. independence between the counterfactual outcome under exposure level Z ¼ z 0 and mediator level M ¼ m and the counterfactual mediator under exposure level Z ¼ z, where z 0 and z are possibly different values. Since in reality different values of the exposure cannot be observed simultaneously this assumption is difficult to verify empirically. The cross-world independence assumption is violated in cases where there is intermediate confounding, i.e. mediator-outcome confounders affected by the exposure.
If these assumptions are fulfilled the natural direct and indirect effects conditional on the covariates are identified by (Pearl 2001) For continuous mediators we replace the sums and probabilities in (1) and (2) with integrals and densities. By marginalizing (1) and (2) over x we obtain the NDE z 1 ;z 0 z ð Þ and NIE z 1 ;z 0 z ð Þ, the natural direct and indirect effects at the population level.
Here we use a parametric approach where the natural direct and indirect effects are estimated by specifying parametric regression models for the outcome and mediator and rewriting (1) and (2) as functions of the regression parameters. The resulting estimators are consistent given that the previously outlined assumptions are fulfilled and the regression models are correctly specified. In the following section we derive estimators of the natural direct and indirect effects for combinations of binary and continuous mediators and outcomes (for the combination binary mediator and outcome, see Lindmark et al. (2018)).

Estimators and sensitivity analysis in the absence of censoring or truncation of the outcome
We specify a parametric regression model for the mediator conditional on the exposure and observed covariates. For a continuous mediator we specify a linear regression model where the g i are i.i.d. (independent and identically distributed) with zero mean and standard deviation r g . For a binary mediator we specify a probit regression model with where g Ã i $ i:i:d: Nð0; 1Þ.
We also specify a parametric regression model for the outcome conditional on the exposure, mediator and observed covariates. For a continuous outcome we specify where the n i are i.i.d.with zero mean and standard deviation r n . For a binary outcome we specify with n Ã i $ i:i:d: Nð0; 1Þ. In Table 1 expressions for the natural direct and indirect effects are presented for different model combinations, first when both mediator and outcome are continuous, then when the mediator is binary and the outcome continuous and lastly when the mediator is continuous and the outcome binary. Note that these are more general versions of previously derived expressions, see Imai et al. (2010a), adding interactions between the covariates and exposure and mediator to the regression models used to allow for moderated mediation, i.e. different direct and indirect effects for different covariate levels. The natural direct and indirect effects are estimated by fitting the mediator and outcome models using maximum likelihood (ML) and plugging the estimated parameters into the appropriate expressions in Table 1. Approximate standard errors of the effects can be obtained using the delta method (Oehlert 1992).

Sensitivity analysis
The sensitivity analysis is presented for mediator-outcome confounding (U 2 in Fig. 2) but can be modified to exposure-mediator (U 1 ) or exposure-outcome (U 3 ) confounding by replacing the mediator model with a model for the exposure assignment conditional on the covariates and the outcome model with the mediator model, or by replacing the mediator model with an exposure model, respectively. For details see Lindmark et al. (2018). We assume that the error terms in the mediator and outcome models are bivariate normal with correlation q. If there is unobserved mediator-outcome confounding then q 6 ¼ 0, otherwise q ¼ 0. The sensitivity analysis is performed by deriving the joint likelihood for M and Y as a function of the regression parameters and q. In Table 2 the log-likelihoods for a sample of n units (i ¼ 1; :::; n) derived for different model combinations are presented. We cannot estimate q from the observed data without further assumptions (Imai et al. 2010b) and instead proceed with a modified maximum likelihood (ML) procedure, where the log-likelihood is maximized with regards to the regression parameters for a fixed value of the correlation, q ¼q. The sensmediation package (Lindmark 2019) uses functions from the maxLik (Henningsen and Toomet 2011;Toomet and Henningsen 2015) package for the maximization. The default maximization method is the Newton-Raphson algorithm which utilizes analytical gradients and Hessians of the log-likelihood functions.
The resulting parameter estimatesbq ð Þ orb Ãq ð Þ andĥq ð Þ orĥ Ãq ð Þ (plusr gq ð Þ for a continuous mediator and binary outcome) are then plugged into the expressions for the NDE z 1 ;z 0 ðz; xÞ and NIE z 1 ;z 0 ðz; xÞ (Table 1). This gives estimates of the conditional natural direct and indirect effects under a given level of unobserved mediator-outcome confounding, d NDE z 1 ;z 0 ðz; x;qÞ and d NIE z 1 ;z 0 ðz; x;qÞ. Estimates of the marginal natural direct and indirect effects under given levels of confounding are given by averaging these estimated conditional effects over the study The results of the sensitivity analysis can be presented in different ways. One is to report the results over a range of the sensitivity parameter. This range can be defined using subject matter knowledge about the probable nature of the unobserved confounding, e.g. whether or not an unobserved confounder is expected to affect both the mediator and the outcome in the same directions and thus induce a positive error term correlation. In the absence of such prior knowledge a wide range encompassing both negative and positive correlations can be used. The results can be summarized through plots of point estimates and CIs and/or so called uncertainty intervals (UIs) (Vansteelandt et al. 2006;Genbäck et al. 2015Genbäck et al. , 2018, the union of all 100 Â ð1 À aÞ% CIs over the range of the sensitivity parameter. An alternative, or complement, to these is to report the values of the sensitivity parameter where the 100 Â ð1 À aÞ% CIs include 0, i.e. where the effect is no longer significant at an a level of significance.

Estimation and sensitivity analysis in the presence of censoring or truncation of the outcome
Here we present estimation methods and sensitivity analyses to mediator-outcome confounding when we have a continuous mediator and the outcome is either left censored or left truncated, i.e. where censoring/truncation occurs for values of the (3) and (5) (3) and (6) Á ð Þ denotes the pdf of a bivariate normal distribution with zero mean vector and covariance matrix R ¼ r 2 g qr g r n qr g r n r 2 n ! : b / Á ð Þ and U Á ð Þ denote the standard normal pdf and cumulative distribution function, respectively. c See Appendix A for the derivation of the joint mediator and outcome distribution. d The joint mediator and outcome distribution is derived as in Appendix A outcome variable that are below a certain point. The methods presented here can be used also for right censoring/truncation, i.e. where censoring/truncation occurs for values of the outcome variable that are above a certain point. This can be accomplished by multiplying the right censored/truncated outcome variable by -1, thus transforming it into a left censored/truncated variable.
These methods are not currently implemented in the sensmediation package but analytic gradients of the log-likelihoods are provided in appendices to facilitate implementation of the optimization to obtain ML estimates. We also provide the code used to perform the analyses in the simulation study which may be adapted to other applications. Here, the optimization was performed using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method implemented in the maxLik (Henningsen and Toomet 2011; Toomet and Henningsen 2015) package, as no analytic Hessian was derived.
For comments on adapting the methods presented to sensitivity analyses of unobserved confounding involving the exposure, see Sect. 4.

Sensitivity analysis unobserved mediator-outcome confounding, censored outcome
Assume that we have a continuous mediator and outcome that follow models (3) and (5), but that we observe Y i ¼ maxðY i ; tÞ, i.e. left censoring at t. To estimate direct and indirect effects the mediator model could be fitted using e.g. OLS or ML while the regression parameters in the outcome model could be estimated using e.g. tobit regression (Tobin 1958). To assess the sensitivity of the estimated effects to mediator-outcome confounding we again assume that the error terms g i and n i are bivariate normal with correlation q. The joint distribution of the mediator and outcome is then given by That is, for observations that are not censored the joint distribution is simply a bivariate normal distribution while for observations that are censored we have: with the resulting density derived as in Appendix A. The joint log-likelihood is To obtainĥq ð Þ andbq ð Þ, (7) is maximized for a fixed q ¼q (gradients are provided in Appendix B). These estimates are then plugged into the expressions for the natural direct and indirect effects in Table 1, model combination (3) and (5) to obtain estimates under a given level of unobserved mediator-outcome confounding and censoring.

Sensitivity analysis unobserved mediator-outcome confounding, truncated outcome
Truncation is a more complicated problem than censoring as observations are completely missing, meaning that truncation of the outcome also leads to missing mediator values. To reduce the complexity here we simplify the models (3) and (5) for M and Y to only include main effects but the results can be extended to models including interactions. The mediator and outcome models used here are: The expressions for the natural direct and indirect effects under these simpler models are given by Now assume that we only observe Y i [ t, i.e. truncation at t. Since truncation of the outcome also leads to missing mediator values we simultaneously estimate the parameters in the mediator and outcome regression models. Assume that g y i and n y i are bivariate normal with correlation q. Then, ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r 2 n y þ h y2 2 r 2 g y þ 2h y 2 qr n y r g y q 0 B @ 1 C A; (See Appendix C for derivation). The joint log-likelihood for the mediator and outcome is given by By maximizing (12) (for gradients see Appendix D) for q ¼ 0 we obtainĥ y andb y under truncation. The relevant parameters can then be plugged into (10) and (11) to obtain estimates of the natural direct and indirect effects. For sensitivity analysis we maximize (12) for non-zero q ¼q to obtainĥ yq ð Þ andb yq ð Þ and in turn d NDE z 1 ;z 0 ðz;qÞ and d NIE z 1 ;z 0 ðz;qÞ.

Simulation scenarios and data generation
To demonstrate the performance of the proposed approach a simulation study was performed. For each replicate, observations of an exposure, an outcome, a mediator and an observed confounder affecting the exposure, mediator and outcome were generated (R code is found at by https://github.com/anitalindmark/Sensitivity_ analysis). Five scenarios were investigated (see Table 3 for a summary). In scenarios a, b, d and e the data generating mediator and outcome models contained all interactions involving the exposure and (for the outcome model) mediator. In scenario c where the outcome was truncated data were generated from models containing only main effects.
The regression coefficients used to generate the mediators and outcomes were selected to yield approximately equal effects within each scenario for comparability. For scenarios a, b, d and e the true effects were obtained by using the data generating regression coefficients in the expressions in Table 1. To obtain marginal effects Monte Carlo integration was performed by generating a very large number (n ¼ 1 Â 10 9 ) of values of the observed covariate, calculating the effects conditional on these values, and averaging the effects. For scenario c the true effects were given by (10) and (11), i.e. by the true regression coefficient for the exposure in the outcome model and the product of the true regression coefficient for the exposure in the mediator model and the true coefficient for the mediator in the outcome model, respectively. True values in all scenarios were NIE 1;0 ð1Þ % 0:041 and NDE 1;0 ð0Þ % 0:038.
In each scenario a-e mediator-outcome confounding was induced by correlating the error terms of the data generating models, with q ¼ 0:5. For scenarios a, d and e separate simulations with exposure-mediator and exposure-outcome confounding were performed. For these simulations confounding was induced by correlating the error terms in the model used to generate the exposure and the model to generate the mediator and outcome, respectively. Samples of size n obs ¼ 500; 1000; 5000 were generated 2000 times from each scenario. In each of the 2000 replicates effects and standard errors were estimated based on two values of the sensitivity parameter:q ¼ 0 (assuming no unobserved confounding) andq ¼ 0:5 (the true value). The sensmediation package was used for estimation. For censoring and truncation separate functions for the optimization, log-likelihoods and gradients were implemented (code for these are found at https://github.com/anitalindmark/Sensitivity_analysis). Functions from the sensmediation package were then used to calculate the effects and standard errors. The input outcome model for scenario b (censored outcome) was estimated using the tobit function from the AER package (Kleiber andZeileis 2008, 2020).

Simulation results
Results were summarized using the rsimsum package (Gasparini 2018) and are presented according to recommendations in Morris et al. (2019). The performance measures used are the bias and empirical coverage rate of 95% CIs over the 2000 replicates. In addition, the SEs of the effects estimated using the delta method are compared to empirical SEs over the 2000 replicates using the relative % error:  (Figs. 3,4,5,8,9,10,11 and 12 in Appendix E), where dots represent the estimated performance measure, with a line from the dot to the target value of that performance measure. The 95% CIs are represented by parentheses and thus parentheses not enclosing the target value indicate evidence that the performance measure does not meet the target. Results for mediator-outcome confounding and scenarios a-e are summarized in Figs. 3, 4 and 5. For all scenarios, not taking into account unobserved confounding (i.e. usingq ¼ 0) led to substantial bias (Fig. 3a). Note that since the total effect is not affected by mediator-outcome confounding and is given by summing the natural direct and indirect effects the biases of the d NDE 1;0 ð0Þ and d NIE 1;0 ð1Þ arising from unobserved mediator-outcome confounding are of similar sizes but opposite signs, i.e. cancel each other out. The delta method SEs appeared to target the empirical SEs ( Fig. 4a) but the large bias resulted in very poor coverage of the 95% CIs (Fig. 5a).
Bias over the 2000 replicates for scenarios a-e when using the true valueq ¼ 0:5 for estimation of effects is illustrated in Fig. 3b. The bias is generally small, especially for the larger sample sizes, although a slightly larger bias was observed for d NDE 1;0 ð0;q ¼ 0:5Þ in scenario c with a sample size of 500. The relative % error in delta method SE shown in Fig. 4b indicates that the delta method SEs generally Fig. 4 Relative % error for simulations with mediator-outcome confounding. Relative % error in delta method standard errors compared to empirical standard errors based on 2000 replicates for effects estimated using A:q ¼ 0 and B:q ¼ 0:5. The dotted vertical lines indicate 0% error. Black dots indicate relative % error in delta method SE for d NDE 1;0 ð0;qÞ and gray dots relative % error for d NIE 1;0 ð1;qÞ. Parentheses represent Monte Carlo 95% CIs appear to target empirical SEs. There is a tendency for a slight overestimation by the delta method SE of the d NIE 1;0 ð1;q ¼ 0:5Þ for sample sizes n obs ¼ 500; 5000 in scenario c, truncated outcome. Values of the empirical and delta method SEs are found in Tables S1 and S2 of Online Resource 1. Looking at the empirical coverage of 95% CIs in all scenarios (Fig. 5b) these are generally close to the nominal level.
Results for scenarios a, d and e with unobserved exposure-mediator confounding are found in Figs. 8, 9 and 10 in Appendix E and Tables S3 and S4 in Online Resource 1. Here the sensitivity parameter is the correlation between error terms in the exposure and mediator models,q zm . Corresponding results for unobserved exposure-outcome confounding are found in Figs. 11, 12 and 13 in Appendix E and Tables S5 and S6 in Online Resource 1. Here the sensitivity parameter is the correlation between error terms in the exposure and outcome models,q zy . We see similar results as in the simulations with unobserved mediator-outcome confounding, with small bias when using the correct correlation value (Figs. 8b and 11b) and  (Figs. 10b  and 13b). A notable difference is that unobserved exposure-mediator confounding leads to large bias for d NIE 1;0 ð1;q zm ¼ 0Þ but small bias for d NDE 1;0 ð0;q zm ¼ 0Þ and conversely unobserved exposure-outcome confounding leads to large bias for d NDE 1;0 ð0;q zy ¼ 0Þ but small bias for d NIE 1;0 ð1;q zy ¼ 0Þ.

Empirical example
To illustrate the method we use the publicly available data set UPBdata from the R package medflex (Steen et al. 2020). The data were used to illustrate functions in the medflex package in Steen et al. (2017)  Following the example in Steen et al. (2017) we look at the relationship between attachment style towards the ex-partner prior to the breakup and unwanted pursuit behaviors (UPBs) towards the ex-partner after the breakup and the extent to which this is mediated by level of emotional distress experienced during the breakup. A binary exposure is used, indicating whether or not the individual's self-reported anxious attachment level was higher than the sample mean. The outcome is whether or not the individual reported that they had displayed UPBs towards their ex-partner after the breakup. The mediator is standardized self-reported experienced level of negative affectivity (emotional distress) during the breakup, a continuous variable. We adjust for age, highest attained education level (high, intermediate, low) and gender (male, female). The hypothesized relationships between the variables are illustrated in Fig. 6. All analyses are performed using the sensmediation package (Lindmark 2019), code found at https://github.com/anitalindmark/ Sensitivity_analysis.
We begin with estimation of the natural direct and indirect effects assuming no unmeasured confounding and then proceed with sensitivity analyses to the three types of unmeasured confounding. Since we have a continuous mediator and a binary outcome the analyses are based on models (3) and (6) and the corresponding estimators from Table 1. In this example we investigate effect modification (moderation) by gender. To this end we include an interaction between the exposure and gender in the model for the mediator (Table S7 of Online Resource 1) as well as interactions between gender and both exposure and mediator in the outcome model (Table S8 of Online Resource 1). An interaction term between exposure and mediator is also included in the outcome model, as recommended by VanderWeele (2015) to fully capture the dynamics of mediation.
Estimated effects averaged over all observed confounders (marginal effects) as well as conditional on male and female gender are presented in Table 4. Looking at the estimated total effects we see that anxious attachment increases the risk of UPBs both marginally and conditional on gender, with a larger effect for men than for women. Over half of this total effect is an indirect effect of anxious attachment on UPBs operating through negative affectivity, with a slightly larger proportion for males than for females. To gauge the effect of possible unobserved confounding on the results we perform sensitivity analyses. Here we choose to focus on the natural indirect effect, which was statistically significant in the original analysis. We present results for all three types of confounding, with sensitivity parameters ranging from -0.9 to 0.9 in increments of 0.1. Plots of point estimates of the marginal natural indirect effect with corresponding CIs over the range of the sensitivity parameters are presented in Fig. 7. For both mediator-outcome confounding (Fig. 7a) and exposure-mediator confounding (Fig. 7b) the overall pattern is that the natural indirect effect decreases over the range of the sensitivity parameter. If additional adjustment were made for a confounder inducing an error term correlation (q orq zm ) of 0.3 or higher the CIs of the effect would include 0 and additional adjustment for a confounder inducing aq of at least 0.6 or aq zm of at least 0.5 would lead to CIs entirely below 0. The natural indirect effect is not sensitive to unobserved exposure-outcome confounding (Fig. 7c). Note that as the exposure and outcome are both binary the sensitivity analyses to exposure-outcome confounding were performed using methods presented in Lindmark et al. (2018).
The results in Fig. 7 are summarized in Table 5 which also shows the 95% UIs over the range of the sensitivity parameter. Corresponding results for men and women are also shown, indicating similar results as those seen for the marginal effect. For exposure-outcome confounding the lower bounds of the 95% UIs over the range of the sensitivity parameter all lie above 0, indicating that the effects are not sensitive to unobserved exposure-outcome confounding.

Discussion
In this paper we have extended results from Lindmark et al. (2018) to provide methods for sensitivity analysis of unobserved confounding in mediation analysis for combinations of continuous and binary mediators and outcomes, as well as for censored or truncated outcomes. Where previous methods focus exclusively on mediator-outcome confounding (VanderWeele 2010; Imai et al. 2010a;Hafeman 2011;le Cessie 2016), this approach is flexible due to the ability to take into account not only mediator-outcome confounding but also exposure-mediator and exposureoutcome confounding. The latter two are of particular importance in observational studies, where the exposure has not been randomized, due to the difficulty in guaranteeing that all relevant confounders have been adjusted for. It also has the advantage that sensitivity analyses can be performed also in the presence of interactions involving the exposure, mediator and covariates, allowing more complicated models and facilitating sensitivity analyses also when the interest lies in exploration of effect modification. We performed simulations that showed that the method targets the true effects when the error term correlation induced by unobserved confounding is taken into account in the estimation. Generally this illustrates that the method does indeed capture the effect that would have been observed under a given level of correlation. In reality this correlation will be unknown to the researcher and therefore further simulation studies investigating, e.g. the performance of UIs based on a range of correlation levels is of interest.
The method has some limitations that should be subjects for future development. One such limitation is the reliance on the specification of parametric regression models with distributional assumptions on the error terms. The results may therefore be sensitive to model misspecification and the nature of this sensitivity should be subject to further study. On the other hand, since the method allows inclusion of interactions involving the exposure, mediator and covariates, rich parametric models can be specified which can reduce the risk of model misspecification bias. This under the condition that the data allow such a specification and are large enough to lessen the impact of the increase in variance. In any case, further developments utilizing either semi-parametric techniques (Tchetgen Tchetgen and Shpitser 2012; Huber 2014) or retaining parametric regression models but relaxing the multivariate normality assumption of the error terms upon which the method introduced here relies are warranted. The issue of model misspecification is even more important for the proposed methods for censoring or truncation since maximum likelihood estimators for regression parameters when the outcome is censored or truncated have been shown to be sensitive to violations of distributional assumptions (Vijverberg 1987). Semiparametric estimators that impose fewer assumptions on the error term have been developed both for censoring, e.g. Powell (1986), and truncation, e.g. Powell (1986); Lee (1993); Laitila (2001). Further research into the usefulness of such models in the context of mediation is of interest.
For the cases with a censored/truncated outcome and a continuous mediator we have presented results for sensitivity analyses to mediator-outcome confounding only. Since we assume that only the outcome is censored, not the exposure or mediator, a sensitivity analysis to exposure-mediator confounding is straightforward and can be performed using the methods presented in Sect. 2.2.1. For a truncated outcome this would be more complicated since truncation of the outcome means that values will be missing for the exposure and mediator as well, which would need to be taken into account. For exposure-outcome confounding and a censored outcome, if the exposure is continuous and can be modeled with a linear regression model a sensitivity analysis could be performed by replacing the mediator model with the exposure model in the joint log-likelihood. The situation is again less straight-forward for truncation where the joint exposure-outcome distribution would need to be derived.
The methods presented in this paper evaluate sensitivity to each type of unobserved confounding separately, assuming that the other two kinds are not present. Extending the method to investigate sensitivity to all three types of confounding simultaneously is therefore of interest.
In this paper we present results for natural direct and indirect effects on the mean difference scale. Adapting the methods to other scales is of interest, in particular for cases with a binary outcome where the researchers may be interested in effects on the risk ratio or odds ratio scales (VanderWeele and Vansteelandt 2010;Valeri and VanderWeele 2013;Doretti et al. 2021).
Finally, it is important to note that the natural direct and indirect effects are not identified when the cross-world assumption introduced in Sect. 2.1 is violated. Such violations include the presence of mediator-outcome confounders that are affected by the exposure, regardless of whether these are observed or not. In such cases we either need to make additional parametric assumptions (Robins and Greenland 1992;Petersen et al. 2006;De Stavola et al. 2015) or use different effect definitions, e.g. so called interventional direct and indirect effects (see e.g. VanderWeele et al. 2014;Lok 2016).
To summarize, we have provided sensitivity analysis methods for unobserved confounding that are useful when performing parametric estimation of natural direct and indirect effects even when more complex models including interactions involving the exposure and/or mediator are used. With further developments these methods can be made even more flexible.
Appendix A Derivation of the joint distribution of M and Y under a binary probit mediator model (4) and a linear outcome model (5) To obtain the joint distribution of M i ; Y i given Z i ; X i , with M i following model (4) and Y i following model (5) we see that: where the final equality follows from ðg Ã i ;n i Þ 0 $N 0 0 ! ; 1 qr n qr n r 2 n ! , so that Using the same reasoning we have and finally the joint distribution Appendix B Gradients of the joint log-likelihood, censored outcome (7)  Appendix C Derivation of P(Y i > t), the probability of being included in the sample Assume that M i and Y i can be modeled with (8) and (9) and that only Y i [ t are observed. We have that P Y i [ t ð Þ¼1 À P Y i 6 t ð Þ¼1 À P n i 6 t À h y0 C 4i À Á , and P n y where in the second equality we replace M i with b y0 C 3i þ g y i , thus incorporating the regression parameters of the mediator model in the part of the distribution that takes truncation into account (similarly to the suggestion of Hausman and Wise (1977) for a two equation model). Since n y i þ h y 2 g y i $ i:i:d: Nð0; r 2 n y þ h y2 2 r 2 g y þ 2h y 2 qr n y r g y Þ we have that ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi r 2 n y þ h y2 2 r 2 g y þ 2h y 2 qr n y r g y q 0 B @ 1 C A: Appendix D Gradients of the joint log-likelihood, truncated outcome (12) n y þh y2 2 r 2 g y þ2h y 2 qr n y r g y q . Then, o' b y ; h y ; r g y ; r n y ; q o' b y ; h y ; r g y ; r n y ; q n y þ h y2 2 r 2 g y þ 2h y 2 qr n y r g y q 0 B @ 1 C A; o' b y ; h y ; r g y ; r n y ; q n y þ h y2 2 r 2 g y þ 2h y 2 qr n y r g y q þ t À h y0 n2 C 3i À h y 2 b y0 C 3i h y 2 r 2 g y þ qr n y r g y r 2 n y þ h y2 2 r 2 g y þ 2h y 2 qr n y r g y 3=2 ! ; o' b y ;h y ;r g y ;r n y ;q or g y ¼ À n r g y þ 1 ð1Àq 2 Þr 2 g y X i ( ðm i Àb y0 C 3i Þ 2 r g y À qðm i Àb y0 C 3i Þðy i Àh y0 C 4i Þ r n y !) t Àh y0 n2 C 3i Àh y 2 b y0 C 3i h y 2 h y 2 r g y þqr n y r 2 n y þh y2 2 r 2 g y þ2h y 2 qr n y r g y 3=2 ! ; o' b y ; h y ; r g y ; r n y ; q or n y ¼ À n r n y þ 1 ð1 À q 2 Þr 2 n y X i ( ðy i À h y0 C 4i Þ 2 r n y À qðm i À b y0 C 3i Þðy i À h y0 C 4i Þ r g y ) r n y þ h y 2 qr g y r 2 n y þ h y2 2 r 2 g y þ 2h y 2 qr n y r g y 3=2

!
: Appendix E Simulation results exposure-mediator and exposure outcome confounding  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain