Risk Estimation with Epidemiologic Data When Response Attenuates at High-Exposure Levels

Briefly, we previously conducted a cohort analysis, with a focus on breast cancer, of 7,576 women exposed to EtO while sterilizing medical supplies (see Steenland et al. 2003 for details of this study). Here we focus on the 5,138 women from that cohort who we were able to interview; of these women, 232 had been diagnosed with breast cancer. Estimated daily exposures to EtO (average intensity) across different jobs and time periods ranged from 0.05 ppm to 77 ppm—the current U.S. Occupational Safety and Health Administration standard is 1 ppm (U.S. Department of Labor 2007), with an average duration of exposure of 11 years and a median of 7 years (range, 1–50 years). Cumulative exposure was calculated by multiplying duration of exposure by intensity of exposure for each job or time period and summing across different jobs and time periods, with a 15-year lag period such that exposure within 15 years of case diagnosis (or the corresponding date in noncases) was excluded from analyses. The 15-year lag was optimal (best fitting) in log-linear models (Steenland et al. 2003), and we used it again here for linear RR models for consistency.

Models were fit using Cox regression, in which each case was compared with its matched risk set (with matching based on white or nonwhite race), cumulative exposure was time dependent, and the time variable was age. Analyses included the following potential confounders: parity (continuous), date of birth, and family history of breast cancer in first-degree relatives (yes/no). In the linear RR models, one has the choice of combining additive terms for the main exposure with multiplicative terms for covariates, for example, RR = (1 + β₁X₁) × exp(β₂X₂) × exp(β₃X₃), where X₂ and X₃ represent covariates, or of using a strictly additive model for both covariates and exposure (e.g., RR = 1 + β₁X₁ + β₂X₂ + β₃X₃). We present results for the strictly additive model, because an inspection of model likelihoods indicated that this type of linear RR model fit the data best. SAS NLP (version 9.1; SAS Institute Inc., Cary, NC) was used for analyses of both log RR models and linear RR models, as suggested by Langholz and Richardson (2010). As noted, the SAS procedure provides comparable model likelihoods for both types of models.

