Empirical evidence of mental health risks posed by climate change

Data.

We dichotomize our BRFSS measure of mental health difficulties, as recall of mental health difficulties may be susceptible to memory limitations given the interaction of psychological distress and cognitive functioning (56), adding error to the 30-d measure (48). Supporting this, the distribution of the answers to our mental health variable presents a clustering at weeks of the month as well as at 5-d intervals, although the true underlying distribution is likely continuous (SI Appendix, Fig. S9). Whether or not a person experienced any mental health difficulties over the past 30 d (our dichotomized measure) is likely more reliable. SI Appendix, Fig. S1 plots respondents over time.

We use gridded (at ∼4 km) daily temperature and precipitation from the PRISM Climate Group (57) and average daily cloud cover, relative humidity, and wind speed from the National Centers for Environmental Prediction Reanalysis II project (58). Global circulation model projections use NASA Earth Exchange’s maximum temperature and precipitation projections for 2010, 2050, and 2099 (59) from 21 of the Coupled Model Intercomparison Project Phase 5 (CMIP5) models (60) run on Representative Concentration Pathway 8.5 (RCP8.5) (high-emission scenario) (61). SI Appendix, Figs. S10 and S11 presents the spatial and temporal distributions of these projected changes.

Pooled Cross-Section.

Our relationship of interest is the marginal effect of recent meteorological conditions on the probability of experiencing mental health issues. We model this asYijkts=f(AVG.TMAXijkts)+g(PRECIP.DAYSijkts)+Zη+αj+μt+νks+ϵijkts.[1]In this pooled cross-sectional linear probability model, i indexes individuals, j indexes cities, k indexes states, t indexes calendar days, and s indexes calendar years. Our dependent variable Yijkts is binary and represents whether respondents in city j in state k on calendar day t within calendar year s reported experiencing any days of poor mental health over the 30 d before their interview date.

AVG.TMAXijkst and PRECIP.DAYSijkts represent the 30-d average of daily maximum temperatures and number of precipitation days over the same 30-d window as respondents’ reported mental health, respectively. We control for average diurnal temperature range, percentage of cloud cover, relative humidity, and wind speed, which are represented via Zη, as failure to do so may bias our primary estimates (45). Bold terms in equations denote matrices.

We use separate indicator variables for each 5 °C 30-d average maximum temperature bin and for each 5-d bin of precipitation days, which are represented by f() and g(), respectively. This allows for flexible estimation of a nonlinear relationship between our main meteorological variables and mental health outcomes (50). We omit the 10 °C to 15 °C maximum temperature and 0–5 precipitation days indicator variables and interpret our estimates as the change in probability of reporting recent mental health issues associated with a particular temperature or precipitation day range relative to these baselines.

To ensure that time-invariant city-specific factors do not bias our estimates of the effect of weather on mental health outcomes, we include αj—representing nonparametric city fixed effects—in Eq. 1 (48). Furthermore, there may be unobserved daily considerations or region-specific time trends influencing our mental health outcomes that could spuriously correlate with the weather. To control for these potential confounds, we include μt and νks in Eq. 1, representing calendar date and state-by-year fixed effects, respectively.

Our empirical identifying assumption is that the remaining variation in 30-d maximum temperature and precipitation is as good as random after conditioning on these fixed effects (45, 62) (SI Appendix, Fig. S2). The estimated coefficients from f(AVG.TMAXijkst) and g(PRECIP.DAYSijkts) from this estimation can thus be interpreted as the causal effects of maximum temperature and precipitation days on the 30-d occurrence of mental health issues (45, 62). We adjust for spatial and serial correlation in ϵijkts by using heteroskedasticity-robust SEs clustered at the state level (48). We omit nonclimatic control variables from Eq. 1 because of their potential to generate posttreatment bias in our parameters of interest (18, 45).

Our estimates are robust to varying functional form specification of our main meteorological variables. SI Appendix, Fig. S3 shows that results are robust to varying the size of our temperature and precipitation bins. SI Appendix, Fig. S4 and Table S2 shows that results are robust to modeling average maximum temperatures with a quadratic and precipitation days linearly. SI Appendix, Fig. S5 and Table S3 shows that results are robust to using the full distribution of daily maximum temperatures over the period (63). Coefficient estimates in this model represent the effect on probability of reporting mental health issues of having 1 additional day of the 30-d window falling in the specified maximum temperature bin.

Our estimates are also robust to varying control variable inclusion. SI Appendix, Table S1 shows that results are robust to progressively dropping meteorological controls. SI Appendix, Table S4 shows that our functional forms persist even when choosing less stringent sets of fixed effects (although we prefer our primary fixed effects specification, as it controls for the most confounding variation). SI Appendix, Table S5 shows that our results are robust to varying inclusion of demographic controls. We model our demographic controls with nonparametric bins and introduce them as fixed effects into the estimation. Of note, our sample size in these regressions varies, as not all respondents answered each question.

Furthermore, SI Appendix, Fig. S6 shows that our results are not driven by any one particular state’s respondents, SI Appendix, Fig. S7 indicates that our results are not entirely driven by any one census division, and SI Appendix, Table S6 shows that our results are robust to aggregating to the city-day unit of analysis and weighting regressions by the number of respondents on each city-day. Finally, one might be concerned that meteorological factors produce differential rates of response to the BRFSS survey in ways that limit the external validity of our findings. SI Appendix, Table S7 presents minimal evidence of the effect of meteorology on sample collection, indicating a low threat of sample selection bias driven by the meteorological variables that we use. Finally, our pooled cross-sectional results are robust to using the count of mental health days outcome variable in the context of both an ordinary least squares and a negative binomial model (SI Appendix, Table S19). In the negative binomial model, we substitute year-month, day-of-week, and city-level fixed effects for computational tractability (to our knowledge, no package exists in R to estimate the size of our full fixed effects matrix in a negative binomial model).

