Forecasting COVID-19 recovered cases with Artificial Neural Networks to enable designing an effective blood supply chain.

This study introduces a forecasting model to help design an effective blood supply chain mechanism for tackling the COVID-19 pandemic. In doing so, first, the number of people recovered from COVID-19 is forecasted using the Artificial Neural Networks (ANNs) to determine potential donors for convalescent (immune) plasma (CIP) treatment of COVID-19. This is performed explicitly to show the applicability of ANNs in forecasting the daily number of patients recovered from COVID-19. Second, the ANNs-based approach is further applied to the data from Italy to confirm its robustness in other geographical contexts. Finally, to evaluate its forecasting accuracy, the proposed Multi-Layer Perceptron (MLP) approach is compared with other traditional models, including Autoregressive Integrated Moving Average (ARIMA), Long Short-term Memory (LSTM), and Nonlinear Autoregressive Network with Exogenous Inputs (NARX). Compared to the ARIMA, LSTM, and NARX, the MLP-based model is found to perform better in forecasting the number of people recovered from COVID-19. Overall, the findings suggest that the proposed model is robust and can be widely applied in other parts of the world in forecasting the patients recovered from COVID-19.

Introduction

After being seen for the first time in December 2019 in Wuhan, China, COVID-19 became a pandemic within a matter of months [1,2]. Reportedly for China, 81% of the patients showed mild symptoms during the recovery period and the mortality rate was 2.3%. Moreover, approximately 5% of the patients suffered more serious symptoms, including respiratory failure, septic shock, and multi-organ failure; in this group, the mortality rate reached 5% [3]. It is widely found that the spread of COVID-19 varies among men and women, different age groups, and those with underlying medical conditions [4]. In response to the pandemic, governments all over the world including Turkey have quickly taken a series of drastic measures to slow down its spread.

A challenging task in developing effective infection prevention and control strategies to stop the spread of the virus is being able to forecast new cases, deaths, and recoveries. Many studies have already proven the effectiveness of various forecasting approaches in estimating daily new infections and deaths. However, the efforts dedicated to forecasting the recovery cases are relatively limited. Some patients can be treated using the plasma received from recovered patients who are reportedly able to donate plasma after 14 days [[5], [6], [7]]. To facilitate the process of promptly transferring the plasma from recovered patients to new COVID-19 patients, developing an effective supply chain system is crucial. Therefore, accurately forecasting the number of recovered people as potential donors is key in developing such a system. As such, this study introduces a forecasting model to estimate the number of recovery cases for effective plasma donation in CIP therapy that is a treatment of COVID-19. The proposed model is based on the Artificial Neural Networks (ANN) and hypothesized to achieve a high prediction rate.

Research to develop vaccines and drugs for COVID-19 is ongoing and no effective therapy or antiviral drug has been found yet. At this point, a passive antibody therapy CIP that is currently applied to COVID-19 patients is an option for COVID-19 therapy. In producing CIP, the liquid component (plasma) of whole blood is taken from recovered COVID-19 patients [8]. Better planning and execution of supply chain processes with plasma donation to patients of CIP therapy is crucial for effective management of the fight against COVID-19. It is aimed to develop a forecasting model in this study that will play a key role in designing an effective blood supply chain mechanism in the struggle with the treatment of this disease. Therefore, we focus on the problem of forecasting the number of people diagnosed with COVID-19 who have completed successful treatment (who have a negative COVID-19 test result) and now have the potential to donate plasma.

Blood, plasma, or stem cell donations were studied in many dimensions before the pandemic. This issue is addressed from many different aspects, such as the need for immediate use, supply and storage of blood and other tissues, as well as the psychological aspects of people's reluctance to donate blood, misdirection, and anxiety about getting sick. The severe progress of COVID-19 in Turkey, as in all other countries, accelerates the studies focusing on the elimination of this disease as soon as possible; and at this point, CIP therapy stands out as one of the most effective treatments to tackle the disease [[8], [9], [10], [11]]. Recovered patients can donate plasma through the Red Crescent and pandemic hospitals in Turkey. At this point, the number of recovered patients must be already determined. The donation process of blood and other tissues is a vital problem even in pre-pandemic periods. In this context, it clearly shows that proper planning of the plasma donation process under pandemic conditions will contribute significantly to the literature.

Related works

The impacts of the COVID-19 pandemic on the status quo have been attracting researchers' attention. The researchers focus specifically on forecasting cases of COVID-19 for death and infections. In this section, a systematic literature review is conducted, which adopts the Preferred Reporting Items for Systematic Reviews and Meta-Analyses (PRISMA) to gather literature by reviewing the studies focusing on forecasting of COVID-19 cases. PRISMA is a systematic approach developed by David Moher to review literature [12]. The PRISMA approach has five steps in literature search: defining criteria, identifying sources, selecting literature, collecting data, and selecting data items [13].

We adopt the PRISMA approach as Systematic Literature Review Method (SLR) to minimize bias in the review process and make it more systematic [14]. The literature search was performed from January 9, 2020 to August 1, 2021, with the keywords used shown in Table 1 .

Table 1

Studies found in literature search.

Database	Details of the Search	Number of studies
SCOPUS	(TITLE-ABS-KEY ( covid-19) AND TITLE-ABS-KEY (prediction) AND TITLE-ABS-KEY ("artificial neural network") )	145
	(TITLE-ABS-KEY ( covid-19) AND TITLE-ABS-KEY (estimation) AND TITLE-ABS-KEY ("artificial neural network") )	21
	(TITLE-ABS-KEY ( covid-19) AND TITLE-ABS-KEY (forecasting) AND TITLE-ABS-KEY ("artificial neural network") )	82

A total of 248 papers is found as a result of the search in the Scopus database with the keywords shown in Table 1. Later, papers directly related to this study are determined by specifying both inclusion and exclusion criteria for the studies found. Thus, it is planned to reveal similar studies that can form a background for the approach adopted in the study and compare the inferences. Table 2 shows the inclusion and exclusion criteria for the papers.

Table 2

Inclusion and exclusion criteria for the papers in the literature search.

Inclusion Criteria	Exclusion Criteria
The studies include forecasting implementation for COVID-19 cases	Studies whose full text could not be reached
The studies include forecasting implementation for COVID-19 deaths	Studies that do not explicitly mention the method used and the results
The studies include forecasting implementation for COVID-19 recoveries	Studies are adopting different artificial intelligence approaches in the prediction of COVID-19 cases
Studies using up-to-date forecast data on COVID-19 and taking advantage of the ANN method	Studies written in languages other than English

When the studies in the literature search are filtered with the criteria in Table 2, it is seen that 26 studies are similar to the methodology adopted in this paper. To reveal the differences between these studies with our paper and to determine their contributions, a summary table has been presented in Table 3 .

Table 3

The literature that adopted the ANN approach for COVID-19 forecasts.

