Evaluating the population-level effects of overdose prevention sites and supervised consumption sites in British Columbia, Canada: Controlled interrupted time series

Background On 14 April 2016, British Columbia’s Provincial Medical Health Officer declared the overdose crisis a public health emergency, sanctioning the implementation of new overdose prevention sites (OPS) and supervised consumption sites (SCS) across the province. Methods We used the BC Centre for Disease Control’s Provincial Overdose Cohort of all overdose events between 1 January 2015 and 31 December 2017 to evaluate the population-level effects of OPSs and SCSs on acute health service use and mortality. We matched local health areas (LHA) that implemented any site with propensity score matched controls and conducted controlled interrupted time series analysis. Results During the study period, twenty-five OPSs and SCSs opened across fourteen of British Columbia’s 89 LHAs. Results from analysis of LHAs with matched controls (i.e. excluding Vancouver DTES) were mixed. Significant declines in reported overdose events, paramedic attendance, and emergency department visits were observed. However, there were no changes to trends in monthly hospitalization or mortality rates. Extensive sensitivity analyses found these results persisted. Conclusions We found OPSs and SCSs reduce opioid-related paramedic attendance and emergency department visit rates but no evidence that they reduce local hospitalization or mortality rates.


Study design
This is a controlled interrupted time series study using data from the retrospective BC overdose cohort for the period 1 January 2015 to 31 December 2017.

Setting
British Columbia is Canada's most western province and has a population of approximately 4.8 million residents. The province provides single payer coverage of inpatient and outpatient health services through its Medical Services Plan. Residents excluded from the insurance program include newly landed immigrants and people covered under federal insurance programs including refugees, asylum seekers, military personnel and First Nations' members (representing less than 4% of the population) [19].
In 2012, fentanyl was first detected in the illicit drug supply, and 4% of the province's 270 overdose deaths were fentanyl-related. By 2019, fentanyl was detected in over 85% of the province's 984 drug overdose deaths [1]. Here, illicit overdose events include indication of street drugs (controlled and illegal: heroin, cocaine, MDMA, methamphetamine, illicit fentanyl), and medications not prescribed to the decedent but obtained/purchased on the street, from unknown means, or where origin of drug not known.
To test our hypotheses, we compared outcomes at the local health area (LHA) level. LHAs are a mutually exclusive and exhaustive classification of the land area [20]. Although they have no administrative functions, until 2019, LHAs were the smallest geographic boundaries used for health services planning and delivery [21]. There are eighty-nine LHAs ranging from 2,000 residents in sparsely populated remote communities to 485,000 in heavily dense urban centres. BC's largest city, Vancouver, is comprised of six LHAs ranging from 63,000 to 145,000 residents. During the study period, fourteen LHAs implemented the intervention, and seventy-five did not.

Study population and data sources
We used the BC Provincial Overdose Cohort [22][23][24][25][26] which captures all coroner confirmed drug-related mortality and non-fatal opioid-related overdose events involving health service provider interactions that occurred in British Columbia. Because the province's coroner services are required to investigate and determine cause of death for all "unnatural, sudden and unexpected, unexplained, or unattended deaths," [27] there is a two-year lag between confirmed opioid-related mortality and access to data for research purposes.
The cohort included all events treated with naloxone administered by paramedics; calls to the Drug and Poison Information Centre about an opioid-related event; coroner-determined illegal drug overdose deaths; hospital admissions with ICD-10 code T40.0, T40.1, T40.2, T40.3, T40.4, or T40.6; emergency department visits with ICD-10 code T40.0 or T40.6; and outpatient physician visits with ICD-9 code 965.0 or E850.0. Non-fatal opioid-related overdose events treated with naloxone by a peer in the community or staff at an OPS or SCS were not included unless there was follow-up with a health services provider. Although BC Coroners Service data include deaths involving any street drugs (heroin, cocaine, MDMA, methamphetamine, illicit fentanyl etc.), as well as medications not prescribed to the deceased, combinations of the preceding with prescribed medications, and those overdoses where the origin of the drug was not known, we restricted our analysis to deaths involving opioids (heroin, codeine, oxycodone, morphine, hydromorphine, methadone, fentanyl and analogues, etc.).
For individuals with multiple opioid-related overdose events during the study period, each event was included separately in the data set. Event location was assigned as the LHA where the event initial occurred or where the deceased was found if no previous contact with health services was made. Detailed description of the cohort is available elsewhere [22].

