Data

Our empirical analysis was developed using a panel of data from the Spanish provinces from 2002 to 2017, the year when the latest data on provincial GDP were published. For cigarette consumption, we used the official annual tobacco sales and the average price of a pack of 20 cigarettes in euros, as published by the Commission for the Trade of Tobacco. The real gross domestic product (GDP) is available from the National Institute of Statistics in Spain. All series employed here are per capita (for those aged ≥18 years) and expressed in real terms using the consumer price index (CPI base 2016). The descriptive statistics of the dataset are shown in Table 1.

Empirical methodology

Data-driven anomaly detection systems have been discussed in the literature as distortion detection systems in many fields of application^20-23. Such systems aim to detect any abnormal deviations from the normal observations in any given data set. Therefore, these methodologies provide a good opportunity to detect anomalies in tobacco sales. Furthermore, given its abovementioned characteristics, Spain seems to be a reasonable candidate country in which to quantify anomalies.

The aim of this work is to detect tobacco sales anomaly at the provincial level. The prediction of the upper and lower limits of tobacco sales at the provincial level is proposed as a methodology to identify any abnormal deviations in tobacco sales behavior. The proposed methodology is a supervised learning method, which adjusts the model of the relationship between tobacco sales (the dependent variable) and price and GDP (the independent variables). On the other hand, the detection and estimation of anomalies is performed through an unsupervised method, as mentioned before, by means of the computation of upper and lower intervals. Several statistical and machine learning models were compared to find the best model for predicting the tobacco sales in each province (these methods/models are presented in this section).

The main methodology consists of splitting the data into training and test sets for all the available provinces in the Spanish territory, where the training set consists of all the provinces available with the exception of the province to be predicted, which is in the test set. In other words, all provinces are used to predict the sales in a given province without including that province. As is commonly used to explain tobacco consumption behavior in Spain^24-26, the dependent variable is per capita tobacco sales for each province, and the independent variables are per capita GDP and price:

To model the relationship between the dependent variable and the independent variables (the characteristic vector x), two supervised learning methods have been used. In addition, to estimate the upper and lower limits of the prediction interval, quantile predictions are used to generate intervals following the methods explained in this section.

The first method used to model the relationship between the variables is the quantile regression (QR) method. This method was introduced²⁷ for the estimation of models in which the quantiles of the response variable are modeled to depend on the independent variables. The τth quantile for a population is the sample where the 100/τ percent proportion of the population lies. This model of the relationships between different quantile predictors and the dependent variable, in this case, provides clear interpretability for the anomaly detection results, as it is possible to identify an anomaly within a given range²⁸.

The conditional α-quantile q of a scalar variable Y is P(Y≤q|I)=α, where the probability 0<α<1 is given and I denotes an information set generated by independent variables X. A complete justification of this method is given elsewhere²⁷.

For the purposes of this work, two models were combined to build the intervals for detecting anomalies. That is, the conditional 0.1-quantile was set as a lower bound and the 0.9-quantile as an upper bound for each province. By construction, the probability that a value belongs to the interval between the upper and lower bounds is:

In contrast to ordinary least squares, which estimates the conditional mean, this method is based primarily on choosing a model for estimating the conditional quantile. Depending on the strength of the assumptions imposed, a range of parametric or nonparametric options are available²⁹.

For assessing the models, the conditional median response (the 0.5 quantile) for each province was estimated. Not only are the models evaluated for the quality of their punctual predictions, but the ability of their intervals to select the best model with good performance in both tasks is also evaluated, which is discussed later in this section.

In this work, a bagging method is proposed to estimate the conditional quantiles. A combination of random forest (RF) and QR was proposed³⁰, resulting in the quantile regression forest (QRF) approach. One of the main differences between RF and QRF is that in QRF, each node of each tree maintains the values of all observations in that node, but RF only maintains the mean of the observations found in the node³⁰. Ranger is a high-speed implementation of RF or recursive partitioning^31,32 and is particularly well-suited for high-dimensional data^a.

To detect anomalies, two methods were selected to build the prediction intervals (PIs) through conditional quantiles; for every new observation of the response variable, there is a high probability that it lies within the PI²⁴. Furthermore, an anomaly detection and quantification system is proposed that uses an upper bound and lower bound computed using the fitted models.

