# Capturing the Data Uncertainty Change in the Cocaine Consumption in Spain Using an Epidemiologically Based Model

Authors: Christopher Anaya; Clara Burgos; Juan-Carlos Cortés; Rafael-J. Villanueva
Journal: Abstract and Applied Analysis
 (2016)

---

## Abstract

A probabilistic model is proposed to study the transmission dynamics of the cocaine consumption in Spain during the period of 1995–2011. Using the so-called probabilistic fitting technique, we study if the model is able to capture the data uncertainty coming from surveys. The proposed model is formulated through a nonlinear system of difference equations whose coefficients are treated as stochastic processes. A discussion regarding the usefulness and limitations of probabilistic fitting technique in order to capture the data uncertainty of the proposed model is presented.


---

## Body

## 1. Introduction

In [1], the authors presented an epidemiologically based model to study the transmission dynamics of the cocaine consumption in Spain. This model worked well in the presented scenario, 1995–2005, where the population of cocaine consumers had an increasing trend. However, in 2008, the economic crisis began and the Spanish Health Ministry started to apply the National Drug Plan. As a consequence, the nonconsumer population increased and the consumer populations started to decrease. This sudden change could not be captured by the model as it was shown in [2].In this work, we propose to use a variation of the model stated in [1] and assume that the model parameters are variable over the time in order to see if the variable model parameters are able to adapt and help the model to capture the observed changes in the cocaine consumption.With this aim, we use the probabilistic fitting technique recently developed in [3, 4] by some of the authors. With this technique, we expect to capture the data uncertainty (data coming from surveys), despite the change of trend, and provide an estimation to the probability distribution of the model parameters. The latter will help us study if the observed data variation also produces variations in some or all model parameters as well as their quantification.These types of mathematical models also have been used in the study of epidemics and other drug addictions, such as alcohol, tobacco, ecstasy, or heroin addiction [5–9], and in the approach to other sociological topics that are spread by social contact as obesity or extreme ideological behaviors [10, 11].The work is organized as follows. In Section2, we present the data we are going to use about cocaine consumers in Spain. In Section 3, we build a type-epidemiological model describing the cocaine consumption dynamics in Spain. Also, the model is scaled in order to match data and model magnitudes. In Section 4, we recall and adapt theprobabilistic fitting technique introduced in [3]. This technique is applied to the model in order to obtain both data estimation with uncertainty and an estimation of the probability distribution of the model parameters. In Section 5, we present and discuss the results. In Section 6, conclusions are drawn.


## 2. Data

In this study, we use data from the Survey on Alcohol and Drugs in Spain (EDADES), part of the Spanish National Drug Plan [12]. This survey is published every two years. Specifically, we have focused on the following key question quoted in this survey:how often have you consumed cocaine? The available responses were as follows: “never consumed (nonconsumers),” “at least once in your life (occasional consumers),” “at least once in the past year (regular consumers),” and “at least once in the past 30 days (habitual consumers).” This survey was launched in 1995 and we have data until 2011. The data are collected in Table 1.Table 1
Survey results.% of people who belong to nonconsumers, occasional consumers, regular consumers, or habitual consumers of cocaine in Spain from 1995 to 2011.
Survey dates Nonconsumer Occasionaconsumer Regular consumer Habitual consumer t 1 = 1995 0.94 4 0.03 4 0.018 0.004 t 2 = 1997 0.948 0.032 0.015 0.005 t 3 = 1999 0.948 0.031 0.015 0.006 t 4 = 2001 0.911 0.049 0.026 0.014 t 5 = 2003 0.903 0.059 0.027 0.011 t 6 = 2005 0.884 0.070 0.030 0.016 t 7 = 2007 0.874 0.080 0.030 0.016 t 8 = 2009 0.860 0.102 0.026 0.012 t 9 = 2011 0.879 0.088 0.022 0.011In order to simplify the notation, we are going to define the set of time instants where survey data are available:(1)Ω=t1=1995,t2=1997,t3=1999,t4=2001,t5=2003,t6=2005,t7=2007,t8=2009,t9=2011.


## 3. Model Building

