Article Text

Concurrent prescriptions for opioids and benzodiazepines and risk of opioid overdose: protocol for a retrospective cohort study using linked administrative data
  1. Erin Y Liu1,2,
  2. Robyn Tamblyn2,3,
  3. Kristian B Filion3,4,
  4. David L Buckeridge2,3
  1. 1Department of Epidemiology, Biostatistics and Occupational Health, McGill University, Montréal, Quebec, Canada
  2. 2McGill Clinical and Health Informatics, McGill University, Montréal, Quebec, Canada
  3. 3Departments of Medicine and of Epidemiology, Biostatistics and Occupational Health, McGill University, Montréal, Quebec, Canada
  4. 4Centre for Clinical Epidemiology, Lady Davis Institute, Jewish General Hospital, Montréal, Quebec, Canada
  1. Correspondence to Erin Y Liu; erin.y.liu{at}


Introduction Opioid overdoses have increased substantially over the last 20 years, with over 400 000 deaths in North America. While opioid prescribing has been a target of research, benzodiazepine and opioid co-intoxication has emerged as a potential risk factor. Our aim was to assess the risk of opioid overdose associated with concurrent use of opioids and benzodiazepines relative to opioids alone.

Methods and analysis A retrospective cohort study will be conducted using medical claims data from adult residents of Montréal, Canada. We will create a cohort of new users of opioids (ie, no opioid dispensations in prior year) in 2000–2014 from people with at least 2 years of continuous health insurance. Those with any diagnosis or hospitalisation for cancer or palliative care in the 2 years before their first opioid dispensation will be excluded. On each person-day of follow-up, exposure status will be classified into one of four mutually exclusive categories: (1) opioid-only, (2) benzodiazepine-only, (3) both opioid and benzodiazepine (concurrent use) or (4) neither. Opioid overdose will be measured using diagnostic codes documented in the hospital discharge abstract database, physician billing claims from emergency department visits and death records. Using a marginal structural Cox proportional hazards model, we will compare the hazard of overdose during intervals of concurrent opioid and benzodiazepine use to intervals of opioid use alone, adjusted for sociodemographics, medical and psychiatric comorbidities, and substance use disorders.

Ethics and dissemination This study is approved by the McGill Faculty of Medicine Institutional Review Board and the Commission d’access à l’information (Québec privacy commission). Results will be relevant to clinicians, policymakers and other researchers interested in co-prescribing practices of opioids and benzodiazepines. Study findings will be disseminated at relevant conferences and published in biomedical and epidemiological peer-reviewed journals.

  • epidemiology
  • adverse events
  • public health

This is an open access article distributed in accordance with the Creative Commons Attribution Non Commercial (CC BY-NC 4.0) license, which permits others to distribute, remix, adapt, build upon this work non-commercially, and license their derivative works on different terms, provided the original work is properly cited, appropriate credit is given, any changes made indicated, and the use is non-commercial. See:

Statistics from

Request Permissions

If you wish to reuse any or all of this article please use the link below which will take you to the Copyright Clearance Center’s RightsLink service. You will be able to get a quick price and instant permission to reuse the content in many different ways.

Strengths and limitations of this study

  • Our choice to conduct our study in new opioid users reduces the chance for selection bias because prevalent opioid users will be more tolerant of the medication.

  • Other strengths in this study include the use of a validated outcome definition and modelling exposure as time-varying, which will reduce misclassification bias.

  • The use of a marginal structural Cox proportional hazards model can control both time-fixed and time-varying confounders, therefore reducing confounding bias.

  • Our study in a Canadian sample can inform whether previous results are robust across different healthcare systems, types of prescribing behaviour and drug insurance policies.

  • A potential limitation is that our results may not be generalisable to all opioid users, as we are restricted to individuals insured through the public drug plan.