Long Differences.

To model the relationship between multiyear changes in average maximum temperatures and multiyear changes in average mental health outcomes, we take the mean of our dependent and independent variables over two separate periods, a and b, with period a falling temporally before period b (45, 47). Our main model uses data from period a (2002–2006) to period b (2007–2011) for each city. For example, for our mental health issues outcome measure, Yijkts, we calculate Yjka¯ as 1n∑t∈aYijkts and similarly, TMAXjka¯ for the TMAXjkts and Zη over the same period. This results in the below equation for period a:Yjka¯=βTMAXjka¯+Za¯η+αj+ϵjka¯.[2]In Eq. 2, j indexes cities, k indexes states, and αj represents city-level fixed effects that control for time-invariant city-specific factors. Calculating the same metrics for period b results in the below equation:Yjkb¯=βTMAXjkb¯+Zb¯η+αj+ϵjkb¯.[3]Then, differencing Eqs. 2 and 3 results inYjkb¯−Yjka¯=β(TMAXjkb¯−TMAXjka¯)+(Zb¯−Za¯)η+(αj−αj)+(ϵjkb¯−ϵjka¯).[4]In Eq. 4, the city-level, time-invariant fixed effects (αj) terms drop out (and are thus not potential confounds). Thus, the long-differences approach implicitly controls for time-invariant city factors in the way that we explicitly control for them in the pooled cross-sectional regressions. We are left with our long-differences model:ΔYjk=βΔTMAXjk+ΔZη+Δϵjk.[5]β from Eq. 5 provides an estimate of multiyear warming on multiyear change in rates of individuals reporting mental health issues. In Eq. 5, ΔZη represents changes in other meteorological variables. We weight the regression by the number of respondents in a city for statistical efficiency and again, cluster our errors Δϵjk at the state level. Unbiasedness of β requires that city-level changes in mean maximum temperature between our two periods are uncorrelated with unobserved time-varying variables that alter our mental health outcomes of interest (47).

SI Appendix, Table S9 shows that results are robust to excluding those cities with missing observations in either period a or period b. SI Appendix, Table S10 shows that results are robust to using percentage change in maximum temperature from baseline rather than absolute level of temperature change.

Furthermore, given that unobserved time-varying variables could confound our estimation, SI Appendix, Table S11 shows that our long-differences results persist across a varied set of demographic controls calculated using the BRFSS demographics questions. Adding to these results is SI Appendix, Table S12, where we use demographic variables drawn from the US Census Bureau to externally measure changes in demographic composition of cities. Our results are robust to the inclusion and varied specification of these controls as well. However, as these demographic variables may themselves be driven by longer-term changes in meteorological variables and thus, serve as “bad controls,” we prefer the specifications that exclude demographic controls (45). Finally, our long-differences results are also robust to using the count of mental health days outcome variable (SI Appendix, Table S19).

Difference-in-Differences.

We leverage the landfall of Hurricane Katrina to estimate the relationship between disaster exposure and mental health outcomes using a difference-in-differences design (48):Yijkts=γPOSTijkts+κTREATijkts+ψPOSTijkts×TREATijkts+ϵijkts.[6]In Eq. 6, Yijkts represents an individual’s mental health as in Eq. 1. We retain indices from earlier equations and again cluster our errors at the state level. POSTijkts is equal to one for the post-Katrina period in our data (from August 24, 2005 to February 1, 2012) and zero otherwise. TREATijkts is equal to one if city j received a federal disaster declaration in response to Hurricane Katrina and zero otherwise. POSTijkts×TREATijkts represents the interaction of these two terms. Our estimate of interest, ψ, represents the relationship between a city’s experience with Hurricane Katrina and residents’ subsequent mental health outcomes. Unbiasedness of ψ requires the assumption of no contemporaneous shocks and that disaster affected areas were trending similarly to nondisaster areas before hurricane landfall (48).

SI Appendix, Table S14 indicates that results are robust to using a Katrina exposure indicator variable coupled with city and date fixed effects. This specification more flexibly controls for time-invariant city-specific factors as well as any idiosyncratic temporal factors. Furthermore, SI Appendix, Table S15 shows that results are robust to the inclusion of a varied set of respondent-level demographic controls to attempt to control for possible bias induced by unobserved heterogeneity (48). However, as these demographic variables may be affected by Katrina exposure, we prefer the specifications that exclude demographic controls (45). SI Appendix, Table S16 displays that our results are robust to aggregating to the city-day and weighting the regression by number of respondents per city-day.

We also explore alternative selection rules for Katrina-exposed and non–Katrina-exposed groups. SI Appendix, Table S17 shows that results are robust to using only neighboring nondisaster declaration cities from states that geographically border the states that had cities with disaster declarations as the control group. Furthermore, not all cities in the sample are represented in both pre- and post-Katrina periods. SI Appendix, Table S18 shows that results are robust to selecting respondents from only those cities with representation in the sample both before and after landfall. Finally, our difference-in-differences results are robust to using the count of mental health days outcome variable (SI Appendix, Table S19).