#	Authors	Aim	Input variables	Contributions	Limits
1	Ahmad and Asad [15]	Predicting of coronavirus COVID-19 cases	Actual confirmed casesActual DeathsActual Recoveries	A new ANN design is proposed to estimate the number of deaths, recovered and confirmed COVID-19 cases	Insufficient number of input parameters collected from real life
2	Bodapati et al. [16]	Forecasting of Daily Cases, Deaths Caused, and Recovered Cases	Number of casesNumber of DeathsThe total number of people recovered from the virus	Ability to integrate with an application that streams live data from government sites to view real-time graphs of data	Exposure of the model to limited data
3	Saba and Elsheikh [17]	Predicting the different number of COVID-19 cases at the end of the epidemic	Confirmed cases	Applying ARIMA and NARANN methods to forecast COVID-19 cases	Not making comparisons using different time series methods
4	Istaiteh et al. [18]	Predicting the COVID-19 cases in each country all over the world	Confirmed cases	Using spatio-temporal forecasting for 189 countries worldwide to predict COVID-19 cases	Encountering high error rates with the methods adopted for some countries
5	Al-qaness et al. [19]	Forecasting the number of confirmed cases of COVID-19 in Brazil and Russia	Confirmed cases	Using an enhanced version of marine predators algorithm (MPA), called chaotic MPA to improve ANFIS performance and avoid the shortcomings of traditional ANFIS	Having longer computation time of the proposed method than the compared methods
6	Abbasimehr and Paki [20]	Predicting COVID-19 confirmed cases	Confirmed cases	Utilizing Bayesian Optimization in determining estimation parameters and adopting a multi-output modeling approach	Exposure of the model to limited data
7	Elsheikh et al. [21]	Forecasting the number of total confirmed cases, total recovered cases, and total deaths in Saudi Arabia	Total number of confirmed casesTotal recoveriesTotal deaths	Using seven different statistical evaluation criteria in the estimation accuracy of the model (RMSE, R-squared, MAE, efficiency coefficiency, overall index, coefficient of variation, and coefficiency of residual mass)	Using a limited number of data
8	Hamadneh et al. [22]	Forecasting the number of confirmed and recovered cases of COVID-19	The requested date	Determining the parameters of the MLPNN using the prey-predator algorithm (PPA)	There is no comparative analysis to compare the results
9	Moftakhar et al. [23]	Predicting the number of patients	The observed number of newly infected cases	Applying ANN and ARIMA models for prediction	Models cannot be trained well due to few observations and the Inability to evaluate any risk factor for this disease due to insufficient data on demographic information and social networks of patients.
10	Ünlü and Namlı [24]	Predicting COVID-19 confirmed cases and deaths in seven countries	Confirmed casesDeaths	Applying e Support Vector Machines (SVM), Holt-Winters, Facebook's Prophet, and Long-Short Term Memory (LSTM) for forecasting	Using only RMSE in the interpretation of results
11	Niazkar and Niazkar [25]	Estimating the confirmed cases of COVID-19 in China, Japan, Singapore, Iran, Italy, South Africa, and the United States of America.	Chronological data of confirmed and death cases	Applying fourteen ANN-based models to predict the COVID-19 outbreak	Using a limited number of data
12	Hamadneh et al. [26]	Forecasting the number of total confirmed cases, total recovered cases, and total deaths in Brazil and Mexico	The requested data	Using ANN to estimate the number of cases of COVID-19 with prey predator algorithm (PPA)	There is no comparative analysis to compare the results
13	Toga et al. [27]	Predicting the infected cases, the number of deaths, and the recovered cases with ARIMA and ANN in Turkey	Susceptible casesDaysCurfewsLaboratory tests	Determining the susceptible case number for each day after the first recovered case
14	Kumari and Toshniwal [28]	Forecasting the COVID-19 outbreak in India with ANN	Cumulative confirmed, new, and cumulative deceased cases recorded daily	Utilizing a mathematical curve fitting model to understand the performance of the proposed model	Not conducting a comparative analysis
15	Conde-Gutiérrez et al. [29]	Estimating the cumulative number of deaths from COVID-19 in México	Cumulative number of deaths	Comparing ANN with Gompertz model in estimating the number of deaths
16	Safi and Sanusi [30]	Predicting the total confirmed cases and the recovery or death rate worldwide	Total casesDeath rates	Applying the ARIMA, ETS, ANN models, and the hybrid combination of the three models.	Using data around the world and not separating regions or countries by region (??)
17	Wieczorek et al. [31]	Forecasting number of cases each day worldwide	Information from the last 12 days plus geolocation coordinates of latitude and longitude have an impact on the cases correlations between neighboring countries	Proposing a model, which can work as a part of an online system as a real-time predictor to help in the estimation of COVID-19 spread worldwide	Using a limited number of data
18	de Barros Braga et al. [32]	Estimating the daily and cumulative cases and deaths caused by COVID-19 and demand for hospital beds	Cumulative casesCumulative deathsMunicipal demography Occurrence dateName of State's municipalities Names of the health regions	Training ANN with data from 6 different moments for providing the ability to evaluate the forecasting quality	Using a limited number of data
19	Shetty and Pai [33]	Forecasting the number of infected cases	Daily reported cases	Applying a fast training algorithm that is Extreme Learning machine to reduce the training time and using cuckoo search (CS) algorithm to select the parameters	Not conducting a comparative analysis
20	Tamang et al. [34]	Estimating the number of rising cases and deaths in India, the USA, France, and the UK	The number of days	Presenting intelligent based optimum curve fitting and forecasting for different non-linear models	Not conducting a comparative analysis
21	Melin et al. [35]	Predicting confirmed and COVID-19 deaths for 26 countries	The confirmed and deaths	Adopting the firefly algorithm for ensemble neural network optimization to COVID-19 time series prediction with type-2 fuzzy logic in a weighted average integration method	Not conducting a comparison with different fuzzy extensions to measure the performance of the adopted model
22	Ardabili et al. [36]	Estimating the cumulative number of cases for five countries	Cumulative number of cases	Demonstrating a comparative analysis of machine learning and soft computing models in forecasting the COVID-19 outbreak	Not applying comparative studies on various machine learning models for individual countries.
23	Ahmar and Boj [37]	Predicting infection fatality rate of COVID-19 in Brazil using NNAR and ARIMA	Total data of confirmed casesTotal data of death of COVID-19	Presenting Neural Network Time Series (NNAR) and ARIMA to Forecast Infection Fatality Rate (IFR)	Using a limited number of dataUsing only MSE as a performance measureNot specifying which approach is used in parameter selection
24	Ardabili et al. [38]	Estimating the cumulative number of cases for five countries total and daily cases worldwide	Total casesDaily cases	Proposing artificial neural network-integrated grey wolf optimizer for COVID-19 outbreak estimations	Not conducting a comparative analysis
25	Alsmadi et al. [39]	Forecasting the cumulative cases of COVID-19 for the four Canadian provinces	The cumulative number of infected cases	Applying three models, smooth transition autoregressive (STAR) models, neural network (NN) models, and susceptible-infected-removed (SIR) models for forecasting the cases	Data reliability problem
26	de Oliveira et al. [40]	Estimating the number of confirmed cases and deaths of COVID-19 for Brazil, Portugal, and the United States	The daily cumulative number of cases and deaths	Comparing different training functions	Not conducting a comparative analysis

In addition, the input data used in these studies and the limitations of these studies, as well as their contributions are summarized and a detailed perspective that can be used in a glance at the relevant literature is presented.

When we examine the studies in Table 3, it has been observed that no analysis has been conducted to validate the results in many forecasting studies. Generally, many of them use much less data than the number of data in our study. In addition, no study has carried out as much detailed input research as in this study, and the problem has not been examined in detail from a multi-faceted perspective. In most of the studies, forecasting has been applied using only time-series data. However, in this study, many input candidates that can effectively forecast the number of recovered patients are identified and the inputs to be used in the proposed model are selected through analytical methods. The dataset in different countries is used to determine the effectiveness of the proposed method and it is concluded that the results of the method are reliable and useful for other countries. Besides, it is investigated whether the proposed model gives the best results by applying to ARIMA, LSTM, and NARX artificial intelligence tools for comparative analysis. As a result of the analysis carried out, it has been seen that the model proposed by adopting multi-layer perceptron (MLP) gives the best results according to the select key performance indicators (KPIs).

Apart from the reviewed studies, our model using MLP forecasts the number of recovered patients to ensure the efficient plasma donation used for treating COVID-19 patients. Unlike existing studies using MLP, a large number of potential input variables are investigated. MLP is a valuable network structure to model complex nonlinear problems [41]. MLPs can handle large datasets for inputs. They also work well even with small datasets. MLP provides faster decision-making after network training [42]. Unlike other statistical models, MLP does not need probabilistic information about the dataset or assumptions regarding probability density functions. MLP also can approximate unknown functions. MLP is implemented in many studies to solve complex real-world problems due to its high self-adapting and self-learning ability [42]. In this paper, the forecasting process using the MLP model has been conducted in Turkey for COVID-19 patients surviving the disease. By estimating patients who recover daily, the number of people who can become donors in the treatment of plasma can be determined for people who have the disease. At this point, when the patients recovering are determined, important information such as the number of patients to donate for CIP treatment, which is currently used as an effective method in the treatment of the disease, will be obtained, and thus effective planning can be realized. This study is thought to contribute to the literature with the subjects can be listed below:

The daily number of recovered patients is forecasted, and a new perspective that will contribute to the treatment process is presented for COVID-19.

