Statistical, machine learning, and deep learning models for COVID-19 forecasting in Kenya

ARIMA is a statistical model for forecasting univariate time series with autocorrelation and non-stationarity, integrating autoregressive (AR), differencing (I), and moving average (MA) components [13,19]. The AR component links current observations to past values, while the MA component models dependencies on past forecast errors. Differencing achieves stationarity. The ARIMA( p , d , q p,d,q ) model, where p p is the AR order, d d is the differencing order, and q q is the MA order, is expressed as follows: (1) Δ d Y t = μ + ∑ i = 1 p ϕ i Δ d Y t − i + ε t + ∑ j = 1 q θ j ε t − j , {\Delta }^{d}{Y}_{t}=\mu +\mathop{\sum }\limits_{i=1}^{p}{\phi }_{i}{\Delta }^{d}{Y}_{t-i}+{\varepsilon }_{t}+\mathop{\sum }\limits_{j=1}^{q}{\theta }_{j}{\varepsilon }_{t-j}, where μ \mu is a constant, ϕ i {\phi }_{i} are the AR coefficients, θ j {\theta }_{j} are the MA coefficients, and ε t {\varepsilon }_{t} is the white noise. The autocorrelation function (ACF) and partial autocorrelation function (PACF) guide p p and q q selection, with PACF cutoffs indicating AR processes and ACF cutoffs suggesting MA processes.

SVR is an extension of the SVM framework for regression tasks. It aims to find a function f ( x ) f\left(x) that deviates from the true target values y i {y}_{i} by no more than a small margin ε \varepsilon , while keeping the model as flat as possible [3,11]. The function is defined as follows: (2) f ( x ) = ⟨ w , x ⟩ + b , f\left(x)=\langle w,x\rangle +b, where w w is the weight vector, x x is the input vector, and b b is the bias term. To enforce flatness, the norm ‖ w ‖ 2 \Vert w{\Vert }^{2} is minimized. This leads to the following primal optimization problem: (3) min w , b 1 2 ‖ w ‖ 2 . \mathop{\min }\limits_{w,b}\frac{1}{2}\Vert w{\Vert }^{2}. subject to the constraint: (4) ∣ y i − ⟨ w , x i ⟩ − b ∣ ≤ ε , ∀ i = 1 , … , n . | {y}_{i}-\langle w,{x}_{i}\rangle -b| \le \varepsilon ,\hspace{1.0em}\forall i=1,\ldots ,n.

To allow for infeasible constraints due to noise or outliers, slack variables ξ i , ξ i * ≥ 0 {\xi }_{i},{\xi }_{i}^{* }\ge 0 are introduced. The problem becomes (5) min w , b , ξ i , ξ i * 1 2 ‖ w ‖ 2 + C ∑ i = 1 n ( ξ i + ξ i * ) \mathop{\min }\limits_{w,b,{\xi }_{i},{\xi }_{i}^{* }}\frac{1}{2}\Vert w{\Vert }^{2}+C\mathop{\sum }\limits_{i=1}^{n}\left({\xi }_{i}+{\xi }_{i}^{* }) subject to (6) y i − ⟨ w , x i ⟩ − b ≤ ε + ξ i , {y}_{i}-\langle w,{x}_{i}\rangle -b\le \varepsilon +{\xi }_{i}, (7) ⟨ w , x i ⟩ + b − y i ≤ ε + ξ i * , \langle w,{x}_{i}\rangle +b-{y}_{i}\le \varepsilon +{\xi }_{i}^{* }, (8) ξ i , ξ i * ≥ 0 . {\xi }_{i},{\xi }_{i}^{* }\ge 0. Here, C > 0 C\gt 0 is a regularization parameter that balances the model complexity and the amount up to which deviations larger than ε \varepsilon are tolerated.

Several commonly used kernels in SVR include the following:

Linear kernel: (9) K ( x , x ′ ) = ⟨ x , x ′ ⟩ , K\left(x,x^{\prime} )=\langle x,x^{\prime} \rangle ,

