Meteorological and human mobility data on predicting COVID-19 cases by a novel hybrid decomposition method with anomaly detection analysis: A case study in the capitals of Brazil

2.1. Human mobility in the prediction of COVID-19 cases

According to Nayak et al. (2021), one of the primary impacts on predicting the COVID-19 cases consists of the variations in engagement, i.e. how committed people are to taking measures to reduce the number of COVID-19 cases. These measures include washing hands, wearing face masks, and maintaining social distancing. Concerning social distancing, one way to estimate the level of commitment is by analyzing the rates of human mobility. In line with this, Oztig and Askin (2020) considered the flow of people at airports as a human mobility measure, and observed that the greater the number of airports in a country, the more likely it is for the country to have a higher number of COVID-19 cases. This conclusion was drawn by the use of negative binomial regression analysis.

Badr et al. (2020) studied the correlation between social distancing and COVID-19 cases, where social distancing was quantified by mobility patterns. To model the mobility data, the authors considered changes in commuting patterns between and within counties in the USA. The data to model the mobility patterns were obtained by Teralytics (Zürich, Switzerland). Wang et al. (2020) coupled the data of confirmed COVID-19 cases with the Google mobility data in Australia. The authors concluded that the social restriction policies imposed in the country at the emergence of the first COVID-19 case were effective in curbing the spread of the virus. Moreover, they observed that the correlation between human mobility and the spread of COVID-19 varies according to the type of mobility.

Shao et al. (2021) used human mobility data from 47 countries in 6 continents collected from Mobility Trends Reports (from Apple Inc.), and showed that human mobility is strongly related to the COVID-19 transmission rate. Zhu et al. (2020) demonstrated a positive link between human mobility and the number of people infected by COVID-19, considering data from 120 cities in China. These studies show a clear influence of human mobility on the spread of COVID-19. However, cities present different human mobility patterns depending on factors such as how technological the cities are, the conditions of public and private transportation systems, among others. Therefore, the relationship/correlation between human mobility and dissemination of COVID-19 must be evaluated considering the cities’ particularities. Bearing this in mind, in this study we focus on the relationship between human mobility and the spread of COVID-19 cases in each of the 27 Brazilian capitals.

2.2. Meteorological variables in the prediction of COVID-19 cases

Sahin, 2020, Sharma and Gupta, 2021 state that meteorological features should be used to improve the accuracy of COVID-19 predictions. Such variables are crucial factors affecting infectious diseases, whether in terms of changes in the transmission dynamics, regarding host susceptibility, or the survival of the virus in the environment (McClymont & Hu, 2021). In line with this, Gupta et al. (2020) studied the relationship among new COVID-19 cases, absolute humidity, and temperatures in the USA. The authors observed that the spread of COVID-19 was majorly influenced by the absolute humidity in a narrow range of 4 to 6 g/m³.

Islam et al. (2020) investigated the link between some environmental factors and COVID-19 cases in 206 countries/regions (until April 20, 2020). The relationship between the spread of COVID-19 and humidity, and UV index were inconclusive. Their investigation suggested a negative relationship between wind speed and COVID-19 cases. Moreover, a higher rate of COVID-19 cases was observed in environments with an absolute humidity between 5 and 10 g/m³. In Singapore, Pani et al. (2020) revealed that temperature, absolute humidity, and dew point have a positive correlation with the number of daily COVID-19 cases. Wind speed, atmospheric boundary layer height, and ventilation coefficient, on the other hand, showed a negative correlation with the number of COVID-19 cases. In Turkey, Sahin (2020) showed that wind speed has a positive correlation with COVID-19 cases, and that temperature and COVID-19 cases are negatively correlated.

In Brazil, Neto and Melo (2020) concluded that only the average air humidity was significantly correlated with the number of COVID-19 cases (considering data from Brazilian capitals, and data available from April 2020 to May 2020). The study revealed a positive correlation, in contrast with the results obtained by others studies performed in cities in China, Spain, and the United States. The authors also demonstrated that population density presented a strong positive correlation with the number of COVID-19 cases in the Brazilian capitals. They emphasize that population density, which is linked with higher human mobility, and poorer social-economic environments that have deficient sanitary conditions contribute to the spread of the virus.

