Risk factors for the development of opioid use disorder after first opioid prescription: a Swedish national study

Abstract

Background

We need to better understand the frequency and predictors of opioid use disorder (OUD) after first opioid prescription (OP).

Methods

We followed 1 516 392 individuals from the Swedish population born 1980–2000, from 1 July 2007, until 31 Dec 2017. We examined putative risk predictors with univariable and multivariable Cox Models and the potential causal effects of predictors by propensity score and co-sibling analyses.

Result

Of the individuals in our cohort, 24.8% (375 404) received a first OP, of whom 3034 (0.90%) developed a subsequent first OUD. The hazard ratio (HR) (± 95% CIs) for OUD after OP equaled 7.10 (6.75–7.46), with a mean time to onset of 3.41 (2.39) years. The strongest putative risk factors for development of OUD after OP were prior psychiatric and substance use disorders, criminal behavior, parental divorce/death, poor school performance, current community deprivation, divorce, and male sex. Few predictors differed across sexes. OP renewal was associated with a HR of 3.66 (3.41–3.93) for OUD. Co-sibling and propensity score analyses suggested that at least a moderate proportion of the risk factor-OUD association was likely causal. A risk score to predict OUD after OP had an AUC of 0.85, where nearly 60% of cases scoring in the top decile.

Conclusions

In a general population sample, an OP represents a substantial risk factor for subsequent OUD. Many of the risk factors for OUD after OP can be readily assessed at the time of potential OP, permitting clinicians to evaluate the risk of iatrogenic OUD.

Opioid use disorder (OUD) is a public health problem of major proportions (Beauchamp, Winstanley, Ryan, & Lyons, 2014). Because a significant proportion of OUD develops after receipt of an opioid prescription (OP) (Edlund et al., 2014; Fishbain, Cole, Lewis, Rosomoff, & Rosomoff, 2008; Kolodny et al., 2015; Minozzi, Amato, & Davoli, 2013; Voon, Karamouzian, & Kerr, 2017), the indications for OP for treatment of pain has come under increased scrutiny with a need to balance the importance of pain control with the long-term risk of iatrogenic OUD (Sehgal, Manchikanti, & Smith, 2012). This public health importance of these issues has grown along with the large rise in rates of opioid use and opioid-related deaths in recent decades in most countries of Europe and North America including Sweden (Andersson, Håkansson, Krantz, & Johnson, 2020; Fugelstad, Thiblin, Johansson, Ågren, & Sidorchuk, 2019; Hastie, Gilson, Maurer, & Cleary, 2014).

In response to this urgent question, empirical studies have sought to predict OUD risk after OP (Boscarino et al., 2010; Cochran et al., 2014; Klimas et al., 2019; Reps, Cepeda, & Ryan, 2020), most of which have had modest sample sizes, samples of convenience, limited predictors and/or short follow-up periods (Beauchamp et al., 2014) (Minozzi et al., 2013). In this article, we attempt to address many of the methodological concerns of these earlier studies, examining five specific questions:

1. 1. How much does receipt of an OP increase risk for subsequent OUD?

2. 2. Explore, in univariable and multivariable analyses, a wide range of risk factors assessed at the time of OP that might predict development of OUD.

3. 3. Determine whether duration of OP alters risk for subsequent OUD.

4. 4. Perform propensity-score matching and co-sibling analyses to gain insight into the causal nature of the association between our predictors and OUD risk.

5. 5. From multivariable analyses, develop, in a training split-half our sample, an aggregate risk score for OUD onset after OP and then, in the test split-half, evaluate its performance.

Methods

We linked nationwide Swedish registers via the unique 10-digit identification number assigned at birth or immigration to all Swedish residents. The identification number was replaced by a serial number to ensure anonymity. We secured ethical approval for this study from the Regional Ethical Review Board of Lund University (No. 2008/409, 2012/795, and 2016/679).