Morbidity and mortality attributable to opioid use have increased sharply in North America over the last 20 years. In the USA, deaths from drug overdose increased from 6.1 per 100 000 in 1999 to 20.7 in 2018,1 with 446 032 of these deaths involving opioid use.2 In Canada, hospitalisations related to opioid overdose rose by 53% from 2007 to 2017,3 with a national rate of 15.5 per 100 000 in 2019.4 Deaths related to opioids also remains high, with Canada’s public health agency recording 1634 deaths from January 2016 to March 2020.5 Since harm from opioid overdoses persist, and they affect younger age groups disproportionately, these alarming trends demand effective and immediate public health actions.

Substantial research conducted in North America has focused on the role of prescription opioids in these overdoses and deaths.6–9 Increasingly though, the co-prescribing of opioids and benzodiazepines and their potential role in causing overdose is receiving attention.10–13 Studies in Canada and the USA showed that 27.6% and 30% of opioid overdose deaths involved benzodiazepines, respectively.14 15 This has paralleled high rates of concurrent use of both medications. A 2017 Alberta study showed that 17.6% of opioid users had overlapping dispensations with a benzodiazepine,16 and a 2014 study using USA pharmacy claims estimated the prevalence to be 10%.17 What is concerning is that the prevalence of concurrent use increased by 41% from 2002 to 2014, despite guidelines from both countries cautioning against combining their use.18 19 In this context, the extent of the role of benzodiazepines in the opioid overdose epidemic demands further research.

Benzodiazepines are a class of anticonvulsant/anxiolytic medications used to treat anxiety, depression, panic disorders, insomnia, seizure disorders, alcohol dependence and musculoskeletal pain.20 When consumed alone, benzodiazepines do not cause respiratory depression, the main consequence of drug overdose. But animal studies have shown that the risk of respiratory depression is elevated when benzodiazepines are taken concurrently with opioids.21 This interaction is due to the presence of receptors for both opioids (mu and delta) and benzodiazepines (GABA) in the brain region responsible for respiratory control.22 Since benzodiazepines act on the inhibitory pathway for respiration and opioids inhibit the excitatory receptors, the activation of both pathways can reduce the respiratory drive more than the activation of either pathway alone. Whether this reduction is a consequence of an additive effect of both drugs inhibiting the respiratory system or an amplification of benzodiazepines on the effects of opioids is unknown. Nevertheless, the biological mechanism suggests greater harm from their combined use, and there is a need to understand if these effects are sufficient to cause respiratory depression at therapeutic doses.

To date, three epidemiological studies conducted in population-based samples of adults have examined the risk of overdose associated with concurrent use of opioids and benzodiazepines. We considered other studies on concurrent use and overdose risk, but due to their use of distinct subpopulations such as veterans,23 24 we precluded them from our literature review as their generalisability is unclear.

The three relevant studies were all retrospective analyses based in the USA that used health insurance claims data. They reported effects ranging from HR=1.2 (95% CI 1.16 to 1.34; Cho et al),25 OR=2.14 (95% CI 2.05 to 2.24; Sun et al),26 to HR=5.05 (95% CI 3.68 to 6.93; Hernandez et al),27 when comparing concurrent opioid and benzodiazepine use to opioid use alone. However, potential methodological limitations make these results difficult to interpret. In Hernandez et al, the exclusion of cohort members based on non-benzodiazepine sedatives use after cohort entry could induce selection bias, with an unknown effect on the point estimate if concurrent users are more, or less, likely to be dispensed these medications. A potential source of selection bias in Sun et al is restriction to prevalent opioid users, since this would limit the sample to those who tolerated their opioid medication long enough to continue its use. Bias from exposure misclassification is also a concern. In the retrospective cohort study conducted by Hernandez et al, comparing concurrent use in the 1–90 days before overdose to opioid-only use in the 1 day before means concurrent users had more opportunity (up to 90 days) for opioid exposure. Longer opportunity for opioid use means concurrent users are more likely to be chronic opioid users,16 17 which independently increases the risk for overdose. The resultant effect would be an overestimation of the association. In Sun et al, misclassification bias from inclusion of person-years of non-opioid use in the reference group could have also overestimated the association. Periods of non-opioid use are less likely to lead to opioid overdose, which would artificially underestimate the rate of overdose in the reference group. Immortal time bias28 due to misclassification of person-time could have manifested in the retrospective cohort study conducted by Cho et al. In their main analysis, use of a time-fixed exposure means patients who initiated opioids and later added a benzodiazepine to their opioid regimen would have their opioid-only (and, by definition, event-free) person-time misclassified as concurrent drug use. This misclassification would bias the estimate downward, underestimating the risk.

