Selection Bias When Using Instrumental Variable Methods to Compare Two Treatments But More Than Two Treatments Are Available

1 Introduction

Confounder adjustment is critical in the estimation of treatment effect in observational studies. In practice, however, there is no guarantee that all the confounders (i.e., baseline variables that affect both the treatment and outcome) have been measured. These unmeasured confounders can bias the treatment effect estimation.

Instrumental variable (IV) methods are designed to adjust for the bias caused by unmeasured confounders in estimating a difference measure of effect such as the risk difference. IVs are predictors of treatment allocation that are independent of unmeasured confounders and do not have a direct effect on the outcome (i.e., only affect the outcome through affecting the treatment). Intuitively, the IV method seeks to extract variation in treatment that is free of unmeasured confounders and uses this variation to estimate the treatment effect. For more information on IV methods, see Angrist et al.[1], Newhouse and McClellan [2], Greenland [3], Hernán and Robins [4], Cheng et al. [5], Baiocchi et al.[6] and Imbens [7]. This paper is motivated by provider preference IV (PP IV) which is commonly used as an IV in health studies. PP IV utilizes the natural variation in medical practices to construct an IV. For example, in studying the effect of cyclooxygenase-2 (COX-2) inhibitors vs. nonselective nonsteroidal anti-inflammatory drugs (nonselective NSAIDs) on gastrointestinal complications Brookhart et al. [8] used whether a patient’s physician last prescription for a patient prescribed a nonsteroidal anti-inflammatory drug (NSAID) was for a Cox-2 inhibitor or a nonselective NSAID as an IV for whether the patient was prescribed a Cox-2 inhibitor or a nonselective NSAID. Other examples of PP IV can be found in Brooks et al. [9] and Brookhart and Schneeweiss [10].

For a binary IV and binary treatment, Angrist et al. [1] proposed a nonparametric estimator which identifies the treatment effect among subjects who would take the treatment if encouraged to do so by the IV and not take the treatment if not encouraged to do so by the IV (compliers). This estimator is equivalent to a 2 stage least square (2SLS) estimator where at the first stage the treatment is regressed on IV and at the second stage the outcome is regressed on fitted values of the first stage regression (i.e., predicted value of the treatment variable given the IV). The coefficient of the independent variable of the second stage regression is equivalent to the estimator proposed by Angrist et al. [1]. The 2SLS estimator is also called the Wald estimator after Wald [11]. More efficient versions of the 2SLS estimator has been developed by Imbens and Rubin [12, 13] and Cheng et al. [14].

Although treatment effect estimation using IV approaches for binary treatments is straightforward, it is challenging when there are more than two treatment arms. In fact, Cheng and Small [15] shows that the treatment effects are not point identifiable and proposes bounds on treatment effects using method of moments. Long et al. [16] studies the same problem and proposes a Bayesian approach for finding credibility intervals for the identification region.

Sometimes investigators are interested in comparing the effect of two treatments while the data consist of more than two treatment options for patients. In this situation a common approach is to select patients who are assigned to the treatments of interest and perform the standard IV approach for binary treatment and IV to estimate the treatment effect [17–21].

The objective of this manuscript is to provide empirical and theoretical evidence that, in general, excluding patients based on their assigned treatment can result in selection bias. This problem has also been discussed in an independent work by Swanson et al. [22] and Swanson [23]. In the current manuscript, we formalize the mechanism of this selection bias within a principal stratification framework and provide a procedure that identifies the treatment effect of interest as a function of a vector of sensitivity parameters. Specifically, the sensitivity parameters are used to estimate the probability of being assigned to one of the treatments of interest as a function of the observed outcome and the IV value [24]. We also list assumptions under which selecting patients based on their received treatment does not lead to a biased treatment effect estimate.

The manuscript is organized as follows. Section 2 provides the setup and defines the causal estimand. It also lists the required assumptions for our procedure. In Section 3, we discuss the selection bias issue and present the proposed sensitivity analysis method. We report on the results of a simulation study in Section 4, where we examine the performance of our proposed method. In Section 5, we study the comparative effect of metformin and sulfonylureas on weight gain using The Health Improvement Network (THIN) database.

