Use of Multiple Imputation to Correct for Bias in Lung Cancer Incidence Trends by Histologic Subtype

Lung cancer is the leading cause of cancer death in women and men in the United States. On average, only approximately 15% of newly diagnosed cases survive for 5 years or longer (1). Histologically, lung cancers are classified as small cell and non–small cell (NSC) carcinoma (2). The latter is usually further divided into squamous cell carcinoma, adenocarcinoma, and large cell carcinoma. Within the NSC category, etiologic and morphologic differences by histology have been recognized, but in the past, treatment and prognosis were considered relatively homogeneous for different histologies of the same stage. Emerging data now increasingly identify subsets of adenocarcinoma (3) and squamous histologies (4) with specific genetic alterations. For example, the epidermal growth factor receptor (EGFR) protein overexpression and activating EGFR mutations, associated with responsiveness to EGFR therapies (tyrosine kinase inhibitors; refs. 5 and 6), are almost exclusively found in adenocarcinoma histology. Similarly, echinoderm microtubule-associated protein-like 4 (EML4)-anaplastic lymphoma kinase (ALK) rearrangements are also more common in adenocarcinoma and these mutations indicate responsiveness to another therapeutic agent, crizotinib (7). As we move into the future, clinical strategy for tumor management will be determined by molecular studies of the tumors and their underlying mutations (8). Inherited variation in lung cancer that has been identified may eventually have therapeutic implications in terms of efficacy and side effects. The recent results from the National Lung Screening Trial further suggested that histology might be attributable to the differential computed tomography (CT) screening efficiency (9). As the broader implications of histologic classification are becoming increasingly relevant to screening, treatment, prognosis, and etiology, so will the examination of temporal trends separately for each subtype.

Cancer registry data collected by the National Cancer Institute (NCI)'s Surveillance Epidemiology and End Results (SEER) Program have been a primary source of data for providing national trends of lung cancer incidence and mortality (10). SEER registries have been coding cancer histology according to the International Classification of Diseases for Oncology (ICD-O). In the 1990s, pathologists tended not to report NSC carcinomas with specificity because their treatments and prognoses were considered similar, thus an increasing number of cases are coded with 8010 (carcinoma, NOS) since 1980. In recognition of this trend, 8046 (NSC carcinoma) was added into ICD-O-3 in 2001 to group cases that could not be classified beyond the exclusion of small cell. Collectively, the percentage of cases coded with 8010 or 8046 increased dramatically, from 5% in 1982 to more than 22% in 2005 (11). Some of these cases could have been derived from one of the specific histologic subtypes, which would have subsequently reduced their incidence rates. However, this increasing use of nonspecific codes did not continue. In light of the advances in cancer research and therapy, increasingly NSC cases have been diagnosed with more histologic specificity (12) over the last few years, which may have driven up the rates for squamous or adenocarcinoma. Such differential use of nonspecific morphology codes could bias the estimated temporal trends of histologic subtypes and complicate interpretations. Appropriate statistical adjustments are necessary to improve the quality of inferences using the authoritative cancer registry data, which otherwise has been compromised by the unavoidable limitations imposed by the imperfect earlier classification system.

Multiple imputation (MI) has been shown to be a useful approach for handling measurement or coding changes for settings both in the presence (13–16) and absence (17) of calibration data (observations that are measured in all measurement scales or coding systems). When calibration data (usually on a random subsample) are available, one can generate plausible values in all measurement scales from an imputation model and analyze the imputed data using the preferred scale. For the issue associated with the change in the use of nonspecific morphology codes, that is 8010 and 8046, 2 types of calibration data could be useful for correcting coding inconsistency. The first type comprises cancer cases that are originally assigned to a nonspecific code, but are updated with a specific code through reexamination. Such data provide information about the association between nonspecific and specific histologies that one can use to recover the missing histology for all cases with a nonspecific code. Because nonspecific codes no longer exist in the imputed data, the trend analysis of incidence by histology is valid (provided that the imputation model is correct). The second type consists of cancer cases coded in multiple classification systems. Using these data as a bridge, one can convert data from one system to another. Although nonspecific codes still exist, temporal comparisons of imputed histology in any classification system is valid because coding consistency is maintained. However, neither type of calibration data could be easily obtained because of practical reasons, such as budget constraints and the lack of diagnostic data sources. Thus, this problem becomes a missing data issue where the specific histologies for cases with a nonspecific morphology code are missing and an assumption about the association between the missing specific histology and observed data (18–20) is required. We make a reasonable assumption that for cancer cases with similar tumor, treatment, survival, patients' demographic characteristics, the distribution of nonspecific and specific histology is similar. Based on this assumption, we developed an MI approach using the sequential regression imputation method (SRMI; ref. 21) to redistribute cases without specific histology to one of specific subtypes, thus correcting the biased estimates of incidence rates.