Outcomes
The primary outcome was coroner confirmed mortality. Secondary outcomes were opioidrelated health service encounters defined as hospitalizations using admissions records in the Discharge Abstract Data base, emergency department visits captured in the National Ambulatory Care Reporting System, and paramedic attendance in BC Emergency Health Services data holdings.

Data preparation and analysis
We estimated the population-level effects of the newly implemented OPS/SCSs using propensity score matched interrupted time-series analysis. This quasi-experimental study design enabled us to estimate the effects of the intervention despite the lack of randomization by geography [28].
For each LHA that implemented the intervention we selected a matched counterfactual from the pool of LHAs that did not begin operating any OPS or SCS within twelve months of the exposed LHA's implementation month to account for confounding by indication. Owing to constraints with the data, we worked with LHA characteristics plausibly predicting the implementation of OPS or SCS (i.e. population-level age, sex, and income demographics and opioidrelated overdose mortality rate). Low income residents were those covered by PharmaCare Plan C-the provincial drug plan available only to individuals and families receiving income assistance. We applied matching with replacement using a maximum-likelihood logit model to calculate each LHA's propensity score [29]. The matching algorithm used these scores and a caliper of width equal to 0.15 of the standard deviation of the logit of the propensity score to create two evenly matched groups while accounting for the variance-bias trade-off [30].
For each outcome, we organized exposed-control LHA pairs by study time and created twelve weighted monthly outcome rates per 100,000 population by exposure group pre-and post-implementation, censoring on the month of implementation. For LHA-level population size, we used Statistics Canada Census data (periods 2011, 2016) and linear interpolation and extrapolation to estimate monthly values during the study period. Where LHAs implemented more than one OPS or SCS during the study period (e.g. Vancouver's Downtown Eastside, Victoria, and Surrey), we set study time = 0 as the operation month of the first new OPS and considered subsequent OPS/SCS openings within the LHA as 'scale-up'. As part of our sensitivity analysis, we assigned the LHA of the first point of medical contact (e.g. hospital) where opioid-related overdose event location was missing.
We inspected time trends for autocorrelation using the 2-sided Durbin-Watson test and visual plots of the autocorrelation and partial autocorrelation functions [31,32], adjusting where necessary. We tested our hypotheses using ordinary least squares and segmented regression, by fitting the following regression model, per outcome: where j was the intervention, t was the study time in monthly intervals pre-(negative time) and post-intervention (positive time), and k distinguished between intervention and control group. Significant values for coefficients β 6 and β 7 indicate an effect of the intervention on exposed LHA after controlling for level and trend changes in the controls, respectively. Where a matched control could not be found for an exposed LHA, we analyzed this LHA separately using a traditional interrupted time series analysis. Data were prepared using SAS 9.4 [33]. All statistical analyses and graphics were conducted in R 3.6.1 [34] using the nlme [35], car [36] and ggplot2 [37] packages.
McGill University's Institutional Review Board (Certificate Number: A12-E79-18A) and the University of British Columbia's Research Ethics Board (Certificate Number: H18-03361) approved this study. This study only uses secondary data and did not require consent from participants included in the datasets.

Results
Between January 2015 and December 2017, there were 36,576 unique overdose events in British Columbia. Of these, 26,223 (71.7%) were attended by paramedics, 24,171 (66.1%) included a visit to an emergency department, 3356 (9.2%) events resulted in a hospitalization, and 3604 (9.9%) resulted in mortality (not mutually exclusive events). Of the 36,576 overdose events, 5836 (16.0%) were missing the event location. The majority (76.3%) of events missing location information were first captured in emergency department (ED) visit records, and are likely indicative of patients who self-transported to the ED. When we restricted our analysis to overdose events with complete location information, we were left with 30 During the study period, fourteen LHAs implemented at least one OPS-some of which have since transitioned to SCSs (Table 1). When we matched with replacement, we found eleven control LHAs for thirteen exposed LHAs (Tables 1 and 2). Propensity score matching did not identify a suitable control for Vancouver's DTES and we did not include data from this LHA as part of the aggregate analyses.
Focusing on the year before and after implementation in exposed and matched control LHAs captured 19,119 opioid-related overdose events (14,141 in exposed, and 4,978 in control LHAs) across the thirteen exposed LHAs (n~1,859,477 residents) and eleven matched LHAs (n~1,203,260 residents). Of these unique overdose events, 16,551 were attended by paramedics (12,277 in exposed, and 4,274 in control); 13,569 included a visit to the emergency department (10,704 in exposed, and 2,865 in control); 1,760 included a hospitalization (1,273 in exposed, and 487 in control); and 2,129 resulted in death (1,499 in exposed, and 630 in control).
After accounting for observed changes in levels and trends in matched controls, our analysis found no significant changes in monthly overdose mortality or hospitalization rates. However, we did observe significant decreases in trends of paramedic attended event rates (1.14 fewer events per 100,000 population per month; 95% CI: -1.99 to -0.28) after an initial (level) increase of 7.43 events per 100,000 population (95% CI: 1.27 to 13.58) immediately following implementation (Table 3 and Fig 1). In other words, paramedic events increased 30.1% after implementation of OPS/SCS but declined 3.0% per month thereafter, for an absolute difference of 6.19 fewer paramedic attended events per 100,000 (23.5% relative decrease) by twelve months post-implementation compared with expected. Similarly, we found emergency department visits declined (1.25 fewer events per 100,000 population; 95% CI: -1.95 to -0.55) after an increase of 3.93 events per 100,000 population (95% CI: -1.14 to 9.00) post-implementation, for a 22.8% initial increase followed by a 3.6% decline per month ultimately resulting in 11.11 fewer emergency department visits per 100,000 (39.0% relative decrease) than expected at twelve months post-implementation.

PLOS ONE
To test the robustness of the results we ran a series of sensitivity analyses. For previously excluded events with missing location information, we assigned the location of the first point of medical contact. Doing so captured 34,256 (93.7%) opioid-related overdose events, 26,220 (100.0%) of which were attended by paramedics, 23,266 (96.3%) included emergency department visits, 3217 (95.6%) included hospitalization, and 3528 (97.9%) resulted in death. Restricting analysis to twelve periods before and after implementation, we had 20,960 opioidrelated overdose events (+9.6%; 16,701 in early adopters-LHA that implemented OPS/SCS commencing December 2016 and matched controls, and 4,259 in late exposed-LHA that implemented their first site after December 2016 and matched controls); no change in the number of paramedic attended events (13,277 early and 3,324 late adopter matched pairs); 15,320 emergency department visits (+12.9%; 12,625 early and 1,658 late adopter matched pairs); 2,082 hospitalizations (+18.3%); and no change in deaths (1,669 early and 460 late adopter matched pairs). When we repeated our analysis using assigned event location, and again separating early adopters and their controls, and late adopters and their controls, the results remained unchanged (Table 3 and Fig 1).

Discussion
The results from our analyses were mixed. We found no statistically significant changes in opioid-related overdose mortality and hospitalization rates' levels or trends following OPS/ SCS implementation across matched LHA pairs. However, we did observe significant declines in monthly rates of paramedic attendance and emergency department visits following the introduction of OPS/SCSs. In the case of paramedic attended events, the overall decline in trend was preceded by an initial spike in level of events post-implementation. The observed level change was not surprising given original guidelines for OPSs developed with the BC CDC and The Portland Hotel Society recommended calling emergency health services and transferring care to paramedics for all overdose events reversed with naloxone [51][52][53]. With time, protocols for when to call emergency health services changed and allowed more discretion by the attending peer or staff (e.g. Fraser Health's 2018 Policy and Protocol Recommendations for Service Providers included text on determining a priori with clients when to call 9-1-1) [54].
Our results echo findings reported in other evaluations of OPS/SCSs locally and internationally. A 2007 Australian report by the National Centre for HIV Epidemiology and Clinical  Table 3. Interrupted time series for exposed and propensity score matched controls. All output is reported per 100,000 population. Research for the New South Wales Department of Health found "no statistically significant differences in the rates of decrease of opioid-related deaths between Kings Cross [where the SCS was implemented] and the rest of [New South Wales]" but did observe a significant decrease in ambulance attendances during the six years of follow-up [55]. Similarly, a report by British Columbia's Island Health Authority examining the early effects of recently implemented OPSs found some effect on local ambulance dispatches and emergency department visits but did not examine effects on mortality [43]. Finally, a mathematical model of the individual and combined impacts of BC's recently implemented or otherwise expanded harm reduction interventions (i.e. Take Home Naloxone program, OPSs, and opioid agonist treatment) estimated that 5% (95% credible interval (CRI) = 3-7%) of deaths were averted by OPSs alone, with 1.3 (95% CRI = 0.9-1.7) deaths averted per site per month between April 2016 and December 2017 [56]. The model's estimates may appear at odds with the reported 3476 overdose events reversed by OPS/SCSs in the first year post-implementation [13,56,57]. However, although overdoses account for substantial morbidity, not all events are fatal in the absence of medical or peer intervention (e.g. naloxone or oxygen); and some clients experience multiple overdose events [55,[58][59][60][61].

