Associations between air pollution and psychiatric symptoms in the Normative Aging Study

We studied men enrolled in the VA NAS. In brief, the NAS study is an ongoing longitudinal closed-cohort study of aging established by the US Department of Veterans Affairs. The participants were recruited initially between 1961 and 1970, and consisted of 2280 community-dwelling men who lived in the greater Boston, Massachusetts area, aged 21–81 years at the time of enrollment (Bell et al 1972). The majority of them were Caucasian men. At enrollment, men with known chronic medical conditions (i.e. history of cardiovascular diseases, hypertension (HT), cancer and respiratory diseases) were excluded. Participants underwent examinations at approximately three to fiveyear intervals, including routine physical examination and laboratory tests, and self-reported information on medical history, drinking or smoking history, physical activities and other demographic characteristics.

We studied members of the cohort who completed the BSI. Details of the administration of BSI in the study can be found in the online-only text (available online at stacks.iop.org/ERL/17/034004/mmedia). We focused on the study period of 2000–2014 when both BSI records and geographical fine-resolution ambient air pollution predictions were available, leaving us with 764 NAS individuals. We further excluded 145 people without sufficient response to BSI items (<41 items answered, as recommended by the BSI scoring manual (Derogatis 1993)). Another 49 individuals who were on psychiatric medication use during the visits (including anti-anxiety or antidepressant drug use) were also excluded, limiting the bias from the influence of psychiatric medication use on the response of psychiatric symptoms. All participants were sent BSI questionnaires, however some failed to complete it. The final study sample consisted of 570 men with 1114 visits collected over the period of 2000–2014.

