Lung cancer incidence decreases with elevation: evidence for oxygen as an inhaled carcinogen

Data collection & preparation

From 11 publicly-available databases, we compiled US county data on cancer risk factors, environmental features, demographics, and quality control metrics, while avoiding redundancy. When selecting resources, we balanced several considerations including coverage, precision, collection period, and accessibility. To capture the long latency period of cancers and minimize observational error, we preferred collection periods preceding the cancer incidence timeframe and spanning multiple years. Resources were integrated using county FIPS (Federal Information Processing Standards) codes if available and name matching otherwise. Variables collected as cancer incidence predictors are displayed in Table 1. Variables were averaged over their entire collection periods. Unless otherwise noted, averaging was performed by the source databases.

County filtering

We restricted the analysis to states in the contiguous United States with elevation spans exceeding 3,000 m. The selected states—AZ, CA, CO, ID, MT, NV, NM, OR, UT, WA, WY—contained 414 counties composing the Western United States. Next, counties were filtered for quality control. Counties with populations below 10,000 were excluded due to high missingness (values were missing for many of the variables) and observational error (values were present but subject to large margins of error, evidenced by source-reported confidence intervals). Counties with high Native American composition or immigration rates were considered potentially problematic: cancer rates among Native Americans are prone to misestimation (Puukka, Stehr-Green & Becker, 2005); and immigrants accumulate cancer risk prior to migration, outside of their destination county. Accordingly, we found that predictions of all-site cancer incidence, based on eight general demographic and health predictors, diverged from reported incidence for Native American and immigration-rich counties (Fig. S1). Selecting exclusion thresholds corresponding to this divergence, we omitted counties with five-year immigration rates exceeding 40% or Native American population exceeding 25%. After filtering, 260 counties remained.

Regression analysis

We evaluated the association between elevation and cancer incidence using multivariate linear regression. Counties were weighted by their population square root up to a maximum population of 250,000 where measurement uncertainty leveled off to minimal levels. The weighting scheme accounted for increasing measurement uncertainty among low population counties without granting heavily populated counties an overwhelming influence.

To minimize confounding effects, we selected well-established factors to include as covariates with elevation. We avoided excessive collinearity (Dormann et al., 2013) by carefully identifying major cancer-specific risk factors with available county-level estimates. Smoking, radon, fine particulate matter, and percent male were identified for lung cancer; female smoking, mammogram, and drinking for breast cancer; smoking, drinking, diabetes, meat consumption, and percent male for colorectal cancer; and meat consumption for prostate cancer. We included six additional covariates—metro, white, black, education, income, and obesity—for all cancers to indirectly account for unknown or immeasurable risk factors or biases. Since a large degree of risk is shared between cancers (Danaei et al., 2005), for each cancer we included the incidence of all other cancer. In addition to elevation, a total of 11 covariates were included for lung, 10 for breast, 12 for colorectal, and 8 for prostate (Table 1). We created cancer-specific datasets by removing counties with any missing data for included variables. Standardized versions of each dataset were created by converting cancer incidence and all predictors to weighted z-scores.

We employed two regression methods, best subset and lasso, in parallel. The best subset approach allowed us to force elevation into the model and exhaustively evaluate all possible models while remaining amenable to statistical interpretation. However, evaluating all possible subsets creates the potential for overfitting. Lasso addresses this concern by introducing coefficient shrinkage and variable selection (Tibshirani, 1996). Despite efforts to exclude redundant predictors, variables were characterized by moderate yet pervasive collinearity (Fig. 1). In the best subset method, severe collinearity could cause unstable and unreliable coefficient estimates with inflated standard errors. In addition to scrutinizing best subset results, we adopted the lasso, whose variable selection mechanism tends to include a single member from a group of correlated predictors. This characteristic makes the lasso effective at identifying truly associated predictors in the presence of high collinearity (Dormann et al., 2013).

County stratifications

