A cautionary note on the use of attributable fractions in cohort studies

1 Introduction

The attributable fraction (AF; also known as the attributable risk) is a widely used measure to quantify the public health impact of an exposure on an outcome. It was originally defined for binary outcomes as the proportion of outcomes that would be eliminated in the population if the exposure was hypothetically eliminated. ¹ In observational studies, consistent estimation of the AF typically requires appropriate confounder adjustment. Benichou^2,3 gave an overview of various adjusted estimators, including the Mantel–Haenszel estimator, the weighted-sum estimator, and the model-based estimator. The latter is based on the following algebraic rearrangement of the AF

AF = 1 - ∑ x , z ρ xz RR x | z

Benichou^2,3 claimed that the model-based estimator can be adapted to cohort studies with right censored time-to-event outcomes, by replacing RR x | z with an estimate of the adjusted (hazard) rate ratio obtained from a Poisson or Cox regression. However, he provided no formal justification for this claim. Furthermore, he did not clarify what to replace ρ xz with when the outcome is no longer binary. We made a survey of recent literature and found that the cohort adaption of the model-based estimator proposed by Benichou^2,3 has been used in many studies.^4–11 In all these studies, ρ xz was replaced with the sample proportion with levels x and z among those who were observed to experience the event before end of follow-up; we refer to the resulting estimator of the AF as the “cohort-adapted estimator.” Without exception, the authors of these studies appear to have interpreted the cohort-adapted estimator as estimating the proportion of events that would have been eliminated in the cohort during follow-up, if the exposure had been eliminated; we refer to this proportion as the “cohort attributable fraction” (CAF). In this note, we show that the cohort-adapted estimator only estimates the CAF under the following two assumptions: (1) the cumulative incidence is small during the whole follow-up, and (2) there is no censoring, apart from administrative censoring at the end of follow-up. When either of these assumptions is violated, the cohort adapted estimator does not estimate the CAF and may have no meaningful interpretation. We further argue that even if both these assumptions hold, the CAF may not be of substantial scientific interest. We propose a potentially more relevant measure of attributable fraction in cohort studies, the baseline attributable fraction (BAF). We show how the BAF can be conveniently estimated in cohort studies, using Cox proportional hazards models. ¹² Our estimator does not require the two aforementioned assumptions required by the cohort-adapted estimator.

The paper is organized as follows. In Section 2, we show that the cohort-adapted estimator only converges to the CAF under the aforementioned two assumptions. In Section 3, we argue that the CAF may not be a relevant measure of attributable fraction in the first place, and we propose the BAF as a potentially more relevant measure. In Section 4, we show how the BAF can be conveniently estimated in Cox proportional hazards models. In Section 5, we carry out a small simulation study to demonstrate the consistency of our proposed estimator. In Section 6, we present an application to real data.

2 The underlying assumptions of the cohort-adapted estimator

We first introduce some notation. We consider a binary exposure X, with levels 0 and 1 for “unexposed” and “exposed,” respectively. We let T denote the random time to event. We assume that the event of interest is “absorbing,” so that when the event occurs the individual is no longer part of the cohort, e.g. when the subject dies. We let τ denote the end of follow-up. To begin with, we assume that there is no censoring apart from administrative censoring at the end of follow-up.

We will throughout assume that the conditional log hazard function, given X and adjustment covariates Z, follows the Cox proportional hazards (PH) model

λ ( t | X , Z ) = μ ( t ) exp ( β X + γ Z )

(1)

1 - ∑ x , z Pr ( X = x , Z = z | T ≤ τ ) exp ( - β x ) = Pr ( X = 1 | T ≤ τ ) { 1 - exp ( - β ) }

To relate the cohort-adapted estimator to the CAF, we need a formal definition of the latter. Using standard potential outcome notation^13,14 we let T₀ denote the counterfactual time to event that would have been observed for a given subject, had that subject been unexposed. We let F ( t ) = Pr ( T ≤ t ) denote the factual proportion of events that occurs at or before time t in the cohort (i.e. the cumulative distribution function at T = t ), and we let F 0 ( t ) = P ( T 0 ≤ t ) denote the counterfactual proportion of events that would occur at or before time t in the cohort, had everybody in the cohort been unexposed. In this notation, the CAF is equal to

CAF = 1 - F 0 ( τ ) F ( τ )

