2.1. Multiple linear regression

The multiple linear regression (MLR) model is a simple and widely used modeling technique. The MLR model is given as:

y_j = β₀ + β₁x_1j + β₂x_2j + … + β_px_pj + u_j, j = 1,2, . . . , n (1)

where y_j is the dependent (output) variable, β_j are the regression coefficients, x_j are the independent (input) variable, n is the number of observations, p is the number of independent variables, and u_j is the residual term.

2.2. Least absolute shrinkage and selection operator

Tibshirani (1996) introduced the least absolute shrinkage and selection operator (LASSO) as a variable-selection approach for regression models. The method minimizes the residual sum of squares subject to the absolute values of the regression coefficients. LASSO performs variable selection and regularization simultaneously to enhance the interpretability and precision of statistical models (Tibshirani 1996). This study applies LASSO to determine important meteorological variables for predicting river discharge.

Assuming a sample contains M events where each event has p number of independent variables and one dependent variable, let y_i be the dependent (output) variable, and x_i = (x₁, x₂, … , x_p)^T be the vector of i^th independent (input) variables, then the objective function of LASSO is:

where λ is a pre-determined parameter that determines the regularization degree and β = (β₁, β₂, … , β_p) is the vector of regression coefficients. Let X be the matrix of independent variables, i.e. X_ij = (x_i)_j, where i = 1, 2, … , M, j = 1, 2, … , p and x_i^T is the i^th row of X. Then, the above formula in a compact form can be written as:

where βp=∑i=1Mβip1/p is the standard L_p norm, 1_M is a column vector of M dimensions with entries 1. In this study, LASSO is employed using the optimal glmnet library in R language and the optimal value of the LASSO parameter using this library is obtained using a 5-fold cross-validation approach.

2.3. Artificial neural network

The artificial neural network (ANN) is a robust modeling tool in which information processing is a representation of biological systems (Kachrimanis et al. 2003). The network is constructed from interconnected neurons, which can determine values from the inputs through network processing. The neuron receives input signals and provides the output signal that mainly depends on the neuron processing function. The ANN architecture consists of a series of interlinked neuron layers. Every layer is linked with another layer through neurons, which transfer information between these layers. Through this processing, the information reaches the output (dependent variable) layer. The ANN mechanism follows four assumptions:

Figure 1a shows the architecture of the ANN model, and Figure 1b presents the structure of a neuron where every input (independent variable) comes from other neurons and are multiplied by their weights (w_j; j = 1,2, … , n) respectively and then aggregated with the bias (b) vector. This aggregated input (s) is passed using the transfer or activation function (f) to obtain the output (a) of a specific neuron. Letting x be the vector of independent (input) variables, the neural network maps into another output vector a through:

a= f(x. w +b) (4)

The mean squared error (MSE) is computed and using the back-propagation process, the weights of the entire network are modified in the training process. The accuracy of the ANN depends on the quality and amount of data in training.

Fig. 1.

The mathematical model of ANN (a) and systematic representation of a neuron (b).

In this study, the ANN algorithm is trained by a back-propagation technique where the output and input variables are applied in the network. The Broyden-Fletcher-Goldfarb-Shanno (BFGS) optimization is employed in a three-hidden-layer network. In the input layer of the ANN algorithm, the activation function is applied with 1000 iterations in the hidden layers. In this study, the ANN algorithm is applied using the validant library in the R programming language.

2.4. Empirical mode decomposition

Huang et al. (1998) introduced empirical mode decomposition (EMD) as an adaptive method for signal analysis. The EMD is designed to analyze non-linear series. The EMD approach assumes that a signal contains different intrinsic mode functions (IMFs) of oscillations. Every mode has the same number of extrema and zero-crossings. There is a single extremum between successive zero-crossings. In this way, the signal is decomposed into different IMFs and residuals. A component is an IMF if it satisfies two conditions: (i) the number of extrema and the number of zero-crossings must be equal to one or differ at most by one, and (ii) at any point, the average of the envelope is zero (Huang et al. 1998). Any original signal y(t) can be decomposed using the EMD algorithm as follows (Lei et al. 2003; Jungsheng et al. 2006):

a) Find the local minima and maximum through the cubic spline line as the upper envelope and lower envelope, respectively.

b) Find the mean (m₁) of upper and lower envelopes.

c) The difference between the y(t) and the 1^st component m₁ is the first component denoted as h₁ i.e. h₁ = y(t) − m₁. If h₁ is an IMF, then it is said to be the first IMF component of y(t).