In summary, the current evidence base could be affected by important bias. As such, there remains uncertainty about the magnitude of the risk of overdose due to the concurrent use of opioids and benzodiazepines. To reduce this uncertainty, we propose a retrospective cohort study using population-based data sources.

Prescription opioid use has more than doubled in North America from 2001 to 2013, and combined with Europe and Oceania, these three continents account for 95.7% of worldwide use.29 Benzodiazepine use is also high in these regions, with the USA reporting 5.2% of adults filling a prescription in 2008,30 4.9% among those 45 years and older in Canada in 2002,31 17.7% prevalence of use in France in 201232 and 2.4% of Australians from 2002 to 2007.33 Evidently, the use of prescription opioids and benzodiazepines is highly prevalent in many countries. An accurate estimate of the risk of overdose associated with their concurrent use is important for setting public health policy in this area.

Methods and analysis

Data source

This study will use anonymised population-based administrative data from Montréal, Québec, to define our cohort and measure all study variables. The source population consists of 1.4 million people representing a 25% random sample of health-insured Montréal residents sampled from 1 January 1998 to 31 December 2014. For each individual, data on their physician services and community pharmacy drug dispensations come from the Régie de l’assurance maladie du Québec (RAMQ), the provincial health authority that reimburses all publicly insured medical and pharmacy services. Data on hospitalisations are provided by the Ministère de la santé et des services sociaux (MSSS), and mortality data come from the Institut de la statistique du Québec (ISQ). As such, we have linked information on each person’s demographics, inpatient and outpatient medical services, hospitalisations, emergency department (ED) visits, and the date and causes of death from the time they were sampled until the time that an individual loses health insurance through death or emigration from province. Half of residents are also insured for prescription drugs through the RAMQ if they are welfare recipients, aged 65 years or older, or lack other drug insurance coverage (eg, through their employer). Thus, data on dispensed drugs in the community (including the date, chemical name, dose and duration), the prescribing physician and dispensing pharmacist are available for approximately 50% of Montréal residents. The data sources have been validated previously,34 and used extensively for research.35–38

Study population

New opioid users with RAMQ drug insurance will form our study population. These are drug-insured Montréal residents aged ≥18 who started using opioids between 2000 and 2014. We chose to include new users starting on 1 January 2000 to allow a 2-year lookback period for exclusion criteria. A new opioid user will be defined as a person with no opioid dispensations for at least 1 year before their first opioid dispensation in this interval. A new user design was chosen to reduce the chance for selection bias, as prevalent drug users are limited to those that tolerated the drug well enough to continue its use.39 The date of their first opioid dispensation will define the cohort entry date, and opioid users will be followed up from this date until the occurrence of the outcome (defined below) or censoring due to (1) death, (2) cancer diagnosis, (3) palliative care admission, (4) loss of drug or health insurance or (5) end of study period (31 December 2014). Since patients with cancer and those undergoing palliative care will have a different risk–benefit ratio when using opioids, cohort members will be excluded if they have at least one cancer or palliative care-related diagnostic code in the 2 years prior to cohort entry. Patients whose first opioid dispensation is methadone or buprenorphine likely have a history of illicit opioid use, as these medications are indicated for an opioid use disorder. Given that we will be adjusting for history of opioid substance abuse, these patients will remain in the study.


