We will use health care encounter data drawn from KPWA EHR data and health care insurance claims to construct initial states, transient states, and terminal states. Data elements include ICD-10 (International Classification of Diseases, Tenth Revision, Clinical Modification) diagnosis codes, health care service procedure codes, encounter type classifications (eg, inpatient, ambulatory, virtual visit, or ED), encounter subtypes (eg, urgent care), departments (eg, mental health), clinical specialty (eg, mental health, social worker, or psychiatry), revenue codes, and place of service codes. Visit timestamps are used in models to determine follow-up time and transition times.

Patient self-reported substance use frequency, patient self-reported SUD symptoms, and clinician-documented SUD diagnoses will be used to construct states. In the creation of initial states, diagnosis codes, service procedure codes, and revenue codes help determine whether a given health care visit was related to SUD. For instance, some initial states are defined by a new clinician-documented SUD diagnosis, while others are also defined to capture clinical workflows where patients are first identified as having substance use by the nature of a health screening (eg, substance use screening in primary care). Care states require an accompanying clinician-documented SUD diagnosis, unless the care state involves a procedure or medication that is solely indicated for SUD (eg, outpatient addiction treatment or disulfiram medication). The SUD remission state is captured using self-reported substance use frequency and SUD symptoms, and clinician-reported SUD remission diagnoses.

To obtain patient self-reported substance use frequency data, clinicians are prompted by the EHR to conduct annual screening with patients in primary care. During in-person visits, screenings are administered on paper or tablet by medical assistants, then entered into the EHR. In telephone or video visits, patients complete a questionnaire via an online patient portal. Validated substance use screening measures assess past-year alcohol, cannabis, and other drug use [32,39,40]. Patients with screening scores that are deemed to be at high risk for SUD are asked to complete a follow-up assessment for symptoms of SUD using the Alcohol Symptom Checklist and the Substance Use Symptom Checklist, respectively. These 11-item checklists capture symptoms of SUD that align with the 11 DSM-5 diagnostic criteria [41]; the cutpoint of 2 or more symptoms is consistent with SUD. The checklists have demonstrated strong psychometric properties [42-44].

SUD diagnoses include clinically documented ICD-10 diagnoses for alcohol, cannabis, opioid, stimulant, sedative, and other drug use disorder (eg, inhalants or hallucinogens), which are recorded in the EHR and claims data when clinicians conduct visits. SUD diagnoses will include those identified in the alcohol, opioid, and other drug “abuse and dependence” value sets from the HEDIS (Healthcare Effectiveness Data and Information Set) developed by the National Committee for Quality Assurance [45,46].

Death, classified as a terminal state in the MSM, will be ascertained from the Washington State Department of Health death data. These data are available for all patients, regardless of whether they remain enrolled in the health system. Following standard practice [47], death data are triangulated from several data sources. Reliability of death date (excellent, fair, or poor) is determined based on the source, the number of reporting sources, and the strength of the match in the matching process (eg, no discrepancies in identifiers).

Patient sociodemographic and clinical characteristics include both time-invariant and time-varying measures and will be used to describe the analytic sample. In addition, they may be used as covariates in MSMs. Planned demographic characteristics include age, race or ethnicity, sex, and insurance history (eg, Medicaid, Medicare, employer-based, or individual or family plans purchased through WA state’s health plan marketplace). Demographics are typically collected at enrollment and confirmed or updated at health care visits over time [48]. Socioeconomic status will be approximated from census block-level measures from the AHRQ SDOH Database [49]. Additionally, clinician-documented housing instability ICD (International Classification of Diseases) codes are from patients’ EHR data [50], and patients’ inability to pay for care is measured by a recent request to the KP Medical Financial Assistance program for low-income or uninsured individuals [51]. Medical morbidity [52,53], depression [54], anxiety [55], and serious mental illness [56] will be defined by ICD codes. AHRQ SDOH database measures will be linked to patient-level records at the county level based on patient addresses.

SUD-related health conditions will be used to describe the analytic sample and as time-varying covariates in the MSMs. These include documented diagnoses for health conditions in which substance use is always the cause (eg, alcoholic gastritis, alcohol neuropathy, accidental drug poisoning, and nonfatal alcohol and illicit drug overdoses). These are measured using ICD codes specified by the Centers for Disease Control and Prevention’s Alcohol-Related Disease Impact codes [46,57] and the HEDIS diagnosis grouping of “unintentional drug overdose.”