d) If h₁ is not an IMF, then it is treated as an original signal, and the steps (a)-(c) are repeated, then h₁ − m₁₁ = h₁₁.

After repeating the sifting process k times, h_1k becomes an IMF, i.e. h_1(k−1) − m_1k = h_1k, then it is termed as:

c₁ = h_1k (5)

The first IMF component from the data.

e) Next, subtract c₁ from y(t) to obtain u₁ = y(t) − c₁ where u₁ denotes the treated data, and the process is repeated n times to get n IMFs of y(t). Then,

At the end of the process, we have IMFs (c_j; j = 1,2, … , n) and residual (u_j). By summation of all the components, the original signal y(t) can be obtained as:

The EMD method is implemented using the EMD library in R language in this study.

2.5. Support vector machine

Support vector machine (SVM) is a popular modeling technique for classification and regression problems. The SVM algorithm maps complex high-dimensional data into high-feature space (Vapnik 1995). We assume a training set with n observations, {x_d, y_d}, d = 1,2, … , n, x_d ε R, y_d ε R, where y_d denotes the estimated value of the dependent (output) variable, x_d is the corresponding lagged values of the dependent variable, and n is the sample size. Then, the SVM is developed as:

f(x) =ω^Tφ(x) + b (8)

where f(x) is the estimated dependent variable, b ∈ R is the bias, and ω ∈ R represents the vector of weights. The transfer function φ(x) maps input data into high-dimensional space. The Eq. (8) is solved by risk minimization as follows:

where c > 0 represents the penalty parameter, ξ and ξ^* are slack variables that show the upper and lower constraint of f(x), and ε denotes the insensitive loss function. Further, the Lagrangian function is used as the non-linear regression function, which replaces φ(x) and ω in Eq. (8) as:

where k(x, x_d) = ⟨φ(x), φ(x_d)⟩ is the kernel function. The α_d^* and α_d represents the Lagrange coefficients.

In this study, SVM is applied to capture the features of the error series using the radial basis function (RBF) kernel, i.e. k(x, x_d) = ex−xd2g2, where g is the width of RBF (Baydaroğlu et al. 2018). The SVM algorithm is applied in this study using the R language e1071 library with unit cost and g = 1/m where m is the number of input variables.

3.1. Limitations of LAES hybrid framework

The efficiency of the LAES hybrid framework depends on the optimal choice of parameters of the LASSO approach. This framework works efficiently when the independent variables are selected using the optimal value of the LASSO parameter and the information loss by dropping variables is minimal. A high value of the LASSO parameter can contribute toward a loss of information, which may result in poor model fit. Secondly, the performance of the proposed hybrid method depends on the availability of data variables that may vary in different regions of the world due to differences in weather characteristics. The performance of the LAES hybrid model may vary with respect to changes in region (or location) of study and climatic conditions.

4.1. Description of data

The Khyber Pakhtunkhwa province is a mountainous region, including the Tirich Mir, Lalazar, Hindu Kush, and some other mountain ranges. The changing climate of this region affects air temperature, water flows, precipitation, and groundwater resources for irrigation systems and domestic use. These conditions make the northern area of Pakistan prone to drought or flooding due to changing environment and weather conditions.

The Kabul River begins at the Unai pass base from the Hindu Kush mountains in Afghanistan, flowing toward the east and spanning 700 km to drain into the Indus River of Pakistan (Mehmood et al. 2021). The Kabul River at Nowshera station is located at a latitude of 34°0'25''N and longitude of 71°58'50''E. The hydrometeorological regime is characterized by rain in the spring and snow in the winter. The melting of glaciers in summer is increasing each year due to high temperatures, leading to rising water levels in the river (Rasouli 2022). In addition, rainfall in the monsoon season also affects water levels in the river. The Kabul River is influenced by varying climatic conditions, which may lead to hydrometeorological hazards (i.e., heatwaves, floods or drought).

Figure 3 shows the location of the Kabul River in Pakistan. Kabul River data was collected from the Surface Water Hydrology Project (SWHP) Department of the Water and Power Development Authority of Pakistan (WAPDA) from 1st January 2005 to 31st December 2017. The data contain river discharge and meteorological variables. The meteorological variables include air temperature (minimum and maximum), pan water (minimum and maximum), relative humidity (8 AM and 5 PM), dew point (8 AM and 5 PM), evapotranspiration, and wind speed. Average temperature and precipitation have high variability across the basin. River flow has been high during the monsoon period in Pakistan, particularly in July and August. In the midst of 2005, 2010, and 2015, there was extensive flooding due to high temperatures and heavy rainfall in the region. The discharge had some missing values, which were replaced with the monthly average (mean) value. Outliers present in the data were also replaced by median of the respective month. The number of observations for each variable is 4748, approximately 365 daily values for 13 years.