In randomized trials, exposed and unexposed are exchangeable so that F 0 ( τ ) is equal to the proportion of subject who would experience the event before end of follow-up among those factually unexposed, i.e. F 0 ( τ ) = F ( τ | X = 0 ) . In observational studies, F 0 ( τ ) can be estimated by adjusting for covariates, if these are appropriately selected and sufficient for confounding control. To determine whether a particular set of covariates is sufficient or not is a difficult problem, and one that typically requires strong subject matter knowledge. This problem of covariate selection is beyond the scope of this paper; we proceed by assuming that we have measured a set of covariates Z which is sufficient for confounding control. Under this assumption, F 0 ( τ ) is equal to E Z { F ( τ | X = 0 , Z ) } , where the outer expectation is taken over the marginal distribution of Z.^15,16 It follows that the CAF can be rewritten as

1 - E Z { F ( τ | X = 0 , Z ) } F ( τ ) = 1 - E X , Z { F ( τ | X = 0 , Z ) F ( τ ) F ( τ | X , Z ) F ( τ | X , Z ) } = 1 - E X , Z { F ( τ | X = 0 , Z ) F ( τ | X , Z ) | T ≤ τ } = 1 - E X , Z { 1 - exp { - ∫ 0 τ λ ( u | X = 0 , Z ) d u } 1 - exp { - ∫ 0 τ λ ( u | X , Z ) d u } | T ≤ τ } ≈ 1 - E X { exp ( - β X ) | T ≤ τ } = Pr ( X = 1 | T ≤ τ ) { 1 - exp ( - β ) }

To see that the cohort-adapted estimator can be severely biased if the cumulative incidence is large, consider the case when the event of interest is death and the follow-up is infinite; τ = ∞ . In this case, the CAF is equal to 0, since F 0 ( τ ) = F ( τ ) = 1 (eventually, everyone dies). However, the cohort-adapted estimator equals Pr ( X = 1 | T ≤ ∞ ) { 1 - exp ( - β ) } = Pr ( X = 1 ) { 1 - exp ( - β ) } , which may be far from 0, unless β = 0 (i.e. unless the exposure has no effect).

We next consider the impact of random right censoring. Let C denote the random time to right censoring, so that for each subject we observe the event if T < min ( C , τ ) . ¹⁷ In the presence of random censoring, we cannot observe Pr ( X = 1 | T ≤ τ ) , but have to replace this proportion in the cohort-adapted estimator with Pr { X = 1 | T ≤ min ( C , τ ) } , that is, with the proportion in the cohort who are exposed among those who are observed to experience the event during follow-up. To see that this modification may induce bias, suppose that the censoring is heavy, so that a large fraction is censored early during follow-up. Those who contribute to Pr { X = 1 | T ≤ min ( C , τ ) } will then mainly be those who die early. If there is an association between the exposure and the time to death, then there will be an overrepresentation of exposed subjects in this group, as compared to those who die before end of follow-up. This implies that Pr { X = 1 | T ≤ min ( C , τ ) } > Pr ( X = 1 | T ≤ τ ) , so that the cohort-adapted estimator becomes biased upwards. We note that this bias will be present even if C is completely independent of both X, Z, and T.

3 The baseline attributable fraction

Many cohort studies have relatively long follow-up, or relatively high baseline hazard, so that the assumption of small cumulative incidence is not tenable. Censoring, either random or nonrandom, is also commonly present. Thus, the cohort-adapted estimator may often be biased, as an estimator of the CAF. Recently, methods have been proposed which can be used to consistently estimate the CAF even if the cumulative incidence is not small and random censoring is present.^18,19 However, it can be questioned whether the CAF is a relevant measure of attributable fraction in the first place. The problem is that the CAF depends on the actual follow-up time τ in a rather arbitrary fashion. In particular, as already noted, the CAF is 0 when the event of interest is death and τ is sufficiently large, since eventually everyone dies regardless of whether we are able to eliminate the exposure or not. To bypass this problem, one could instead present (an estimate of) the whole AF function

AF ( t ) = 1 - F 0 ( t ) F ( t )

We will focus on the limit value lim t → 0 AF ( t ) , which measures the proportion of events that would have been eliminated during a short (infinitesimal) period of time immediately after baseline in the cohort, if the exposure had been eliminated; we refer to this measure as the “baseline attributable fraction” (BAF). The BAF coincides with the maximum value of the AF function when the AF function is monotonically decreasing. An attractive feature of the BAF is that it does not only apply to the cohort but also to the population from which the cohort was drawn. This is not generally true for other summaries of the AF function, since the characteristics of a (closed) cohort often change and diverge from the characteristics of the source population over time, as subjects gradually drop out from the cohort due to death or for other reasons. At baseline though, the cohort is still guaranteed to be representative for the source population, and thus the BAF can be interpreted as the proportion of events that would have been eliminated during a short period of time in the source population, if the exposure had been eliminated. We emphasize that the BAF may not always be a relevant measure of attributable fraction, e.g. when the exposure has no short-term effect, so that the BAF is equal to 0 even though the exposure may have a long-term effect. However, in many scenarios, in particular when the effect of the exposure is relatively constant over time, the BAF summary measure may offer a useful compromise between relevance and parsimony.

Under Cox PH models, the BAF has a simple and instructive expression. Consider the following Cox PH model, which generalizes the model in (1) in that it allows for both X and Z to vary over time

λ { t | X ( t ) , Z ( t ) } = μ ( t ) exp { β X ( t ) + γ Z ( t ) }

(3)

BAF = 1 - exp ( - β ) 1 + exp { - ( β + Ω + α 1 - α 0 ) }

(4)

4 Estimation of the BAF under Cox PH models

The Cox PH model in (3) can be fitted to data with standard software to obtain estimates β ^ and γ ^ . Subsequently, the log of the sample means of exp { γ ^ Z ( 0 ) | X ( 0 ) = 1 } and exp { γ ^ Z ( 0 ) | X ( 0 ) = 0 } can be used as estimates α 1 ^ and α 0 ^ , respectively. The sample log odds of being exposed at baseline can be used as an estimate Ω ^ . Finally, the estimates ( β ^ , Ω ^ , α 1 ^ , α 0 ^ ) can be plugged into (4), to obtain an estimate BAF ^ .

To assess the statistical uncertainty in BAF ^ , it is desirable to derive its asymptotic distribution. Using the multivariate delta method, ²⁰ we have that BAF ^ has an asymptotic normal distribution with variance var ( BAF ^ ) = g T Σ g , where g is the column vector of derivatives of BAF with respect to ( β , Ω , α 1 , α 0 ) , and Σ is the asymptotic variance–covariance matrix of ( β ^ , Ω ^ , α 1 ^ , α 0 ^ ) . By differentiating the BAF with respect to ( β , Ω , α 1 , α 0 ) , we have that

∂ BAF ∂ β = exp ( - β ) [ 1 + exp { - ( Ω + α 1 - α 0 ) } ] [ 1 + exp { - ( β + Ω + α 1 - α 0 ) } ] 2

∂ BAF ∂ Ω = ∂ BAF ∂ α 1 = - ∂ BAF ∂ α 0 = { 1 - exp ( - β ) } exp { - ( β + α 1 - α 0 + Ω ) } [ 1 + exp { - ( β + α 1 - α 0 + Ω ) } ] 2

U i ( θ ) = ( S i ( β , γ ) X i ( 0 ) - expit ( Ω ) X i ( 0 ) [ exp { γ Z i ( 0 ) } - exp { α 1 } ] { 1 - X i ( 0 ) } [ exp { γ Z i ( 0 ) } - exp { α 0 } ] )

The derivation of var ( BAF ^ ) outlined above assumes that the subjects are independent. Often, this assumption may be violated due to clustering, e.g. in studies of siblings or twins. When the clusters are independent, it is easy to correct var ( BAF ^ ) for dependencies within clusters. Let U ij ( θ ) denote the estimating equation contribution from subject j in cluster i, and define U i ( θ ) = ∑ j U ij ( θ ) . With this minor modification, the U i ( θ ) ’s are now independent and var ( BAF ^ ) is equal to g T ( h T Σ ∼ h ) θ * g as before. ²¹

5 Simulation study

