Accounting for biases in survey-based estimates of population attributable fractions

Background This paper discusses best practices for estimating fractions of mortality attributable to health exposures in survey data that are biased by observed confounders and unobserved endogenous selection. Extant research has shown that estimates of population attributable fractions (PAF) from the formula using the proportion of deceased that is exposed (PAFpd) can attend to confounders, whereas the formula using the proportion of the entire sample exposed (PAFpe) is biased by confounders. Research has not explored how PAFpd and PAFpe equations perform when both confounding and selection bias are present. Methods We review equations for calculating PAF based on either the proportion of deceased (pd) or the proportion of the entire sample (pe) that receives the exposure. We explore how estimates from each equation are affected by confounding bias and selection bias using hypothetical data and real-world survey data from the National Health Interview Survey–Linked Mortality Files, 1987–2011. We examine the association between cigarette smoking and all-cause mortality risk in the US adult population as an example. Results We show that both PAFpd and PAFpe calculate the true PAF in the presence of confounding bias if one uses the “weighted-sum” approach. We further show that both the PAFpd and PAFpe calculate biased PAFs in the presence of collider bias, but that the bias is more severe in the PAFpd formula. Conclusion We recommend that researchers use the PAFpe formula with the weighted-sum approach when estimates of the exposure-outcome relationship are biased by endogenous selection.


Background
This paper discusses best practices for estimating the fraction of mortality attributable to health exposures in survey-based data that are biased by both observed confounders and unobserved endogenous selection. Much extant work has reviewed errors in computing population attributable fractions (PAFs) in the presence of confounders [1][2][3][4][5][6], but little work has considered how different formulae for computing PAFs are affected by endogenous selection biases (e.g., collider bias).
Endogenous selection bias can affect estimates of statistical associations in many ways. Conditioning on a collider variable-that is, a variable caused by two other variables that are associated with the exposure and the outcome-can occur through statistical control, stratification of the sample into different groups, or the selection of participants into a study [7][8][9][10][11]. Introducing collider variables through any of these mechanisms can bias estimates of associations between exposure and outcome. In this study, we focus on unobserved endogenous selection-a problem that commonly occurs in health studies through the sampling process of recruiting study participants. Simply put, the likelihood of participation in a health study can be affected by both the exposure and outcome, which can bias estimates of the true association between them.
The most common PAF formulae are based on either the proportion of deceased (pd) in the sample that receives the exposure or the proportion of the entire sample (pe) that receives the exposure [1]. The two main aims of this investigation are to examine the performance of these model-based methods for calculating PAF in the presence of (1) known and observable confounders of the exposure-mortality association and (2) collider bias. We focus on the association between cigarette smoking and all-cause mortality risk in the US adult population, which is confounded by other variables and also a likely contributor to unobserved endogenous selection bias in survey-based data of smoking and mortality risk [7].

Methods
We use hypothetical data and real-world survey data to calculate PAF in the presence of confounding and unobserved endogenous selection. In all of our exercises, non-exposed cases are respondents who have never smoked cigarettes and exposed cases are respondents who are current or former smokers. The association of interest is how smoking affects all-cause mortality risk. For each exercise, we estimate the fraction of US mortality attributable to cigarette smoking using the PAF pd formula: where pd is the prevalence of a health exposure among the deceased cases and RR is the mortality risk ratio between the exposed and non-exposed subjects [12]. We also estimate this fraction using the PAF pe formula: where pe is the prevalence of the exposure among all cases in the sample [6,13]. For each formula, we adopt a "weighted-sum" approach [1,3,14,15], which uses model-based adjusted estimators of PAF separately for each adjustment level i as well as the distribution of cases by the adjustment levels: where i indicates the adjustment level (i.e., confounder) and W i indicates the proportion of deaths in adjustment level i. The weighted-sum approach is mathematically equivalent to the PAF pd [1,14]; combining the PAF pd with the weighted-sum approach is therefore redundant. Nevertheless, we apply it in all of our exercises to maintain consistency.

Exercise 1: Observed confounding bias in hypothetical data
In our first exercise, we examine PAF estimates from Eqs. (1)(2)(3)(4) in the presence of a single confounder, race/ ethnicity. For simplicity, we consider race/ethnicity using only two categories, non-Hispanic black and non-Hispanic white (hereafter black and white). The hypothetical data are composed of 1000 black respondents and 4000 white respondents. Both smoking prevalence (i.e., pe) and mortality risk are higher among black respondents than among white respondents, which confound the smoking-mortality association. In these data, black pe is 0.35 compared to white pe of 0.2, and overall mortality risk for black respondents is 0.3 compared to 0.2 for white respondents.