Separate MSMs will be run for each of the 5 subcohorts specified in Table 1. The 5 subcohorts are mutually exclusive. All prespecified states that are anticipated will be considered. States may be collapsed to address states with a small number of events, and new states may emerge if they occur with sufficient frequency (eg, concurrent receipt of pharmacotherapy and behavioral treatment and counseling). When collapsing or expanding states, we will strike a balance between clinical interpretability and reliable parameter estimation. As a general guideline, ≥40‐50 transition events from S_j to S_j per covariate helps ensure stable hazard ratio estimates and corresponding SEs, but if collapsing states compromises the clinical interpretability of the results, we will adopt a more permissive threshold of at least 10 transition events per covariate for modeling [58]. This criterion will first be applied to define states; specifically, a minimum of 40 (or 10 under the more permissive threshold) transition events is required to define a state, allowing inclusion of 1 covariate in the transition model. When transitions are sufficient to include covariates in the model, at least 40 × m transition events are required. Covariate selection for each transition model, given the number of transitions, is described below. Follow-up time or care episode duration is defined as the interval from entry into the initial state (time 0) to entry into a terminal state.

Figure 1 provides a graphical representation of a simplified MSM illustrating potential care pathways. In this example, patients who meet the inclusion criteria and are identified as having an incident SUD in primary care enter the initial state C₁. Patients then transit into 1 of 15 care states (S₁, ... ,S₁₅) defined by a service type and setting. From these states, patients may transition to temporary discontinuation of care, D, or an SUD remission state, R. Transitions between the transient states (care states, temporary discontinuation, and SUD remission) are bidirectional, whereas any transition to a terminal state (T₁, ..., T₄) is unidirectional, reflecting that patients cannot leave terminal states once they are entered in MSMs. At the end of follow-up, patients occupy 1 of the 4 terminal states.

We will adopt a semi-Markov assumption [59] under which transition intensities depend on only the duration that patients spend time in the current state before transitioning. This framework allows transition intensities to vary over time rather than assuming constant transition intensities as in the Markov framework. Time is defined as the duration since entry into a state, with the clock reset to 0 upon entry into the cohort and each subsequent state. At any time since entry into the initial state, a patient occupies 1 of the 21 states (17 transient states and 4 terminal states). Movements between states j and k (j,k=S1,…,S15,D,R,T1,…,T4), are governed by a 21-by-21 matrix of transition intensities qjk(t,z(t)), which is a function of time t and covariates z(t). Intensities capture the instantaneous probability of moving from state j to k.

Under the semi-Markov assumption, we will fit separate flexible parametric survival models [60,61] for each transition:

log⁡{qjk(t,z(t))}=s(log⁡(t)∣η,k0)+βTzi(t).

For simplicity, we omit transition-specific subscripts for the parameters (η,k0,β). In the model, s(log⁡(t)∣η,k0)=∑i=1k0ηiBi(log⁡(t)) represents the log cumulative baseline intensity function, which we model flexibly using restricted cubic spline basis functions (also known as natural cubic splines) Bi(⋅), with corresponding coefficients ηi. The number of knots, k0, determines the number of basis functions and thus the flexibility of the baseline intensity functions. For covariates z_i(t), the coefficients β represent the log hazard ratios at time t. The unknown parameters to be estimated or selected from data are (η,k0,β). We will use the Akaike information criterion based on data to select the number of knots, with lower Akaike information criterion values indicating better model fit. The regression coefficient parameters (η,β) will be estimated by using the maximum likelihood approach for statistical inference [62,63]. For each transition, we will include prespecified covariates for adjustment when the number of transitions allows (at least 10 transition events per covariate), and if additional covariate selection is required, variables that are statistically significant at the 5% level based on Wald tests will be included [64]. The estimated transition intensities are transformed to a transition probability from states j to k at time t since entry into state j at time s (s<t, s and t>0), by integrating, pjk(s,t)=exp⁡{−∫stqjk(u,z(u)) du}.

Multiple incident SUD events may be observed within the same patient (see Quantitative Study Population). For example, a patient may have an incident SUD event, eventually experience a terminal state such as ceasing SUD care (eg, one who receives no SUD care for 194 days), and then experience a second incident SUD event. Including multiple incident events from the same patient within the same subcohort in the MSM framework may introduce within-patient correlation and reflect care pathways that differ from those of the initial episode. We will account for within-patient correlation using a sandwich variance estimator for the parameter estimates. The number of prior incident events will be included as a covariate or used as a stratification factor, depending on its distribution.