An opioid overdose occurs when opioids, possibly together with other substances, cause sedation, leading to loss of consciousness and respiratory depression.40 We chose to study events where the use of opioids alone or in combination with benzodiazepines induces respiratory depression severe enough to require medical attention. To identify these events, we will use International Classification of Diseases (ICD) codes for opioid poisoning (ICD-9 codes E8500-E8502, 96500-96509 and ICD-10 codes T400-T404, T4060, T4069 and X42) and acute effects of opioids (E9350-E9352 and Y450) (online supplemental table E1, appendix I). The presence of any one of these codes in medical service claims, hospital discharge abstracts or vital statistics data during follow-up will define an event, and follow-up will be censored on this date. We chose to use ICD codes related to opioid poisoning and acute effects of opioids because there are no codes specific to drug overdose. However, validation studies comparing opioid poisoning codes to physician chart review have found positive predictive values ranging from 70.0% to 84.6%.41–43 As such, this outcome definition will provide a high specificity in capturing the mechanism of interest (eg, respiratory depression caused by medication).


To ascertain each person’s opioid and benzodiazepine status, we will calculate their average daily dose using information on the drug type (online supplemental table E2, appendix I), the date the prescription was filled, the dose, the quantity dispensed and the duration of the prescription. To account for different combinations of use and non-use, we will model our exposure using a time-varying approach. On each day of follow-up, exposure status will be classified into one of four mutually exclusive categories: (1) opioid-only, (2) benzodiazepine-only, (3) both opioid and benzodiazepine (concurrent use) or (4) neither (figure 1). Since opioid overdose is an acute event, we will consider exposure status on the day before the event as the primary aetiological window of interest.

Figure 1

The purpose of this figure is to show how exposure will be assigned using three fictitious patients to illustrate. Patient (A) enters the cohort with no opioid and benzodiazepine use, as shown by the section in white. Then, they receive their opioid prescription (blue), which ends right before they begin their benzodiazepine prescription (orange). During their period of benzodiazepine use, however, they receive another opioid prescription that overlaps this benzodiazepine period (red). Since they were using both drugs the day before an overdose, this patient’s event would be exposed to concurrent drug use. Patient (B) enters our study already using benzodiazepines (orange), then stops for a period of time (white) before receiving a short prescription for opioids (blue) and then continues their period of non-use (white). Patient (C) enters the study using neither drug (white), but then begins a long duration of opioid-only use (blue) before experiencing the event. Since they were exposed to opioids only in the day before overdose, their event would be considered unexposed to concurrent drug use.

To ensure equivalency across opioid types, we will use published conversion factors to calculate each drug’s oral morphine milligram equivalent (MME).3 44 To calculate MME, we will use the standard method of multiplying the average daily dose by its conversion factor:

Embedded Image

To standardise benzodiazepine potencies, we will use the defined daily dose (DDD) methodology as recommended by the WHO Collaborating Centre for Drug Statistics Methodology (WHOCC).45 The DDD is the ‘the assumed average maintenance dose per day for a drug used for its main indication in adults’ and is calculated as follows:

Embedded Image

Each drug’s Anatomical Therapeutic Chemical (ATC) category and route of administration has its own recommended daily dose, and values will be obtained from the WHOCC website.


Confounding can occur in our study if we fail to adjust for factors that differ between opioid users and concurrent users, and these factors independently affect one’s risk for overdose. Opioid users additionally prescribed a benzodiazepine are often older,17 26 46 more likely to be women,17 23 26 46 require social assistance,27 47 be a chronic opioid user,16 17 have more pain conditions,46 use more potent opioids,46 48 use other central nervous system medications,23 25 and are more likely to have comorbidities,25 26 mental health23 25 27 46 and substance use disorders.23 26 46 Concurrent users are also more likely to seek out multiple prescribers and pharmacies in a given time period.27 48 49 Since these characteristics are also associated with overdose in studies on opioid prescribing patterns,50 it is important that they are adjusted for in this study.

