Forecasting infection fatality rate of COVID-19: measuring the efficiency of several hybrid models

Full Text

Introduction

A highly contagious respiratory illnes, COVID-19 caused by the SARS-CoV-2 virus was initially discovered in China, in December 2019, and since then, it has spread over the world.

This pandemic has had significant impacts on many aspects of life, including public health, the economy, education, and social interactions. For this reason, forecasting COVID-19 is an important tool in managing the pandemic, helping to minimize its impact and inform public health interventions.

Numerous researchers extensively explored this pandemic. For instance Alzahrani et al. [4] used ARIMA model to predict the spread of the pandemic, Dahesh et al. [8] treated the new cases using ARIMA model, Roy et al. [16] focused on spatial prediction. Note-worthy contributions also include the research conducted by Rath et al. [15], Chen [9], Lukman et al. [13], Yousaf et al. [23].

Numerous researchers have examinated this phenomena through the application of genetic algorithm, such as the works of Deif et al. [10], Salgorta et al. [17], Acosta et al. [1].

The deep learning tools are also used to predict the new cases such as the work of Alazab et al. [3], Tamang et al. [19] Kapoor et al. [12] and Namasudra et al. [14] use neural networks, Zeroual et al. [24] make a compartive study between different deep learning models, ArunKumar et al. [5] compared between statistical models ARIMA, seasonal ARIMA model and machine learning models Gated Recurrent unit (GRU), Long-Short term memory (LSTM).

IFR stands for infection fatality rate, which is the proportion of people who die from an infectious disease among all those who have been infected, was also estimated and forecasted in many works such as Singh et al. [18], Vattay et al. [22] Forecast the outcome and estimating the epidemic model parameters from the fatality time series. Ahmar et al. [2] use ARIMA and nonlinear AR model.

In the remainder of this paper we deal with forecasting daily IFR using hybrid models then we evaluate their effectiveness. In the second section we present both the data and the descriptive statistics, which provide insights into the behavior of the phenomenon. Moving on to the third section we elaborate the methodology employed in our study, The final section encompasses the implementation of our approach which includes statistical models (ARIMA, BATS) and hybrid models (ARIMA-SVM, BATS-SVM, BATS-RVM, ARIMA-Random Forest and BATS-Random Forest), These models are subsequently subjected to comparison using performance metrics such as RMSE, MSE, and MAE.

Materials and methods

Forecasting Models

ARIMA and ARFIMA model

Autoregressive integrated moving average (ARIMA) models predict future values based on past values, it gauges the strength of one dependent variable relative to other changing variables.

A stochastic process (X_t)_{t ≥ 0} is said to be an ARIMA(p, d, q) an integrated mixture autoregressive moving average model if it satisfies the following equation:

ϕ(L)(1 − L)^dX_t = θ(L)ε_t ∀t ≥ 0 (1)

where d ∈ N, L is lag operator, ε_t ∼ N (0, σ ²) i.i.d. errors, with σ ² < ∞.

ϕ(L) = (1 − ϕ₁L − ... − ϕ_pL^p) with ϕ_p≠ 0

θ(L) = (1 − θ₁L − ... − θ_qL^q) with θ_q 0

• In the case of d = 0, we obtain ARMA(p, q) process;
• In the case of d ∈ R, coincide with Fractional ARIMA(p, d, q) process.

BATS model

The BATS (Exponential smoothing state space model with Box-Cox transformation, ARMA errors, Trend and Seasonal components) model is a time series forecasting model that was proposed by De Livera et al. [11].

Box-Cox Transformation component is used to transform the data to achieve normality and stabilize the variance. The ARMA (Autoregressive Moving Average) Errors component is used to model the residuals of the time series data, which are assumed to be independent and identically distributed. Finally, the Seasonal component is used to model the seasonal patterns in the data.

Random Forest (RF)

The random forest for regression algorithm is a machine learning algorithm that combines multiple decision trees to predict continuous target variables. The algorithm works as follows: Select a random subset of the training data, with replacement.

Construct a decision tree for the subset of data by recursively partitioning the data into subsets based on the values of the input features. At each node, randomly select a subset of features to consider for splitting.

Repeat the previous steps to create multiple decision trees.

For prediction, pass the input data through all the decision trees and obtain the predicted target variable for each tree.

