Assessment of the Legionnaires’ disease outbreak in Flint, Michigan

Deidentified data on LD cases from 2010 to 2016 were obtained from the MDHHS by Data Use and Confidentiality Agreements following approval from the Institutional Review Board for the Protection of Human Research Subjects (MDHHS IRB 201608-01-EA, Wayne State University IRB 067016B3E). Data represent a complete inventory of LD cases in Genesee, Oakland, and Wayne Counties over this 6-y period. Each case is time-stamped with dates of symptom onset, patient diagnosis, and referral to the Michigan Disease Surveillance System (MDSS). The precise date of referral to the MDSS is available for all LD cases. Only 623 of the 833 (25.2% missing) LD cases in Genesee, Oakland, and Wayne Counties have a recorded diagnosis date. Although MDHHS staff generously provided enhanced onset date data collected from case investigations, only 694 cases (16.7% missing) had a verified date of symptom onset. In consultation with scientific personnel at MDHHS, the date of the referral is the most reliable and valid indication of timing. Median difference in elapsed time between referral and diagnosis dates is 1 d, and 5 d between referral and symptom onset dates. These lags inform our use of referral timing data in chlorine analyses. Given the very high correlation between referral and onset date (R² = 0.9973), using referral date with a time lag adjustment in chlorine models resolves the timing error and preserves maximum information (limiting the missing information bias that arises with use of onset date). MDHHS data are also referenced geographically by the residence of the LD case at the census tract scale. Of the 833 LD cases observed over this time period, all but 27 cases included a verifiable address (or census tract residential indicator), including 3 in Genesee, 2 in Oakland, and 22 in Wayne County.

To test the water regime hypotheses, we exploited the temporal and spatial properties of MDHHS case data to develop an outcome variable, LD incidence, that is observable in time (before and after the switch in water regime) and space (in and outside regime and chlorine-treated neighborhoods). LD incidence is a binary variable equal to 1 if a confirmed case of LD is observed in census tract $i$ in week $t$ and 0 otherwise.

To estimate LD effects from the switch in water source and treatment, and the ensuing variability in chlorine residuals, we also collected a suite of demographic and meteorological control variables. With respect to neighborhood (census tract) demography, two variables from the US Census Bureau are used: percent of population receiving public assistance and percent of population $\ge$50 y of age. Both socioeconomic status and age ($\ge$50 y) are known correlates of LD risk. Three county-level meteorological variables are used: average weekly temperature, average weekly humidity, and average weekly precipitation. The thermal forces of temperature, humidity, and precipitation are known to govern the observed seasonality of LD incidence rates through growth effects on legionellae bacteria (25–27).

To capture effects of changes to source water and treatment, or water regime, we deployed a quasi-experimental method called difference-in-differences. An illustration of the method is summarized in Table S2. Our first difference is spatial, corresponding to whether a census tract is located inside Flint (F), and therefore treated by the shift in water regime, or not in Flint (NF). Our second difference is period-based, corresponding to whether a census tract is observed before (A) or after (B) the switch in water regime. The difference-in-differences method infers a causal effect by comparing the difference between A and B on the outcome of interest (LD incidence) for treated census tracts (F) relative to not treated census tracts (NF).

A key assumption of the difference-in-differences method, known as the parallel paths requirement, posits that the average period difference (A − B) in control group (NF) census tracts constitutes the counterfactual average difference between A and B in F census tracts if not for the treatment or switch in water regime. In this analysis, the preperiod parallel paths requirement is satisfied (Supporting Information). The differential behavior of LD incidence rates in post-water regime switch period of Fig. S1 is indicative of a powerful place-specific period effect. Note in Fig. S1 the extraordinary increase in the LD incidence rate in Genesee County. Although it is analytically tempting to conclude that the switch in water regime governs the LD spike in Genesee, our analysis plan aimed to rule out forces coincidental with the regime switch, evaluate plausible alternative explanations, and identify a potential causal mechanism for the observed increase in LD in Genesee County.

