8  Sensitivity analysis for mediation analysis

This session is inspired on Tyler VanderWeele’s book: Explanation in Causal Inference: Methods for Mediation and Interaction 1st Edition
Learning outcomes

After this session you will be able to:

  • Explain the assumptions for conducting sensitivity analysis for unmeasured confounding

  • Describe sensitivity analysis techniques to address unmeasured confounding

  • Perform sensitivity analysis using the CMAverse package and interpret the results

Unmeasured or uncontrolled confounding is a common problem in observational studies. This is a challenge to observational research even in the analysis of total effects.

When we are interested in direct and indirect effects, the assumptions about confounding that are needed to identify these effects are even stronger than for total effects.

We might be worried that these assumptions are violated and that our estimates are biased.

Sensitivity analysis techniques can help assess HOW ROBUST results are to violations in the assumptions being made.

These techniques assess the extent to which an unmeasured variable (or variables) would have to affect both the exposure and the outcome in order for the observed associations between the two to be attributable solely to confounding rather than a causal effect of the exposure on the outcome.

It can also be useful in assessing a plausible range of values for the causal effect of the exposure on the outcome corresponding to a plausible range of assumptions concerning the relationship between the unmeasured confounder and the exposure and outcome.

8.1 Sensitivity analysis for unmeasured confounding for total effects

Consider the following figure in which U represents an unmeasured confounder, C measured covariables, A the exposure and Y the outcome.

The basic idea of sensitivity analysis is to specify parameters corresponding to the relationships between U and Y and between U and A and from these, along with the observed data, to obtain “corrected” effect estimates corresponding to what would have been obtained had control been made for U and not only C.

The results essentially compare:

What we obtain adjusting only for measured covariables C with what we would have obtained had it been possible to adjust for measured covariables C and unmeasured covariable(s) U.

If it is thought that adjusting for C and U together would suffice to control for confounding, then we may also interpret the results as comparing the effect estimate that is obtained adjusting only for measured covariables C versus the true causal effect.

8.1.1 Continuous outcomes

Suppose then we have obtained an estimate of the effect of the exposure A on the outcome Y conditional on measured covariables C using regression analysis.

We will define the bias factor Badd(c) on the additive scale as the difference between the expected differences in outcomes comparing A = a and A = a* conditional on covariables C = c and what we would have obtained had we been able to adjust for U as well.

If the exposure is binary, then we simply have a = 1 and a* = 0.

A simple approach to sensitivity analysis is possible if we assume that (A8.1.1) U is binary and (A8.1.2) that the effect of U (on the additive scale) is the same for those with exposure level A = a and exposure level A = a*(no U × A interaction).

If these assumptions hold, let γ be the effect of U on Y conditional on A and C, that is:

\(γ = E(Y|a,c,U = 1)\) − \(E(Y|a,c,U = 0)\)

Note that by assumption (A8.1.2),

\(γ = E(Y|a,c,U = 1)\) − \(E(Y|a,c,U = 0)\)

is the same for both levels of the exposure of interest.

Note also that γ is the effect of U on Y already having adjusted for C; that is, in some sense the effect of U on Y not through C

Now let δ denote the difference in the prevalence of the unmeasured confounder U for those with A=a versus those with A = a*, that is:

\(δ = P(U = 1|a,c)\) − \(P(U = 1|a^*,c)\)

Under assumptions (A8.1.1) and (A8.1.2), the bias factor is simply given by the product of these two sensitivity analysis parameters:

\(B_{add}(c) = γδ\) (8.1)

Thus to calculate the bias factor we only need to specify the effect of U on Y and the prevalence difference of U between the two exposure groups and then take the product of these two parameters.

Once we have calculated the bias term Badd(c), we can simply estimate our causal effect conditional on C and then subtract the bias factor to get the “corrected estimate” - that is, what we would have obtained if we had controlled for C and U.

Under these simplifying assumptions (A8.1.1) and (A8.1.2), we can also get adjusted confidence intervals by simply subtracting γδ from both limits of the estimated confidence intervals.

We may not believe any particular specification of the parameters γ and δ, but we could vary these parameters (based on expert knowledge or previous reported estimates of the associations of the C and Y) over a range of plausible values to obtain what were thought to be a plausible range of corrected estimates.

Using this technique, we could also examine how substantial the confounding would have to be to explain away an effect (we could do this for the estimate and confidence interval).

8.1.2 Continuous outcome with different sensitivity analysis parameters for different covariable values

