Sensitivity analysis of selection bias: a graphical display by bias-correction index

Selection bias due to of selected participants’ characteristics

Suppose that we are interested in the association between a binary exposure, X, and a binary outcome, Y, where Y =1 is the event of interest. Z denotes a potential confounder of the associations between X and Y. We assume the probability of interested outcome in a population given X and Z as logit  P(Y  = 1|X, Z) =  β₀ + β₁X + β₂Z, where β₁ is a parameter of interest. To describe how selection bias influences parameter estimation, we introduce the indicator S, which equals 1 when an individual is selected for the survey sample and 0 otherwise. The probability of an individual being willing to participate (or being selected) in a study can be expressed as a linear logistic form logit P(S = 1|Y , X, Z) =  α₀ + α₁Y  + α₂X + α₃Z. The linear logistic transformation of outcome prevalence with selection bias among selected participants can be written as $\text{logit}\mathrm{P}\left(\mathrm{Y}=1|\mathrm{X},\mathrm{Z},\mathrm{S}=1\right)=\varnothing \left(\alpha ;\mathrm{Y},\mathrm{X}\right)+{\beta }_0+{\beta }_1\mathrm{X}+{\beta }_2\mathrm{Z}.\varnothing \left(\alpha ;\mathrm{Y},\mathrm{X}\right)$ is a function of coefficients from the equation of $\text{logit}\mathrm{P}\left(\mathrm{S}=1|\mathrm{Y},\mathrm{X}\right).\varnothing \left(\alpha ;\mathrm{Y},\mathrm{X}\right)\text{is equal to}{\alpha }_1+log\frac{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}}}{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_1+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}}}$, which is called BCI. As seen in the equation: logit P(Y  = 1|X, Z, S = 1), BCI is a term that corresponds to the difference between the estimate in the outcome model for the target and that from the sample. Ergo, the BCI represents the bias.

When Z is not only a confounder but also an effect measure modifier in the association between X and Y, it adds complexity to selection bias estimation. We also present another linear logistic form with an interaction term: logit  P(Y  = 1|X, Z, XZ) =  β₀ + β₁X + β₂Z + β₃XZ, where β₁ is a parameter of interest. Z denotes effect measure modifier. The probability of an individual being willing to participate in a study is written as logit P(S = 1|Y , X, Z, XZ) = α₀ + α₁Y  + α₂X + α₃Z + α₄XZ. The outcome prevalence model with selection bias among the selected participants can be written as $\text{logit}\mathrm{P}\left(\mathrm{Y}=1|\mathrm{X},\mathrm{Z},\mathrm{XZ},\mathrm{S}=1\right)=\varnothing \left(\alpha ;\mathrm{Y},\mathrm{X}\right)+{\beta }_0+{\beta }_1\mathrm{X}+{\beta }_2\mathrm{Z}+{\beta }_3\mathrm{XZ}.\varnothing \left(\alpha ;\mathrm{Y},\mathrm{X}\right)$ is equal to ${\alpha }_1+log\frac{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}+{\alpha }_4\mathrm{XZ}}}{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_1+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}+{\alpha }_4\mathrm{XZ}}}$.

There may be many factors influencing participants’ willingness to participate in a study. ∅(α; Y , X) can take different forms, depending on selection probabilities, and these differences directly affect the outcome model. Six scenarios of participant selection models are listed below. Scenarios 1-5 are based on the outcome model without interaction in the population, and Scenario 6 represents a general form of BCI with selection probability depending on Y, X, Z, and XZ.

• Scenario 1: The selection of participants is affected by neither Y nor X.

Selection model: logit P(S = 1|Y , X, Z) = P(S = 1) = α₀

Outcome model with selection bias among selected participants: logit P(Y  = 1|X, Z, S = 1) = β₀ + β₁X + β₂Z

• Scenario 2: The selection of participants is affected only by Y.

Selection model: logit P(S = 1|Y , X, Z) = P(S = 1|Y ) = α₀ + α₁Y

Outcome model with selection bias among selected participant: logit P(Y  = 1|X, Z, S = 1) = ${\beta }_0+{\beta }_1\mathrm{X}+{\beta }_2\mathrm{Z}+\left.\left({\alpha }_1+log\frac{1+{\mathrm{e}}^{{\alpha }_0}}{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_1}}\right.\right)$

• Scenario 3: The selection of participants is affected only by X.

Selection model: logit P(S = 1|Y , X, Z) = P(S = 1|X) = α₀ + α₂X

Outcome model with selection bias among selected participant: logit P(Y  = 1|X, Z, S = 1) = β₀ + β₁X + β₂Z