Our analysis begins by identifying regime effects. In the regime effect equations detailed below, our first difference is always geographic, corresponding to whether a census tract is located in Flint (and is therefore a recipient of Flint water) or is not in Flint (either located in Oakland, Wayne, or parts of Genesee County not in Flint). Our second difference is a period indicator that is variously defined by whether or not parameters from a given census tract are observed before the first switch in water source (phase A), after the switch to the Flint River but before the issuance of water quality alerts (phase B), after the switch and after the issuance of water quality alerts (phase C), and after the switch back to the DSWD water supply (phase D). Fig. 2A summarizes the precise timing of each phase. Expectations of parameter behavior involving various statistical comparisons by phase are described below.

We begin with a more global test of a water regime effect, estimating a baseline census tract random effects logistic equation for the probability of census $i$ in week $t$ presenting with a case of LD (1 = yes, and 0 = no):

where $\mathrm{\Lambda }[\cdot ]$ is the cumulative distribution function (CDF) of the logistic distribution; $R_{it}$ is a vector of temperature, precipitation, and humidity measures (from Weather Underground); $X_{it}$ is a vector of census tract control variables including percentage of population $\ge$50 y of age and percent of households receiving public assistance (from US Census Bureau); ${\zeta }_i$ is the random effect of census tract $i$; $F_i$ is an indicator variable = 1 if the census tract is in Flint; and $P_t$ = 1 if the census tract is observed in the postswitch period (combining phases B and C detailed in Fig. 2A), with the treatment effect of the regime switch captured by the estimated coefficient $\left(\delta \right)$, constituting our difference-in-differences of $F$ and $P$. In the presentation of logit model results below we exponentiate the estimated coefficient $\delta$ to derive an odds ratio, with the expectation that ${\text{exp}}^{\delta }>1$ indicating that the switch from Detroit to Flint River water caused an increase the risk of a census tract in Flint presenting with an LD case.

Next, we test whether boil water advisories issued by authorities in Flint attenuated the risk of LD, providing an additional probe of whether the switch in water regime caused the observed spike in LD incidence in Flint. We divide the postswitch period into two phases, B and C (as detailed in Fig. 2A); compare LD outcomes in census tracts in phase B (after switch, before advisories) and phase C (after switch, after advisories) to phase A (before the switch in water regime); and then compare outcomes in phase C to B. These comparisons assume that the issuance of boil water advisories induced a meaningful reduction in water use by residents. Boil advisories were not issued to address the presence of legionellae in the Flint water supply: authorities issued advisories after positive tests for E. coli contamination. However, the advisories may have confirmed suspicion among residents that the drinking water was unsafe, thereby unintentionally increasing water avoidance by residents. Two sources indicate that the advisories meaningfully affected residential water exposure. First, Google Trend search interest data on water contamination in the Flint–Saginaw–Bay City metropolitan area increased measurably around boil advisory dates, indicating awareness of the official warnings among the local population (see ref. 19). Second, a large, statistically significant, and sustained increase in sales of bottled water in Genesee County corresponded with the issuance of advisories (18).

To examine whether advisories reduced LD risk (through a water avoidance pathway), we estimate a census tract random effects logistic equation for the probability of census $i$ in week $t$ presenting with an incidence of LD (1 = yes, and 0 = no):

where all terms carry from Eq. 1, with the exception of $P{B}_t$, which is equal to 1 if the census tract is observed in the postperiod but before boil water advisories, and $P{C}_t$, which assumes a value of 1 if the census tract is observed in the postperiod and after boil water advisories. The comparison of estimated coefficients ${\delta }_1$ and ${\delta }_2$ indicates whether water avoidance behavior of residents in Flint helped attenuate the LD outbreak and provides support for the water regime hypothesis. Insofar as waterborne exposure to legionellae in Flint is linked to LD risk in a given census tract in time, and boil water advisories helped to reduce LD risk in Flint by inducing water avoidance in resident population, it is expected that ${\text{exp}}^{{\delta }_1}>1$, ${\text{exp}}^{{\delta }_2}>1$, and ${\text{exp}}^{{\delta }_1}>{\text{exp}}^{{\delta }_2}$.

We analyze whether the switch back in water supply on October 16, 2015, caused a reduction in LD incidence by estimating a census tract random effects logistic equation for the probability of census $i$ in week $t$ presenting with an incidence of LD (1 = yes, and 0 = no):