From the available survey data, we then introduce an epidemiological model to describe the dynamics of the four subpopulations over the time. While a traditional epidemiological study considers the transmission of a biological contagion, for our study, we will apply asocial contagion approach. The underlying assumption is the following: we can consider cocaine consumption to be a contagion insofar as the probability that one individual becomes a cocaine consumer depends on the interaction with people who already exhibit this behavior, that is, fellow cocaine consumers [1, 13]. Let us assumehomogeneous population mixing; that is, each individual can contact with any other individual [14, 15]. With this assumption, we can now build a system of difference equations based on the model introduced in [1] to further describe the transmission dynamics.As with the survey results, we divide the Spanish 15–65-year-old population (the age group considered in the surveys) into four subpopulations:(1)
Nt, the amount of people who have never consumed cocaine at year t(2)
Ot, the amount of people who occasionally consume cocaine (at least once in their lives) at year t(3)
Rt, the amount of people who regularly consume cocaine (at least once in the past year) at year t(4)
Ht, the amount of people who habitually consume cocaine (at least once in the past month) at year tFurthermore, we have the total population(2)Tt=Nt+Ot+Rt+Ht,which is not going to be a constant value over the time.The transitions between the different subpopulations are determined as follows:(i)
Let us consider that the newly recruited 15-year-old individuals become members ofNt subpopulation at rate μ; that is, we consider that they have never consumed cocaine before.(ii)
An individual in the nonconsumer population can become an occasional consumer through contact with occasional, regular, or habitual consumers. Then, this social contagion term can be modeled withβtNt/TtOt+Rt+Ht.(iii)
An occasional consumer can increase their consumption and become a regular consumer. This transition term can be modeled byγtOt.(iv)
A regular consumer can increase their consumption and become a habitual consumer. We assume that this can be modeled by the termσtRt.(v)
A habitual consumer can decide to go to drug-treatment therapy. If he/she stays in therapy for 6 months, according to experts, he/she can be considered a nonconsumer [16]. This term is modeled by ϵtHt.Using the above assumptions, a dynamic cocaine consumption model for the 15–65-year-old Spanish population is given by the following nonlinear system of difference equations:(3)Nt+1=Nt+μTt-dNNt-βtNtTtOt+Rt+Ht+ϵtHt,(4)Ot+1=Ot-dOOt+βtNtTtOt+Rt+Ht-γtOt,(5)Rt+1=Rt-dRRt+γtOt-σtRt,(6)Ht+1=Ht-dHHt+σtRt-ϵtHt,where(i)
μ is the birth rate (the mortality between 0 years old and 14 years old is so small that it can be discarded),(ii)
dN is the death rate for nonconsumers,(iii)
dO is the death rate for occasional consumers,(iv)
dR is the death rate for regular consumers,(v)
dH is the death rate for habitual consumers,(vi)
βt is the time-dependent transmission rate for cocaine consumption,(vii)
γt is the time-dependent transition rate between occasional and regular cocaine consumers,(viii)
σt is the time-dependent transition rate between regular and habitual cocaine consumers,(ix)
ϵt is the time-dependent rate at which habitual consumers enter and complete drug-treatment therapy.Figure1 shows the flow diagram of the dynamic cocaine consumption model. The boxes represent the subpopulations and the arrows represent the transitions between the subpopulations. Arrows are labeled by their corresponding parameters of the model.Figure 1
Flow diagram of the mathematical model (3)–(6) for the dynamics of cocaine consumption in Spain. The boxes represent the subpopulations and the arrows represent the transitions between the subpopulations. Arrows are labeled by the parameters of the model.Following the idea proposed in [1], we are going to determine bounds for the model parameters. The birth rate μ∈0.009168052,0.011275930; both ends of this interval represent the minimum and maximum values for the birth rate during the period 1995–2011 [17]. Analogously, the death rate dN∈0.008166102,0.009228115. Taking into account the fact that the death rate of the consumers may be up to 6.8% higher than dN [12] and the fact that the more an individual consumes drugs the higher probability to die because of drug abuse is, we have that dN<dO<dR<dH≤0.009228115×1.068.Also, we have the time-dependent parametersβt, γt, σt, and ϵt for t∈Ω for which we consider that all of them lie between 0 and the double of the values given in [1] (maximum likelihood). Then, βt∈0,2×0.09614, γt∈0,2×0.0596, σt∈[0,2×0.0579], and ϵt∈0,2×0.0000456, for t∈Ω.The time step for simulations will be of two years according to the time instants inΩ. Then, the total number of model parameters is 1 birth rate, 4 death rates, 9βt, 9γt, 9σt, and 9ϵt, that is, 41 model parameters.


### 3.1. Model Scaling

