The Iranian Drug-related Deaths with Time Series Modeling and Forecasting: 2014–2016

Background: Investigating the temporal variations and forecasting the trends in drug-related deaths can help prevent health problems and develop intervention programs. Iran's recent policy is strongly focused on deterring drug use and replacing illicit drugs with legal drugs. This study was conducted to investigate drug-related deaths in 2014-2016 in Iran and forecast the death toll by 2019. Methods: This longitudinal study used a Box-Jenkins time series analysis to forecast referrals with a diagnosis of drug-related death. The number of referrals was extracted from March 2014 to March 2017 by month. After using appropriate data transformation methods to create a stationary time series and carrying out a more accurate examination of the stationary assumption by the Dickey-Fuller test, the ARIMA model parameters were determined using Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) graphs. Comparing the Akaike statistics of the proposed models led to the selection of ARIMA (0,1,2) as the best model for t to the data. In the nal stage, the number of deaths was forecast in Iran until 2019 using the ARIMA (0,1,2) model. The nal extracted data were analyzed in R software, Minitab and SPSS-23. Results: The death toll in Iran during the entire study period of three years was 8883 according to the Iranian Ministry of Health and the Legal Medicine Organization. The number of deaths varied by year and was 2840 in 2014, 2810 in 2015 and 3233 in 2016. According to the time-series ndings, data on the number of deaths showed an increasing trend and was not under any regular seasonal effects. In addition, the mean number of deaths per month in Iran until 2019 was forecast as 245.8. Conclusion: This study shows that the trend of drug-related deaths rose during the study period, and the modeling process for their forecasting suggests that this trend shall continue until 2019 if proper interventions are not instituted.


1-Background
Illicit drug abuse is a serious concern with grave consequences for public health in Iran [1]. The prevalence of drug abuse in Iran and other Muslim countries, such as Afghanistan and Pakistan, is much higher than that in western societies [ 2]. Iran's long history of drug use and neighboring Afghanistan as the largest drug producer in the world facilitate Iranian's access to drugs [3 ]. Alcohol and drug abuse disorders make up about 2% of the burden of diseases in Iran, and although similar to global statistics, they comprise one of the parameters placing a heavy burden on the Iranian health system [4]. Many policies have been implemented in Iran to reduce the harm caused by drug abuse, including the Harm Reduction Program by Methadone Maintenance Treatment since 2005. Iran is on the drug transit route, and border surveillance activities as well as the execution of drug tra ckers since January 2018 have been some other policies taken to reduce the drug supply, demand and harms [5][6][7][8]. Moreover, the ageadjusted rate has decreased in overdose and drug-related deaths compared to the past decade and various patterns of smoking, alcohol consumption and drug use are now observed among adolescents in Iran. Government o cials should therefore revise their policies on harm, supply and demand reduction [3,9 ].
In the long run, the mortality rate is six to 54 times higher than expected among drug users, and the rate is eight times higher in patients who have recently been treated for drug poisoning compared to the general population [10]. Histories of severe psychiatric illness, drug poisoning and suicidal drug poisoning are risk factors for mortality in this group of the population [10,11 ]. In Iran, despite the presence of harm reduction programs, the rate of fatal drug poisoning is still 152% higher than that in the past two decades [12]. This study was therefore conducted to examine temporal changes such as increasing or decreasing trends of fatal poisoning (drugs, stimulants and alcohol) in Iran in 2014-2016 and forecast its future trends with Autoregressive Integrated Moving Averages (ARIMA). This method is the most widely-used model for short-term forecasting and investigating the stationarity and seasonality of a phenomenon, and is a powerful analytical tool with simple interpretations; however, it remains rarely used in public health programs [13]. This forecast process can provide valuable information and credible evidence for policymaking. This study also aimed to examine the trend and make forecasts for the future in order to enable better decision-making and ultimately prevent the phenomenon in question.

