Association Between Automotive Assembly Plant Closures and Opioid Overdose Mortality in the United States

Key Points Question Are closures of US automobile assembly plants associated with increases in opioid overdose mortality rates among working-age adults? Findings In this difference-in-differences study, US manufacturing counties that experienced an automotive assembly plant closure were compared with counties in which automotive plants remained open from 1999 to 2016. Automotive assembly plant closures were associated with a statistically significant increase in county-level opioid overdose mortality rates among adults aged 18 to 65 years. Meaning Automotive assembly plant closures were associated with increases in opioid overdose mortality, highlighting the potential importance of the role of declining economic opportunity in the US opioid overdose crisis.


Section 1. Identification of automotive assembly plants
Our strategy to identify automotive assembly plant closures followed the approach used in the academic literature and federal government reports, which involves triangulating data from industry trade publications, automotive company websites, and newspaper articles. [1][2][3] We first obtained an initial listing of automotive assembly plants in operation in North America as of 2005-2006 from Automotive News (AutoNews.com), 4 a leading industry trade publication. This was the earliest publicly-available and industry-verified census of automotive assembly plants during our study period. We considered all plants that were involved in manufacturing of automobiles, trucks, and heavy trucks. To identify plants that were in operation as of 1999, but had closed by 2005-2006, we conducted a series of Google and Google News searches using the each of the following search terms: [ The field [Company name] refers to each of the following automotive manufacturers: AM General, Autoalliance, BMW, FCA US LLC (or "DaimlerChrysler" or "Chrysler"), Ford, Freightliner, General Motors (or "GM" or "G.M"), Honda, Hyundai, International, Kenworth, Mack, Mitsubishi, Nissan, Nummi, Peterbilt, Sterling, Subaru, Toyota, Volvo, and Western Star. The field [Year] refers to individual years between 1999 and 2005. Using this strategy, we found six additional automotive assembly plants that were operational as of 1999 but closed before 2005.
To identify automobile assembly plant closures from 2005 onwards, we used a combination of company websites (e.g. https://corporate.ford.com/company/operation-list.html#s1f0) and Google news searches using either specific plant names or company names with plant locations. Specifically, we conducted searches using the following terms: The list of automotive assembly plants identified through our search algorithm, along with their location, and, if relevant, year of closure are provided eTable 1. The list is divided by plants that were represented in the main study sample versus those that were not, based on the sample restrictions discussed in Section 2 of this Supplement.

Section 2. Sample and exposure
To define the study sample, we first restricted the study sample to all counties within commuting zones that contained at least one operational automotive assembly plant as of 1999. Commuting zones represent groups of counties that capture predominant residential and commuting patterns. The use of commuting zones to define exposure therefore accounts for the possibility that individuals likely to be affected by automotive assembly plant closures may live in a county other than the one in which the plant was located. Commuting zones have been widely used to define local labor markets. [11][12][13] Because our study period ranged from 1999-2016, we used commuting zone definitions from the year 2000. Counties were considered exposed to an automotive assembly plant closure if they were located in a commuting zone that experienced a plant closure during the study period. In the 4 commuting zones where more than one automotive plant closure occurred, we assigned exposure based on the year of the first closure.
We defined counties that were most likely to be affected by automobile assembly plant closures using data on the share of county residents employed in manufacturing at the beginning of the study period. Specifically, we restricted our sample to counties in the top quintile nationwide with respect to the share of employed residents aged 16 years and above who were working in the manufacturing sector, based on data from the 2000 Decennial Census 14 ("manufacturing counties," eFigure 1). The use of the top quintile versus bottom four quintiles of the share of the workforce employed in manufacturing to define manufacturing counties follows from prior work, which (1) demonstrated that top quintile manufacturing counties experienced distinct socioeconomic trends during the study period, 15 and (2) defined the U.S. manufacturing sector as consisting of ~600 counties (i.e., approximately 1/5 th of all U.S. counties). 16 Imposing this sample restriction meant that 18 of the 48 commuting zones with at least one automotive assembly plant in operation as of 1999 were excluded, as they did not contain any counties with a share of manufacturing workers in the highest national quintile. As we discuss in Section 5C, our substantive findings were unchanged if we employed alternate definitions of exposure that included counties from all 48 of these commuting zones.
Our approach to identify sample counties was motivated in part by limitations in vital statistics and employment record data. Vital statistics data do not include information on previous occupation or industry of employment, and data on relevant proxies such as education are often missing or measured with error. 17 Employment databases, such as the U.S. Census Bureau's Quarterly Workforce Indicators, 18 allow for an exact enumeration of individuals employed in the automotive industry. However, the assignment of counties in these data are based on place of employment and not place of residence, which creates similar potential for measurement error as in the vital statistics data (which identify only the county of residence or death).
Our approach was also motivated by the fact that the consequences of an automotive assembly plant closure likely extend beyond just those workers employed in the automotive industry. These spillovers are likely present because individuals employed in other manufacturing industries have long sought employment in the historically higher-paying automotive sector. 19 Moreover, contraction of opportunities within the automotive sector likely affected other workers in other manufacturing industries through reduced demand for automotive parts from component suppliers 1,20 and potentially through changing expectations about the potential for socioeconomic mobility through manufacturing work.
Our final study sample consisted of 112 manufacturing counties situated in 30 commuting zones, including 29 exposed counties located in 10 commuting zones in which an automotive assembly plant closed during the study period, and 83 unexposed counties located in 20 commuting zones in which all automotive assembly plants remained open during the study period.

Section 3. Outcome variables
The primary outcome was the annual age-adjusted opioid overdose rate for 18-65-year-old adults in the county. We calculated these rates using restricted-access, individual-level multiple cause of death vital statistics data from the U.S. National Center for Health Statistics (NCHS) from January 1, 1999 to December 31, 2016. 5 We followed the strategy recommended by the U.S. Centers for Disease Control and Prevention by first using ICD-10 underlying cause codes X40-X44, X60-X64, X85, Y10-Y14 to identify drug overdose deaths and then using contributing cause codes T40.0-T40.4 to identify deaths specific to opioid overdoses. 6,7 Doing so, we identified a total of 10,759 opioid overdose deaths of 18-65 year-old adults in the 112 sample counties during the study period. We then used these data, along with population counts from the U.S. Census Bureau, 8 to compute age-adjusted rates per 100,000 individuals (based on a direct standardization method). Counties were assigned based on reported county of residence for each decedent in the NCHS data. We used a similar procedure to construct our secondary outcomes: rates of overall drug overdose mortality (ICD-10 codes X40-X44, X60-X64, X85, Y10-Y14), prescription opioid overdose mortality (subset of drug overdose deaths assigned contributing cause codes, T40.2, T40.3, and T40.4), and illicit opioid overdose mortality (subset of drug overdose death contributing cause codes, T40.0 and T40.1). 7 We examined overall drug overdose mortality as a secondary outcome in order to address potential bias from underreporting of opioids as a specific cause of drug overdose mortality during the study period. 9 This potential bias includes differential changes in the likelihood of identifying opioids as the causal agent specifically related to plant closures (e.g., due to rising overall drug overdose mortality rates raising concerns among medical examiners). The possibility of such bias is suggested by geographic variation in the identification of different contributing causes of drug overdose deaths. 10 The concordant results for both opioid and overall drug overdose mortality ( Figure 2) suggests that any bias from differential reporting was unlikely to have been substantial.
We used the same procedures to construct subgroup-specific opioid overdose mortality rates reflecting each combination of age (18-34 versus 35-65), sex (women versus men), and race/ethnicity (non-Hispanic white versus all other).

Section 4. Statistical model
For each primary and secondary outcome, we estimated the following least squares multivariable regression model: where i indexes the county, j the commuting zone within which the county is located, and t the calendar year. The subscript p refers to event years, which are annual intervals relative to the calendar year of automotive assembly plant closure. For example, period -3 refers to 3 years before plant closure. Event periods 5 or more years prior to plant closure were combined into a single bin, as were event periods 7 or more years after plant closure. This aggregation follows prior literature 21 as a means of avoiding difficulties in interpreting coefficients due to sample size imbalances created by differences in the timing of plant closures. The variable Closurejp denotes a series of binary indicators = 1 if the closure of an automotive assembly plant in operation since 1999 had occurred in commuting zone j during or before the calendar year associated with eventtime period p. Yijt denotes the outcome of interest. County fixed effects (denoted by i) captured time-invariant differences in socioeconomic and cultural characteristics across sample counties. Nationwide secular trends in the outcomes, for example, due to the Great Recession or differential changes in availability of prescription or illicit opioids were captured by calendar year fixed effects (denoted by t).
This model estimates the difference in outcomes for leads and lags of plant closure relative to a reference year (event time -1) and relative to all manufacturing counties in commuting zones where a plant had not closed during the study period (i.e., for whom Closurejp is equal to zero for all event periods). These relative differences are captured by the coefficients p, estimates for which are plotted in all of the main and supplement figures.
By allowing associations between exposure and the outcome to vary over time, our specification represents a generalization of the method of difference-in-differences -known as the "event study specification" in modern econometrics. [22][23][24] Specifically, in cases where the timing of exposure varies across units and where effects vary with time, event study models are preferred to standard difference-in-differences models because the single coefficient reported in difference-in-differences analyses represents a weighted average of all combinations of comparisons between sample units. 22 However, in some of these comparisons, units exposed earlier may serve as controls for units exposed later. If differences in outcomes are increasing with event time (as is the case for the present study), the estimates from these particular comparisons will be attenuated toward the null. Hence, the single difference-in-differences estimate may provide a misleading estimate of the true association between exposure and outcome. The event study specifications, which also provide a more transparent test of violations of the parallel trends assumption required for causal inference, 26 avoid this problem by indexing the reference point to time since the event (as opposed to calendar time) and by allowing associations to vary over time.
For all models, we computed 95% confidence intervals adjusting for serial correlation in the outcome at the commuting zone level. 27 We weighted all regressions by county population size at baseline (1999).

A. Modelling overdose deaths as a count variable
We estimated our main regression model (Equation 1) using least squares given this estimator is less prone to finite sample bias from the inclusion of numerous fixed effects (a well-known problem that arises with maximum likelihood estimators such Poisson or negative binomial regression 28 ). We nevertheless assessed the sensitivity of our findings to this choice by modelling the number of opioid overdose deaths, using a generalized linear model with a negative binomial distribution and a log-link (and with county-year population specified as the exposure). 29 The resulting coefficients can be interpreted as relative changes in mortality rates.
The substantive findings were robust to using this alternate estimate strategy (eFigure 6). In particular, the estimates suggest that at 5 years after a plant closure, opioid overdose mortality rates increased by 85% (95% CI: 41%, 128%, p<0.001) in exposed counties relative to unexposed counties. (The corresponding estimate from our main specification implies a 75% increase relative to unexposed counties by the same time point.)

B. Wild-cluster bootstrap-t procedure for variance estimation
The Huber-White method to account for geographic clustering may result in standard errors that are biased downwards in cases where there are few clusters. 30 What constitutes a small number of clusters depends on the application. Our main study sample specified 30 clusters (commuting zones), which is generally considered a large enough number of clusters according to conventional rules of thumb. Nevertheless, we tested the robustness of our statistical inference using the wild cluster bootstrap method. 31 To implement this procedure, we used the Stata command -bootwildct-developed by Bansi Malde and Molly Scott (available at https://www.ifs.org.uk/publications/6231). As shown in eTable 2, the p-values on the event study coefficients were substantively unchanged when using this alternate procedure.

C. Including counties from commuting zones excluded from main sample
As discussed in Section 2, we restricted our sample to counties in the top quintile nationwide with respect to the share of employed residents working in the manufacturing sector, in order to identify areas of the country most likely to be economically affected by automotive assembly plant closures. However, doing so excluded counties in 18 of the 48 commuting zones with at least one operational automotive assembly plant as of 1999. To assess whether our results were sensitive to the inclusion of these commuting zones, we added to our main sample the county in each excluded commuting zone with the highest share of residents working in the manufacturing industry. Results for this alternate sample, which included 130 counties (37 exposed and 93 unexposed), were similar to those for the main study sample (eFigure 7).

D. Alternate control group
Our primary regression models defined unexposed manufacturing counties as those located in commuting zones without an automotive assembly plant closure. We hypothesized that these unexposed counties were likely to follow similar pre-existing trends in the outcomes compared with exposed counties, given that they would also be subject to similar automotive industryspecific social, economic, and cultural factors that may be relevant to population health. However, plant closures in one commuting zone may cause automotive workers in commuting zones where plants have remained open to experience economic uncertainty, which may be consequential for health. 32 Another possibility is that automotive manufacturers preferentially closed plants in areas where their workers were experiencing downward trends in health. However, both the known rationales for specific plant closure decisions 19 and the absence of differential pre-existing trends in the outcome between exposed and unexposed counties ( Figure  2) suggests that bias from selective plant closure is unlikely.
These potential spillovers could violate the Stable Unit Treatment Value Assumption (SUTVA). 26 To assess the robustness of our findings to this potential violation, we repeated the analysis with an alternate control group: manufacturing counties located in non-automotive commuting zones that were situated in the same states as commuting zones with automotive plant closures. The restriction of counties to the same states helps to ensure similar exposure to broad, regional trends in social or economic factors that could influence health outcomes.
However, this control group should be less prone to cross-commuting zone spillovers (though the problem is not fully eliminated given that firms supplying automobile component parts, who would face the downstream consequences of plant closures, tend to locate in similar regions as automotive plants 1 ). The estimates using this alternate control group were substantively similar to our main findings (eFigure 8), suggesting that cross-commuting zone spillovers are unlikely to bias our findings.

E. Repeating the analysis under conditions expected to produce a null result
We did not expect to find strong associations between automotive assembly plant closures and opioid overdose mortality in non-manufacturing counties (i.e., counties in the bottom four quintiles of the county workforce employed in manufacturing). Consistent with this expectation, the estimated associations between plant closure and opioid overdose mortality in nonmanufacturing counties were smaller in magnitude and not statistically significant (eFigure 9).

F. Selective migration
The association between automotive assembly plant closures and rising opioid overdose mortality rates may reflect, in part, selective outmigration of individuals at lower risk of developing substance use disorders. While recent work in labor economics suggests that differential outmigration after adverse economic events may be minimal, 12,33 we nevertheless examined whether manufacturing counties exposed to plant closures experienced greater (net) outmigration relative to unexposed manufacturing counties.
We obtained intercensal, county-level data on net migration rates from the University of Wisconsin-Madison, Net Migration Patterns for US Counties database. 34 We used these data to create net migration rates by each combination of age group (18-34 years versus 35-64 years) and sex. We investigated these groups separately given potential differences in the propensity to migrate, along with differential exposure to automotive assembly plant closures (Figure 3). Because these data are only available at intercensal (i.e., decadal) intervals, we could not estimate the event study specification in Equation 1. Instead, we estimated the following least squares regression model, for each age-sex group: (2) ∆ if the manufacturing county was located in a commuting zone where an automotive assembly plant in operation at the start of the study period closed during the specified time frame. We split closures to take into account any potential lagged effects of economic decline on migration. 35 Equation 2 assesses whether net migration rates changed more in manufacturing counties located in a commuting zone where an automotive assembly plant closed versus manufacturing counties not exposed to plant closure . Because the dependent variable is the change in decadal net migration rates, and the independent variables effectively capture a commuting zone change in the operational status of one or more automotive assembly plants, Equation 2 effectively represents a first-differences model. Like our event study specifications, which include countyspecific fixed effects, first-differences models adjust for potential confounding by time-invariant, county-specific variables, whether observed or unobserved. 36 This specification is the closest analog to the model in Equation 1 that was possible to fit given the constraints of the migration data.
Consistent with the recent labor economics literature, 12,33 we found no evidence that automotive assembly plant closures were associated with net migration rates for any era of plant closure (eFigure 10). Across specifications, even the estimates with the largest magnitudes represent only a 0.13 standard deviation decrease in net in-migration (not outmigration) rates.
There are two caveats to these findings. First, the confidence intervals are relatively wide, so we cannot exclude the possibility that there were substantively meaningful differences in net outmigration rates attributable to plant closures. Second, the findings only apply to a time frame of 8 years post plant closure (because the last decadal time point to assess net migration was 2010 and the first closure in our sample was 2002). Thus, we cannot rule out differences in out-migration over a longer time frame. That said, the null findings here are useful to support the veracity of our event study estimates through at least 6 years of follow-up after plant closure. Notes: Coefficient estimates from the main event study model for the primary outcome of opioid overdose mortality rates (per 100,000 individuals aged 18-65 years). All models include county and calendar year fixed effects (see main text and Section 4 of this Supplemental Appendix for further details on the event study specification). Point estimates for each event time binary variable are bolded. 95% CI, adjusting for clustering at the commuting zone level, are in square brackets. Below these are p-values based on clustercorrected Huber-White standard errors. Wild cluster bootstrap p-values are denoted by "wild p." Section 5B of this Supplement provides further details on statistical inference using the wild-cluster bootstrap-t method.