Study area

We conducted this study among Geisinger’s adult patients58, located in over 40 counties in central and northeast Pennsylvania in a region with a range of UNGD activity. Geisinger’s primary care population is representative of the general population of the region based on distributions of age, sex, and race/ethnicity11. Geisinger has had a fully-operational EHR installed since 2005. All Geisinger patients had the option to opt out of all research, but less than 0.1% did so at the time of the study. Those that did not opt out were informed that their EHR data could be used for research.

Survey design and study population

Our study population came from a cohort originally designed to study the epidemiology of chronic rhinosinusitis and related nasal and sinus symptoms18,58. The survey design oversampled racial/ethnic minorities and people more likely to have nasal and sinus symptoms. In April 2014, a baseline questionnaire, a cover letter that explained the study, and a $1 bill as incentive was sent to 23,700 adults 18 years of age and older. The cover letter explained that study participation was voluntary and that if participants did not return the questionnaire they would not be included; by returning the questionnaire participants provided informed consent. Of the 23,700 letters sent, 7,847 participants responded (response rate = 33.1%)58. Six months later, a follow-up questionnaire, which included the eight-item Patient Health Questionnaire-8 (PHQ-8), was sent to all responders of the baseline questionnaire, of whom 4,966 responded (secondary response rate = 63.3%). Follow-up questionnaires were received from November 2014 to May 2015 (median of November 12, 2014). Using previously described methods70, we geocoded study subjects to their residential address listed in the EHR, 89.1% to street address, 3.1% to ZIP + 4, and 7.7% to ZIP code centroid. After excluding respondents living outside Pennsylvania (n = 34), the analysis consisted of 4,932 participants. All study protocols were reviewed and carried out in accordance with guidelines approved by the Geisinger Institutional Review Board.

Outcome ascertainment

Depression symptoms

The PHQ-8 asked individuals to report symptoms in the prior two weeks, for example, “how often have you been bothered by feeling down, depressed, hopeless?”56. Each question on the PHQ-8 has response options: “not at all”, “several days”, “more than half the days”, or “nearly every day”, scored as 0–3 respectively. Over 90% of participants answered all eight questions, for whom we calculated their total score by summing scores from the eight questions57. For participants who answered fewer than eight questions, we calculated their total score as a pro-rated sum using the formula: (sum of answered questions × 8)/(number of questions answered). Of participants who answered 1–7 questions, 81% answered 7 questions. We defined current depression symptoms using the total PHQ-8 score based on previously established categories56, but combined the two most severe groups because few participants had a “severe” total score. Scores were categorized into 0 to <5, no significant depression symptoms; 5 to <10, mild depression symptoms; 10 to <15, moderate depression symptoms; and 15 to 24, moderately severe/severe depression symptoms57. We excluded participants who did not answer any PHQ-8 questions (n = 170) from statistical analyses.

Disordered sleep diagnoses

Disordered sleep diagnoses (case-events) among the study population were identified in Geisinger’s EHR from January 2009 to June 2015. We identified encounters (98% outpatient) in the EHR that were accompanied by ICD-9 codes for disordered sleep (see Supplementary Table S1). We also identified orders for disordered sleep medications, using drug class “hypnotics” as well as drug subclass and name. We included all medications in the drug subclass antihistamine hypnotics, selective melatonin receptor agonists, hypnotics – tricyclic agents, and orexin receptor antagonists. In the subclass non-barbiturate hypnotics, we included all medications except midazolam hydrochloride, which is more often used for procedural sedation. We considered either an appropriate medication order or an encounter with a disordered sleep ICD-9 code as a disordered sleep outcome. We only retained disordered sleep diagnoses from when the participant was 18 years of age or older and randomly selected one disordered sleep diagnosis per participant per year so that study subjects with many encounters for sleep disorders would not unduly contribute (see Supplementary Figure S1).

For control dates, we identified all their dates of contact with the health system, excluded contact dates within one year of a disordered sleep diagnosis, and randomly selected one encounter date per year per participant. Control dates were frequency-matched to cases on age category (i.e., 18–44, 45–61, 62–74, 75+ years), sex, and year. We used encounter dates, rather than patients, to match as controls because of the time-varying nature of UNGD and many covariates.

Well data and activity metric assignment

We compiled well data from the Pennsylvania Department of Environmental Protection, the Pennsylvania Department of Conservation and Natural Resources, and SkyTruth, as described previously11,15,18. These data, collected for all unconventional natural gas wells in Pennsylvania from 2005–2015, included: latitude and longitude; dates of well pad construction, drilling, stimulation, and production; total well depth; and volume of natural gas produced biannually or annually.

We assigned UNGD activity for the four phases of well development (pad preparation, drilling, stimulation, and production) to each study subject (in the depression symptom analysis) or index date (in the disordered sleep analysis) using metrics that incorporated distances from participant residence to wells, and the density and size of wells, as in prior studies11,15,18. The metric has the potential to incorporate a variety UNGD-related hazards that exist on different temporal and spatial scales (e.g., regional air pollutants, local noise, truck traffic, activities that may lead to stress)71. We calculated the metric for each phase of well development:

$${\rm{Metric}}\,{\rm{for}}\,{\rm{participant}}\,{\rm{j}}=\sum _{t=t}^{d}\,\sum _{i=1}^{n}\frac{{s}_{i}}{{m}_{ij}^{2}}$$