Exercise 2: Unobserved endogenous selection bias in hypothetical data
Our second example uses the same data as before, but presupposes that estimates of the smoking-mortality association are biased by differential selection into the sample. We assume that current smokers sampled are relatively more select on health than are non-smokers. That is, both the non-smoking and the smoking samples are healthier than the true populations, but the difference between the smoking sample and the smoking population is greater than the difference between the non-smoking sample and the non-smoking population. This unobserved process of health selection biases downward the all-cause mortality RR estimated in the sample data. When these conditions hold, both PAF estimates will be biased due to the central role of the RRs (see Eqs. [1][2][3][4]. Moreover, the distribution of deaths by exposure and by adjustment levels, W i , will also be biased. This is because counts of deaths among the exposure group in the sample will be artificially low and, consequently, W i will be incorrect. Thus, PAF pd and PAF pe estimates will remain biased via W i even if our adjusted RRs account for collider bias. Finally, the estimated PAF from the PAF pd formula will be additionally biased, due to the central role of the pd in the calculation of the PAF. That is, the pd in the observed data, like the RRs in the observed sample data, will be downwardly biased because deaths among the smoking sample are underreported.

Exercise 3: PAF estimation with real-world survey data
Finally, we illustrate the points above by analyzing the smoking-mortality association in the National Health Interview Survey-Linked Mortality Files (NHIS-LMF) for years 1987-2009. These data are composed of NHIS waves from 1987 and 1989-2009 that have been linked to official death records at the National Death Index through December 31, 2011 (the 1988 NHIS survey did not contain information about respondents' smoking behavior). The NHIS-LMF are designed to form a representative sample of non-institutionalized US adults [12]. To simplify the example, we limit the analytic sample to contain only US adult black and white men and women aged 40 through 84 at time of interview and whose survival is followed between ages 50 and 84. We extend the example by considering two levels of smoking exposure, "former smoker" and "current smoker," and by considering three possible confounders of the smoking-mortality association: race/ethnicity (i.e., white and black), gender (i.e., men and women), and age group (i.e., 50-59, 60-69, 70-79, and 80-84).
We fit a series of clog-log discrete-time survival models to estimate smoking-based differences in US adult mortality risk. First, we fit a baseline model that estimates differences in mortality risks between current, former, and never smokers (reference category). Next, we fit a confounder model that estimates age-specific differences in mortality risks between current, former, and never smokers, adjusting for race/ethnicity and gender as categorical confounders of the smoking-mortality association. We also fit models separately for black and white men and women that estimate age-specific RRs for former and current smokers compared to never smokers (i.e., confounder-specific models to be used with the weighted-sum approach to calculate PAFs). Finally, we fit a bias model that refits the confounder model by accounting for cohort-based variation in mortality risk and age-related selection biases in the NHIS-LMF data.
Participants in health surveys like the NHIS are positively selected on survival, health, and non-institutional living arrangements [16]. These selection biases tend to grow stronger with increasing age [17]. Thus, older respondents in NHIS-LMF data are selected on the outcome of interest (i.e., survival) and inclusion in the NHIS sampling frame (i.e., healthy and living in non-institutionalized housing). Combined, the selective nature of the sample results in collider biases via age-related selection into the sampling frame and the selective factors associated with age are likely stronger among respondents with health risk factors such as smoking than among healthy respondents [8].
Survival models fitted separately by cohort of entry into the NHIS sample provide evidence consistent with these assumptions about collider biases. For example, the estimated RR between current smokers and never smokers who died at age 70-80 ranges from 1.51 [1.43-1.58 95%CI] among respondents surveyed at age 70-75 to 4.11 [3.02-5.57 95%CI] among respondents surveyed at age 50-55. The bias model is a shared frailty survival model that estimates random effects variation in mortality risk by NHIS respondents' 5-year age cohorts at the time of sampling. Overall, the model fits age-specific mortality risks separately for current, former, and never smokers, adjusting for gender, race/ethnicity, birth year, and random effects for a 5-year cohort of entry into the data. Mortality differences between US adults self-reported to be current, former, and never smokers between ages 50 and 84 are estimated across these three models. We use the adjusted RRs between (1) current smokers and never smokers and (2) former smokers and never smokers, which are estimated from confounder-specific survival models and the weighted-sum approach to calculate the PAF for smoking as a cause of death in the US adult black and white populations between ages 50 and 84 for years 1987-2011. For all models, we contrast PAFs calculated from PAF pd with PAFs calculated from PAF pe to examine how each formula is affected by (1) confounders in the estimated smoking-mortality association and (2) collider bias.
The NHIS-LMF data analyzed for the current study are public-use files made available by the NCHS (https:// www.cdc.gov/nchs/data-linkage/mortality.htm). The analytic scripts (Additional file 1) and calculations to generate results (Additional file 2) for Exercise 3 are available in the appendix.

Exercise 1: Observed confounding in hypothetical data
The confounding effect of race/ethnicity on the smokingmortality association is illustrated in Table 1. The allcause mortality RR for smoking when unadjusted for the confounding effects of race/ethnicity is (450/1150)/(650/ 3850) = 2.32. Alternatively, the RR adjusted for race/ethnicity is 2.23. That is, when we estimate separate RRs for each race/ethnicity sample, we observe When these race/ethnic-specific RRs for smoking are standardized by the race/ethnic distribution of deaths and the race/ethnic distribution of smoking prevalence, the adjusted RR is 2.23. If one does not account for the confounding effects of race/ethnicity on both mortality risk and the probability of smoking, one would incorrectly estimate the PAF by the following: As a result, estimates of the PAF for smoking, irrespective of the formula used, would be biased by not attending to the confounding effects of race/ethnicity: The actual PAF shown in the counterfactual example above is (1100 − 856)/1100 = 0.222 Thus, by failing to account for (1) the higher prevalence of smoking among black respondents and (2) the higher mortality risks among black respondents, we would incorrectly inflate the RR associated with smoking and misattribute numerous deaths to smoking as a cause of mortality in the population. As such, it is necessary to identify the RR by accounting for confounders in model estimates, and then use this confounder-adjusted RR to calculate PAFs [18]. It has been argued that only the PAF pd formula can accurately estimate the PAF when using confounder-adjusted RRs [3][4][5][6]. Yet, as others have noted, one can use the PAF pe equation with the confounder-adjusted RR to derive the true PAF [1,15]. To do so, one needs to first estimate separate PAFs for each confounder group (i.e., each adjustment level i), and then standardize these confounder-specific PAFs by the distribution of deaths across groups (i.e., W i ).
To illustrate, when we estimate separate PAFs for black and white respondents, we see for black: and for white: To estimate the total PAF, we further attend to the distribution of deaths across groups. That is, we simply weight the confounder-specific PAFs by the proportion of total deaths occurring in the confounder groups (i.e., W i ) [18]. The proportion of the total deaths that occurred among black respondents = (300/ 1100) = 0.273 and the proportion of total deaths that occurred among white respondents = (800/1100) = 0.727. When we weight the confounder-specific PAFs by the proportion of deaths in the two groups, W i , we retrieve the true overall PAF: This shows that the weighted-sum approach can calculate the true PAF regardless if one uses the PAF pd or PAF pe formula. So long as (1) unobservable confounders or unobservable selection do not induce bias, and (2) one attends to observable confounders of the smokingmortality association, one can use adjusted RRs with either PAF pd or PAF pe and the weighted-sum approach to calculate the PAF for smoking-related mortality in the sample [18].

Exercise 2: Unobserved endogenous selection in hypothetical data
In the next exercise, we extend the previous example to consider sample data that are biased by unobserved selection, causing underestimation of mortality risk in the smoking population. To simplify matters, let us assume that the prevalence of smoking is the same in both the sample and population so that the only change pertains to q x for smokers in the sample. The new information about population parameters is presented in Table 2 below.
The mortality probabilities for the non-smoking populations equal those in the sample data (0.231 among blacks and 0.156 among whites). Smoking prevalence is also the same (pe NHB = 0.35 and pe NHW = 0.20). However, we now see discrepancies in the mortality risks for the smoking populations (0.500 in the white population vs. 0.375 in the white sample, and 0.500 in the black population vs. 0.429 in the black sample). These, in turn, affect the RRs for smoking (e.g., 2.17 vs. 1.86 for black and 3.21 vs. 2.40 for white), the pds (0.538 vs. 0.500 for black and 0.444 vs. 0.375 for white), and the W i (e.g., 0.265 of population deaths are among blacks vs. 0.273 of sample deaths).
The confounder-specific PAFs using both the PAF pd and PAF pe formulae are as follows (estimates might be slightly different due to rounding): non-Hispanic black: non-Hispanic white: Standardizing these confounder-specific PAFs by the distribution of deaths, W i , we use the weighted-sum approach to calculate the true PAF: We see that the PAFs in the sample data underestimate the true PAF in the population (0.222 vs. 0.301), and this bias is the same in the PAF pd and PAF pe formulae. The discrepancy arises from one's inattention to (unobservable) endogenous selection bias in the sample data, resulting in biased sample estimates of the mortality RRs associated with smoking as well as biased W i in the sample.
Imagine that we had accounted for unobservable selection bias in our survival models and correctly identified the RRs for smoking for both the black and white samples. Even though the adjusted RRs would be correct in our survival models, the counts of deaths in the sample data would remain biased. Consequently, the pd values in the sample stay at 0.50 and 0.375, and the proportion of deaths occurring among blacks and whites stay at 0.273 and 0.727, respectively. As a result, if we were to calculate the PAF using the adjusted RRs with PAF pd , we would find The confounder-specific PAFs are biased (i.e., 0.270 estimated vs. 0.290 actual for blacks and 0.258 estimated vs. 0.306 actual for whites) even when using the adjusted RRs. Furthermore, when we use the weighted-sum approach and standardize these PAFs by W i , we add another source of bias because the distribution of deaths in each confounder group is biased as well: total PAF = (0.270*0.273) + (0.258*0.727) = 0.263. Yet, were we to follow conventional wisdom [3][4][5][6] and use the adjusted pe proportion smoker in population, pd proportion smoker among deceased in population, q x (NS) probability of death among nonsmokers in population, q x (S) probability of death among smokers in population, RR risk ratio RR with the PAF pd for the entire sample, we would estimate the same biased PAF: Thus, even if we accurately accounted for selection bias in our survival models and estimated an unbiased RR (e.g., by fitting frailty models that account for selection bias in the smoking RR [19]), the PAF calculated from the PAF pd formula will still be biased. In this case, a biased 0.263 is estimated for the sample when the true PAF in the population is 0.301 (a bias on the proportionate scale of 12.6%: (0.263 − 0.301)/0.301).
If we calculate the PAF using the confounder-and selection-adjusted RRs with the PAF pe formula, we find We see that the confounder-specific PAFs are unbiased. Only when we standardize these PAFs by the distribution of deaths, W i , do we introduce slight bias in the total PAF = (0.290*0.273) + (0.306*0.727) = 0.302 (a bias on the proportionate scale of − 0.3%: (0.301 − 0.302)/0.302). Thus, when we account for selection bias in our survival models and estimate unbiased adjusted RRs, the PAF calculated from PAF pe will be biased, but only via W i . By using the PAF pe equation, we avoid bias in estimates from the pd and dramatically reduce the overall bias in the PAF estimate (0.3% vs. 12.6%).
To recap, when sample data are biased by unobserved selection, both the PAF pd formula and the PAF pe formula will calculate a biased PAF-even if researchers adjust for selection bias in the data. However, the PAF pd formula is far more affected by the bias than is the PAF pe formula because bias is introduced in both the pd and W i . Conversely, estimates of the confounder-specific PAF from the PAF pe equation are not biased, but some bias is introduced in the weighted-sum approach via W i . Theoretically, one could completely eliminate bias by identifying the true RR (i.e., attend to both observable confounders and unobservable selection biases) and standardizing the PAFs by the true distribution of deaths for each adjustment level (i.e., use population data to estimate W i ).

Exercise 3: PAF estimation with real-world survey data
For the final exercise, we calculate PAF for smoking as a cause of US adult mortality in the NHIS-LMF data, which are biased by confounding (i.e., age, race/ethnicity, and gender) and likely biased by endogenous selection (i.e., likelihood of sample inclusion depends on health). Table 3 shows age-specific mortality risks between years 1987 and 2011 for NHIS respondents who are current, former, and never smokers. The pd for former smokers (0.352) combined with the pd for current smokers (0.338) indicates that nearly 70% of the deceased NHIS sample had been exposed to smoking.
From the sample data in Table 3, we calculate the unadjusted RRs: Because we are calculating a PAF for two-levels of an exposure, former smokers and current smokers, the PAF formulae change slightly [18,20]: We see that if we did not consider age, race/ethnicity, or gender as confounders of the smoking-mortality association in these NHIS-LMF data, we would estimate about 23% of US black and white adult deaths between ages 50 and 85 for years 1987-2011 were attributable to cigarette smoking.
Average RRs for current smoking estimated from cloglog discrete time hazard models are presented in Table 4, and overall PAFs estimated from the PAF pe and PAF pd formula are included as well.
The baseline model estimates mortality risks for former and current smokers relative to never smokers that match the RRs observed in Table 3 (i.e., 1.49 and 1.52, respectively). Using these RRs, we estimate the same 0.232 PAF for smoking as a cause of US adult mortality, regardless if we estimate the PAF from the PAF pe formula or the PAF pd formula. The confounder model estimates age-specific RRs for former and current smokers relative to never smokers while controlling for confounding by gender and race/ethnicity. The age patterns in the RRs for current smokers suggest that the mortality consequences of smoking significantly decline with age. For example, current smokers are estimated to have about 2.6 to 2.7 times the mortality risk as never smokers in age-groups 50-59 and 60-69, but only about 1.2 times the mortality risk in age-group 80-84. When using these confounder-adjusted and age-specific RRs for smoking, we estimate a 0.247 PAF for smoking as a cause of US adult mortality.
Finally, the estimated age-specific RRs from the bias model are significantly larger than the age-specific RRs from the confounder model, especially at older ages. Although the smoking-mortality relationship attenuates with age, it is substantially less than the attenuation observed in the confounder model. Using these confounder-and selection-adjusted RRs, we calculate a PAF of 0.289 from the PAF pd formula and a PAF of 0.326 from the PAF pe formula. This is the only case in which we observe different PAF values depending on the formula used. This is because the PAF pd formula remains biased by pd and likely underestimates the amount of mortality attributable to cigarette smoking in the US adult population. In this case, the PAF estimated from the PAF pd formula is likely additionally biased by − 11.3% over the PAF pe (0.289 − 0.326)/0.326) because it does not fully account for collider bias in estimates of the smoking-mortality association in the NHIS-LMF data.

Discussion
Between-group differences in mortality (e.g., smokers and non-smokers) estimated from survey data are often biased by unobserved endogenous selection [8,10]. These biases can distort research findings and lead to incorrect conclusions and misguided policy recommendations. Researchers should therefore be wary of collider biases and, when possible, adjust estimates to account for them. Relatedly, researchers should be wary of how these biases affect PAF calculations. In this paper, we demonstrated that the PAF pd formula is far more sensitive to collider bias than the PAF pe formula. Results from both our hypothetical examples and real-world illustration using the NHIS-LMF show the PAF pd formula calculated severely biased estimates of the PAF for smoking as a cause of mortality. As such, if estimates of the exposure-outcome association are likely biased by endogenous selection, researchers should consider calculating PAFs using the PAF pe formula with the weighted-sum approach. The main challenge to using the weighted-sum approach is the data required to scale estimates by W i , which increase with the number of confounders in the model. In addition, the weighted-sum approach may not be appropriate in small samples because estimates of W i are unreliable [1].
The findings are important for researchers aiming to estimate the mortality burden of exposures that may induce collider bias in sample data. For example, estimates from the NHIS-LMF data indicate that widening educational disparities in US adult mortality have greatly increased deaths attributable to low educational attainment [21]. Yet, estimates of the education-mortality association in the NHIS-LMF data may be biased by mortality and health selection across age [22]. Deaths attributable to low education in the USA may, in fact, be underestimated by not accounting for collider bias in PAF calculations. Also, researchers have reported discrepant PAFs for obesity as a cause of US mortality. For example, Flegal et al. [5] review PAF values indicating 2-15% of adult deaths are attributable to high BMI. The discrepancies likely reflect the extent to which researchers attend to confounder and collider biases in model estimates and how these biases affect PAF calculations. While Flegal et al. ([5] p. 203) consider the PAF pe to be "the invalid formula" and PAF pd to be the "formula appropriate for use with adjusted relative risks when confounding exists," their review did not consider how the PAF formulae were affected by collider bias. Results here indicate that the PAF pd is, in fact, the formula that calculates more biased estimates when relative risks are adjusted for confounding and selection biases.

Conclusion
Many studies have addressed best practices for calculating and interpreting PAFs for causes of mortality [1,3,5,6,20,[23][24][25]. In this paper, we extend these discussions to consider how unobserved endogenous selection bias (e.g., collider bias) distorts calculations of PAFs in the PAF pd and PAF pe formulae. Prior research has highlighted the importance of confounding bias in PAF calculations, but it has not considered how collider bias may affect PAF calculations. We used both hypothetical and real-world data on the smoking-mortality relationship to explore these considerations. Results from our examples demonstrate that both the PAF pd and PAF pe formulae can equally attend to observable confounders and accurately calculate PAFs via the weighted-sum approach [1,3,18]. Yet, the PAF pe formula via the weighted-sum approach is preferred to the PAF pd formula if RR estimates for the exposure are biased from endogenous selection. In contrast to conventional wisdom that recommends using the PAF pd formula with adjusted RRs [3,5,6], we conclude by recommending the use of the PAF pe formula with the weighted-sum approach when using RRs adjusted for both confounding bias and selection bias.