Aggregate the predictions of all trees to obtain the final prediction. This can be done by taking the average of the predicted values or using weighted averaging.

SVM model

Support vector machine (SVM) analysis is a popular machine learning tool for classification and regression, it is considered a nonparametric technique because it relies on kernel functions.

Given deviation data of training x_i(i = 1, 2, , m) where x_i Rⁿ is the input vector with n-dimension, y_i Rⁿ is the associated desired output value of x_i. Then the SVMs model is formulated as follows:

f (x) = wϕ(x) + b (2)

Where (φ(x)) is called the feature that is non linearly mapping from the input space x. The w and b are coefficients that are estimated by minimizing the regularized risk function shown in formula:

$min\frac{1}{2}{w}^{T}w+C\sum _{}_{i=1}^{m}L\left(y-f\left(x\right),x\right)$ (3)

where C is the regularized constant determining the trade-off between the empirical error and the regularization term. The larger the constant C is, the more the minimum experience risk is emphasized, and the lower the generalization of function f.

Using the Lagrange function and duality theory, and with the kernel function k(x, x) introduced, the function given in 3 can be transformed into a quadratic programming problem as follows:

$min\sum _{}_{i,j=1}^m(\alpha i - \alpha i*)(\alpha j - \alpha j*)k(xi, xj) + ϵ\sum _{}_{i=1}^m\left(\alpha i + \alpha i*\right) \sum _{}_{i=1}^{m}yi(\alpha i - \alpha i*)$(4)

0 ≤ α_i,α_i∗ ≤ C,

where α_i and α_i∗ are Lagrange multipliers. They are obtained by solving this quadratic programming problem, and the input vector x_i corresponding the nonzero α_i and α_i∗ ia is the support vector. Thus, we transform 6 the following equation:

$f (x, \alpha i, \alpha i*) =\sum _{}_{i=1}^m\left(\alpha i - \alpha i*\right)k(xi, xj) + b$ (5)

i = 1, 2, , m

Empirical results and discussion

Infection fatality rate (IFR) is a measure used to assess the proportion of infected individuals with fatal outcomes. Here is the formula used to calculate daily IFR for COVID-19:

$\frac{Numbe of daily new confirmed deaths}{Number of daily new confirmed cases }\times 100$ (6)

Source of data: World health organization.

We have dealt with new cases and new deaths over the world from January, 3rd 2020 to March 16th 2023 Using (6), we calculate IFR index.

As shown in Fig. 1 the IFR is may appear higher in the first three months but it subsequently decreased significantly. This is due to the social awareness and the effective implementation of public health measures, such as mask mandates, social distancing, and lockdowns.

 

Figure 1. Daily IFR of COVID-19 from January, 3rd 2020 to March 16th 2023

 

Two years into the pandemic, the IFR dropped significantly, approaching close to zero, largely due to widespread vaccination among the population, especially those at higher risk of severe illness, the overall death rate can decrease.

In the early stages of the pandemic, testing availability was limited, and many mild or asymptomatic cases went unreported. As testing capacity increases and more people get tested, health authorities can identify a larger proportion of mild cases.

Methodology

The methodology used in our work is to implement a statistical approach to model the data.

X_t = Y_t + U_t (7)

X_t is decomposed into two parts linear component Y_t and nonlinear component U_t, we apply statistical models (ARIMA and BATS) which are more suited to linear pattern Ŷ(t), then the difference ε_t = X_t – Ŷ(t) is the residual.

The residual series ε_t contains the nonlinear parts, thus we use SVM and RVM to fit residuals gives the predicted values Ȗ(t), finally we combine the two predictive results to get:

$\overbrace{X}{\left(t\right)}_t = \overbrace{Y}{\left(t\right)}_t + \overbrace{U}\left(t\right)$ (8)

Results

We decompose the data into training data and test data, we use test data as horizon of forecasting to validate the results. We treat two cases, in the first case we use horizon of 60 days (short term) in other words the last 60 observations are test data. the second case horizon is 120 days (long term).

In Table 1 We provide a summary of how our data is described, we detect the missing values (Na’s). In this case the mean is close to the median it suggests that the data has a relatively symmetrical distribution, which can be a useful insight for understanding the central tendency and overall shape of the dataset.

 

Table 1. Descriptive statistics for Daily IFR