Despite the fact that (3)–(6) describe the cocaine consumption dynamics in Spain, there is a mismatch between the input for these difference equations, which is formulated for total population counts, and the survey data available, which is stated in percentages. We therefore need to scale the equations. To do this, we add together (3)–(6). Then, taking into account (2), the left-hand side yields(7)Nt+1+Ot+1+Rt+1+Ht+1=Tt+1,while the right-hand side yields(8)Tt+μTt-dNNt-dOOt-dRRt-dHHt.Now, we shall define(9)nt=NtTt,ot=OtTt,rt=RtTt,ht=HtTt,from which we get nt+ot+rt+ht=1. From here, we will scale the equations, starting with (3). Here, using (8)-(9), one obtains(10)nt+1=Nt+1Tt+1=Nt+μTt-dNNt-βtNt/TtOt+Rt+Ht+ϵtHtTt+μTt-dNNt-dOOt-dRRt-dHHt,and by dividing the numerator and denominator by Tt, one gets(11)nt+1=nt+μ-dNnt-βtntot+rt+ht+ϵtht1+μ-dNnt-dOot-dRrt-dHht.Applying the same method to the other equations, the following expressions are obtained:(12)ot+1=Ot+1Tt+1=Ot-dOOt+βtNt/TtOt+Rt+Ht-γtOtTt+μTt-dNNt-dOOt-dRRt-dHHt=ot-dOot+βtntot+rt+ht-γtot1+μ-dNnt-dOot-dRrt-dHht,rt+1=Rt+1Tt+1=Rt-dRRt+γtOt-σtRtTt+μTt-dNNt-dOOt-dRRt-dHHt=rt-dRrt+γtot-σtrt1+μ-dNnt-dOot-dRrt-dHht,ht+1=Ht+1Tt+1=Ht-dHHt+σtRt-ϵtHtTt+μTt-dNNt-dOOt-dRRt-dHHt=ht-dHht+σtrt-ϵtht1+μ-dNnt-dOot-dRrt-dHht.And the scaled equations can be written in the following form: (13)nt+1=nt+μ-dNnt-βtntot+rt+ht+ϵtht1+μ-dNnt-dOot-dRrt-dHht,ot+1=ot-dOot+βtntot+rt+ht-γtot1+μ-dNnt-dOot-dRrt-dHht,rt+1=rt-dRrt+γtot-σtrt1+μ-dNnt-dOot-dRrt-dHht,ht+1=ht-dHht+σtrt-ϵtht1+μ-dNnt-dOot-dRrt-dHht.


## 3.1. Model Scaling

Despite the fact that (3)–(6) describe the cocaine consumption dynamics in Spain, there is a mismatch between the input for these difference equations, which is formulated for total population counts, and the survey data available, which is stated in percentages. We therefore need to scale the equations. To do this, we add together (3)–(6). Then, taking into account (2), the left-hand side yields(7)Nt+1+Ot+1+Rt+1+Ht+1=Tt+1,while the right-hand side yields(8)Tt+μTt-dNNt-dOOt-dRRt-dHHt.Now, we shall define(9)nt=NtTt,ot=OtTt,rt=RtTt,ht=HtTt,from which we get nt+ot+rt+ht=1. From here, we will scale the equations, starting with (3). Here, using (8)-(9), one obtains(10)nt+1=Nt+1Tt+1=Nt+μTt-dNNt-βtNt/TtOt+Rt+Ht+ϵtHtTt+μTt-dNNt-dOOt-dRRt-dHHt,and by dividing the numerator and denominator by Tt, one gets(11)nt+1=nt+μ-dNnt-βtntot+rt+ht+ϵtht1+μ-dNnt-dOot-dRrt-dHht.Applying the same method to the other equations, the following expressions are obtained:(12)ot+1=Ot+1Tt+1=Ot-dOOt+βtNt/TtOt+Rt+Ht-γtOtTt+μTt-dNNt-dOOt-dRRt-dHHt=ot-dOot+βtntot+rt+ht-γtot1+μ-dNnt-dOot-dRrt-dHht,rt+1=Rt+1Tt+1=Rt-dRRt+γtOt-σtRtTt+μTt-dNNt-dOOt-dRRt-dHHt=rt-dRrt+γtot-σtrt1+μ-dNnt-dOot-dRrt-dHht,ht+1=Ht+1Tt+1=Ht-dHHt+σtRt-ϵtHtTt+μTt-dNNt-dOOt-dRRt-dHHt=ht-dHht+σtrt-ϵtht1+μ-dNnt-dOot-dRrt-dHht.And the scaled equations can be written in the following form: (13)nt+1=nt+μ-dNnt-βtntot+rt+ht+ϵtht1+μ-dNnt-dOot-dRrt-dHht,ot+1=ot-dOot+βtntot+rt+ht-γtot1+μ-dNnt-dOot-dRrt-dHht,rt+1=rt-dRrt+γtot-σtrt1+μ-dNnt-dOot-dRrt-dHht,ht+1=ht-dHht+σtrt-ϵtht1+μ-dNnt-dOot-dRrt-dHht.