For the first time, apart from the number of COVID-19 patients or deaths, the number of recoveries, which can be used in planning the treatment process, is estimated.

Forecasting of the number of recoveries from COVID-19 is modeled using an MLP structure with multiple inputs rather than a time series problem.

The most influential factors to forecast the number of recoveries are determined using correlation analysis and causal relationships.

The rest of the paper is organized as follows: The proposed methodology is presented in Section 3. Section 4 presents the data cleaning process and real case application of the proposed methodology for Turkey. The comparative analysis is presented in Section 5. Finally, the conclusions and future recommendations are given in the last section.

Proposed methodology

This study proposes a novel strategy to develop more accurate forecast models. The proposed strategy can be divided into three steps: (i) candidate input variables are determined, (ii) predictors are selected by correlation analysis and causal analysis, (iii) different hyperparameters are tested to find the best model.

First, the relevant literature is reviewed, and various variables are defined as candidate input variables to determine the daily number of recovered COVID-19 cases. Then, the Pearson correlation coefficients for the candidate variables are calculated, and some candidate variables are eliminated. Thus, the predictors of the number of people who have recovered from COVID-19 are determined. After the predictors are determined, the values of the predictors are compiled, and a dataset is created. Finally, the most appropriate artificial neural network structure is determined using these data, and the daily number of patients recovering from COVID-19 is forecasted. The steps of the proposed method are given in Fig. 1 .

Fig. 1

The proposed methodology.

Determining the predictors

Some predictors have little or no effect on forecast performance or may cause overfitting. Therefore, a predictor selection/reduction process is required for effective forecasting. Correlation analysis is one of the data pre-processing techniques that provides model simplicity and prevents overfitting [43].

Correlation analysis is a statistical method used to determine whether there is a linear relationship between two numerical measurements; the direction and severity of this relationship, if any [43]. In short, it shows whether the changes in variables affect each other (x and y). The correlation coefficient (r) is a value used to establish the relationship between variables, and this coefficient (Pearson's correlation) takes a value between −1 and 1 [44]. Pearson correlation coefficient can determine the level of the linear relationship between the variables [45]. The Pearson correlation coefficient (r) for two independent variables, x and y, is calculated in Equation (1).where represents the number of observations. If the coefficient is close to ±1, it can be said that there is a perfect correlation, which means as one variable increases, the other variable tends to increase (if positive) or decrease (if negative). If the coefficient value lies between ±0.50 and ± 1, it is said to be a strong correlation. When the value lies below + 0.29, then it is said to be a weak correlation. When the value is zero, it can be said that there is no linear correlation between the variables [46].

The predictor selection step can be applied to improve the accuracy of machine learning models such as support vector machine (SVM), NARX. It provides a more relevant, consistent, and robust input dataset to forecast target variable(s).

Artificial Neural Networks (ANNs)

ANNs are computer systems designed to automatically perform new capabilities such as producing and discovering new information through learning [47]. The ANNs can be trained with examples as human intelligence [48]. The more examples are shown in the network training, the higher the success of solving the problem [49]. The training of the network occurs with the acceptance of certain error values. Learning in a network can be considered as a change in the weight matrix [50]. ANNs are used to estimate output values considering input values for estimation. In this context, ANNs can be used in tasks such as forecasting, classification, data association, data interpretation, data filtering [51].

ANNs can simulate the performance of especially complex systems when available experimental data is limited [52]. ANN are nonlinear information and data processing systems. They are developed from processing units called neurons and the connections between neurons [53]. ANN can be trained by changing the weights of these connections to deal with nonlinear problems that cannot be handled with traditional analysis methods [54].

The three main components of the ANNs structure are the neurons, connections established between these neurons, and functions [55]. ANNs consist of layers; the input layer is the layer where data is entered into the network, and the hidden layer(s) is the layer that processes the data received by the network in the input layer. The number of these layers can be one or more. The output layer is the layer where the result is expressed as output after the data is processed in hidden layers.

Multi-layer perceptron

ANN can be either feed-back (recurrent) networks or feed-forward networks. Recurrent neural networks contain a self-connecting hidden layer, and the connections between nodes form a directed graph along a temporal sequence [56]. Recurrent neural networks contain feed-forward loops and feedback loops, which can store information. Increasing the number of weights increases the complexity of the structure for recurrent neural networks [57]. Perceptrons are arranged in layers in the feed-forward networks so that the first layer receives inputs, and the last layer produces outputs. In contrast, the middle layers have no connection with the outside, therefore, they are called the hidden layers [58]. The feed-forward network adopted in this study is MLP. MLP can solve complex nonlinear problems effectively and is generally applied for classification and forecasting [59]. Backpropagation is a supervised learning technique that is used in training MLP [60]. Backpropagation is a version of the Widrow–Hoff learning rule with multiple layers and differentiable nonlinear activation functions [61]. The backpropagation algorithm consists of the feed-forward phase, where the output of the network is determined, and the feed-back phase, in which the weights are updated by spreading the error back to reduce the gradient [62]. In this algorithm, the difference between the output and the target, the error, is gradually reduced to the minimum level by propagating it on all weights [63]. In the backpropagation algorithm, training starts with a random set of weights [64]. In the feed-forward phase, the inputs of the training set are presented to the input layer of the network. The input layer contains the neurons that receive these inputs. The neurons in the input layer transmit the input values directly to the hidden layer. Each neuron in the hidden layer calculates the total value by adding the threshold value to the weighted input values, and processes them with an activation function and transmits them to the next layer or directly to the output layer. After the net input of each neuron in the output layer is calculated by adding the threshold value to the weighted input values, this value is again processed with the activation function to determine the output values. The error value is calculated by comparing the output values of the network with the expected output values [54]. MLP has an input layer, one or more hidden layers, and an output layer composed of neurons.

In this study, an MLP based methodology, which is one of the most known ANN models, is proposed. MLP is developed to express the nonlinear relationship between independent and dependent variables [65]. A basic MLP architecture contains layers of consecutive neurons. MLP is a feed-forward ANN model. Input data is transferred only from the input layer to the output layer. In MLP models, the backpropagation training method is generally used to minimize the error between the output of the network and the actual results of the training set [66]. The main function of backpropagation is to determine the most appropriate weight values using the gradient descent method to identify the most accurate weights used in ANN [67]. The feed-forward step generates a set of weights, then the weights are updated by the backpropagation of the error in each iteration when training an MLP [68].

Each neuron has incoming signals () into neurons through weighted () input connections. ANNs are trained with these weights. The signals, modified according to the weight changes, are summed in the neuron. The results proceed through the activation or transfer function used to modify the results before the signal becomes the output signal. Generally, there are multiple inputs () that enter the neuron. Input weights () are used to calculate the weighted sum () of inputs as given below [69]:

The value becomes the input of the nonlinear activation function, and it is modified according to function. The bias is another input to the activation function. In this way, the bias and weighted sum form the inputs of the activation function [54]. Then, an activation function is applied to produce the output signal (), as follows:where denotes the activation function. Different activation functions can be used in this stage. Sigmoid function is one of the most used activation functions in feed-forward networks [70]. Sigmoid function is represented with [71]:

The linear activation function is described as follows [72]:

Hyperbolic tangent activation function is given in Equation (6) [73]:

To calculate output ) of the MLP structure, the weight value between the hidden neuron and the output neuron () is multiplied by the value of the hidden neuron (). This value is used as in the determined activation function. In this study, the sigmoid function is used as the transfer function.

The supervised learning method is one of the most used techniques to train networks in ANN. In supervised training, both input and output are provided to the network. Then the network processes the inputs and compares the outputs with the targets (desired outputs) [74]. A comparison is performed between the calculated output of the network and the target to determine the network error. The errors are propagated back to the network, and the weights controlling the neural network are updated [75]. This happens over and over again as the weights are constantly being changed. These weights are determined carefully to reach a minimum error or tolerance, as explained earlier. Besides, errors can be used to change network parameters to improve network performance. In other words, the inputs are assumed to be at the beginning and the outputs to be at the end of the causal sequence. Supervised learning is the most common form of neural network training [76].

The mean squared error (MSE) function is used as a loss function for the training and testing of the proposed ANN, whereas to improve the performance of ANN, the correlation coefficient (R) is used. In this study, the early stopping method is used as a regulation method to avoid overfitting [[77], [78], [79]]. Early stopping consists of breaking or stopping the network training when the loss function increases. Two basic rules are considered stopping conditions (i.e., loss function value and loss function change) [80]. The loss function, MSE, is evaluated for each iteration in training. When the value of the loss function begins to increase without decreasing for the next iterations, training stops in the next iteration with the lowest loss function value in the dataset [81]. Accordingly, if the generalization capacity of the model decreases, the training stops. The general flow of the training state of ANN can be shown in Fig. 2 .

Fig. 2

The training of ANN.

While designing the ANN model, the number of hidden layers, the number of neurons in the layers, and the activation functions used can be different. So, determining the number of hidden layers, fixing the number of neurons in each hidden layer, applying a suitable training algorithm and activation functions play a vital role in achieving better network performance. All these values create the neural network. The number of neurons of each hidden layer affects network performance. The learning rate and momentum coefficient are also used to achieve minimum error, and similarly, they can be different to minimize error.

In this study, trials with different network structures are created, and the best results are given in the following sections.

Forecast key performance indicators (KPIs)

Different error models are applied to evaluate the performances of the forecasting models. Through calculated metrics, the accuracy of the forecasting models is measured. The smaller the metric value means the higher the forecast accuracy. In the literature, MSE, root means square error (RMSE), mean absolute error (MAE) and mean absolute percentage error (MAPE) are used as forecasting performance evaluation metrics (KPIs) [[82], [83], [84], [85]]. Different KPIs are used in this study to evaluate the accuracy of the proposed methodology from various perspectives [86]. MSE can be adopted to determine the performance of our proposed model. The MSE can be calculated as in the following [87]:where shows the predicted value and y shows the true value where j = 1,2, …,m.

RMSE is another KPI that can evaluate forecasting models. The mathematical expression for RMSE is given [88]:

MAPE is one of the most common KPIs used to determine forecasting accuracy [89] due to its features such as scale independence and interpretability [90,91]. The MAPE is calculated as [92]:

MAE determines the variance between two continuous variables [93]. Equation (10) is used to calculate MAE [94]:

The KPIs define the difference between the output and target data [95,96]. Many forecasting networks proposed are based on the minimization of the MSE [[96], [97], [98], [99], [100], [101]]. Using the MSE as KPI, the performance of neural networks is measured using a qualitative approach based on the adjustment ratio of weights and MSE [95]. MAPE is also applied in many forecasting problems to examine the accuracy of the models [90,[102], [103], [104]]. For these reasons, we used MSE, RMSE, MAPE, and MAE as the KPIs for the forecast models in this paper.

A real case application for Turkey

In this study, the daily number of patients recovering from COVID-19 in Turkey is forecasted using ANNs. In this context, the variables predicted to affect the number of patients recovered have been determined by interviews with experts from the health sector and data published by the Ministry of Health in Turkey. The input variables are given in Table 4 .

Table 4

Input variables.

	Candidate Input Variables
I-1	Total tests
I-2	Total cases
I-3	Total death
I-4	Intensive cares
I-5	Intubated patients
I-6	Total recovered
I-7	Daily tests
I-8	Daily cases
I-9	Daily deaths
I-10	Daily change in the intensive cares
I-11	Daily change in the intubated patients

Predictor selection

After the experts determine the variables, the data of the variables are collected from the official sources of the Ministry of Health [105]. On March 10, 2020, the Ministry of Health of Turkey announced that a Turkish man who caught the virus while traveling to Europe was the first coronavirus case in the country. Then, on March 13, the Health Minister of Turkey explained that the coronavirus was detected in the relatives of the first patient, who was taken under observation. Making a statement after the Coronavirus Scientific Committee convened on March 27, 2020, the Minister of Health announced that 42 people recovered, 341 people were in intensive care, 2069 positive cases were detected in the last 24 h, and 17 people died due to COVID-19. The number of recovered patients was first announced on this date. Therefore, relevant data for the next 245 days (between March 27 and November 26) are collected and analyzed [105].

Firstly, the correlation coefficients of candidate input variables with the number of daily recovered patients are calculated using Equation (1) and given in Table 5 .

Table 5

Correlation coefficients with the output variable.

	Input Variables	Correlation Coefficient
I-1	Total tests	−0.021
I-2	Total cases	0.178
I-3	Total death	0.129
I-4	Intensive cares	0.602
I-5	Intubated patients	0.402
I-6	Total recovered	−0.159
I-7	Daily tests	0.414
I-8	Daily cases	0.255
I-9	Daily deaths	0.561
I-10	Daily change in the intensive cares	−0.513
I-11	Daily change in the intubated patients	−0.516

After determining correlation coefficients between all variables with the number of daily recovered patients, we try to decide which variables would be used in the forecasting process. Using highly correlated and the most influential input variables can give better forecast results [106]. If the correlation coefficient takes a value between −0.5 and 0.5, that means there is a low or negligible correlation between variables [43]. Therefore, variables that have a negligible effect (correlation coefficient between −0.5 and 0.5) on the number of forecasted daily recovered patients are excluded. So “I-4 Intensive cares”, “I-9 Daily deaths”, “I-10 Daily change in the intensive cares”, and “I-11 Daily change in the intubated patients” are determined as potential input variables in forecasting calculations. Then, the correlation coefficients among the remaining four input variables are calculated as given in Table 6 .

Table 6

Correlation coefficients among the determined input variables.

Variable	I-4	I-9	I-10	I-11
I-4	1	0.8586	0.1038	0.4403
I-9	0.8586	1	0.1402	0.4684
I-10	0.1038	0.1402	1	0.2333
I-11	0.4403	0.4684	0.2333	1

When we examine Table 6, mostly weak relationships appear between the variables. The absolute value of the correlation coefficient between two variables being greater than 0.9 indicates that the values of these two variables increase or decrease together. Hence, it can be reasonable to use one of these as the predictor. The correlation coefficients between individual input variables shown in Table 6 do not take values close to ±1 to eliminate one pair of the highly correlated variables [45]. In this analysis, we do not need to remove any input variables. After all data cleaning processes, we use four different input variables as predictors such as “I-4 Total intensive cares”, “I-9 Daily deaths”, “I-10 Daily change in the intensive cares”, and “I-11 Daily change in the intubated patients”.

After correlation calculations, causal analysis is also performed to determine whether input variables had a causal relationship with the target variable and whether the model fits well. For this aim, we apply structural equation modeling (SEM) to analyze the relations. In this SEM, it is checked whether the inputs have a causal relationship with the target variable. For this aim, we construct the following hypothesis:

I-4 variable has a causal relationship with the target variable

I-9 variable has a causal relationship with the target variable

I-10 variable has a causal relationship with the target variable

I-11 variable has a causal relationship with the target variable

Relationship analysis is conducted with the IBM AMOS program. The diagram obtained from the IBM AMOS program and results are given in Fig. 3 and Table 7 as follows:

Fig. 3

Relationship diagram between variables.

Table 7

SEM results.

			Estimate (β1)	Estimate (β2)	S.E.	C.R.	P
Target Variable	<---	I-11	0.038	1
Target Variable	<---	I-10	−0.253	−5.665	0.927	−6.113	<0.001
Target Variable	<---	I-9	0.948	29.746	2.156	13.795	<0.001
Target Variable	<---	I-4	−0.326	−0.721	0.152	−4.748	<0.001

Fig. 3 shows the inputs and outputs, arrows extending from inputs to outputs to test the hypothesis of whether inputs have an effect on output, error variables denoted by "e", which indicates the error of each factor. In addition, the relations of the errors with each other are also shown by arrows as a result of the analysis.

In Table 7, β1 shows the standardized path coefficients, β2 shows the path coefficients, S.E. shows the standard estimates, C.R. shows the critical ratio, and p shows the marginal significance level in a hypothesis test representing the probability of a particular event occurring.

In the structural model, 4 hypotheses identified as H1, H2, H3, and H4 are examined. The goodness of fit of the SEM model is also evaluated by the chi-square test and the goodness of fit indices, such as Root Mean Square Error of Approximation (RMSEA), goodness-of-fit index (GFI), and comparative fit index (CFI). According to the results obtained in the constructed structural model, it can be interpreted that the model is compatible since all fit indices are within the desired limits (RMSEA<0.08, CFI>0.9, GFI>0.9 [107]). It has also been observed that all inputs are related to the output variable. By checking the p-values in Table 5, it can be interpreted that all the inputs have a significant effect on the target variable.

Forecasting of the daily number of patients recovering from COVID-19 in Turkey

After determining the predictors to be used and the values of these predictors, these values are used in ANN as inputs for forecasting the number of recovered patients in Turkey for COVID-19 cases. In our network structure, the data of the “I-4 Total intensive cares”, “I-9 Daily deaths”, “I-10 Daily change in the intensive cares” and “-11 Daily change in the intubated patients” are adopted as input variables, and the actual number of “the Daily Number of Patients Recovering” is evaluated as the output variable. Fig. 4 shows the values taken by our input variables for the time interval specified for the forecasting process.

Fig. 4

The values of input variables.

As shown in Fig. 4, the number of patients receiving intensive care varies between 500 and 2500 and shows an increasing trend. The number of daily deaths shows a steady trend between May and August, then increases after August. The daily change in intensive care and daily change in intubated patients do not have a trend and, the numbers take both positive and negative values. There are lots of spikes for these two input variables.

The input and output variables are normalized according to Equation (11) before being used in the neural network. The normalized values of all variables take the value in the range [0,1].where and represent values before and after normalization of the variable, respectively. and represent the minimum value and maximum value of the variable, respectively.

In our proposed methodology, n-day lag is used to make the n-day forecast. Therefore, the network is constructed using (245-n) samples. n is assumed to be 7, 10, and 14 to include medically justified recovery times [[108], [109], [110], [111], [112]]. The data of each input variable with an n day delay are used to forecast daily recovered patients. Once the data is compiled, the partition for training and test is applied. Furthermore, the first n-days period is excluded from data due to n-day lags. Subsequently, the remaining data is divided into two parts as training and test data. To achieve the best performance, one and two layers are tried as the number of hidden layers. The number of hidden neurons is also tested between 4 and 12 as the most commonly used values. The learning rate and momentum coefficient are tested starting from 0.1 increased by 0.05 up to 1. In this network structure, TRAINLM (Levenberg-Marquardt backpropagation) is used as a training function. The tangent sigmoid function and logarithmic sigmoid function are tried as transfer functions at the hidden layer(s) and pure linear function in the output layer. All combinations are made for 1000 epochs and run 100 times to test different starting weights. Finally, the most appropriate ANN network structure consists of two layers with eight and five hidden neurons in the first and second layers. The sigmoid function is used as a transfer function for each layer, and the pure linear function is used to determine the value of the output layer. The learning rate is determined as 0.7, and the momentum coefficient is determined as 0.45. The network structure used in this study is shown in Fig. 5 .

Fig. 5

The neural network structure used in this study.

MLP models in this paper have been developed in MATLAB R2020b software environment with different parameter values. The k-fold cross-validation is used to avoid overfitting with the dataset. For this purpose, we randomly divide the training dataset into k-folds. One of these folds is used to test the network, and the remaining folds are used to train the network. This process is repeated k times for each fold. The MSE metric for each fold is calculated, and the cross-validation error is computed by averaging the MSE values. The proposed model with predetermined parameters is validated in this way, and 5-fold cross-validation is also used. To improve the accuracy of the proposed network structure, we performed the re-forecasting process by changing the rate of training data. We compare the forecast performance with respect to MSEs at different training and test percentages to determine the best combination of parameters.

While calculating the MSE values, the network was run 100 times for each distribution, and randomness was kept by calculating the average values. The 80/20 training/test split gave the best forecasting performance in our case. In this way, the split of test and training data that gives the lowest error rate were successfully determined. We adopted MSE as a performance measurement technique in different studies based on ANN [[113], [114], [115], [116], [117]]. Then, the dataset is divided into training and test sets. This division is performed considering the first 80% of the initial data for the training step and 20% of the last data for the testing step.

After determining the best network structure for the forecasting problem under consideration, we also used different time lags between the predictor data and the target variable. Additionally, we have added the last value of the daily number of patients recovering from COVID-19 as network inputs to better forecast performance. Table 8 summarizes the results of different input variable sets for the determined neural network structure. While calculating the forecast KPIs, the network was run 100 times for each model, and randomness was preserved by calculating average values.

Table 8

The results of MLP models.

Inputs	MSE	RMSE	MAPE	MAE
7 Day Lag	27118	165	9.063	131.554
7 Day Lag with Last Value	25876	160	7.139	112.993
10 Day Lag	60195	245	13.245	187.785
10 Day Lag with Last Value	225267	474	20.342	307.185
14 Day Lag	331131	575	35.681	484.2
14 Day Lag with Last Value	567360	753	38.376	497.279

After the forecasting process was performed, the best model was determined as 7-Day Lag with Last Value with respect to all KPIs. This model uses both 7-day lag and the last value of output to make the forecast. Fig. 6 shows the R-squared (R2) values of the proposed network structure and the predicted values of daily recovered patients; in other words, the ANN model outputs against the desired targets.

Fig. 6

R2 values for the training and tests.

R2 is one of the most used statistical measures in forecasting models, which determines the proportion of the variance for a predicted output explained by input variables in a forecasting model. R2 evaluates the fit quality of the proposed model [94]. The formula to calculate R2 is expressed as follows:where shows the predicted value, and yj shows the true value (j = 1,2, …,m).

When considering the relationship of all four predictors, the ANN model accurately explained 95.89% of the variation in the training data and 97.50% on the test data. The plots of dispersion of actual values and predicted values of daily recovered people are presented in Fig. 6 for training and test data sets. In Fig. 6, the two plots represent the training and testing data. The dashed lines show the perfect results, namely outputs = targets. The solid blue and red lines show the best fit linear regression line between outputs and targets for training and testing data, respectively. For this study, both the training and testing data indicate good fits. A small quantity number for outliers in the plots means that a perfect matching of the proposed model to the actual data. This model uses predictors and generates more complex equations to interpret them simply. Fig. 7 shows a comparison of the predicted and actual data for the test stage.

Fig. 7

The predicted and actual data for daily recovered patients in Turkey.

As can be seen from Fig. 7, there is not much difference between forecast and actual data. The biggest error was calculated for November 21, 2020. The forecast value was 3816 on the day November 21, 2020, when 4485 people were recovered. As shown in Fig. 7, there was a sudden increase in the number of recovered patients after November 6, 2020. The proposed network model has forecasted this change well, as shown, with an error rate of about 5%. Eventually, it can be inferred that the recovered patients forecasted for each day will be ready to donate plasma.

After that, we forecasted the next 7 days (December 5 to 11) using the past 7 days' data in the proposed MLP structure. For example, to forecast the number of daily recovered people on December 5, we used the value of input variables on November 28. Thus, a 7-day forecast was performed. Fig. 8 demonstrates both the actual data for the last 3 days (from December 2 to December 4) and the forecasted data for the next 7 days from December 5 to December 11, 2020.

Fig. 8

Forecasting of next 7 days.

In this study, apart from forecasting only those who died or became infected by COVID-19, the daily number of recovered patients has been forecasted, and a perspective that would contribute to the treatment process has been presented. Therefore, it has an original and different perspective as compared to other COVID-19 forecasting studies.

CIP therapy is currently the only known treatment method that contributes to treating severe COVID-19 cases. As a result of the proposed approach, appropriate strategies for CIP therapy will be developed with the estimated number of recovered patients, and more people can survive until the vaccine is available for COVID-19. With estimating the number of recovered patients, the number of people who can donate plasma will be determined to realize the correct planning and timely implementation of plasma donation. This proposed forecasting mechanism is verified for Turkey. More planned steps can be taken regarding the treatment of the disease in all countries of the world struggling with this disease, and thus an effective blood supply chain system can be developed.

Comparison and testing procedure for the proposed algorithm

For testing and comparing the proposed network structure, the number of deaths from COVID-19 across the country was forecasted. For this aim, the adopted methodology had been run from scratch. Firstly, the correlation coefficients of the candidate input variables with the number of daily death from COVID-19 were calculated using Equation (1) and given in Table 9 .

Table 9

Correlation coefficients.

	Input Variables	Correlation Coefficient
I-1	Total tests	0.677
I-2	Total cases	0.558
I-3	Total deaths	0.590
I-4	Intensive cares	0.830
I-5	Intubated patients	0.922
I-6	Total recovered	0.436
I-7	Daily tests	0.692
I-8	Daily cases	0.943
I-9	Daily recovered	0.489
I-10	Daily change in the intensive cares	0.249
I-11	Daily change in the intubated patients	0.605

After the determination of correlation coefficients between all variables with the number of daily death from COVID-19, “I-1 Total tests”, “I-2 Total cases”, “I-3 Total deaths”, “I-4 Intensive cares”, “I-5 Intubated patients”, “I-7 Daily tests”, “I-8 Daily cases” and “I-11 Daily change in the intubated patients” were determined as potential input variables in forecasting calculations. After that, correlation coefficients among the remaining 8 input variables were calculated as given in Table 10 .

Table 10

Correlation coefficients among the determined input variables.

Variable	I-1	I-2	I-3	I-4	I-5	I-7	I-8	I-11
I-1	1	0.9729	0.9759	0.7968	0.8666	0.9737	0.5164	0.6023
I-2	0.9729	1	0.9945	0.7138	0.8056	0.9503	0.4146	0.5127
I-3	0.9759	0.9945	1	0.7174	0.8330	0.9438	0.4506	0.5311
I-4	0.7968	0.7138	0.7174	1	0.8428	0.8127	0.7199	0.4403
I-5	0.8666	0.8056	0.8330	0.8428	1	0.8588	0.8379	0.6453
I-7	0.9737	0.9503	0.9438	0.8127	0.8588	1	0.5360	0.5724
I-8	0.5164	0.4146	0.4506	0.7199	0.8379	0.5360	1	0.5208
I-11	0.6023	0.5127	0.5311	0.4403	0.6453	0.5724	0.5208	1

When we analyze Table 10, it is noticed that there is a solid relationship between the “I-1 Total tests”, “I-2 Total cases”, “I-3 Total deaths” and I-7 Daily tests”. Having such a high correlation coefficient means that the values of these four variables increase or decrease together. Hence, it can be reasonable to use one of them as the predictor. So we used only “I-1 Total tests” as a predictor among these four variables. After all the data cleaning processes, we determined five different input variables as predictors such as “I-1 Total tests”, “I-4 Intensive cares”, “I-5 Intubated patients”, “I-8 Daily cases”, and “I-11 Daily change in the intubated patients”.

Once the predictors to be used and the values of these predictors were determined, these values were used in the same ANN structure used to forecast the number of daily recovered patients as inputs in forecasting the number of daily death from COVID-19 in Turkey.

When calculating the MSE values, the network structures were run 100 times for each combination, and the average MSE values were calculated to eliminate the chance factor to ensure that the lowest MSE is achieved. Fig. 9 shows the comparison of the forecasted and actual data for the test stage.

Fig. 9

The forecasted and actual data for daily deaths.

As shown in Fig. 9, there is not much difference between the forecasted and actual data. The biggest error is calculated for November 20. The forecasted value is 185, and 177 people died on November 20. The MAPE is calculated as 5.967 for the test stage.

As a result of this analysis, it was seen that our forecasting model could be used with different data, and successful results were achieved by testing it with the number of deaths. Our proposed network structure also gave a very low MAPE value in forecasting the number of deaths from COVID-19.

Robustness analysis

To confirm the robustness of the model presented, we analyzed the COVID-19 data of Italy and forecasted the number of recovered people. For this aim, we used data from February 24, 2020, to January 10, 2021, meaning 322 data points were used in the forecasting process. 77 of this data were used to test the network. The proposed ANN model had made very successful forecasting results with the data determined for Italy. The MAPE was calculated to be only 12.8. MSE and RMSE values were determined to be 1500 and 38.73, respectively. Fig. 10 shows the actual and predicted data for the test stages.

Fig. 10

The forecasted and actual data of recovered people from COVID-19 in Italy.

As can be seen from Fig. 10, the proposed forecasting model can be used for different countries. The biggest error was calculated for November 10. The forecast value was 12941, while 17734 people were recovered. As shown in Fig. 10, there were sudden changes, up and down, in the number of patients recovering after the first day of November. The proposed network model had forecasted this change well, as it is shown.

As a result of the robustness analysis, the proposed forecasting model has again yielded successful results. It has been proven to be a model that can be used to determine strategies to combat COVID-19 and the flawless operation of the blood supply chain mechanism.

Comparative analysis