To investigate the potential of a smoking-elevation interaction affecting lung cancer, we partitioned counties into smoking prevalence terciles as follows: high-smoking (49.9–61.9), mid-smoking (44.9–49.9), low-smoking (28.2–44.9). Within each tercile, we regressed lung cancer incidence against elevation for visual examination. For statistical evaluation, we refitted the optimal best subset model with an added interaction term (standardized smoking prevalence × standardized elevation).

To mitigate uncontrolled confounding (i.e., omitted-variable bias), we stratified the lung cancer dataset by state. Health policy and data collection are often enacted at the state level, making stratification by state a sensible choice for maximizing within-strata homogeneity. Within each stratum, lung cancer incidence was regressed against elevation and smoking prevalence. The elevation effect size was estimated across the eleven state-specific models using a fixed-effects meta-analysis, which averaged elevation coefficients weighted by their inverse variances.

Population subgroupings

We evaluated the association between elevation and lung cancer incidence measured for the following population subgroups: under 65 years old, 65 or older, males, and females. These subgroups were chosen because exposure to risk factors, such as occupational or lifestyle hazards, often segregates by sex or age. For this analysis, we created a dataset with counties that had no missing data for the four subgroup outcomes and the predictors from the optimal best subset model. Proceeding with the best subset covariates, we fit a separate regression model for each population subgroup. We used sex-specific smoking prevalence for the male and female models.

Elevation substitutions

Environmental variables without established cancer risk were not evaluated in the previous regression analyses to avoid problems of collinearity. For lung and breast cancer, we investigated whether substituting elevation with each of seven environmental variables produced a more likely model, which could indicate an indirect elevation-cancer association. For each environmental variable, we performed best subset regression to find the BIC-minimizing set of predictors. The covariates evaluated for each cancer matched those from before (Table 1) with the exception of radon and fine particulate, which were excluded as lung cancer covariates and instead included as elevation replacements. For both cancers, the increase in minimum-BIC for each substitution compared to elevation was recorded. Change in BIC was converted to a Bayes factor, $K\approx e^{-\frac{1}{2}\left(\Delta BIC\right)}$ (Raftery, 1995, p. 139). K > 1 provides evidence favoring replacement whereas K < 1 provides evidence against.

Software

Analyses were performed using the statistical-computing language R. County elevation computation relied on the GIS packages raster and rgdal. The best-subset regression analysis used the leaps package, which efficiently identifies top performing models from the complete search space. The glmnet package implemented the lasso (Friedman, Hastie & Tibshirani, 2010). The state-specific lung cancer elevation coefficients were meta-analyzed using the metafor package (Viechtbauer, 2010). Tables were exported using the Hmisc package. Plots were created with the ggplot2 package. Correlation plots were ordered using Ward’s hierarchical clustering.

Data availability