To investigate the finite sample properties of the proposed methods, we carried out a small simulation study. In all simulated scenarios, we generated a time-stationary covariate Z from a normal distribution with mean 0 and standard deviation σ. We generated a time-stationary exposure X from the logistic model logit ( p ) = ω + Z . For each simulated scenario, ω was set to achieve a desired value of p, as described below. We generated survival times T from the Cox proportional hazards model in (1), with μ ( t ) = constant = 1 . We generated random censoring times C from the same conditional distribution as T, given X and Z. Under this scheme, C is associated with X and Z, but conditionally independent of T, given X and Z. The marginal (over X and Z) censoring rate is 50%. From this model, we generated 1000 samples of n = 500 independent subjects, for each of the eight possible parameter combinations of σ = ( 0 . 5 , 1 . 5 ) , p = ( 0 . 05 , 0 . 15 ) , and β = ( 0 . 5 , 1 . 5 ) . For each sample, we estimated the BAF and its standard error using the method described in Section 4. Table 1 displays the mean (over the 10,000 samples) estimated BAF ( BAF ^ ) together with the true value ( BAF ), and the mean estimated standard error ( SE ) together with the empirical standard error ( ESE ) for each scenario. We observe that the mean estimated BAF is close to the true value and that the mean estimated standard error agrees well with the empirical standard error, for all scenarios.

OPEN IN VIEWER

Table 1.

Simulation results for estimation of the BAF, for 1000 samples of n = 500 independent subjects.

σ	p	β	BAF	BAF ^	SE	ESE
0.5	0.05	0.5	0.04	0.04	0.03	0.03
0.5	0.05	1.5	0.18	0.18	0.06	0.07
0.5	0.15	0.5	0.11	0.11	0.04	0.05
0.5	0.15	1.5	0.39	0.39	0.07	0.06
1.5	0.05	0.5	0.14	0.14	0.08	0.10
1.5	0.05	1.5	0.46	0.47	0.10	0.13
1.5	0.15	0.5	0.24	0.23	0.07	0.09
1.5	0.15	1.5	0.63	0.62	0.06	0.08

6 Real data illustration

In this section, we present an application to real data, borrowed from Carlsson et al. ²³ These authors aimed to study the association between body mass index (BMI) and mortality. Their cohort data comprise 44,258 Swedish same-sex twins, of which 16,793 are monozygotic (MZ) and 27,465 are dizygotic (DZ). All twins filled in a questionnaire in 1967 or 1972 on lifestyle factors, health, height, and weight. Death records were obtained by linkage to the National Causes of Death Registry for the years 1972–2004. The authors stratified the analysis on zygosity and sex, and controlled for smoking in the statistical models.

We reanalyzed the Carlsson et al data as follows. Baseline (i.e. t = 0 ) was defined as the date when the questionnaire was filled in. BMI was dichotomized into BMI ≤ 25 (unexposed) and BMI > 25 (exposed). The BAF was estimated for the DZ and MZ twins separately, controlling for smoking, sex, and age at baseline. Clustered standard errors were computed to account for the correlation between twins in the same pair, as described in Section 4. These standard errors were subsequently used to construct 95% Wald confidence intervals (CIs). For the DZ twins we obtained BAF ^  = 0.06, with 95% CI ranging from 0.02 to 0.10. For the MZ twins, we obtained BAF ^  = 0.07, with 95% CI ranging from −0.11 to 0.25. These results indicate that around 6–7% of all deaths would have been eliminated in this cohort of DZ/MZ twins, immediately after baseline, if all twins would have had BMI ≤ 25 at baseline.

Before leaving this example, we take the opportunity to note that attributable fractions are conceptually difficult when the exposure of interest is BMI. The reason is that it may not be entirely clear what it means to be “unexposed,” e.g. to have BMI ≤ 25. Clearly, there are many ways to reduce BMI; we may for instance drastically reduce the amount of body fat or the amount of muscles on everybody. These two hypothetical interventions may both render everybody “unexposed,” but the mortality rates under these interventions may be very different. Thus, unless we make the exposure definition more precise, we may consider the attributable fraction as somewhat ill-defined. We therefore conclude by recommending that researchers who wish to study the public health impact of obesity should make an effort to make their definition of obesity/nonobesity as precise as possible, in order to reduce vagueness in the research question being posed.

7 Discussion

We have studied the commonly used cohort-adapted estimator, which is supposed to estimate the CAF. We have shown that this estimator relies on two assumptions, which are often not tenable in practice. We have argued that the CAF may not be a relevant measure of public health impact in the first place, and we have proposed to use the BAF instead. We have shown how the BAF can be conveniently estimated in Cox proportional hazards models.

A major obstacle for practitioners who wish to use new methodology is often the lack of implementation in standard software. To facilitate the proposed methods, we have written a convenient R-function, which implements the BAF estimator proposed in Section 4. This R-function is available at http://www.meb.ki.se/personal/arvsjo/, and can be obtained from the author upon request. We describe its usage in the Appendix 1.