Unveiling the future: Wavelet- ARIMAX analysis of climate and diarrhea dynamics in Bangladesh’s Urban centers | BMC Public Health

Study area

The study area focuses divisions in Bangladesh, namely Dhaka, Chattogram, Rajshahi and Khulna due to variation in climatic conditions, population density. Dhaka is the capital of Bangladesh. Another study location is Chattogram division. Chattogram is one of the largest cities in the country, and it is the Port City of Bangladesh, the next study locations are Khulna and Rajshahi, which are the second and 3rd largest among eight divisions in Bangladesh (see Fig. 1).

Data source

Our study collected data for different divisions for two types of data: meteorological and virological (diarrhea). From the Directorate of General Health in Bangladesh^{Footnote 1} (DGHs), we collected this diarrhea dataset. There are 704 daily entries for each division consists of daily recorded data of diarrheal cases, spanning January 31, 2022, through December 12, 2023. Due to limitations on the data availability we only have 704 entries. Data on temperatures (maximum and minimum), humidity, and rainfall between January 31, 2022 and December 12, 2023 are acquired from NASA’s Power^{Footnote 2} and the Bangladesh Meteorological Department^{Footnote 3}.

Data preprocessing

There were few missing values in diarrhea cases. We address these missing values via linear interpolation methods. Linear interpolation is a regularly used method in numerous domains for guessing values between two known data points. In the context of time series data, linear interpolation acts as a basic yet efficient technique for filling in missing values or providing a smoother representation of the data. In the domain of geo-statistics, linear interpolation has been compared with other approaches like quantile kriging, showing its role as a linear estimator for interpolating variables of interest to unknown locations [44].

The formula of Linear Interpolation is:

In this formula, we have the following terms: \(\:{x}_{1}\) and \(\:{y}_{1}\) are the first coordinates, \(\:{x}_{2}\) and \(\:{y}_{2}\) are the second coordinates, x is the place where we execute the interpolation, y is the interpolated value.

Stationarity test

Time series data frequently show seasonality; however, some may not. When dealing with non-stationary data, accurate forecasting becomes difficult and requires a check for stationarity. In this work, we employed the ADF test to evaluate the stationarity of a time series and establish the confidence level for further analysis and forecasting. The ADF test is referred to as a “unit root test,” which is the right way for determining time series stationarity. Furthermore, in the context of the Augmented Dickey-Fuller (ADF) test, the null hypothesis (\(\:{H}_{0}\)) argues for the presence of a unit root in the time series, implying non-stationarity. The alternative hypothesis (\(\:{H}_{1}\)) claims that the time series is stationary.

Wavelet transformation

Wavelet transform is a mathematical tool that decomposes a signal into different frequency components, allowing for a multi-resolution study of the data. Wavelet has advantages in modeling seasonal and nonlinear climate variations, which are some of the deficiencies in traditional models, such as ARIMA and SARIMA. Incorporating wavelets into the model allows for the decomposition of climate data and representation of different frequency components, which enables the modeling of irregular climate changes such as fluctuations in monsoon seasons quite accurately. There are various types of wavelet functions are employed in wavelet transformations, such as the Haar wavelet, the Daubechies wavelet, coifflets, symlets, the biorthogonal wavelet, and the discrete Meyer wavelet [45]. In this study, we use Daubechies wavelet functions with level 4, due to the best compromise between localization and smoothness. These wavelets serve a significant role in the decomposition process, enabling the extraction of both short-term fluctuations and long-term trends in the data. Wavelet transforms have been intensively explored for their applications in several domains, including signal processing, picture compression, and feature extraction [46]. The numerous types of wavelets employed in wavelet transformations give versatility in analyzing and extracting information from signals, making the wavelet transform an important approach in current data analysis and processing. Wavelets effortlessly handle stationary or non-stationary time series data and deliver better outcomes for future forecasting. For this reason, we apply wavelet decomposition to extract the high- and low-frequency signals. There are two types of wavelet transformation, namely continuous wavelet transformation (CWT) and discrete wavelet transformation (DWT). Wavelets are defined by the wavelet function \(\:\psi\:\left(t\right)\) (mother wavelet) and scaling function φ(t) (father wavelet) in the time domain. Wavelets are normalized, finite, short-duration, zero mean functions.

The wavelet function, typically represented as \(\:\psi\:\left(t\right)\), acts as the “Mother wavelet” from which “Child wavelets” are formed by dilation and translation operations. These Child wavelets are then used to evaluate the input signal and produce wavelet coefficients, which offer information about the signal’s frequency content at multiple scales and places. The wavelets created by scaling and translating the mother wavelet are known as “children” wavelets. They are indicated as \(\:{\psi\:}_{k,s}\left(t\right)\:\), where \(\:k\) indicates the translation factor, and s represents the scaling factor. These parameters control how the Wavelet is positioned and scaled relative to the original (mother) wavelet.