Although some studies have addressed the studies involving the relationship between meteorological variables and COVID-19, some results appear to be inconsistent. On the one hand, temperature and humidity, for example, were reported as having a significant impact in the majority of the studies. On the other hand, the correlation was positive in some cases and negative in others. These observations suggest that the link between meteorological features and the number of COVID-19 cases is complex and hard to generalize. There is evidence that meteorological variables contribute to the increase in the transmission of COVID-19, but the effect of these relationships should be studied locally, since other factors such as human mobility and public health measures (lockdown, for example) also have a strong influence on the number of COVID-19 cases.

2.3. COVID-19 and forecasting

From a methodological point of view, several studies attempt to understand the spread of COVID-19 using artificial intelligence. Albahri et al. (2020) provided an exhaustive overview of integrated artificial intelligence based on data mining and machine learning algorithms. The authors pointed to a need for integrated sensor technologies for outdoor scenarios to control the spread of the coronavirus. This process is only possible when there is an interconnection with IoT technologies. Nayak et al. (2021) present an overview of the applicability of intelligent systems such as machine learning and deep learning to solve COVID-19 outbreak-related issues. Sharma and Gupta (2021) reported and summarized the research performed on COVID-19 with machine learning and big data.

The literature presents few studies that address the problem of predicting new COVID-19 cases through decomposition methods. To our knowledge, only Silva et al., 2020, Mousavi et al., 2020 used decomposition methods in their predictions considering independent variables. Both proposed strategies were based on the variational mode decomposition (VMD). Silva et al. (2020) used VMD and some prediction techniques including deep learning and machine learning, to predict COVID-19 cumulative confirmed cases in five Brazilian states and five American states with high daily incidences. The authors used temperature and precipitation as exogenous variables. They pointed that the VMD coupled with cubist regression achieved the best results among the tested techniques. Mousavi et al. (2020) proposed a model based on the combination of VMD with Long Short Term Memory considering the daily temperature, humidity, and transmission rates in the prediction of new COVID-19 daily cases in Maharashtra, Tamil Nadu, and Gujarat, India. Among these works, only Mousavi et al. (2020) addressed the prediction of daily COVID-19 cases, since such prediction is more difficult because of the accumulated cases. In this study, we address the prediction of new daily cases of COVID-19.

3.1. Data

COVID-19 data were obtained from Brasil.io (Justen et al., 2020), which compiles newsletters from the State Health Secretariats of Brazil. Meteorological data were obtained from the Centro de Previsño de Tempo e Estudos Climáticos located at the Instituto Nacional de Pesquisas Espaciais (CPTEC, 2020). The meteorological data considered in this study are:

• •

  Minimum Temperature (Min Temp): refers to the daily minimum temperature in degrees Celsius;

• •

  Maximum Temperature (Max Temp): refers to the daily maximum temperature in degrees Celsius;

• •

  Humidity (Hum): refers to the daily air humidity in percentage;

• •

  Rainfall (Rain): refers to the daily total precipitation in millimeters.

Human mobility data were obtained from the COVID-19 Community Mobility Reports (GOOGLE, 2020) prepared by Google. These reports point to geographical movement trends over time, across different categories of places. The place categories are:

• •

  Retail and recreation (RR): refers to mobility trends to places like restaurants, shopping centers, theme parks, etc;

• •

  Grocery and pharmacy (GP): refers to mobility trends to places like grocery markets, farmers markets, pharmacies, etc;

• •

  Parks (PA): refers to mobility trends to places like local parks, public beaches, public gardens, etc;

• •

  Transit stations (TS): refers to mobility trends to places like subway, bus, train stations, etc;

• •

  Workplaces (WO): refers to mobility to places of work;