## 4. Probabilistic Fitting

This technique consists of using information of the data survey to assign probability distributions to the data. Then, we sample data values from these probability distributions and fit the model to the sampled data. Thus, we find model parameters that fit not only the data but also the uncertainty contained into the intrinsic survey error. Thus, these model parameters will allow the model to capture the data uncertainty (with95% confidence intervals).


### 4.1. Data 95% Confidence Intervals (95% CI)

Data in Table1 correspond to the mean percentage obtained from the EDADES survey conducted between 1995 and 2011 every two years. In the technical specifications of each survey, we can see sample sizes of 15000, 15000, 15000, 15000, 15000, 27934, 23715, 20109, and 22128 interviews, respectively.Taking into account the fact that the sample is not the same for each survey, let us assume that the survey outputs are independent. For each one of the 9 available surveys, let us denote byXj=X1j,X2j,X3j,X4j, 0≤Xij≤nj, i=1,2,3,4, and j=1,…,9, a random vector whose entries are X1j=# nonconsumers, X2j=# occasional consumers, X3j=# regular consumers, and X4j=# habitual consumers and n1=15000, n2=15000, n3=15000, n4=15000, n5=15000, n6=27934, n7=23715, n8=20109, and n9=22128 are the sample sizes of surveys. These components represent exclusive selections (events) with probabilities(14)PjX1j=x1=θ1j,PjX2j=x2=θ2j,PjX3j=x3=θ3j,PjX4j=x4=θ4j,j=1,…,9,where θ1j, θ2j, θ3j, and θ4j are the percentages collected in Table 1 for each survey j: j=1,…,9. Thus, each random vector Xj follows a multinomial probability distribution. Therefore, the probability that X1j occurs x1 times, X2j occurs x2 times, X3j occurs x3 times, and X4j occurs x4 times is given by(15)Pnjjx1,x2,x3,x4=nj!x1!x2!x3!x4!θ1jx1θ2jx2θ3jx3θ4jx4,j=1,…,9,where x1+x2+x3+x4=nj. The resulting multinomials for each EDADES survey can be seen in Table 2.Table 2
Date and data joint multinomial probability function of each survey.
Survey dates Joint multinomial probability functions t 1 = 1995 P 1 = P 15000 1 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.94 4 x 1 0.03 4 x 2 0.01 8 x 3 0.00 4 x 4  t 2 = 1997 P 2 = P 15000 2 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.94 8 x 1 0.03 2 x 2 0.01 5 x 3 0.00 5 x 4  t 3 = 1999 P 3 = P 15000 3 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.94 8 x 1 0.03 1 x 2 0.01 5 x 3 0.00 6 x 4  t 4 = 2001 P 4 = P 15000 4 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.91 1 x 1 0.04 9 x 2 0.02 6 x 3 0.01 4 x 4  t 5 = 2003 P 5 = P 15000 5 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.90 3 x 1 0.05 9 x 2 0.02 7 x 3 0.01 1 x 4  t 6 = 2005 P 6 = P 27934 6 x 1 , x 2 , x 3 , x 4 = 27934 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.88 4 x 1 0.0 7 x 2 0.0 3 x 3 0.01 6 x 4  t 7 = 2007 P 7 = P 23715 7 x 1 , x 2 , x 3 , x 4 = 23715 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.87 4 x 1 0.0 8 x 2 0.0 3 x 3 0.01 6 x 4  t 8 = 2009 P 8 = P 20109 8 x 1 , x 2 , x 3 , x 4 = 20109 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.8 6 x 1 0.10 2 x 2 0.02 6 x 3 0.01 2 x 4  t 9 = 2011 P 9 = P 22128 9 x 1 , x 2 , x 3 , x 4 = 22128 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.87 9 x 1 0.08 8 x 2 0.02 2 x 3 0.01 1 x 4Now, we compute the quantiles 2.5 and 97.5 (95% CI) of each one of the joint multinomial distributions in Table2: j=1,2,…,9. The obtained 95% CI are collected in Table 3.Table 3
95% CI of the EDADES surveys data using the joint multinomial probability function of each survey.
Survey dates Nonconsumer Occasional consumer Regular consumer Habitual consumer t 1 = 1995 [ 0.940,0.947 ] [ 0.031,0.037 ] [ 0.016,0.020 ] [ 0.003,0.005 ] t 2 = 1997 [ 0.944,0.952 ] [ 0.029,0.035 ] [ 0.013,0.017 ] [ 0.004,0.006 ] t 3 = 1999 [ 0.944,0.952 ] [ 0.028,0.034 ] [ 0.013,0.017 ] [ 0.005,0.007 ] t 4 = 2001 [ 0.906,0.916 ] [ 0.046,0.053 ] [ 0.023,0.029 ] [ 0.012,0.016 ] t 5 = 2003 [ 0.898,0.908 ] [ 0.055,0.063 ] [ 0.024,0.030 ] [ 0.009,0.013 ] t 6 = 2005 [ 0.880,0.888 ] [ 0.067,0.073 ] [ 0.028,0.032 ] [ 0.015,0.018 ] t 7 = 2007 [ 0.870,0.878 ] [ 0.076,0.083 ] [ 0.028,0.032 ] [ 0.014,0.018 ] t 8 = 2009 [ 0.855,0.865 ] [ 0.098,0.106 ] [ 0.024,0.028 ] [ 0.011,0.014 ] t 9 = 2011 [ 0.875,0.883 ] [ 0.084,0.092 ] [ 0.020,0.024 ] [ 0.010,0.012 ]


