Sensitivity analysis for causal effects with generalized linear models

1 Introduction

Bias due to residual confounding is one of the main threats to the validity of causal conclusions in observational studies. A common way to deal with confounding bias is to carry out a sensitivity analysis, where the analyst is required to speculate about certain sensitivity parameters that quantify the magnitude of residual confounding. These parameters are used, together with the observed data, to reason about the magnitude of the causal effect.

The literature on such sensitivity analyses can be roughly classified along two dimensions. The first of these refers to the type of sensitivity parameters used. In particular, one may start by considering an unmeasured (set of) confounder(s) U that is sufficient for confounding control, and parametrize its effect on the exposure and the outcome. The resulting parameters are then used as sensitivity parameters. Alternatively, one may use sensitivity parameters that quantify the deviation from (conditional, given measured confounders) exchangeability. Formally, exchangeability is defined as statistical independence between the exposure and the potential outcomes that define the target causal effect [1,2]. Residual confounding leads to nonexchangeability; the stronger the confounding, the larger the deviation from exchangeability. This approach has the advantage of not requiring the analyst to explicitly consider a set of unmeasured confounders, but the resulting sensitivity analysis parameters may be more abstract and difficult to interpret. The second dimension refers to the type of inference that the analysis produces about the target causal effect. In particular, the sensitivity analysis may allow point estimation of the target causal effect, given the sensitivity parameters, or it may only produce bounds for it. How published sensitivity analysis methods fall into the four categories defined by the cross-tabulation of these two dimensions is illustrated in Table 1.

In the analyses listed in the top-left quadrant of Table 1, the analyst parametrizes the effect of a pre-specified set of unmeasured confounders that they deem sufficient for confounding control. The estimation problem is formulated so that the target causal effect can be point estimated given assumed values of the sensitivity parameters. Arguably, this is the most conceptually straight-forward way of carrying out a sensitivity analysis, and perhaps for this reason, it was also the earliest to appear in the literature [3]. Later authors considered various scenarios and model assumptions for the exposure, the outcome and the unmeasured confounders. In contrast, in the analysis by Rosenbaum [4] (top-right quadrant of Table 1), the analyst does not assume fixed values of the sensitivity parameters; rather, they specify a threshold for the maximal effect that U may have on the exposure, when contrasting all possible levels of U . Rosenbaum [4] showed how this threshold value translates into bounds for the target causal effect. Yadlowsky et al. [5] developed several sophisticated estimators of these bounds, based on semi-parametric theory and machine learning methods. Ding and Vander Weele [6] and Vander Weele and Ding [7] developed a similar sensitivity analysis and bounds, but based on specified thresholds for both the maximal effect of U on the exposure, and the maximal effect of U on the outcome. They coined the term “ E -value” for the common threshold value of these effects required to “nullify” an observed exposure–outcome association.

The analyses by Cinelli and Hazlett [8] and Chernozhukov et al. [9] stand out from other analyses on the top row of Table 1 in an important aspect. Although the other analyses typically define the sensitivity parameters as relatively standard regression coefficients (e.g., mean differences or risk ratios), Cinelli and Hazlett [8] and Chernozhukov et al. [9] defined them in terms of the variation in the exposure and the outcome, respectively, that is explained by U , e.g., as partial R 2 coefficients. Although these measures may be perceived as somewhat more abstract than standard regression coefficients, they may also enable a more general analysis. For instance, analyses by Cinelli and Hazlett [8] and Chernozhukov et al. [9] allow for arbitrary exposures, whereas many of the other analyses require the exposure to be either binary or normally distributed.

For binary exposures, Brumback et al. [10] (bottom-left quadrant of Table 1) considered a “confounding function” that describes how the mean potential outcome depends on the exposure. They developed an analysis with sensitivity parameters that index this confounding function and are defined as differences in mean potential outcomes. Chiba [11] extended this analysis to allow for sensitivity parameters defined as ratios of means of potential outcomes. Both Brumback et al. [10] and Chiba [11] focused on causal effects defined through marginal structural models and estimated with inverse probability weighting, but Blackwell [12] showed how his sensitivity analyses can be used in combination with other models and estimation methods as well. Ciocănea-Teodorescu et al. [13] developed a related sensitivity analysis, but using a parsimonious set of sensitivity parameters defined marginally over the measured confounders, instead of conditionally as Brumback et al. [10] and Chiba [11]. Franks et al. [14] and Scharfstein et al. [15] used a somewhat different way of parametrizing the deviation from exchangeability. Instead of the confounding function, they considered a “tilting function” that describes how the probability of exposure depends on the potential outcomes. They developed an analysis with sensitivity parameters that index this tilting function and are defined as odds ratios. For both the analyses using the confounding function and the tilting function, the target causal effect can be point estimated given assumed values of their respective sensitivity parameters.

Finally, all cited authors in the bottom-right quadrant of Table 1 used the tilting function (or some transformation thereof) like Franks et al. [14] and Scharfstein et al. [15]. In the spirit of Rosenbaum [4], they developed sensitivity analyses where the analyst specifies a threshold for the tilting function, which translates into bounds for the target causal effect. They developed various parametric and semi-parametric estimators for these bounds.