To assess model adequacy, we will compare spline-based estimates of the cumulative baseline hazard with nonparametric estimates derived from the Kaplan-Meier method for visual agreement. To evaluate the proportional hazards (time-independent effects) assumption, we will test for time-dependent effects by including interactions between covariates and spline functions of time. If these interaction terms are not statistically significant at the 5% level, they will be excluded from the final model.

Additionally, we assume that transition times are exactly observed. Lastly, it is assumed that the distribution of health care visit times is independent of transition probabilities for simplicity (noninformative visiting process). We will primarily use the R (R Foundation) package flexsurv to implement the MSMs and develop our own modules to transform the fitted results to answer research questions.

We will consider sensitivity analyses if sample size and other feasibility constraints permit. This includes stratified analyses by subgroups defined by the number of outpatient visits (eg, patients with more visits may have poorer health) [65,66] or clinically meaningful periods. For instance, impactful operational changes occurred due to COVID-19 mitigation restrictions (eg, reduced in-person care), so we will consider estimating models before, during, and after these restrictions were in place. Additionally, stratified analyses by SUD type will be considered to acknowledge that each substance may sometimes require unique forms of care, such as pharmacotherapy usage.

We will apply the estimated MSM parameters from aim 1 to generate inferences that yield actionable data-driven insights to improve care delivery for SUD. As shown in Table 3, these inferences fall into three general categories: (1) examine SUD care states and terminal event rates, (2) identify empirically common care pathways, and (3) identify the extent to which guideline-concordant care is achieved.

Table 4 summarizes estimates that can be obtained from fitted MSMs to generate our inferences and answer additional research questions. We will calculate point estimates and CIs and further characterize variation in output quantities across states graphically.

How does the accessibility of each state vary? We will use the models to describe how the accessibility of each state varies by estimating the dynamics of transitioning between different types of care for SUD. We will report and characterize variation in (1) transition probabilities, (2) hitting probabilities, (3) mean sojourn time, (4) expected first passage time, and (5) expected number of visits (definitions in Table 4).

How does the probability of occupying a given care state or terminal state vary over time? We will describe how the probability of occupying a state varies over time since cohort entry. To address this question, we will construct a stacked plot in which the y-axis represents the probability of being in each state, the x-axis shows time since cohort entry, and the stacked bands correspond to different states (Figure 2). Depicting the timing and probability of occupying states can help inform quality improvement interventions. For instance, by visualizing the time at which the probability of ceasing SUD care increases, it can help guide the timing of interventions to check for SUD care engagement and outreach to re-engage patients when needed.

What are inefficient sequences of SUD care following each initial state over time? Several inferences from the MSM may be a signal for inefficient SUD care. First, prolonged time in some states, such as inpatient care, may reflect challenges in coordinating care with services following an inpatient stay. We will thus identify states with relatively long mean sojourn time [67]. Second, care states with a low probability of being reached, or that take a long time before they are reached, may reflect barriers to care such as long wait times. To identify these care types, we will explore states with long expected first passage times [68] and low hitting probabilities [67]. We will also construct survival plots to display how the probability of starting in an initial access point (eg, index SUD event in primary care) and not reaching a target care state changes over time. For example, Figure 3 displays the probability of not reaching a target care state (eg, behavioral treatment and counseling in an outpatient setting) from different initial states. In this example, the probability of not reaching behavioral treatment and counseling in an outpatient setting following entry into the ED cohort is markedly higher than the corresponding probabilities following entry into the primary care or MHAC cohorts, suggestive of referral inefficiencies in ED services.

What is the contribution of care states in achieving terminal states? We will calculate restricted mean survival times (RMSTs) [69] between a given SUD care state and each terminal state. An example of RMST quantities between select care states and a terminal state is presented in Table 5. The RMST captures the average terminal-event-free time, which in this study reflects time before reaching a terminal state, calculated as the area under the terminal-event-free curve over a prespecified time interval since the entry to the given care state. Figure 4 presents an example of terminal-event-free curve estimates using states from Figure 1. In this example, care types that produce lower RMSTs are those more likely to end in ceasing SUD care. Identifying these care types can inform targeted efforts, such as enhancing care coordination or improving staff levels in specific settings or for patients who receive certain treatment modalities to increase retention in SUD care.