We examined three air pollutants in this study: PM_2.5, ozone, NO₂. The selected air pollutants are the three major air pollutants in the US with fine resolution prediction data available at daily and 1 km × 1 km grid cell level, which is one of the strengths of this study. National fine-resolution ambient air pollution spatiotemporal predictions were estimated by combined machine learning algorithms using a geographically weighted regression. The algorithm fused ground monitoring data (from the Air Quality System operated by the US Environmental Protection Agency and other monitoring data sources), satellite-derived measurements, 32 km × 32 km meteorological conditions retrieved from the North American regional reanalysis data set (air temperature, humidity, wind speed, etc), chemical transport model simulations, and 1 km × 1 km land-use data (land-use coverage types, road intensity, elevation and Normalized Difference Vegetation Index (NDVI), etc over 2000–2016 (Di et al 2019, 2020, Requia et al 2020). Cross-validation comparing ambient predictions and held out monitored values indicated good performance of the prediction models (ten-fold cross validation: PM_2.5: R² = 0.86; O₃: R² = 0.91; NO₂: R² = 0.79). The predictions were spatially available at each 1 km × 1 km grid cell covering the contiguous US on a daily basis from 2000 to 2016. Daily ambient surface temperature and relative humidity were obtained from the gridMET project predictions using NASA's Land Data Assimilation System (NLDAS-2) with approximately 4 km × 4 km spatial resolution (ClimatologyLab 2021). The gridMET project blends spatial attributes of gridded climate data from the parameter-elevation relationships on independent slopes model data with desirable temporal attributes from regional reanalysis (NLDAS-2) to provide high-resolution (1/24th degree ∼4 km) gridded datasets of surface meteorological variables in the contiguous US. Validation of the resulting data was conducted against an extensive network of weather stations and achieved good performance (Abatzoglou 2013). We matched predicted ambient air pollutants, temperature and relative humidity based on the exact residential address at the visit dates, as well as different time periods before the visit dates, by averaging the daily predictions for the various time windows of interest. The exposure time windows were averages of one week (1 wk), four weeks (4 wk), eight weeks (8 wk) and one year (1 year) prior to the visit. Residential history was used to assign the appropriate address for each of the prior days before computing the different time averages.

The BSI is a brief self-report psychiatric symptom scale (Derogatis and Melisaratos 1983) applied to both clinical psychiatric patients and non-patient general population. It consists of 53 items covering nine symptom dimensions: somatization, obsession-compulsion, interpersonal sensitivity, depression, anxiety, hostility, phobic anxiety, paranoid ideation and psychoticism; as well as three global indices of distress: global severity index (GSI), positive symptom distress index (PSDI), and positive symptom total (PST) (Derogatis and Melisaratos 1983). Respondents rank each feeling item out of the total 53 items (e.g. 'your feelings being easily hurt') on a five-point scale ranging from 0 (not at all) to 4 (extremely). Rankings characterize the intensity of distress during the past 7 d. Higher scores indicated worse symptoms experienced. The GSI was calculated using the sums for the nine symptom dimensions plus the four additional items not included in any of the dimension scores, and dividing by the total number of items to which the individual responded including answers of zero score. It represented the respondent's global distress level. The PSDI was the sum of the values of the items receiving non-zero responses divided by the total number of non-zero response positive symptoms. This index provided information about the average intensity of distress the respondent experienced for the positive symptoms he had. The PST was the total number of non-zero responses. Computation of the scores were conducted following the BSI scoring manual (Derogatis 1993). Standardized T scores for GSI and PSDI were additionally computed based on their raw scores following the scoring manual. The BSI was previously proved to work as a general indicator of psychopathology with good validity and reliability (Boulet and Boss 1991).

In this study, the outcomes of interest are the three BSI global indices. Research with the global indices has confirmed that the three scales reflect distinct aspects of a psychiatric disorder. They measure level of symptomatology (GSI), intensity of symptoms (PSDI), and number of reported symptoms (PST), respectively (Derogatis et al 1975).

For each study visit, we obtained basic demographic information (age, race, ethnicity, education, marital status, etc), physical health indicators body mass index (BMI, any major illness, cardiovascular disease, respiratory disease, HT, diabetes, total and high-density cholesterol levels, etc), physical activity (via information on total metabolic task equivalent level), smoking, drinking, etc. Area-level contextual characteristics, such as census tract-level population density, median value of house and median household income, were extracted from the US census (US census decennial data 2000/2010, American Community Survey for years after 2010) and further linked with the participant's address over the study period. Continuous variables include age (years), BMI (kg m⁻²), census tract median household income (USD/household, year), census tract median value of house (USD/unit), census tract population density (persons/mile²) and physical activity as measured by total metabolic equivalent of task (hours/week). Categorical variables include calendar year of visit, season of visit (spring = March–May, summer = June–August, autumn = September–November, winter = December–February), race (Caucasian or not), education (none, grade school, some high school/college, graduate, professional), marital status (married, single, separated, divorced, widowed, other), smoking (never, current or former cigarette smoker) and drinking (usually have more than two drinks per day or not).

We applied generalized additive mixed effects regression with subject-specific random intercepts to examine the associations between air pollution exposure and the levels of psychiatric symptoms. The mixed model takes into account the correlation between repeated measures collected for the same individual at different study visits. BSI GSI T scores, PSDI T scores as well as PST were analyzed as the outcome variables in separate models. Air pollution measures were fitted as predictors in the mixed effects models using averages of one week, four weeks, eight weeks and 1 year prior to the visit date to explore the associations at different exposure windows. For each exposure window and outcome of interest, we included all three air pollutants (PM_2.5, ozone, NO₂) in a multi-pollutant model simultaneously.

We used substantive knowledge (expert knowledge, literature evidence and directed acyclic graph (figure S1)) to identify the confounding structure, and hence the covariates to adjust for in the models. Details of the discussion on the confounding structure are provided in online-only text in the supplement. Two parallel sets of multi-pollutant models were conducted as the main analyses (a) adjusting for a key set of covariates that we determined via the variable selection procedure (temperature, relative humidity, calendar year, age, BMI, race (Caucasian or not), education, marital status, season, census tract median household income, census tract median value of house, census tract population density) (b) adjusting for the same covariates and three behavioral factors (smoking, drinking and physical activity status). The additional adjustment for behavioral factors is discussed in online-only text in the supplement. Nonlinearity of the exposure-outcome or covariate-outcome relationships were examined by fitting penalized splines of the continuous variables in the mixed effects model. If covariates or exposures followed an approximately linear relationship with the outcome, we modeled these associations as linear terms in the final models. Otherwise, the penalized spline term was used in the final models. Existing literature supported short-term temperature associations with adverse mental outcomes (Basu et al 2018, Mullins and White 2019, Li et al 2020). Therefore, we adjusted for temperature in those models. Temperature followed a close-to-linear relationship with the outcome in our study, and was thus modeled as a linear term for each window (figure S2). However, due to concerns regarding potential over-control when adjusting for season and temperature simultaneously for the window of time one year prior to a visit, we did not control for temperature for this exposure window in the main analysis.

The key effect estimates reported were different in each of the three global BSI scores (as well as 95% confidence interval (CI)) per interquartile (IQR) increase in PM_2.5, ozone, NO₂ over different time windows prior to visit. In addition, we conducted effect modification analyses by area-level contextual factors using multiplicative interaction terms between contextual factors and exposures. Specifically, we created an area-level income variable (high-income vs low-income areas) using the 75th percentile of median household income ($85 457). We dichotomized median house value (high house value vs low house value areas) using the 25th percentile (∼$213 070). Census tract population density was used to create the high vs low population density area indicator using the 75th percentile distribution cutoff of 1917 persons/mile². Sensitivity analyses were conducted for (a) single-pollutant models; (b) using a different community reference sample to create the T scores; (c) further including people on psychiatric medication; (d) additionally adjusting for temperature in the one year window prior to a visit. All the analyses were conducted in R (version 3.6.1).

A total of 570 participants with at least one BSI psychiatric symptom assessment during the study period entered into our final analyses (with a total of 1114 visits). Of them, 301 provided a second assessment, 176 of them provided a third, 66 of them provided a fourth and one provided a fifth visit. The median time interval between visits with BSI assessment is 3.96 years (IQR: 3.03–6.64 years) table 1 presents the overall characteristics for our study sample (included) as well as the characteristics for those who were still active during the study period but were excluded based on the restriction criteria (excluded). Our participants had an average age of 72.4 years. Most of them were married (74.4%) with an education level equivalent to some college at baseline. About 65% were former cigarette smokers. Since these are relatively older men, some of them had physical health morbidities, such as coronary heart diseases (CHDs), HT, stroke, etc. There is no strong evidence of substantial differences between the study sample and those who did not enter our study. However, men selected in our study sample had a lower prevalence of HT and CHD, compared with those who did not enter.

The distributional characteristics of air pollutants are shown in table 2. Our study participants had an overall average annual residential exposure level of about 9.0 µg m⁻³, 34.8 and 22.3 ppb to PM_2.5, O₃, and NO₂ respectively. We also examined the correlation between the ambient exposures in this study and observed low to medium correlation among them (figure S3).

Table 3 presents the estimates for the associations between air pollutants at different time windows and PSDI T score based on two sets of models (model I, model II). Based on model I, for each IQR increase in the exposure to the prior one or four weeks average of NO₂, PSDI T score increased by 1.60 (95% CI: 0.31, 2.89) and 1.71 (95% CI: 0.18, 3.23) respectively, controlling for covariates. Similarly, for each IQR increase in the exposure to prior one or four weeks average of ozone, PSDI T score increased by 1.66 (95% CI: 0.68, 2.65) and 1.36 (95% CI: 0.23, 2.49) respectively, adjusting for the same key set of covariates. The estimates remained robust when additionally controlling for smoking, drinking and physical activity (model II). We did not see associations between elevated exposure to air pollution and the increase in GSI T score or PST (tables 4 and 5). Results remained robust in single-pollutant models (table S1). Results also remained robust when we used different community reference samples to compute the T scores (online-only text, table S2). Further, including people on psychiatric medication did not substantially change the findings, but overall tended to attenuate the ozone and NO₂ effect estimates (table S3). We did not observe a substantial difference in the results for the average one year prior to visit window after additional adjustment of temperature (table S4).

In the modification analyses by community contextual factors on PSDI T score (table 6), we found that for the exposure windows of ozone at an average of one week and four weeks prior to the visit, area-level household income and house value modified the associations between ozone and PSDI T score reported. However, we did not observe the same modification in NO₂—PSDI relationships. For each IQR increase in one week ozone, PSDI T score increased by 1.79 (95% CI: 0.8, 2.78) in low-income areas vs −0.49 (95% CI: −2.48, 1.49) in high income areas (interaction P-value: 0.010). Additionally, for each IQR increase in one week ozone, PSDI T score increased by 3.12 (95% CI: 1.26, 4.98) in low house value areas vs 1.02 (95% CI: −0.06, 2.10) in high house value areas (interaction P-value: 0.007). Similar differential associations were found for average four weeks exposures of ozone.