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

Environmental risk factors for psychiatric health are poorly identified. We examined the association between air pollution and psychiatric symptoms, which are often precursors to the development of psychiatric disorders. This study included 570 participants in the US Veterans Administration (VA) Normative Aging Study (NAS) and 1114 visits (defined as an onsite follow-up at the VA with physical examination and questionnaires) from 2000 to 2014 with information on the brief symptom inventory (BSI) to assess their psychiatric symptom levels. Differences in the three BSI global measures (global severity index (GSI), positive symptom distress index (PSDI) and positive symptom total (PST)) were reported per interquartile (IQR) increase of residential address-specific air pollutants levels (fine particulate matter—PM2.5, ozone—O3, nitrogen dioxide—NO2) at averages of one week, four weeks, eight weeks and one year prior to the visit, using generalized additive mixed effects models. We also evaluated modification by neighborhood factors. On average, among the NAS sample (average age: 72.4 years (standard deviation: 6.7 years)), an IQR increase in one and four week averages of NO2 before a visit was associated with a PSDI T score (indicator for psychiatric symptom intensity) increase of 1.60 (95% confidence interval (CI): 0.31, 2.89), 1.71 (95% CI: 0.18, 3.23), respectively. Similarly, for each IQR increase in one and four week averages of ozone before a visit, the PSDI T-score increased by 1.66 (95% CI: 0.68, 2.65), and 1.36 (95% CI: 0.23, 2.49), respectively. Stronger associations were observed for ozone and PSDI in low house-value and low household income areas. No associations were found for PM2.5. Exposure to gaseous air pollutants was associated with a higher intensity of psychiatric symptoms among a cohort of older men, particularly in communities with lower socio-economic or housing conditions.


Introduction
An increasing proportion of people are being affected by psychiatric problems. In the US, about 7.1% adults are living with diagnosed major depressive disorders (Center TNIoMHIR 2019). Studying the risk of development of psychiatric disorders is complicated by the uncertainty of disease onset as well as diagnosis error. In contrast, studying psychiatric symptoms in a designated time window can provide valuable information on people's mental state. Acute symptoms, particularly among people with psychiatric disorders who have not received proper treatment, usually correlate with longer-term experience (Valiengo Lda et al 2016). Mood alterations in the elderly population have been associated with a disruption of their physical, mental and social wellbeing (Penninx et al 1998, Geerlings et al 2000, Abrams et al 2002, Lenze et al 2005, Zhang et al 2005, Thurston et al 2006. An increase in the frequency or intensity of any psychiatric symptom could mean a higher risk of developing clinically relevant psychiatric disorders and often serve as the underlying criteria for the diagnosis manuals of mental diseases. Growing evidence has shown that air pollution may be harmful for psychiatric outcomes. Toxicological studies based on animal models and human brain neuro-imaging also support such a link (Peterson et al 2015, Wilker et al 2015. However, most previous epidemiology studies focused on particulate matter with supporting meta-analysis evidence (Braithwaite et al 2019). In the National Social Life, Health and Aging Project study, increased exposure to fine particulate matter (PM 2.5 ) was related to depressive and anxiety symptoms in US older adults (Pun et al 2017). Similar harmful associations were found for children between acute exposure to PM 2.5 and increased risk of psychiatric emergency department visits (Brokamp et al 2019). There is much less evidence on whether gaseous pollutants (such as ozone and nitrogen oxides) can affect psychiatric symptoms or the risk of developing psychiatric disorders.
In this study, we investigated the associations between exposure to particulate matter and gaseous air pollutants with reported psychiatric symptoms, using brief symptom inventory (BSI) measures from the Veterans Affairs (VA) Normative Aging Study (NAS). We aimed to provide high quality evidence based on a longitudinal cohort with fine resolution exposure data. We hypothesized that exposure to increased levels of air pollution would be associated with increased psychiatric symptoms. We also expected the above associations to be dependent on the exposure window as well as area-level contextual factors. If these associations can be demonstrated, this study could provide new insights into the role of air pollution in the development of psychiatric symptoms. Further, if air pollution is a risk factor for psychiatric diseases, then intervention can be applied at the population scale to reduce the related disease burden.

Sample
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 fineresolution 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.

Air pollution and meteorology assessment
We examined three air pollutants in this study: PM 2.5 , ozone, NO 2 . 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), satellitederived 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, 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 2 = 0.86; O 3 : R 2 = 0.91; NO 2 : R 2 = 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 highresolution (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.

Outcome measurement
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, obsessioncompulsion, 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).

Covariates information
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 −2 ), census tract median household income (USD/household, year), census tract median value of house (USD/unit), census tract population density (persons/mile 2 ) 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).

Statistical analysis
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 2 ) in a multipollutant 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 2 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 (highincome 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 2 . 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).

Sample characteristics
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.