• •

  Residential (RE): refers to mobility to places of residence.

The Residential category shows a change in the permanence of people in their homes, while the other categories measure changes in the total number of visitors. Changes in mobility patterns each day were compared with a baseline corresponding to the same day of the week. This baseline corresponds to the median of the corresponding day of the week, during the five weeks from January 3 to February 6, 2020.

The number of observations in human mobility data and meteorological data varies according to the number of data in the variable relative to daily COVID-19 cases in each city. In each city, all reported data start at the day they confirmed the first COVID-19 case in the city (column “first case” in Table 1 of the Supplementary Material) and end at the final compiled day: November 6, 2020. A more descriptive analysis of the data is presented in Section 1 of the Supplementary Material.

3.2. Ensemble Empirical Mode Decomposition

Time series decomposition techniques have the goal of extracting simple periodic signals from the original time series, which can be used as inputs to machine learning approaches or other statistical models. Our study focuses on the use of the Ensemble Empirical Mode Decomposition (EEMD) technique (Huang & Wu, 2008), an adaptive data analysis method based on local characteristics of the data. EEMD catches nonlinear, non-stationary oscillations effectively. EEMD has been successfully used in several types of datasets (Lin et al., 2021, Niu et al., 2019), mainly in meteorological data, such as temperature (Liu et al., 2019), precipitation (Alizadeh et al., 2019), and wind speed (Santhosh et al., 2019).

EEMD is an improvement of the empirical mode decomposition (EMD) method (Huang et al., 1998, Huang and Wu, 2008). It aims at decomposing the original data into a series of modes, called finite intrinsic mode functions (IMFs) and a residual, identifying the oscillatory modes that coexist. EEMD overcomes the so-called mode-mixing problem found in EMD. The mode-mixing occurs when different oscillation components coexist in a single IMF and very similar oscillations reside in different IMFs (Huang & Wu, 2008).

EEMD uses an ensemble of IMFs obtained by applying EMD to several different series of the original time series obtained by adding white Gaussian noise. Adding a white Gaussian noise reduces the mode-mixing problem by occupying the whole time–frequency space (Huang & Wu, 2008). In summary, EEMD has the following steps.

• 1.

  Let $W_t,m$ and s be the input data corresponding to, respectively, the original time series that will be decomposed, the number of ensembles, and the number of IMFs to be extracted from $W_t$.

• 2.

  Make $k=1$, a control variable that indicates the ensemble to be generated in the iteration.

• 3.

  Generate a new time series $Z_t$, obtained from $W_t$ for the ensemble k, adding to it a white noise with a standard deviation ${\sigma }_{\mathit{noise}}$ proportional to the standard deviation of $W_t$, called ${\sigma }_{\mathit{original}}$. Therefore, ${\sigma }_{\mathit{noise}}=\mu {\sigma }_{\mathit{original}}$, where $\mu$ is a relatively small number which must be empirically determined.

• 4.

  Make $j=1$, a control variable related to the index of the IMF of the k-th ensemble to be defined in the following steps, referred to as ${\text{IMF}}_j^k$.

• 5.

  Identify all the local extreme values of $Z_t$ – a combination of high and low values of the series. After that, interpolate all this values by a cubic spline interpolation as the upper (high values) and lower envelopes (low values), respectively $e_{\mathit{max}}^k$ and $e_{\mathit{min}}^k$.

• 6.

  Calculate the point-to-point arithmetic mean between the envelopes – $m_t^k=(e_{\mathit{min}}^k+e_{\mathit{max}}^k)/2$ – and subtract this “average time series” from time series $Z_t$, obtaining the time series $d_t^k$ – $d_t^k=Z_t-m_t^k$.

• 7.

  If $j⩽s$, then ${\text{IMF}}_j^k=d_t^k,j=j+1,Z_t=Z_t-d_t^k$, and repeat steps 5 and 6. If $j>s$, assign $Z_t$ to the residual time series, called ${\text{Res}}_{W_t}^{}$.