In this article, we build on the work by Brumback et al. [10] and Chiba [11]. Like these authors, we use sensitivity parameters that quantify the deviation from exchangeability in terms of the confounding function. However, in our generalized sensitivity analysis, we allow for an arbitrary exposure, and sensitivity parameters defined as any contrast between counterfactual means, e.g., odds ratios. We further suggest a simple way to perform the sensitivity analysis in practice, which can be easily implemented with standard software for generalized linear models. Both Brumback et al. [10] and Chiba [11] considered scenarios where the exposure and outcome are measured at a single occasion, and scenarios where they are measured repeatedly during follow-up; we restrict attention to the former.

This article is organized as follows. We first introduce notation, definitions and assumptions that we will use throughout the article in Section 2. We next briefly review the sensitivity analyses previously proposed by Brumback et al. [10] and Chiba [11] in Section 3. Then, we present our generalized sensitivity analysis in Section 4 and an estimation method based on generalized linear models in Section 5. To help practitioners select plausible values for the sensitivity parameters, we propose a method of calibrating them against measured confounders in Section 6. We present a small simulation study, which empirically confirms the theoretical properties of our estimation and calibration methods in Section 7. Finally, we demonstrate our sensitivity analysis by re-analyzing data from a cohort study of behavior patterns and coronary heart disease in Section 8. These data are publicly available, which enables the reader to replicate and elaborate on our analysis.

2 Notation, definitions and assumptions

Let X and Y be the exposure and outcome of interest, and let Z be a set of measured confounders that we wish to control for in the analysis. We allow for these variables to be either continuous, discrete or (possibly multilevel) categorical. To simplify notation, we use p ( ⋅ ) for both probabilities and probability densities. We assume that data consist of n independent and identically distributed observations from p ( Z , X , Y ) . We let ( Z i , X i , Y i ) be the observation for subject i , for i = 1 , … , n , and suppress the index i when not needed.

Let Y ( x ) be the potential outcome if X were, possibly contrary to fact, set to level x [1,2,26]. We define the target causal effect as some contrast between the mean potential outcomes Y ( x ) for different values of x . For binary exposures such contrasts include, for instance, the conditional (on Z ) mean difference E { Y ( 1 ) ∣ Z } − E { Y ( 0 ) ∣ Z } or the marginal (over Z ) odds ratio [ p { Y ( 1 ) = 1 } / p { Y ( 1 ) = 0 } ] / [ p { Y ( 0 ) = 1 } / p { Y ( 0 ) = 0 } ] . This formulation of the target causal effect is very general and entirely nonparametric. In Section 5, we suggest an estimation method for the target effect based on parametric models.

We make the usual consistency and positivity assumptions. The consistency assumption states that the potential outcome Y ( x ) is equal to the factual outcome Y whenever the exposure X is factually equal to x :

The positivity assumption states that, for each Z = z with nonzero probability, there is a nonzero probability of X = x :

When there is no residual confounding, the potential outcome Y ( x ) is conditionally independent of the exposure X , given Z :

which is referred to as conditional exchangeability, given Z [1,2]. Under (1)–(3), both conditional and marginal causal effects are identifiable from the joint distribution p ( Z , X , Y ) . Specifically, E { Y ( x ) ∣ Z } is identifiable through

where the first equality follows from (2) and (3), and the second follows from (1). It follows from the law of total probability that E { Y ( x ) } is identifiable through

where the outer expectation is taken over the marginal distribution of Z . When there is residual confounding, the relation in (3) does not generally hold, and neither conditional nor marginal causal effects are generally identifiable.

3 Previously proposed sensitivity analyses

To assess the sensitivity of observed exposure–outcome associations to residual confounding, in settings where the exposure is binary, Brumback et al. [10] proposed to use the “confounding function”

with η being the identity link. Brumback et al. [10] further proposed to parametrize the confounding function with a (scalar or vector) sensitivity parameter α that measures the deviation from conditional exchangeability; the larger the magnitude of α , the stronger the deviation from (3). Conversely, under conditional exchangeability, we require that α = 0 . A simple example of such a parametrization is the function c ( x , Z ; α ) = α ( 2 x − 1 ) , for which E { Y ( 1 ) ∣ X = 1 , Z } − E { Y ( 1 ) ∣ X = 0 , Z } = − [ E { Y ( 0 ) ∣ X = 0 , Z } − E { Y ( 0 ) ∣ X = 1 , Z } ] = α . When α > 0 , both mean potential outcomes are higher for those factually exposed than for those factually nonexposed; thus, for this choice of confounding function, the residual confounding biases the observed exposure-outcome association upwards. Conversely, when α < 0 , the residual confounding biases the observed association downwards. A slightly more complex function allows for the deviation from exchangeability to depend on the levels of the measured confounders, e.g., by including an interaction term between ( 2 x − 1 ) and Z as c ( x , Z ; α ) = α 0 ( 2 x − 1 ) + α 1 ( 2 x − 1 ) Z . We refer readers to the study by Brumback et al. [10] for more examples and discussion of possible confounding functions.

Given c ( x , Z ; α ) , Brumback et al. [10] proposed to construct a “bias-corrected” outcome

for each subject i . They showed that this correction removes the confounding bias, in the sense that E ( Y i , α ∣ X i = x , Z i ) = E { Y ( x ) ∣ Z i } . This suggests that a sensitivity analysis can be carried out by first using subject matter knowledge to specify the confounding function c ( x , Z ; α ) and a set of plausible values for α ; then by using the conditional (given Z ) association between Y i , α and X i as an estimate of the target causal effect for each value of α . Chiba [11] showed that an analogous sensitivity analysis can be carried out by letting η in (5) be the log link and defining the bias-corrected outcome Y i , α as follows:

for which we again have that E ( Y i , α ∣ X i = x , Z i ) = E { Y ( x ) ∣ Z i } . Computing the bias-corrected outcomes proposed by Brumback et al. [10] and Chiba [11] requires the conditional distribution p ( X ∣ Z ) . This is usually not known but can be estimated from data, e.g., with a logistic regression model. Both Brumback et al. [10] and Chiba [11] focused on causal effects defined through marginal structural models and estimated with inverse probability weighting, but Blackwell [12] noted that the sensitivity analysis can also be carried out with other models and methods, such as ordinary outcome regression or matching.

4 Generalized sensitivity analysis

To accommodate arbitrary exposures and link functions, we propose to generalize the confounding function in (5) as follows:

where η is an arbitrary smooth link function, e.g., the identity, log or logit link. Like the simpler confounding function in (5), we propose to parametrize it with a (scalar or vector) sensitivity parameter α that measures the deviation from conditional exchangeability (3) and is equal to 0 when conditional exchangeability holds. In addition, the parametrization is such that c ( x , x , Z ; α ) = 0 for all ( x , Z ) when α = 0 . In contrast to the simpler confounding function in (5), the generalized confounding function in (8) allows for any measure of deviation that can be defined as a difference between transformed (by the link function η ) means, and it allows for the exposure to be (possibly multilevel) categorical or continuous.

A straightforward generalization of the simple confounding function c ( x , Z ; α ) = α ( 2 x − 1 ) for binary exposures is the function c ( x , x ′ , Z ; α ) = α ( x − x ′ ) . When α > 0 , this function assumes that the (transformed by η ) mean potential outcome increases linearly with the factual exposure level, and it reduces to c ( x , Z ; α ) = α ( 2 x − 1 ) when X is binary. For instance, suppose that Y is binary, η is the logit link and α = 0.1 . Then, the confounding function c ( x , x ′ , Z ; α ) = α ( x − x ′ ) implies that the log odds difference logit [ p { Y ( x ) = 1 ∣ X = x , Z } ] − logit [ p { Y ( x ) = 1 ∣ X = x ′ , Z } ] increases with 0.1 units for each unit increase in the exposure difference x − x ′ . More complex confounding functions can be modeled with, for instance, splines or polynomials. One may also let the deviation from conditional exchangeability depend on levels of the measured confounders, e.g., by including an interaction term between ( x − x ′ ) and Z as c ( x , Z ; α ) = α 0 ( x − x ′ ) + α 1 ( x − x ′ ) Z . When the exposure is categorical with ordered levels x 1 , … , x J , one may, for instance, parametrize the confounding function as c ( x j , x j ′ , Z ; α ) = α ( j − j ′ ) . When the exposure levels are not ordered, it may be difficult to defend such a parsimonious model for the confounding function, and a realistic model may need to include (at least) one parameter for each combination ( x j , x j ′ ) . As the dimension of the sensitivity parameter increases, it also becomes increasingly more difficult to suggest plausible values for it. Thus, we conjecture that, for numeric (e.g., binary and continuous) exposures, many analysts will in practice be content with the simple confounding function c ( x , x ′ , Z ; α ) = α ( x − x ′ ) , which only requires specification of a scalar sensitivity parameter α .

Given c ( x , x ′ , Z ; α ) , E { Y ( x ) ∣ Z } is identifiable through

where the first equality follows from the law of total probability, the second from positivity (2) and the definition of c ( x , X , Z ; α ) and the third from consistency (1). The outer expectation on the right-hand side of (9) is taken over the conditional distribution p ( X ∣ Z ) . The relation in (9) suggests that a sensitivity analysis can be conducted as shown in the studies by Brumback et al. [10] and Chiba [11], by first using subject matter knowledge to specify the confounding function c ( x , x ′ , Z ; α ) and a set of plausible parameter values for α , and then estimating the target causal effect for each value of α in this set.

Allowing for an arbitrary link function η is important, as it enables the analyst to define the confounding function on the same scale as the target causal effect, which in most situations would be perceived as natural. Furthermore, we show in Appendix A that, for binary outcomes, the logit link has an important advantage over the identity and log links, in that it makes the confounding function variation independent of the observed data distribution p ( Z , X , Y ) ; that is, all values of c ( x , x ′ , Z ; α ) are logically compatible with all values of p ( Z , X , Y ) . This has two important consequences [14,27]. First, it enables separation of the analysis into the modeling of p ( Z , X , Y ) and the choice of confounding function, thus effectively separating the assessment of uncertainty due to sampling variability and the uncertainty due to unmeasured confounding. Second, regardless of whether the value for α is determined by subject matter knowledge or guided by calibration against measured confounders (see Section 6), the analyst does not need to worry about potential conflicts between the chosen value for α and the observed data distribution. We show in Section 8 how this property of the logit link plays out with real data, and how the absence of this property for the log link makes inference more challenging.

5 Estimation and asymptotic inference

Suppose that the confounding function c ( x , x ′ , Z ; α ) and the sensitivity parameter α are fixed. If both X and Z are low dimensional (e.g., binary), then the identifying relation in (9) can be used to estimate E { Y ( x ) ∣ Z } nonparametrically by replacing the true mean E ( Y ∣ X = x , Z ) with the corresponding sample mean and taking the outer expectation over the conditional distribution of X given Z in the sample. To handle more general situations, we propose an estimation method based on parametric models.