We selected 522,416 malignant lung cancer cases diagnosed from 1975 to 2010 from the SEER 9 registries database (including Atlanta, Connecticut, Detroit, Hawaii, Iowa, New Mexico, San Francisco–Oakland, Seattle–Puget Sound, and Utah). We created 6 histologic categories according to the most recent NCI's SEER Cancer Statistics Review (1) and they are small cell carcinoma (8041–8045), squamous and transitional cell carcinoma (8051–8052, 8070–8084, 8120–8131), adenocarcinoma (8050, 8140–8149, 8160–8162, 8190–8221, 8250–8263, 8270–8280, 8290–8337, 8350–8390, 9400–8560, 8570–8576, 8940–8941), large cell carcinoma (8011–8015), other NSC carcinoma (8020–8022, 8030–8040, 8090–8110, 8150–8156, 8170–8175, 8180, 8230–8231, 8240–8249, 8340–8347, 8561–8562, 8580–8671), and other specified and unspecified types (8680–8713, 8800–8912, 8990–8991, 9040–9044, 9120–9136, 9150–9252, 9370–9373, 9540–9582, 8720–8790, 8930–8936, 8950–8983, 9000–9030, 9060–9110, 9260–9365, 9380–9539, 8000–8005). We singled out 8010 and 8046 from these categories, for which we performed statistical adjustments. We excluded cases with other specified and unspecified types because their incidence is not likely to be affected by the recent change in coding system. We also excluded cases that were not histologically confirmed or with unknown histologic confirmation status, because their diagnoses tended to be inaccurate and lacked specificity. The final sample size for this analysis is 470,326.

We treated the cases with 8010 or 8046 as missing data that we dealt with by MI (22). This MI approach took each case with missing histology and imputed it with a specific histologic subtype. Cases coded with 8010 were imputed with one of the 5 carcinoma subtypes, that is small cell, squamous, adenocarcinoma, large cell, and other NSC. For cases coded with 8046, the imputation was limited to one of the NSC subtypes, that is excluding small cell. This process was repeated independently 10 times to create 10 completed datasets to account for imputation uncertainty. Age-adjusted incidence rates (using the 2000 U.S. standard population in 19 age groups) were estimated from each completed dataset in the same way as using the original dataset, thus producing 10 sets of estimates. We then combined these estimates to produce MI estimates. For a single incidence rate, the MI point estimate was the average of 10 imputed data estimates. The associated standard error was calculated by combining the average of the squared standard errors of the 10 estimates and the variance of the 10 rate estimates (22). Joinpoint linear regression models (23) were used to fit connected linear trends on a log scale with up to 4 joinpoints using the Joinpoint regression program version 3.5.0 developed by the NCI. Annual percentage change (APC) with a corresponding 95% confidence interval (CI) was calculated to describe each joined trend.

The nonspecific histologic diagnoses are highly likely to have nonrandom characteristics. For example, patients may not merit further histologic diagnostic procedures because they have diseases too advanced to permit curative surgery (i.e., stage IIIB or greater) or because their medical status preclude surgery or other modalities with curative intent. When surgery is not a clinical option, obtaining adequate tissue to establish a histologic subtype may be impossible and, in this circumstance, clinicians may elect to forgo further histologic classification. Therefore, we considered using the information that is predictive of histology and the missingness of specific histology to recover the incomplete specific histology. We assumed the missingness is random conditional on this information, and this assumption has been shown to be reasonable in most practical situations (24, 25).