• 8.

  Make $k=k+1$ and repeat Steps 3 to 7 until $k>m$, i.e. until the method obtains the m ensembles.

The values of $\mu$ and m were empirically chosen after several computational tests. These tests indicated that an ensemble number $m=125$ and the $\mu$ value equals $0.01$ presented better outcomes. Furthermore, because of the proposed EEMD-ARIMAX method in Section 3.4, the number of IMFs into which the time series is decomposed was fixed in advance and, after tests, we found that a decomposition into 5 IMFs plus a residual was the most appropriate, i.e., $s=5$.

Fig. 1 shows the IMFs extracted from the data of São Paulo COVID-19 cases, by applying the EEMD algorithm. The IMFs were plotted from the first to the last component extracted from the series, where the last plot corresponds to the residual. The x-axis indicates the days, whereas the y-axis represents the values of the decomposed time series.

3.3. Autoregressive Integrated Moving Average Exogenous inputs (ARIMAX)

The Autoregressive Integrated Moving Average (ARIMA) model proposed by Box and Jenkins (1990) is the most general class of models for forecasting time series due to its simplicity of application and capability of handling non-stationary data. The AR part of ARIMA indicates that the variable of interest is regressed on its own lagged values. The MA part indicates that the regression error is a linear combination of error values that occurred in the past. Finally, the I (for “integrated”) part represents the order of differencing to turn the time series into a stationary series (if necessary). Differencing means replacing the original series by the difference between their values and the previous values (Box & Jenkins, 1990). The ARIMA model that includes other time series as input variables (exogenous variables) is referred to as Autoregressive Integrated Moving Average Exogenous inputs (ARIMAX) model.

The parameters of ARIMAX($p,d,q,n$) model are: p, the number of autoregressive terms; d, the number of nonseasonal differences needed for stationary; q, the number of lagged forecast errors in the prediction equation; n, the number of exogenous variables; $\eta$, a constant; and, ${\varphi }_i$, for $i=1,\dots ,p,{\theta }_j$, for $j=1,\dots q$, and ${\zeta }_l$, for $l=1,\dots ,n$, the model parameters. Mathematically, this model can be formulated as in Eq. (1).

where $W_t$ and $W_{t-i}$, for $i=1,\dots ,p$, are the predicted values of the time series; $Y_l$, for $l=1,\dots ,n$, are the exogenous variables; and $e_{t-j}$, for $j=1,\dots ,q$, represent the error terms.

3.4. EEMD-ARIMAX

To our knowledge, EEMD has not yet been used to predict time series with independent variables. The main idea behind EEMD-ARIMAX is to predict time series of independent and dependent variables. For this, we first decompose each time series of the independent variables ($Y_1,Y_2,\dots ,Y_n$) and dependent variables by applying the EEMD method, creating s levels of decomposition for each variable. Then, in each level of the decomposition, we use the ARIMAX method to predict the IMFs related to the dependent variables, by considering the IMFs of the variables $Y_1,Y_2,\dots ,Y_n$ as the exogenous variables. We employ the same procedure to predict the time series of the residual values. Finally, by summing the predicted time series, we obtain the prediction for the original time series of the daily number of COVID-19 cases. The algorithm of the proposed EEMD-ARIMAX method can be described by steps 1–5:

• 1.

  Let $X_t$ be dependent variable under study, $Y_1,Y_2,\dots ,Y_n$ the independent/predictor variables, m the number of ensembles, and s the number of IMFs that will be extracted of each time series;

• 2.

  Apply EEMD to decompose the time series of the dependent and independent variables individually, to obtain a set of s IMFs and a time series Res, in each decomposition;