### 4.2. Probabilistic Estimation

LetM(t;α) be a short representation of the scaled model (13), where α denotes the list of 41 model parameters; that is,(16)α=μ,dN,dO,dR,dH,βt1,…,βt9,γt1,…,γt9,σt1,…,σt9,ϵt1,…,ϵt9.Also, we have the data95% confidence intervals (95% CI) of Table 3 obtained by sampling the joint probability distributions of Table 2 and calculating the percentiles 2.5 and 97.5.Given the probability distributionsPi, i=1,…,9, in Table 2, we take a sample dtij∗, i=1,…,9 and j=1,…,4, and we look for the values to the parameters α∗ so that(17)Mt1;α∗Mt2;α∗⋮Mt9;α∗-dt11∗,dt12∗,dt13∗,dt14∗dt21∗,dt22∗,dt23∗,dt24∗⋮dt91∗,dt92∗,dt93∗,dt94∗Fis as minimum as possible, with ·F being the Frobenius norm [18]:(18)aijF=∑i∑jaij2.This procedure is a classic optimization problem that can be carried out through genetic algorithms, PSO, Nelder-Mead, and so forth [19, 20] but having the following key difference: now we will fit the data sampled from probability distributions instead of the raw data. We shall perform the fitting N times (N being a large number), storing both the parameter values α∗ and the calculated errors e∗, ordered from smallest to largest errors. The result of this procedure is a list of model parameters fitted to a sample of the data with their corresponding errors represented in Table 4.Table 4
Model fitting toN samples of the data’s probability distributions.
Error Parameters Model e 1 ∗ α 1 M ( t ; α 1 ) e 2 ∗ α 2 M ( t ; α 2 ) ⋮ ⋮ ⋮ e N ∗ α N M ( t ; α N )The valueN should be a large number in order to capture as much data uncertainty as possible during the sampling process and this uncertainty could be fitted by the model. In our case, we take N=25000.Now, we takeMt;α1 and Mt;α2 in Table 4 and calculate the outputs for times t1,…,t9 in Ω, the time instants where data are available. For each time instant, we shall calculate percentiles 2.5 and 97.5, one for each one of the 4 subpopulations. Hence, we will name m2 the sum of the following:(i)
The Frobenius norm of the difference between the percentiles2.5 from the model output and from the data percentiles 2.5 in Table 3(ii)
The Frobenius norm of the difference between the percentiles 97.5 from the model output and from the data percentiles 97.5 in Table3We will repeat the above process with the outputs fromM(t;α1), M(t;α2), and M(t;α3), obtaining m3, the measure between the confidence bands from the outputs and the data. The same shall be done for m4 with M(t;α1), M(t;α2), M(t;α3), and M(t;α4) and so on, until Mt;α1,Mt;α2,…,M(t;αN), obtaining mN as the measure between the confidence bands from the outputs and the data.Takingmk=min⁡m1,…,mN, we ensure that the 95% confidence band of the outputs from Mt;α1,Mt;α2,…,Mt;αk is the closest to the 95% confidence band from the data, with which our model will capture the maximum uncertainty of the data from the output of the models M(t;α1),M(t;α2),…,M(t;αk).