We assume that the target causal effect can be expressed through the conditional (on Z ) causal model

or through the marginal (over Z ) causal model

In these models, η is an appropriate link function, e.g., the identity, log or logit link; ψ is a column vector of parameters and g is a regression function with the same dimension as ψ . In principle, we may allow for the link function to be different in the causal models (10) and (11) than in the confounding function (8), but we suspect that these will usually be the same in practice. A simple example is when g in model (10) is linear in x and Z : η [ E { Y ( x ) ∣ Z } ] = ψ 0 + ψ 1 x + ψ 2 Z . In this model, the target causal effect is the mean counterfactual contrast ψ 1 = η [ E { Y ( x + 1 ) ∣ Z } ] − η [ E { Y ( x ) ∣ Z } ] , which measures the conditional (on Z ) effect of increasing the exposure with 1 unit on the (transformed by η ) mean potential outcome. The formulations in (10) and (11) cover more complex models as well, including, for instance, splines, polynomials and interaction terms. Neither the conditional model (10) nor the marginal model (11) makes any distributional assumptions, apart from their respective mean specifications.

To address the fact that the potential outcome in regression models (10) and (11) is incompletely observed, we make use of relation (9). For this, we first need to estimate the mean E ( Y ∣ X = x , Z ) . Together with our choice of c ( x , x ′ , Z ; α ) and α , this will allow us to compute a “bias-corrected” outcome and fit the causal model with standard software for generalized linear models. To estimate the mean E ( Y ∣ X = x , Z ) in relation (9), we use a “nuisance model” of the form

where β is a column vector of parameters and h ( X , Z ) is a regression function with the same dimension as β . This model can be fitted with standard methods (e.g., maximum likelihood) to produce an estimate β ˆ . Since the outcome nuisance model is not of interest per se, there is no need to keep this model simple to have a transparent interpretation of β . Rather, it would be typically advantageous to use as flexible model as possible, to reduce the risk of bias in the estimated target causal effect due to model misspecification. When both X and Z are low dimensional (e.g., binary), the outcome nuisance model may even be saturated.

The relation in (9) contains the exposure as both a fixed constant x and a random variable X ; the former corresponds to the counterfactual scenario where everybody receives exposure level x , whereas the latter corresponds to the factual scenario where the exposure varies randomly across subjects. In practice, these scenarios cannot be observed simultaneously. However, we can “simulate” the counterfactual scenario by artificially assigning a fixed exposure level to each subject and use the postulated confounding function to compute the potential outcome under this scenario, thereby being able to “observe” both scenarios simultaneously. This intuition motivates us to “augment” the data used to fit the models in (10) and (11), whereby, for each subject i , we generate multiple realizations of the outcome, corresponding to different values of x . More precisely, for each subject i , we generate a random draw of an auxiliary multivariate variable V i = ( V i 1 , … , V i K ) T for each subject i from a pre-specified distribution independently of X . We discuss the choice of distribution below. For each subject i and observation V i k , we define the bias-corrected outcome

Thus, for each subject i , we generate K bias-corrected outcomes. The term β ˆ T h ( V i k , Z i ) is a prediction of the transformed mean η { E ( Y ∣ X = V i k , Z i ) } and can easily be obtained from standard software after having fitted the outcome nuisance model (12). A similar data augmentation technique was used by ref. [28], for estimation of natural direct and indirect effects. We show in Appendix B that the bias-correction in (13) removes the confounding bias, in the sense that E ( Y i k , α , β ∣ V i k = x , Z i ) = E { Y ( x ) ∣ Z i } .

For fixed α and β ˆ , we estimate the parameter vector ψ in the conditional model (10) as the solution to the equation system:

where H i ( α , β ˆ , ψ ) = ∑ k = 1 K H i k ( α , β ˆ , ψ ) and

Estimation of ψ in the marginal model (11) is carried out with the same sets of equations, but with g ( V i k , Z i ) replaced by g ( V i k ) . We show in Appendix B that the estimating equation system is unbiased, i.e., that E { H i ( α , β , ψ ) } = 0 at the true value of ( α , β , ψ ) . It follows from the standard theory for M-estimators [29] that ψ ˆ is consistent and asymptotically normal, provided that the confounding function c ( x , x ′ , Z ; α ) is correctly specified with correct value of the sensitivity parameter α and that both the causal model (10) (or (11)) and the outcome nuisance model (12) are correct.

The equation system in (14) is generally nonlinear and has one equation for each element in ψ . Thus, if ψ is high dimensional, it may be hard and/or time consuming to find a solution by attempting to solve the equation system explicitly. Fortunately, a solution can be easily obtained by noting that the estimating function H i k ( α , β ˆ , ψ ) is identical to the maximum likelihood score function for a generalized linear model of the form

with canonical link function η . Thus, once the bias-corrected outcomes have been computed, ψ can be estimated by fitting the model in (16) with standard software for generalized linear models. Such software typically maximize a scalar likelihood function rather than solving a potentially large set of score equations, and are usually highly optimized. Thus, this estimation method is expected to be quick and stable even for large data sets and high dimensional parameter vectors ψ .