2-Methods
This time-series study was conducted in 2019 following publication of a paper which contained parts of studies by the same authors in 2019, titled "An estimation of drug-related deaths in Iran using the capturerecapture method" with the ethics code IR.MAZUMS.REC.1398.445. To collect the data, a letter of permission was sent by Mazandaran University of Medical Sciences to the research settings (LMO and MOHME), which issued an authorization letter for data collection. Drug-related death data for the years 2014-2016 were collected from the Iranian MOHME and LMO with the relevant codes [14] using the capture-recapture method. The data include records from March 2014 until March 2017 in the form of the number of referrals by month with relevant diagnostic codes in the International Classi cation of Diseases (ICD-10). The trend of deaths was evaluated over a period of 36 months.

a. Modeling
This study used time series analyses of Box-Jenkins and Box-Cox models to forecast records of drugrelated deaths. The Kolmogorov-Smirnov test was used to test the normality and Bartlett's test to examine the equality of variances in drug-related death data. In order to analyze the observations of the whole country's death toll in different months from 2014 to 2016, the time-series graph of the observations, Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) dependency quantity graphs and series structural decomposition graphs were drawn. The graphs showed non-stationarity in the mean and variance. The Cox-Box transformation was used to stabilize the variance for a more detailed analysis. The estimated value of the λ parameter in the Cox-Box transformation was 2.27. The Cox-Box transformation with a parameter of 2.27 and then a rst-order differential operator were therefore applied to the main observations. For further investigation, the Cox-Box transformation was used again on the transformation observations and the parameter value was obtained as 1. In the time-series theory, when the λ parameter value is 1 in the Cox-Box transformation, there is no more need to apply the transformation to the observations. It should be noted that the λ parameter value is estimated by the maximum likelihood method in Cox-Box transformation. The non-stationarity of the variance was therefore removed. In addition, the Dickey-Fuller test was used to examine the time series stationarity. The results of this test con rmed the stationarity of the time series of the transformed observations. Subsequently, the aforementioned graphs were again drawn for the transformed observations, and it was clearly observed that the minor trend in structural decomposition was removed and stationarity had occurred in the mean and variance. According to the ACF and PACF graphs, the proposed model in the ARIMA family is ARIMA (0,1,1) for the transformed observations. To enable evaluation and comparison with models close to ARIMA (0,1,1), the ARIMA (0,1,2) and ARIMA (1,1,1) models were also tted to the initial time-series observations, and the AIC model selection criterion was used for the excellence of these models. According to the AIC criterion, a model with a lower AIC is better. The results presented in Table 2 suggest that the ARIMA (0,1,2) model provides a better t. The nal model of the observations was therefore considered to be ARIMA (0,1,2) ( Table 2).
The residuals corresponding to this model were also evaluated for their dependence and normality, and the ACF, PACF and QQ-Plot graphs for the residuals showed uncorrelated and normal observation. The Ljung-Box test was also used to con rm the validity of the residual in this section, and its results con rmed that the residuals were not correlated. In the nal section, forecasts were made for the next 36 months using the ARIMA (0,1,2) model. The same con dence intervals applied to all the months, except for the rst month. The graph of this forecast has also been presented until 2019 using the ARIMA (0,1,1) model. Model tting and forecasting of the intended time series data were carried out using Time Series and Forecast Package in R software. In addition, some analyses were performed in MINITAB software. The nal data extracted were analyzed in R software, Minitab and SPSS-23.

Inclusion criteria
A search was carried out of the codes for drug-related deaths by cause of poisoning based on the ICD-10. The leading cause of drug-related death was drug overdose, including acute alcohol intoxication with ethanol and methanol, poisoning as a result of opium overdose, tranquilizer overdose (barbiturates and benzodiazepines) and stimulant overdose (cocaine and amphetamines) [14].