Polynomial kernel: (10) K ( x , x ′ ) = ( ⟨ x , x ′ ⟩ + c ) d , K\left(x,x^{\prime} )={\left(\langle x,x^{\prime} \rangle +c)}^{d},

Radial basis function (RBF) kernel: (11) K ( x , x ′ ) = exp ( − γ ‖ x − x ′ ‖ 2 ) , K\left(x,x^{\prime} )=\exp \left(-\gamma \Vert x-x^{\prime} {\Vert }^{2}), where c c is a constant, d d is the degree of the polynomial, and γ \gamma controls the width of the RBF kernel.

Regularization through the parameter C C controls the trade-off between model complexity and tolerance to errors. A large C C assigns a higher penalty to errors, leading to less margin violation, while a small C C allows for a more generalized model by tolerating larger deviations. The choice of kernel and regularization parameters significantly affects the SVR model’s ability to capture nonlinear relationships and its generalization performance.

RF is an ensemble learning algorithm that builds multiple decision trees and aggregates their outputs to perform regression. It was introduced by Breiman [5] and is particularly effective in handling high-dimensional data, nonlinear relationships, and reducing overfitting.

In regression settings, the goal is to predict a continuous target variable y ∈ R y\in {\mathbb{R}} based on an input vector x ∈ R d {\bf{x}}\in {{\mathbb{R}}}^{d} . RF achieves this by constructing a collection of T T decision trees, each trained independently on a different subset of the training data. The final prediction for a new input x {\bf{x}} is obtained by averaging the predictions from all individual trees: (12) y ˆ = 1 T ∑ t = 1 T h t ( x ) , \hat{y}=\frac{1}{T}\mathop{\sum }\limits_{t=1}^{T}{h}_{t}\left({\bf{x}}), where h t ( x ) {h}_{t}\left({\bf{x}}) denotes the prediction made by the t t th decision tree.

Each tree is trained on a bootstrap sample, which is created by randomly sampling (with replacement) from the original training dataset. This introduces variability among the trees. Additionally, during the tree-building process, RF introduces further randomness by selecting a random subset of features at each split rather than considering all features. This decorrelates the trees and enhances generalization.

At each node of a decision tree, the algorithm selects the best split point among the randomly chosen features based on a splitting criterion that minimizes prediction error. For regression, the most commonly used criterion is the mean-squared error between predicted and actual values in the resulting child nodes. This local optimization helps guide the recursive partitioning of the input space.

Unlike a single decision tree, which may overfit the training data, RF reduces the variance of predictions by averaging across many trees. The ensemble effect ensures that individual overfitting errors are averaged out, leading to more robust and accurate predictions.

RNNs are a class of neural architectures tailored for sequential data modeling. They maintain a hidden state h t {h}_{t} that evolves over time by incorporating information from both the current input x t {x}_{t} and the previous hidden state h t − 1 {h}_{t-1} , thereby capturing temporal dependencies [26]. The hidden state is updated using the recurrence relation: (13) h t = σ ( W h h t − 1 + W x x t + b h ) , {h}_{t}=\sigma \left({W}_{h}{h}_{t-1}+{W}_{x}{x}_{t}+{b}_{h}), where W h {W}_{h} and W x {W}_{x} are the weight matrices corresponding to the hidden state and the input, respectively, b h {b}_{h} is the bias term, and σ \sigma is a nonlinear activation function such as tanh or ReLU. The output at each time step y t {y}_{t} is computed as follows: (14) y t = W y h t + b y , {y}_{t}={W}_{y}{h}_{t}+{b}_{y}, where W y {W}_{y} is the weight matrix for the output, and b y {b}_{y} is the output bias term. This recurrent structure, combined with the output transformation, allows the network to retain and propagate information through time steps, making it suitable for tasks involving time series or sequences.

LSTM is a specialized type of RNN designed to capture long-range dependencies in sequential data by using gates to control the flow of information. Each LSTM cell contains three primary gates: the forget gate, input gate, and output gate, which regulate how information is retained, updated, and outputted at each time step [6].