• 3.

  Fit IMFs of the same “levels” using ARIMAX – meaning that the j-th IMF of the time series represented by the “Daily number of COVID-19 cases”, denoted here by ${\text{IMF}}_{X_t}^j$, will be fitted by the j-th IMFs of the same level j of the time series related to meteorological/mobility variables $Y_l$, denoted by ${\text{IMF}}_{X_l}^j$, for all $l=1,\dots ,n$. The estimated j-th IMF is denoted by ${\widehat{\mathrm{IMF}}}^j$;

• 4.

  Denote the residual values obtained by applying EEMD in $Y_1,\dots Y_n$ by ${{\text{Res}}_{Y_1}}_{}$,...,${{\text{Res}}_Y}_n$, respectively. Denote the residual value found by applying EEMD in $X_t$ by ${{\text{Res}}_X}_t$. Let $\widehat{\mathrm{Res}}$ be the estimated time series of Res${{}_X}_t$ through ARIMAX using ${{\text{Res}}_Y}_1$,...,${{\text{Res}}_Y}_n$ as exogenous variables;

• 5.

  Denote the fitted values of variable $X_t$ by ${\widehat{X}}_t$. Thereby, ${\widehat{X}}_t={\widehat{\mathrm{IMF}}}^1+\dots +{\widehat{\mathrm{IMF}}}^s+\widehat{\mathrm{Res}}$.

A flowchart of the proposed EEMD-ARIMAX method is presented in Fig. 2 .

4.1. Correlation analysis

We evaluate the pairwise correlation between the number of COVID-19 cases and meteorological/mobility variables using Spearman correlation. For more details about this measure, see Section 2.1 of the Supplementary Material.

Tables 1 and 2 show the correlation values – columns “$\rho$” – between the number of COVID-19 cases and the meteorological and human mobility variables, respectively, for all Brazilian capitals, in the period considered. In addition, these tables present the p-value regarding the statistical significance of the corresponding variables at a significance level of $\alpha =0.01$. Therefore, if the p-value of the indicated correlation is less than or equal to $0.01$, the correlation is said to be statistically significant.

Table 1.

Spearman correlation between the number of COVID-19 cases and the meteorological data.

Open in a new tab

Table 2.

Spearman correlation between the number of COVID-19 cases and the human mobility data.

Open in a new tab

We consider that two variables are correlated if $\rho ⩾0.3$ or $\rho ⩽-0.3$. As stated before, if $\rho$ is positive, the variables are directly proportional, otherwise, they are inversely proportional. Therefore, on the one hand, we say that there is a positive correlation between a pair of variables when $\rho ⩾0.3$, meaning that there is evidence that the variables grow together. On the other, when the correlation is negative, i.e., $\rho ⩽-0.3$, it means that the analyzed pair of variables has an opposite behavior: the greater the values of one variable, the smaller the values of the other variable. For better visualization, we highlighted the positive correlations in dark gray, and negative correlations in light gray.

According to the results, the number of COVID-19 cases and meteorological variables were correlated in 16 cities. In 11 of them, the correlated meteorological variable was the minimum temperature. The number of COVID-19 cases and meteorological variables were not correlated in any of the cities in the South region. The maximum temperature and the number of COVID-19 cases were correlated in all cities in the Midwest region.

The correlations between the number of COVID-19 cases and minimum temperature were negative, indicating that the number of cases increases when the minimum temperature decreases. The same behavior was observed between the number of COVID-19 cases and maximum temperature, except in Teresina-PI and Palmas-TO. In these two cities, the relationship between the number of COVID-19 cases and maximum temperature was inversely proportional.

Humidity and the number of daily COVID-19 cases are correlated in the following cities: Palmas-TO, Fortaleza-CE, Recife-PE, Teresina-PI, Brasília-DF, and Goiânia-GO. Particularly in Palmas-TO and Teresina-PI, the humidity and number of COVID-19 cases showed a strong correlation. The average humidity of Palmas-TO was the lowest among the capitals of the North region. The average humidity of Teresina-PI was the second lowest average of the capitals of the Northeast region. Since humidity is directly linked to temperature, these facts could explain the inversely proportional correlations between the number of COVID-19 cases and the maximum temperature in both cities.