Specifically, we selected the covariates to be included in the imputation model following the principle of reducing missing data bias in a statistical analysis (26). Sociodemographic covariates include age, gender, race, Hispanic origin, nativity, and marital status. Covariates describing tumor characteristics and treatment include tumor size (27), grade, stage, survival time, and receipt of cancer-directed surgery. Certain therapies have shown to be more responsive in some histologic subtypes, thus making them important predictors. However, such information can only be made available for patients 65 years and older through the linked SEER-Medicare database (28) for 1991 and later. Considering the lack of analytics tools to handle the dynamics of the availability and access to particular regimen over time and patients' age, we did not include more detailed treatment variables in the model. We also did not include lymph node involvement in the final model because it is highly collinear with stage. We included a nominal variable of 9 SEER registries to reflect the variability among registries in the use of nonspecific morphology codes. Cancer diagnosis year was entered into the models as a nominal variable (instead of a continuous variable) to relax the temporal assumption about the intervariable relationships. Smoking and socioeconomic deprivation are also strongly predictive of histology (29), but they are not routinely collected in SEER. To substitute, we used county-level smoking prevalence estimates obtained from the Model-based Small Area Estimates Projects of NCI (http://sae.cancer.gov/; ref. 30), and poverty prevalence estimates from the 2000 U.S. Census Bureau (31).