where d was the date of return of the questionnaire, n was the number of wells in the given phase, m ij 2 was the squared-distance (meters) between well i and participant j, and s i was 1 for the pad production and drilling phases, total well depth (meters) of well i for the stimulation phase, and daily natural gas production volume (m3) of well i for the production phase. For the depression symptom analysis, for each phase of development, the metric was summed for the 14 days prior (t = d − 14) to the date of the returned follow-up questionnaire (d). We chose 14 days prior to the survey return because the PHQ-8 ascertained depression symptoms over the past two weeks. For the disordered sleep analysis, we summed the UNGD activity for the three months prior to the date of the sleep disorder diagnosis (t = d − 90). In both analyses, we z-transformed the activity metrics for each of the four phases of development, summed the transformed values, and calculated quartiles of the sums to create a composite UNGD metric that represented very low, low, medium, and high residential proximity to more and bigger UNGD wells.

Covariates

Using the EHR and questionnaire, we created covariates for potential confounding variables: race/ethnicity; sex; Medical Assistance (a means-tested program used as a surrogate for family socioeconomic status)72; age at time of questionnaire return, disordered sleep diagnosis, or control date; smoking and alcohol use status; and body mass index. Antidepressant medication use in the month prior to survey return was ascertained with medication orders based on drug group (e.g., antidepressants), class (e.g., selective serotonin reuptake inhibitors [SSRIs]), sub-class, and name. We evaluated antidepressant use as an effect modifier; we hypothesized that antidepressant medication could attenuate associations of UNGD with depression symptoms. Time-varying-covariates (all but race/ethnicity and sex) were assigned before the date of questionnaire return (for the depression symptom analysis) or before the disordered sleep diagnosis or comparison date (for the disordered sleep analysis). Based on the participants’ geocoded coordinates, we assigned them to a community using a mixed definition of place (township, borough, or census tract in cities)70. For each community, we used the 2006–2010 American Community Survey to calculate community socioeconomic deprivation73. We used public water supplier service areas from the Pennsylvania Department of Environmental Protection to assign residential water source (municipal water or ground water)74. Patients residing in homes outside the public water supplier service area were assumed to use ground water.

Statistical analysis

We employed sampling weights to estimate unbiased measures of association while accounting for the survey stratified sampling design, the response rate to the baseline questionnaire, and loss to follow-up from the baseline to the follow-up questionnaires (see Supplementary Table S2). Because one weight was much larger than the others, we truncated the largest weight to the next largest for our primary analyses75.

To build models, we first included the UNGD variable representing residential proximity to more and bigger wells (quartiles: very low, low, medium, high), and then added potential confounding variables identified a priori: race/ethnicity (non-Hispanic White, non-Hispanic Black, Hispanic), sex (male, female), Medical Assistance (no, yes), age (years), smoking status (never, former, current), alcohol use status (heavy [based on the Centers for Disease Control definition of heavy drinking as 8 or more drinks per week for females and 15 or more drinks for males76]; vs. not heavy [which included no alcohol use]), body mass index from the EHR at the visit closest to questionnaire return (BMI, kg/m2), community socioeconomic deprivation, and water source (municipal water, well water). We centered the continuous covariates (age, BMI, and community socioeconomic deprivation) and included them as linear and quadratic terms to allow for non-linearity. We did not include community type (i.e., township, borough, or city) in final models because it may lie on the causal pathway between UNGD and sleep and mental health. We used a 2-sided type 1 error rate of 0.05 for significance testing and Stata version 11.2 (StataCorp Inc.) and R version 3.2.2 (R Foundation for Statistical Computing) for analyses.

We fit multinomial logistic models to estimate the association UNGD with each level of depression symptoms (mild, moderate, moderately severe/severe) compared to no depression symptoms (reference outcome). We also evaluated the association of UNGD with depression symptoms using negative binomial regression, which treated the PHQ-8 score as a continuous outcome, allowing us to evaluate associations between UNGD and the continuous burden of depression symptoms, rather than with the screening tool’s categories77. To assess the association of proximity to more and bigger UNGD wells with disordered sleep diagnoses, we fit a survey-weighted generalized estimating equations model, to account for multiple diagnoses within participants. Because sleep diagnoses spanned 2009–2015 but most UNGD began after 2010 and attitudes towards UNGD may have changed over time, in a secondary analysis we restricted the sleep analysis to diagnoses between January 2014 and June 2015. To test a definition with a higher positive predicative value for sleep disorder, we also ran a model where a case had to have both a diagnostic code for a sleep disorder and a sleep-related medication order.

We hypothesized reduced susceptibility to stressors like UNGD among participants taking antidepressants78,79. To test this, we evaluated effect modification by antidepressant use by including cross-products of the UNGD variables and antidepressant medication use to our final depression symptom models. We used a Wald test to evaluate the significance of the cross-products.

In a sensitivity analysis, we evaluated the influence of our sample weights by examining associations of UNGD with depression symptoms among all subjects using the final multinomial logistic and negative binomial models without weights and with full and truncated weights. In theory, weighted models will provide less precise, but more unbiased estimates than unweighted models80, our rationale for using truncated weighted models as the primary analysis.

Data Availability

Data on unconventional natural gas development in Pennsylvania are publicly available from the Pennsylvania Department of Environmental Protection, the Pennsylvania Department of Conservation and Natural Resources at http://www.dcnr.state.pa.us/topogeo/econresource/oilandgas/resrefs/wis_home/. The health data that were used in this study are protected health information and subject to many restrictions. Data may be available from the corresponding author upon reasonable request and with specific required agreements in place.