Suppose now that instead of focusing on effects conditional on a particular covariable value C = c or specifying the sensitivity analysis parameters γ and δ to be the same for each covariable C, we were interested in the overall marginal effect averaged over the covariables and we wanted to specify different sensitivity analysis parameters for different covariable levels.

Suppose then for each level of the covariables of interest C = c we specified a value for the effect of U on Y

\(γ(c) = E(Y|a, c,U = 1) − E(Y|a, c,U = 0)\)

and also a value for the prevalence difference of U between those with exposure status A = a and A = a* and covariables C = c

\(δ(c) = P(U = 1|a, c)−P(U = 1|a^*, c)\)

We could then obtain an overall bias factor, Badd, by taking the product of the bias factors in each strata of C and then averaging these over C, weighting each strata of C according to what proportion of the sample was in that strata. The overall bias factor is then

\(B_{add}=\sum~c~\{γ(c)δ(c)\}P(C=c)\)

We could then subtract this overall bias factor from our estimate adjusted only for C to obtain a corrected estimate.

In this case, however, we can no longer simply subtract the bias factor from both limits of the confidence intervals because this does not take into account the variability in our estimates of the proportion of the sample in each strata of the covariables P(C = c).

Corrected confidence intervals could instead be obtained by bootstrapping.

8.1.3 Binary outcomes

Sensitivity analysis approach for odds ratios or risk ratio estimates of the effect of the exposure A on outcome Y, conditional on the covariables C obtained from a logistic regression model.

The bias factor is defined as Bmult(c) on the multiplicative scale as the ratio of:

1- The risk ratio (or odds ratio, with a rare outcome) comparing A = a and A = a* conditional on covariables C= c and

2- What we would have obtained as the risk ratio (or odds ratio) had we been able to condition on both C and U.

We now make the simplifying assumptions that (A8.1.3) U is binary and that (A8.1.2b) the effect of U (on the risk ratio scale) is the same for those with exposure level A = a and exposure level A = a*(no U × A interaction on the risk ratio scale).

If these assumptions hold, we will let γ be the effect of U on Y conditional on A and C on the risk ratio scale, that is:

\(γ = \frac{ P(Y = 1|a, c,U = 1)}{P(Y = 1|a, c,U = 0)}\)

By assumption (A8.1.2b)

\(γ = \frac{P(Y=1|a,c,U=1)}{P(Y=1|a,c,U=0)}\)

is the same for both levels of the exposure.

This is the effect of U on Y adjusted for C; this is the effect of U on Y not through C.

Under assumptions (A8.1.1) and (A8.1.2b), the bias factor on the multiplicative scale is given by:

\(B_{mult}(c) = \frac{1+(γ −1)P(U = 1|a,c)}{1+(γ −1)P(U = 1|a∗, c)}\) (8.2)

We can use the bias formula by specifying the effect of U on Y on the risk ratio scale and the prevalence of U among those with exposure levels A = a and A = a*.

Once we have calculated the bias term Bmult(c), we can estimate our risk ratio controlling only for C (if the outcome is rare, fit a logistic regression) and we divide our estimate by Bmult(c) to get the corrected estimate for risk ratio—that is, what we would have obtained if we had adjusted for U as well.

Under the simplifying assumptions of (A8.1.1) and (A8.1.2b), we can also obtain corrected confidence intervals by dividing both limits of the confidence interval by Bmult(c).

Note that to use the bias factor in (8.2), we must specify the prevalence of the unmeasured confounder in both exposure groups P(U= 1|a, c) and P(U = 1|a*, c), not just the difference between these two prevalences as in (8.1) for outcomes on the additive scale.

At a minimum, it may be useful to present:

1- the sensitivity analysis parameters that would suffice to completely explain away an effect and also

2- the sensitivity analysis parameters that would be required to shift the confidence interval to just include the null.

8.2 Sensitivity analysis for controled direct effects

8.2.1 Continuous outcomes

Assume that controlling for (C,U) would suffice to control for exposure–outcome and mediator–outcome confounding but that no data are available on U and that U confounds the mediator–outcome relationship.

If we have not adjusted for U, then our estimates controlling only for C will be biased.

We will consider estimating the controlled direct effect, CDE(m), with the mediator fixed to m conditional on the covariables C = c.

Let \(B^{CDE}_{add}(m|c)\) denote the difference between:

1- the estimate of the CDE conditional on C

2- what would have been obtained had adjustment been made for U as well.

As with total effects, we will be able to use a simple formula for sensitivity analysis for CDE under some simplifying assumptions.