The default standard errors obtained from the software are generally not correct, since they ignore the uncertainty in β ˆ and the fact that the K score contributions from each subject are correlated. This problem can be solved by using bootstrap standard errors, provided that the outcome nuisance model in (12) is re-fitted and the bias-corrected outcomes are re-computed for each bootstrap sample. Alternatively, analytic standard errors can be obtained with theory for M-estimators. Specifically, define θ = ( β T , ψ T ) T and the stacked estimating function

where S i ( β ) is the maximum likelihood score contribution from the outcome nuisance model for subject i . It follows from the standard theory for M-estimators [29] that the asymptotic variance for θ ˆ is given by

An estimate of the asymptotic variance is obtained by replacing the true value of θ in (17) with the estimate θ ˆ and the population moments with their sample counterparts.

Consistency of ψ ˆ does not require any particular auxiliary variable distribution p ( V ∣ Z ) ; in particular, it does not require that p ( V i k = x ∣ Z ) = p ( X = x ∣ Z ) . However, the asymptotic variance of ψ ˆ may depend on the auxiliary variable distribution, with some distributions giving more efficient estimators. For instance, in analogy with ordinary linear regression, we would expect that the estimated “regression slope” ψ ˆ is less stable if the “regressors” V i 1 , … , V i K are close to each other than if they are spread out. At the extreme end the “regressors” are all equal, V i 1 = , … , = V i K = x ∗ , say, in which case, it is clearly not possible to estimate E { Y ( x ) ∣ Z } for any other value than x ∗ , so that the estimated “regression slope” ψ ˆ has infinite variance. When X is categorical with levels x 1 , … , x J , we conjecture that a reasonable efficiency may be obtained by setting K = J and V i k = x k for k = 1 , … , J . When X is continuous with quartiles q 1 , q 2 , q 3 , we may set K = 3 and V i k = q k for k = 1 , … , 3 .

6 Calibration of the sensitivity parameter α

The technical definition of conditional exchangeability in (3) is quite abstract, and it may be hard, even for a subject matter expert, to determine how the impact of unmeasured confounders translates into particular values of the sensitivity parameter α , for a chosen parametric form of the confounding function (8). To aid in this process, we propose a method for calibration of α , which uses information about the measured confounders Z , and is applicable when the confounding function c ( x , x ′ , Z ; α ) of interest is assumed not to depend on Z , e.g., c ( x , x ′ , Z ; α ) = α ( x − x ′ ) .

To this end, for any disjoint sets of confounders W 1 and W 2 , we first define the “extended” confounding function for W 2 , given W 1 , as follows:

where the outer expectation in the second term on the right-hand side is taken over the conditional distribution of W 2 , given ( X = x ′ , W 1 ) . The purpose of this function is to quantity the potential for W 2 to further confound the association between X and Y , once W 1 has been controlled for. If either W 2 is conditionally independent of X , given W 1 , or conditionally independent of Y , given ( X , W 1 ) , then we would expect that W 2 has no confounding potential, given W 1 . Indeed, it can be shown (see Appendix C) that

and

Conversely, the stronger the deviations from the conditional independencies in (19) and (20), the larger the magnitude of d W 2 ( x , x ′ , W 1 ) we would generally expect. The extended confounding function has a simple relation to the ordinary confounding function in (8). If W 1 and W 2 are two disjoint sets of confounders requiring control to achieve exchangeability, then it can be shown (see Appendix C) that the extended confounding function for W 2 , given W 1 , is identical to the ordinary confounding function given W 1 :

These relations motivate the following calibration strategy. Let U be an additional set of confounders requiring control to achieve exchangeability, given Z , so that d U ( x , x ′ , Z ) = c ( x , x ′ , Z ) . We first select a subset Z ∗ of Z for which the need for control is fairly certain; this subset can be the full set of all measured confounders, i.e., Z ∗ = Z . We then consider the extended confounding function d Z ∗ ( x , x ′ , ∅ ) for Z ∗ given the empty set; this function measures the “marginal” confounding potential of Z ∗ . We further assume that this function has the same parametric form as the ordinary confounding function c ( x , x ′ , Z ; α ) for U , given Z , that we are ultimately interested in, and that c ( x , x ′ , Z ; α ) does not depend on Z. For example, for a confounding function given by c ( x , x ′ , Z ; α ) = α ( x − x ′ ) , we would consider the extended confounding function d Z ∗ ( x , x ′ , ∅ ; α ∗ ) = α ∗ ( x − x ′ ) . Since Z ∗ is in fact measured, we can estimate α ∗ , thereby obtaining a concrete reference point for the sensitivity parameter α of main interest. For instance, if we believe that the confounding potential of U , given Z , is no larger in magnitude than the marginal confounding potential of Z ∗ , then we would consider values of α in the range ( − ∣ α ∗ ∣ , ∣ α ∗ ∣ ) as plausible.

To estimate α ∗ , we use the fact that

where the first equality follows from (18) and the second from the law of total probability. This suggests that α ∗ can be estimated as the value that minimizes the discrepancy between the left-hand side and the right-hand side of (21). Specifically, we propose the following five-step estimation strategy:

1. Fit regression models E ( Y ∣ X , Z ∗ ; ξ ) and E ( Y ∣ X ; δ ) are ideally as flexible as possible, e.g., saturated for binary X and Z ∗ .