The following sources were used to create our dataset: Total Population Register, containing information about year of birth, and sex; Multi-Generation Register, linking individuals born after 1932 to their parents; the Longitudinal Integration Database for Health Insurance and Labor Market Studies (LISA) containing information about education and income from 1990 to 2018; the Hospital Discharge Register, containing hospitalizations for Swedish inhabitants from 1964 to 2018; the Day Surgery Register, containing diagnosis from 1997 to 2000; the Prescribed Drug Register, containing all prescriptions in Sweden picked up by patients from July 2005 to 2018; Outpatient Care Register, containing information from all outpatient clinics from 2001 to 2018; Crime Register that included national complete data on all convictions in lower court from 1973–2018; Swedish Suspicion Register that included national data on individuals strongly suspected of crime from 1998–2018; and the Mortality Register with dates and causes of death from 1952 until 2018. In addition, we had medical diagnosis from Primary Health Care clinics from counties in Sweden which counties included 87% of the Swedish population [for details see (Sundquist et al., 2017)]. We included everyone born between 1980 and 2000 without a registration of OUD prior to July 2007.

Measures

The outcome measure was OUD which we assessed using a medical registration of the ICD-10 code of F11 from any of the above registers, or from the Prescribed Drug Register using the ATC code representing prescribed drugs for opioid dependence; N07BC. In addition, we used criminal convictions from 1996 to 2006 and in 2009 when the type of drug used was registered. OPs were identified from the Prescribed Drug Register using ATC code for opioids; N02A.The outcome measure was OUD which we assessed using a medical registration of the ICD-10 code of F11 from any of the above registers, or prescribed drugs for opioid dependence; N07BC. In addition, we used criminal convictions from 1996 to 2006 and in 2009 when the type of drug used was registered. OPs were identified from the prescribed drug register using ATC code N02A. As we want to assess the first OP, we censored the first two years from the prescribed drug register to reduce the probability that we would count a repeat OP as a first OP. Thus, we only counted prescriptions after July 2007 when no prior opioids are prescribed. Using the Swedish Register Data, we assessed a number of possible predictors which are all outlined in online Supplementary Appendix Table 1.

Statistical methods

We began by assessing the overall risk of OUD as a consequence of OP utilizing a Cox proportional Hazard model with first OP as a time dependent covariate. We follow individuals from July 2007 until first OUD diagnosis, death, emigration, or end of follow-up, whichever comes first. We estimated both the raw association and adjusted for fixed covariates including sex, parental education, parental divorce, educational achievement, and Familial Genetic Risk Score for drug use disorder (DUD) (FGRS_DUD).

In the second set of analyses, we focused on individuals who have received an OP. Using Cox proportional hazard models we follow individuals from time of first OP until first OUD diagnosis, death, emigration, or end of follow-up.

We first estimate the overall raw associations for each of the predictors in the whole study sample, adjusting for sex. Next, we investigated whether the association differed between males and females by allowing for an interaction term between each predictor and sex. The results from this model were presented as HRs for males and females separately and in addition, the ratios between these ratios, representing the interactions, are presented. Finally, we included all predictors in a full, multivariable model. We investigated the proportionality assumption for the categorical variables visually by plotting the Kaplan–Maier curves. For the continuous variables we tested whether the HR was stable during follow-up time.

In a third set of analyses, we addressed the question of causality using two different approaches. First, we use propensity score matching (Austin, 2011) to construct groups that are as equal as possible with respect to the probability of having the exposure (an OP) and then we compare the risk of OUD for the group exposed compared to the ones unexposed, for each one of the predictors. We used the matchIt package in R and the nearest neighbor method to match individuals. Next, we use a co-sibling design that account for genetic and environmental factors shared within sibling pairs (Ohlsson & Kendler, 2020). Siblings discordant on the exposure and outcome, or the time of the outcome contribute to the analysis. To quantify the amount of the raw associations that can be explained be confounding factors, we compare the raw HRs with the HRs obtained from propensity score matching and co-sibling analysis. We compare the linear form of the parameters from the Cox regressions and present, as a percentage, the amount of the association that is explained by cofounding. For further background on these analyses, see online Supplementary Appendix Table 2.

In the final set of analysis to evaluate the predictive effect of our model we split the analytic sample into two random halves and run a logistic model, including all predictors, on the first or training half. Than we evaluate our model using the second, testing half of the sample by estimating the predictive probabilities of OUD given the parameter estimated from the training half and compare with the observed values. The results are presented as a ROC curve with the corresponding AUC value.