Suppose that (A8.1.1) U is binary and (A8.2.2b) the effect of U on Y on the additive scale, conditional on exposure, mediator, and covariables, (A,M,C), is the same for both exposure levels A = a and A = a*.

Let γm be the effect of U on Y conditional on A, C, and M = m, that is:

\(γm = E(Y|a,c,m,U = 1)−E(Y|a,c,m,U = 0)\)

Note that by assumption (A8.2.2b) is the same for both levels of the exposure.

Let δm be the difference in the prevalence of the unmeasured confounder for those with A = a versus those with A = a* conditional on M = m and C = c, that is:

\(δm = P(U = 1|a,m,c) − P(U = 1|a^*,m,c)\)

Under assumptions (A8.1.1) and (A8.2.2b), the bias factor is simply given by the product of these two sensitivity-analysis parameters (VanderWeele, 2010a):

\(B^{CDE}_{add}(m|c) = δmγm\) 8.3

This formula states that under assumptions (A8.1.1) and (A8.2.2b) the bias factor \(B^{CDE}_{add}(m|c)\) for the CDE(m) is simply given by the product δmγm.

Under these simplifying assumptions, this gives rise to a particularly simple sensitivity analysis technique for assessing the sensitivity of estimates of a CDE to an unmeasured mediator–outcome confounder.

We can hypothesize a binary unmeasured mediator–outcome confounding variable U such that the difference in expected outcome Y comparing U = 1 and U = 0 is γm across strata of A conditional on M = m, C = c and such that the difference in the prevalence of U, comparing exposure levels a and a* (comparing the exposed and unexposed), is δm conditional on M =m, C = c.

For such an unmeasured mediator–outcome confounding variable, the bias of our estimate of the CDE controlling just for C is given simply by δmγm.

We can assess sensitivity to the presence of such an unmeasured confounding variable by varying γm (which is the direct effect of U on Y) and by varying δm, interpreted as the prevalence difference of U, comparing exposure levels a and a* conditional on M = m and C = c.

We can subtract the bias factor \(B^{CDE}_{add}(m|c)\) from the observed estimate to obtain a corrected estimate of the effect (what we would have obtained had it been possible to adjust for U as well).

Under the simplifying assumptions (A8.1.1) and (A8.2.2b), we could also subtract this bias factor from both limits of a confidence interval to obtain a corrected confidence interval.

Note that the CDE(m), may vary with m, and so for different values of m we will likely want to consider different specifications of the values δm and γm in the sensitivity analysis.

If there is no interaction between the effects of A and M on Y, then this simple sensitivity analysis technique based on using formula above will also be applicable to natural direct effects as well.

8.2.2 Binary outcomes

We will consider estimating the controlled direct effect odds ratio as \(OR^{CDE}(m)\), with the mediator fixed at level m, conditional on the covariates C = c.

This approach will assume a rare outcome but can also be used for risk ratios with a common outcome. Let \(B^{CDE}_{mult}(m|c)\) denote the ratio of

1- The estimate of the controlled direct effect conditional on C

2- What would have been obtained had adjustment been made for U as well.

Suppose that (A8.1.1) U is binary and that (A8.1.2d) the effect of U on Y on the ratio scale, conditional on exposure, mediator, and covariables (A,M,C), is the same for both exposure levels A = a and a*.

Let γm be the effect of U on Y conditional on A, C, and M = m, that is:

\(γm = \frac{P(Y = 1|a, c,m,U = 1)}{P(Y = 1|a, c,m,U = 0)}\)

Note that by (A8.1.2), γm is the same for both levels of the exposure of interest.

Under assumptions (A8.1.1) and (A8.1.2d), the bias factor on the multiplicative scale is given by:

\(B^{CDE}_{mult}(m|c) = \frac{1+(γm −1)P(U= 1|a,m, c)}{1+(γm−1)P(U = 1|a∗,m, c)}\) (8.4)

Once we have calculated the bias term \(B^{CDE}_{mult}(m|c)\), we can estimate the CDE risk ratio controlling only for C (if the outcome is rare), we fit a logistic regression) and we divide our estimate and confidence intervals by the bias factor \(B^{CDE}_{mult}(m|c)\) to get the corrected estimate for CDE risk ratio and its confidence interval—that is, what we would have obtained if we had adjusted for U a well.

We have to specify the two prevalences of U, namely \(P(U = 1|a,m, c)\) and \(P(U = 1|a∗,m, c)\), in the different exposure groups conditional on M and C.

As with CDE on an additive scale, the issue of conditioning on M in the interpretation of these prevalences is important