Fig. 3.

Location of Kabul River in Pakistan.

Table 1 shows summary descriptions of all the variables of the Kabul River data. The air temperature (maximum), air temperature (minimum), pan water (maximum), pan water (minimum), dew point (8 AM and 5 PM), relative humidity (8 AM and 5 PM) have negatively skewed distributions, while river discharge, wind speed, evapotranspiration, precipitation and rainfall have positively skewed distributions. The average discharge in the Kabul River is 871.8 m 3/s. Figure 4 shows the Kabul River discharge series. It shows that there are non-linear relationships between river discharge and all meteorological variables.

Fig. 4.

Kabul River discharge series.

The data variables were normalized using the following (Duan et al. 2021):

where z is the original data variable, z_normal is the normalized data variable, z_min is the minimum value, and z_max is the maximum value of the original data variable. After normalization, the dataset is divided into two parts, where 80% of the data are used for training and the remaining 20% for testing (Kisi et al. 2021; Shabbir et al. 2022). The performance of models is evaluated by 5-fold cross-validation using different performance evaluation measures and the average results of these indicators for training and testing data.

4.2. Performance evaluation measures

The prediction performance of the proposed hybrid framework is evaluated on both training and testing datasets. A 5-fold cross-validation approach and different goodness-of-fit measures are selected to assess the performance of models. These measures include root mean square error (RMSE), mean absolute percentage error (MAPE), root-relative square error (RRSE), mean absolute error (MAE) and coefficient of determination (R ²). These measures are given as follows (Zeinali et al. 2020; Shabbir et al. 2023):

where n denotes the total number of observations, y_j denotes the actual observation and ŷ_j denotes the predicted values. The terms ȳ and y^¯ denote the average of observed and predicted values, respectively.

To compare the performance of the different models for river discharge prediction, the improvement percentages of RMSE, MAPE, RRSE, and MAE are also used and are given as:

where subscript i denotes the competing model and subscript j indicates the proposed LAES hybrid model. These quantities indicate the degree of improvement in the prediction performance of one model relative to another model (Duan et al. 2021).

The Diebold-Mariano (DM) test has been widely used in literature to compare the forecast accuracy of two models (Silva et al. 2021; Shabbir et al. 2022). The null and alternative hypotheses are:

H₀: E[d_t] ≥ 0 (22)

H₁: E[d_t] < 0

where d_t is the difference loss function, i.e., d_t = e_1t − e_2t, e_1t and e_2t denotes the set of prediction errors of two competing models. The test statistic is DM=d¯2πf^d0/m1/2, where m is the length of prediction errors, d¯=1m∑t=1mdt is the average loss differential between two predictions, The DM statistic follows the standard normal distribution (i.e. N(0,1)) and f^d0 is the spectral density. The 2πf^d0 is the consistent estimator of the asymptotic variance. The null hypothesis (H₀) is rejected if DM<−Z_α, where Z is the standardized normal percentile with probability α.

In this study, a one-sided DM test is used to compare the prediction accuracy of the LAES model with six models. This test uses subscript 1 for the proposed LAES model and subscript 2 for the competing models. This test is applied using the squared differences loss function to compare models at a 1% significance level. If DM < −2.326, we will reject the null hypothesis. The proposed LAES hybrid model is compared with MLR, SVM, ANN, LASSO-MLR, LASSO-SVM and LASSO-ANN models in this study.

5.1. Comparison of model accuracy

The daily river discharge was estimated against various meteorological variables. Table 2 presents the training and testing phase results of daily river discharge prediction. In the training phase, the MLR model is the worst performer among all models (RMSE = 533.822 m ³/s, MAE = 378.003 m ³/s, RRSE = 0.711, MAPE = 66.786% and R² = 49.4%). However, the SVM and ANN models performed relatively better than the MLR model. For example, in the training phase, the RMSE for MLR, SVM and ANN models is 533.822 m ³/s, 511.262 m ³/s and 507.015 m ³/s, respectively. Similar to this study, Zhang et al. (2018b) found that the MLR model is the worst performer for predicting river discharge in the East River basin of China. Some other studies found that the non-linear features of river discharge are captured well by SVM and ANN models (see Poul et al. 2019 and Meng et al. 2021).