## 4.1. Data 95% Confidence Intervals (95% CI)

Data in Table1 correspond to the mean percentage obtained from the EDADES survey conducted between 1995 and 2011 every two years. In the technical specifications of each survey, we can see sample sizes of 15000, 15000, 15000, 15000, 15000, 27934, 23715, 20109, and 22128 interviews, respectively.Taking into account the fact that the sample is not the same for each survey, let us assume that the survey outputs are independent. For each one of the 9 available surveys, let us denote byXj=X1j,X2j,X3j,X4j, 0≤Xij≤nj, i=1,2,3,4, and j=1,…,9, a random vector whose entries are X1j=# nonconsumers, X2j=# occasional consumers, X3j=# regular consumers, and X4j=# habitual consumers and n1=15000, n2=15000, n3=15000, n4=15000, n5=15000, n6=27934, n7=23715, n8=20109, and n9=22128 are the sample sizes of surveys. These components represent exclusive selections (events) with probabilities(14)PjX1j=x1=θ1j,PjX2j=x2=θ2j,PjX3j=x3=θ3j,PjX4j=x4=θ4j,j=1,…,9,where θ1j, θ2j, θ3j, and θ4j are the percentages collected in Table 1 for each survey j: j=1,…,9. Thus, each random vector Xj follows a multinomial probability distribution. Therefore, the probability that X1j occurs x1 times, X2j occurs x2 times, X3j occurs x3 times, and X4j occurs x4 times is given by(15)Pnjjx1,x2,x3,x4=nj!x1!x2!x3!x4!θ1jx1θ2jx2θ3jx3θ4jx4,j=1,…,9,where x1+x2+x3+x4=nj. The resulting multinomials for each EDADES survey can be seen in Table 2.Table 2
Date and data joint multinomial probability function of each survey.
Survey dates Joint multinomial probability functions t 1 = 1995 P 1 = P 15000 1 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.94 4 x 1 0.03 4 x 2 0.01 8 x 3 0.00 4 x 4  t 2 = 1997 P 2 = P 15000 2 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.94 8 x 1 0.03 2 x 2 0.01 5 x 3 0.00 5 x 4  t 3 = 1999 P 3 = P 15000 3 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.94 8 x 1 0.03 1 x 2 0.01 5 x 3 0.00 6 x 4  t 4 = 2001 P 4 = P 15000 4 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.91 1 x 1 0.04 9 x 2 0.02 6 x 3 0.01 4 x 4  t 5 = 2003 P 5 = P 15000 5 x 1 , x 2 , x 3 , x 4 = 15000 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.90 3 x 1 0.05 9 x 2 0.02 7 x 3 0.01 1 x 4  t 6 = 2005 P 6 = P 27934 6 x 1 , x 2 , x 3 , x 4 = 27934 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.88 4 x 1 0.0 7 x 2 0.0 3 x 3 0.01 6 x 4  t 7 = 2007 P 7 = P 23715 7 x 1 , x 2 , x 3 , x 4 = 23715 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.87 4 x 1 0.0 8 x 2 0.0 3 x 3 0.01 6 x 4  t 8 = 2009 P 8 = P 20109 8 x 1 , x 2 , x 3 , x 4 = 20109 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.8 6 x 1 0.10 2 x 2 0.02 6 x 3 0.01 2 x 4  t 9 = 2011 P 9 = P 22128 9 x 1 , x 2 , x 3 , x 4 = 22128 ! x 1 ! x 2 ! x 3 ! x 4 ! 0.87 9 x 1 0.08 8 x 2 0.02 2 x 3 0.01 1 x 4Now, we compute the quantiles 2.5 and 97.5 (95% CI) of each one of the joint multinomial distributions in Table2: j=1,2,…,9. The obtained 95% CI are collected in Table 3.Table 3
95% CI of the EDADES surveys data using the joint multinomial probability function of each survey.
Survey dates Nonconsumer Occasional consumer Regular consumer Habitual consumer t 1 = 1995 [ 0.940,0.947 ] [ 0.031,0.037 ] [ 0.016,0.020 ] [ 0.003,0.005 ] t 2 = 1997 [ 0.944,0.952 ] [ 0.029,0.035 ] [ 0.013,0.017 ] [ 0.004,0.006 ] t 3 = 1999 [ 0.944,0.952 ] [ 0.028,0.034 ] [ 0.013,0.017 ] [ 0.005,0.007 ] t 4 = 2001 [ 0.906,0.916 ] [ 0.046,0.053 ] [ 0.023,0.029 ] [ 0.012,0.016 ] t 5 = 2003 [ 0.898,0.908 ] [ 0.055,0.063 ] [ 0.024,0.030 ] [ 0.009,0.013 ] t 6 = 2005 [ 0.880,0.888 ] [ 0.067,0.073 ] [ 0.028,0.032 ] [ 0.015,0.018 ] t 7 = 2007 [ 0.870,0.878 ] [ 0.076,0.083 ] [ 0.028,0.032 ] [ 0.014,0.018 ] t 8 = 2009 [ 0.855,0.865 ] [ 0.098,0.106 ] [ 0.024,0.028 ] [ 0.011,0.014 ] t 9 = 2011 [ 0.875,0.883 ] [ 0.084,0.092 ] [ 0.020,0.024 ] [ 0.010,0.012 ]