The county-level dataset compiled for this study is available (Data S1). The project GitHub repository (https://github.com/dhimmel/elevcan) contains the code used to perform analyses as well as all intermediate files.

Strong, negative association between elevation & lung cancer incidence

Performing best subset regression for each cancer, we found a highly significant, strong negative association between elevation and lung cancer incidence with a standardized coefficient (β_z) of −0.35 99% CI [−0.46, −0.25] (p < 10⁻¹⁶, one-tailed t-test) (Table 2). Lung cancer incidence decreased by 7.23 [5.67–8.80] cases (per 100,000 individuals) per kilometer rise in elevation, equating to 12.7% [9.1%–16.4%] of the mean lung cancer incidence. For other cancers, we found a weak, negative association with breast cancer (β_z = −0.15, p < 10⁻²) but not with colorectal (p = 0.88) or prostate (p = 0.97) cancer.

The optimal (BIC-minimizing) models contained five predictors for lung and colorectal cancers, six predictors for breast, and four predictors for prostate cancer (Table 2). Within each cancer, we compared the elevation coefficients across a range of model sizes (Fig. 2). Unique to lung cancer, elevation confidence intervals were consistent and wholly negative, indicating robustness to collinearity as well as to confounding by included covariates. Other cancers displayed greater coefficient variability and uncertainty, possibly due to covariate collinearity with elevation, which led us to implement lasso regression.

Lasso regression performs variable selection that operates well under moderate collinearity and coefficient shrinkage that prevents overfitting. Using a conservative setup of the lasso, we again observed a strong, negative association between elevation and lung cancer incidence with a standardized coefficient of −0.33, changing minimally from the best subset estimate (β_z = −0.35) despite the strong regularization of the lasso (Table 3). For breast cancer, where the best subset model yielded an elevation coefficient with high uncertainty, the lasso reduced the estimate to a trivial level (β_z = −0.02) indicating that overfitting contributed to best subset negativity. Meanwhile, the elevation term was absent in the colorectal and prostate lasso models. Together the regression methods indicated a negative association with elevation that was unique to lung cancer in terms of strength, significance, and statistical robustness.

Models accurately assess known cancer associations

Models produced for each cancer by best subset (Fig. 3A) and lasso (Fig. 3B) regression corresponded with the literature. The lasso (and best subset) models explained 67% (70%) of variation in lung cancer incidence, 51% (57%) in breast, 29% (34%) in colorectal, and 9% (19%) in prostate, (Tables 3 and 2) mirroring a previously described trend in fraction of risk attributable to modifiable factors for each of the four cancers (Danaei et al., 2005).

For lung cancer, both regression methods found previously characterized positive associations with smoking prevalence, percent of black residents (Greenlee et al., 2000, p. 10), and rate of other cancer (Ahlbom et al., 1997), as well as large negative associations with elevation and education. The lasso also found a small negative association with radon (β_z = −0.01), attributable to radon’s strong positive correlation with elevation (Fig. 1). Including covariates sharpened the association between elevation and lung cancer, as evidenced by elevation’s higher partial R² = 0.252 in the multivariate model versus the bivariate R² = 0.202 (Fig. 4). Of note, both best subset and lasso regression attributed the two largest effect sizes to smoking (best subset: 0.59 and lasso: 0.52) and elevation (−0.35 and −0.33). The best subset model also found that smoking (p < 10⁻³⁵) and elevation (p < 10⁻¹⁶) were the two most significant associations with lung cancer incidence (Table S1).

For breast cancer, the best subset regression model captured known positive associations with education and income (Devesa & Diamond, 1980), other cancer (Ahlbom et al., 1997), metropolitan status (Hall et al., 2005), and percent white (Greenlee et al., 2000, p. 10) (Howlader et al., 2014). The lasso corroborated the top three positive associations found by best subset regression, and added other known factors, mammography frequency and income (Devesa & Diamond, 1980), in place of percent white. As for lung and breast cancer, colorectal models were sensible, finding positive associations with other cancer, meat consumption (Norat et al., 2002), and percent of black residents (Greenlee et al., 2000, p. 10), as well as a negative association with education (Jemal et al., 2008)—all previously reported. In line with past ecological analyses (Danaei et al., 2005), prostate models were incomplete and inconclusive, failing to find the known positive association with percent black (Greenlee et al., 2000, p. 10). Overall, the plausible and well-fitting lung, breast, and colorectal models reflected our ability to recapitulate known associations and therefore characterize the prospective association of elevation with lung cancer.

Elevation’s association with lung cancer is robust to stratification & subgrouping

Given the respiratory intersection of oxygen inhalation and smoking, we asked whether elevation associated differently across smoking prevalences. Stratifying counties into smoking terciles, the strong effect of smoking on lung cancer was illustrated by the non-overlapping confidence bands (Fig. 5A). Tercile slopes were approximately parallel providing no evidence for an interaction between smoking and elevation. Corroborating the lack of interaction, an added smoking × elevation term was not significant (p = 0.47) when refitting the best subset model.

While focusing on the elevation-varying Western United States minimized the risk of regional confounding, possible differences in health policy or practice between higher and lower elevation states could still exist. State-specific models that accounted for smoking found negative elevation coefficients for ten of eleven states (Fig. 5B). Moreover, while small intrastate sample sizes created coefficient uncertainty, meta-analysis estimated that lung cancer incidence decreases by 10.8 99% CI [7.0–14.6] cases (per 100,000 individuals) per kilometer, matching the best subset regression confidence estimate of 7.2 99% CI [5.7–8.8] and indicating that state-based confounding did not underlie the negativity of the region-wide elevation association. To evaluate the consistency of state-specific elevation effects, we calculated I², which equaled 71% 95% CI [35–90%]. Despite its uncertainty, the observed I² suggests moderate to high heterogeneity across states (Higgins et al., 2003), which is unsurprising given the small sample sizes and limited covariate adjustment of the state-specific models.

Disparate habits and lifestyles exist across age groups and sexes. Using subgroup-specific lung cancer incidences, we refit the best subset model asking whether elevation association would change (Fig. 6). All four subgroupings showed significant, negative association with elevation: under 65 years old (p < 10⁻¹²), 65 and older (p < 10⁻¹⁷), males (p < 10⁻¹⁴), and females (p < 10⁻¹⁸). Standardized coefficients were large, ranging from −4.39 to −3.58. Subgroup confidence intervals all overlapped the confidence interval from the global model, showing no evidence for subgroup-specific effect mediation or confounding.

Lung cancer associates with elevation over environmental correlates

Rising elevation leads to lower atmospheric pressure, which helps drive a repertoire of interconnected climatic changes, including perturbations in sun exposure, temperature, and precipitation. This trend was apparent in our data, as many environmental variables correlated with elevation, and thus also with lung and breast cancer (Fig. 7A). If either cancer’s elevation association was indirect of atmospheric pressure but rather a product of secondary climatic changes, we expected that environmental correlates could outperform elevation in best subset regression. For example, vitamin D synthesis is stimulated by sunlight and UV exposure (Gilchrest, 2008). The hormonally active form of vitamin D, calcitriol, potentially possesses anti-cancer properties (Krishnan & Feldman, 2011). Sunlight and UVB exposure correlate positively with elevation in our data. Hence, a reasonable proposition would be that increased vitamin D synthesis is driving elevation’s association with lower cancer rates (as posited by Hayes, 2010) and that replacing elevation with UVB would improve model likelihood.

We approximated the likelihood (as a Bayes factor, K) that any correlated environmental variable could replace elevation in our models of lung and breast cancer (Fig. 7B). Remarkably, for lung cancer, all of the variables tested produced models between 10⁸ and 10¹³ times less likely than the elevation-including model. Among these, fine particulate matter was over one trillion times (K < 10⁻¹²) less likely than elevation. In contrast, breast cancer was more effectively explained by environmental correlates or secondary climatic changes, such as precipitation, UVB, and sunlight, than by elevation. For example, precipitation was 17 times more likely than elevation. Together, these results indicated that lung cancer but not breast cancer was directly associated with atmospheric pressure.

Radon and UVB associations with lung cancer confounded by elevation

As a consequence of elevation’s numerous environmental correlates, we speculated that previous ecological studies of lung cancer may have fallen prey to uncontrolled elevation confounding. Radon (ρ = 0.71) and UVB (ρ = 0.50) correlated highly with elevation in our data (Fig. 7A) and both had previously been reported to associate negatively with lung cancer (Cohen, 1995; Hayes, 2010). To test for confounding, we fit three models of lung cancer incidence with either radon or UVB exposure plus the following predictors: (1) smoking alone; (2) smoking and elevation; and (3) the predictors from the optimal best subset model. Similar to previous studies, model 1, which excluded elevation, identified a negative association for both radon (p < 10⁻¹⁴, one-tailed) and UVB (p < 10⁻³). However, models 2 and 3, which both included elevation, erased these associations (Table 4), indicating high potential for erroneous lung cancer associations when not accounting for elevation.

Confounding effect of elevation

Since elevation commonly covaries and its effect size on lung cancer is large, the potential for confounding is high. We identified two reported lung cancer associations—radon (Cohen, 1995) and UVB (Hayes, 2010)—that we attributed wholly to elevation in our analyses. Previously, Lagarde & Pershagen (1999) implicated ecological fallacy in the appearance of weak inverse associations between radon and lung cancer in Sweden. However, the plausibility of ecological fallacy resulting in the strong inverse association Cohen (1995) observed across American counties is less clear. Alternatively, Van Pelt (2003) attributed “some, but not all” of the Cohen (1995) radon association to elevation. Follow-up correspondences by each author revolved around the difficulty in assigning the effect wholly to elevation or radon when both of these highly-correlated predictors remained significant (Cohen, 2004; Van Pelt, 2004). We believe that our data quality improvements, including county-specific smoking prevalences and population-weighted elevations, were responsible for wholly attributing the effect to elevation.

Studies where lung cancer is incorporated as a predictor rather than outcome may also be susceptible to confounding. Since the collection and availability of tobacco data has historically lagged behind lung cancer data, lung cancer rates have often been adopted as a proxy for smoking prevalence (Peto et al., 1992). Ezzati et al. (2012) evaluated whether elevation was associated with several mortality outcomes across US counties. Their study relied on “lung cancer as the indicator of accumulated population exposure to smoking. This adjustment for lung cancer in multivariable regressions may have over-adjusted, if altitude has a beneficial effect on lung cancer.”

The confounding potential of elevation extends to any analysis of lung cancer across an elevation-varying region. To protect against this uncontrolled confounding, we urge future lung cancer studies to strongly consider adjusting for elevation. Fortunately, elevation is well-documented across the globe, and many existing epidemiological datasets contain locality information.

Limitations & future directions

Cross-sectional study designs are susceptible to uncontrolled confounding where associations arise due to an unmeasured confounding factor. We designed the study to minimize this risk by focusing on a homogeneous and elevation-varying region, filtering error-prone counties, including established covariates, and performing multiple stratifications and subgroupings. Furthermore, while environmental correlates of elevation represented likely confounders, all seven pervasive environmental factors we investigated could not replace elevation in models of lung cancer. Therefore, our findings gave no indication of uncontrolled confounding for lung cancer. As relevant data becomes available, follow-up across different regions will provide additional assessment of uncontrolled confounding.

Since we evaluated counties rather than individuals, ecological fallacy was also a concern. However, several of our methods and findings limit this possibility. By focusing on US counties, the smallest population grouping with systematic data available for the target region, we inherently reduced the risk of ecological fallacy. Moreover, further increases in group specificity through population subgrouping did not alter the elevation-lung association. Additionally, we reason that an exposure affecting only a portion of the population must confer a very large risk to produce the strong association observed. We find it unlikely that an extremely damaging risk factor on the individual level, such as smoking, would have evaded detection until now. Therefore, we speculate the causal factor is likely mild in carcinogenicity but universal in exposure and thus amenable to translation from the individual to population level. Follow-up biological and experimental analyses will be critical to understanding the causal factor and potential mechanisms underlying the observed elevation association. If future research confirms oxygen-driven tumorigenesis in the human lung, the present study will join the substantial list of ecological analyses that spurred new insights into cancer etiology (Pearce, 2000).

Open data

This study was made possible by excellent county level resources, many of which have only recently become available. While the trend towards available, accessible, and reusable data is encouraging, barriers still remain. For example, the 2006–2010 State Cancer Profiles release is missing incidence for two thirds of the counties in Washington due to “state legislation and regulations which prohibit the release of county level data to outside entities.” As the world comes online and the number of people with access to informatics tools expands, we see the spread of open data as a vital catalyst for progress.