min	1stQu	Median	Mean	3rdQu	max	NA’s
0	0.373 9	1.5395	1.817 1	2.167 6	29.867 7	8

 

From Table 2 we can conclude some results about our data:

 

Table 2. Some characteristics of daily Infection Fatality Rate

Tests	Daily IFR	Comment
KPSS test	0.01	Non stationary
Kolmogorov-Smirnov Test	2.2e-16	Non normal
Terasvirta test	2.2e-16	Non linear
Hurst Exponent	0.99908	Short term memory

 

The primary purpose of the KPSS test is to check for the presence of a unit root in the time series data. A unit root indicates nonstationarity. If the test statistic is greater than the critical values at a chosen significance level, you fail to reject the null hypothesis, suggesting that the data is nonstationary and this is our case 0.01 < 0.05.

Identifying nonstationarity is important because many time series forecasting models assume stationarity, and addressing nonstationarity may involve transformations or differencing to make the data suitable for modeling.

Kolmogorov–Smirnov test is a valuable statistical tool for assessing goodness of fit between a sample distribution and a theoretical distribution (Normal distribution in our case).

he primary utility of the Terasvirta test is to detect nonlinearities in time series data, instead of using linear models, you may choose to employ nonlinear modeling techniques and that is explain our choice to hybridize statistical models and machine learning models which handle nonlinear patterns.

The Hurst exponent serves as a valuable tool for evaluating whether long memory models are suitable candidates and for estimating the memory parameter in an ARFIMA model.

d = H – ½, thus d = 0.49908 which is close to 0.5, This proximity to 0.5 implies that ARFIMA model non-invertible suggesting that short memory process such as ARIMA more suitable to model our data.

From Figure 2 the ACF function did not decrease hyperbolically which confirm that short term memory model is a good candidate, and we remark clearly from the PACF function has a periodic component.

 

Figure 2. ACF and PACF for Daily IFR

 

Handling missing values

We remove missing values (NA’s), it appears when the number of the new cases is 0, thus it is more accurate to replace missing values with zero’s.

ARIMA:

Due to the ARIMA model’s efficiency in analyzing time series, we first apply it. By using the Box-Jenkins approach and the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) to select the the best ARIMA model, which is in our case ARIMA (5,1,3), we obtain the following equation.

For the horizon of forecasting h = 60, 120.

(1+0.2104L−0.1261L²−0.0137L³−0.3397L⁴−0.3157L⁵)X_t = (1+0.7631L−0.0849L²+0.269L³)ε_t (9)

We have used ARFIMA function to fit the model, we get a Fractional AutoRegressive model FAR(1,0.4998,2) the integrated fraction d is close to 0.5 which make the model not invertible.

We have used Maximum likelihood method to estimate the parameters of ARIMA and ARFIMA models.

BATS:

For the horizon h = 60, 120.

We perform BATS model, the output in R language is BATS (1, 1,4, –, –). The first is the ω parameter of the Box-Cox transform, the second is ARMA order of the errors, the third is φ trend damping, the fourth is the seasonal periods, which in your case are none. ω = 1, meaning that indeed, there is no Box-Cox transformation.

ARMA order of errors is ARMA(1,4) with autoregressive coefficients 0.608 and moving average coefficients 0.658, 0.415, 0.369635, 0.52285 for the horizon h = 120.

SVM:

For the horizon h = 60.

SVM models have various hyperparameters that can be tuned using cross validation method, we have used epsilon regression which is tolerance error, we use also radial kernel K(x_j, x_k) exp(x_j x_k²) it allows to capture nonlinear patterns. ∈ = 0.2, gamma = 2, cost = 512.

For the horizon h = 120.

The best hyperparameters ∈ = 0.2, cost = 4, gamma = 1

Random Forest:

For the horizon h = 60.

As Random forest is an ensemble of decision trees, thus the number of trees is 500 and mean of squared residuals is 3.001598.

For the horizon h = 120.

The number of trees is 500 and mean of squared residuals is 3.162188.

Forecasting Hybrid models

We treat the residuals of ARIMA with SVM then with Random Forest model in order to improve the forecasting results, we do the same thing with BATS-SVM and BAATS RVM.

We tune the hyperparameters using cross validation technique.

Illustration

Performance Measures