In this study, we forecasted daily recovered people from COVID-19 by MLP structure, unlike the current studies in the literature that use time series based forecasting methodologies, such as ARIMA, LSTM, NARX. In this section, the proposed MLP based forecasting methodology is compared with different forecasting methods to demonstrate its robustness. For this purpose, ARIMA (a time series-based method, and LSTM and NARX, which are two machine learning methods, are utilized in this study to forecast daily recovered people from COVID-19 in Turkey. The same dataset, designed after the predictor selection step, is used to compare the results with the proposed MLP application. The results were compared according to their KPIs. This section gives brief information and results about the methods employed in this study.

ARIMA

ARIMA is one of the most used techniques for parametric univariate time series modeling. ARIMA models are applied to non-stationary series but can also be converted to stationary series by difference-taking. ARIMA strongly assumes that future data values are linearly dependent on current and past data values, similar to other linear methods [118]. In this way, ARIMA provides highly accurate results in stationary time series forecasting. The ARIMA method uses autoregressive (AR) and moving average (MA) models. AR includes lagged terms, and MA includes lagged terms on the residuals or noise [119]. ARIMA is used stationary time series data with no missing values. So time series data can be modeled as stationary or can be transformed into stationary by differencing. The Dickey-Fuller statistics are used to test whether the series is stationary [120]. Thus, the letter ‘I’ (Integrated) in ARIMA means that the first-order difference is applied to transform time series to stationary. The general representation of the models is ARIMA (p, d, q). Here, p and q are the degrees of the AR model and MA model, respectively, and d is the degree of difference. The equation representing the ARIMA model for the time sequence is given in Equation (13). is a normal random variable white noise sequence with zero mean and variance σ2 and B is the backshift operator whose effect on the can be represented as .

In this study, the Box-Jenkins approach is utilized in ARIMA models to find the most suitable time series model for the past values of the time series. The preliminary analysis of the daily recovered people data is performed using time plots of both the original time series and 1-differencing, as shown by Fig. 11 .

Fig. 11

Original series and 1-differencing.

The time series used in the study was non-stationary according to the result of the Dickey-Fuller test with a test result of a p-value of 0.256, at the significance level of 0.05. The data reached stationarity with one order of differencing, as can be seen in Fig. 11. To construct the ARIMA model, the non-stationary data was transformed stationary by the 1-differencing method. Before further analysis using the Box-Jenkins approach, the data were transformed to achieve stationary time series data. After applying the 1-differencing method, the time series data became stationary. So d was determined as “1” in this study.

The next step in the Box Jenkins approach was to identify whether the model needs any AR terms or not. A partial Autocorrelation plot can be helpful to determine the required number of AR terms, namely the value of p. Any autocorrelation in a stationary time series can be revised by adding enough p. So, the order of p was taken to be equal to many lags that cross the significance limit in the partial autocorrelation plot. Fig. 12 shows a partial autocorrelation plot of stationary time series.

Fig. 12

Partial autocorrelation plot.

Lag 1 crossed the significance limit (blue region), as can be seen in Fig. 12. So lag 1 was quite significant. Therefore, p could be fixed 1 in this time series.

An autocorrelation plot can be used to determine the number of MA terms, namely q. Fig. 13 shows the autocorrelation plot of stationary time series.

Fig. 13

Autocorrelation plot.

Again, lag 1 crossed the significance limit (blue region), so q could be fixed as 1. Therefore, according to the Box Jenkins approach, the most suitable model was determined as ARIMA (1, 1, 1). The results of an ARIMA (1,1,1) model were estimated by maximum likelihood estimation and given in Table 11 .

Table 11

Parameter estimates of ARIMA (1, 1, 1) model.

Variable	Coefficient	Standard Error	z-statistic	p-value
AR(1)	0.4613	0.884	0.522	0.602
MA(1)	−0.4295	0.899	−0.478	0.633
Constant	15.7149	20.679	0.760	0.447

Akaike's Information Criterion (AIC) and Bayesian Information Criterion (BIC) values are determined as 3492 and 3506, respectively.

To achieve the best performance, different values of “p” and “q” are tried. In this study, the ARIMA models are implemented in Python ARIMA. The results of the ARIMA for forecasting the daily recovered patients are presented in Table 12 . The best model with respect to MAPE is determined as (1,1,1) model. To make a fair comparison of the different ARIMA models, MSE, RMSE, MAPE, and MAE are calculated.

Table 12

The results of ARIMA models.

(p,d,q)	MSE	RMSE	MAPE	MAE	(p,d,q)	MSE	RMSE	MAPE	MAE
(0,1,0)	92848	305	12.579	185.736	(2,1,0)	92654	304	12.501	185.595
(0,1,1)	92795	305	12.595	185.790	(2,1,1)	92643	304	12.267	185.449
(0,1,2)	92631	304	12.572	185.772	(2,1,2)	85867	293	14.051	187.560
(0,1,3)	92488	304	12.556	185.317	(2,1,3)	89114	299	12.696	185.313
(1,1,0)	92791	305	12.578	185.787	(3,1,0)	92548	304	12.445	185.159
(1,1,1)	92724	305	12.223	185.498	(3,1,1)	92259	304	12.633	185.857
(1,1,2)	92612	304	12.292	185.644	(3,1,2)	88995	298	12.764	185.432
(1,1,3)	92200	304	12.731	186.211	(3,1,3)	88747	298	12.762	184.982

LSTM

Long short-term memory networks, often referred to as LSTM, are a special type of RNN that can learn long-term dependencies [121]. It is aimed to store and transfer the state information of the ANN while processing data in RNNs. However, as a result of continuous processing and transfer of state information, it is unlikely to be transferred without breaking long-term dependencies [122]. In other words, while short-term dependencies are transferred very successfully, sometimes long-term dependencies are not transferred effectively. LSTMs are designed to address long-term dependency issues. LSTM provides faster convergence for the training data and can identify long-term dependencies in the input data [123]. The LSTM module consists of 3 separate gates: forget gate, input gate, and output gate. Forget gate determines how much of the information is forgotten and how much is passed on to the next stage. In the next step, the input gate controls which information should be stored. Finally, the output gate determines which information to read and output. The equations for the forget gate, input gate, and output gate are given in Equations (14), (15), (16), respectively [124]. Where be the sigmoid function, are relevant weights in forget, input and output gate, denotes the previous output at time , is the input vector at time , and are bias neurons at the respective gate associated with each LSTM block. In this study, the sigmoid function is used as the transfer function. Then, Equation (17) is given for cell state:

and represent the new and previous cell states. means the element-wise multiplication of the vectors. The equations for candidate cell state are given:

In this study, tanh is used as an activation function. The LSTM is implemented in Python by using the Keras library with Tensorflow backend. We have the data of daily recovered patients in Turkey from March 28, 2020, to November 27, 2020. The first 80% of the time-series data is used to train the network, and the remaining data is used to test the network. Namely, 196 data is used for training, and 49 data is used for testing. Fig. 14 shows the training and test data used in this study.

Fig. 14

Input data of LSTM network.

Several parameters were considered in this study. 5, 10, 15, and 20 hidden neurons were tested in the hidden layer, and the time step of the input sequence (lag) was analyzed between 1 and 14 to achieve better performance, that means the number of daily recovered people in the next day was forecasted using past data of n time step. The “Adam” optimizer was used, and batch size was set to 4. The number of epochs was determined as 100 in the LSTM structure. The loss function was MSE in the LSTM application. Then forecast KPIs of each trial were calculated. The best performance was determined as 11 lags with respect to all performance measures. The results of LSTM implementation are presented in Table 13 .

Table 13

The results of LSTM models.

Neuron	lag	MSE	RMSE	MAPE	MAE	Neuron	lag	MSE	RMSE	MAPE	MAE
5	1	56572	238	7.971	177.638	15	1	55754	236	7.894	175.740
5	2	51744	227	7.349	168.283	15	2	48239	220	7.286	166.117
5	3	82729	288	9.619	220.673	15	3	59089	243	7.987	181.683
5	4	197461	444	17.096	377.443	15	4	78560	280	9.289	215.793
5	5	49644	223	7.235	164.462	15	5	238811	489	17.360	401.676
5	6	56641	238	7.667	179.532	15	6	52072	228	7.207	165.488
5	7	64499	254	8.208	193.347	15	7	91028	302	9.800	234.754
5	8	45494	213	7.423	156.764	15	8	150036	387	12.193	299.410
5	9	100120	316	10.768	232.969	15	9	270676	520	19.087	410.643
5	10	358982	599	14.415	374.843	15	10	248754	499	16.302	403.890
5	11	168621	411	14.747	339.102	15	11	366927	606	19.189	470.035
5	12	65727	256	8.143	200.108	15	12	233278	483	13.260	345.633
5	13	371776	610	19.705	502.980	15	13	292434	541	16.486	422.905
5	14	121581	349	10.931	276.037	15	14	369737	608	16.818	449.302
10	1	45461	213	8.856	161.969	20	1	57070	239	8.020	178.793
10	2	46816	216	8.952	167.084	20	2	50316	224	7.289	167.332
10	3	122655	350	14.59	279.596	20	3	70514	266	9.224	204.683
10	4	48323	220	8.822	167.67	20	4	83512	289	9.943	227.505
10	5	57832	240	9.505	177.153	20	5	81765	286	9.722	218.089
10	6	62834	251	9.351	177.218	20	6	99130	315	9.958	221.806
10	7	41645	204	8.136	154.659	20	7	108665	330	9.484	224.204
10	8	49748	223	8.626	162.764	20	8	125848	355	13.068	289.266
10	9	47921	219	8.776	170.753	20	9	151110	389	13.340	314.699
10	10	63552	252	10.018	197.963	20	10	84775	291	9.522	218.877
10	11	42519	206	7.807	149.97	20	11	128943	359	10.230	262.864
10	12	55855	236	9.031	180.647	20	12	80277	283	9.397	205.864
10	13	53383	231	9.031	181.973	20	13	167281	409	12.355	317.922
10	14	69337	263	9.152	184.405	20	14	301992	550	18.122	456.975

NARX

NARX networks are a kind of dynamic neural network model that gives successful results in nonlinear system modeling and time-series forecasting applications. Compared to traditional neural networks, NARX networks converge faster and demonstrate more effective learning. The values of the output signals depend on both the input signals and the historical values of output signals in a dynamic system [125]. Thus, NARX can provide more effective results than traditional neural networks [126]. NARX model for the time sequence is given in Equation (20). Let be the inputs, f be the nonlinear ambiguity function, and lastly d be the feedback.

The NARX was implemented in MATLAB R2020b using the “narx net” structure in this study. Daily recovered patients were used as the target variable, “I-4 Intensive cares”, “I-9 Daily deaths”, “I-10 Daily change in the intensive cares”, and “I-11 Daily change in the intubated patients” were used as the input variables. The tangent sigmoid transfer function and linear transfer functions were used for the hidden layer and the output neuron. The number of epochs was determined as 1000 in the NARX structure. The dataset was divided into training and test sets. This division was performed considering 80% of the initial data for the training step and the remaining 20% of the data for the testing step. Additionally, the delay was analyzed between 1 and 7 with different numbers of neurons to achieve better performance. The best performance is achieved for 10 neurons with 5 delays. The results of NARX are presented in Table 14 .

Table 14

The results of NARX models.

Neuron	Delay	MSE	RMSE	MAPE	MAE	Neuron	Delay	MSE	RMSE	MAPE	MAE
3	1	67942	261	13.334	164.800	7	1	86169	294	27.684	193.713
3	2	82851	288	17.843	188.910	7	2	115212	339	33.117	244.976
3	3	84870	291	14.715	187.677	7	3	58795	242	14.722	167.026
3	4	54929	234	13.890	147.664	7	4	63246	251	14.444	162.860
3	5	66299	257	14.155	173.820	7	5	76667	277	11.512	152.969
3	6	79779	282	17.251	200.794	7	6	119339	345	21.878	240.404
3	7	246981	497	23.515	313.483	7	7	108467	329	20.154	243.278
4	1	76025	276	18.575	192.396	8	1	71621	268	15.176	207.234
4	2	70606	266	13.824	175.397	8	2	66490	258	14.280	161.866
4	3	79649	282	14.052	178.730	8	3	60521	246	13.271	166.725
4	4	75539	275	15.828	189.546	8	4	80732	284	19.862	206.159
4	5	66019	257	16.383	186.433	8	5	46325	215	14.685	140.594
4	6	76218	276	16.056	198.327	8	6	90021	300	21.696	230.516
4	7	82728	288	14.551	197.595	8	7	42322	206	10.293	139.212
5	1	75272	274	14.607	192.396	9	1	74267	273	20.446	175.982
5	2	97525	312	26.906	204.123	9	2	64569	254	14.494	150.948
5	3	59589	244	13.540	163.058	9	3	61212	247	13.261	160.839
5	4	77601	279	18.307	194.039	9	4	114996	339	19.474	239.338
5	5	65230	255	15.994	174.969	9	5	119551	346	25.606	260.760
5	6	61764	249	14.511	162.376	9	6	178826	423	25.254	323.668
5	7	86714	294	13.433	191.343	9	7	64268	254	10.874	157.449
6	1	80379	284	19.120	178.416	10	1	62522	250	11.805	175.982
6	2	75049	274	18.837	184.113	10	2	68830	262	16.206	180.345
6	3	120091	347	32.092	264.169	10	3	68581	262	15.155	170.314
6	4	59897	245	14.039	160.584	10	4	73782	272	15.579	178.332
6	5	111323	334	26.986	257.046	10	5	51834	228	10.986	146.103
6	6	70788	266	16.536	187.946	10	6	68087	261	14.778	170.311
6	7	161004	401	21.611	286.878	10	7	134816	367	20.565	258.958

Table 15 summarizes the investigated parameters for each method. 16 different ARIMA models, 56 different LSTM structures, and 56 different NARX models were utilized to compare the proposed methodology results.

Table 15

Investigated parameters for each method.

Method	Parameter	Values
ARIMA	p	0, 1, 2, 3
	d	1
	q	0, 1, 2, 3
LSTM	# of neurons	5, 10, 15, 20
	lag	1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14
NARX	# of neurons	3, 4, 5, 6, 7, 8, 9, 10
	lag	1, 2, 3, 4, 5, 6, 7
MLP	# of layers	1, 2
	# of neurons	4, 5, 6, 7, 8, 9, 10, 11, 12
	learning rate	starting from 0.1 increased by 0.05 up to 1
	momentum coefficient	starting from 0.1 increased by 0.05 up to 1
	transfer function	tangent sigmoid, logarithmic sigmoid

After four different models were applied to forecast daily recovered patients, their best performance results were compared and evaluated according to forecasting KPIs as given in Table 16 . For the LSTM and ARIMA applications, three different models could be determined as the best models in terms of different KPIs. In the NARX application, one model was determined as the best structure with respect to each KPI.

Table 16

The comparison of forecasting methods.

Method	Structure	MSE	RMSE	MAPE	MAE
MLP	Two layers: 8 neurons, 5 neurons	25876	160	7.139	112.993
LSTM	10 neurons, 11 lags	42519	206	7.807	149.970
LSTM	10 neurons, 7 lags	41645	204	8.136	154.659
LSTM	15 neurons, 6 lags	52072	228	7.207	165.488
NARX	8 neurons, 7 delays	42322	206	10.293	139.212
ARIMA	(1,1,1)	92724	305	12.223	186.651
ARIMA	(2,1,2)	85867	293	14.051	187.560
ARIMA	(3,1,3)	88747	298	12.762	184.982

As shown in Table 16, the proposed MLP model provides the best performance with respect to all KPIs. LSTM model with 15 hidden neurons and uses 7-day lag is determined as the second-best performance with the MAPE value of 7.207. However, ARIMA shows the worst performance with a 12.223 MAPE value as expected and all other KPIs. It can be said that the proposed forecasting method is very effective for forecasting daily recovered patients.

Conclusion

Nowadays, where COVID-19 still seriously affects the whole world, making forecasts about the number of cases and deaths is one of the most critical steps to take the necessary precautions and implement the right strategies to fight against the disease. However, besides these forecasts, it is crucial to estimate the number of patients recovering for CIP treatment, which is applied by taking plasma from recovered patients and transferring it to infected people. It is currently considered an effective treatment method in the fight against the disease. In this way, accurate planning can be done to accomplish blood supply between the infected people and donors, and it can create a smooth supply chain system for a perishable product like blood. The daily number of recovered people, who play a key role in planning plasma supply and distribution required for the implementation of CIP therapy is forecasted. This study introduces a dynamic forecasting model to the literature, which can be used daily, especially during emergencies. The proposed forecasting model is used with different datasets. With using this model, the number of people, who have been diagnosed with COVID-19, received treatment, and tested negative for COVID-19, is forecasted on a daily basis.

For this aim, ANN approach is used to forecast the daily number of people recovered from COVID-19 to determine the number of donors in the treatment of the disease. The variables used in the forecasting are determined as a result of expert opinions and a literature review. Subsequently, the variables are analyzed and eliminated to be used as predictors by correlation analysis and causal analysis. Thus, a novel strategy is proposed to increase forecast accuracy by systematically determining the predictors in forecasting models. After the predictors are determined, their data is compiled, and the best ANN network structure is determined. The most suitable structure is found by trying many different structures. Besides, the robustness of the proposed model is evaluated by analyzing data from Italy for the number of patients recovered from COVID-19 and the number of deaths. This helped to validate and test the proposed model. As a result of the study, it is shown that ANN is an effective approach to forecast the daily number of people recovered from COVID-19.

This study has a different point of view than other forecasting studies on COVID-19, and no similar study has been found in the literature to the best of the author's knowledge. In this regard, this study is thought to have a high contribution to the literature and practice, and it is believed that it would be a leading paper for the researchers and practitioners who want to study in this area. This paper can be a guide to establish a smooth and secure collection/distribution infrastructure on supplying plasma for giving a more effective response to the pandemic. It can also serve as a reference point for the uninterrupted flow of COVID-19 treatment and solutions that enable matching the highest number of donors with patients.

In the future, the approach adopted in this paper can be used with updated data or different forecasting methods can be applied to compare the current results. Other forecasting methods may be applied to find even better forecasting results, or easier to use methodologies would be developed apart from ANN. A comparison of different ANN-based forecasting models would also be interesting. Support vector regression may also improve the accuracy of the proposed methodology as used in many studies [[127], [128], [129], [130], [131]]. With the increase in available data, time series-based forecasting models can be developed.

The proposed forecasting approach could eventually help the decision-makers or managers in planning new pandemic strategies. Studying with a network trained on historical data of other pandemics can be helpful in the long term. To improve the developed predictor selection methodology, data collection and data selection methods based on big data using deep-learning technologies can be applied. This study can be improved with model selection for each input setting for time lags. Additionally, the training stage can be enhanced with a calculation of entropy in the input dataset. The forecasting methodology can be extended to other countries because of the same pandemic process and the treatment mechanism they need, including guidelines for predictor selection.

Funding

No funds, grants, or other support was received.

Ethical approval

Ethics committee approval is not required.

Availability of data and material

Not applicable.

Code availability

Not applicable.

Consent to participate

Not applicable.

Consent to publish

The authors confirm that the final version of the manuscript has been reviewed, approved, and consented to for publication by all authors.

Authors' contributions

Ertugrul Ayyildiz and Melike Erdogan developed the theoretical formulation, performed the analytic calculations, and performed the numerical analysis. Ertugrul Ayyildiz and Alev Taskin contributed to the sensitivity analysis and final version of the manuscript. Ertugrul Ayyildiz and Melike Erdogan drafted the manuscript. Alev Taskin supervised the work.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.