Among the 6 capitals that showed a correlation between the number of COVID-19 cases and humidity, Fortaleza-CE and Recife-PE presented correlations of 0.329 and 0.316 respectively. The other four capitals showed negative correlations. One can observe that the rainfall variable and the number of COVID-19 cases are not correlated in Fortaleza-CE and Recife-PE. In Palmas-TO, Teresina-PI, Brasília-DF, and Goiânia-GO, on the other hand, it is possible to see that they were negatively correlated.

It is known that meteorological data regarding temperature, humidity, and rainfall are related and, therefore, influence one another. In this study, however, we will only consider the meteorological variables of each capital that had a correlation with the number of COVID-19 cases greater than 0.3, in absolute value. These values are summarized in Table 3 .

As mentioned before, Table 2 shows the correlation between the mobility variables and the number of COVID-19 cases. The mobility variables and the corresponding cities with which the number of COVID-19 cases have a positive correlation are:

• •

  Retail and recreation: Boa Vista-RR, Palmas-TO, Teresina-PI, Campo Grande-MS, Goiânia-GO, Florianópolis-SC, Porto Alegre-RS;

• •

  Grocery and pharmacy: Boa Vista-RR, Palmas-TO, Teresina-PI, Brasília-DF, Campo Grande-MS, Goiânia-GO, Vitória-ES, Florianópolis-SC, Porto Alegre-RS;

• •

  Parks: Palmas-TO, Teresina-PI, Campo Grande-MS, Goiânia-GO, Florianópolis-SC;

• •

  Transit stations: Boa Vista-RR, Palmas-TO, Teresina-PI, Campo Grande-MS, Goiânia-GO, Florianópolis-SC, Porto Alegre-RS;

• •

  Workplaces: Boa Vista-RR, Palmas-TO, Teresina-PI, Campo Grande-MS, Goiânia-GO, Belo Horizonte-MG, Vitória-ES, Florianópolis-SC, Porto Alegre-RS;

• •

  Residential: Fortaleza-CE, Maceió-AL, Recife-PE, São Luis-MA, Cuiabá-MT.

On the one hand, the positive correlation between the mobility parameters and the number of COVID-19 cases, except for the Residential variable, shows that the increase in the number of COVID-19 cases is directly proportional to the rise in the populations’ mobility trends in traffic, pharmacies, work, parks, and retail. This means that the higher the mobility rate, the greater the number of cases. On the other hand, some cities showed a negative correlation between mobility variables and the number of COVID-19 cases. They are:

• •

  Retail and recreation: Fortaleza-CE, Cuiabá-MT;

• •

  Grocery and pharmacy: Cuiabá-MT;

• •

  Parks: Fortaleza-CE, Recife-PE, Cuiabá-MT;

• •

  Transit stations: Fortaleza-CE, São Luis-MA, Cuiabá-MT;

• •

  Workplaces: Cuiabá-MT;

• •

  Residential: Palmas-TO, Teresina-PI, Goiânia-GO, Florianópolis-SC, Porto Alegre-RS.

Therefore, for example, the negative correlation between Residential and daily COVID-19 cases means that the fewer people stayed at home, the greater the number of COVID-19 cases.

The negative correlations between COVID-19 cases and meteorological parameters in Cuiabá-MT are due to a sequence of null values at the end of the series describing the number of COVID-19 cases. If these values were excluded, the correlation coefficient between these variables would be positive.

4.2. Analysis of the number of predicted cases

The EEMD-ARIMAX method was implemented in the R software using the “Rlibeemd” and “forecast”. We generated 125 new time series for each variable considering that the standard deviation of Gaussian noise was 1% of the standard deviation of the corresponding original time series.