Because missing histology cannot be imputed for cases that are associated with missing covariates using simple regression-based imputation approaches, we developed an algorithm using SRMI technique to deal with multivariate missing data with arbitrary missing patterns. Specifically, SRMI fits a conditional model for each variable at a time on the remaining variables sequentially for multiple rounds to achieve convergence. The form of conditional model depends on the type of variable imputed. Our algorithm offers 2 new capacities beyond what is available in existing SRMI-based imputation packages, such as IVEware (http://www.isr.umich.edu/src/smp/ive/) and MICE (http://cran.r-project.org/web/packages/mice/index.html). First, for imputing binary data (categorical variables with more than 2 levels can be expressed as a series of nested dummy variables), we used ridge-penalized logistic regressions (32, 33) to improve imputation precision in the presence of binary outcome with skewed distribution and highly correlated covariates (34). The standard approach for imputing missing binary data is usually based on a logistic regression model (21, 35). However, the adequacy of logistic models could highly depend upon the extent to which the binary outcome is balanced and there is an absence of collinearity. In the presence of either condition or both at the same time, logistic regression coefficients may still be unbiased, but the precision could be very low, which could lead to poorly imputed data. The proposed approach improves the imputation by estimating a penalized log likelihood to obtain coefficients estimates with minimum prediction errors. Optimizing the penalty parameters is critical and usually requires intensive cross-validation studies (36). We follow the simplified approach proposed by Yu (34) and obtain the optimized parameters directly from the data by estimating the unrestricted log likelihood. The remaining steps are similar to those when standard logistic models are used (21). Second, we added a module to impute discrete right-censored survival data. For the data we chose for this study, more than 25% of survival time was censored because the patient was still alive at the end of study or died from other causes. Because both survival and censoring are highly correlated with histology as well as other covariates such as age, stage, tumor size, and grade, it is problematic to use relatively simple approaches, such as the indicator method where censoring is taken care of by including a censoring indicator (37, 38). The proposed method applies the MI principle to impute the censored time with a plausible future survival time. Specifically, to generate the imputed values, we first aggregate continuous survival time (in month) into several meaningful categories and sort them in an increasing order of survival. We then define an imputing risk set for each censored case as the cases with observed survivals no shorter than the censoring time. Using data from this imputing risk set, we finally estimate the predictive conditional distributions of survival categories, from which we randomly draw a value to be the imputed survival. Note that the possible value of an imputed survival is always equal to or longer than the censoring time category. This is reasonable because a censored case could only die at a later time in its own survival category or be still alive and die at a future category, but not die at a past category. This imputation process starts with censored cases in the first survival category and cycles through all categories to complete one imputed survival data. Because the survival is now a discrete variable, we estimate its predictive conditional distribution using nested ridge-penalized logistic models similar to what we have outlined for categorical data. Furthermore, to deal with the inconsistency in stage definitions over time, we conducted the imputation separately for 1975 to 1982, 1983 to 1987, and 1988 to 2010, so that staging is comparable within each period.

To explore information recovery from the MI in estimating the distribution of histology, we generated a simulated dataset from the analysis data with only complete observations included (n = 10,659). We considered a situation similar to the main analysis where histology is missing at random and the probability of the induced missingness is determined by a logistic regression model with the coefficients estimated using the analysis data. The rate of induced missing data was 8.4% (the observed missing rate was 10.0% for the portion of data with all covariates observed). Twenty imputed datasets were generated using the proposed approach and the standard logistic regression method, respectively.

The ridge-penalized logistic regression model outperformed the standard logistic regression model in recovering the missing information based on the Akaike information criterion (AIC; ridge-penalized method: AIC = 31,008 and standard method: AIC = 31,065). The imputed distributions of histology obtained using the proposed method were similar to the complete data distribution (with absolute difference less than 2% in estimating the percentage of cases in each histology and gender group). We also calculated the overlap probability (39) to evaluate how much the associated 95% CI estimated from the imputed and complete data overlap. Suppose (L_imp, U_imp) and (L_com, U_com) are the 95% CIs for estimating P, the percentage of adenocarcinoma among men, using the imputed and complete data, respectively. The probability overlap in the CIs for P is |$I} = \frac{1}{2}\left[ {\int_{L_{imp}}^{U_{imp}}\ {f_{com} (t)dt} + \int_{L_{com}}^{U_{com}}\ {f_{imp} (t)dt}}\right]$|⁠, where f_imp and f_com are the distributions of P computed under the imputed and complete data, respectively. Note that f_com could take a different form of distribution depending on the type of statistics for which one wish to obtain estimates, but f_imp is always t-distributed according to Rubin's rules (22). I takes value 0.95 if 2 CIs overlap perfectly and 0 if they do not overlap at all. A large value in I suggests that the imputed data highly maintains the analytical properties of the complete data. This measure provides more information than a simple comparison of 2 point estimates by also considering the standard errors. Estimates with large standard errors might still have a high CI overlap even if their point estimates differ considerably from each other because the CI will increase with the standard error of the estimate. In this simulation study, most overlap probabilities (for estimating the distributions of cases by histology and gender) were more than 0.8, which suggested a very strong agreement, with a few exceptions in which the probabilities were around 0.75, which still suggested a strong agreement. These evaluation results provided strong evidences for model adequacy in the proposed method.

Table 1 shows the distribution of histologic categories by histology confirmation status. Ninety percent of cases are histologically confirmed. Among the cases that are not confirmed and the cases for which the confirmation status is unknown, 8010 accounts for about 50% of the total whereas 8046 only accounts for less than 2%. Possible explanation for the differential use of 8010 and 8046 could be that the latter is mainly used when histologic diagnosis, although not quite specific, exists, and the former is also used when the diagnosis is not available.

Table 2 shows the distributions of lung cancer cases by histology and selected covariates. All covariates are closely associated with histology. Men and older patients were more likely to be diagnosed with squamous type. Squamous and adenocarcinoma tumors tended to be more well-differentiated than large cell and other specific NSC tumors. Squamous and large cell tumors tended to be larger at diagnosis. Small cell tumors were likely detected at a later stage (61.6%) as compared with other types. In contrast, tumors of squamous and adenocarcinoma types tended to be detected at early stage. There are also a few notable differences in the use of nonspecific codes across registries. For example, a lower use of 8046 (15.2% in 8046 compared with the overall percentage of 20.8%) is observed in Detroit, and a higher use of both 8010 (16.9% compared with the overall percentage of 15.0%) and 8046 (19.9%) is observed in Seattle. The use of nonspecific code is also slightly higher for cases not reported by a hospital (2.8% in 8010 and 2.9% in 8046 compared with the overall percentage of 1.8%). These variables are also predictive to the use of nonspecific morphology codes. As we expected, tumors without specific histologic diagnosis tended to be less well differentiated, diagnosed at a late stage, had shorter survivals, and were less likely to be candidates for surgery.

