Fluorosurfactants in groundwater increase the incidence of chronic health conditions among California Medicare beneficiaries

Health Outcome Data

We obtained individual-level information on chronic health conditions among Medicare beneficiaries from the Centers for Medicare and Medicaid Services’ Chronic Conditions Warehouse. For each Medicare beneficiary, these annual files provide indicators for 27 chronic health conditions as well as their annually updated ZIP-code of residence (US Department of Health and Human Services 2023). Since Medicare is exceptional for those under 65 years, we restricted our analysis to Medicare beneficiaries aged 65 years and older for whom coverage is nearly universal.

The 27 chronic health conditions were separated into nine common conditions (cataract, chronic kidney disease, hypertension, congestive heart failure, hyperlipidemia, anemia, diabetes, rheumatoid/osteo arthritis, ischemic heart disease), eight less common conditions (chronic obstructive pulmonary disease (COPD), stroke, atrial fibrillations, glaucoma, hyperplasia, acquired hypothyroidism, osteoporosis, Alzheimer’s, and dementia), and three negative controls (lung cancer, depression, and hip fracture). We excluded from the analysis seven outcomes that averaged less than one incident case per year in a ZIP code: prostate cancer, breast cancer, endometrial cancer, colorectal cancer, asthma, Alzheimer’s, and acute myocardial infarction. We used indicators of each chronic health condition to calculate the annual incidence of newly diagnosed conditions over the seven-year period (2011-2017) for each ZIP code of residence.

For each ZIP code of residence, each unique beneficiary with an incident chronic health condition at any time during the seven-year study period contributed one count to the incidence of the corresponding condition. Person-years at risk were calculated in a similar manner; for each ZIP code of residence, each unique beneficiary contributed a beneficiary-year at risk for every year that they did not have the condition as well as the year in which they were first diagnosed.

Medicare beneficiary residential ZIP codes (postal delivery routes assigned by the US Postal Service) were mapped to both ZIP code population weighted centroids (US Department of Housing and Urban Development 2023) and ZIP Code Tabulated Areas (ZCTAs) (US Bureau of the Census 2023). In this process, multiple ZIP codes may have been mapped to the same ZCTA (2,444 ZIP codes were mapped to 1,756 ZCTAs).

PFAS Groundwater Concentrations

PFAS concentrations in groundwater were obtained from California’s Groundwater Ambient Monitoring and Assessment (GAMA) program, a collaboration between the U.S. Geological Survey and California’s State Water Resources Control Board, which began testing domestic and public supply wells for PFAS in 2018 (US Geological Survey 2023). Well location data and PFAS sampling results were downloaded from the GAMA website for the period 2018 until 2021, however most samples were analyzed in recent periods (California Water Boards 2023).

Individual wells may have been tested several times for multiple PFAS species and all results were reported in nanograms per liter (ng/L). Values reported as less than the detection limit were interpreted as non-detects and converted to zero before calculating exposure; however, minimum detection limits varied by PFAS testing method. Maximum PFAS concentrations at individual wells were aggregated to the ZCTA-level as an inverse squared distance weighted concentration averaged across all wells within 20-km of the ZCTA population centroid, with wells closest to the centroid contributing the most weight.

Sociodemographic Data

Population characteristics, including age, race, and income, were downloaded from the US Census Bureau’s 5-year American Community Survey (ACS-5) for all ZCTAs and census tracts in California for the year 2015 (US Bureau of the Census 2023). Additional sociodemographic data were retrieved from California’s Office of Environmental Health Hazard Assessment’s CalEnviroScreen 4.0 (CES) dataset which is an interactive mapping tool and database of environmental, socioeconomic, and public health data (California Office of Environmental Health Hazard Assessment 2023). The CES data were available by census tract, and for the purpose of the analysis they were averaged for all census tracts within each ZCTA as well as linked to the wells found within the census tract.

PFAS Source

Locations of industries such as airports, fire training facilities, chemical manufacturing sites, and waste management sites registered through USEPA’s Enforcement and Compliance History Online (ECHO) database and included in EPA’s PFAS Analytics Tools were treated as potential PFAS sources in this study (‘PFAS Sources’, 23 variables, Table S1) (US Environmental Protection Agency 2023). Although, except for most military installations, this dataset does not include any actual data or confirmation about each facility’s PFAS use or discharge, these industry points are considered to have some relevance to PFAS use or discharge. Potential PFAS sources were counted within well site buffers (1- and 5-km) and ZCTAs by industry sector, and distances to the nearest PFAS industry source were measured from well sites. The number of PFAS industry sources were included as predictors of the exposure and were not considered to be on the pathway between PFAS concentration and health outcomes.