As mentioned before, the PIs are computed using the calculation of the chosen conditional 0.1-quantile and 0.9-quantile for the lower bound and upper bound, respectively. One requirement for the choice of α for the intervals is the use of a symmetric range (i.e. one cannot use the 0.1 quantile as the lower bound and the 0.7 quantile as the upper). It is not of interest for this study to find the best intervals within a model but to provide a methodology for computing the share of abnormalities, as shown in this section.

The proposed method assumes a uniform distribution with endpoints at the lower and upper limits of the computed PIs. Every point outside this interval range is considered abnormal, and the intervals are also used to quantify the share of abnormalities in the response variable. Specifically, the observed response for province p is abnormal if either case is true:

where η (0 ≤ t ≤ 100%) represents the chosen quantile level given that the interval is symmetrical, and the limits y(100-t) percent and y(t) percent represent the upper and lower conditional quantiles, respectively. The choice of a smaller value of t will lead to a larger number of provinces being predicted to be abnormal.

To train the model, a data partition was performed, as explained in the above sections, and the predictive accuracy of the models was measured by splitting the data into training and test sets.

The error assessment was performed either by using a 0.5-quantile prediction for the quantile analyses or the interval prediction was used to determine the magnitude of the abnormality of tobacco sales that were evaluated with the results of surveys, and certain metrics were used to assess the quality of these intervals.

The performance of the predicted responses (ŷ_i) in relation to the observed responses (y_i) of the training and test sets were assessed by computing the error metrics given in the Supplementary file.

In addition to evaluating the punctual predictions, the estimated prediction interval has also been evaluated in this work. The academic literature has placed a special emphasis on point prediction relative to interval predictions and predictive densities; consequently, there has been little work on the evaluation of PIs³³.

A review of the methods for evaluating point predictions used in this work is presented in the Supplementary file, where the selected metrics assess

the accuracy of the training and test sets. However, as the main idea in forecasting is to decrease uncertainty, an interval prediction evaluation is performed as well. Supplementary file Table 3 summarizes the metrics used to evaluate the prediction interval.

Although anomaly detection has been used in many previous works in other disciplines, the novelty of this method is that it complements abnormality detection with abnormality quantification. As explained previously, the computed intervals are used to quantify the share of abnormalities, and Table 2 shows the formulas that are applied to carry out this task.

As the final product of this work, the model should be able to discern between a province with abnormal tobacco sales and a province without abnormal tobacco sales but also, when an abnormality is detected, it should also be able to determine whether this abnormality is due to a quantity below the lower bound or to a quantity above the upper bound.

National temporal evolution of the anomalies detected

As indicated, to quantify the anomalies in the Spanish provinces, we use the upper anomaly ratio (UAR) and the lower anomaly ratio (LAR). The average UAR and LAR for the Spanish territories are calculated by averaging the ratio of tobacco sales per capita observed to be below the lower limit (the lower prediction bound) and the ratio of tobacco sales per capita observed to be above the upper bound (the upper prediction bound). Figure 2 shows the aforementioned index, and the answers to two important questions can be observed: 1) the magnitude of the average upper anomaly exceeds 40%, while the average lower anomaly is <15%; and 2) the lower anomalies are trending upward over time and the upper anomalies are trending downward.

Figure 2

Average UAR and LAR for the Spanish territory

Regarding the magnitude of the LAR and UAR, although, as will be seen in this section, geographical patterns play a key role, it seems that there is a significant difference between the indices. The average LAR represents the average percentage of those provinces that have observed sales below the estimated values. Therefore, given that this index is conditioned by cross-border and illegal trade, it seems that both activities have remained constant during the period studied. However, the average UAR increased notably until 2005 and has decreased since then. The UAR index represents the average anomaly in the provinces in which the observed sales exceed those estimated by the model. For this reason, both the effect of tourism and that of the cross-border trade with countries where tobacco is more expensive than in Spain determine the magnitude of the UAR index. The provinces that border France have been considered to exhibit a cross-border effect, while the rest of the provinces are highly influenced by tourism. Supplementary file Figure 1 shows the part of the average UAR that represents cross-border effects and the part that represents tourism effects. Both effects seem to be decreasing, with the trend in the cross-border effect being stronger. This last highlight is consistent with recent evidence on cross-border transactions between Spain and France¹³.