What are common sequences of care following each initial state? Care pathways consist of a sequence of transitions across care types and settings. To characterize these sequences, we will calculate the probabilities of all identified pathways starting from a given initial state in Table 1, proceeding through multiple care states, and ending in a specific terminal state. We will estimate probabilities of sequences of varying length. The length of sequences will be informed by descriptive analysis examining the mean and median number of transitions among subgroups of patients achieving terminal outcomes in Table 2. Sequence probabilities will be estimated using the fitted MSMs under hypothetical scenarios, such as repeated cycling between care and care-discontinuation and specified transition timings (see Multimedia Appendix 2).

How often and where is guideline-concordant care achieved? We will analyze several health care quality measures, including those published by the National Committee on Quality Assurance [36,70,71] and others [72,73]. These monitor initiation of care, or the receipt of prespecified health services following an incident SUD event; engagement in care, defined as having a minimum number of qualifying services after initiating care; and adequacy of treatment, which captures whether quantity or duration of treatment over time is sufficient. Achievement of these measures can be ascertained from the states in the MSM model. Multimedia Appendix 1 contains details about how we will operationalize these quality measures using estimates from the fitted MSMs. For example, a National Committee on Quality Assurance quality measure monitors follow-up for patients who have an acute care encounter for SUD. This measure monitors the percentage of patients who obtain follow-up SUD care after an ED visit over 7-day and 30-day timeframes. Thus, among patients in the ED and hospitalizations subcohort whose incident SUD event is in an ED setting, we will estimate transition probabilities capturing the likelihood of transitioning from the initial state to any SUD care state within 7-day or 30-day. We will estimate measures reflecting care engagement and adequacy of treatment by estimating the expected length of time in a pharmacotherapy state or the expected number of visits to qualifying states over the duration of time specified in Multimedia Appendix 1.

How does the probability of occupying terminal states differ across patient subgroups? Within study subcohorts, we will conduct descriptive analysis comparing occurrences of terminal states across patient subgroups. For instance, we will compare terminal states at selected time points (eg, 1-year after cohort entry) across mutually exclusive subgroups defined by patient SUD type (eg, sole opioid, alcohol, or cannabis use disorder vs multiple SUDs). Comparisons of terminal states will help provide insights into whether outcomes differ for specific groups of patients who have incident SUD events in each setting [28].

How do answers for the above key questions differ across patient characteristics? First, overall quantities will be estimated to address the key questions in Table 4. Then, Cox proportional hazards models will be fit for each transition intensity to assess associations between patient characteristics, such as demographics, SUD type and severity, mental health diagnosis, socioeconomic status, and transition intensities [9]. Using the fitted model, the quantities required to address the questions in Table 4 will be calculated as a function of patient characteristics included in the transition models.

Aim 3 includes observations and interviews to explore the context for SUD care decisions and experiences with SUD care transitions. We will perform direct observations [74,75] of patient-clinician encounters during salient SUD visits, such as behavioral health intake visits. Witnessing and documenting these interactions will provide key information on decision-making regarding various clinical scenarios. We will also conduct semistructured in-depth interviews with a broad sample of patients and clinicians to collect views, preferences, and needs of individuals engaged in different care pathways. Observations and interviews will be concurrent to inform ongoing data collection and confirm data sufficiency [76]. We will draw on phenomenological [77] and ethnographic theory [78] to explore emic [79], or insider perspectives on care transitions.

The quantitative and qualitative components synergize in 2 ways: first, by using an explanatory sequential sampling mixed methods design [80] whereby qualitative interview participants are purposely sampled from potentially meaningfully rich subgroups of patients who have embarked on specific care pathways identified in the quantitative population. Second, by integrating quantitative and qualitative data in the analysis phase. Specifically, we will reapply the inclusion criteria for aims 1 and 2 to identify patients who were recently identified to start an episode of SUD care at the time of qualitative data collection. We will then match the patient’s observed sequence of care to model-identified pathways estimated in aim 1. We will then purposively sample [81] patients from six pathways: (1) two of the overall most common care pathways, (2) two of the most common care pathways that contain guideline-concordant care, and (3) two of the most common care pathways that have a high probability of reaching negative terminal states (eg, death). We will be open to the emergence of other important care patterns to guide our sample selection, and we may consider other important factors, such as SUD type and/or severity, when defining these groups. We will sample approximately 3‐5 patients from each pathway (estimated to be 35 patients) and ensure representation across demographic groups and those with different SUD types.