2 Settings and notations

For simplicity we focus on observational studies with three treatment options A, B, and C. We assume that we are only interested in comparing the effect of treatment A with B. We use small letters to refer to the possible values of the corresponding capital letter random variable. Let Z be a binary IV where Z∈{A,B} such as provider preference. Z=A means that the IV predisposes the patient to receive treatment level A among the options A and B, and Z=B means that the IV predisposes the patient to receive treatment level B among the options A and B; e.g., in Brookhart et al. [8] study, if A denotes Cox-2 inhibitor and B denotes nonselective NSAID, Z=A if the patient’s physician last prescription to a patient prescribed an NSAID was a Cox-2 inhibotor and Z=B if it was a nonselective NSAID. Let D(z) be the potential treatment assigned given Z=z, where D(z)∈{A,B,C}, i.e., D(A) is the treatment that would be received if the IV was A. The potential outcomes are Y(z,d) for Z=z and D=d, i.e., Y(A,B) is the outcome that would be observed if the IV was A and the treatment level was B. The observed treatment and outcome are D=D(Z) and Y=Y(Z,D(Z)), respectively.

Frangakis and Rubin [25] developed a basic principal strata (PS) framework which stratifies subjects based on their potential assigned treatments. In our setting, we can classify subjects into the following 9 principal strata: S1) always B taker, S2) always A taker, S3) always C taker, S4) subjects who would comply with their value of IV, i.e., D(A)=A, D(B)=B, S5) subjects who would take B if Z=B and would take C if Z=A, S6) subjects who would take C if Z=B and would take A if Z=A, S7) subjects who would take A if Z=B and would take B if Z=A, S8) subjects who would take B if Z=A and would take C if Z=B, and S9) subjects who would take A if Z=B and would take C if Z=A. We denote the proportions in the basic principal strata as πS..

We assume that the IV satisfies the following standard assumptions [1]: (1) stable unit treatment value assumption, (2) The IV is independent of unmeasured confounders, potential outcomes and potential treatment received, (3) the IV affects the treatment but has no direct effect on the outcome. This assumption is known to as the exclusion restriction (ER) assumption. Under the exclusion restriction, we let Y(d) be the potential outcome if treatment d is received which equal Y(z=0,d)=Y(z=1,d), and (4) the IV is independent of the pretreatment covariates, potential outcome and potential treatments. Furthermore, we assume two types of monotonicity assumptions: (5) monotonicity I: there is no defiers (i.e., p(D(A)≠A,D(B)≠B)=0) and (6) monotonicity II: there is no one who would take B if Z=A but take C if Z=B and no one who would take C if Z=A but take A if Z=B. Note that monotonicity assumptions 5 & 6 reduce the number of principal strata to 6. The parameter of interest is defined as

3 IV selection bias

We show in this manuscript that IV analysis methods may not result in an unbiased treatment effect if applied on a data set in which units are preselected based on their received treatments. A standard IV analysis of patients on the subset of patients for whom D=A or B is essentially a comparison of the Z=A vs. Z=B patients on this subset. If Z=A has a different effect on whether D=C than Z=B, then there will be selection bias. For example, patients who are in a more severe stages of the disease may be treated with treatment D=C and therefore have less chance to be observed in the sample; if doctors who prefer drug A(Z=A) are more likely to use treatment C with severe patients than doctors who prefer drug B(Z=B), then there will be selection bias. Another scenario where such selection bias might arise is that treatment C is more likely to be offered if some complications are present and doctors who prefer drug A are more likely to treat a patient with complications with drug C than doctors who prefer drug B.