Results

Descriptive results

We began with a population sample of 1 516 392 individuals (51.7% males and 48.3% female) born from 1980–2000, followed from 1 July 2007, until 31 Dec 2017. During this follow-up period, 375 404 of them (24.8%) received a first OP, of whom 335 833 had complete data for analysis. After their first OP, 3034 of these individuals (0.90%) were registered for the first time for an OUD. The mean time (s.d.) from OP to OUD diagnosis in this sample was 3.41 (2.39) years and was modestly shorter in males [3.34 (2.35)] than in females [3.53 (2.46)] [p = 0.04]. The risk for OUD over the first 1, 2 and 5 years after first OP equaled 0.17, 0.32 and 0.68%, respectively.

Impact of OP on OUD risk

The raw HR (± 95% CIs) for OUD conditional upon receipt of an OP was estimated (± 95% CIs) at 7.10 (6.75–7.46). When including key fixed covariates including sex, parental education, educational achievement and FGRS_DUD, this HR decreased to 5.81 (5.49–6.15).

Prediction of OUD – Univariable models

The frequency of our putative predictors of OUD and their inter-correlations are seen in online Supplementary Appendix Tables 3 and 4, respectively. Although the mean (s.d.) inter-correlation between them was modest [+0.16 (0.15)], some inter-correlations – for example between prior DUD, alcohol use disorder (AUD) and criminal behavior (CB) –were substantial. Table 1 provides the univariable analysis of our predictors in our entire sample in the left-hand column. Our putative predictors were divided into eight categories in rough chronological order of their occurrence: (i) demographic features, (ii) parental characteristics, (iii) genetic risk, (iv) educational performance, (v) prior psychiatric, and substance use disorders and CB, (vi) prior injuries and pain diagnoses, (vii) current marital status and (viii) community characteristics.

Table 1. Univariable, sex-specific and multivariable analyses of putative predictors of opiate use disorder after opiate prescription

Only one putative predictor – immigration status – was not significantly associated with OUD risk. The strongest risk factors were all in the ‘prior disorder’ category with HRs ranging from 5.17 (4.80 5.56) for major depression or anxiety disorders (DAD) to 12.59 (11.69–13.55) for prior DUD. Three predictors produced HRs in the range of 2–3: low v. high community deprivation and parental education and parental death/divorce. HRs between 1.5 and 2 were seen with male sex and poor grades in high school. The genetic risk had a relatively modest impact with a HR of 1.31 per s.d..

Table 1 also presents results for males and females separately and then, in the right-hand column, the difference between the two HRs. Only 4 predictors were significantly different, and all were more strongly associated in males: immigrant status, prior CB, DAD and bipolar or non-affective psychotic disorders.

While our registry data on the nature of the OP is limited and did not include the specific form of opioid and was often missing information about dosage, we could determine its renewal status in our entire sample. In a univariable analysis, renewal of the OP within six months was a robustly associated with future OUD with a HR of 3.66 (3.41–3.93). This HR, whose impact is over and above the 7-fold increased risk for OUD from the first OP, did not differ significantly in males and females, with a M/F ratio of 0.96 (0.83–1.11).

Prediction of OUD – Multivariable models

The final column in Table 1 presents the multivariable results. As would be expected given the positive correlations between most of our predictors, all the HRs declined with four losing significance (low v. high parental education, unmarried v. married, divorced v. married and % of DUD in neighborhood). The five strongest multivariable predictors were, in order: prior DUD, prior DAD, prior CB, prior AUD and male sex.

Causal effects

While casual effects of putative predictors cannot be evaluated definitively in observational data, we apply two distinct methods to provide insight into the degree to which our observed associations arise from causal effects v. confounders. Propensity Score Matching can be applied to all of our putative risk factors. Our co-sibling method, by contrast, cannot be utilized for risk factors for which full siblings are obligatorily concordant such as parental education or history of divorce. Furthermore, our co-sibling analyses can only be performed on sibling pairs concordant for OP exposure but discordant for the particular risk factor being evaluated. Even within our large sample size, such pairs were sometimes rare, which resulted in parameters estimated with poor precision.

