Association of State Medicaid Expansion Status With Rates of Suicide Among US Adults

Key Points Question What is the association of Medicaid expansion under the Patient Protection and Affordable Care Act with suicide rates? Findings Although suicide rates increased from 2000 to 2018 across all states in this cross-sectional study, a statistically significant attenuation of −0.40 suicides per 100 000 people was found in Medicaid expansion states compared with nonexpansion states, translating to 1818 fewer suicides from 2015 to 2018. Meaning These findings suggest that Medicaid expansion, which increases mental health coverage, is associated with a decrease in suicide mortality.


Data Format
We utilized mortality rate data, which captures the number of suicides within a given population at a given point in time. For our analyses, suicide mortality data were obtained at the state level for all states and D.C. for each year from 2000 -2018 and for each age group (20-29, 30-44, 45-64). This results in 19 years * 51 states (including D.C.) * 3 age groups = 2907 "observations," each of which includes the number of suicides, the number of the population, and the age-adjusted suicide rate and associated standard error for state i, year t, and age group j.
The population included in the denominator of the rate (# suicides / # population) represents the entire U.S. population ages 20-64 years, and the suicide mortality data encompass all reported suicides (based on ICD codes on death certificates) within this population. Due to the nature of this data format, conceptualizing the concept of "sample size" is difficult. One could view all individuals ages 20-64 years living in the U.S. from 2000 -2018 as the sample given their inclusion in the suicide rate calculations. Alternatively, since data are reported at the state level, one could also view the sample size as 51, one for each state, with repeated measures over time.
We prefer to report sample size in terms of the number of unique observational units (at the state-age-year level), which in our case is 2907.

Covariates
As above, data were obtained by state, year, and age group. These factors were all included in our model as covariates. Year was treated as a linear continuous variable (rescaled such that the year 2014 was 0). Due to concerns of heterogeneous trends in suicide rates over time by state, we also included a state*year interaction effects. Further information regarding the selection of the state-level trends is presented in the statistical analysis section.
We obtained time-varying (including both pre-and post-expansion values) state-level covariates, which were selected a priori due to suspected associations with suicide mortality and/or access to care. Census estimates from the United States Census Bureau were used to define race/ethnicity (% non-Hispanic white) data in 2010 (pre-expansion) and 2020 (postexpansion). 1 4 State-level data was then merged with the mortality data by state and year (or nearest matching year). The number of state opioid prescribing laws and firearm laws were obtained from recent publications. [5][6][7]
It can be shown that the DID estimate is equivalent to the interaction effect between state Medicaid expansion status (expansion vs. non-expansion) and time period (pre-vs. postexpansion) in a linear regression model, which also enables adjustment for other factors and modeling of correlation between observations. 8 We use a flexible hierarchical Bayesian regression framework for our linear regression approach to difference-in-differences to overcome some of the challenges associated with these data. Our modeling approach is largely based on a previously reported method. 9 For the sake of completeness of our methodology, we include that information here: With observed vector of mortality rates Yi and vector of associated estimated variances si 2 for state i, age group j, and year group t, diagonal matrix Si containing the standard deviations si, correlation matrix R, mortality rate average for the national population of interest m, and covariates k=1, …, K, we developed a hierarchical model such that: k~ Cauchy(0, 2.5) T(0, Infinity)  ~ Uniform(l,u) Suicide mortality rates were assumed to follow a normal distribution by the central limit theorem, as the rate can be viewed as a mean. 10 However, the variance for the mortality rates were not be equal across state-age-year units given the differences in population and number of observed events. A prior was not placed for these values given that the estimated variance is simply a function of other values. We utilize a multivariate normal distribution to obtain all draws for a given state (across all years t and age groups j) to be able to explicitly model the correlation from the repeated measures from the state. The correlation matrix R was a 57 * 57 (note 57 = 19 years * 3 age groups) symmetric matrix with diagonal elements equal to 1 and offdiagonal elements equal to . Note that such a structure is considered an exchangeable correlation matrix. Such a specification was based on empiric estimates of the correlation matrix, which showed remarkably similar correlations across years that did not diminish over time or across age groups. The values l and u for the flat prior for  were also based on the empiric estimates of the correlation matrix and were determined to 0.4 and 0.9, respectively, corresponding to the range of values within the estimated correlation matrix. Note that specifying an independent variance matrix (i.e. no correlation component) gave in very similar results (data not shown), perhaps due to the use of the year*state interaction terms to allow for state-specific temporal trends (which, in the Bayesian framework where all parameters are treated as random variables, function similarly to random intercepts and random slopes in mixed effects models, though without the distinction between fixed and random effects). Note that the use of the multivariate normal distribution provides two distinct advantages over potential binomial approaches: (1) the use of linear regression, which is often preferable in DID analyses, 11 and (2) the ability to directly incorporate correlation structures accounting for autocorrelation into the likelihood for the response variable.
The prior for  was selected to be centered roughly around the corresponding suicide mortality rate for the nation over the study period, however with much greater variance. We specified m=16 for the overall and age subgroup analyses, and used m=19, 10, 12, 28, and 8 for subgroup analyses for White, Black, Other, males, and females, respectively. The variance for the prior for  was set to a=100 (corresponding to standard deviation of 10) for all analyses. The prior for  is typical for regression analyses. The prior for k was a half Cauchy distribution, as recommended by Gelman et al., 12 given its precision with standard deviations closer to 0 and its ability to accommodate much larger standard deviations should the need arise. Note that the prior distribution for the k (assuming the null, as specified incorporating the k) has ~75% of its density between -5 and 5, over 10% of its density with values more extreme than 10 (or -10), and still over 1% of its density for values more extreme than 100 (or -100), which was deemed appropriate for the effect sizes expected with the present data.
Draws from the posterior distribution were obtained with Gibbs sampling via R2jags. We obtained 20,000 draws and discarded the first 2,000 as burn in. Convergence of the draws were assessed visually and with the Geweke, Raftery and Lewis, and Heidelberger and Welch diagnostics. Additional draws were discarded as burn in as necessary (note eTable 1 footnotes) based on chain diagnostic criteria. A summary of those tests for the posterior distribution of the DID estimator (the interaction between time period and Medicaid expansion status) are given in eTable 1. .77 a The DID estimate, just as in other analyses, represents the interaction between state expansion status and time period. b The group trend differential is equal to the coefficient for the Year * state expansion status interaction term. CrI = Credible Interval. Note that the Pr(Est>0) is a 1-tailed probability, so Pr < 0.025 is required for statistical significance.