Table 4 shows the results of the EEMD-ARIMAX method for all Brazilian capitals. We compared EEMD-ARIMAX with the ARIMAX method. The objective was to analyze the effect of EEMD on the prediction. In both methods and for each city, we present the widely employed mean error (ME), root-mean-square deviation (RMSE), and mean absolute error (MAE) measures to describe the results of the predictions. For details about these measures, see Section 2.1 of the Supplementary Material. Column “City-Federative unit” shows the pair city-federative unit and the parameters used by ARIMAX to forecast the number of COVID-19 cases in this capital. These parameters were calibrated for each city using auto.arima() function in R. The independent variables that were considered to predict the number of cases of COVID-19 in each corresponding city are shown in Table 3 and follow the Spearman correlation coefficients shown in Table 1, Table 2.

In all cities, the proposed decomposition method improved the predictions of the time series in terms of RMSE values. The average RMSE of the predictions considering only the ARIMAX method was 211.987 with a standard deviation of 186.335. Using the EEMD-ARIMAX method, the average RMSE was 155.330 with a standard deviation of 145.645. EEMD-ARIMAX showed an improvement of 26.73% over ARIMAX. Section 2.2 of the Supplementary material presents some graphics comparing the original time series with the predicted values by EEMD-ARIMAX in all Brazilian regions.

5.1. Analyzing model errors

We analyzed the days for which the model significantly missed the predicted number of cases for each city. The error made by the model was quantified by the difference between the observed and predicted number of cases, as shown in Eq. (2).

where $c_i^t$ and ${\widehat{c}}_i^t$ are, respectively, the observed and predicted number of COVID-19 cases on day t in city i. An error is considered significant when $e_i^t>\mathit{TD}\left(E_i\right)$, where $E_i$ is the vector formed by the elements $e_i^t$ $\forall t$, and TD is defined by Eq. (3), where $\mathit{Mean}\left(E_i\right)$ and $\mathit{STD}\left(E_i\right)$ are the arithmetic mean and standard deviation of $E_i$, respectively.

Fig. 3 illustrates the values of $e_i^t$ considering the city of Goiânia - GO. The threshold value $\mathit{TD}\left(E_i\right)$ is highlighted in red. Therefore, every day t whose $e_i^t$ is above the red line corresponds to a significantly mispredicted day by the model.

5.2. Analyzing data anomalies

A spectral anomaly detection strategy was adopted to detect days when the recorded number of daily COVID-19 cases was potentially anomalous. While the model errors are identified by comparing the predicted values with the observed values, the anomaly detection strategy analyzes the daily variation in the number of cases considering the distance between cities to identify potentially anomalous variations.

For example, if a city has a slight variation over two days in the number of COVID-19 cases, we expect nearby cities to have a similar variation. Similarly, if the number of cases in a city suffers a significant increase from one day to the other, we expect nearby cities also to have a relative increase in the number of cases.

To perform this analysis, we model a complete and weighted network where a node $v_i\in V$ represents a city, and the weight of the edge $w_{\mathit{ij}}$ is the Euclidean distance between cities i and j. Each node $v_i$ carries the daily variation in the number of cases in city i, with the daily variation $s_i^t$ defined by Eq. (4), $\forall t>1$.

A signal $S^t$ contains the values $s_i^t$ for every city i in the dataset. We calculate the spectra ${\widehat{S}}^t$ of $S^t$ signal, $\forall t>1$, using the Fourier transform for graphs (Sandryhaila & Moura, 2013). In graph Fourier analysis, the graph Laplacian eigenvectors associated with small eigenvalues ${\lambda }_l$ vary slowly across the graph, whereas eigenvectors associated with larger eigenvalues oscillate more rapidly (Ortega, Frossard, Kovacevic, Moura, & Vandergheynst, 2018). It means that if two vertices are connected by an edge with a large weight, the values of the eigenvector at those locations are likely to be similar. This concept is then used to define low and high frequencies for signals indexed by graphs.