The forget gate f t {f}_{t} determines what proportion of the previous cell state C t − 1 {C}_{t-1} should be forgotten. It outputs values between 0 and 1, where 0 means “completely forget” and 1 means “completely retain.” The forget gate is calculated as follows: (15) f t = σ ( W f x t + U f h t − 1 + b f ) , {f}_{t}=\sigma \left({W}_{f}{x}_{t}+{U}_{f}{h}_{t-1}+{b}_{f}), where f t {f}_{t} is the forget gate at time t t , W f {W}_{f} , U f {U}_{f} , and b f {b}_{f} are the weights and bias for the forget gate, and σ \sigma is the sigmoid activation function.

The input gate i t {i}_{t} determines how much new information should be added to the cell state C t {C}_{t} . It also uses a sigmoid function to output values between 0 and 1. The input gate is computed as follows: (16) i t = σ ( W i x t + U i h t − 1 + b i ) , {i}_{t}=\sigma \left({W}_{i}{x}_{t}+{U}_{i}{h}_{t-1}+{b}_{i}), where i t {i}_{t} is the input gate at time t t , W i {W}_{i} , U i {U}_{i} , and b i {b}_{i} are the weights and bias for the input gate.

The cell candidate C t ˜ \tilde{{C}_{t}} represents the new information that is to be added to the cell state. This candidate value is then filtered by the input gate i t {i}_{t} . The cell state C t {C}_{t} is updated as follows: (17) C t = f t C t − 1 + i t C t ˜ , {C}_{t}={f}_{t}{C}_{t-1}+{i}_{t}\tilde{{C}_{t}}, where C t − 1 {C}_{t-1} is the previous cell state, f t {f}_{t} is the forget gate, i t {i}_{t} is the input gate, and C t ˜ \tilde{{C}_{t}} is the cell candidate.

The output gate o t {o}_{t} determines what part of the cell state should be output as the hidden state h t {h}_{t} . It acts as a filter to decide how much of the cell state should contribute to the next hidden state. The output gate is calculated as follows: (18) o t = σ ( W o x t + U o h t − 1 + b o ) . {o}_{t}=\sigma \left({W}_{o}{x}_{t}+{U}_{o}{h}_{t-1}+{b}_{o}).

Finally, the hidden state h t {h}_{t} is updated as (19) h t = o t tanh ( C t ) , {h}_{t}={o}_{t}\tanh \left({C}_{t}), where o t {o}_{t} is the output gate, C t {C}_{t} is the cell state, and tanh \tanh is the hyperbolic tangent activation function, applied to the cell state.

Figure 1 presents the LSTM network, which incorporates memory cells to handle long-term dependencies effectively.

Figure 1

Architecture of the LSTM network. Source: From ref. [39].

The GRU is a simpler version of LSTM with fewer parameters. It combines the forget and input gates of LSTM into a single update gate. The GRU has two main gates: the update gate and the reset gate [39].

The update gate z t {z}_{t} controls the extent to which the previous hidden state h t − 1 {h}_{t-1} should be retained. A higher value means more of the previous state is retained, while a lower value means more of the current candidate hidden state h t ˜ \tilde{{h}_{t}} is used. The update gate is computed as (20) z t = σ ( W z x t + U z h t − 1 + b z ) , {z}_{t}=\sigma \left({W}_{z}{x}_{t}+{U}_{z}{h}_{t-1}+{b}_{z}), where z t {z}_{t} is the update gate, and W z {W}_{z} , U z {U}_{z} , and b z {b}_{z} are the weights and bias for the update gate.

The reset gate r t {r}_{t} determines how much of the previous hidden state should be used when computing the candidate hidden state. It is calculated as follows: (21) r t = σ ( W r x t + U r h t − 1 + b r ) , {r}_{t}=\sigma \left({W}_{r}{x}_{t}+{U}_{r}{h}_{t-1}+{b}_{r}), where r t {r}_{t} is the reset gate, and W r {W}_{r} , U r {U}_{r} , and b r {b}_{r} are the weights and bias for the reset gate.