For the patient interviews, a semistructured interview guide will include three main topics: (1) care choices and constraints at critical junctures (eg, referral to specialty services), (2) care alignment with patient needs and preferences (eg, care frequency and location), and (3) clinical and care quality outcomes (eg, treatment adequacy). Interviews will use a modified form of free listing [82-84], an ethnographic technique increasingly used in health research [85-87] to elicit insider views and conceptions of a cultural or experiential domain (eg, process of obtaining SUD care). We will ask patients to rapidly say or write words that come to mind when thinking of an experience (eg, patient decision-making). Their words are assumed to indicate salient constructs within a domain (eg, treatment options) common to people who share the experience and help reveal the conceptual boundaries of the domain of interest. This provides a foundation for further inquiry into health-related experiences, attitudes, beliefs, and behaviors [88]. We will query how patients conceive of and define care journeys, potential choices, and instructions from clinicians. We will draw from the Chronic Care Model [29] to inquire how the free-listed words affected the ability to achieve outcomes, and how clinicians and systems can better support care transitions.

For clinician interviews, we will interview employees of KPWA or external organizations in the contracted care network who help patients make SUD care decisions in primary care, ED, urgent care, behavioral health, and specialty SUD care. We will oversample clinicians who provide care in roles that are consistent with the model-identified common care pathways and ensure representation across clinician disciplines and care settings. We will recruit clinicians over email to arrange a video interview [89]. We anticipate recruiting approximately 25 clinicians. Guided by the Chronic Care Model, interviews will focus on decision-making and system- and practice-level factors affecting quality and coordination of SUD services. Questions will consider practice guidelines and information systems used to determine treatment plans, and how delivery system design affects referral choices. Per relevant literature, we will probe barriers and facilitators to care coordination (eg, program capacity), patient engagement, and meeting care needs. We will query beliefs and experiences tied to the role clinicians play in enhancing care continuity, and how roles can be strengthened [90]. All interviews will be conducted by phone or videoconference platform, recorded, and professionally transcribed.

We will also conduct patient-clinician observations during visits where SUD treatment planning takes place, such as scheduled SUD intake evaluations, primary care visits that follow ED or inpatient discharge, addiction psychiatry consultation, and buprenorphine induction. We will conduct up to 20 observations (≈5 per visit type), monitoring for data saturation, to capture nuances of the patient experience at critical transition points [91]. We will contact clinicians via email to invite them to have their visits observed. When a clinician agrees, a researcher will be available when the clinician begins the scheduled video or phone-based session. The clinician will ask the patient if a researcher could join to observe. If the patient agrees, the researcher will join the video call, invite the patient to participate in study observations, and obtain verbal informed consent. A semistructured observation template will be developed for taking field notes to document the context in which care options are considered and decided upon, and assess multilevel factors affecting care planning and pathways pursued. The template will be organized by the Chronic Care Model constructs and allow for unstructured notes and reflections. After each data collection (interview or observation), qualitative researchers complete a structured memo to identify data gaps and challenge assumptions. Prompts such as “what surprised you” and “what is unanswered” help the researcher reflect on initial assumptions and reveal data gaps to explore. After 5 data collections, memos are summarized and discussed to monitor data quality and completeness. Physician and social work investigators, partners, and consultants will join debriefs as data collection unfolds to identify new questions to be pursued.

A total of 2 researchers will analyze all qualitative data using thematic analysis [92]. We will code transcripts and observation notes using deductive and inductive codes [93]. The initial code list will be developed from interview and observation guide questions and prompts and Chronic Care Model domains. Inductive codes will be developed from reading transcripts and notes. Further, 2 researchers will read the same subset of transcripts to derive emergent codes, discussing and resolving discrepancies by consensus. We will combine all the codes, refine the code list via an iterative process, and rely on the Atlas Ti qualitative platform [94] to apply the final codes to transcripts and observation notes. Lastly, we will examine the coded text to highlight key themes and subthemes and use coding memos to summarize and further refine themes. Coding memos will form the basis for disseminating study results. We estimate that our patient and clinician samples will ensure the analytic depth, rigor, and trustworthiness of the findings, but we will monitor thematic saturation during analysis and recruit additional participants if needed [91,95]. We will triangulate interview and observation data to interrogate qualitative findings [96,97] and reach robust conclusions. Then we will integrate qualitative and quantitative findings by constructing a joint display of key qualitative themes and salient estimated quantities [80]. The joint display will highlight patterns, reveal concordant or discordant findings, and expand overall interpretation [98].