Regional temporal evolution of the anomalies detected

Once the magnitudes have been estimated at the national level, it is interesting to analyze the temporal evolution of the UAR and LAR in the Spanish provinces. First, regarding the LAR, there are 6 provinces that stand out for their behavior. As shown in Figure 3, there are three provinces in southern Spain (Sevilla, Cádiz and Córdoba) in which the LAR is estimated for the first time in 2010 and grows notably until 2017, reaching values close to 40% in some cases. On the other hand, in three other provinces (Orense, Pontevedra and Lugo), the LAR is decreasing; in addition, it rarely takes values close to 20%. In addition, Figure 3 also shows the values given by the EPSs that were conducted by the TTCs from 2012 to 2017. As the figure shows, in line with the previous literature, it seems that the use of independent data provides estimates of the illicit market that are lower than those provided by the EPSs. Specifically, in Cádiz, the EPS is slightly overestimated, and in Córdoba, there are substantial differences.

Figure 3

Temporal evolution of UAR and LAR in the Spanish provinces

Six provinces also stand out in the temporal evolution of the UAR in the Spanish territory. On the one hand, Malaga, Alicante and Islas Baleares provinces that are highly influenced by tourism^a present a similar trend. On the other hand, in the provinces that border France, such as Gerona, Huesca and Guipúzcoa, the UAR exhibits a decreasing trend, with the anomalies detected in Gerona being of a much higher magnitude. As indicated in the introduction, Gerona is the Spanish region in which sales fell the most due to border closures related to the public health crisis of COVID-19. Therefore, these results are consistent with what the raw data indicated.

Regional distribution of the anomalies detected

Although the provincial evolution of the UAR and LAR provides important information, the geographical distribution of the anomalies helps us to understand the contagion effect of the UAR and the LAR, as well as tobacco consumption behavior at the borders with other countries. In this sense, Figure 4 shows the geographical distribution of the LAR and UAR according to the QR model, respectively. As can be seen, the anomalies in the provinces in which lower-than-estimated sales have been observed reach 35% and are concentrated in the Northwest in 2002 and in the South in 2017, which is consistent with the previous literature. In addition, regarding the UAR, the anomalies in the provinces with sales above those estimated reach values of 190%. Finally, while these high UAR values were concentrated in tourist provinces and provinces bordering Portugal in 2002, in 2017 the concentration in tourist provinces remained the same, but it is now the border with France where the anomalies are located.

Figure 4

Geographical distribution of LAR and UAR in the Spanish provinces (QR model)

One of the main contributions of the present study is the quantification of the LAR and UAR. Along these lines, the closure of Spain’s borders with other countries due to the health crisis caused by COVID-19 during the months of April and May 2020 has made it possible to analyze the robustness of the results of this work. As shown in the Supplementary file Figures 2 and 3, both the geographical pattern and the magnitude of the UAR estimated in this work are robust to the decreases in tobacco sales in Spain in the months of April and May 2020. On the one hand, tobacco sales fell in April and May by up to 180% and 160%, respectively. This magnitude is consistent with the anomalies shown in Figure 4. Furthermore, the geographical patterns are also consistent, with the greatest decreases in tobacco sales being observed in the areas near the border with France and in the tourist provinces. In addition, the areas near the border with Gibraltar are those in which tobacco sales have decreased the least during the border closures. It is precisely in these provinces that the highest UARs are observed in the latest years analyzed.

Limitations

The results shown are not without limitations. On the one hand, cigarette sales have been used. Although cigarettes account for almost 90% of total tobacco consumption in Spain, part of the anomalies detected may be motivated by the consumption of substitute products. In addition, the use of aggregated data prevents knowing individual behaviors. By using individual data, it could be known whether individual sociodemographic characteristics are influencing the abnormalities detected. Another line of future research consists of separating which part of the detected anomalies come from illegal activities. Finally, the used method assumes that tobacco control measures have the same compliance rate in all the provinces studied.