2. Let X be the set of possible values for X . If X is continuous, we suggest to provisionally categorize X to obtain a finite set X . For each x in X , use the regression model E ( Y ∣ X , Z ∗ ; ξ ) to create an estimate of E { E ( Y ∣ X = x , Z ∗ ) } , which we denote by μ ˆ 1 ( x ; ξ ˆ ) . This is done with regression standardization [30], by using the model to create a prediction E ( Y ∣ X = x , Z i ∗ ; ξ ˆ ) for each subject i and then averaging these predictions:

   μ ˆ 1 ( x ; ξ ˆ ) = ∑ i = 1 n E ( Y ∣ X = x , Z i ∗ ; ξ ˆ ) / n .

3. For each x in X , and a fixed value α ∗ , use the regression model E ( Y ∣ X ; δ ) to create another estimate of E { E ( Y ∣ X = x , Z ∗ ) } , which we denote by μ ˆ 2 ( x ; δ ˆ , α ∗ ) . This is done by using the model to create an estimate E ( Y ∣ X = x ; δ ˆ ) and then combining this with the extended confounding function d Z ∗ ( x , X , ∅ ; α ∗ ) to create a prediction E { E ( Y ∣ X = x , Z ∗ ) ∣ X i ; δ ˆ , α ∗ } = η − 1 [ η { E ( Y ∣ X = x ; δ ˆ ) } − d Z ∗ ( x , X i , ∅ ; α ∗ ) ] for each subject i , and then finally averaging these predictions:

   μ ˆ 2 ( x ; δ ˆ , α ∗ ) = ∑ i = 1 n E { E ( Y ∣ X = x , Z ∗ ) ∣ X i ; δ ˆ , α ∗ } / n .

   Note that μ ˆ 2 ( x ; δ ˆ , α ∗ ) depends on the value of α ∗ through the extended confounding function d Z ∗ ( x , X i , ∅ ; α ∗ ) .

4. Define the estimated squared distance D ˆ ( x ; ξ ˆ , δ ˆ , α ∗ ) = { μ ˆ 1 ( x ; ξ ˆ ) − μ ˆ 2 ( x ; δ ˆ , α ∗ ) } 2 and compute the estimated sum of squared distance:

   D ˆ ( ξ ˆ , δ ˆ , α ∗ ) = ∑ x ∈ X D ˆ ( x ; ξ ˆ , δ ˆ , α ∗ ) .

5. Estimate α ∗ as the value that minimizes the estimated sum of squared distance:

   α ˆ ∗ = argmin α ∗ D ˆ ( ξ ˆ , δ ˆ , α ∗ ) .

We note that, in contrast to the ordinary confounding function c ( x , x ′ , Z ) , the extended confounding function d Z ∗ ( x , x ′ , ∅ ) is nonparametrically identifiable from the observed data. However, it is only nonparametrically estimable if Z ∗ is categorical. Thus, if Z ∗ is categorical, we may informally verify the parametric model d Z ∗ ( x , x ′ , ∅ ; α ∗ ) by comparing with a fully saturated model for d Z ∗ ( x , x ′ , ∅ ) ; a large agreement between these indicates that the parametric model fits well to the data.

7 Simulation

We carry out a small simulation study to confirm the theoretical properties of our estimation and calibration methods, in the ideal scenario when all model assumptions are correct. We generate data from the model

where U and Z represent an unmeasured and measured confounder, respectively. We assume conditional exchangeability given ( U , Z ) , so that the last row of (22) can be interpreted as the conditional causal model

Thus, ψ is the conditional causal odds ratio, given ( U , Z ) , and the extended confounding function d U ( x , x ′ , Z ) is equal to the ordinary confounding function c ( x , x ′ , Z ) . We set ψ = 0 , thus simulating under the causal null. For a fixed value of α , we generate data so that d Z ( x , x ′ , ∅ ; α ∗ ) = α ∗ ( x − x ′ ) and d U ( x , x ′ , Z ) = c ( x , x ′ , Z ; α ) = α ( x − x ′ ) for Z ∈ { 0 , 1 } , with α = α ∗ , so that the confounding potential of U , given Z , is equal to the marginal confounding potential of Z . These relations imply the following six constraints on the observed data distribution:

For a fixed value of α , we use these constraints to determine the six parameters ( ν 0 , ν 1 , ν 2 , τ 0 , τ 1 , and τ 2 ); we provide details in Appendix E. We note that the obtained parameter values may not be unique; there could possibly exist different sets of parameter values that satisfy the constraints for a given value of α .

For each α ∈ { 0.05 , 01 , 0.25 } , we generate 300 samples of size n = 1,000 from the model mentioned earlier. For each sample, we fit the conditional causal model logit [ p { Y ( x ) = 1 ∣ Z } ] = ψ 0 + ψ 1 x + ψ 2 Z with the estimation methods outlined in Section 5. This model is correctly specified and saturated under our data-generating mechanism, with ψ 1 = ψ = 0 . We fit the model twice: first using the correct value of α , then estimating α according to the calibration scheme outlined in Section 6 and using the estimated value of α in the estimation of ψ . In this calibration, we used Z ∗ = Z and saturated models for E ( Y ∣ X ) and E ( Y ∣ X , Z ) , and correct models for the (extended) confounding functions. We estimate the standard error with both the sandwich formula and a nonparametric bootstrap, using 100 bootstrap samples for each original sample. In the bootstrap, we re-estimate α for each bootstrap sample.