Comparing the performance of models based on meteorological variables selected by LASSO, we found that the performance of all models is improved in most of the instances. The performance of the LASSO-MLR model is better than the MLR model in the testing phase (RMSE = 543.559 m ³/s, MAE = 381.889 m³/s, RRSE = 0.725, MAPE = 67.758% and R² = 47.4%). However contrary results are obtained in the training phase, in which the LASSO-MLR model has a similar fit to the MLR model. The prediction ability of LASSO-ANN and LASSO-SVM is better than ANN and SVM models respectively. Mehr and Gandomi (2021) found that LASSO improved the predictive ability of a multi-stage genetic programming model by reducing the number of genes for predicting river discharge in the Sedre River of Turkey. In the training phase, the proposed LAES hybrid model has the best fit for river discharge data based on various performance criteria (RMSE = 302.952 m ³/s, MAE = 201.022 m ³/s, RRSE = 0.404, MAPE = 30.494% and R² = 83.7%).

Comparing the results in the testing phase, the MLR model has the poorest performance when all the meteorological variables were used as inputs (RMSE = 554.277 m ³/s, MAE = 383.541 m ³/s, RRSE = 0.739, MAPE = 68.134% and R² = 45.3%). The use of LASSO for dimension reduction enhanced the performance of MLR, SVM, and ANN models in the testing phase. Judging by RMSE, RRSE and R ², the LASSO-ANN model is a better performer than the LASSO-SVM and LASSO-MLR models. However, comparing MAE and MAPE, the LASSO-SVM model performs better than the LASSO-MLR and LASSO-ANN hybrid models (MAE = 307.124 m ³/s and MAPE = 39.394%). The proposed LAES model outperforms all competing models in the testing phase (i.e., RMSE = 337.143 m ³/s, MAE = 218.353 m ³/s, RRSE = 0.449, MAPE = 32.354% and R² = 79.8%). Overall, the proposed LAES hybrid model has higher prediction accuracy than single and LASSO-based ANN, SVM, and MLR models.

Figure 8a presents the goodness-of-fit measure values of all the models considered in the study in both training and testing data. It shows that the proposed LAES hybrid model has the highest accuracy among all models considered in the study. The Taylor diagram in Figure 8b shows that the proposed LAES model is the most efficient among all models considered in predicting daily river discharge based on its relationship with meteorological variables.

Fig. 8.

Prediction results of models in training and testing phase (a) and Taylor diagram (b) for the Kabul River.

The improvements of the proposed LAES hybrid model are shown in Table 3 in terms of P_RMSE, P_MAE, P_RRSE, P_MAPE and P_R² for both training and testing phases. The proposed LAES hybrid model has 43.3%, 40.7% and 40.3% lower RMSE than the MLR, SVM, and ANN models, respectively, in the training phase. The findings indicate that the MLR model is least efficient for non-linear data, consistent with the findings of Zhang et al. (2018b).

Comparing the LAES model to LASSO-based models, we found that their promoting improvements were lower compared to single MLR, SVM, and ANN models in the majority of the scenarios. During testing, the reduction in RMSE by LASSO-MLR and MLR models is 38% and 39.2%, respectively. Similarly, the improvements by the LAES model vs. the SVM model (36.1%) are higher than the LAES model vs. the SVM model (32.6%). The proposed LAES hybrid model has 68.2%, 43.6%, and 42.7% better prediction accuracy than the LASSO-MLR, LASSO-SVM, and LASSO-ANN models. Kang et al. (2023) also stated that LASSO helps enhance the predictive performance of monthly run-off, which is influenced by meteorological events.

Generally, the proposed LAES hybrid model has promising predictions compared to all six models. During the training phase, the MAE of LAES compared to MLR, SVM, ANN, LASSO-MLR, LASSO-SVM, and LASSO-ANN decreased by 46.8%, 35.1%, 39.9%, 46.9%, 28.4%, and 33.6% respectively. These results are in agreement with the findings of Duan et al. (2021). They reported that the decomposition-based error correction approach significantly improves the accuracy of models.

The DM test results on the testing data of Kabul River discharge are given in Table 4. The null hypothesis for all competing models is rejected at a 1% significance level. Thus, the prediction accuracy of the proposed hybrid LAES model is higher than the six benchmark models. Therefore, the DM test confirms that the proposed LAES hybrid model has higher prediction accuracy than the competing models in predicting river discharge.