• Scenario 4: The selection of participants is affected by both Y and X.

Selection model: logit P(S = 1|Y , X, Z) = P(S = 1|Y , X) = α₀ + α₁Y  + α₂X

Outcome model with selection bias among selected participant: $\text{logit}\mathrm{P}\left(\mathrm{Y}=1|\mathrm{X},\mathrm{Z},\mathrm{S}=1\right)={\beta }_0+{\beta }_1\mathrm{X}+{\beta }_2\mathrm{Z}+\left.\left({\alpha }_1+log\frac{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_2\mathrm{X}}}{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_1+{\alpha }_2\mathrm{X}}}\right.\right)$

• Scenario 5:The selection of participants is affected by Y, X and Z.

Z is a confounder variable.

Selection model: P(S = 1|Y , X, Z) = P(S = 1|Y , X, Z)  = α₀ + α₁Y  + α₂X + α₃Z

Outcome model with selection bias among selected participant: $\text{logit}\mathrm{P}\left(\mathrm{Y}=1|\mathrm{X},\mathrm{Z},\mathrm{S}=1\right)={\beta }_0+{\beta }_1\mathrm{X}+{\beta }_2\mathrm{Z}+\left.\left({\alpha }_1+log\frac{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}}}{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_1+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}}}\right.\right)$

• Scenario 6:The selection of participants is affected by Y, X and Z.

Z is an effect measure modifier.

Selection model: P (S = 1 | Y,X,Z,XZ) = P(S = 1 | Y, X,Z,XZ) = α₀ + α₁Y  + α₂X + α₃Z + α₄XZ

Outcome model with selection bias among selected subjects: $\text{logit}\mathrm{P}\left(\mathrm{Y}=1|\mathrm{X},\mathrm{Z},\mathrm{XZ},\mathrm{S}=1\right)={\beta }_0+{\beta }_1\mathrm{X}+{\beta }_2\mathrm{Z}+{\beta }_3\mathrm{XZ}+\left.\left({\alpha }_1+log\frac{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}+{\alpha }_4\mathrm{XZ}}}{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_1+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}+{\alpha }_4\mathrm{XZ}}}\right.\right)$

A detailed derivation is shown in Appendix S1. In Appendix S1, it is shown that the probability P (Y = 1 | X, Z, S = 1) can be expressed as a function of P (Y | X, Z) and P (S | Y, X, Z). This observation makes it possible to express P (Y = 1 | X, Z, S = 1) as an equation containing both ∅(α;Y,X) and β₀ + β₁X + β₂Z.

Correction for selection bias and sensitivity analysis

Scenarios 4–6 indicate that when selection probability is influenced by both X and Y, the coefficient of X (the parameter of interest) is biased. For the outcome model without interactions between X and Z (Scenario 5), where selection probability is affected by X, Y, and Z, a BCI is as follows. Bias correction index (BCI) = ${\alpha }_1+log\frac{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}}}{1+{\mathrm{e}}^{{\alpha }_0+{\alpha }_1+{\alpha }_2\mathrm{X}+{\alpha }_3\mathrm{Z}}}$. When external data are available to estimate the BCI, the selection bias may be corrected by including the BCI as an offset in the outcome model. However, when external data are not available, the conditional probability of selection P (S = 1 | Y, X, Z) with the selection parameters (α₀,  α₁,  α₂,  α₃) is unknown. To quantify the magnitude of bias for the odds ratio (OR) of interest, selection parameters (α) can be determined based on the researchers’ understanding of the participants. Alternatively, α can be set within a range of values. We performed sensitivity analyses by assuming selection probability model with different parameters (α₁, α₂) and presents the results using a graphical display. In the next section, we will illustrate how to construct a bias-correction sensitivity plot, which shows varying possible bias-corrected ORs under different magnitudes of selection parameters.

Application: analysis of survey data with potential selection bias

National Health Interview Survey (NHIS)

The NHIS in Taiwan involves a national representative sample and has been conducted every four years since 2001. The survey is conducted using a multistage stratified systematic sampling design, with sample selection based on probability proportional to size. Goodness-of-fit tests revealed no significant differences of urbanization level, sex, and age between the survey sample and household registration data (National Health Research Institutes & Health Promotion Administration (HPA) of the Ministry of Health and Welfare (MOHW) of Taiwan, 2019).

The NHIS database can be linked to Taiwan’s National Health Insurance Research database (NHIRD) to explore the associations between disease, health utilization outcomes, and health behaviors. The National Health Insurance is a mandatory program for all citizens in Taiwan with coverage rate over 99% (Hsieh et al., 2019). For NHIS data to be linked with the NHIRD, participants in the NHIS must provide written informed consent. The following section will illustrate the practice of a sensitivity analysis using the NHIS data of only those who consented to the data linkage.