Landcover Data

Land use and landcover data were obtained from the Multi-Resolution Land Characteristics Consortium’s 2016 National Land Cover Dataset (NLCD) (Dewitz and Survey 2021, Wickham, Stehman et al. 2021). The NLCD uses Landsat imagery to classify landcover into 16 categories at a 30m resolution for the entire US. The NLCD Land Cover Change Index raster dataset, which measures changes in landcover types between 2001-2019, was also included in this analysis. The NLCD also provides a raster of urban imperviousness that assigns each 30m pixel a value of percent developed impervious cover, for example, pavement and rooftops. Landcover data were summarized within California’s ZCTAs as well as within 1- and 5-km buffers around each well site. Local soil characteristics, including percent organic matter, percent clay, and pH, were summarized at each well site with data from the POLARIS database (Chaney, Wood et al. 2016, Moro Rosso, de Borja Reis et al. 2021). Elevation (meters above sea level) was estimated at each well site using the “elevatr” package in R (Hollister, Shah et al. 2022).

Statistical approach

We used Poisson regression (“glm” package in R) (R Core Team 2022) to model the incidence of chronic health conditions as a function of PFAS concentrations (ng/L) in the groundwater. We first restricted study area to regions where at least 25 percent of drinking water was sourced from groundwater wells. This restriction was not based on PFAS groundwater concentrations nor incidence of chronic health conditions. However, to control for potential selection bias and for confounding by common causes introduced by choice of drinking water source we developed a propensity score (‘GW PS’). Second, we consider the possibility that wells were not tested at random and may have been differentially sampled with respect to unknown factors. To adjust for selection bias due to preferential sampling of wells, we also developed a well testing propensity score (‘Testing PS’) defining the probability of a well being tested based on measured covariates. Both propensity scores were included in the model as natural splines with three degrees of freedom.

Additionally, three population-level variables that are known determinants of health were added to the regression (‘Regional Determinants of Health’; percent of workforce employed in agriculture, percent Hispanic ethnicity, and percent with a disability). The prevalence of smoking (‘Smoking’) was unknown, both for the general population and for Medicare beneficiaries. Smoking is known to be the overwhelming cause of lung cancer and any contribution of PFAS must be orders of magnitude less, so we used lung cancer incidence as a surrogate to account for the unmeasured variation in the prevalence of smoking across the ZCTAs in the study area.

Finally, we adjusted for two sets of covariates that may predict the exposure: covariates related to the known and potential sources of PFAS (‘PFAS Sources’, 23 variables, Table S2) and covariates that were indicative of PFAS concentrations (‘PFAS Covariates’, 10 variables, Table S2).

We also considered depression and hip fracture, as negative controls to assess residual confounding by unknown or omitted sources of confounding. These outcomes have no known biological mechanisms that link the incidence of disease to PFAS exposure. However, we expect there may be other factors associated with negative controls and PFAS that are also likely to confound the association with other outcomes of interest. Therefore, the common causes of exposure to PFAS and outcome are hypothesized to be like the common causes of exposure to PFAS and incidence of chronic conditions. These factors may be directly or indirectly measured through variables describing the residential location, type of employment, proximity to sources, demographic characteristics, etc. Since social isolation is the major risk factor for depression, we believe that there exists no biological relationship between PFAS and depression. Any residual association between PFAS and depression would be evidence of residual confounding. We did not find literature reporting a direct link between hip fracture with PFAS exposure. The major risk factor for hip fracture is fall and age, but factors such as arthritis and osteoporosis may contribute to the severity of outcomes associated with hip fracture. Those outcomes are analyzed separately.

Propensity score for PFAS testing