PLOS ONE
In a separate set of analyses examining the effects of newly implemented OPS/SCSs in Vancouver's DTES using traditional interrupted time series without a matched control, the positive effects observed therein may be a result of intervention-specific and community-level features unique to this LHA. For example, higher volume facilities, longer hours of operation, and no client restrictions (Table 1) may overcome barriers to access. Further, OPS/SCSs in the DTES may benefit from a legacy of strong, grass roots activism that reduce drug use-related stigma, thereby improving acceptability of this intervention locally. However, without an appropriate control, we cannot dismiss the possibility that regression to the mean or some uncontrolled historical bias explain the observed effect. Difference in trend between exposed and matched controls, preimplementation ( Over 69% of people who died from an illicit drug overdose in British Columbia between 2016 and 2017 consumed alone at the time of the fatal event [1]. This population is plausibly different from the target population of OPS/SCSs including Vancouver's Insite whose goal was to "attract the [individuals] who were at highest risk of health-related harms and those responsible for public order problems (e.g., public injection drug use)" [5,62]. Persistent stigma and police presence may impede the social acceptability and uptake of OPSs by at-risk and vulnerable populations [63]. Meanwhile, hours of operation, facility capacity (see Table 1), residence requirements and absence of safe inhalation rooms limit the effectiveness of these services [43,64]. Together, these factors may explain why there was limited uptake in some communities (e.g. 13 visits daily at Duncan site) [44] and no observed effect on overdose mortality rates.
Our study had several limitations. For starters, to ensure the privacy of individuals captured in the cohort, the smallest geographic unit available for analysis was the LHA. This restricted our ability to test the effects of OPSs at a more granular level, particular given the distancedecay effect observed in Marshall et al.'s 2011 seminal paper. This work demonstrated that Vancouver's Insite had a significant effect on overdose mortality within 500m of the supervised injection facility but minimal impact outside this radius [14]. More recently, the Island Health Authority's evaluation of OPSs in Campbell River, Courtenay, Cowichan Valley, and Port Alberni found their effects on paramedic dispatches waned outside a 1km radius [43]. While these studies suggest our unit of analysis was not fine enough to isolate truly localized effects on overdose mortality and hospitalization rates, they also indicate that it is unlikely population mobility across LHA boundaries diluted the effects of the new OPS/SCSs. Importantly, Marshall et al.'s work was specific to the very concentrated DTES community, aggregated over