The candidate hidden state h t ˜ \tilde{{h}_{t}} is computed using the reset gate to decide how much of the previous hidden state should be incorporated. It is calculated as follows: (22) h t ˜ = tanh ( W h x t + U h ( r t ⋅ h t − 1 ) + b h ) , \tilde{{h}_{t}}=\tanh \left({W}_{h}{x}_{t}+{U}_{h}\left({r}_{t}\cdot {h}_{t-1})+{b}_{h}), where h t ˜ \tilde{{h}_{t}} is the candidate hidden state and r t ⋅ h t − 1 {r}_{t}\cdot {h}_{t-1} is the element-wise multiplication of the reset gate and the previous hidden state.

Finally, the hidden state h t {h}_{t} is updated as a linear interpolation between the previous hidden state h t − 1 {h}_{t-1} and the candidate hidden state h t ˜ \tilde{{h}_{t}} , controlled by the update gate z t {z}_{t} : (23) h t = ( 1 − z t ) ⋅ h t − 1 + z t ⋅ h t ˜ , {h}_{t}=\left(1-{z}_{t})\cdot {h}_{t-1}+{z}_{t}\cdot \tilde{{h}_{t}}, where h t − 1 {h}_{t-1} is the previous hidden state, h t ˜ \tilde{{h}_{t}} is the candidate hidden state, and z t {z}_{t} is the update gate.

The structure of the GRU is illustrated in Figure 2, where gating mechanisms help regulate information flow and improve efficiency.

Figure 2

Architecture of the GRU. Source: From ref. [39].

The predictive performance of the models is evaluated using four standard metrics: RMSE, MAE, MAPE, and R 2 {R}^{2} . These metrics assess the accuracy of point forecasts in both absolute and relative terms. Additionally, to quantify the improvement of the best-performing model (RF) compared to other models, the percentage improvement for RMSE, MAE, and MAPE is calculated, providing a relative measure of performance enhancement: (24) RMSE = 1 n ∑ i = 1 n ( y i − y ˆ i ) 2 , \hspace{0.1em}\text{RMSE}\hspace{0.1em}=\sqrt{\frac{1}{n}\mathop{\sum }\limits_{i=1}^{n}{\left({y}_{i}-{\hat{y}}_{i})}^{2}}, (25) MAE = 1 n ∑ i = 1 n ∣ y i − y ˆ i ∣ , \hspace{0.1em}\text{MAE}\hspace{0.1em}=\frac{1}{n}\mathop{\sum }\limits_{i=1}^{n}| {y}_{i}-{\hat{y}}_{i}| , (26) MAPE = 100 n ∑ i = 1 n y i − y ˆ i y i , \hspace{0.1em}\text{MAPE}\hspace{0.1em}=\frac{100}{n}\mathop{\sum }\limits_{i=1}^{n}\left|\frac{{y}_{i}-{\hat{y}}_{i}}{{y}_{i}}\right|, (27) R 2 = 1 − ∑ i = 1 n ( y i − y ˆ i ) 2 ∑ i = 1 n ( y i − y ¯ ) 2 , {R}^{2}=1-\frac{{\sum }_{i=1}^{n}{\left({y}_{i}-{\hat{y}}_{i})}^{2}}{{\sum }_{i=1}^{n}{\left({y}_{i}-\bar{y})}^{2}}, (28) Percentage improvement = M other − M RF M other × 100 , \hspace{0.1em}\text{Percentage improvement}=\frac{{M}_{\text{other}}-{M}_{\text{RF}}}{{M}_{\text{other}}}\times 100, where M other {M}_{\text{other}} represents the RMSE, MAE, or MAPE of a competing model (ARIMA, SVR, RNN, LSTM, or GRU), and M RF {M}_{\text{RF}} is the corresponding metric for RF, the best-performing model.

To statistically evaluate whether the forecasting performance of two competing models differs significantly, the DM test is applied [8]. The test is based on the loss differential series d t {d}_{t} , defined as the difference in forecast losses from two models: (29) d t = L ( e 1 t ) − L ( e 2 t ) , {d}_{t}=L\left({e}_{1t})-L\left({e}_{2t}), where L ( ⋅ ) L\left(\cdot ) is a loss function such as MSE, MAE, and MAPE, and e 1 t {e}_{1t} , e 2 t {e}_{2t} are the forecast errors at time t t from the two models being compared.