Figure 1 shows the percentages of cases coded with 8046 and 8010 by year of diagnosis for men and women separately. The temporal distributions are similar for both genders. The percentage of cases coded with 8010 had increased from 1982 until the introduction of 8046 into ICD-O-3 in 2001, when it dropped to around 3%. There seems to be a smooth compensation between 8010 and 8046 in 2001, which suggests that 8010 and 8046 are probably used interexchangebly in practice.

Figure 2 shows the rates of incidence by imputed histology among cases coded with 8010 or 8046. Overall, the amount of imputed histology differs by histologic subtype and year. For both 8010 and 8046, the rates of incidence raised by imputation were greatest for adenocarcinoma and squamous. For both histologic subtypes, the rates followed an n-shaped pattern over the most recent 15 years. Small cell was the third most raised category, although only contributed from imputing 8010 cases, and the amount of increases was relatively stable over time.

Figure 3 compares the before and after imputation temporal trends in age-adjusted incidence rate of lung cancer by histology for men and women separately (see Table 3, for detailed results of the joinpoint trends analysis) The numbers listed over (imputed) or under (original) each segment represents the APC for that portion of the trend and an asterisk indicates a statistically significant trend at 0.05 level. The rates for 8010 (small cell type) and 8010 and 8046 combined (NSC subtypes) are also included in these plots to help examine how cases are distributed by the imputation procedure.

The imputation adjustment affected the incidence trends differently for each histologic subtype. For small cell in both genders, the original and imputed trends are similar. For squamous cell cancer in both genders and adenocarcinoma in men, the trends showed a similar pattern overall from 1970 to early 1990s before and after imputation. From early 1990s to 2005, the decreasing trends also remained unchanged after imputation, but the pace of decline slowed. After 2005, the increasing trends based on the original data had been replaced by the steady continuations of earlier decreasing trends for squamous and adenocarcinoma in men, a constant trend for squamous in women, after imputation. For adenocarcinoma in women, the trends, before and after imputation, exhibited similar patterns overall before early 1990s. From 1992 to 2007, the plateau followed by an increasing trend started in 2004 changed to a continuously increasing trend after imputation. It is also worth noting that the imputed rates showed a nonsignificant decreasing tendency during the most recent 3 years starting in 2007. For large cell cancer and cancer in other specified NSC type, the imputed rates were similar to the original rates and the imputation did not change the overall trends.

To rule out the possibility that changes in trends may be because of the absence of cases that are not histologically confirmed or have missing confirmation status, we conducted a sensitivity analysis on all cases. The imputation affected the trends similarly (see Supplementary Fig. S1 and Table S1, for detailed results on the rates and jointpoint analysis), which suggests that excluding these cases does not affect the overall findings and conclusions.

In cancer surveillance data collections, it is common for the morphological classification systems to change to reflect the contemporary pathology practice. Hence, the data often comprise cancer cases coded one way at one time and others a different way at another time. When classification systems differ in coding histology, temporal inferences by histologic subtype can be misleading and difficult to interpret. Without access to calibration data to inform the underlying distribution of histology among cases coded without specificity or the association in histology between editions of classification systems, we carefully developed an MI approach to correct for biases in statistical inferences about temporal trends of lung cancer incidence based on the MAR assumption.

