Research Article 106 Ved ABSTRACT ABSTRACT Objective: Little is known about the population-level impact of recreational cannabis legalization on trends in opioid-related mortality. Increased access to cannabis due to legalization has been hypothesized to reduce opioid-related deaths because of the potential opioid-sparing effects of cannabis. The objective of this study was to examine the relations between national retail sales of recreational (non-medical) cannabis and opioid overdose deaths in the 5 years following legalization in Canada. Method: Using time-series data, we applied Granger causality methods to evaluate the association between trends in legal recreational cannabis sales and opioid-related deaths over time. Both sales and opioid deaths grew over time, with the latter exhibiting significant increases following the onset of the COVID-19 pandemic. Results: We found no support for the hypothesis that increasing post-legalization sales Granger caused changes in opioid-related deaths in British Columbia, Ontario, or at the national level. Conclusions: These findings suggest that increases in legal recreational cannabis sales following legalization were not meaningfully associated with changes in opioid-related mortality. Further examination with longer follow-up periods will be needed as the legal cannabis market becomes more entrenched in Canada, but these findings converge with previous work suggesting legalization is not related to opioid overdose mortality and further undermine that hypothesized link as a basis for legalization in other jurisdictions. Key words: = cannabis; legalization; opioid; mortality; behavioural economics As more jurisdictions legalize recreational (non-medical) cannabis, there has been obvious interest regarding how legalization impacts cannabis use (Kelsall, 2018), but also its collateral impact on use of other substances including opioids (Mathur & Ruhm, 2023). In Canada, more than 30,000 people have died from opioid-related overdoses since 2016, which is more than all other accidental death causes combined (Fischer et al., 2023). The opioid overdose epidemic has been recognized as a public health crisis in Canada (Public Health Agency of Canada, 2023), and, André J. McDonald1-3*, Alysha Cooper1-3*, Amanda Doggett1-3, Kyla Belisario1-3, James MacKillop1-3 1Peter Boris Centre for Addictions Research, St. Joseph's Healthcare Hamilton 2Department of Psychiatry and Behavioural Neurosciences, McMaster University 3Michael G. DeGroote Centre for Medicinal Cannabis Research, McMaster University *Co-first (equal contributions) Cannabis 2024 © Author(s) 2024 researchmj.org 10.26828/cannabis/2024/000261 Volume 7, Special Issue 3 Legal Recreational Cannabis Sales and Opioid-Related Mortality in the 5 Years Following Cannabis Legalization in Canada: A Granger Causality Analysis Corresponding Author: André Mcdonald, PhD, McMaster University. 100 West 5th Street, Hamilton, Ontario, Canada, L8P 3P2. Email: mcdona36@mcmaster.ca Cannabis, A Publication of the Research Society on Marijuana 107 based on potential substitutability relationships —so-called “opioid-sparing effects” where cannabis reduces the need for opioids without loss of analgesic efficacy (Nielsen et al., 2017)—there is interest in whether regulatory changes in cannabis access may have a favorable impact (Mathur & Ruhm, 2023). Previous research—mostly from the U.S— examining the relationship between cannabis legalization and opioid-related outcomes is mixed. Some studies have found a substitution effect between cannabis legalization and opioid-related mortality (e.g., Chan et al., 2019; Powell et al., 2018), with one major study suggesting that states with medical cannabis laws experienced slower increases in opioid-related deaths from 1999 to 2010 (Bachhuber et al., 2014). However, a subsequent study extended the observation period and found that this trend reversed direction, suggesting that the association may be spurious (Shover et al., 2019). Another large study also found no significant association between state cannabis laws (recreational and medical) and opioid-related mortality (Nguyen et al., 2024). Meanwhile, a recent study found that legalized cannabis in the U.S. was associated with increased opioid mortality rates (Mathur & Ruhm, 2023), though the authors acknowledge that the association may have been driven by the emergence and proliferation of fentanyl, a particularly potent opioid in the drug supply. Few national jurisdictions have legalized cannabis for recreational purposes, and those that have, such as Canada in 2018, did so relatively recently, meaning little is known about the population-level impact of legal cannabis sales on trends in opioid-related mortality. This represents a significant research gap that could help inform the continuing policy response in Canada and provide guidance to other jurisdictions considering cannabis policy reform. National retail sales data provide an opportunity to objectively study trends and dynamics between population-level purchasing patterns of legal commodities, including cannabis (MacKillop et al., 2021), yet few studies have leveraged these data. Therefore, the objective of this study was to examine the relationship between national retail sales of legal recreational cannabis in Canada from its inception in relation to opioid deaths. Using time-series data, we used Granger causality—i.e., one variable forecasting or “Granger causing” another variable—to evaluate the interrelationships between trends over time. This approach evaluates whether, in two parallel time series, indicators earlier in one time series systematically predict an indicator at later points in a second time series (Granger, 1969). METHODS Data Sources This study used national level data as well as provincial level data from Ontario and British Columbia, the two provinces with the greatest number of opioid-related deaths (Public Health Agency of Canada); monthly data were not available from other provinces. At the national level, data included cannabis sales ($ CAD) and number of opioid-related deaths from October 2018 to June 2023. Given that the national data were measured quarterly for opioid-related deaths, the total monthly cannabis sales were aggregated into quarterly sales. Given the small number of data points of the quarterly time series, only a single time lag was considered when fitting the model. In BC, data were available from October 2018 to September 2023, and in Ontario from October 2018 to March 2023. For the sales data, we adjusted for inflation and price changes (cannabis prices declined markedly over the study period) by dividing the series by the consumer price index series and multiplying by 100 to revert to the original scale. The number of opioid-related deaths were divided by the estimated jurisdictional population for the given time frequency and reported as number of deaths per 1 million. National and provincial sales data were collected from Table 20-10-0056-01 in Statistics Canada (2023a), consumer price index (CPI) data were collected from Table 18-10-0004-10 (Statistics Canada, 2023b), the number of opioid-related deaths in BC was collected from the British Columbia Centre for Disease Control (BCCDC, 2024), the quarterly estimates of the population in BC were collected from Table 17-10-0009-01 in Statistics Canada (2023c), and both the number of opioid-related deaths in Ontario along with the estimated population were collected from the Interactive Opioid Tool available through Public Health Ontario (2023). We note that for BC, monthly data were only available for unregulated drug deaths at the aggregate level, which included Legal Cannabis Sales and Opioid Deaths in Canada 108 some deaths that were related to drugs other than opioids; however, on average, over 85% of unregulated drug deaths involved fentanyl, and 24% involved other opioids for the years studied (BCCDC, 2024). National level opioid-related death data were obtained from the Public Health Agency of Canada (2023). All data used in this study are publicly available (see eTable 1 in supplementary materials for all time-series data used for this study). Vector Auroregression and Granger Causality Vector autoregression (VAR) is a series of regression equations where each variable at time t is predicted by time lags of itself and the other variable(s) in the series. The number of time lags, p, included in the model can be chosen using information criteria and/or final prediction error (FPE). In this study, we chose the number of time lags based on the majority vote between the Akaike Information Criterion (AIC; Akaike, 1969), Schwarz Criterion (SC; Schwarz, 1978), Hannan-Quinn criterion (HQ; Hannan & Quinn, 1979), and the FPE. After fitting the VAR models, the Granger causality tests were conducted to determine whether one variable could be used to forecast another variable. Prior to conducting Granger causality, it is important to ensure that all model assumptions of the VAR model are met. The assumptions of the VAR model include non-serially correlated residuals, homoscedastic residuals, normally distributed residuals, and stability of parameter estimates. The asymptotic Portmanteau test was conducted to evaluate whether model residuals were serially correlated. It is recommended to conduct the Portmanteau test at lags considerably larger than p, but not too large such that the power of the test is reduced (Kilian & Lütkepohl, 2017). Therefore, we considered lags from p+2 to f (i.e., the frequency of the series) in the Portmanteau test to ensure the results were robust across lags. In addition, the ARCH-LM test (Engle, 1982) was conducted to evaluate whether or not there was evidence of heteroskedasticity of model residuals from 1 to f lags. Finally, the Jarque-Bera (JB) test (Jarque & Bera, 1987) was used to determine whether there was evidence of deviation from the normal distribution in the model residuals, and the OLS-CUSUM test (Ploberger & Krämer, 1992) was applied for evaluating the stability of model parameters. VAR models can also include other variables which may change over time. Therefore, to ensure that the onset of COVID-19 did not change our results, all analyses were re-conducted with a dummy variable equal to 1 during and after March 2020 and 0 otherwise. The BC government implemented a safer supply program in March 2020 (Nguyen et al., 2024); therefore, for BC, the dummy variable represented both COVID-19 and the safer supply program. No conclusions changed when fitting the monthly models with the dummy variable; therefore, this variable was not included in the reported analyses. All models were fit with the vars package in R (Pfaff, 2008a; R Core Team, 2024). Pre-Processing Data A requirement of VAR models is that the series are stationary (i.e., its statistical properties such as mean and variance do not depend on time). Therefore, the series were evaluated for seasonal patterns and trends prior to fitting the models. In terms of the trend, if the trend was deterministic (i.e., predictable and described as a function of time), then it would be accounted for within the VAR models, whereas if the trend was stochastic (i.e., unpredictable movements of the series), then it would be removed from the series prior to fitting the model. In addition, if seasonality was present in the series, it would be removed prior to model fitting as well, to ensure the series in the models were white noise. A combination of the Kruskal-Wallis (KW) test (Kruskal & Wallis, 1952) and the modified QS test (Maravall, 2011) was used to detect seasonality, such that the data were deemed to be seasonal if both tests were significant at the 5% level of significance. The tests were conducted within the seastests package in R (Ollech, 2021). If the time series had a significant seasonal component, it was estimated with “Seasonal and Trend decomposition using Loess” (STL; Cleveland et al., 1990) and subtracted from the series. The Kwiatkowski–Phillips–Schmidt–Shin (KPSS) test (Kwiatkowski et al., 1992) and the Augmented Dickey-Fuller (ADF) test (Dickey & Fuller, 1979) were used for testing whether there was a unit root in the series. Each series was first evaluated with the KPSS and ADF tests for whether they were stationary around a trend. If Cannabis, A Publication of the Research Society on Marijuana 109 the series were found to be stationary around a trend, they were evaluated for whether they are stationary around a constant. Finally, if the series were stationary around a constant, then a final ADF test was conducted for whether they were stationary without a trend or constant. If a series was found to be stationary only after accounting for a trend/constant, then a trend/constant was included in the respective VAR model. The KPSS tests were conducted with the kpss.test function within the tseries package in R (Trapletti et al., 2007), and the ADF tests were conducted with the ur.df function in the urca package in R (Pfaff, 2008b). If the series were non-stationary at the 5% level of significance according to either the KPSS test or the ADF test (i.e., p < 0.05 for KPSS or p > 0.05 for ADF), they were first differenced until stationarity was achieved. If a time series appeared to have a structural break (i.e., sudden change in data generating process), then the structural break point would be estimated via the breakpoints function within the strucchange package in R (Zeileis, Kleiber, Krämer, & Hornik, 2003; Zeileis, Leisch, Hornik, & Kleiber, 2002). In the presence of a structural break, the Zivot- Andrews test was conducted in addition to the KPSS and ADF tests. The Zivot-Andrews test accounts for an unknown single structural break when evaluating the stationarity of a time series. A rejection of the null hypothesis in the Zivot- Andrews test implies the series is trend stationary with a break in either the a) intercept, b) trend, or c) the intercept and the trend, depending on the test conducted. For the purpose of our analyses, we conducted the Zivot-Andrews test accounting for a potential change in the intercept and trend. Toda-Yamamoto Approach Granger causality in differences can potentially be problematic if there is co- integration between the variables in the model. To ensure robustness of our results, we additionally applied the Toda-Yamamoto (1995) approach for testing Granger causality. The Toda-Yamamoto approach is less sensitive to incorrect identification of the order of integration or co- integration of the series. To ensure the proper distribution of the test statistic, this approach involves augmenting a VAR model with additional lags up to the “max integration order”. The max integration order is defined as the order of integration required to make both series stationary. For example, if we initially fit a VAR model with p = 2 lags on the first-differenced data (i.e., integration order of 1), then we would fit the VAR model with p + 1 = 3 lags on the non- stationary data for the Toda-Yamamoto approach. However, the test of causality is only evaluated for the first p lags of the model and the augmented lags are not included in the test. In addition, the Fourier Toda Yamamoto method (Nazlioglu, Gormus & Soytas, 2016) was applied to confirm reliability of our findings when treating the structural break from COVID-19 as unknown (see eResults 1 in supplemental materials for further details). Interrupted Time Series Analysis We additionally conducted interrupted time series analysis via autoregressive integrated moving average (ARIMA) modelling to evaluate the impact of legalization on opioid deaths in BC and Ontario (there were not enough quarterly data at the national level). For both Ontario and BC, we used monthly opioid-related deaths spanning from January 2013 (pre-pandemic) to October 2023. These analyses included a variable that allows for a gradual change in opioid-related mortality following the onset of the COVID-19 pandemic (and the safer supply program in BC). We conducted interrupted time series analysis with ARIMA models to account for any auto-correlation in the series and/or seasonality when required. We included two variables for the potential impact of legalization implementation: a variable for a step/level change immediately following legalization and a variable for a change in the outcome over time following legalization. The breakpoint of interest for both analyses was October 2018, the date of legalization implementation. However, we also looked to see whether the passage of the Cannabis Act in November 2017 (Cannabis Act, 2018; Library of Parliament, 2018) had an effect or improved model fit. We also examined naloxone policy changes in Ontario (naloxone provided for free in pharmacies June 2016) and BC (naloxone unscheduled in Sept 2016) as intervention points. As passage of the Cannabis Act and naloxone policy changes did not have an effect on opioid- related mortality in either province, we excluded these intervention points from the final models. Legal Cannabis Sales and Opioid Deaths in Canada 110 Cannabis was not legalized until October 2018 in Canada (Cannabis Act, 2018; Library of Parliament, 2018), and therefore, the data in our analyses were restricted in size, thereby limiting the methods that could be used for seasonal adjustment. We chose to use the STL approach for seasonal adjustment, as it preserves the sample size of the data and can be estimated with a few years worth of data. However, the moving average filters used in STL may reduce the power of unit root tests (Ghysels & Perron, 1993), and in doing so may impact the perceived order of integration. To overcome this, we could have adjusted for seasonality inside the model, as opposed to outside the model, through the inclusion of eleven seasonal dummy variables in the model. In considering the ratio of sample size to number of coefficients to estimate, we opted to seasonally adjust the data outside the model via STL. However, a sensitivity analysis confirmed that all conclusions remained the same to those reported below when adjusting for seasonal patterns inside the model as opposed to outside the model. RESULTS Visualizing Changes Over Time The time series of cannabis sales and national and provincial opioid-related mortality rates from October 2018 up to September 2023 are presented in Figure 1. A steady upward trend can be observed for legal cannabis sales since legalization, increasing by more than 10 times from $42 million per month in October 2018 to $444 million per month five years later in September 2023. In terms of the opioid mortality rates, both nationally and provincially, the trend and intercept appear to change with the onset of COVID-19 in March 2020, such that the mean mortality rate increased and the trend becomes more positive with the onset of the pandemic. Interrupted time series analyses were conducted to examine whether the structural break in these time series were significant. Nationwide, there was a significant immediate increase in the mortality rate at the onset of COVID-19 (p < .001) but not a significant change in the trend (p = .466). The same findings were found for Ontario, where there was a significant immediate increase in the mortality rate (p < .001), but not a change in the trend (p = .142). On the other hand, both a significant increase in the level (p < .001) and trend (p < .001) of the mortality rate were found in BC following the onset of COVID-19. This structural break in the data was accounted for in the pre-processing and analysis stages of our methods. Figure 1A. Quarterly Legal Cannabis Sales ($) and Opioid-Related Mortality Rate in Canada Between October 2018 and September 2023 Cannabis, A Publication of the Research Society on Marijuana 111 Figure 1B. Monthly Legal Cannabis Sales ($) Between October 2018 and September 2023 and Opioid-Related Mortality Rate in Ontario Between October 2018 and March 2023 Figure 1C. Monthly Legal Cannabis Sales ($) and Opioid-Related Mortality Rate in British Columbia Between October 2018 and September 2023 Pre-Processing Data Table 1 presents the results of the KW and QS tests for seasonality of each time series. Due to an apparent structural break in the opioid-related mortality rates, these series were first split in two (pre- and post- COVID-19 onset) for assessment of seasonal patterns to ensure the seasonal patterns did not change with the structural break. Significant seasonality was detected in the BC cannabis sales series. The seasonality of the series was estimated with STL decomposition (STL decompositions are presented in eFigure 1 in supplemental materials) and subtracted from the series. After seasonal adjustment, the tests for seasonality did not find evidence of seasonal patterns remaining in the data. All series were first differenced to achieve stationarity (see Table Legal Cannabis Sales and Opioid Deaths in Canada 112 2). The Zivot-Andrews test suggested that first- differencing was necessary for the national opioid mortality rate data when accounting for an unknown structural break (p > .05) in addition to being necessary for both the BC and Ontario opioid mortality rate data (p > .05 for each). After first-differencing, a structural break was no longer detected in the national, ON, or BC time series. Table 1. Seasonality Tests for Unadjusted and Seasonally Adjusted Time Series QS KW Seasonal? CA cannabis sales Unadjusted 4.48 9.23** No Seasonally adjusted N/A N/A N/A CA opioid mortality rate Unadjusted 0.00 3.40 No Seasonally adjusted N/A N/A N/A BC cannabis sales Unadjusted 8.80** 21.46** Yes Seasonally adjusted 0.00 5.48 No BC opioid mortality rate Unadjusted 1.80 17.24 No Seasonally adjusted N/A N/A N/A ON cannabis sales Unadjusted 1.01 17.84 No Seasonally adjusted N/A N/A N/A ON opioid mortality rate Unadjusted 6.44** 19.27* No Seasonally adjusted N/A N/A N/A Note. *, **, and *** denote significance at the 10%, 5%, and 1% level, respectively. N/A = not applicable; QS = modified QS test (Maravall, 2011); KW = Kruskal-Wallis test (Kruskal & Wallis, 1952). Table 2. KPSS and ADF Unit Root Test Statistics for Seasonally Adjusted Cannabis Sales ($) in Canada from October 2018-September 2023, and Opioid Mortality Rate per Million in British Columbia from October 2018-September 2023, Ontario from October 2018-June 2022, and Canada from October 2018-June 2023 KPSS ADF CA cannabis sales Level 0.143 (t) * -2.214 (t) First difference 0.102 (c) -4.317 (c) *** CA opioid mortality rate Level 0.088 (t) -2.342 (t) First difference 0.061 (c) -3.050 (n) *** BC cannabis sales Level 0.248 (t) *** -2.199 (t) First difference 0.075 (c) -7.939 (c) *** BC opioid mortality rate Level 0.188 (t) ** -1.885 (t) First difference 0.102 (c) -5.737 (n) *** ON cannabis sales Level 0.159 (t) ** -2.478 (t) First difference 0.102 (c) -4.264 (n) *** ON opioid mortality rate Level 0.197 (t) ** -3.376 (t) * First difference 0.061 (c) -5.844 (n) *** Note. The letter in parentheses is for whether the KPSS and ADF tests are testing for stationarity with trend (t), stationary with constant (c), or stationary without trend/constant (n). *, **, and *** denote significance at the 10%, 5%, and 1% level, respectively. Vector Autoregression Models and Granger Causality of Cannabis Sales and National Opioid Mortality Rates. A VAR model with a single time lag was fit to the first-differenced, quarterly cannabis sales and national opioid-related mortality rate per million. Cannabis, A Publication of the Research Society on Marijuana 113 The Portmanteau test did not find evidence to suggest serial correlation among model residuals, (p > .05 for lags 3 and 4) and the ARCH-LM test verified that the residuals were homoscedastic up to 4 lags, respectively (p > .05 for all tests). The residuals did not deviate from normality, 𝜒2(4) = 6.310, p = 0.177, and the parameters of the model were stable over time. Table 3 presents the summary of results for the VAR model. Table 3. Results of VAR Models for Cannabis Sales ($) and Rate of Death per 1 Million due to Opioid Overdose in Ontario, British Columbia, and Canada with Rate of Death due to Opioids as the Outcome Opioid deaths Jurisdiction Time series Predictor �̂� (SE) t p-value Canada Quarterly Oct 2018 to Jun 2023 Constant -1.208 (2.107) -0.573 .576 Cannabis(t-1) 0.193 (0.242) 0.798 .438 Opioid(t-1) 2.768e-05 (1.954e-05) 1.417 .178 Ontario Monthly Oct 2018 to Mar 2023 Constant 0.163 (0.248) 0.656 .515 Cannabis(t-1) -1.394e-05 (2.065e-05) -0.675 .503 Opioid(t-1) -0.152 (0.141) -1.079 .286 British Columbia Monthly Oct 2018 to Sept 2023 Constant 0.454 (0.617) 0.730 .465 Cannabis(t-1) -1.777e-04 (1.755e-04) -1.010 .317 Opioid(t-1) -0.327 (0.127) -2.568 .013 First-differenced quarterly cannabis sales at time t-1 were not associated with the national opioid-related mortality rate at time t. The Granger causality test did not find evidence to suggest that cannabis sales in Canada Granger caused changes in the national opioid-related mortality rate, F(1,28) = 2.007, p = .168. Vector Autoregression Models and Granger Causality of Cannabis Sales and Provincial Opioid Mortality Rates. Two separate VAR models were fit: one for the relation between first-differenced cannabis sales and opioid mortality rate per million in BC, and the other for the relation between first-differenced cannabis sales and the opioid mortality rate per million in ON. A constant was included in the ON model to fulfill the stability assumption, and it was also incorporated into the BC model, given that BC cannabis sales exhibited stationarity around a constant. One lag was selected for each model. The Portmanteau test and the ARCH-LM test verified that the residuals of each model were not serially correlated between 3 and 12 lags. However, the ON model had significant heteroskedasticity of model residuals at lag 1 (p < .01). To account for the heteroskedasticity of residuals, heteroskedastic-consistent standard errors were used when evaluating Granger causality. The residuals did not deviate from normality for the BC model, 𝜒2(4) = 1.671, p = .796, or the ON model, 𝜒2(4) = 3.746, p = .442, and the parameters of each model were stable over time. Cannabis sales at time t-1 were not significantly associated with deaths caused by opioid overdose in BC, nor opioid overdose in ON at time t. The Granger causality test for the BC model did not find evidence of cannabis sales Granger causing the opioid mortality rate in BC, F(1,110) = 1.019, p = .315, nor did the Granger causality test for the ON model, F(1,98) = 0.753, p = .388. The Toda-Yamamoto method did not find evidence that cannabis sales Granger caused changes in the opioid mortality rate in BC, 𝜒2(4) = 4.839, p = .304, or ON, 𝜒2(1) =1.196, p = .274. The Toda-Yamamoto method could not be conducted for the opioid-related mortality rate at the national level due to the small sample size. Interrupted Time Series Analysis via Autoregressive Integrated Moving Average (ARIMA) Modelling Like the Granger causality analyses, we did not find a significant association between legalization and opioid deaths in either Ontario or Legal Cannabis Sales and Opioid Deaths in Canada 114 BC. We found a significant outlier in Ontario in October 2017, when the mortality rate appeared to start decreasing. An ARIMA(1,1,0) model with drift was selected, and the data was found to be non-seasonal. There was not a significant level change in mortality rate following legalization implementation (estimated level change: 0.060, SE = 1.185, Z = 0.051, p = .960), nor a significant change over time (estimated slope change: -0.130, SE = 0.177, Z = -0.732, p = .464). Figure 2A below shows the observed opioid-related mortality rate versus the counterfactual rate, which is the prediction based on what would have happened had legalization implementation not occurred. COVID-19 did have a significant impact on the outcome, with a gradual increase to a new level of mortality rate. Figure 2A. Observed Opioid-Related Mortality Rate Versus the Counterfactual Rate in Ontario In BC, we found similar results. We found a significant outlier in November 2016, where it appears that the mortality rate started to increase. An ARIMA (0,1,1) model was selected and no seasonality was detected. There was not a significant level change in mortality rate following legalization implementation (estimated level change: -3.101, SE = 2.479, Z = -1.251, p = .211) and no evidence to suggest that the rate changed over time after this point (estimated slope change: -0.182, SE = 0.175, Z = -1.040, p = .298). Figure 2B shows the observed opioid- related mortality rate versus the counterfactual rate had legalization implementation not occurred [counterfactual based on ARIMA(0,1,0) model]. We note that the confidence interval was wide for BC, likely because the time series was noisy pre- legalization and did not follow a stable pattern. Moreover, the forecast is based on the ARIMA(0,1,0), which forecasts just a single value, and the error accumulates linearly because of this over time. As shown, COVID-19 led to a significant, gradual increase to a new level of mortality rate. Cannabis, A Publication of the Research Society on Marijuana 115 Figure 2B. Observed Opioid-Related Mortality Rate Versus the Counterfactual Rate in British Columbia DISCUSSION This study found no evidence that increasing access to cannabis following legalization, as measured by increasing cannabis sales, Granger caused changes in opioid-related deaths over time. When examining each trend in isolation, cannabis sales had a pronounced and consistent increase, while the rate of opioid-related deaths increased sharply at the beginning of the COVID-19 pandemic and did not return to pre-pandemic levels. This qualitative step-like increase has been recognized as a major public health problem by the Government of Canada (2023). Potential reasons for the increase in opioid-related harms during the pandemic include border and travel restrictions contributing to unsafe and interrupted drug supply, physical distancing and isolation increasing the likelihood of solitary consumption, increased stress and worsening mental health, and reduced capacity in harm reduction services, among others (Government of Canada, 2023). The findings of this study add to a relatively scarce literature on the relationship between various forms of cannabis legalization and opioid- related mortality. Previous research is mixed, with some studies finding that cannabis legalization led to a reduction in opioid-related mortality (Chan et al., 2020; Powell et al., 2018; Bachhuber et al., 2014), and others an increase (Mathur & Ruhm, 2023) or no effect at all (Nguyen et al., 2024; Shover et al., 2019). Our study is largely consistent with recent studies examining recreational cannabis legalization suggesting that legalization is not associated with meaningful changes in opioid-related mortality (Nguyen et al., 2024; Mathur & Ruhm, 2023; Shover et al., 2019). As such, arguments for cannabis legalization based on the potential benefits for opioid overdoses are not rooted in evidence, nor are arguments against cannabis legalization on the grounds that it exacerbates the opioid crisis. This is one of the first studies to examine the impact of cannabis legalization on opioid-related deaths in a jurisdiction outside the US. The policy Legal Cannabis Sales and Opioid Deaths in Canada 116 landscape is different in Canada compared to the US, as cannabis is federally legal in Canada, whereas only individual states have legalized cannabis for medical and/or non-medical purposes in the US. Cannabis use was already highly prevalent in Canada prior to recreational cannabis legalization, with approximately 47% of people 15 years of age and over reporting lifetime cannabis use and 15% using in any given year (Statistics Canada, 2018). Meanwhile, problematic use of opioids was reported by only 0.4% of Canadians in the year prior to cannabis legalization (Statistics Canada, 2018); therefore, the potential for legalization to effect changes in opioid-related mortality was limited to a small subgroup of Canadians. This study had several strengths. We used national and provincial data that are inherently representative and generalizable to Canadians using legal cannabis. We used a time-series design, which is superior to cross-sectional designs used by most previous studies because it allowed for the establishment of temporality and followed best practices for Granger causality analysis. Moreover, previous research is largely based on self-reported data, whereas we analyzed sales data, which are objective and do not suffer from social desirability bias, recall bias, or underreporting. The study also had limitations. Most importantly, this study did not have a non- legalization control jurisdiction to compare to. It did not have pre-legalization data on purchasing patterns of cannabis because cannabis was still illegal, which prevented a pre-post comparison that would have allowed us to observe whether pre-legalization trends in our indicators changed or stayed the same following cannabis legalization. Some of our analyses were likely underpowered due to the limited availability of post-legalization time points, especially for opioid- related deaths. We were unable to control for the impact of naloxone distribution on opioid mortality, as no monthly or quarterly data were available over the study period, though we conducted interrupted time series analysis, which found no association between naloxone policy changes and opioid mortality. Legal cannabis sales data did not measure quantity used, although we adjusted for the consumer price index of cannabis, which declined markedly over time. Legal sales are also an imperfect measure of people’s access to cannabis, as many people acquired cannabis from illegal sources before and after legalization. The data we used were ecological; further studies are needed that examine the impact of legalization on subgroups of interest, such as by age, sex, and gender. Finally, the COVID-19 pandemic was a major event that impacted substance use and opioid- related deaths at the population level. While we adjusted for the pandemic using a dummy variable, it is impossible to know what the counterfactual would have been had the pandemic not occurred. To our knowledge, this is the first study to provide statistical evidence of significant increases in opioid deaths following the onset of the pandemic, as well as a steeper escalation in BC. Conclusions This study found that increases in legal recreational cannabis sales following legalization were not meaningfully associated with changes in opioid-related deaths. Although it has been speculated that cannabis legalization might support a greater cannabis-opioid substitution effect, thereby reducing harms related to opioids, this study extends previous findings suggesting that this does not happen in practice. Further studies with longer follow-up periods will be needed to confirm whether these findings hold as the legal cannabis market becomes more entrenched in Canada. REFERENCES Bachhuber, M. A., Saloner, B., Cunningham, C. O., & Barry, C. L. (2014). Medical cannabis laws and opioid analgesic overdose mortality in the United States, 1999-2010. JAMA Internal Medicine, 174(10), 1668–1673. https://doi.org/10.1001/jamainternmed.2014.4 005 British Columbia Centre for Disease Control. (2024). Unregulated Drug Poisoning Emergency Dashboard. https://public.tableau.com/app/profile/bccdc/vi z/UnregulatedDrugPoisoningEmergencyDash board/Introduction Cannabis Act, S.C. 2018, c. 16. (2018). https://laws-lois.justice.gc.ca/eng/acts/C-24.5/ Chan, N. W., Burkhardt, J. & Flyr, M. (2020). The effects of recreational marijuana legalization Cannabis, A Publication of the Research Society on Marijuana 117 and dispensing on opioid mortality. Economic Inquiry, 58, 589- 606. https://doi.org/10.1111/ecin.12819 Cleveland, R. B., Cleveland, W. S., McRae, J. E., & Terpenning, I. (1990). STL: A seasonal- trend decomposition. Journal of Official Statistics, 6(1), 3-73. Dickey, D. A., & Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association, 74(366a), 427-431. https://doi.org/10.2307/2286348 Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987-1007. https://doi.org/10.2307/1912773 Ghysels, E., & Perron, P. (1993). The effect of seasonal adjustment filters on tests for a unit root. Journal of Econometrics, 55(1-2), 57-98. https://doi.org/10.1016/0304-4076(93)90004-O Government of Canada. (2023). Data on cannabis for medical purposes. https://www.canada.ca/en/health- canada/services/drugs- medication/cannabis/research-data/medical- purpose.html Granger, C. W. J. (1969). Investigating Causal Relations by Econometric Models and Cross- spectral Methods. Econometrica, 37(3), 424– 438. https://doi.org/10.2307/1912791 Jarque, C. M., & Bera, A. K. (1987). A test for normality of observations and regression residuals. International Statistical Review/Revue Internationale de Statistique, 55(2), 163-172. https://doi.org/10.2307/1403192 Kilian, L., & Lütkepohl, H. (2017). Structural vector autoregressive analysis. Cambridge University Press. https://doi.org/10.1017/9781108164818 Kruskal, W. H., & Wallis, W. A. (1952). Use of ranks in one-criterion variance analysis. Journal of the American Statistical Association, 47(260), 583-621. https://doi.org/10.1080/01621459.1952.104834 41 Kwiatkowski, D., Phillips, P. C., Schmidt, P., & Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root: How sure are we that economic time series have a unit root?. Journal of Econometrics, 54(1-3), 159-178. https://doi.org/10.1016/0304-4076(92)90104-Y Library of Parliament. (2018). Legislative Summary of Bill C-45: An Act respecting cannabis and to amend the Controlled Drugs and Substances Act, the Criminal Code and other Acts. https://lop.parl.ca/sites/PublicWebsite/default/ en_CA/ResearchPublications/LegislativeSum maries/421C45E MacKillop, J., Cooper, A., & Costello, J. (2021). National Retail Sales of Alcohol and Cannabis During the COVID-19 Pandemic in Canada. JAMA Network Open, 4(11), e2133076. https://doi.org/10.1001/jamanetworkopen.202 1.33076 Maravall, A. (2011). Seasonality Tests and Automatic Model Identification in TRAMO- SEATS. Bank of Spain: Madrid, Spain. Mathur, N. K., & Ruhm, C. J. (2023). Marijuana legalization and opioid deaths. Journal of Health Economics, 88, 102728–102728. https://doi.org/10.1016/j.jhealeco.2023.102728 Nazlioglu, S., Gormus, N. A., & Soytas, U. (2016). Oil prices and real estate investment trusts (REITs): Gradual-shift causality and volatility transmission analysis. Energy Economics, 60, 168-175. https://doi.org/10.1016/j.eneco.2016.09.009 Nguyen, H. V., McGinty, E. E., Mital, S., & Alexander, G. C. (2024). Recreational and Medical Cannabis Legalization and Opioid Prescriptions and Mortality. JAMA Health Forum, 5(1), e234897. https://doi.org/10.1001/jamahealthforum.2023 .4897 Nielsen, S., Sabioni, P., Trigo, J. M., Ware, M. A., Betz-Stablein, B. D., Murnion, B., Lintzeris, N., Khor, K. E., Farrell, M., Smith, A., & Le Foll, B. (2017). Opioid-sparing effect of cannabinoids: A systematic review and meta- analysis. Neuropsychopharmacology, 42(9), 1752–1765. https://doi.org/10.1038/npp.2017.51 Ollech, D. (2021). seastests: Seasonality Tests. R package version 0.15.4. [Software]. https://CRAN.R- project.org/package=seastests Pfaff, B. (2008a). VAR, SVAR and SVEC models: Implementation within R package Legal Cannabis Sales and Opioid Deaths in Canada 118 vars. Journal of Statistical Software, 27, 1-32. https://doi.org/10.18637/jss.v027.i04 Pfaff, B. (2008b). Analysis of integrated and cointegrated time series with R. Springer Science & Business Media. Ploberger, W., & Krämer, W. (1992). The CUSUM test with OLS residuals. Econometrica: Journal of the Econometric Society, 60(2), 271- 285. https://doi.org/10.2307/2951597 Powell, D., Pacula, R. L., & Jacobson, M. (2018). Do medical marijuana laws reduce addictions and deaths related to pain killers?. Journal of Health Economics, 58, 29–42. https://doi.org/10.1016/j.jhealeco.2017.12.007 Public Health Agency of Canada. (28 June, 2023). Modelling opioid-related deaths during the overdose crisis. Ottawa: Public Health Agency of Canada. https://www.canada.ca/en/health- canada/services/opioids/data-surveillance- research/modelling.html Public Health Ontario. (2023, March). Interactive Opioid Tool. https://www.publichealthontario.ca/en/Data- and-Analysis/Substance-Use/Interactive- Opioid-Tool#/trends R Core Team (2024). R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. https://www.R-project.org Shover, C. L., Davis, C. S., Gordon, S. C., & Humphreys, K. (2019). Association between medical cannabis laws and opioid overdose mortality has reversed over time. Proceedings of the National Academy of Sciences, 116(26), 12624–12626. https://doi.org/10.1073/pnas.1903434116 Statistics Canada. (2018, June). Canadian Tobacco, Alcohol and Drugs Survey (CTADS): Summary of results for 2017. https://www.canada.ca/en/health- canada/services/canadian-alcohol-drugs- survey/2017-summary.html Statistics Canada. (22 November, 2023a). Table 20-10-0056-01 Monthly retail trade sales by province and territory (x 1,000) [Detailed Tables]. https://www150.statcan.gc.ca/t1/tbl1/en/tv.act ion?pid=2010005601 Statistics Canada. (2023b). Table 18-10-0004-10 Consumer Price Index, monthly, percentage change, not seasonally adjusted, Canada, provinces, Whitehorse and Yellowknife - Alcoholic beverages, tobacco products and recreational cannabis. [Detailed Tables]. https://www150.statcan.gc.ca/t1/tbl1/en/tv.act ion?pid=1810000410 Statistics Canada. (2023c). Table 17-10-0009-01 Population estimates, quarterly. [Detailed Tables]. https://www150.statcan.gc.ca/t1/tbl1/en/tv.act ion?pid=1710000901 Toda, H. Y., & Yamamoto, T. (1995). Statistical inference in vector autoregressions with possibly integrated processes. Journal of econometrics, 66(1-2), 225-250. https://doi.org/10.1016/0304-4076(94)01616-8 Trapletti, A., Hornik, K., & LeBaron, B. (2007). tseries: Time series analysis and computational finance. [Software]. R package version 0.10-11. Zeileis A, Kleiber C, Krämer W, Hornik K (2003). “Testing and dating of structural changes in practice.” Computational Statistics & Data Analysis, 44(1-2), 109-123. https://doi.org/10.1016/S0167-9473(03)00030- 6 Zeileis A, Leisch F, Hornik K, Kleiber C (2002). “strucchange: An R Package for Testing for Structural Change in Linear Regression Models.” Journal of Statistical Software, 7(2), 1-38. https://doi.org/10.18637/jss.v007.i02 Funding and Acknowledgements: This project was supported by the Peter Boris Chair in Addictions Research and Canada Research Chair in Translational Addiction Research (CRC-2020- 00170). AJM (MFE-193995) is supported by a Postdoctoral Fellowship from the Canadian Institutes of Health Research (CIHR). JM is a principal in Beam Diagnostics, Inc and has consulted to Clairvoyant Therapeutics, Inc. All other authors have no conflicts of interest to declare. Copyright: © 2024 Authors et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction, provided the original author and source are credited, the original sources is not modified, and the source is not used for commercial purposes. https://creativecommons.org/licenses/by/4.0/ Cannabis, A Publication of the Research Society on Marijuana 119 Citation: McDonald, A. J., Cooper, A., Doggett, A., Belisario, K., & MacKillop, J. (2024). Legal recreational cannabis sales and opioid-related mortality in the 5 years following cannabis legalization in Canada: A Granger causality analysis. Cannabis, 7(3), 106–119. https://doi.org/10.26828/cannabis/2024/000261 Issue Date: December 12, 2024