Table 2 shows the regression coefficients ( ν 0 , ν 1 , ν 2 , τ 0 , τ 1 , and τ 2 ) for each value of α . Table 3 shows the mean (over the 300 samples) estimate ψ ˆ together with its mean sandwich standard error, mean bootstrap standard error and its empirical (over the 300 samples) standard deviation, when using the true value of α in the estimation. Table 4 shows the same statistics when using the estimated values of α ; this table also shows the mean estimate α ˆ .

We observe that ψ ˆ is virtually unbiased ( ≈ 0 ) for all values of α , both when using the true value and the estimated value of α in the estimation of ψ . We further observe that α ˆ is virtually unbiased as well. The sandwich and bootstrap standard errors of ψ ˆ agree well with the empirical standard deviation of ψ ˆ when using the true value of α . When using the estimated value of α the sandwich standard error is slightly too small; this is because it does not account for the sampling variability in α ˆ .

8 Illustration

Chiba [11] analyzed data from the Western Collaborative Group Study (WCGS); a cohort study that aimed to assess risk factors for CHD. In this section, we illustrate our sensitivity analysis by applying it to the WCGS data. As a comparison, we provide a sensitivity analysis conducted using the method presented by Franks et al. [14] in Appendix F. The R-code for the analysis is provided in Appendix G. The WCGS data are publicly available in the R package faraway, which enables the reader to replicate and elaborate on our analyses.

The cohort consists of 3,154 middle-aged men employed at 10 California companies in 1960. Several covariates were measured at baseline, including behavior pattern classified as type A (aggressive and ambitious) or type B (calm and relaxed), age, body mass index, systolic and diastolic blood pressure, and smoking (yes/no). The subjects were followed until death or end of 1969, whichever came first. During follow-up, all diagnoses of CHD were recorded. In line with Chiba [11], we ignore complications due to censoring by death and consider CHD at end of follow-up as a binary outcome: Y = 1 for CHD and Y = 0 for no CHD.

8.1 Analysis with a binary exposure

Following Chiba [11], we start by considering the behavior pattern as the exposure ( X = 1 for type A and X = 0 for type B) and the other measured covariates as potential confounders Z . Chiba [11] aimed to estimate the marginal causal risk ratio for behavior pattern on CHD, whereas we focus on the marginal causal odds ratio. Specifically, we consider the marginal causal logistic model

(25) logit { E ( Y ( x ) ) } = ψ 0 + ψ 1 x ,

where the causal log odds ratio ψ 1 is of main interest. Since X is binary, this model is saturated.

To estimate ψ while assessing its sensitivity to unmeasured confounding, we use the confounding function c ( x , x ′ , Z ; α ) = α ( x − x ′ ) , with η being the logit link. Chiba [11] used the same confounding function, but with a log link. We use the outcome nuisance model

(26) logit { E ( Y ∣ X , Z ) } = β 0 + β 1 X + β 2 T Z + β 3 T X Z ,

which includes both main effects of X and each element of Z , and interaction terms between these.

To appreciate what values of α can be considered plausible, Chiba [11] computed the crude log risk ratio between smoking and CHD, which was found to be equal to 0.58. He then considered values of α in the range ( − 0.58 , 0.58 ) as plausible, thus presumably thinking of the unmeasured confounders as being about “equally strong” as smoking. However, while the sensitivity parameter α depends on the smoking-CHD risk (or odds) ratio, it does this in a complex fashion. Thus, there is no reason why the parameters should be similar, even if smoking was sufficient for confounding control. In particular, α also depends on the association between smoking and behavior pattern, and equals 0 if these are independent, even if smoking and CHD are strongly associated.

Instead, using our proposed method to calibrate α , with Z ∗ being smoking, gives a logically consistent range for α . Since both behavior pattern and smoking are binary, we use saturated logistic models for E ( Y ∣ X , Z ∗ ) and E ( Y ∣ X ) in the first step of the estimation, obtaining α ˆ ∗ = 0.03 . This indicates that, if we believe that the confounding potential of the unmeasured confounders, given the measured confounders, is no larger than the marginal confounding potential of smoking, then we may consider values of α in the range ( − 0.03 , 0.03 ) , which is an order of magnitude narrower than the range ( − 0.58 , 0.58 ) considered by Chiba [11]. As an informal model verification, we fit the saturated model d Z ∗ ( x , x ′ , ∅ ; α ∗ ) = α 1 ∗ I ( x = 1 , x ′ = 0 ) − α 0 ∗ I ( x = 0 , x ′ = 1 ) for the extended confounding function, obtaining α ˆ 1 ∗ = 0.028 and α ˆ 0 ∗ = 0.051 . Neither of these values are very far off from the estimate α ˆ ∗ = 0.03 obtained from the nonsaturated model d Z ∗ ( x , x ′ , ∅ ; α ∗ ) = α ∗ ( x − x ′ ) , which gives some credibility to this model.

We next repeat the calibration, now letting Z ∗ be the set of all measured confounders, Z ∗ = Z . For this, we use the logistic model logit { E ( Y ∣ X ) } = δ 0 + δ 1 X and the logistic model for E ( Y ∣ X , Z ) = E ( Y ∣ X , Z ∗ ) in (26). The resulting estimate of α ∗ is 0.19, which indicates that, if we believe that the confounding potential of the unmeasured confounders, given the measured confounders, is no larger than the marginal confounding potential of the whole set of measured confounders, then we may consider values of α in the range ( − 0.19 , 0.19 ) .