where all terms carry from Eq. 1, with the exception of $P{D}_t$ which is equal to 1 if the census tract is observed in the switch back period (phase D in Fig. 2A) and 0 if observed in the preswitch period (phase A in Fig. 2A). The causal effect of the switch back in water regime is captured by the estimated coefficient $\left(\delta \right)$. Insofar as the rise and fall of LD incidence in Flint was caused by a city-wide failure in water treatment, this test is expected to yield an ${\text{exp}}^{\delta }\approx 1$, indicating that the switch back to Detroit water returned the risk of a census tract in Flint presenting with an LD case to precrisis levels.

Free chlorine (mg/L as Cl₂) measurements at eight monitoring locations in Flint from 2013 to 2016 were obtained from the Monthly Operating Reports provided by the Flint Water Department (Supporting Information). Chlorine was measured at each location two to three times per week. Fig. 1B illustrates the behavior of average free chlorine at the weekly time step at each monitor site. Vertical lines bisecting the space correspond to water regime switch moments. Note the high between-monitor agreement in the level of free chlorine in the preswitch period, indicating negligible spatial variation in water quality across the City of Flint. In the postswitch period, we observe extraordinary temporal (or within monitor) and spatial (or between monitor) variation in the level of free chlorine.

Table S1 summarizes the statistical behavior of free chlorine within and between monitors in time illustrated in Fig. 1B. Analytically, the unprecedented exogenous variation in free chlorine levels observed in Fig. 1B and Table S1 is what we exploit to identify statistically the effect of changes in water quality on LD incidence in Flint, MI.

To test the chlorine residual hypothesis, we needed to select a method for associating a location in the City of Flint map with chlorine residual monitoring points. Although a number of approaches could be used, we choose to develop a physically relevant monitor-to-parcel assignment algorithm that leverages best available information on parcel occupancy/vacancy, residence time of water, and the Flint water distribution system pipe network. We tried a number of common approaches, including various proximity and Thiessen polygon-based methods that delivered quantitatively similar results. A hydraulic-based chemical transport model of the water distribution system could be used to estimate free chlorine residual. However, it is unclear if this will have a meaningful effect on the results given the spatial and epidemiologic limitations of surveillance data. The algorithm begins by finding the shortest path (or spine) from the centroid of each parcel to the Flint Water Plant (FWP) via the pipe network obeying the water age gradient. The water age gradient used is correlated with spatial variation in blood lead levels during the switch in water supply (28), demonstrating the utility of this metric to account for physical variability within the water distribution system. Water age is dynamic and likely varied spatially during the study period. Utilizing the water age gradient helps to incorporate major hydraulic constraints that proximity-based methods fail to accommodate (e.g., flow restrictions due to pipe size). This results in 41,286 parcel to FWP spines. Next, the algorithm finds the shortest path of each monitor to each spine, again obeying the water age gradient. Each parcel is then assigned the (weekly average) chlorine value of the monitor with a spine juncture nearest to the parcel. Because LD incidence data are organized at the census tract scale, we average parcel chlorine to the census tract to generate 8,000 fully observed census tract ($i$) by week ($t$) observations in Flint from 2013 to 2016. The outcome of the monitor-to-parcel assignment algorithm is provided in Fig. S2.

At adequate levels (commonly assumed to be concentrations $\ge$0.2 mg/L as Cl₂), chlorine effectively suppresses legionellae growth. Under normal circumstances, it is near impossible to identify statistically a chlorine residual → legionellae → LD incidence pathway because of insufficient time and space variation in chlorine within water distribution systems. As observed in Fig. 1B, the switch in water regime in Flint induced striking temporal and spatial variation of free chlorine in the Flint water distribution system. Our estimation strategy analytically leverages this quasi-random behavior in free chlorine throughout the city. We estimate a census tract random effects logistic equation for the probability of a census $i$ in week $t$ presenting with an incidence of LD (1 = yes, and 0 = no):