Across all statistical analyses, we will adjust for the following confounders (table 1): demographic factors such as age, sex, drug insurance type and neighbourhood income quintile, medical comorbidities such as history of myocardial infarction, dementia, and chronic pulmonary disease, mental health comorbidities such as diagnoses and prescriptions for anxiety, depression, and mood disorders, opioid use disorders, non-opioid substance abuse/misuse such as alcohol or nicotine dependence and measures of drug-seeking behaviour such as the number of distinct prescribers and pharmacies used in the past 30 days. Receipt of concurrent muscle relaxants, gabapentinoids or Z-drugs will also be adjusted for since these medications can also depress the respiratory system.

Table 1

Confounders under study

To ascertain confounder status, we will check for its associated diagnostic or service code(s) in the defined time period. Values of covariates will be measured in the baseline period (1 year before cohort entry) and every 30 days during follow-up. For age, sex, drug insurance type, neighbourhood income quintile and history of opioid poisoning, the baseline values will remain fixed throughout follow-up. For medical comorbidities, mental health conditions, opioid use disorder and substance use disorders, once a patient is diagnosed, we will assume these conditions persist. All other covariate values will be updated every 30 days during follow-up.

Statistical analyses

Descriptive analysis

We will first perform descriptive analyses to better understand our study participants. This will include comparing concurrent and non-concurrent users with respect to baseline covariates such as those listed in table 1. Continuous variables will be compared using means and medians, and categorical variables using frequency distributions and proportions.

Primary analysis

Time-varying confounding

Time-varying confounding is a concern when both exposure and confounder values are likely to change over the study period and affect each other.51 Given that opioids and benzodiazepines are often used intermittently, bias could occur if there are factors that both predict a change in benzodiazepine status (eg, development of drug misuse disorders or mental health exacerbation) and independently increase the risk for overdose. The three previous epidemiological studies assessing concurrent use and risk of overdose measured confounders at baseline only.25–27 This approach assumes, often implicitly, that comorbidities and drugs used at baseline remain the same throughout follow-up and that changes in comorbidities and concomitant drugs that occur during follow-up do not impact subsequent exposure to concurrent opioid and benzodiazepines. However, for studies with long follow-up periods and transient exposures, such as our study, this assumption will likely be violated.

To adjust for changes in confounder status without adjusting for variables on the causal pathway, we will employ a marginal structural Cox model with inverse probability of treatment weights (IPTW).52 These models reweight study participants based on their inverse probability of exposure to each treatment group, conditional on past use and all confounders.

Logistic regression for calculation of weights

To construct these weights, we will first develop a multinomial logistic regression (exposure model) for the four exposure categories. Given A=exposure and W=Covariates, the goal is to estimate P(A|W), which is the predicted probability of being exposed to (1) opioid-only, (2) benzodiazepine-only, (3) both opioid and benzodiazepine (concurrent use) or (4) neither, as a function of covariates. Covariates used to generate the exposure model will include all time-varying confounders listed in table 1. These weights will be re-estimated during study follow-up using updated covariate values from the past 30 days. For values that do not change, the previous value will be carried forward. Indicator variables will be used for all categorical covariates and continuous covariates will be tested for non-linear effects using Akaike information criterion to determine the optimal form. Once we convert the log odds into predicted probabilities, P(A|W), we can calculate stabilised weights (P(A)/P(A|W)). This is achieved by dividing P(A), the proportion of the cohort exposed to each opioid and/or benzodiazepine category, by P(A|W), the predicted probability estimated from the exposure model. Use of stabilised weights is recommended because it leads to smaller variance as compared with unstabilised weights, which are simply the inverse of P(A|W).52 53 This ensures that results do not depend on a few individuals with extremely large weights. Similarly, we will truncate any extreme values (eg, 99th or 95th percentiles) and verify that the mean distribution of weights is close to 1 to ensure the robustness of our weights.

Cox proportional hazards regression to estimate effect on overdose

The final Cox proportional hazards regression model will estimate the HR for time to first overdose event, when comparing events during concurrent opioid and benzodiazepine person-time to opioid-only person time, adjusted for baseline confounders, and weighted according to each person’s vector of IPTWs. We will use bootstrapping to construct 95% CIs.

Subgroup analyses