3-Results
The number of deaths varied from 2840 in 2014 and 2810 in 2015 to 3233 in 2016. The total death toll in all three years was 8883 according to MOHME and LMO. Figure 1 shows the time-series graph of the death toll.
According to Fig. 1, the series is almost stationary in its mean. There is also a slight trend conceived in the series. There is no seasonal or periodic trend for the time series.
In order to more precisely examine the time series components, including trend, seasonal effect and random component, the time series is decomposed using the moving average method. The time series decomposition graphs are plotted as in Fig. 2.
Structural time series decomposition shows a slight increasing trend in the time series. An increasing trend can affect the assumption of stationarity in the mean of time series. A rst-order differential operator was therefore used for stationarity in the mean of time series. In addition, stationarity in the variance of these observations also had to be examined. Cox-Box transformation was used for this purpose. To determine the λ parameter for Cox-Box transformation, the observations were transformed and λ = 2.27 was obtained (Fig. 3). It is important to note that, in order to apply the intended transformations for variance and mean stationarity, the required transformations for variance are rst applied, and then the differential operator is used for stationarity in the mean. Applying the noted transformations yielded the time-series graph for the transformed observations as shown in Fig. 4. Figure 5 shows that the trend has been eliminated and the time series of the transformed observations is stationary in the mean and variance. Cox-Box transformation was used again for a more precise examination and λ = 1 was obtained, which means that variance has become stabilized and there is no more need for a transformation. Figure 6 shows the structural time series decomposition of the transformed data under Cox-Box and the rst-order differential operator.
After the Cox-Box transformation and the rst-order differential operator, the data trend was eliminated and the time series of the transformed data became stationary. The Dickey-Fuller test was used to further investigate the stationarity of the transformed data time series. Table 1 shows the Dickey-Fuller test results (Table 1). 8 present the ACF and PACF graphs of the transformed data, respectively.
As can be observed, the ACF graph has exceeded the mean in the rst observations. MA (1) and MA (2) models are therefore suitable for t to the data.
According to the graph, PACF did not exceed the mean for any of the delays. In general, Figs. 7 and 8 show the suitability of the ARIMA (0,1,1) and ARIMA (0,1,2) models for t to the data.

a. Model selection
According to the results presented in Figs. 7 and 8, the ARIMA (0,1,1), ARIMA (0,1,2) and ARIMA (1,1,1) models appear suitable for t to the data. The Akaike coe cients were compared to select the best model. Table 2 presents the proposed models for t to the data. A comparison of the Akaike coe cients of the proposed models for t to the data shows that the ARIMA (0,1,2) model has the lowest Akaike coe cient and is best for t to the data. The model parameters were estimated using the maximum likelihood method, and in order to test the goodness of t of the proposed model, the residuals were analyzed intuitively using the Ljung-Box test. Figure 9 presents the results of the model residual analysis.
As shown in Fig. 9, in the ACF graph, after the initial delay of the ACF value, the other residuals do not exceed the mean value; also, the PACF graph con rms that the correlation is zero for all the delays. In addition to these two graphs, the signi cance levels associated with the Ljung-Box statistic are all above 0.05, which means that the residuals are uncorrelated for all the delays. The proposed tted model is therefore valid for checking the residuals.

b. Forecasting
The forecast death toll in all the months of 2017 to 2019 was estimated by the model as 245.8501 and Table 3 presents the interval estimation of the months in each year with a probability of 95% (Table 3).