Table 2 presents the results of these causal inference methods where we show the raw univariable findings which we then correct for sex. Next, we present the results from propensity score matching, from which we estimate the proportion of the association that is likely causal. Across all our predictors the propensity score matching, after eliminating the potentially protective effect of being non-Swedish born, the mean (s.d.) of the proportion of the association that is causal equaled 30.5% (25.6). Then we present the parallel results for the co-sibling analyses which estimates a considerably higher proportion of causal effects for our predictors: 62.9% (28.0). Because discordant co-siblings do not control completely for genetic confounding (sharing 50% of their genes identical by descent), the causal estimates from this method are likely biased upward.

Table 2. Analysis of potential causal effect of putative predictors by propensity score and co-sibling analyses

^a Results missing from co-sibling analyses for predictors that do not differ among siblings.

Both methods suggest that single marital status and prior DUD have strong causal effects on risk for OUD after OP with moderate to substantial causal effects seen for prior CB and DAD and at least modest causal impact for prior AUD, bipolar or non-affective psychotic disorders and pain diagnosis. We also examined the potential causal effect of renewal of the OP within 6 months which had a sex-adjusted HR of 3.66 (3.41, 2.93). The HR estimated from propensity score and co-sib analyses were 2.55 (2.32–2.81) and 1.56 (1.30–1.88), respectively, suggesting that 72% and 34% of its impact on OUD risk was likely causal.

Development of predictive model

We utilized 21 predictors which we examined in a multivariable logistic regression model predicting OUD in a random half of our sample to estimate beta parameters weights. We then applied those beta- parameters to the second half of our sample with the results displayed in Table 3. The pattern of findings results was, as expected, similar to that seen in the multivariable HR analyses in Table 1 with the strongest associations seen for prior DUD, DAD, CB, AUD and sex. As seen in Fig. 1a, 60% of all cases of OUD were found in the top decile of the risk score which has an OR per decile of 1.79 (1.73–1.84). For the AUC curve, equal to 0.84, see Fig. 1b.

Fig. 1. (a) The performance of our risk prediction model – The Y axis reflects the observed risk for OUD while the X-axis represents the decile of the risk score generated from our logistic regression model on the training split half of our sample and then applied to the testing split-half sample. (b) An AUC curve demonstrating the performance of our risk prediction model as applied to our testing split-half sample. The Y-axis represents the specificity of the model in predicting OUD while the X-axis represents 1- specificity.

Table 3. Prediction of OUD from a multivariable logistic regression confirmed on a second split-half of the sample

Discussion

We addressed five questions in this paper all related to the delineation of risk predictors of the onset of OUD after a first OP. We review these results in turn.

First, we found that the increased risk for OUD association with an OP in the Swedish population was large, with a raw HR of 7.10 which was reduced to 5.86 with standard covariates. This is a much larger increased risk than recently seen in one-year follow-up of nearly 50 000 patients who received an OP as part of an ER visit (HR, 1.1; 95% CI 0.9–1.4) (Punches, Ancona, Freiermuth, Brown, & Lyons, 2021), and more modest than the nearly 15-fold increased risk (5.8 v. 0.4%) risk of an ‘opioid-abuse related diagnosis’ in a one year follow-up of a sample of over 750 000 adolescents and young adults from the Optum Research Database receiving a first time OP during a dental visit (Schroeder, Dehghan, Newman, Bentley, & Park, 2019). The magnitude of the OP-OUD association found in the Swedish general population highlights the public health importance of the risk for iatrogenic OUD. It also suggests that the risk for OUD after OP is not restricted to individuals whose OP is for the treatment of chronic pain.

