Designing and validating a Markov model for hospital-based addiction consult service impact on 12-month drug and non-drug related mortality

ACS Referral

Once patients are admitted to the hospital, they can be referred to ACS care. ACSs exist across a growing number of North American hospitals. Typically, they include care from an interprofessional team that may include of medical providers, social workers, nurses, and alcohol and drug counselors (17). Some intentionally include people with lived experience in recovery (18-20). ACSs typically address needs of people use any substance (for example, stimulant, alcohol, and opioids) and care includes comprehensive assessments, withdrawal management, medication treatment, psychosocial and harm reduction interventions, and efforts to support patient engagement and linkage to care across settings. ACSs commonly also provide staff education and patient advocacy (17, 21, 22). For this model, we used an intention-to-treat approach; all patients referred to ACS were included regardless of level of care engagement or specific services received.

Post-discharge OUD treatment engagement

We used a modified Healthcare Effectiveness Data and Information Set (HEDIS) measure of engagement to stratify for post-discharge OUD treatment engagement. The original measure requires that patients initiate treatment and have two or more additional alcohol or drug services or medication for OUD within 34 days of initiation (23). Recent research has shown that evidence-based MOUD has superior outcomes in preventing mortality and decreasing opioid use (7). For this reason, we defined post-discharge OUD treatment engagement as: 1) at least two filled prescriptions for buprenorphine, extended-release naltrexone, or methadone from an Opioid Treatment Program in the 30 days following hospital discharge, or 2) a prescription for extended-release naltrexone or buprenorphine that covered 28 of the 30 days post-hospital discharge (24).

12-month mortality

At twelve months, deaths are classified as drug related versus non-drug related (including as circulatory, neoplasm, infectious, digestive (including alcohol-related liver disease), external (including suicide and unintentional injury), respiratory, endocrine, and other) by ICD-10 mortality codes described by Hser et al (2).

Cohorts for transition points

We defined three cohorts for our analyses utilizing Oregon Medicaid data. First, we included all patients who met eligibility criteria in analysis for our first transition, referral to ACS. Then, we used a matched cohort of three controls to one IMPACT patient for our post-discharge OUD care engagement and mortality analyses. We matched without replacement on hospital admission quarter and admission number, including one admission per person.

Bayesian priors via expert elicitation

We used expert elicitation to capture prior information for our models. We identified important covariates at each transition point, including age (in years), gender (female/male), race (White/not White/unknown), ethnicity (Hispanic/Not Hispanic), concurrent alcohol use disorder (yes/no), concurrent stimulant use disorder (yes/no), hospital length of stay (in days), rural residence (yes/no), filled at least one prescription for medication for OUD in the month before hospital admission (yes/no), previously admitted to the hospital (yes/no), and Chronic Illness and Disability Payment System (CDPS) Score (continuous). The engagement model also included referral to an ACS (yes/no). The mortality models included engagement in care after discharge (yes/no) and filled a naloxone prescription in the 30 days after hospital discharge (yes/no).

We used a clinical-vignette design to ask providers about the relevance of covariates on patient outcomes. To do this, participants provided a probability estimate for different events: referral to an ACS, post-discharge engagement, and mortality.

For example, a vignette could read:

“The patient is a young White man with OUD and AUD. He was in the hospital for several days. He was on medication for OUD at admission. He had never previously been admitted to the hospital. He has many comorbidities. He is not from a rural area. What is the probability he engaged in post-discharge treatment for OUD within 30 days of discharge?”

Experts evaluated 16 (referral to ACS), 17 (engagement) and 18 (mortality) vignettes selected from an optimal experimental design generated for each model (29). From the optimal design, we chose a subset of the vignettes that were substantially different from one other for ease of interpretability and to maximize the information gathered about each covariate.

As part of our IRB-approved research, study authors (HE, PTK) generated lists of experts in addiction consult services and hospital-based addiction treatment in general in the United States. Each participant took only one survey. We aimed to recruit at least five participants for each survey, with a goal of at least three responses per survey. For the referral to ACS survey, we also asked participants to refer hospitalists at their institutions to complete the survey, as hospitalists are frequently providers who refer patients to ACS. Ultimately, six participants took the ACS survey (6 of 11, 54.5%), four took the engagement survey (4 of 5, 80%), and three took the mortality survey (3 of 8, 37.5%).

Bayesian logistic regression models

We used the transformed prior information from expert surveys and Oregon Medicaid cohort data to fit Bayesian logistic regression models at each transition point. Models were fit using Markov Chain Monte Carlo methods (30). We sampled each parameter 10,000 times with 2000 burn-in chains. We used multiple metrics to assess model convergence. First, we used Gelman and Rubin’s potential scale reduction factor; all values in all models equal 1.0. Values close to 1.0 are suggestive of convergence. Effective sample sizes all approximated the number of posterior draws requested. All model trace plots appear to have a caterpillar-like distribution, and there were no divergent transitions. Autocorrelation plots for all parameters suggest low autocorrelation. We used the package Shiny Stan to evaluate Bayesian model fit (31).

We tested different prior information strengths: first, using a cohort sample size method, where the prior information equivalates a percent of the study sample size (0.1%, 1%, 5% and 10%); second, using a confidence interval method, where we fit beta distributions to the range of survey responses, and then used the maximum and minimum values as borders for 80%, 85%, 90%, and 95% confidence intervals. We picked the best-fit model using Pareto smoothed importance-sampling leave-one-out cross validation using the loo package in R where lower expected log predictive density values indicate a better model fit (32). We also prioritized models where Pareto k diagnostic values had at least good reliability for all estimates.