4-Discussion
ARIMA time-series models generate data for the future, which, given the current data conditions, are useful for making e cient and useful decisions for the future, because these periodic autoregressive models have the most accurate forecasting power [15].
The results of this study con rm that, during the three years of this survey, mortality rates increased slightly among illicit drug users in Iran. According to the WHO, during the same period, the number of drug-related deaths in the US and UK also increased, which could be attributed to people's increased access to drugs and the ine ciency of harm reduction approaches. Iran's criminal justice system has taken many steps after the Islamic Revolution to combat drug use, but none of them have managed to reduce this phenomenon [14]. Since 2005, when the National Association of Medical Examiners passed the bill that all drug-and alcohol-related deaths should be autopsied by a forensic pathologist, the development and spread of anatomical and toxicology research led to the discovery of a large number of drug-related deaths [15 ].
Therefore, in recent years, the Iranian government has focused on abuse harm reduction approaches to at least prevent harm to the users and the society. Iran has the most expanded level of harm reduction infrastructure in the Middle East, which re ects the good access to methadone maintenance treatment [16]. Although methadone is a viable alternative to heroin and reduces harm to the users, surveys show that many clients of these centers continue to use other illicit drugs, such as methamphetamine, despite having access to methadone. The main reasons for this tendency include compulsory drug treatment programs, social stigma, police encounters, di culties in obtaining governmental identi cation documents and other problems, which socially marginalize the clients of harm reduction centers [16]. The Forensic Toxicology Center has also identi ed methadone as a major cause of death in Iran [17]. It therefore appears that the application of methadone maintenance treatment has not been effective, given the increased death toll in recent years, and experience shows the superiority of consequence prevention over simple focus on quitting [14].
First, a wide range of measures can be taken for the general health of at-risk groups by physicians as the rst target group. Physicians should be adequately educated about the risks and side-effects of abused drugs and be equipped with multiple doses of Naloxone. These efforts should, however, be part of a comprehensive plan to reduce illicit drug supply and expand access to pharmaceutical treatments for drug use [18]. This program can help reduce the number of referrals due to drug abuse in hospitals and its continuation in Iran.
The second goal of the program is to protect young adults' and adolescents' health as a correct indicator of health equity [19]. Due to the global shift of the burden of poverty-related epidemiological diseases to non-communicable diseases [20], the burden of health risks due to diseases has also changed in Iranian young adults and adolescents due to the undeniable role of drug abuse 21 . The analysis of health indicators suggests the need for intervention priorities for the age group of 10-19 years due to the decreasing age of onset of smoking and the increasing trend of use of new drugs among children and adolescents in Iran [19,22]. There is therefore a paramount need to teach harm reduction policies in schools and such measures might indeed prove helpful [14,21]. One of the most successful programs in this area is providing skill training to parents and their participation in school programs for preventing their adolescents' rst time of drug use [23].
To accomplish the third objective, the Iranian criminal law should be examined. Although the criminalization of drug use has accelerated since the Islamic Revolution [14], it has had no effect on reducing drug use [16]. Iran executed 509 drug tra ckers in 2011, triggering a human rights outcry and ending its cooperation with the EU on reducing drug use [8]. According to UN documents, if drug dependence is diagnosed as a complex multidimensional health-disrupting problem 24 , and if, by this view, the behavior of all drug users is decriminalized, or if health and social care is provided to these individuals before executing criminal justice, access to these people will be facilitated and therapists can freely implement harm reduction programs [14,25]. For the rst time, in 2001, Portugal implemented a policy of decriminalization of all drugs, and subsequently gained the second lowest rate of drug-related deaths compared to the world average among EU countries [26]. This project has also been successfully implemented in many Asian countries and has led to lower rates of HIV among drug users. Nevertheless, these countries are now faced with a wave of new drugs, such as methamphetamine and amphetamine stimulants (ATS), which lack alternative drugs, and only long-term psychological interventions can help reduce their harm, which could undermine the public opinion about these policies [24].

5-Limitations
There were two limitations in this study.
-Although a time series model requires short intervals for analysis [27], the overall number of years surveyed in this study was small, and the number of points required for ARIMA analysis thus decreased.
-This study was the rst to examine the trend of drug-related deaths (i.e. by opioids, stimulants and alcohol) in Iran in 2014-2016 and forecast the death toll until 2019. Since there is no dedicated national record in Iran for registering drug-related deaths or drug poisoning, and since different organizations and sources report different statistics, there are many undercounts in this regard. To address this problem, this study used statistics from two major death registries in Iran to correct the error caused by undercounts. One of the major limitations of the present study was the Iranian Legal Medicine Organization's lack of cooperation in providing accurate statistics on the type of poisoning substance. Understanding the pattern and type of poisoning substance in each region of Iran can help design suicide prevention strategies and reduce the risks of accidental poisoning.

6-Conclusion
This study showed that the trend of drug-related deaths rose during the study period, and the modeling process for forecasting showed that this trend will continue until 2019 if proper interventions are not implemented. The discussed policies are therefore recommended to be adopted; otherwise, larger expenditures will be required for combating addiction and suicide or their combination.

7-Declarations
Ethics approval and consent to participate: We did not have participants for this study, so obtaining consent was not necessary. However, all activities performed in this study was in accordance with the ethical standards at institutional and/or The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.
Permission to collect data: A letter was sent for permission to collect data by Mazandaran University of Medical Sciences to the research settings (LMO and MOHME) and they issued an authorization letter for the researcher to collect the data.