Exposure distribution and correlation structure
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 −3 , 34.8 and 22.3 ppb to PM 2.5 , O 3 , and NO 2 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 2 , 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 2 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).

Effect measure modification analysis
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 2 -PSDI relationships. For each IQR increase in one week ozone, PSDI T

Discussion
In this longitudinal cohort study of NAS participants during the period 2000-2014, we observed positive associations between short-medium term (average one and four weeks prior to visit) ambient exposure to gaseous pollutants (ozone and NO 2 ) and elevated psychiatric symptom intensity. No associations were found for PM 2.5 , nor on reported frequency of symptoms. Area-level household income and house value modified the associations between ozone exposure and psychiatric symptoms, with residents of lower income or house-value neighborhoods having an increased response to ozone. Previous literature looking at psychiatric healthrelated outcomes and air pollution has not shown consistent results and direct comparison is difficult as methods of exposure and outcome measurement vary widely between studies. A previous study looking into the same NAS cohort found that increased levels of air pollutants such as PM 2.5 and NO 2 measured from local stationary monitoring sites were associated with an increase in perceived stress score (Mehta et al 2015). A time-series study conducted in California looking at emergency department visits and psychiatric hospital admissions, found that, per 10 ppb increase in the 7 d average daily 8 h ozone from monitors increased the number of overall psychiatric admissions by 0.64% (95% CI: 0.21%-1.07%), depression admissions by 1.87% (95% CI: 0.62%, 3.15%), and bipolar admissions by 2.83% (95% CI: 1.53%, 4.15%) (Nguyen et al 2021). Researchers in Italy found same-day ozone levels, but not NO 2 or CO, to be associated with an increase in admissions for psychiatric conditions in a multipollutant model (Bernardini et al 2019). In that study, the exposure data were averaged across up to 16 monitor stations dispersed in the Umbria region of Italy. More evidence is needed looking at psychiatric symptoms and ambient exposures using finer geographical resolution for exposures, and in a longitudinal setting to help establish more solid temporal links between them.
The main hypothesized biological mechanism by which air pollution could lead to adverse psychiatric health outcomes is through increased inflammation and oxidative stress in the central nervous system (Ventriglio et al 2021). Animal models have shown that exposure to ozone can increase inflammatory biomarkers not only in the lungs but also in the cerebral cortex (González-Guevara et al 2014). Pollutants may move across biological membranes in the lungs to enter the systemic circulation and spread throughout the body, including the nervous system (Thomson 2019). Our results showed elevated levels of symptoms reported only on the PSDI and not the GSI or PST under increased ozone and NO 2 levels. One reason could be the specific sample that we studied. Our cohort consists of elderly, mostly Caucasian veterans. They are mostly better educated than their peers, and their psychiatric health status may differ from those whose health profiles were used to develop the BSI. In general, their BSI scores are lower than average. It might be that the GSI or PST are less sensitive in this sample as compared to the PSDI. It is also possible that air pollution only affected the intensity of the reported symptoms, but not the overall symptomatology or number of symptoms.
Previously, in the Women's Health Initiative Observational Study looking at long-term PM 2.5 and risk of Cardiovascular Disease (CVD), higher risk estimates were observed among participants living in low-SES (socio-economic status) neighborhoods, indicating people with lower neighborhood socioeconomic status (nSES) may be more susceptible to air pollution-related health effects (Chi et al 2016). That susceptibility could stem from higher levels of stress, weaker community support and higher rates of unemployment when living in a poorer neighborhood, all of which could also predispose the residents to a higher risk of developing psychiatric symptoms when exposed to community ambient air pollution (Pun et al 2017). In this study, we found neighborhood SES or housing conditions (such as area-level household income and house value) modified the associations between ozone and a global symptoms severity measure for psychiatric symptoms. This suggests that living in a low SES area may enhance the effect of ozone on the intensity of psychiatric symptoms. This further suggests that an intervention to lower ozone levels in low SES areas would be protective for psychiatric symptoms. The higher effect estimates for ozone in these neighborhoods may also be due to those communities having more open windows and less air filtration/conditioning, and hence more ozone penetration. Therefore, people living in poorer communities may also be more exposed to ozone compared to richer communities.
In the Derogatis 1983 US community sample looking at predominantly Caucasian nonpsychiatric-patient adults of both sexes with a mean age of 46 years (SD: 14.7) the population had an average GSI and PSDI raw score of 0.30 (SD: 0.31) and 1.29 (SD: 0.41) respectively (Derogatis and Melisaratos 1983). The average BSI scores for male populations are expected to be lower than for females. This was confirmed in both Israel and British community samples (Francis et al 1990, Gilbar andBen-Zur 2002). In the British community norms sample (mean age: 44 years), among male adults, the average GSI raw score was 0.36 (SD: 0.40) and the average PSDI raw score was 1.32 (SD: 0.51) compared with 0.50 (SD: 0.50) and 1.42 (SD: 0.58) in females (Francis et al 1990). In our analytical sample, we had an average GSI raw score of 0.24 (SD: 0.32) and an average PSDI raw score of 1.22 (SD: 0.34). This indicates that we studied a group of older men with, on average, less distress compared with general adults. However, harmful effects of ambient exposures on psychiatric symptoms were still observed among them.
We acknowledge that this study has some limitations: first, our sample is a subset of NAS participants, who are mostly older Caucasian men with mean BSI scores lower than the general adults in community samples. This may limit the generalizability of our findings to other race-ethnicity, age and sex populations as well as to other countries. Second, although the participants were followed up at each scheduled visit for routine physical exams and the collection of many covariates, the BSI questionnaire was mailed to the participants about a month before his scheduled visit and the participants were asked to fill in the survey and bring it to the visit. Thus, there is some uncertainty in identifying the outcome assessment date which affects the assignment of prior exposure windows. However, we believe the date difference would be relatively small for most of the assessments, since people most likely fill out the questionnaire shortly before the visit. Exposure misclassification will be larger for the average prior one week window, but less important for the average exposure in the previous four weeks or longer windows. This kind of exposure measurement error should be non-differential, so the significant findings with an average one week before should still be valid. Third, residual confounding is still possible. We did not have information on individuals' stress, sleep or trauma, which could be risk factors for psychiatric symptoms, particularly over the short-term. However, for those to be confounders they would need to be associated with air pollution. We were able to control for a relatively comprehensive list of covariates, including basic demographic, physical activity (metabolic levels), behavioral factors and nSES. In fact, research has shown that controlling for area-level factors will remove most of the confounding bias for ambient-environment-based studies since ambient exposures are related to area-level, rather than individual level factors (Sheppard et al 2012, Weisskopf andWebster 2017). We were cautious about controlling for additional physical health variables, such as CVD, respiratory diseases, stroke or HT status, etc, because they could be instead on the causal pathway. Our area-level ambient exposure level itself would not be affected by any individual behaviors, thus the chance for smoking, drinking and physical activity to behave as confounders is low. We presented results under both adjustment types in this paper, but would give preference to the results generated from model I. Lastly, the ambient exposures were matched to each individual based on his residential addresses history. Because address history was obtained at each visit, we were able to assign exposure at each visit, for the different averaging times, based on how long they had lived in each location. Consequently, we were able to match the nationwide predictions for air pollution and temperature to their addresses in other states as well as in the Boston area. While the exposure prediction models we fit were nationwide, using over 2000 monitoring locations for calibration, we did not see any evidence of more exposure error in New England, where the large majority of the participants lived. Exposure measurement error in the predictions themselves remains present, however, this type of error tends to be nondifferential and bias the effect estimates towards the null. Similarly, recall bias is also more likely to be nondifferential for the covariates and outcome information collection in NAS questionnaire administration. Selection bias is likely if we had non-random missing/selection from the inclusion of people with sufficient BSI response and without psychiatric medication use. To test whether this induced selection bias (which requires selection being associated with exposure), we fit logistic regressions predicting inclusion (being selected) in the study by each of the exposures. We found no evidence of selection being associated with the exposures, indicating no evidence of related selection bias due to the inclusion criteria (table S5).
In this longitudinal study looking at men in the VA NAS, we found evidence linking ambient gaseous pollutant exposure, such as ozone and NO 2 , to increased psychiatric symptom intensity. The associations were for short-medium term exposure windows. Area-level household income and house value modified the associations.