## 4.2. Probabilistic Estimation

LetM(t;α) be a short representation of the scaled model (13), where α denotes the list of 41 model parameters; that is,(16)α=μ,dN,dO,dR,dH,βt1,…,βt9,γt1,…,γt9,σt1,…,σt9,ϵt1,…,ϵt9.Also, we have the data95% confidence intervals (95% CI) of Table 3 obtained by sampling the joint probability distributions of Table 2 and calculating the percentiles 2.5 and 97.5.Given the probability distributionsPi, i=1,…,9, in Table 2, we take a sample dtij∗, i=1,…,9 and j=1,…,4, and we look for the values to the parameters α∗ so that(17)Mt1;α∗Mt2;α∗⋮Mt9;α∗-dt11∗,dt12∗,dt13∗,dt14∗dt21∗,dt22∗,dt23∗,dt24∗⋮dt91∗,dt92∗,dt93∗,dt94∗Fis as minimum as possible, with ·F being the Frobenius norm [18]:(18)aijF=∑i∑jaij2.This procedure is a classic optimization problem that can be carried out through genetic algorithms, PSO, Nelder-Mead, and so forth [19, 20] but having the following key difference: now we will fit the data sampled from probability distributions instead of the raw data. We shall perform the fitting N times (N being a large number), storing both the parameter values α∗ and the calculated errors e∗, ordered from smallest to largest errors. The result of this procedure is a list of model parameters fitted to a sample of the data with their corresponding errors represented in Table 4.Table 4
Model fitting toN samples of the data’s probability distributions.
Error Parameters Model e 1 ∗ α 1 M ( t ; α 1 ) e 2 ∗ α 2 M ( t ; α 2 ) ⋮ ⋮ ⋮ e N ∗ α N M ( t ; α N )The valueN should be a large number in order to capture as much data uncertainty as possible during the sampling process and this uncertainty could be fitted by the model. In our case, we take N=25000.Now, we takeMt;α1 and Mt;α2 in Table 4 and calculate the outputs for times t1,…,t9 in Ω, the time instants where data are available. For each time instant, we shall calculate percentiles 2.5 and 97.5, one for each one of the 4 subpopulations. Hence, we will name m2 the sum of the following:(i)
The Frobenius norm of the difference between the percentiles2.5 from the model output and from the data percentiles 2.5 in Table 3(ii)
The Frobenius norm of the difference between the percentiles 97.5 from the model output and from the data percentiles 97.5 in Table3We will repeat the above process with the outputs fromM(t;α1), M(t;α2), and M(t;α3), obtaining m3, the measure between the confidence bands from the outputs and the data. The same shall be done for m4 with M(t;α1), M(t;α2), M(t;α3), and M(t;α4) and so on, until Mt;α1,Mt;α2,…,M(t;αN), obtaining mN as the measure between the confidence bands from the outputs and the data.Takingmk=min⁡m1,…,mN, we ensure that the 95% confidence band of the outputs from Mt;α1,Mt;α2,…,Mt;αk is the closest to the 95% confidence band from the data, with which our model will capture the maximum uncertainty of the data from the output of the models M(t;α1),M(t;α2),…,M(t;αk).


## 5. Results