The challenge is that θ is not identifiable using the observed data when there is selection bias. This is due to the fact that the observed outcomes of Y|Z,D is a mixture of potential outcomes from principal strata. Let fA(y) be the density of outcome for a patient who received D=A when assigned to Z=A (i.e., fA(y)=f(y|D=A,Z=A). Then

where γAS.=πS.πS2+πS4+πS6 and fAS.(y) is the density function of the outcome for subjects in principal stratum S. with D=A. The density fAS2(y) and fAS6(y) can be written as

Thus,

Similarly, let fB(y) be the density of outcome for a patient who received D=B when assigned to Z=B. Then

where γBS.=πS.πS1+πS4+πS5 and fBS.(y) is the density function of the outcome for subjects in principal stratum S. with D=B. Also,

Thus,

The mixture proportions γAS. and γBS. are unknown and need to be estimated using the observed data. However, the strata proportions πS. are not identifiable using the proportions in the observable (D,Z) strata, because

does not have a unique solution. Also, the density functions fAS4(y) and fBS4(y) are not identifiable using the observed data because p(D(B)=A|D(A)=A,Y(A)=y), p(D(B)=C|D(A)=A,Y(A)=y), p(D(A)=B|D(B)=B,Y(B)=y) and p(D(A)=C|D(B)=B,Y(B)=y) are unknown. So we assume logit models for these probabilities and perform the sensitivity analysis based on these models;

Note that for any fixed (α1A,α1B), (α0A,α0B) can be determined since ∫fAS4(y)dy=1 and ∫fBS4(y)dy=1. Therefore, the parameter of interest θ can be estimated using the four dimensional sensitivity parameter (γAS4,γBS4,α1A,α1B).

Fixing α1A<0 makes the conditional probability monotone decreasing in y, which means the odds of being a complier among patients who have taken D=A if Z=A (e.g., seeing a physician with A preference) is lower for a larger Y(A)=y. Similarly, α1A>0 makes the conditional probability monotone increasing in y, which means the odds of being a complier among patients who have taken D=A if Z=A is higher for a larger Y(A)=y. The parameter α1B has a similar interpretation. The parameter γAS4 is the probability of being a complier (i.e., take D=A if Z=A and D=B if Z=B) for patients who have taken D=A when Z=A (i.e., γAS4=p(D(B)=B|D(A)=A). The parameter γBS4 is defined similarly.

3.1 No selection bias

Let R be a binary variable which is 1 if a patient is assigned to one of the treatments of interest and 0 otherwise. Under either of the following two conditions, the treatment effect is identifiable when we have access to the entire data (i.e., data includes patients who have received treatment C as well as A and B):

1. p(R=0|Z,U,X)=p(R=0|U,X) where X and U are measured and unmeasured confounders, respectively.

2. p(R=0|Z,U,X)=p(R=0|Z,X).

Assumption A.1 implies that the IV is unrelated to the patient’s selection. In other words, a patient would have the same chance of being selected for both values of the IV. A.2 means that the selection probability is independent of the unmeasured confounders given the IV and measured confounders. Figure 1 depicts the association between different variables. The left plot presents a scenario in which analyzing the preselected data will cause bias because conditioning on R will open the path Z→R→U→Y which violates one of the IV assumptions. This path is closed under either of assumptions A.1 or A.2 as shown in the center and right plots, respectively.

Figure 1:

Causal diagram. Conditioning on subjects with R=1 leads to a bias treatment effect estimate in the left figure. Assumptions A.1 and A.2 hold for the center and right figure, respectively.

Under either of assumptions A.1 or A.2, the treatment effect can be identified. Specifically, under assumption A.1, the principal strata S5 and S6 does not exist. Thus

E[Y|Z=A,D=A]=γS2E[Y(A)|PS=S2]+γS4E[Y(A)|PS=S4]

E[Y|Z=B,D=B]=γS1E[Y(B)|PS=S1]+γS4E[Y(B)|PS=S4]

E[Y|Z=A,D=B]=γS1E[Y(B)|PS=S1]

E[Y|Z=B,D=A]=γS2E[Y(A)|PS=S2].

These equations can be used to identify θ=E[Y(A)−Y(B)|PS=S4]. Also, under assumption A.2, θ can be estimated using the following three-step procedure:

1. Using the entire data, estimate ϖ=p(R=1|X,Z) by postulating a logit model on R.

2. Include patients who have been assigned to either treatments A or B.

3. Estimate the treatment effect θ using 2SLS IV approach. At each stage fit a weighted least square where weights are ϖˆ.

Fitting weighted least squares accounts for the selection procedure and hence avoids the bias by creating a pseudo-population in which there is no association between (X,Z) and R.

4 Simulation studies

Through simulations we examine the performance of our proposed sensitivity analysis assuming known and unknown (γAS4,γBS4). We assume two scenarios with the following principal strata proportions: Scenario 1:

1. Scenario 1: πS1=πS2=πS3=πS5=0.1, πS4=0.3 and πS6=0.3 which gives us (γAS4,γBS4)=(0.30.7,0.30.5).

2. Scenario 2: πS1=πS2=πS3=0.1, πS4=0.3, πS5=0.4 and πS6=0.0 which gives us (γAS4,γBS4)=(0.20.3,0.20.7).

Table 1 presents the statistical power of detecting the true treatment effect θ=0.0, 0.50, and 0.80 based on 500 datasets of sizes n=500 with (α1A,α1B)=(1,2), (1,1) and (0,0). The nominal type-I error rate is 5%. For both scenarios when the true amount of selection bias is presumed, the empirical and nominal sizes are close. In general, the empirical sizes are inflated when the selection bias parameters are not set correctly. For some of the combinations of the true and presumed sensitivity parameter values the power and type-I error rate are much higher than the nominal one which is caused by the overestimation of the treatment effect. Also, in scenario 2, when there is no selection bias (i.e., true (α1A,α1B)=(0,0)) and we presumed (α1A,α1B)=(1,1) the power is very low because of the underestimation of the treatment effect. Table 2 shows the bias (0.80−θˆ), standard error and mean squared error of the treatment effect estimation for different combinations of the true and presumed sensitivity parameters when the true treatment effect θ=0.80. This table highlights the over- and underestimation phenomenon discussed earlier. Specifically, when the true sensitivity parameters are (α1A,α1B)=(1,2), presuming (α1A,α1B)=(0,2) will result in underestimating the counterfactual mean outcome E[Y(A)|PS=S4]. This is because when we set the presumed α1A less than the true value, the larger outcomes are less likely to be a complier among patients with D=Z=A. In other words it will shift the distribution of this compliers to left. Similarly when the presumed α1A more than the true value, E[Y(A)|PS=S4] will be overestimated (See rows 3 and 4 in Table 2).

Figure 2 shows how the treatment effect estimate changes by the sensitivity parameters when the true (α1A,α1B)=(1,2) and θ=0.80. The darker area in left panel of Figure 2 indicates the smaller bias and the darker area in the right panel shows the region in which the bias is not significant at 5% (i.e., the confidence interval for (0.80−θˆ) contains zero). See also Table 2.

Figure 2:

Simulation study: sensitivity analysis of θ=0.80 when γAS4 and γBS4 are known. The true sensitivity parameter values are (α1A,α1B)=(1,2). The darker area in left panel indicates the smaller bias and the darker area in the right panel shows the region in which the bias is not significant at 5%. The first and second rows represent Scenarios 1 and 2, respectively.

4.1 Unknown (γAS4,γBS4)

In this case our sensitivity parameter space is (γAS4,γBS4,α1A,α1B). Figures 3 and 4 present the bias in the treatment effect estimation for different values of the presumed sensitivity parameters based on scenarios 1 and 2, respectively. In these figures the true (γAS4,γBS4,α1A,α1B)=37,35,1,2 and (γAS4,γBS4,α1A,α1B)=23,27,1,2 in scenarios 1 and 2, respectively. Also the true treatment effect is θ=0.80. The darker area in the left panel shows the region in the sensitivity parameter space in which the bias is smaller and right panel shows the area in which the bias is not significant based on 500 samples (i.e., the confidence interval for (0.80−θˆ) contains zero). The shape of the area which leads to an unbiased estimate changes by the value of the presumed (γAS4,γBS4) and, in general, when the presumed values of (γAS4,γBS4) are far from the true values, the unbiased area is small which means it would be less likely to find an unbiased treatment effect estimate using sensitivity analysis (see also Tables 4 and 5 in the Appendix). We have tabulated the power and type-I error rate for different values of the presumed sensitivity parameters in Table 3. This Table summarizes the result for three different values of the true treatment effect θ=0.00,0.50, and 0.80. The nominal type-I error rate is 5%.

Figure 3:

Simulation study (Scenario 1): Sensitivity analysis of θ=0.80. Sensitivity parameters are (γAS4,γBS4,α1A,α1B). The darker area in left panel indicates the smaller bias and the darker area in the right panel shows the region in which the bias is not significant at 5%. The true vector of sensitivity parameters are (γAS4,γBS4,α1A,α1B)=37,35,1,2.

Figure 4:

Simulation study (Scenario 2): Sensitivity analysis of θ=0.80. Sensitivity parameters are (γAS4,γBS4,α1A,α1B). The darker area in left panel indicates the smaller bias and the darker area in the right panel shows the region in which the bias is not significant at 5%. The true vector of sensitivity parameters are (γAS4,γBS4,α1A,α1B)=23,27,1,2.

Table 3:

Simulation study: Power analysis of θ=0.00,0.50 and 0.80. The true vector of sensitivity parameters in scenarios 1 and 2 are (γAS4,γBS4,α1A,α1B)=(37,35,1,2), and (23,27,1,2), respectively.

Presumed	Scenario 1			Scenario 2
(γAS4,γBS4,α1A,α1B)	0.00	0.50	0.80	0.00	0.50	0.80
(0.3,0.6,1,2)	0.38	0.90	0.98	1.00	1.00	1.00
(0.3,0.6,0,2)	1.00	0.00	0.00	0.07	0.73	0.96
(0.3,0.6,1,0)	1.00	1.00	1.00	1.00	1.00	1.00
(0.3,0.6,0,0)	0.12	0.35	0.80	1.00	1.00	1.00
(0.5,0.5,1,2)	0.40	0.08	0.38	0.93	1.00	1.00
(0.5,0.5,0,2)	1.00	0.00	0.00	0.04	0.40	0.75
(0.5,0.5,1,0)	0.96	1.00	1.00	1.00	1.00	1.00
(0.5,0.5,0,0)	0.10	0.29	0.76	1.00	1.00	1.00
(0.2,0.8,1,2)	0.98	1.00	1.00	1.00	1.00	1.00
(0.2,0.8,0,2)	0.94	0.00	0.05	0.61	1.00	1.00
(0.2,0.8,1,0)	1.00	1.00	1.00	1.00	1.00	1.00
(0.2,0.8,0,0)	0.06	0.24	0.65	1.00	1.00	1.00
(0.7,0.2,1,2)	1.00	0.00	0.00	0.30	0.08	0.32
(0.7,0.2,0,2)	1.00	0.00	0.00	0.98	0.00	0.00
(0.7,0.2,1,0)	0.65	0.99	1.00	1.00	1.00	1.00
(0.7,0.2,0,0)	0.12	0.39	0.87	1.00	1.00	1.00

5 Application to real data

The Health Improvement Network (THIN) is a large database, which contains the electronic medical record of more than 11 million patients. The data contains longitudinal measurements of diagnostic and prescription data and baseline characteristics collected from over 500 general practices in the UK. We are interested in assessing the effect of metformin and sulfonylureas on body mass index (BMI, calculated as mass in kilograms divided by the square of height in meters) among diabetic patients. Our focus is on patients who are taking these treatments as an initial treatment and the BMI is measured two years after treatment initiation. In our analysis we coded metformin and sulfonylureas as A and B, respectively.

In the modern era, metformin is universally acknowledged as the appropriate first-line medication for diabetes type 2 except where contraindicated [26]. While this dominance was being established in the late 1990s and early 2000s, sulfonylureas were the competing first line agent. During that time the prevalence of metformin use in clinical practice rose very quickly, at the expense of sulfonylureas.

The dynamic aspect of the preference justifies the use of provider preference as an IV. We define the preference as a time-varying quantity such that for each practice ID, the IV is defined as an average of metformin use during each two year timeframe from 1998 to 2012. We assign Z=A if the average is more than 50% and B otherwise. Thus, a particular practice ID may have Z=A for one period of time and B for the other.

Our data includes patients who have been assigned to medications other than metformin or sulfonylureas as an initial therapy. Figure 5 presents the sensitivity analysis for the estimation of the treatment effect of interest based on a four dimensional sensitivity parameter (γAS4,γBS4,α1A,α1B).

Figure 5:

THIN data: Sensitivity analysis for comparing the effects on BMI of sulfonylureas vs. metformin. Numbers on the plot present treatment effect estimates. A: metformin; B: sulfonylureas.

The parameter α1A is the log odds coefficient of BMI in predicting the probability that a patient would take D=B (sulfonylureas) if if she saw a physician with sulfonylurea preference (Z=B) given that a patient would take D=A if she saw a physician with metformin preference (Z=A). The parameter γAS4 is the probability of being a complier (i.e., always taking the same medication as the physician preference) for a patient who would have taken D=A if she saw a physician with metformin preference (Z=A). The parameters α1B and γBS4 are defined similarly.

Since there is a clinical sense that sulfonylurea is associated with weight gain [27, 28], physicians may be less likely to prescribe sulfonylurea to high weight patients. This would imply that among patients with D(A)=A, the odds of being a complier should decrease with y. Analogous reasoning suggests that among patients with D(B)=B, the odds of being a complier will increase with y. If this reasoning holds, then we can restrict our sensitivity parameter domain to α1A<0 and α1B>0. Also, on average patients who have taken a sulfonylurea if seeing a physician with a sulfonylurea preference are more likely to be compliers compared with those who have taken metformin given seeing a physician with metformin preference. This means that it is very likely that γAS4<γBS4. In the Appendix, we summarize the results based on the extreme values of the sensitivity parameters (Table 6).

The range of the treatment effect estimate varies by different choices of γAS4 and γBS4. Specifically, for γAS4=0.2 and γBS4=0.6, the treatment effect lies in [−0.5,1]; while for γAS4=0.5 and γBS4=0.3, it lies in [−0.7,0.5]. Based on our sensitivity analysis, it is unlikely that metformin be associated with increase in BMI and a large negative effect (i.e., E[Y(A)−Y(B)|PS=S4]) is also unlikely because it requires very small and large values of α1A and α1B, respectively. Hence, if anything, the treatment effect of interest should lie in [−0.6,−0.1].

6 Discussion

This manuscript shows that IV analysis methods may fail to provide an unbiased estimate of treatment effects when analyzing data in which subjects are selected based on their received treatment. Specifically, this selection bias happens if the chance of including patients in the preselected data varies based on their IV value and/or some unmeasured confounders. For example, suppose we are interested in studying the comparative effect of treatment A and B while there is a third choice of treatment C. Suppose patients with more severe stages of the disease are more likely to be treated with C if seeing a doctor with A preference. Then if we analyzed a data that only includes patients who have received treatment A or B, severe cases will be under represented among patients who have seen doctors with A preference. This means that in the preselected data, the IV is associated with an unmeasured confounder which makes the IV invalid and results in a biased treatment effect estimate.

Within the principal strata framework, we develop a sensitivity analysis that can be used to estimate the treatment effect among compliers as a function of a vector of sensitivity parameters. Specifically, the sensitivity parameters are used to identify the probability of being a complier given the IV value, received treatment and the observed outcome. The dimension of the sensitivity parameter can be reduced if it is plausible to assume that there is no always C takers. This is because the proportions in the principal strata are identifiable which means γAS4 and γBS4 can be estimated using the observed data.