RMSE, MAE, and MSE are commonly used evaluation metrics in machine learning and statistics to assess the performance of models. They are used to measure the accuracy of predicted continuous values compared to the actual values.

$RMSE=\sqrt{\frac{1}{n}\sum _{}_{i=1}^n{\left(y_i-z_i\right)}^2}$ (10)

$MSE= \frac{1}{n}\sum _{}_{i=1}^n{\left(y_i-z_i\right)}^2$ (11)

$MAE=\frac{1}{n}\sum _{}_{i=1}^n\left(\left|y_i-z_i\right|\right)$ (12)

Discussion

A small MAE means that our model is excellent at predictions, while a large MAE suggests that our model does not perform well at predictions. Unlike MSE, we do not square the residuals, thus MAE is more robust to outliers.

A higher MSE indicates that the model will be penalized for making predictions that significantly differ from the actual value. This means that a large difference between predicted and actual values will be more heavily penalized in MSE than in MAE. RMSE and MSE are sensitive to outliers.

We remark in Fig. 3 and 4 (see cover III) Hybrid models in two cases for h = 60 and h = 120 perform better than ARIMA and BATS because these models are linear and can not handle nonlinear pattern, contrarily SVM can treat nonlinear patterns due to its form of kernel, Random Forest has a good performance because it aggregate the results of many decision tree.

We are unable to make predictions using ARFIMA due to the parameter d being extremely close to 0.5 which makes the model non-invertible The advantage of hybrid models is that we can deal with linear components (ARIMA and BATS) and nonlinear component (SVM and RF).

In Fig. 3 (see cover III), for h = 120, BATS has a good performance then ARIMA because it deals with short term memory phenomenon, same thing for BATS but this model can deal with complex pattern in time series such seasonal component and it assumes that there is a correlation between the errors.

 

Figure 3. Daily IFR Forecasting with horizon h = 60

 

In Fig. 4 (see cover III), for h = 120, ARIMA-RF has a good performance for the 60 days, but thereafter, its performance declines, primarily attributed to the limited memory capacity within Random Forest (RF) models.

 

Figure 4. Daily IFR Forecasting with horizon h = 120

 

Based on the performance measures in Table 3, it is evident that ARIMA-RF consistently demonstrates the superior performance in both scenarios.

 

Table 3. Performance Measures

h = 60	RMSE	MSE	MAE	h = 120	RMSE	MSE	MAE
ARIMA	0.367	0.135	0.248	ARIMA	0.369	0.136	0.272
ARFIMA	NaN	NaN	NaN	ARFIMA	NaN	NaN	NaN
BATS	0.388	0.151	0.265	BATS	0.363	0.132	0.289
SVM	0.362	0.131	0.243	SVM	0.345	0.129	0.266
RF	0.308	0.095	0.234	RF	0.377	0.142	0.278
BATS-SVM	0.345	0.119	0.251	BATS-SVM	0.345	0.129	0.264
BATS-RF	0.301	0.094	0.219	BATS-RF	0.360	0.094	0.289
ARIMA-SVM	0.324	0.105	0.216	ARIMA-SVM	0.342	0.117	0.266
ARIMA-RF	0.290	0.084	0.215	ARIMA-RF	0.268	0.071	0.183

 

Conclusion

To sum up, we have studied the behavior of daily IFR of COVID-19 using statistical models and machine learning models, we have fitted the values of daily IFR which help us to understand the phenomena. To improve the forecasting results we have put into practice hybrid models such as ARIMA-SVM, BATS-SVM, ARIMA-RF and BATS-RF, basing on performance measures ARIMA-RF is the best model.

Acknowledgement

We acknowledge the support of “Direction Générale de la Recherche Scientifique et du Développement Technologique DGRSDT”.MESRS ALGERIA.

About the authors

D. Seba

Higher School of Informatics; University of Bejaia

Author for correspondence.
Email: d.seba@esi-sba.dz

Doctor in Mathematics, Assistant Professor, Laboratory of Applied Mathematics, Department of Mathematics

Алжир, Sidi Bel Abbes; Bejaia

K. Belaide

Higher School of Informatics; University of Bejaia

Email: d.seba@esi-sba.dz

Doctor in Mathematics, Full Professor, Laboratory of Applied Mathematics, Department of Mathematics

Алжир, Sidi Bel Abbes; Bejaia

References

Supplementary files