According to this definition, abrupt oscillations are concentrated at the high frequencies of the signal spectrum. To highlight abrupt variations and expose anomalies, we accentuate the magnitude of the high frequencies of ${\widehat{S}}^t$ spectrum to make anomalous variations more evident, generating a new spectrum ${\widehat{R}}^t$. We apply the inverse Fourier transform to ${\widehat{R}}^t$ to get a new $R^t$ signal, that contains the accentuated variation in the number of cases in each city. The intuition behind this operation is that, if the variation in the number of cases in a city i is normal, then $r_i^t<s_i^t$ probably holds, where $r_i^t$ is the i-th element in vector $R^t$. On the other hand, if the city has an anomalous variation, $r_i^t>s_i^t$ probably holds. Fig. 4 shows a graphic visualization of the normal variation and the accentuated variation.

The threshold used to determine whether a variation is anomalous or not is calculated in the same way for errors, as defined in Eq. (3), for $R^t$. Fig. 5 illustrates the values of $R_i$, which is a vector with the $r_i^t$ of a given city i, and the threshold $\mathit{TD}\left(R_i\right)$. It is worth noting the similarity between Fig. 3, Fig. 5, which points out that there is a direct relation between the cases in which EEMD-ARIMAX significantly missed the prediction and the days when the attenuator pointed out potentially anomalous variations.

5.3. Comparing errors and anomalies

As presented in Section 5.2, Fig. 3, Fig. 5 indicate that there is a direct relationship between the days when EEMD-ARIMAX made a significant error and the days whose variation in the number of cases was interpreted as potentially anomalous. We compared the model’s errors with anomalous variations to establish a quantifier that indicates whether there is, in fact, a direct relationship between them.

For each city, two sets were defined: set CE, containing the days on which the model made a significant error; and set CA, with the days whose variation was detected as potentially anomalous. To quantify the relationship between CE and CA, we adopted the following criterion: if a day $t\in \mathit{CE},1<t<\mathit{nc}$ and $\{t-1,t,t+1\}\cap \mathit{CA}=\varnothing$, then the error made by the model on day t and the anomalous variation that occurred on the days adjacent to t are directly related.

Fig. 6 shows the percentage of days the model made significant mispredictions and which are directly related to a day with a potentially anomalous variation. On average, more than 60% of the days when the model was wrong were detected as anomalous, as indicated by the red line, which represents the average.

This result indicates that EEMD-ARIMAX was affected by errors in the data and the results point that the model’s errors are directly related to anomalous variations. To overcome this problem and correct the anomalies, we adopted a spectral strategy for removing anomalies, also based on the Fourier transform. Section 5.4 presents the methodology employed to correct the data anomalies.

5.4. Normalizing Data

As discussed before, the abrupt variations are concentrated in the high frequencies. To detect the anomalies, the presented anomaly detection strategy accentuated the magnitude of the high frequencies. Then, to correct the anomalies, we adopted a strategy that does the opposite, using a low-pass filter that attenuates high frequencies. Unlike the high-frequency accentuator, the low pass filter decreases the magnitude of the high frequencies, attenuating abrupt variations.

While the accentuator was applied to the signal that carried the daily variation in the number of cases, the low-pass filter was applied to the signal formed by the number of cases on each day, that is, the C signal, where $c_i^t$ is the number of cases in city i at day t, to generate a filtered signal $\widehat{C}$. Fig. 7 compares an original signal and a filtered signal. It is possible to note that, in general, the signal oscillation is mitigated, ensuring a more reliable signal.

By applying both the ARIMAX and the EEMD-ARIMAX methods to the normalized data, we obtained the results shown in Table 5 , which are presented as in Table 4. The average RMSE for the forecasting considering only the ARIMAX method was 142.981 with a standard deviation of 122.703. The average RMSE for the EEMD-ARIMAX method was 107.664 with a standard deviation of 99.917. Therefore, EEMD-ARIMAX was 24.70% better than ARIMAX.

There was an improvement of 30.69% in the prediction by EEMD-ARIMAX when normalized data were used. Fig. 8 shows all the RMSEs obtained by the EEMD-ARIMAX method using non-normalized (black) and normalized (red) data.