The DM test evaluates the following hypotheses: (30) H 0 : E [ d t ] = 0 (Equal predictive accuracy) , {H}_{0}:{\mathbb{E}}\left[{d}_{t}]=0\hspace{1.0em}\hspace{0.1em}\text{(Equal predictive accuracy)}\hspace{0.1em}, (31) H 1 : E [ d t ] ≠ 0 (Unequal predictive accuracy) . {H}_{1}:{\mathbb{E}}\left[{d}_{t}]\ne 0\hspace{1.0em}\hspace{0.1em}\text{(Unequal predictive accuracy)}\hspace{0.1em}.

The test statistic is given by (32) DM = d ¯ σ ˆ d 2 ∕ T , {\rm{DM}}=\frac{\bar{d}}{\sqrt{{\hat{\sigma }}_{d}^{2}/T}}, where d ¯ = 1 T ∑ t = 1 T d t \bar{d}=\frac{1}{T}{\sum }_{t=1}^{T}{d}_{t} is the sample mean of the loss differential, T T is the forecast horizon, and σ ˆ d 2 {\hat{\sigma }}_{d}^{2} is a consistent estimate of the variance of d t {d}_{t} .

Under the null hypothesis, the DM statistic asymptotically follows a standard normal distribution: (33) DM ∼ N ( 0 , 1 ) . {\rm{DM}}\hspace{0.33em} \sim \hspace{0.33em}{\mathcal{N}}\left(0,1).

If the absolute value of the DM statistic exceeds the critical value from the standard normal distribution, the null hypothesis is rejected, indicating a statistically significant difference in predictive performance between the two models. In this study, the DM test is applied using MSE, MAE, and MAPE as loss functions to ensure comprehensive evaluation.

The experimental setup follows an 80:20 data split, where 80% of the time series is used for training and 20% for testing. We evaluate the forecasting performance of six models: ARIMA, SVR, RF, RNN, GRU, and LSTM. These models are applied to four key variables, with all learning carried out under a consistent framework.

To prepare the time series data for supervised learning for SVR, RF, and the RNNs, we apply a sliding window approach. 30-day sliding window is used to prepare the data. Given a univariate sequence y 1 , y 2 , … , y t {y}_{1},{y}_{2},\ldots ,{y}_{t} , each input-output pair is constructed as: (34) X t = ( y t − w , y t − w + 1 , … , y t − 1 ) , y t = y t , {X}_{t}=\left({y}_{t-w},{y}_{t-w+1},\ldots ,{y}_{t-1}),\hspace{1.0em}{y}_{t}={y}_{t}, where w w is the window size. This transformation allows the models to capture and learn temporal dependencies by treating past observations as features for predicting the next value.

For the DL models (RNN, GRU, LSTM), feature scaling is performed using Min-Max normalization to stabilize and accelerate the training process. The normalization formula is (35) X scaled = X − X min X max − X min . {X}_{\text{scaled}}=\frac{X-{X}_{\text{min}}}{{X}_{\text{max}}-{X}_{\text{min}}}.

After prediction, the inverse transformation is applied to recover the original scale: (36) X original = X scaled × ( X max − X min ) + X min . {X}_{\text{original}}={X}_{\text{scaled}}\times \left({X}_{\text{max}}-{X}_{\text{min}})+{X}_{\text{min}}.

Rescaling ensures that input features are within a uniform range, preventing dominance by features with larger numerical values. This enhances the convergence behavior of DL algorithms and contributes to better generalization performance.

Hyperparameter tuning was carried out for RF, ARIMA, and SVR using grid search or automated routines, while DL models (RNN, LSTM, GRU) employed fixed configurations.

For ARIMA, the model parameters were automatically selected using the auto_arima package in Python, resulting in an optimal configuration of ARIMA(4,1,0) for all datasets. The differencing order d = 1 d=1 was used to make the time series stationary.

For SVR, hyperparameter tuning was performed using grid search, which explored the following parameter space:

• C ∈ [ 1 , 10 , 100 ] C\in \left[1,10,100] ,