We present a variety of log RR and linear RR models, as described below, including a categorical exposure model (deciles), a linear (untransformed) exposure model (e.g., RR = 1 + cumexp), a model using the log transformation of exposure [(e.g., RR = 1 + log(cumexp)], a model using a square-root transformation (e.g., RR = 1 + square root (cumexp), a two-piece spline model, a linear-exponential model (used only in the linear RR case), and a Michaelis–Menten model. The square-root transformation fits curves that attenuate at high exposures and was chosen to represent the class of polynomial models.

We restricted ourselves to relatively simple statistical models and did not use any biologically based models (e.g., two-stage models), which are not commonly used for risk assessment using human data.

In categorical analyses, exposure was classified into 11 groups, with the reference group including cases and noncases who had no EtO exposure > 15 years before the case was diagnosed with breast cancer (62 breast cancers), and 10 groups with cumulative exposures > 0 with cut points selected such that an approximately equal number of breast cancer cases (n = 17) occurred in each decile: 1–355, 356–842, 843–1,361, 1,362–2,187, 2,188–3,772, 3,773–5,522, 5,523–7,891, 7,892–14,483, 14,484–25,112, and > 25,112 ppm-days.

The linear spline model had the form (linear RR version) of RR = 1 + β₁ × cumexp + β_2z × maximum (0, cumexp – knot), where the value corresponding to maximum (0, cumexp – knot) is either the cumulative exposure value less the exposure level used to define the knot, or 0 if the cumulative exposure is less than the value of the knot. The slope of the first exposure–response segment based on the model is represented by β₁, whereas the slope of the second segment is β₁ – β₂. Knots for the two-piece spline model were chosen to maximize the log likelihood after trying knots at 100 ppm-day intervals. In our case, there was a clear optimal knot at 5,800 ppm-days (linear RR model). [See Supplemental Material, Figure 1 (doi:10.1289/ehp.1002521).]

The linear-exponential model had the form RR = 1 + (β₁ × cumexp) × exp(β₂× cumexp). This model has been used in radiation studies (Schneider and Walsh 2008) as well as in modeling other human carcinogens (e.g., arsenic) (Lubin et al. 2008) and allows for the linear exposure–response in the usual linear RR model to be modified to exhibit various degrees of curvature.

We also fit the Michaelis–Menten model of enzyme kinetics, in which the rate of reaction tails off as substrate concentration increases (Lehninger 1982). The model has been used rarely in epidemiology, but it is simple to implement, because it merely involves the transformation of exposure and the estimation of one additional parameter to describe the exposure–response relationship. For the linear RR model, the model is RR = 1 + β₁ × [cumexp/(cumexp + k)], where k represents a constant, analogous to the Km used in biochemistry for enzyme-substrate kinetics, and is estimated from the data.

Another option, not often used in etiologic research but used in studies of the effects of radiation (e.g., Little 2008), is not an RR model at all, but an excess absolute risk (EAR) model, also referred to as an additive Poisson model. This model has the form R_exp = R_nonexp + β × cumexp. We fit this model following the presentation of Boshuizen and Feskens (2010), who present SAS code for the GENMOD procedure to fit this model via Poisson regression. Poisson regression requires categorized data for rates, unlike the Cox regression models used to fit the linear RR and log RR models here. We cross-categorized events and person-time by age, cumulative exposure (lagged 15 years), and calendar time for these analyses and compared the EAR model with the log RR model. We also used a scale-parameter statistic that adjusted for over- or underdispersion.

Excess risk for environmental exposure to EtO according to each model was calculated using life table techniques accounting for all-cause mortality and applying an adaptation of the method described in the 1988 Biological Effects of Ionizing Radiation (BEIR) report (BEIR 1988). Because the exposure of interest was environmental (i.e., risk to the general public due to EtO in ambient air) and was assumed to begin at birth but with a 15-year lag, risk was effectively considered to begin at age 15 years. More details on the excess risk calculation can be found in Supplemental Material (doi:10.1289/ehp.1002521) and also in Appendix C of the 2006 draft U.S. EPA cancer risk assessment for EtO (U.S. EPA 2006). Briefly, background mortality rates from U.S. vital statistics for 2000 (Centers for Disease Control and Prevention 2010) were stratified by 5-year age groups, and the probability of surviving each age interval was calculated. For incidence calculations, we calculated survival during the interval so that it was the product of the probability of surviving to the interval without a diagnosis of breast cancer and the probability of not getting incident breast cancer during the interval. Surveillance, Epidemiology and End Results Program (SEER; National Cancer Institute 2010) rates for incident breast cancer for 1997–2001 were similarly stratified by 5-year age groups. Then, for each age interval, we calculated the cumulative probability of getting incident breast cancer during the interval, given survival without breast cancer up to that interval. These age-specific probabilities of getting incident breast cancer were then summed across age groups to get the background lifetime (≤ age 85 years) risk of developing breast cancer (result: 0.147). We repeated this same procedure for the exposed group, except that the age-specific background incidence rates of breast cancer were multiplied by the rate ratio predicted by the exposure–response model for the cumulative exposure for that age group. The age-specific cumulative exposures were themselves derived by assuming a constant exposure intensity (set at a given level), which then accumulated daily with a 15-year lag (i.e., starting at age 15 years). Lifetime risk for the exposed population was again the sum of the conditional probabilities of getting breast cancer in each age-group interval, and excess risk was derived by subtracting the risk in the nonexposed population from the risk in the exposed (divided by 1 – lifetime risk in nonexposed). The EC₀₁ (i.e., the estimated constant lifetime exposure level associated with a 1% excess lifetime risk of developing breast cancer) is determined by iteratively varying the exposure level until a lifetime excess risk of 0.01 is obtained. Similarly, a 95% one-sided lower confidence limit (LEC₀₁) is estimated using the upper one-sided confidence limit of the rate ratio (β + 1.64 × standard error of β) and iterating until the exposure resulting in a 1% excess risk is found.

As noted by Langholz and Richardson (2010), in linear RR models the profile likelihood may be less likely to be approximately normal than in log-linear models, so that p-values based on the Wald statistic and the change in log likelihood may differ, and Wald-type confidence bounds may diverge from bounds derived from the profile likelihood. We found such divergence in our linear RR models and used bounds from the profile likelihood to estimate the LEC₀₁ for all models except the two-parameter linear-exponential model, where we used the delta method to derive the variance of the RR and the upper one-sided 95% confidence intervals (CIs).

Table 1 gives the model fit results for the log RR and linear RR models, respectively. Overall, the linear RR models fit the data better than the log RR models, based on the corresponding model likelihoods. Based on the Akaike information criterion (AIC; Akaike 1974), which is a measure of model fit, where a lower value signifies a better fit, the linear RR model using the two-piece linear spline fits better than all other models except the square-root model, which provided the best global fit among both log-linear and linear RR models.

Results for Michaelis-Menten model are not shown. Graphically, this model looked almost identical to the square-root model, but the goodness of fit was not as good as the fit for the square-root model. The AIC was 1951.3; the best estimated k was about 13,000 ppm-days.

The EAR (additive Poisson) model did not fit as well as the RR models. Because the Poisson model requires categorized data and has a different likelihood than do the Cox models used elsewhere here, it is not possible to directly compare the goodness of fit between the Poisson EAR model and the Cox RR models. Therefore, we compared the AIC from an EAR model with the AIC from a log RR Poisson model. The EAR model fit considerably less well than the log RR Poisson model (AIC 1918 vs. 1867; smaller is better). The scale parameter was close to 1 with both models. Consequently, we do not present quantitative risk estimates based on the EAR model.

Figures 1 and 2 show the results for the log RR and linear RR models graphically. The shapes of the curves for corresponding log RR and linear RR models are similar. RR estimates based on the categorical models do not follow a monotonic trend with increasing exposure, but the overall pattern suggests increasing RRs with increasing exposure, with some attenuation at the highest exposure decile (i.e., a line fit through the other deciles would pass well above the highest decile). Note that although the cut points for forming decile categories were the same in Figures 1 and 2, the linear RR model gives different categorical results than the log-linear model. Modeling exposure as a simple continuous term appears to underestimate the exposure–response relationship suggested by the categorical point estimates from corresponding log RR and linear RR models. Modeling log-transformed exposure produces log RR and linear RR model estimates that conform reasonably well to the categorical point estimates, but the exposure–response curve is very steep in the low-dose region. Estimates from the two-piece spline models appear to conform well to the categorical results and also suggest a pronounced exposure–response slope for exposures below the knot (< 5,600 ppm-days). The square-root log RR and linear RR models also fit the data well but generate a supralinear shape (although less pronounced than the log-transform model) in the low-dose region.

Table 2 presents the resulting exposure levels that correspond to a 1% lifetime extra risk of breast cancer (EC₀₁) based on the five different linear RR models, as well as the LEC₀₁ based on using the upper one-sided 95% bound for the exposure–response coefficient. The linear exposure model, with its low slope, had the highest EC₀₁, whereas the log-transformed exposure model had the lowest, and the two-piece spline had an EC₀₁ between the two. The EC₀₁ for the two-piece spline resulted entirely from the first segment of the spline model, because the 1% excess lifetime risk resulted from cumulative exposures < 5,800 ppm-days, that is, within the range of exposures corresponding to the first segment of the estimated exposure–response curve.

Our goal here has been to illustrate key issues in risk assessment related to the choice of a statistical model when risk attenuates at high exposures. We have chosen to rely on relatively simple, easily interpretable models that may be appropriate to risk assessors operating in a public health setting, without entering into a comprehensive discussion of statistical modeling of exposure–response trends. Outside the realm of risk assessment, a wide variety of options are available for choosing the best model, but we argue that these choices are somewhat more limited if they are to conform to the needs of the risk assessor, especially if the goal is estimation of permissible environmental levels rather than occupational levels. Then the shape of the curve in the low-dose region is key. When there is a choice of different models that fit the data well, models that are approximately linear in the low-dose region are preferred, in part because they will not be highly sensitive to the relatively arbitrary choice of a point of departure for low-dose linear extrapolation. The point of departure is particularly sensitive to the shape of the curve in the low-dose region, whereas overall model fit pertains to the validity of the model across the entire exposure range.

Occupational studies for a wide variety of chemicals, often used for risk assessment for both occupational and environmental settings, have frequently shown exposure–response relationships in which the risk of cancer plateaus or attenuates at high exposures (Stayner et al. 2003). A number of reasons have been advanced to explain such attenuation, including depletion of susceptible individuals in the population, saturation of biological processes, misclassification of exposure at high levels, and the healthy worker survivor effect. One might argue that the principal reason for such attenuation is increased misclassification and mismeasurement of exposure at high levels, in which case it might make more sense to attempt to correct the exposure–response data for measurement error than to model mismeasured data. There are several arguments against this approach. First, from a practical standpoint, it is rare that a gold standard of well-measured data is available to estimate the direction and amount of mismeasurement. Second, there is no reason a priori to assume that measurement error is the principal reason for observed attenuation, and it is not necessarily clear that attenuation will result from this error. Measurement error that is nondifferential and of the Berkson type (Armstrong 1998; as in occupational studies where the mean job-specific exposure is assigned to all workers in that job) with constant error variance across exposure levels has been shown in RR models to result in a downward bias on the exposure–response relationships in RR models (Armstrong 1990). However, others have shown in log-linear RR models that increased mismeasurement of high exposures versus low exposures does not necessarily lead to estimates of exposure–response relationships biased to the null and can, in fact, lead to overestimates of exposure–response relationships (Steenland and Deddens 2000). Whether one might expect attenuation due to mismeasurement in linear RR models under this scenario is not clear.

We argue that a two-piece spline model is a good candidate for risk assessment when there is risk attenuation at high exposure levels. We have used data on breast cancer incidence and EtO exposure, which exhibit attenuation, to illustrate the two-piece spline model; however, as noted above, there are many other examples where this approach would be suitable. We also argue that additive models (linear RR) offer some advantages for risk assessors, namely, easy interpretability and a strictly linear form in the low-dose region of interest.

When there is attenuation, the risk assessor has a number of choices. For example, if the exposure–response curve plateaus as exposure increases, the risk assessor may face a dilemma of whether to consider all the data or make some relatively arbitrary decision to exclude high-exposure data altogether from the exposure–response model to focus on the low-exposure region that is the most informative and relevant for estimating low-dose risk (e.g., from environmental exposures). We have argued that the two-piece linear model provides one solution to this dilemma that avoids the need to make a subjective decision about what to define as high exposure and has the added benefit of a linear exposure–response curve in the low-exposure region.

One potential disadvantage to the two-piece spline model can arise in choosing the knot based on the likelihood, particularly when data are sparse. Sometimes the best likelihood for choosing the knot is not much better than another, that is, the profile likelihood is fairly flat or has several local maxima. This will also occur when the data are not sparse but the two-piece spline model does not provide a good fit to the data.

The variety of models illustrated here provide a statistically reasonable fit to the overall breast cancer morbidity data considered, but they have very different shapes in the low-exposure region, resulting in very different EC₀₁ estimates. The two-piece spline model yields estimates between the more extreme results obtained from the log-transformed model (very supralinear in the low-exposure region) and the untransformed (linear in exposure) model (too sublinear in the low-exposure region). The square-root model also is somewhat supralinear in the low-dose region, resulting in a low EC₀₁.

An alternative to choosing the best of a series of models that fit the data reasonably well is model averaging, in which the exposure–response results from different models are averaged using some sort of weighting scheme. When these models have the same parametric form, a weighted average of the parameters can be calculated using weights determined from the likelihood in a frequentist setting or from Bayesian posterior probabilities for the parameters (assumed to be random variables). Schwartz et al. (2008) provide an example from air pollution where a five-knot linear spline model was used with different lags on pollution exposure. A series of models was run in which the change in slope at each knot was constrained in different ways (e.g., a threshold model in which the slope of the first piece of the linear model was constrained to be 0). The parameters for slope in each piece were then averaged using as weights the Bayesian posterior probabilities, with noninformative priors or priors favoring a nonthreshold model over a threshold model. Wang et al. (2004) provide a good discussion of this same type of model averaging in a slightly different context. They show via simulation how Bayesian model averaging can outperform backward and forward selection techniques to pick the best model from a large set of candidate predictors, then show how parameter estimates for a given set of predictors from the best models can be averaged using posterior probabilities to allow inference about these parameters while reflecting the uncertainty of model selection.

When the models have different parametric forms, as in our case (e.g., square root–transformed model, two-piece model), model averaging might be done via calculation of different EC₀₁ values and taking a weighted average across them using model likelihoods or Bayesian posterior probabilities. Noble et al. (2009) have used model averaging to calculate a benchmark concentration, analogous to an EC, in a study of lung function and coal dust, using as weights Bayesian posterior probabilities with noninformative priors, across linear, square root, and quadratic models, with the selection of a set of covariates, such as height and weight, from among six possible ones. The subset of best models used in the model averaging, out of all possible models, was chosen using a type of selection rule called Occam’s window (Noble et al. 2009), based on the posterior likelihood for the model. One difficulty which arises is deciding what the necessarily arbitrary cut point will be for excluding some range of poorly fitting models from the average.

Model averaging makes more sense, in general, in situations where the results of different models are similar, for example, using a set of best models chosen via Occam’s window. Assuming a relatively broad inclusion rule, a model result that is an outlier on one side of the central tendency may have undue influence on the model average if the fit of the outlier model is reasonable (thus increasing the weight of its contribution to the average across all models) and if it is not balanced by an outlier on the other side of the central tendency. In our data, we found that different models, all of which fit the data reasonably well, produced very different estimates of EC₀₁. Although an AIC-weighted average of EC₀₁ estimates from Table 2 yielded a model-average EC₀₁ of 0.013, which was not far from the estimates based on the two-piece linear and linear-exponential models, this occurred because the two outliers (from the linear and log-transformed exposure models) were on opposite sides of 0.013—a situation that might not have occurred if we had overlooked one model form or another. Furthermore, models that are not linear in the low-dose region (which is approximately 0–1,000 ppm in our data) will lead to very different estimates of permissible exposure levels, depending on the somewhat arbitrary choice of EC₀₁ or EC₀₅ or EC₀₀₁, for example, as the point of departure for low-dose extrapolation. In our case, limiting ourselves to models that were approximately linear in the low-dose region implies limiting ourselves to the linear exposure model, the two-piece linear model, and perhaps the linear-exponential model. In our view the linear model is an outlier that did not fit the data as well as the other two models, and we would not advise averaging these three models. In sum, from both a public health and risk assessment perspective, we believe it best to choose a best model based on a combination of model fit and a desire that the resulting model be linear or close to linear in the low-dose region.

Supplemental Material is available online (doi:10.1289/ehp.1002521 via http://dx.doi.org/).

This study was supported by the U.S. Environmental Protection Agency.