Data availability statements
Ms Qiu has full access to all of the de-identified analytical data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis. Data is under the regulation of Massachusetts Veterans Epidemiology Research and Information Center. Access to the data requires approval from the center and the US Department of Veterans Affairs first.
The data generated and/or analyzed during the current study are not publicly available for legal/ ethical reasons but are available from the corresponding author on reasonable request.

Fundings
This work was supported by National Institute of Environmental Health Sciences (NIEHS) Grants (P30 ES000002, R01ES015172) and US Environmental Protection Agency (EPA) Grant (RD83587201). It was also supported by AG002207 and AG018436 grants from the National Institute on Aging (NIA). Its contents are solely the responsibility of the grantee and do not necessarily represent the official views of the USEPA. Further, USEPA does not endorse the purchase of any commercial products or services mentioned in the publication. The funders had no role in the design and conduct of the study; collection, management, analysis, and interpretation of the data; preparation, review, or approval of the manuscript; and decision to submit the manuscript for publication. The funding sources support open access publishing.

Ethical statement
This study was approved by the Human Research Ethnics Committees of the Harvard School of Public Health and the VA Boston Healthcare System. The research was conducted in accordance with the principles embodied in the Declaration of Helsinki and in accordance with local statutory requirements. All the participants gave written informed consent to participate in the study.