PLOS ONE
five years of pre-and post-data, and did not explore changes in monthly rates. Conversely, our work covered geographically dispersed clients, evaluated the effects of the intervention in the year following implementation, and examined changes in outcomes by level (immediate) and trends (over time). To replicate Marshall et al.'s analysis across our exposed regions while protecting individual's privacy, accounting for temporal trends, and correcting for autoregression, we would need a much longer observation period with study intervals of quarters or years instead of months-negatively impacting the timeliness and relevance of the findings [28,65].
With respect to the propensity score matching approach applied, we were confined to using individual-level administrative data to create population-level summary health statistics and had to assume that these, combined with LHA-level 2015 overdose mortality rates, sufficiently predicted the probability of treatment assignment. This approach may not adequately capture the qualitative differences in local political leadership and culture which foster the successful implementation of harm reduction interventions, as well as differences in illicit drug market toxicity (all unobservable) which affects event counts. In other words, the LHAs that implemented SCSs may differ from the pool of potential controls such that residual confounding biases the effect estimates. This may be less of a problem for OPSs. These historically unsanctioned, nimble, grass-roots, and peer-initiated responses to the neglected needs of people who use illicit substances generally precede local political will and support; with several implemented during the observation period despite local back lash [66,67]. Given most sites began as OPSs, the concern that treated LHAs were sufficiently different from controls with respect to local support and thereby probability of treatment initiation, while non-negligible, should not diminish our findings. Meanwhile, BC Coroners Services' publicly available data of fentanyl-related and overall illicit overdose deaths, show consistent trends of overdose mortalities where fentanyl, W-18 and carfentanil were detected in treated vs. control regions [68][69][70]. The observed outcome rate trends pre-implementation (see β 3 estimates in Table 3 for paramedic attendance, hospitalization and mortality rates) met the parallel trends assumption critical in many econometric study designs corroborating the suitability of exposed-control pairings and further mitigating concerns that unobservable differences in political leadership, culture or illicit drug supply should bias our findings. Finally, segmented regression can compensate for imperfect controls (see β 3 estimates in Table 3 for ED visits) in ways that other study designs, such as difference-in-differences, cannot [71]. That said, local differences in the acceptability of the intervention may have facilitated or impeded use by at-risk individuals, affecting their observed effectiveness.
The Overdose Cohort does not include events that occurred at OPS/SCSs unless additional health services were sought. One interpretation of the results then may be that the effects of OPS/SCSs on paramedic attended events and ED visits is due to a substitution effect-these harm reduction sites reversed overdose events that would otherwise have utilized health services; but did not necessarily reduce the number of overdose events. Results from the BC CDC's publicly reported Overdose Response Indicators [72] combined with reports published by the health authorities support this interpretation of the decline in overdose events. For example, between December 2016 and November 2017 there were 235,466 visits and 1,253 overdoses reversed in the newly implemented OPS/SCSs (excluding Vancouver's Insite with 119,395 visits and 1432 overdoses reversed) [43,73].
Further, we were unable to directly account for other factors that may have influenced the rates of overdose mortality and health services used. Our data did not allow us to control for the time-varying distribution of Take Home Naloxone (THN) kits, potency of drugs on the illicit market, or changes in access to opioid agonist treatment. However, publicly available Overdose Response Indicators produced by the BC CDC and the Ministry of Mental Health and Addictions suggest no corresponding changes to the number of THN kits distributed, or clients dispensed opioid agonist treatment that better account for the changes in health service use and overdose mortality rates observed post-implementation [57]. Meanwhile, BC Coroners Services' data suggest fentanyl contamination trends (both relative and absolute) were more aligned between propensity score matched treated-control pairs than between treated and neighboring geographic units, helping rule out that differences in street drug toxicity explain the lack of effect observed for hospitalizations and mortalities post-implementation [68][69][70].
Lastly, we were unable to find a suitable control for Vancouver's DTES given the LHA's unique population profile and needs. At the time of intervention implementation, the DTES' outcome rates were at least one order of magnitude greater than those observed in any other LHA (deaths = 50.86; hospitalizations = 66.07; emergency department visits = 578.33; and paramedic attendance = 664.24 per 100,000 population) with unparalleled growth in rates. Similarly, creating a synthetic control using Abadie et al's. method was infeasible, since the outcomes fell outside the convex hull of available controls [74]. As such, we interpret these findings separately and with caution.
Despite these limitations, this is the largest evaluation of the population-level effects of OPS/SCSs across a variety of settings, available to date. We explored the impact of the intervention on opioid-related overdose mortality and health service use using the most comprehensive data set available, across a variety of settings (i.e. rural, remote and urban communities), and over time. Using matched controls and applying a multiple baseline approach allowed us to better account for potential historical bias [75]. Our analyses included approximately two thirds of BC's population, and over 80% of opioid-related overdose events that occurred in the year before and after implementation of OPS/SCSs in matched pair LHAs. Overall, our analyses showed new OPS/SCSs reduced health service use substantially but we were unable to detect an effect on overdose mortality rate a year after implementation. To our knowledge, this is the first study to comprehensively explore the effects of OPS/SCSs in contexts where the target population is not concentrated in a small geographic area.