Although this assumption is not empirically testable, we argue that MAR is reasonable in our setting because we have identified and included into the imputation models an extensive set of auxiliary variables that can explain the missingness of specific histology, for example receipt of cancer-directed surgery, and that are correlates of histology, for example the stage, grade, and size of a tumor, as well as patient survival. Other important variables that could enhance the MAR assumption plausibility are patients' smoking status and socioeconomic status (40, 41), for which we substituted county level estimates at 2000 (pooled estimates from 2000 to 2003 for smoking) from the decennial census because they are not routinely collected in SEER. Although such estimates are not available for every diagnosis year, we believe the ranking of a county in smoking prevalence or poverty level relative to the rest of the country remains relatively unchanged over time. The potential confounding between smoking status and poverty (40) is not likely a cause for concern in our analysis because both are aggregate measures and neither is a strong predictor to histology after conditional on other patient-level information.

Ensuring the plausibility of MAR assumption imposed 2 modeling challenges of handling a large number of variables with missing data and a general missing data pattern, which often cannot be adequately addressed by simple imputation methods (21). The proposed MI approach based on SRMI is particularly suitable to this complex situation because of its flexibility in specifying and fitting conditional distributions. The search for refined ridge-penalized logistic regression imputation models is necessary because the standard SRMI approach (based on logistic regressions) might be inadequate in handling a categorical outcome with a skewed distribution (e.g., certain histology categories only contain 3%–6% of cases) and correlated covariates (e.g., stage and survival). The simulation study demonstrated the adequacy and prediction benefits of the proposed semiparametric models.

The amount of lung cancer cases lacking specific histologic subtypes was predominantly associated with the year of diagnosis, which reflected the evolution of SEER coding algorithms and recent changes in diagnostic practice. The imputation raised the incidence rates across the entire study period for both genders and histology subgroups. However, the magnitudes of the elevations varied. Of the various histologic subtypes, the most impacted were squamous and adenocarcinoma, on which the most pronounced impacts occurred during the last decade. This result further supports our hypothesis that 8010 and 8046 are mainly used to group cases, which could have been coded as either adenocarcinoma or squamous type if more coding information were extracted and available to support detailed histologic coding. For both subtypes, the decreasing trends from early or mid-1990s to 2005, had persisted, although at a slower pace. The increasing trends after 2005 are apparently an artifact of this coding change and imprecision in histopathologic classification, which, after imputation, became a continuation of earlier decreasing trends. The sensitivity analyses including cases that are not histologically confirmed or have missing histologic confirmation information showed similar results.

We classified lung cancers according to a schema developed based on Travis and colleagues (42) and earlier versions of ICD-Os. WHO recently published a revised version of the histologic grouping for lung cancers (43). Different histologic classification systems have been used in practice, for example, the recently published classification schema by the International Agency for Research on Cancer of the WHO (43) in 2007. The differences between this new classification and the one used in this research are summarized in Supplementary Table S2. Because the groupings of the most frequently used morphologic codes are consistent between the 2 schemas, we suspect that the effect of using this alternative schema on the inferences of incidence trends is noticeable for the histologic subtypes that we investigated in this research.

In summary, molecular, genetic, and etiologic features are increasingly associated with histology distinctions (3, 4, 44). Progress in linking molecular features to morphology will facilitate mechanistic understanding and further characterization of the molecular and genetic features specific to histologic subtypes in lung cancer. These considerations, along with the emergence of targeted therapies within specific histologic subtypes especially adenocarcinoma, clearly indicates that accurate population tracking of trends by lung cancer histology will be increasingly important in the future, and that the MI technique applied in this study can help refine these trends. Planned data collections for bridge data in the future will further enhance the quality of data augmented by MI.

No potential conflicts of interest were disclosed.

Conception and design: M. Yu, E.J. Feuer, K.A. Cronin, N.E. Caporaso

Development of methodology: M. Yu, E.J. Feuer, K.A. Cronin

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): M. Yu

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): M. Yu, E.J. Feuer, K.A. Cronin, N.E. Caporaso

Writing, review, and/or revision of the manuscript: M. Yu, E.J. Feuer, K.A. Cronin, N.E. Caporaso

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): M. Yu

Study supervision: M. Yu, K.A. Cronin

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.