We found a relatively low absolute risk for OUD after OP of 0.90%, considerably below the 4.7% obtained from a recent large-scale meta-analysis by Higgins et al. (Higgins, Smith, & Matthews, 2018). However, the rates varied widely across the 12 studies they reviewed, from 0.2% to 34.2% and most of the follow-up periods short, under one-year in duration (Higgins et al., 2018). Only one study they reviewed was a one-year prospective cohort design which found rates of 0.5% progression to OUD (Cepeda, Fife, Ma, & Ryan, 2013), about three-times higher than we observed over our first year of follow-up. Many of the studies reviewed by Higgins et al. (Higgins et al. 2018) focused on patients with chronic pain and this likely represents a group at considerably higher risk for OUD than our general population sample. As expected, our rates are considerably lower than those reported by Parker and Anthony for adolescents developing OUD within 12 months of starting to use prescription opioids extra-medically which ranged across ages and cohorts from 1.8% to 9.8% (Parker & Anthony, 2015).

Second, our univariable analyses revealed a number of risk factors for the development of OUD after OP from a range of different domains. Consistent with prior evidence, we found, in our univariable analyses, that the strongest predictors of development of OUD after OP was prior psychopathology (Boscarino et al., 2010; Burcher, Suprun, & Smith, 2018; Cochran et al., 2014; Hassan, Le Foll, Imtiaz, & Rehm, 2017; Rice et al., 2012) and substance use disorders (Burcher et al., 2018; Rice et al., 2012). Many of the other predictors of OUD development after OP that we observed in the Swedish population were supported by previous studies including being single (Burcher et al., 2018), low SES (Burcher et al., 2018), high familial risk for substance abuse (Burcher et al., 2018; Rice et al., 2012), being male (Burcher et al., 2018; Cochran et al., 2014), a young age at first OP (Burcher et al., 2018; Cochran et al., 2014), and having a pain disorder (Blanco et al., 2016; Cochran et al., 2014). When our analyses were repeated with a multivariable design, most HRs declined substantially, although the strongest predictors remained prior psychopathology and substance use disorders. These results highlight the multifactorial nature of the development of OUD after OP, consistent with our prior studies of all forms of DUD in Sweden (Kendler, Ohlsson, Edwards, Sundquist, and Sundquist, 2017).

We can particularly compare our results with those of a one year follow-up study of OUD after first OP in four large US data bases (Reps et al., 2020). Rates of OUD onset varied from 0.04–0.26% across cohorts, compatible with the one-year rate of 0.17% found in our sample. Furthermore, they generated an 8-item OUD predictor, six of which we replicated in our sample: young age, prior DUD, prior depression, prior anxiety disorder, and two medical conditions associated with chronic pain (Reps et al., 2020).

Third, in accord with prior studies, we found that longer opioid therapy was correlated with a substantially higher subsequent OUD risk (Cochran et al., 2014; Higgins et al., 2018).

Fourth, we utilized two different designs – co-sibling and propensity-score matching – to gain insight into the nature of the associations between our predictors of OUD risk after OP. These designs are quite different in their assumptions and potential limitations so that the joint use can be best seen as an example of triangulation (Munafo & Davey-Smith, 2018).These analyses are of importance for planning interventions based on putative risk factors as the reduction in caseness that might arise from such interventions would likely be closely related to the proportion of causal association between the risk factor and disease. We found reasonable but not complete agreement across our two methods with associations declining substantially across virtually all our predictors. In aggregate, across our entire range of predictors, as well as the strongest category of prior psychiatric and substance use disorders, approximately one third of the association appeared to be causal in its effects. We are not aware of prior studies attempting to predict rates of OUD after OP which have performed similar analyses.

Fifth, using a split-half method, we generated an aggregate risk score from our 21 predictor variables all of which could, theoretically, be assessed at the time at which a clinician is considering prescribing an opioid. The risk score had good sensitivity. Nearly 60% of all cases had scores in the highest decile and the area under the AUC curve equaled 0.85. However, the score had low specificity. In the highest risk decile, only 7% developed OUD. Many other assessment tools have been developed to predict risk for OUD after OP (Klimas et al., 2019). Only further empirical studies will enable a rigorous comparison of the utility and predictive ability of our measure.