To identify patient characteristics that may modify the opioid and benzodiazepine association with overdose, we will undertake the following secondary analyses. We will stratify the primary analysis by sex (males vs female), age (≥76 vs 65–75 years among everyone on the 65+ drug insurance plan), social deprivation (drug insurance for welfare recipients vs all other), presence of an opioid use disorder or substance use disorder (ever/never), presence of a mental health condition (ever/never), and number and type of distinct opioid prescribers (eg, surgeon, dentist or family physician). To assess differential effects by duration or potency of drug use, we will also stratify the analysis by cumulative days, and standardised dose of opioid and benzodiazepines. To address the possibility that prevalent benzodiazepine users may be less susceptible to respiratory depression due to increased tolerance, we will repeat the primary analysis restricted to those with no use of benzodiazepines at baseline.

In these secondary analyses, the stratification variables will be omitted from the multinomial regression model used to calculate IPTWs. Instead, they will be used to separate the cohort into subgroups and each subgroup will be fit with their own Cox regression model to obtain strata-specific estimates. We intend these analyses to be hypothesis generating for future researchers to identify high-risk patient groups.

Sensitivity analyses

Exposure measurement error

To account for the biological half-lives of opioids and benzodiazepines,54 and the fact that many patients may take their medications on an ‘as needed’ (or ‘PRN’) basis, we will introduce a grace period into our exposure assessment, where an additional 7 and 30 days of drug duration will be added to the end of each prescription. This is to account for situations where some patients may dispense their drug but keep it in their medicine cabinets until it is required. Since chronic opioid users have more medication in supply, this also increases their chance for overdose if they ingest more than the recommended dose at once.

Timing of exposure

Although opioid and benzodiazepine-induced overdose is an acute effect, we will assess the 30-day period before the event date as a secondary aetiological window of interest.

Outcome measurement error

Given the challenges of attributing causality in drug poisoning, we will repeat our primary analysis using the broader definition of ‘any drug poisoning’ to account for outcome measurement error. Additionally, we will analyse fatal and non-fatal overdoses separately, and given the potential for overdose deaths to be misclassified, we will also assess all-cause mortality as a secondary outcome.

Strengths and limitations

Our study will contribute evidence to an important public health question. Given the high prescribing rates of opioids and benzodiazepines, any risk of overdose from their concurrent use would have a large impact at a population level.

We believe our study has several strengths. Restricting our cohort to new opioid users reduces the chance for selection bias because prevalent opioid users will be more tolerant of the medication. Our use of a validated outcome measure will reduce misclassification bias and the sensitivity analysis using the broadened definition can test robustness of our results. By modelling use of opioids and benzodiazepines as time-varying, we will reduce exposure misclassification bias. This is particularly important since prescriptions for these two medications tend to be intermittent in nature. By employing a marginal structural Cox proportional hazards model, we can reduce bias from both time-fixed and time-varying confounders. Conducting observational studies with varied statistical techniques in diverse populations will strengthen the evidence base regarding this issue. Finally, our study in a Canadian sample can inform whether previous results are robust across different healthcare systems, types of prescribing behaviour and drug insurance policies.

A limitation is that our results may not be generalisable to all opioid users, as we are restricted to individuals insured through the public drug plan. The drug plan is limited to all persons 65 years or older, those on social assistance, and those without employer-based access to drug insurance. However, given that we are measuring a biological drug effect and population-level variation in drug metabolism is unlikely to differ substantially, our results should be generalisable to most patients prescribed these medications. Another limitation is that exposure to non-prescribed sources of opioids and benzodiazepines, such as those obtained illicitly, or from friends and family, will be missing in our analysis. We believe any misclassification will be minimal, though, since our exposure contrast includes prescription opioid use in both the exposed (concurrent opioid and benzodiazepine) and reference (opioid-only) groups. Moreover, population-based surveys such as the Canadian Tobacco, Alcohol and Drugs Survey show that in 2015, 2.3% of Canadians aged 15 and older reported any non-cannabis, illicit drug use in the past year, and 0.3% had abused pharmaceutical pain relievers.55 This small percentage of illicit opioid use in the population, in addition to the fact that our study period (2000–2014) predates the rise of synthetic, illicit opioids (eg, fentanyl) and their involvement in overdose deaths,56 means it should not impact our study findings greatly.