8.3 Sensitivity analysis for natural direct and indirect effects

8.3.1 Sensitivity analysis for natural direct and indirect effects in the abscence of exposure-mediator interaction

A simple setting in which we can employ sensitivity analysis for natural direct (NDE) and indirect effects (NIE) is when the NDE and CDE coincide.

This occurs when the four confounding assumptions are satisfied and there is no exposure–mediator interaction in the statistical model.

If this is the case, we can use the sensitivity analysis techniques for CDE.

If we assume an unmeasured mediator–outcome confounder U as in the next figure, we can use these same techniques and the same parameters to do sensitivity analysis for NDE as well.

For NIE, we can use of the decomposition property of the total effect.

Let’s assume that we have unmeasured mediator–outcome, even if this present and our estimates of the NDE and NIE are biased, the NDE and NIE themselves will still combine to the correct total effect.

Because a mediator–outcome confounder does not confound the exposure-outcome relationship, we can still obtain valid estimates of the total effect.

And, it turns out that the combination of the DE and IE do constitute a consistent estimator of the total effect, even though the DE and IE estimators will themselves be biased for the true NDE and NIE.

Knowing that the DE and IE estimates combine to a valid estimate of the total effect then allows us to employ the sensitivity analysis techniques for CDE for NIE as well.

To do so, we use the negation (on the additive scale) or the inverse (on the multiplicative ratio scale) of the bias formulas used for CDE (and NDE). Thus on the additive scale, for a continuous outcome, our bias factor for the NDE would simply be:

−δmγm

and we could subtract this from the estimate and both limits of the confidence interval to obtain a corrected estimate and confidence interval for the NIE.

For a binary outcome, on the odds ratio scale with rare outcome or risk ratio scale with common outcome, our bias factor for the NIE would be the inverse of that in (8.4):

\(\frac{1+(γm−1)P(U=1|a∗,m,c)}{1+(γm−1)P(U=1|a,m,c)}\)

and we could divide our NIE estimates and its confidence interval by this bias factor to obtain a corrected estimate and confidence interval.

8.4 Example using the NHANES data

Let’s first load the nhanes dataset:

# first load the dataset
nhanes <- read.csv(here::here("data/nhanes_dataset.csv"))

nhanes <- nhanes %>%
  select(
    id = seqn,
    w1 = age,
    w2 = gender,
    w3 = education_clean,
    w4 = smoke, 
    a = total_redmeat,            #this is the exposure
    m = magic_biomarker,          #this is the mediator
    y = blood_glucose) %>%        #this is the outcome
    na.omit()

We will use the red meat, inflammation and blood glucose example

res_rb_confounders <- cmest(
  data = nhanes, model = "rb", outcome = "y", exposure = "a",
  mediator = "m", basec = c("w1", "w2", "w3"), EMint = TRUE,
  mreg = list("linear"), yreg = "linear",
  astar = 0, a = 1, mval = list(2.5),
  estimation = "paramfunc", inference = "delta"
)