The left panel of Figure 1 shows the estimate of ψ 1 (solid line) together with pointwise 95% confidence limits (dashed lines), as a function of α in the range ( − 0.19 , 0.19 ) . At α = 0 (no deviation from conditional exchangeability, given the measured confounders), ψ ˆ 1 = 0.67 with 95% confidence interval (0.39, 0.94). As α goes from − 0.19 to 0.19, ψ ˆ 1 goes almost linearly from 0.85 to 0.48. Everywhere in this range, the lower confidence limit for ψ 1 is above 0. Thus, if we do not believe that the unmeasured confounders have a stronger potential than the measured confounders, then we may declare the behavior pattern to have a causal effect on CHD. The dotted line shows the estimated marginal log risk ratio, as obtained with the method given by Chiba [11]. In the range ( − 0.19 , 0.19) for α , this is very similar to the marginal log odds ratio.

Figure 1

The estimated marginal causal log odds ratio (solid lines) and log risk ratio (dotted lines) for behavior pattern and CHD together with pointwise 95% confidence limits (dashed lines) for the marginal causal log odds ratio, as functions of the sensitivity parameter α .

The right panel of Figure 1 shows the same quantities as the left panel, but for a much wider range of α . At α = 0.70 and α = 0.41 , the estimate of ψ 1 and its lower confidence limit, respectively, reach 0. Thus, taking statistical uncertainty into account, the unmeasured confounders would have more than twice as strong confounding potential as the measured confounders, on the scale defined by α , to fully explain away the observed association between the behavior pattern and CHD. This type of reasoning, where one derives a condition under which the causal effect is nullified and then evaluates the plausibility of that condition, is analogous to the reasoning that underlies the popular “ E -value” [7], but expressed in terms of our sensitivity parameter. As α becomes large in magnitude, the estimated causal log odds ratio deviates dramatically from the estimated causal log risk ratio. While the former seems to converge to ≈ 3.87 and ≈ − 2.91 as α goes to − ∞ and ∞ , respectively, the latter appears to approach ∞ and − ∞ almost linearly. However, the causal log risk ratio and odds ratio cannot be arbitrarily large; [31] showed that the observed data distribution implies restrictions on causal effects, the so-called bounds. For the WCGS data, these bounds are given by ( − 2.24 , 3.09) and ( − 2.93 , 3.87) for the causal log risk ratio and odds ratio, respectively, ignoring sampling variability. The sensitivity analysis with a logit link respects these bounds, whereas the sensitivity analysis with a log link does not. The reason for this is that the confounding function with a log link is not variation independent of the observed data distribution, as discussed in Section 4.

One may argue that the simple confounding function c ( x , x ′ , Z ; α ) = α ( x − x ′ ) is not entirely realistic in the current example. For instance, CHD is a largely age-dependent disease, which is quite rare at young ages. Thus, one may suspect that the factors and mechanisms that cause CHD are very different for young people than for old people, in which case the potential for unmeasured confounding bias may be different for young and old people as well. To reflect this, we consider the confounding function c ( x , x ′ , Z ; α ) = α 0 ( x − x ′ ) + α 1 ( x − x ′ ) ( age − k ) , where the magnitude of α 1 determines the degree of such an interaction mechanism between the unmeasured confounders and the age. To gain intuition for this confounding function, consider two ages, age 1 and age 2 , such that age 1 > k > age 2 and Δ = ∣ age 1 − k ∣ = ∣ age 2 − k ∣ , and define r = α 1 / α 0 . The ratio of confounding functions at age 1 and age 2 is then given by

R = α 0 ( x − x ′ ) + α 1 ( x − x ′ ) ( age 1 − k ) α 0 ( x − x ′ ) + α 1 ( x − x ′ ) ( age 2 − k ) = 1 + Δ r 1 − Δ r .

It is easy to show that R is increasing in Δ r , and that R > 0 if − 1 < Δ r < 1 . The age range in the WCGS data is 39 to 59 years. We set k = 49 and consider values of r in the range ( − 1 / 30 , 1 / 30 ) . It follows that − 1 / 30 < Δ r < 1 / 30 , which implies that 0.5 < R < 2 . Thus, for this confounding function and choice of parameter values, the deviation from exchangeability may be up to twice as large at any age in the interval (39, 59), as compared to any other age in the interval.

The contour plots in Figure 2 show the estimate of ψ 1 (left panel) and the lower limit of its 95% confidence interval (right panel) as functions of α 0 in the range ( − 0.19 , 0.19) and α 1 in the range ( − ∣ α 0 ∣ / 30 , ∣ α 0 ∣ / 30 ) . We observe that the ranges of ψ ˆ 1 and its lower confidence limit in Figure 2 (0.47, 0.87) and (0.19, 0.59), respectively, are approximately the same as the corresponding ranges in Figure 1. Thus, the inclusion of the interaction term α 1 ( x − x ′ ) ( age − k ) in the confounding function does not give qualitatively different conclusions for this example.

Figure 2

The estimated marginal causal log odds ratio for behavior pattern and CHD (left panel) and its lower pointwise 95% confidence limit (right panel), as functions of the sensitivity parameters α 0 and α 1 . The gray areas correspond to values of α 1 outside the range ( − ∣ α 0 ∣ / 30 , ∣ α 0 ∣ / 30 ) .