• γ ∈ [ 1 × 1 0 − 6 , 1 × 1 0 − 5 , 1 × 1 0 − 4 ] \gamma \in \left[1\times 1{0}^{-6},1\times 1{0}^{-5},1\times 1{0}^{-4}] ,

• Kernel: RBF, linear and polynomial.

C = 10 C=10

γ = 1 × 1 0 − 6 \gamma =1\times 1{0}^{-6}

For the RF model, hyperparameter tuning was performed using a grid search over a comprehensive parameter space to optimize predictive performance. The following hyperparameters were explored:

• Number of trees ( n estimators {n}_{\text{estimators}} ): { 50 , 100 , 200 } \left\{50,100,200\right\} ,

• Maximum tree depth ( max_depth \hspace{0.1em}\text{max\_depth}\hspace{0.1em} ): { None , 10 , 20 , 30 } \left\{\hspace{0.1em}\text{None}\hspace{0.1em},10,20,30\right\} ,

• Minimum samples required to split a node ( min_samples_split \hspace{0.1em}\text{min\_samples\_split}\hspace{0.1em} ): { 2 , 5 , 10 } \left\{2,5,10\right\} ,

• Minimum samples required at a leaf node ( min_samples_leaf \hspace{0.1em}\text{min\_samples\_leaf}\hspace{0.1em} ): { 1 , 2 , 4 } \left\{1,2,4\right\} ,

• Maximum number of features considered for splits ( max_features \hspace{0.1em}\text{max\_features}\hspace{0.1em} ): { auto , max_features } \left\{\hspace{0.1em}\text{auto},\sqrt{\text{max\_features}\hspace{0.1em}}\right\} ,

• Splitting criterion ( criterion \hspace{0.1em}\text{criterion}\hspace{0.1em} ): { mse , mae } \left\{\hspace{0.1em}\text{mse},\text{mae}\hspace{0.1em}\right\} .

For the DL models , we used a fixed architecture with three layers and ReLU activation. The optimizer for all DL models was Adam with a learning rate of 0.001. Dropout was applied with a rate of 0.2 to prevent overfitting, and early stopping was used to halt training if the validation loss did not improve. The models were trained for 100 epochs with a batch size of 16.

The parameter settings for the DL models are summarized in Table 1. These settings ensure consistency across all models, enabling a fair comparison of their performance.

Explanation of parameters for RNN, LSTM, and GRU models:

• Number of layers: All models have three layers of recurrent units (RNN, LSTM, or GRU). More layers allow the models to capture more complex patterns in sequential data.

• Activation function: The activation function used in all models is ReLU (rectified linear unit), which helps mitigate the vanishing gradient problem and accelerates convergence.

• Loss function: The MSE is used as the loss function for all models. It measures the average of the squares of the errors, making it sensitive to larger errors and ensuring that the model minimizes them.

• Optimizer: Adam (Adaptive Moment Estimation) is used for optimization in all models. It adjusts the learning rate during training to improve convergence speed and accuracy.

• Learning rate: Set to 0.001 for all models; this determines the step size during model training. A small learning rate helps with more stable training.

• Dropout rate: Dropout of 0.2 is used to prevent overfitting. It randomly disables 20% of neurons in the network during training, forcing the model to learn more robust representations.

• Epochs: All models are trained for 100 epochs to ensure enough iterations for convergence.

• Batch size: A batch size of 16 is used, meaning 16 samples are processed at once during each training iteration.

• Units per layer: The number of units in each layer is set to 100 for the first layer, 50 for the second, and 25 for the third layer, with the goal of progressively reducing the complexity of the model.

• Early stopping: This mechanism halts training when the validation loss stops improving, preventing overfitting and unnecessary training.

The Python packages used in this study include Numpy for numerical operations and array handling, Pandas for data manipulation and analysis, Seaborn for data visualization, Matplotlib for creating static, animated, and interactive visualizations, Scikit-learn for ML models including SVR, RF, and performance metrics, Statsmodels for statistical models and time series analysis, Pmdarima for automatic ARIMA model selection, and Keras for building and training DL models.