Mathematically, given a mother wavelet \(\:\psi\:\left(t\right)\), the offspring Wavelet is defined as \(\:{\psi\:}_{k,s}\left(t\right)\).

Wavelet Transform \(\:{W}_{f}\:\)of a function \(\:f\left(t\right)\:\) is computed by taking the inner product of \(\:f\left(t\right)\) with the translated and dilated versions of the mother wavelet.

This equation is known as the continuous wavelet transformation (CWT) general form. Where a is the scaling parameter, b is the shift or translation parameter, and \(\:{\psi\:}^{*}\) is the complex conjugate function.

ARIMAX model

The ARIMA model is a frequently-used technique in time series analysis and forecasting. ARIMA stands for Autoregressive Integrated Moving Average. It is a mathematical model that helps capture trends and patterns in time-series data, making it valuable for forecasting future values. By incorporating autoregressive, moving average, and differencing components, the ARIMA model provides a flexible framework for studying and forecasting time series data. Moreover, ARIMA does not account for exogenous variables, which may alter the time series data. In contrast, ARIMAX combines these external aspects into its forecasting process to generate more accurate and comprehensive forecasts. By including important exogenous variables into the model, ARIMAX enables greater accuracy in forecasting future values of time series data while also allowing for a better understanding of their interaction with exogenous elements. Additionally, it can aid in spotting probable outliers within the dataset for increased forecasting precision. Our dataset comprises these external parameters for forecasting the occurrence of diarrhea sickness. We selected the ARIMAX model due to this rationale and the inclusion of exogenous factors in our dataset. And, we used autoarima algorithm for selection the best lags for training our model, which identifies optimal lag structures based on model performance and AIC criteria. Therefore, we advocate adopting the ARIMAX model instead of the usual ARIMA model for our investigation.

The Autoregressive (AR) model is stated as following Eq. (1):

Here \(\:p=\text{1,2},3,\dots\:\dots\:\).

The Moving Average (MA) model is written as following Eq. (2):

Here \(\:q=\text{1,2},3,\dots\:.\)

Exogenous Term as follows Eq. (3):

Here \(\:k=\text{1,2},3,\dots\:.\)

The mathematical equation for the ARIMAX model can be written as Eq. (4):

\(\:{Y}_{t}\) indicates the value of the time series at time \(\:t\) ( \(\:t=\text{1,2},3,\dots\:.\)). \(\:\alpha\:\) is a constant or intercept term. \(\:{\gamma\:}_{1},\:{\:\:\gamma\:}_{2},\:\:\:{\gamma\:}_{3}\:\) are the autoregressive parameters, representing the coefficients of the lagged values of the time series and \(\:{\theta\:}_{1\:},\:\:{\theta\:}_{2\:},\:\:\:{\theta\:}_{3}\) are the moving average parameters, representing the coefficients of the past forecast errors. Here, \(\:{X}_{1,t}\)\(\:{X}_{2,t}\)\(\:{X}_{3,t}\) are exogenous variables at a time, and \(\:{\beta\:}_{1}\)\(\:{\beta\:}_{2}\)\(\:{\beta\:}_{3}\) are the coefficients associated with the exogenous variables. p is the order of the autoregressive (AR) component. q is the order of the moving average (MA) component. k is the number of exogenous variables.

Model evaluation

After finishing our model training, we evaluate model performance using multiple techniques: mean absolute error (MAE), Root Mean Squared Logarithmic Error (RMSLE), Mean squared error (MSE), and Root mean squared error (RMSE).

Mean absolute error

The Mean Absolute Error (MAE) is derived as the average of the absolute discrepancies between the forecasted values \(\:{\widehat{Y}}_{t}\) and the actual observed values \(\:{Y}_{t}\), Mathematically denoted as follows:

Mean squared error

Mean Squared Error (MSE) is calculated as the average of the squared deviations between the forecasted values \(\:{\widehat{Y}}_{t}\) and the actual observed values \(\:{Y}_{t}\), Mathematically expressed as follows:

Root mean squared error

Root Mean Squared Error (RMSE) is calculated as the square root of the average of the squared discrepancies between the forecasted values \(\:{\widehat{Y}}_{t}\) and the actual observed values \(\:{Y}_{t}\:\), Mathematically stated as follows:

Root mean squared logarithmic error

RMSLE stands for Root Mean Squared Logarithmic Error. It is the square root of the mean squared logarithmic error (MSLE). It is another popular statistic used to evaluate model performance, particularly in regression assignments where the target variable spans multiple orders of magnitude. Mathematically, the Root Mean Squared Logarithmic Error (RMSLE) is calculated as: