Statistical Analysis

We aimed to estimate the association between the duration of MOUD and risk of an opioid-related overdose. To do so, we implemented landmark survival analysis⁵¹ which combats a common source of bias in traditional survival analysis, immortal time bias.³⁶ Immortal time bias arises when treatment assignment at T₀ depends on a behavior or exposure after T₀ when the outcomes are assessed. For example, one might compare the hazard of overdose from T₀ through 12 months for a treatment group and a comparison group that are defined by MOUD duration from T₀ forward (e.g., >= 120 days and <120 days). Treatment subjects have a survival advantage, or “immortal time” of 120 days, because by design they cannot develop the outcome before receiving the treatment.

Landmark survival analysis mitigates immortal time bias by identifying treatment or comparison group status independent of the subject’s exposure to MOUD during the outcome assessment period and allowing study group status to vary within person over time. Specifically, for each subject we defined their treatment group (i.e., continuing MOUD), comparison group (i.e., discontinued MOUD) or censored status at each landmark time point, T_n, relative to the index MOUD date. For each landmark time point, we implemented a separate Cox proportional hazard model to assess the relative hazard of overdose among those who continued MOUD beyond T_n relative to those who discontinued by T_n. The analytic sample for each model included only subjects who had not experienced the outcome and remained eligible for treatment or comparison status up to that landmark time point. In other words, the comparison and treatment group subjects at each landmark time point had comparable immortal time. This approach improved the exchangeability of study groups, because all subjects, whether identified as treated or comparison at T_n, were equally likely to continue as participants to T_n.

We used a two-stage procedure to conduct these analyses. First, all states ran identical analytical SAS code (version 9.4, SAS Institute, Inc.,) on their respective data that had been standardized according to the MODRN common data model. We excluded subjects with missing data for rural/urban living area from the analytical sample.

In the second stage, we conducted random effects meta-analyses to pool the state-specific estimates and generate global estimates adopting methods validated in similar settings and previously deployed in the MODRN.^37,51,52 We generated global estimates of the hazard ratio of MOUD continuation relative to discontinuation at each landmark time point by averaging the model estimates from each state weighted by the inverse of their variances to account for differences in population size. We used Cochran’s Q to measure and test the statistical significance of between-state heterogeneity in the estimates. To construct valid 95% confidence intervals for the global estimate we used the Hartung-Knapp-Sidik-Jonkman method to estimate between-state variances.⁵³ We tested whether the global estimate equaled one using a 2-sided adjusted t-test and a significance threshold of 0.05.⁵⁴ We allowed for between-state variation in the effect of MOUD continuation because we had only a selection of states, but wished to extend our inference to similar states outside of our sample. To do so, we constructed 90% prediction intervals.⁵⁵ The prediction intervals convey state-to-state variability denoting the range within which the hazard ratio would fall for 90% of the states if we drew a different sample of states (eAppendix 4 in Supplement). For this reason, prediction intervals remain relatively stable to the number of states included in the pooled analysis, providing an objective measurement of uncertainty, while confidence intervals (i.e., from fixed-effect pooling) would continue to decrease as more states were included.

We conducted a test for trend over the 12-month follow-up period in the global hazard ratio estimates using random-effects meta-analysis with the inclusion of time. The landmark time point was included as a linear fixed effect, and a random effect of state was also included.

We tested the sensitivity of our results to the inclusion of the type of MOUD at index by re-estimating the analyses including the states that offered all MOUD types throughout the study period. Finally, we estimated the E-values for our main results,⁵⁶ the global HR estimates, to explore the degree to which our main estimates could be explained by unobserved confounding variables. Finally, we explored the moderating role of MOUD type in the association between MOUD continuation and overdose risk. We implemented stratified analyses in the four states that covered all MOUD types in all years and had sufficient sample sizes to compare overdose risk among continuers and discontinuers by MOUD type at all time points. The study’s analytical plan was not pre-registered.