Groundwater wells were not randomly selected for PFAS testing which can potentially lead to selection bias in the analysis. To account for potential selection bias, we estimate the probability of each well being selected for PFAS testing (testing propensity score). We take advantage of the large and comprehensive characterization of residential and potential PFAS sources within 1-km and 5-km radii of each well and use extreme gradient boosting, XGBoost algorithm (Chen and Guestrin 2016), to predict the probability of a well being tested for PFAS. More specifically, the following variables were summarized at the well (buffer) level and used in the XGBoost algorithm: the Californian hydrologic region (https://gis.data.ca.gov/datasets/2a572a181e094020bdaeb5203162de15/explore) that the well is in, land use and landcover variables from the NLCD, location of potential industrial PFAS sources (EPA PFAS Analytics Tool), population characteristics (CES and US Census, ACS), soil characteristics (POLARIS), and elevation.

The XGBoost model was fit in two stages using the test-set validation approach. The model in the first stage was fit on a random subset containing 90% of wells, with the optimal number of gradient boosting rounds determined via early stopping with ten rounds, using the other 10% of wells as a validation set. That is, the optimal round is the last round after which performance on the validation set doesn’t improve for ten rounds. The model in the second stage was then fit on all wells, with the optimal number of rounds found in the first stage. The other parameters of the model are set to the defaults as in the “xgboost” package in R (https://xgboost.readthedocs.io/en/stable/parameter.html). This model was then used to predict propensity scores for all the wells. We aggregated propensities to the ZCTA-level as an inverse distance weighted concentration averaged across all PFAS-tested wells within 20-km of the ZCTA population centroid, as with PFAS concentrations.

Propensity score for drinking water from ground sources

We consider groundwater consumption in the epidemiological model for PFAS in groundwater. While we estimated 22% of Californians rely on groundwater as their primary source of drinking water, this proportion varies considerably by county. In the study area, the overall proportion was 54%. Data on public water systems, including population served and whether sourced from groundwater or surface water, were downloaded from California’s Public Water System Information Dataset (PWSID) (California State Water Resources Control Board. 2019). These water-source estimates were combined with estimates of population served by private drinking water wells (Johnson, Belitz et al. 2019) to calculate percent of county population drinking from groundwater. We restricted the study region to ZCTAs with population centroids within the 40 counties where 25% or more of the population rely on groundwater for their drinking water.

To account for potential selection bias from restricting our analysis to the specified ZCTAs, we calculate the second propensity score for each ZCTA (probability of selection in the study region). We use extreme gradient boosting (Chen and Guestrin 2016) to calculate the propensity scores, as earlier. However, for the two-stage test-set validation approach, the model was trained on 80% of the ZCTAs and validated on the other 20%; the validation set has a larger proportion to adjust for the smaller number of ZCTAs. We predicted these propensity scores conditional on ZCTA-level variables: land use and landcover variables (NLCD), location of potential industrial PFAS sources (EPA PFAS Analytics Tool), and population characteristics (CES and US Census, ACS).

Statistical models

For each outcome of interest, we fit seven models (M1-M7) for weighted PFAS concentrations as a treatment variable of interest. M1 (Base Model) used the log count of the total Medicare beneficiaries (population-at-risk) as an offset. The subsequent models were cumulative adjustments to the basic model: M2 (+ Testing PS) adds the PFAS-well testing propensity score to control for potential selection bias; M3 (+ GW PS) adds the propensity score for the proportion of drinking water from groundwater; M4 (+ Regional Determinants of Health) adds three major regional health determinants for the total population, including percent working in agriculture industry, percent Hispanic and percent with disability; M5 (+ Smoking) adds lung cancer incidence as a surrogate for smoking prevalence. Finally, model M6 (+ PFAS Sources) adds known industrial PFAS sources (Supplemental Table S1) and M7 (+ PFAS Covariates) adds additional covariates (not including sources) that were related to high levels of PFAS (Supplemental Table S2).

Sensitivity analysis was conducted with respect to the exposure metric. In addition to the weighted metric used in the main analysis, we also considered unweighted, and propensity weighted metric. The propensity weighted metric was constructed as , where w_i is the distance weight of the well to the population centroid of ZCTA, S_i is the indicator of a sampled location, and p_i is the propensity score of the sampling and the summation is over sampled locations in D_S. Rescaling the original propensity score , where ) provides a property of . As a result, the constructed propensity weighted metric approximates the weighted average of PFAS had all the wells in the ZCTA D_S been sampled . This metric estimates the total exposure for this ZCTA and adjusts for preferential sampling by weighting.

Results are reported as the percent difference in the annual incidence of a chronic condition for a 10.0 ng/L increment in estimated PFAS groundwater concentrations, which corresponds to the 80^th percentile of the PFAS distribution in the study area. To recapitulate, PFAS groundwater concentrations are estimated for each ZCTA as the inverse distance-squared weighted mean of the well-specific maximum PFAS-species concentration for tested wells within 20-km of the ZTCA population centroid.

Strengths and weaknesses

This study has several strengths and weaknesses. For those over the age 65, Medicare coverage is essentially universal. Thus, data collected on new incidence of chronic disease is representative of the population in the same age group in the study region. A unique strength of Medicare beneficiary summaries is that the chronic diseases are tracked over all beneficiaries and the incidence and prevalence rates can be fully enumerated. Nevertheless, to our knowledge there have not been other studies of PFAS effects using beneficiary files.

In this analysis we payed close attention to various sources of potential bias including preferential sampling of PFAS and confounding. We carefully considered selection bias from the non-random sampling of wells for PFAS testing, but the adjustment by a testing propensity score showed minimal impact on the PFAS effect estimates. The restriction of the study area to California counties where more than 25% of drinking water comes from groundwater was also adjusted by a propensity score for groundwater use and several demographic characteristics associated with the selection. The adjustment had minimal impact on the PFAS effect estimates (Table 2), but the reported findings ought not to be extrapolated to areas that were excluded from the analysis. For example, the region around Los Angeles has a number of wells with detected PFAS levels however, the region is estimated to have low reliance on the use of ground water sources for drinking. Additionally, some of the outcomes that had the largest risk estimates, such as cataract, have previously been noted to have both genetic and age dependent risk factors. Although conceivable that an association exists between average age in population and average PFAS concentration and that the outcome of interest is a surrogate for age, the risk estimates did not change after adjusting for age.

PFAS sources had a small suppression effect (MacKinnon, Krull et al. 2000), which when added to the model increased the effect estimates across most outcomes. The small explanatory power of PFAS sources suggest that presence of sources cannot act as an Instrumental variable and therefore cannot be a substitute for the direct exposure measurements. However, the addition of these variables can increase the magnitude and improve the overall model.

This study is a cross sectional analysis of chronic condition incidence and PFAS exposure in ground water used for drinking water, and as such is subject to several sources of potential confounding. The major confounders in the analysis were the three regional determinants of health (Hispanic ethnicity, disabilities, and agricultural employment) and the incidence of lung cancer as a surrogate for tobacco smoking (Table 2). Adjustment for the three regional variables increased the PFAS associations, while the adjustment for lung cancer reduced these associations. Among other limitations, smoking status of beneficiaries was not available and lung cancer incidence is used to characterize population level rates of smoking. Tobacco smoking does have a strong independent association with hypertension and partial adjustment for smoking using lung cancer incidence shows a large change in estimate, reflecting considerable confounding. The adjustment by the propensity scores for groundwater use showed less evidence of confounding. The widths of the confidence intervals for each chronic condition do not change much across the epidemiologic models (Table 2), providing no evidence of overfitting. This lack of variance inflation persists even for meaningful changes in the point estimates.

This study is an ecological study of a shared environmental exposure within the ZIP code of residence. Population level assessments are challenging due the lack of a comprehensive sampling strategy for PFAS. In this study, we used measurements for PFAS in groundwater from GAMA data repository of various sampling campaigns. While possible, measurement errors are unlikely to have varied across ZCTAs differentially with the incidence of chronic health conditions. Similarly, the clinical diagnosis of chronic health conditions supplied by the Medicare data is unlikely to have been influenced by PFAS groundwater concentrations. Finally, temporal misalignment between sampling of PFAS, exposure, and measured outcomes in these early stages of population level investigation is a knowledge gap and should not be overlooked.

Future research will benefit from enhanced measurements of PFAS in drinking water sources that is ongoing in many States. In present work we use the best available data on groundwater concentrations, but measures of drinking water exposure are likely to provide direct linkage. We found that adjustment for area-level lung cancer incidence was an important confounder, which implies that smoking is likely to be an important factor for consideration in future studies. The second most influential confounders included areal-level determinants of population health, which in the study area were % hispanic population, % disability, and % employed in agriculture. Adjustment with propensity scores for the non- random sampling of wells for PFAS testing and for differences in groundwater use for drinking water produced modest changes in the estimates but are important to consider in studies of ground water exposure. The impact of groundwater use was minimized by our restriction to California counties with at least 25% of drinking water derived from groundwater. Finally, the range of chronic conditions associated with PFAS indicated non-specific modes of action. Separation of PFAS associations into direct and indirect effects would be useful in resolving these pathways, as shown by our analysis using hypertension as a mediator.