The NHIS questionnaires and variables

We used the 2001 NHIS database to perform the sensitivity analysis because these data are often used to link to the NHIRD or cancer registries for long-term follow-up studies (Dai et al., 2020; Lo et al., 2022). This study was approved by the Institutional Review Board (IRB) of National Yang-Ming University (YM108041E).

Interviewers asked all participants through standard questionnaires and collected basic demographic data including age, sex, marital status and education level. Because smokers have a higher risk of colitis and colorectal cancer (Giovannucci, 2002; Srivastava et al., 1990), we hypothesized that nonsmokers have greater health consciousness and are more willing to undergo proctoscopic examination, which served as an example to introduce BCI. Participants were asked whether they had ever smoked in their lives (binary; X; 1 = yes, 0 = no), had undergone a proctoscopic examination in the past year (Y; 1 = yes, 0 = no), and had signed the informed consent for NHIRD linkage (S; 1 = yes, 0 = no).

The parameter of interest was the OR of the association between proctoscopic examination (Y) and smoking status (SMOKE) after adjusting for sex. We performed a sensitivity analysis assuming that data were only available for participants who consented to NHIRD linkage (S = 1). Details regarding variable information and model specifications are described in the following sections.

The main model can be presented as follows: logit P (Y = 1 | smoke, sex) = β₀ + β₁smoke + β₂sex.

Based on previous literature (Tolonen, Dobson & Kulathinal, 2005), health behavior may influence the participants’ willingness to participate in a study. We supposed that both proctoscopic examination (Y) and smoking status (SMOKE) were related to study participation (S = 1) and that sex was not. The selection model can be presented as follows: $\text{logit}\mathrm{P}\left(\mathrm{S}=1|\mathrm{Y},\text{smoke},\mathrm{sex}\right)={\alpha }_0+{\alpha }_1\mathrm{y}+{\alpha }_2\text{smoke}+{\alpha }_3\mathrm{sex},\text{where}{\alpha }_3=0$.

The bias-corrected model is rewritten as follows (Scenario 4): $\text{logit}\mathrm{P}\left(\mathrm{Y}=1|\text{smoke},\mathrm{sex},\mathrm{S}=1\right)={\beta }_0+\left.\left(\widehat{{\alpha }_1}+log\frac{1+{\mathrm{e}}^{\widehat{{\alpha }_0}+\widehat{{\alpha }_2}\text{smoke}}}{1+{\mathrm{e}}^{\widehat{{\alpha }_0}+\widehat{{\alpha }_1}+\widehat{{\alpha }_2}\text{smoke}}}\right.\right)+{\beta }_1\text{smoke}+{\beta }_2\mathrm{sex}$, where BCI = $\widehat{{\alpha }_1}+log\frac{1+{\mathrm{e}}^{\widehat{{\alpha }_0}+\widehat{{\alpha }_2}\text{smoke}}}{1+{\mathrm{e}}^{\widehat{{\alpha }_0}+\widehat{{\alpha }_1}+\widehat{{\alpha }_2}\text{smoke}}}$ .

In practice, researchers may determine the range of α based on their prior understanding of the relationship between S and other variables. A range of positive and negative α values are included if no prior information is available. The value of α₀ is set as 1.96 on the basis of the data, which indicates that marginal proportion of individuals that consented to linkage P(S = 1) was approximately 88%. The bias parameters (α₁, α₂) were set to between −2.5 and 2.5, with an increment of 0.3 for each value. The BCIs were calculated using varying combinations of (α₁, α₂). Then, the bias-corrected ORs of SMOKE and Y were calculated using different BCIs. A bias-correction sensitivity plot was performed by BCI to understand how selection bias influenced the estimate of the variable of interest (exposure-outcome relationship).

In addition to the sensitivity analysis, we also calculated the OR without selection bias (i.e., the OR based on the entire data) for demonstration, so that readers can observe the influence of varying BCIs. In practice, the OR without selection bias is unknown. However, the data of individuals who did not provide informed consent for NHIRD linkage (S = 0) are available. Therefore, in the present study, we were able to demonstrate the difference between the target estimate based on the entire data obtained from the main model and that based on the selected sample obtained from the outcome model in the bias-correction sensitivity plot.

All analyses were conducted two-sided with alpha set at 0.05. We conducted all statistical analyses and generated all charts by using lme4 (Bates et al., 2015), jtools (Long, 2020) and ggplot2 (Wickham, 2016) packages of R (version 4.0.0; R Core Team, 2020).