How should the performance of our risk score be interpreted in the context of the results of our causal analyses? To our knowledge, such analyses have not been previously conducted on prior risk measures for OUD post OP (Klimas et al., 2019). Given the level of confounding that likely exists between our putative predictors and OUD, the expected benefits in reductions in OUD for interventions based on these risk factors would be less than expected from the raw HRs. That is, the risk score predictions will accurately predict caseness for OUD in a sample exposed to OP, but not all of that effect will be due to the OP. That is because a proportion of the subjects would have developed OUD in any event, given our evidence for confounding variables that predict both receipt of an OP and risk for OUD.

We found that 24.8% of our population based cohort had received and picked-up at least one OP during out 10 year follow-up period. Locating comparable data was difficult. We found one study based on the United States 2015 National Survey on Drug Use and Health, which estimated that 37.8% of their adult sample had used prescription opioids in the last year (Han et al., 2017). These results suggest that rates of OP in Sweden are likely substantially lower than those seen in recent years in the United States.

Limitations

These results should be interpreted in the context of five potential methodological limitations.

First, our results are directly relevant only to the Swedish population although several factors suggest these findings likely generalize to other Western countries. We found substantial similarity in the risk factors for OUD that we identified and those reported in the prior literature. Furthermore, like the US, Sweden has seen increases in opioid consumption over recent decades (Hastie et al., 2014) as well as marked rises in opioid-related deaths (Andersson et al., 2020; Fugelstad et al., 2019). Like the US (Baldwin, Seth, & Noonan, 2021), there has been a shift in Sweden toward the use of short-acting opioids, especially oxycodone and fentanyl, with their high addictive potential and risk for overdose death (Baldwin et al., 2021; Heilig & Tägil, 2018) (Muller, Clausen, Sjøgren, Odsbu, & Skurtveit, 2019). Sweden, like the US, has substantial rates of non-medical use of prescription opioids (Han et al., 2017; Novak et al., 2016). Rates of OP in Sweden are in the range of those observed in US general population samples (Jeffery et al., 2018).

Second, while we detected subjects with DUD and OUD from medical, legal and pharmacy records which require neither respondent cooperation nor accurate recall, they can produce false negative and false positive diagnoses. While large interview studies of DUD prevalence are not available in Sweden, lifetime prevalence in near-by Norway was similar to the rates we detected in Sweden (Kringlen, Torgersen, & Cramer, 2001). It is not like that we substantially under-ascertain moderate to severe cases of DUD or OUD. The validity of our definition of DUD in Sweden is also supported by the high rates of concordance for registration observed across our different ascertainment methods (Kendler et al., 2015) and heritability estimates similar to those obtained from personal interview based studies (Kendler, Karkowski, Neale, & Prescott, 2000; Kendler, Maes, Sundquist, Ohlsson, & Sundquist, 2013; Tsuang et al., 1996).

Third, our estimates of the effect size of predictors of OUD and OP might be biased by baseline effects in which certain of these risk factors might have impacted on the clinician's decision to provide an OP in the first place.

Fourth, we do not have population level data on tobacco usage in Sweden, so are unable to model its impact on risk for OUD after OP.

Finally, we utilized a two-year period ‘buffer’ period to detect prior OPs which might be too short, resulting in some of our first OPs having previously received OPs. To evaluate this possible bias, we repeated our key analyses expanding the ‘buffer’ period to five years. As seen in online Supplementary Appendix Table 5, the pattern of predictors of OUD with 2 and 5 year buffer periods were very similar.

Conclusion

In a Swedish national sample of individuals who received an OP for any indication, the risk for subsequent OUD was increased seven-fold. For those affected, the mean time from OP to first registration for OUD was 3.5 years. In univariate analyses, the development of OUD post OP was most strongly predicted by prior psychiatric and substance use disorders, a history of parental divorce, current community deprivation, divorce, poor prior school performance and male sex. A renewal of that first OP was associated with a 3.5-fold elevated risk for OUD. Co-sibling and propensity score analyses suggested a moderate to substantial proportion of the risk factor-OUD association was likely causal. Using split-half methods, we developed a predictive risk score with an AUC of 0.85, where nearly 60% of cases scored in the top decile. Our results document the public health importance of the substantial increase in OUD risk after OP and provide predictive models that permit clinicians to gauge, with at least moderate accuracy, individuals at particularly high risk for the development of OUD after OP.