We used mcmcObsProb in the BayesPostEst package (33) to estimate marginal transition probabilities over observed cases with the fitted Bayesian logistic regression models. We created prior-posterior plots using ggplot2 (34).

Face validity

To assess face validity, one researcher (CK) designed the model and received feedback from experts in addiction medicine outside of the study team about the model’s face validity. Experts agreed that the model reflected the path of care for patients admitted to hospitals in Oregon with OUD (structure). Further, the use of Oregon Medicaid data, versus data from the literature, was considered a strength in deriving evidence for the model by outside experts. ACS and their impact on care for patients with OUD is of immense interest to healthcare systems and policymakers, and experts also agreed that the question was timely and important (problem formation). Finally, after data analysis, the model results were presented to researchers who agreed that estimates from the model matched their expectations (results).

Internal validity

We conducted additional checks and analyses to ensure internal validity of our Bayesian approach (also referred to as technical validity, (36)). First, a recent paper used a similar approach and data structure to evaluate the impact of prenatal maternal factors on nonadherence to infant HIV medication in South Africa. After building our Bayesian model, we used the deidentified data from the South Africa analysis to attempt to replicate identical results as were published. The built model exactly replicated the results of the South African analysis. Second, we conducted classic logistic regression models for each transition point in addition to the Bayesian models. We placed a 1/3, 1/3 noninformative prior (Kerman’s prior) on all covariates, which should be roughly approximate to the classic logistic regression results. Our results with non-informative priors were sufficiently similar to classical logistic regression results. Finally, we conducted code “walk throughs” as suggested, where the analyst (CK) walked through code with an expert in these methods (RC).

In addition to the above steps, because we used Bayesian analyses for our transition probabilities, we needed to ensure that our final estimates of confidence intervals around engagement and mortality estimates actually encompassed the observed number of people who engaged, and people who died from drug-related and non-drug related deaths. Of the 6,654 patients with 12 months follow-up time, the model estimates that 1,330.8 patients engage in care (95% CI =1,064.6, 1,597.0). We observed 1,318 patients who engaged in care in the cohort. Additionally, the model estimated 357.2 drug related deaths (95% CI= (98.5, 632.6)); there were 114 observed drug related deaths in the dataset. Similarly, the model predicted 570.8 non-drug related deaths (95% CI= 263.6, 865.0)); there were 408 observed non-drug related deaths in the dataset. Mortality analyses rarely account for all sources of follow-up which may mean that reported mortality estimates in the literature are lower than in reality. Thus, it was not surprising that modeled transition probabilities from Bayesian logistic regression for 12-month mortality may be higher than raw observed proportions.

Cross-validation

Researchers at a separate academic medical center have developed, validated and calibrated the Reducing Infections Related to Drug Use Cost-Effectiveness (REDUCE) model, a Monte Carlo microsimulation model (37). This model has the capacity to answer similar questions to what we post here, using estimates derived from published data and from expert sources. In contrast to our model which uses a cohort defined by opioid use disorder, the REDUCE model simulates data for people who inject drugs. Because model estimates for the REDUCE model are derived from a variety of sources in different parts of the county, we expected outcomes from the REDUCE model to be different from our model; we felt these differences are important to understand.

To support cross-validation of our model, the research team that developed the REDUCE model generated 4,153 simulated patients admitted to the hospital for the first time. Of those, 36 died while in the hospital (0.9%). Of the 4117 still alive at hospital discharge, 96 (2.3%) died within 12 months of hospital discharge (95% CI = 1.9%, 2.8%). This is lower than our estimated 928 (13.9%) deaths from our Markov model (95% CI = 5.4%, 22.5%).

There are several important differences between the REDUCE model and our model. First, as previously mentioned, the REDUCE model simulates data from patients who inject drugs, while ours models patients who have OUD more generally. There are important demographic differences between these two groups, including that our model also includes patients with a primary diagnosis of cancer. Next, the percentage of people seen by an ACS in the REDUCE model was higher than in our model: 25% of patients in REDUCE were seen by an ACS versus 4% in our model. The REDUCE model uses data from Boston, where higher numbers of patients are seen by ACS. This makes it challenging to understand REDUCE estimates in the context of Oregon specifically. Additionally, patients had a higher post-discharge treatment engagement rate in the REDUCE model. In REDUCE, approximately 25.2% of patients receive medication for OUD for at least one week in the month following discharge, versus our model, where 20% of patients not seen by an ACS receive MOUD after discharge. Finally, data from the first simulated admission was used to estimate 12-month mortality from REDUCE; because we matched our cohort controls on the number of previous admissions among patients seen by an ACS, it is possible that our patients were older and sicker than patients who had never previously been admitted to the hospital. While the base model structures are similar, our model is populated with data that provides a focused understanding of addiction consult services in Oregon. Populating our model with different data, including Boston estimates, could provide tailored explorations of ACS in different settings.

External validity

To examine external validity, we used large, high-quality, recent studies of representative populations in independent cohorts of participants to separately validate post-discharge OUD treatment engagement and 12-month drug related and non-drug related mortality. We simulated a cohort of size determined from outside research and looked to see if our simulated confidence interval (cohort simulation/matrix multiplication method, (36)) was different from observed values or confidence intervals from the published estimates (Table 2). Where there was disagreement, we describe potential causes.

Predictive validity

All relevant data was included in building the Markov model described in this paper. We have planned analyses to evaluate our model predictions versus Medicaid claims data for the same cohort of patients seen in through December of 2020, once data is released.