The probabilistic fitting procedure returns, as the best fitting, the one giving fork=40 with an error of mk=0.0502666. A graphical representation of this result can be seen In Figure 2.Figure 2
Result of the probabilistic fitting. The green band corresponds to95% CI model output and the blue line its mean. The red points are the 95% CI data given in Table 3 and the black points their means. We can see how the model tries to capture the data uncertainty, but the goal is not achieved completely.The model95% CI band captures almost all data uncertainty except 1999 and 2011 in nonconsumers, 2009 in occasional consumers, and 2001 and 2011 in habitual consumers (there is no intersection between the model output band and these data 95% CI).The decrease in the nonconsumer subpopulation during the period of 1997–2001 is not properly captured by the model. With regard to the economic crisis and the National Drug Plan proposed in 2008 by the Spanish government, their effect does not appear until 2011, where an increase in nonconsumers and a decrease in occasional consumers arise. As can be seen in Figure2, for these two subpopulations, the model is unable to capture the data uncertainty. Also, in the habitual consumers’ subpopulation, the model still does not realize the trend change.All these comments lead us to say that the proposed model was fairly good to explain what happened when the number of consumers was increasing, mainly in occasional and regular consumers; nevertheless when the trend changed, it could not adapt and explain entirely what was happening.With respect to the model parameters, in Figure3, we can see the estimation of the probability density function (PDF) of βt, γt, σt, and ϵt for every t∈Ω.Figure 3
Estimation of the random processesβt, γt, σt, and ϵt using the estimated PDF of the parameters in the time instants in Ω.Figure3 has been built taking k=40 sets of model parameters selected by the probabilistic fitting procedure, in particular corresponding to βt, γt, σt, and ϵt for t∈Ω, and generate a kernel PDF in each time instant in Ω using theMathematica [21] command KernelMixtureDistribution[].In Figures3 and 4, we can see the variability of the parameters over the time, trying to adapt to the cocaine consumption behavior changes. Note that ϵt has a mean and 95% confidence interval fairly stable, which means that the flow of people going to therapy is almost invariable over the time with independence of the changes in the consumption behavior. Therefore, we could make ϵt=ϵ for all t∈Ω and transform this random process into a random variable.Figure 4
Estimated95% confidence intervals (band) and means of the random processes βt, γt, σt, and ϵt plotted in Figure 3.The use of variable parameters may allow us to identify changes in the studied phenomenon which have not been considered in the hypotheses or which could arise unexpectedly, as it seems to happen in this regard to the economic crisis and the National Drug Plan, in this case. Accordingly, we can see in Figure4 how parameters βt and γt, and with lower changes parameter σt, try to adapt in the rightmost part of the graphs to data changes, but, unfortunately, they do not change enough to capture the data uncertainty when the consumption trend varies. This leads us to think that the model may not be appropriate to explain the decreasing trends in cocaine consumption.


## 6. Conclusion

In this work, we present a model formulated through a system of nonlinear difference equations, where the model parameters are variable over the time. This model intends to describe the cocaine consumption dynamics in Spain in the period of 1995–2011. The model follows fairly well the dynamics until 2009, where the effect of two nonconsidered facts, the effect of the economic crisis and the implantation of the National Drug Plan, introduced a change in the behavior of the individuals, especially in the young people who are the group that is increasing the nonconsumer subpopulation.We expected that the introduction of variable parameters over the time would allow us to capture the changes in the cocaine consumption dynamics; however, it could not. This indicates that the proposed model may not be appropriate when a decreasing trend in cocaine consumption emerges.With regard to the parameters, we must say that it was not necessary to consider the parameterϵt as variable, because it is fairly stable over the time. The remaining parameters, βt, γt, and σt, are trying to adapt to the data trend changes, but, in nonconsumer and habitual consumer subpopulations, they fail to do it accurately.Finally, the three sides that this study, a model with variable parameters using the probabilistic fitting technique, provides should be pointed out:(1)
Consider the model parameters by default as random processes (a family of random variables indexed by time), which is more versatile than considering them as random variables, and then we can detect if they really vary over the time.(2)
Study the estimation of the parameter random processes and obtain estimations for the means, the confidence intervals, the PDFs, and hence other higher statistical moments.(3)
Determine the usefulness of the model in several realistic situations or when data changes because unexpected scenarios emerge.Although it was not mentioned before, a great computational effort has been done to obtain the presented results and this is one of the major drawbacks we want to improve in future works.



---

*Source: 1758459-2016-11-30.xml*