where all terms carry from Eqs. 1 and 3 with the exception $C_{it-1}$ denoting the average weekly free chlorine (mg Cl₂/L) at census tract $i$ in time $t-1$. The 1-wk lag in free chlorine is included to account for the difference between symptom onset and referral date information in MDHSS case data. Recall the use of referral date information was necessitated due to missing and imprecise data for symptom onset date. In addition to a continuous measure of free chlorine, we examine threshold effects of free chlorine, with $C_{it-1}$ = 1 if < 0.5 mg/L or <0.2 mg/L. Insofar as the loss of free chlorine increases the risk of LD, our expectation is that ${\text{exp}}^{{\beta }_1}<1$ in the continuous model and ${\text{exp}}^{{\beta }_1}>1$ in threshold models.

The census tract-specific residual, ${\zeta }_i$, in the random effects model is meant to capture the combined effect of all omitted census tract-specific covariates that cause neighborhood variation in LD susceptibility. Omitted variables may include the underlying health frailty of residents or other sources of water chemistry parameters that affect legionellae growth, such as iron, pH, water temperature, or assimilable organic carbon. The census tract-specific random effect measures the difference in LD risk in a given census tract versus LD risk across all tracts in the City of Flint. Results from a Hausman specification test (${\chi }^2=2.54,p=0.96$) indicate that model coefficients are efficiently estimated by random as opposed to census tract fixed effects.

To test the plausibility of the hospital-based outbreak hypothesis, we recapitulate Eq. 1 through Eq. 4 but limit analysis to non-hospital A-related LD cases. Potential hospital A-related cases are identified by screening (i) all cases in the MDHSS case file indicating admission to hospital A and LD between 2014 and 2016 and/or whether MDHHS staff, on the basis of case analysis, override the reported admission indication and assign the case as a hospital A admission (n = 64) and (ii) all non-hospital A hospitalized cases with epidemiological investigation notes indicating a prior admission to hospital A between 2014 and 2016 (n = 19). All 19 cases in screen ii appear in screen i, giving a total of 64 cases potentially related to hospital A. Granting the hospital outbreak hypothesis maximum explanatory power, we assume that all 19 cases with a prior admission to hospital A contracted LD at hospital A. Of the 45 cases remaining from screen i, all but 6 were returned to our analysis pool because the recorded symptom onset date was before the hospital admission date. Adding 6 and 19 gives 25 cases plausibly resulting from a hospital-based outbreak.

By limiting the analysis to non-hospital A-related cases, we test whether estimated water regime and free chlorine suppression effects appreciably change. If hospital A-related cases fully govern LD outcomes in Flint, then our statistical results pertaining to water regime and loss of free chlorine effects ought to disappear. However, if water regime and chlorine coefficients do not appreciably change with the exclusion of hospital A-related LD cases, then it is highly unlikely that the LD spike in Flint resulted from a hospital-based outbreak only. Although the hospital exposure thesis is not incompatible with our water regime/chlorine hypothesis—the hospital is similarly drawing water from a portion of Flint’s water distribution system where chlorine residual was often very low—our case exclusion tests allow one to rule out an exclusively hospital-based argument.

To test the hypothesis that the sizeable increase in LD incidence in neighborhoods (or census tracts) adjacent to Flint, in Genesee County, were due to water exposure in Flint, we utilized the Longitudinal Employer-Household Dynamics Employment Statistics dataset (https://lehd.ces.census.gov/data/). This dataset estimates that 15,857 workers flow into Flint every day, with 12,843 of them residing in neighboring areas in Genesee County, constituting a remarkable 61.7% of all employed persons working inside Flint. Although commuting data capture inflows for the purposes of work, it is reasonable to assume exposure to the Flint water distribution system through leisure and other activities as well.

Restricting to not-Flint Genesee County census tracts, we test the commuter flow hypothesis by estimating the following census tract random effects logistic equation for the probability of census tract $i$ in week $t$ presenting with a case of LD (1 = yes, and 0 = no):

where all terms carry from Eq. 1, with the exception of $GC$ which is equal to the observed count of daily commuters to Flint originating in non-Flint Genesee County census tract $i$. Insofar as exposure to Flint water is the source of the non-Flint outbreak, it is expected that $\delta$ increases monotonically in $GC$.