Conclusion
The overdose epidemic in BC is unique to that observed elsewhere [76]. In the United States, the epidemic is described as a triple wave of overdose deaths starting with prescription opioids, followed by heroin, and more recently, fentanyl [77,78]. Other regions in Canada demonstrate a similar epidemiological transition from prescription opioids to illicit substances. Despite the early epidemiological differences, regions across North America are now contending with similar illicit opioid overdose epidemics.
In Canada, OPSs and SCSs are increasingly employed as a strategy to reduce overdoserelated morbidities and mortality. Our study shows that this harm reduction intervention reduces paramedic and emergency department use; and the results may be of particular interest to jurisdictions considering mobile units or implementation in less dense communities (e.g. the Appalachian region of the United States). Additional research is needed to understand the effects of operating hours, service volume, residence requirements and police presence. Studies quantifying the impacts of SCSs as low barrier access to other health services for marginalized populations may identify additional benefits.
OPSs including dates of operation and access restrictions. We would also like to thank the anonymous reviewers at the BC CDC for their thoughtful comments on the original draft.
Disclaimer: All inferences, opinions, and conclusions drawn in this paper are those of the author, and do not necessarily reflect the opinions or policies of data stewards including the Ministry of Health, Ministry of Mental Health and Addictions, BC Centre for Disease Control, BC Coroner's Service, or BC Emergency Health Services.