summary(res_rb_confounders)
#> Causal Mediation Analysis
#> 
#> # Outcome regression:
#> 
#> Call:
#> glm(formula = y ~ a + m + a * m + w1 + w2 + w3, family = gaussian(), 
#>     data = getCall(x$reg.output$yreg)$data, weights = getCall(x$reg.output$yreg)$weights)
#> 
#> Coefficients:
#>                      Estimate Std. Error  t value Pr(>|t|)    
#> (Intercept)         17.793668   0.609882   29.176   <2e-16 ***
#> a                   52.238325   0.758400   68.880   <2e-16 ***
#> m                    0.639863   0.052164   12.266   <2e-16 ***
#> w1                   1.499841   0.001149 1305.718   <2e-16 ***
#> w2Male               0.019112   0.045521    0.420    0.675    
#> w3College and above  0.015804   0.048412    0.326    0.744    
#> a:m                  0.057466   0.055168    1.042    0.298    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 4.028118)
#> 
#>     Null deviance: 7017320  on 8568  degrees of freedom
#> Residual deviance:   34489  on 8562  degrees of freedom
#> AIC: 36266
#> 
#> Number of Fisher Scoring iterations: 2
#> 
#> 
#> # Mediator regressions: 
#> 
#> Call:
#> glm(formula = m ~ a + w1 + w2 + w3, family = gaussian(), data = getCall(x$reg.output$mreg[[1L]])$data, 
#>     weights = getCall(x$reg.output$mreg[[1L]])$weights)
#> 
#> Coefficients:
#>                       Estimate Std. Error  t value Pr(>|t|)    
#> (Intercept)         -1.943e+01  1.339e-01 -145.052   <2e-16 ***
#> a                    2.585e+01  1.491e-01  173.390   <2e-16 ***
#> w1                  -4.449e-04  5.716e-04   -0.778    0.436    
#> w2Male               3.090e-02  2.258e-02    1.369    0.171    
#> w3College and above  2.017e-02  2.416e-02    0.835    0.404    
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for gaussian family taken to be 1.003419)
#> 
#>     Null deviance: 42877.0  on 8568  degrees of freedom
#> Residual deviance:  8593.3  on 8564  degrees of freedom
#> AIC: 24354
#> 
#> Number of Fisher Scoring iterations: 2
#> 
#> 
#> # Effect decomposition on the mean difference scale via the regression-based approach
#>  
#> Closed-form parameter function estimation with 
#>  delta method standard errors, confidence intervals and p-values 
#>  
#>               Estimate Std.error   95% CIL 95% CIU  P.val    
#> cde          52.381990  0.692505 51.024705  53.739 <2e-16 ***
#> pnde         51.121905  1.617002 47.952640  54.291 <2e-16 ***
#> tnde         52.607442  0.637358 51.358243  53.857 <2e-16 ***
#> pnie         16.540891  1.351851 13.891312  19.190 <2e-16 ***
#> tnie         18.026428  0.603132 16.844311  19.209 <2e-16 ***
#> te           69.148334  1.325762 66.549888  71.747 <2e-16 ***
#> intref       -1.260085  1.209714 -3.631081   1.111  0.298    
#> intmed        1.485537  1.426153 -1.309671   4.281  0.298    
#> cde(prop)     0.757531  0.013025  0.732002   0.783 <2e-16 ***
#> intref(prop) -0.018223  0.017834 -0.053177   0.017  0.307    
#> intmed(prop)  0.021483  0.021025 -0.019724   0.063  0.307    
#> pnie(prop)    0.239209  0.015627  0.208581   0.270 <2e-16 ***
#> pm            0.260692  0.011312  0.238522   0.283 <2e-16 ***
#> int           0.003260  0.003191 -0.002993   0.010  0.307    
#> pe            0.242469  0.013025  0.216941   0.268 <2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (cde: controlled direct effect; pnde: pure natural direct effect; tnde: total natural direct effect; pnie: pure natural indirect effect; tnie: total natural indirect effect; te: total effect; intref: reference interaction; intmed: mediated interaction; cde(prop): proportion cde; intref(prop): proportion intref; intmed(prop): proportion intmed; pnie(prop): proportion pnie; pm: overall proportion mediated; int: overall proportion attributable to interaction; pe: overall proportion eliminated)
#> 
#> Relevant variable values: 
#> $a
#> [1] 1
#> 
#> $astar
#> [1] 0
#> 
#> $mval
#> $mval[[1]]
#> [1] 2.5
#> 
#> 
#> $basecval
#> $basecval[[1]]
#> [1] 45.5381
#> 
#> $basecval[[2]]
#> [1] 0.4969075
#> 
#> $basecval[[3]]
#> [1] 0.2800794

To perform sensitivity analysis for unmeasured confounding we simply add another line of code.

uc_sens <- cmsens(
  object = res_rb_confounders,
  sens = "uc"
)

uc_sens
#> Sensitivity Analysis For Unmeasured Confounding 
#> 
#> Evalues on the risk or rate ratio scale: 
#>         estRR  lowerRR  upperRR Evalue.estRR Evalue.lowerRR Evalue.upperRR
#> cde  5.288983 5.066011 5.521768    10.051793       9.604564             NA
#> pnde 5.081255 4.595072 5.618878     9.635142       8.659505             NA
#> tnde 5.327035 5.119991 5.542452    10.128105       9.712845             NA
#> pnie 1.692086 1.555629 1.840512     2.774245       2.485335             NA
#> tnie 1.773932 1.708619 1.841742     2.945642       2.808964             NA
#> te   9.013801 8.300349 9.788578    17.512908      16.084657             NA
Note

  • The E-value is the minimum strength of association, on the risk ratio scale, that an unmeasured confounder would need to have with both the treatment and the outcome to fully explain away a specific treatment–outcome association, conditional on the measured covariates (1).

  • A large E-value implies that considerable unmeasured confounding would be needed to explain away an effect estimate.

  • A small E-value implies little unmeasured confounding would be needed to explain away an effect estimate.

8.5 References