Patient and public involvement

This research was done without patient involvement. Patients were not invited to comment on the study design and were not consulted to develop patient relevant outcomes or interpret the results. Patients were not invited to contribute to the writing or editing of this document for readability or accuracy.

Ethics and dissemination

This study is approved by the McGill Faculty of Medicine Institutional Review Board and the Commission d’access à l’information (Québec privacy commission). Since we will be using pseudo-anonymised data, no consent is required.

Our findings will be important to researchers, public health and medical communities in Canada and abroad, and this will be reflected in our dissemination plan. We plan to first present our findings at pharmacoepidemiology (eg, the International Conference on Pharmacoepidemiology & Therapeutic Risk Management) and public health/health services conferences (eg, the Canadian Association for Health Services and Policy Research and Public Health). We hope to elicit feedback from our scientific and clinical colleagues at these venues. As our methods are novel and other epidemiologists may benefit from learning about our experience applying these methods, we will also be submitting abstracts to the Society for Epidemiology’s Annual Meeting and the Canadian Society for Epidemiology and Biostatistics conference. Our partners at the National Institute of Public Health of Québec (INSPQ) and Canadian Institutes of Health Research (CIHR)’s Drug Safety and Effectiveness Network (DSEN) will be important knowledge users. We anticipate our work will inform decisions around implementing prescription monitoring programmes and guidelines on prescribing of opioids and benzodiazepines. We will reach out to our public health partners to present this work in their organisational seminars.

After integrating feedback from conference peers, we will submit our manuscript to peer-reviewed biomedical journals (eg, BMJ, AJPM) as primary care physicians are responsible for the majority of opioid and benzodiazepine prescriptions. All publications will be reported in accordance with the REporting of studies Conducted using Observational Routinely collected health Data specific to pharmacoepidemiological research (RECORD-PE).57


The authors are grateful for the comments and recommendations provided by Dr Ari Nandi, Dr Robert Platt and Dr Michal Abrahamowicz during their review of an earlier draft of this protocol.


Supplementary materials

  • Supplementary Data

    This web only file has been produced by the BMJ Publishing Group from an electronic file supplied by the author(s) and has not been edited for content.


  • Twitter @ErinYiranLiu

  • Contributors This study was conceived by EYL and all authors contributed to the design, analysis and interpretation of the data. EYL drafted the manuscript and DLB, RT and KBF revised it critically for important intellectual content. All authors gave final approval of the version to be published and EYL will be the guarantor of the work.

  • Funding EYL is supported by the McGill R.H Tomlinson Fellowship, the CIHR Frederick Banting and Charles Best Doctoral Award, the CIHR Drug Safety and Effectiveness Cross-Disciplinary Training Program (DSECT) award, and the Fonds de recherche du Québec – santé (FRQS; Quebec Foundation for Health Research) Doctoral Award. KBF is supported by a salary support award from the FRQS and a William Dawson Scholar award from McGill University. DLB is supported by a Canada Research Chair in Health Informatics and Data Science (CPP-137904). RT is supported by a CIHR Operating Grant (IOP-112675).

  • Competing interests None declared.

  • Patient consent for publication Not required.

  • Provenance and peer review Not commissioned; externally peer reviewed.

  • Supplemental material This content has been supplied by the author(s). It has not been vetted by BMJ Publishing Group Limited (BMJ) and may not have been peer-reviewed. Any opinions or recommendations discussed are solely those of the author(s) and are not endorsed by BMJ. BMJ disclaims all liability and responsibility arising from any reliance placed on the content. Where the content includes any translated material, BMJ does not warrant the accuracy and reliability of the translations (including but not limited to local regulations, clinical guidelines, terminology, drug names and drug dosages), and is not responsible for any error and/or omissions arising from translation and adaptation or otherwise.