Predictive modeling of the COVID-19 data using a new version of the flexible Weibull model and machine learning techniques

1.   Introduction

The first infected case of the COVID-19 pandemic appeared in China during the last week of December 2019. Since that, the deadly pandemic has caused a dramatic loss of human life around the globe and presented a serious challenge to all sectors of life. For example, (i) to read about the COVID-19's effect on the economy, we refer to ^[1], (ii) education ^[2], (iii) healthcare ^[3], (iv) Business sectors ^[4], (v) sports sectors ^[5], and (vi) tourism sector ^[6].

Among the mentioned sectors affected by the COVID-19 epidemic, the healthcare sector is the most affected area. Based on the latest updates up to November 18, 2021, 10:36 GMT, around 255.9 million cases were registered worldwide, 5.1422 million people have died, and the total number of recovered cases has reached 231.26 million. The top five countries with the most registered COVID-19 cases are (i) United States of America (USA) with 48.29 million, (ii) India with 34.48 million, (iii) Brazil with 21.978 million, (iv) United Kingdom with 9.675 million, (v) Russia with 9.219 mission. The top five countries with the most death cases due to COVID-19 are (i) USA with 0.787 million, (ii) Brazil with 0.611 million, (iii) India 0.464 million, (iv) Mexico with 0.291 million, and (v) Russia with 0.260 million, https://www.worldometers.info/coronavirus.

The implementation of statistical models for dealing with lifetime data in healthcare sectors is very crucial and an important research area. In the available statistical distributions, the flexible weibull extension (FWE) model has attracted researchers ^[7]. The CDF (cumulative distribution function) $M(y; \boldsymbol{\kappa})$ of the FWE model is

where $\boldsymbol{\kappa} = (\tau_{1}, \tau_{2})$ with $\tau_{1} > 0$ and $\tau_{2} > 0.$

Corresponding to the CDF in Eq (1.1), the PDF (probability density function) $m(y; \boldsymbol{\kappa})$ of the FWE distribution is

El-Gohary et al. ^[8] developed the inverse version of the FWE distribution with CDF, given by

El-Gohary et al. ^[9] further generalized the WFE model by using the exponentiation strategy with CDF, given by

El-Damcese et al. ^[10] investigated another generalization of the FWE model by adopting the Kumaraswamy approach with CDF, given by

This paper introduces a further extension of the FWE model to improve the fitting power of the FWE model. To perform this activity, we combine $M(y; \boldsymbol{\kappa}) = 1-e^{-e^{\tau_{1} y -\frac{\tau_{2}}{y}}}$ with the Z-family approach ^[11]. The CDF $T(y; \beta, \boldsymbol{\kappa})$ of the Z-family is

with SF (survival function) $S(y; \beta, \boldsymbol{\kappa}) = 1-T(y; \beta, \boldsymbol{\kappa}),$ given by

where $\bar{M}(y; \boldsymbol{\kappa}) = 1-M(y; \boldsymbol{\kappa})$.

Corresponding to $T(y; \beta, \boldsymbol{\kappa})$, the PDF $t(y; \beta, \boldsymbol{\kappa}) = \frac{d}{dy}T(y; \beta, \boldsymbol{\kappa})$ and HF (hazard function) $h(y; \beta, \boldsymbol{\kappa}) = \frac{t(y; \beta, \boldsymbol{\kappa})}{1-T(y; \beta, \boldsymbol{\kappa})}$ are given by

and

respectively.

By incorporating $M(y; \boldsymbol{\kappa}) = 1-e^{-e^{\tau_{1} y -\frac{\tau_{2}}{y}}}$ in Eq (1.2), we obtain a new model, namely, the Z-flexible Weibull extension (Z-FWE) distribution. The Z-FWE is a more flexible version of the FWE model. This fact is illustrated by applying the Z-FWE distribution to a data set in the health sector.

2.   The Z-FWE distribution

If Y has the Z-FWE model with parameters $\beta > 0, \tau_{1} > 0,$ and $\tau_{2} > 0$, its CDF is

with PDF

respectively.

Different behaviors for the PDF of the Z-FWE distribution are shown visually in Figures 1 and 2.

Proposition 1. The PDF of the proposed model has a Unimodal (U) shape of bimodal (B) shape.

Proof. The first derivative of the PDF of the Z-FWE distribution model is determined as follows

where

The terms $Z(y; \beta, \boldsymbol{\kappa})$ and $R(y; \beta, \boldsymbol{\kappa})$ are given by

and

respectively.

We can see that $\frac{d}{dy} t(y; \beta, \boldsymbol{\kappa})$ and $L(y; \beta, \boldsymbol{\kappa})$ have the same signs. Also, $Z(y; \beta, \boldsymbol{\kappa})$ and $R(y; \beta, \boldsymbol{\kappa}) > 0$ $\forall {\; } y, \tau_{1} > 0, \tau_{2} > 0,$ and $\beta \geq 0.37$. When the sign of $L(y; \beta, \boldsymbol{\kappa})$ changes from $+$ to $-$, then there is a unique critical point which maximize the PDF of the Z-FWE distribution and give a unimodal shape. When the sign of $L(y; \beta, \boldsymbol{\kappa})$ changes from $+$ to $-$ twice, then there are two critical points which maximize the PDF of the Z-FWE distribution and give a bimodal shape.

The unimodal band biomodal behaviors of the PDF of the Z-FWE distribution are also proved numerically. The numerical results in Table 1 confirm the unimodal and bimodal behaviors of the PDF of the Z-FWE distribution.

Remark 1. The plots in Figure 1 are actually one plot, but we divided them into three plots to show that the PDF of the Z-FWE distribution can be bimodal. If we combine them in one plot, then, the bimodal shape of the Z-FWE distribution doesn't appear due to the range of the function. The actual range of the plot is from 0.00001 to 2. We have divided this range into three pieces, such as (i) 0.00001 to 0.00010 (ii) 0.00010 to 0.010, and (iii) 0.010 to 2.

3.   The heavy-tailed characteristic

This section offers the heavy-tailed behavior and regular variational results of the Z-FWE distribution. Probability distributions that are right-skewed and possess heavy-tailed behavior are very useful in providing the best description of the biomedical data sets. A probability model is called a heavy-tailed distribution, if it satisfies

An important property of the heavy-tailed probability distributions is called the regular variational property (RVP). A probability distribution is called regularly varying, if it satisfies

where $a$ represents an index of regular variation.

Here, we derive the RVP of the Z-FWE distribution. According to Karamata's theorem (Seneta ^[12]), in terms of SF $S\left(y; \beta, \boldsymbol{\kappa}\right)$, we have

Theorem: If $S\left(y; \boldsymbol{\kappa}\right) = 1 - M\left(y; \boldsymbol{\kappa}\right)$ is the SF of the regular varying model, then $S\left(y; \beta, \boldsymbol{\kappa}\right)$ is also a regular varying distribution.

Proof: Suppose $\lim_{y\longrightarrow \infty} \frac{S\left(p y; \beta, \boldsymbol{\kappa}\right)}{S\left(y; \beta, \boldsymbol{\kappa}\right)} = f\left(p\right)$ is finite but non-zero $\forall \; p > 0$. Then, using Eq (1.3), we have

Applying $\lim_{y\longrightarrow \infty}$ on both sides of Eq (3.1), we get

Using Eq (1.1) in Eq (3.2), we get

where $f\left(p\right)$ has the form $p^{a}$. So, the expression in Eq (3.3) is finite and non-zero $\forall \; p > 0.$ thus, $S\left(py; \beta, \boldsymbol{\kappa}\right)$ is a regular varying distribution.

4.   Characterizations of the Z-FWE model

This section offers the characterizations of the Z-FWE distribution via implementing three different approaches. These characterizations are obtained using (i) a simple relationship between two TMs (truncated moments), (ii) the HF, and (iii) the CE (conditional expectation) of a function of the RV (random variable).

4.1.   Characterization using the TMs

The result of the first characterization is obtained due to the theorem in Glänzel ^[13], see the below Theorem 1. It is important to note that the result of the first characterization also holds true when the interval expressed by $W$ is not closed.

Furthermore, it could also be implemented when the CDF $F$ is not in closed form. As Glänzel ^[14] showed that this characterization is stable even under weak convergence.

Theorem 1. Consider a given PS (probability space) $\left(\Omega, \mathcal{F}, \mathbf{P}\right)$ and let $W = \left[\theta_{1}, \theta_{2}\right]$ represents an interval for some $\theta_{1} < \theta_{2} \left(\theta_{1} = -\infty, \theta_{2} = \infty\text{ might as well be allowed}\right).$ Now, let $Y:\Omega\rightarrow W$ be a continuous RV with the CDF $T$ and let $g(.)$ and $h(.)$ be two real functions defined on the interval $W$ such that

is defined with some real function $\vartheta$. Let suppose that $g$, $h\in K^{1}\left(W\right)$, $\vartheta \in K^{2}\left(W\right)$ and $T$ is twice CD (continuously differentiable) and strictly MF (monotone function) on the set $W$. Also, let assume that the equation $\vartheta h = g$ has no real solution in the interior of $W$. Then $T$ is uniquely determined by the functions $g$, $h$ and $\vartheta$, particularly

where $s$ is a solution of the DE (differential equation) given by $s^{\prime} = \frac{\vartheta^{\prime} \ h}{\vartheta \ h \ -\text{}g}$ and the quantity $K$ represents a normalization constant, such that $\int_{W}dT = 1$.

Proposition 2. Let $Y:\Omega\rightarrow \left(0, \infty\right)$ be a continuous RV and let $h\left(y\right) = \left[1+\log\left(\beta\right) e^{-e^{\tau_{1}y-\frac{\tau_{2}}{y}}}\right]^{-1}$ and $g\left(y\right) = h\left(y\right) e^{-e^{\tau_{1} y-\frac{\tau_{2}}{y}}}$ for $y > 0.$ The RV $Y$ has PDF in Eq (2.2) if and only if the function $\vartheta$ defined in Theorem 1 has the form

Proof. Let $Y$ be a RV with PDF in Eq (2.2), then

and

and finally

Conversely, if the function $\kappa$ is given as mentioned above, then

and hence

Now, in view of the result in Theorem 1, $Y$ has PDF in Eq (2.2).

Corollary 1. Consider a RV $Y:\Omega\rightarrow \left(0, \infty\right)$ and let the function $h\left(y\right)$ be as defined in Proposition 2. The PDF of $Y$ is Eq (2.2), if and only if the functions $g$ and $\vartheta$ defined in Theorem 1 satisfying the below DE

Corollary 2. The general solution of the DE in Corollary 1 is

where the term $C$ is a constant quantity. It should be noted that a set of functions obeying the DE provided above is given in Proposition 2 with $C = 0.$

4.2.   Characterization using the HF

It is obvious that the HF, $h_{T}$, of a twice differentiable DF, $T$ with PDF $t$, obeys the first order DE

In terms of the HF, this is the only characterization available for many univariate CMs (continuous models). In terms of the HF, the following characterization establishes a characterization of Z-FWE distribution, which is not of the trivial form proivded above.

Proposition 3. Consider $Y:\Omega\rightarrow \left(0, \infty\right)$ be a continuous RV. Then $Y$ has the PDF in Eq (2.2), if and only if its HF $h_{T}\left(y\right)$ obyes the DE

with the initial condition given by $\lim_{y\rightarrow0}h_{T}\left(y\right) = 0$.

Proof. Is straightforward and hence omitted.

4.3.   Characterization using the CE

The CE is a well-known approach to characterize statistical distributions; see Hamedani ^[15]. Here, we use this approach to characterize the Z-FWE distribution.

Proposition 4. Let $Y:\Omega\rightarrow$ $\left(\lambda_{1}, \lambda_{2}\right)$ be a continuous RV with CDF $T$. Let $\varphi\left(y\right)$ be a differentiable function on $\left(\lambda_{1}, \lambda_{2}\right)$ with $\lim_{y\rightarrow0^{+}}\varphi\left(y\right) = 1$. Then for $\eta\neq1,$

if and only if

Remark 2. For $\left(\lambda_{1}, \lambda_{1} \right) = \left(0, \infty\right)$ $,$ $\varphi\left(y\right) = \frac{e^{-\frac{1}{2}e^{\tau_{1}y-\frac{\tau_{2} }{y}}}}{\beta^{\frac{1}{2}\left(1-e^{-e^{\tau_{1}y-\frac{\tau_{2}}{y}} }\right) }}$ and $\eta = \frac{2}{3}$, Proposition 2 provides a characterization of Z-FWE distribution.

5.   Estimation and simulation

This section is composed of two subsections. The very first subsection offers the derivation of the estimators $\left(\hat{\tau_{1}}, \hat{\tau_{2}}, \hat{\beta}\right)$ of the parameters $\left(\tau_{1}, \tau_{2}, \beta\right)$ of the Z-FWE distribution. Whereas, the second subsection, offers the evaluation of $\hat{\tau_{1}}, \hat{\tau_{2}}$ and $\hat{\beta}$ through a simulation study.

5.1.   Estimation

Here, we obtain the estimators $\hat{\tau_{1}}, \hat{\tau_{2}}$ and $\hat{\beta}$ of the parameters of the Z-FWE model $\tau_{1}, \tau_{2}$ and $\beta$, respectively. The estimators $\hat{\tau_{1}}, \hat{\tau_{2}}$ and $\hat{\beta}$ are obtained by using five different approaches. These methods are the (i) maximum likelihood estimation, (ii) ordinary least-square estimation, (iii) weighted least-square estimation, (iv) maximum product of spacing estimation, and (v) Anderson-Darling estimation.

5.1.1.   Maximum likelihood estimation

Consider a RS (random sample) say $y_{1}, y_{2}, ..., y_{n}$ taken from $t(y; \beta, \boldsymbol{\kappa})$. Then in link to $t(y; \beta, \boldsymbol{\kappa}),$ the LF (likelihood function) $\lambda \left(y; \beta, \boldsymbol{\kappa}\right)$ is

Using Eq (2.2) in Eq (5.1), we have

Corresponding to Eq (5.2), the log LF $\delta \left(y; \beta, \boldsymbol{\kappa}\right)$ is given by

With respect to $\tau_{1}, \tau_{2}$ and $\beta$, the partial derivatives of $\delta \left(y; \beta, \boldsymbol{\kappa}\right)$ are given by

and

respectively.

In solving $\frac{\partial}{\partial \tau_{1}} \delta \left(y; \beta, \boldsymbol{\kappa}\right) = 0$, $\frac{\partial}{\partial \tau_{2}} \delta \left(y; \beta, \boldsymbol{\kappa}\right) = 0$, and $\frac{\partial}{\partial \beta} \delta \left(y; \beta, \boldsymbol{\kappa}\right) = 0,$ we get the maximum likelihood estimators (MLEs) $\hat{\tau_{1}}, \hat{\tau_{2}}$ and $\hat{\beta},$ respectively.

5.1.2.   Ordinary least-square estimation

Let $y_{(1)}, y_{(2)}, \cdots, y_{(n)}$ be the order statistics of a sample of size $n$ from $T\left(y; \beta, \boldsymbol{\kappa}\right)$ in Eq (1.2). The ordinary least-squares estimators (OLSEs) $\widehat{\beta}_{OLSE}$ and $\widehat{\boldsymbol{\kappa}}_{OLSE}$ can be obtained by minimizing

with respect to $\beta$ and $\boldsymbol{\kappa}$. Or equivalently, the OLSEs follow by solving the non-linear equations

where

and

Note that the solution of $\Delta_{m}$ for $m = 1, 2$ can be obtained numerically.

5.1.3.   Weighted least-square estimation

The weighted least-squares estimators (WLSEs) $\widehat{\beta}_{WLSE}$ and $\widehat{\boldsymbol{\kappa}}_{WLSE}$ can be obtained by minimizing

Moreover, the WLSEs can also be obtained by solving the non-linear equations

5.1.4.   Maximum product of spacing estimation

Let $D_{k}\left(\beta, \boldsymbol{\kappa}\right) = T\left(y_{(k)}|\beta, \boldsymbol{\kappa}\right) - T\left(y_{(k-1)}|\beta, \boldsymbol{\kappa}\right)$, for $k = 1, 2, \ldots, n+1,$ be the uniform spacing of a random sample from the Z-FWE distribution, where $T\left(y_{(0)}|\beta, \boldsymbol{\kappa}\right) = 0$, $T\left(y_{(n+1)}|\beta, \boldsymbol{\kappa} \right) = 1,$ and $\sum_{k = 1}^{n+1} D_{k} \left(\beta, \boldsymbol{\kappa}\right) = 1$. Then, the maximum product of spacing estimators (MPSEs) $\widehat{\beta}_{MPSE}$ and $\widehat{\boldsymbol{\kappa}}_{MPSE}$ can be obtained by maximizing the geometric mean of the spacing

with respect to $\beta$ and $\boldsymbol{\kappa}$, or, equivalently, by maximizing the logarithm of the geometric mean of sample spacing

The MPSEs of the Z-FWE parameters can be obtained by solving the nonlinear equations defined by

5.1.5.   The Anderson-Darling estimation

The Anderson-Darling estimators (ADEs) of the parameters of the Z-FWE distribution are obtained by minimizing

with respect to $\beta$ and $\boldsymbol{\kappa}$. These ADEs can also be obtained by solving the non-linear equations

5.2.   Simulation study

Furthermore, we investigate the performances of $\hat{\tau_{1}}, \hat{\tau_{2}},$ and $\hat{\beta}$ via a simulation study. To carry out the evaluation of $\hat{\tau_{1}}, \hat{\tau_{2}},$ and $\hat{\beta}$, random samples, say, $n = 20, 50,150,300$ are generated from $t(y; \beta, \boldsymbol{\kappa})$ using the inverse DF approach.

The evaluation of $\hat{\tau_{1}}, \hat{\tau_{2}},$ and $\hat{\beta}$ has been done for three sets of parameters values such as (a) $\hat{\tau_{1}} = 0.6, \hat{\tau_{2}} = 1.3, \hat{\beta} = 0.9$ (b) $\hat{\tau_{1}} = 0.12, \hat{\tau_{2}} = 0.8, \hat{\beta} = 1.2$, and (c) $\hat{\tau_{1}} = 1.3, \hat{\tau_{2}} = 1.5, \hat{\beta} = 0.7$.

Furthermore, certain statistical methodologies such as (i) MSEs and (ii) absolute biases were selected to check the performances of $\hat{\tau_{1}}, \hat{\tau_{2}},$ and $\hat{\beta}$. The values of the selected measures were respectively obtained using the expressions

and

where $\hat{\boldsymbol{\epsilon}} = \left(\beta, \boldsymbol{\kappa}\right).$

Corresponding to (a) $\hat{\tau_{1}} = 0.6, \hat{\tau_{2}} = 1.3, \hat{\beta} = 0.9$ (b) $\hat{\tau_{1}} = 0.12, \hat{\tau_{2}} = 0.8, \hat{\beta} = 1.2$, and (c) $\hat{\tau_{1}} = 1.3, \hat{\tau_{2}} = 1.5, \hat{\beta} = 0.7$. the simulation results of the Z-FWE model are provided in Table 2.

6.   Data analysis

This section is devoted to data analysis to illustrate the crucial role of the Z-FWE in data modeling. To carry out the illustration of the Z-FWE model, a data set from the medical sector is considered. The data consists of one hundred and six (106) observations and represents the mortality rate of patients during the COVID-19 pandemic in Mexico. This was recorded during the period between March 4, 2020, to July 20, 2020; see Almongy et al. ^[16]. For the simplicity of analysis, each observation is divided by five and is given by: 1.7652, 1.2210, 1.8782, 2.9924, 2.0766, 1.4534, 2.6440, 3.2996, 2.3330, 1.2030, 2.1710, 1.2244, 1.3312, 0.6880, 1.1708, 2.1370, 2.0070, 1.0484, 0.8688, 1.0286, 1.5260, 2.9208, 1.5806, 1.2740, 0.7074, 1.2654, 0.9460, 0.6430, 1.8568, 2.5756, 1.7626, 2.0086, 1.4520, 1.1970, 1.2824, 0.6790, 0.8848, 1.9870, 1.5680, 1.9100, 0.6998, 0.7502, 1.3936, 0.6572, 2.0316, 1.6216, 1.3394, 1.4302, 1.3120, 0.4154, 0.7556, 0.5976, 0.6672, 1.3628, 1.6650, 1.5708, 1.7102, 0.6456, 1.4972, 1.3250, 1.2280, 0.9818, 0.9322, 1.0784, 2.4084, 1.7392, 0.3630, 0.6654, 1.0812, 1.2364, 0.2082, 0.3600, 0.9898, 0.8178, 0.6718, 0.4140, 0.6596, 1.0634, 1.0884, 0.9114, 0.8584, 0.5000, 1.3070, 0.9296, 0.9394, 1.0918, 0.8240, 0.7844, 0.6438, 0.2804, 0.4876, 0.6514, 0.7264, 0.6466, 0.6054, 0.4704, 0.2410, 0.6436, 0.5852, 0.5202, 0.4130, 0.6058, 0.4116, 0.4652, 0.5012, 0.3846.

The summary measures of this data are: smallest observation (minimum) = 0.2082, 1st Quartile = 0.6578, Median = 1.0559, Mean = 1.1645, 3rd Quartile = 1.5188, the largest value (maximum) = 3.2996, variance = 0.4225, and standard deviation = 0.6500. Corresponding to the mortality rate data, the total time test (TTT) plot is obtained in Figure 3. In addition to the TTT plot, the histogram of the mortality rate data is also provided within the same figure.

Using the given mortality rate data, we show the closer-fitting capability of the Z-WFE distribution. To perform the numerical illustration of the Z-FWE distribution, its comparison is done with some well-known competing models. These models include the (i) baseline FWE model, (ii) the exponentiated version of the FWE distribution, namely, exponentiated FWE (E-FWE), (iii) Weibull model, (iv) a generalization of the Weibull, namely, exponentiated Weibull (E-Weibull) distribution, (v) another famous extension of the Weibull is called, the Kumaraswamy Weibull (K-Weibull) distribution, (v) another generalized form of the Weibull model, namely, a modified new flexible Weibull extension (MNFWE) distribution, and (vi) a new modified flexible Weibull extension (NMFWE) distribution.

The above-mentioned competing models are very popular for modeling and describing real-phenomena of nature. The SFs of the selected competing models are

● FWE

where $\tau_{1} > 0$ and $\tau_{2} > 0.$

● E-FWE

where $\tau_{1} > 0$, $\tau_{2} > 0$, and $\theta_{1} > 0$.

● Weibull

where $\tau_{1} > 0$ and $\tau_{2} > 0$.

● E-Weibull

where $\tau_{1} > 0$, $\tau_{2} > 0$, and $\theta_{1} > 0$.

● K-Weibull

where $\tau_{1} > 0$, $\tau_{2} > 0$, $\theta_{1} > 0,$ and $\theta_{2} > 0$.

● MNFWE

where $\tau_{1} > 0$, $\tau_{2} > 0$, $\theta > 0,$ and $\lambda > 0$.

● NMFWE

After selecting the competing models, the very next step is to choose statistical tools to evaluate the performance of the Z-FWE and other models. For the evaluation of these distributions, certain statistical tools and tests were chosen. These tools are given by

● Akaike information criterion

The Akaike information criterion (AIC) is a mathematical approach for testing and evaluating how close/well a probability model fits the given data. In the literature, the AIC is calculated to compare different selected competing models and determine which one better fits the data set under consideration. In this paper, we consider the AIC to compare the fitting results of the Z-FWE distribution with other considered competing distributions. It is calculated as

● Bayesian information criterion

Another criterion that we considered as a comparative tool is called the Bayesian information criterion (BIC). It is also called the Schwarz information criterion. For a given data set, the BIC also determine the best competing model among a finite set of models. It is calculated as

● Corrected Akaike information criterion

The Corrected Akaike information criterion (CAIC) is another useful statistical tool for checking the quality of the fitting of the possible competing distributions. For a given data set, a probability model with the lowest value of CAIC is considered the best competing distribution. The CAIC is calculated as

● Hannan–Quinn information criterion

The Hannan–Quinn information criterion (HQIC) is another statistical criterion for model selection among a set of possible competing statistical distributions. It can be used as an alternative to AIC and BIC for checking the fitting power of the competing distributions. It is given as

● Anderson Darling test

The Anderson Darling (AD) test is a well-known statistical procedure to check whether a sample of the underlined data is drawn from a given statistical distribution. The AD test is most often implemented in contexts where a family of statistical distributions is being tested. The AD test statistic is calculated as

● Cramér-von Mises test

The Cramér-von Mises (CM) test is another useful statistical approach used for judging the goodness of fit of the given probability distributions. The CM test statistic is calculated as

● Kolmogorov-Smirnov test

The Kolmogorov–Smirnov (KS) test is a prominent statistical approach that is used to quantify the distance between the empirical CDF of the sample and the CDF of the reference probability distribution. A lower distance between the empirical CDF and the CDF of the probability indicates the best fitting of the corresponding probability distribution. The KS test statistic is calculated as

Furthermore, another statistical quantity called p-value is also considered a tool for the illustration of the fitted model. A higher p-value associated with any model indicates the best competitor among the competing models.

After performing the analysis, the numerical values of $\hat{\tau_{1}}, \hat{\tau_{2}}, \hat{\beta}, \hat{\theta_{1}}, \hat{\theta_{2}}, \hat{\theta}$, and $\hat{\lambda}$ are presented in Table 3. Whereas, the values of the comparative tools are presented in Tables 4 and 5.

From the numerical illustrations in Tables 4 and 5, we observe that the Z-FWE distribution is the best competitor for modeling the mortality rate data. For the Z-FWE distributions, the numerical values of the selected measures are AIC = 186.07860, CAIC = 186.31390, BIC = 194.06900, HQIC = 189.31720, CM = 0.03228, AD = 0.20038, KS = 0.04683, with p-value = 0.97420. Whereas, the second best model is the NMFWE distribution. For the NMFWE distribution, the values of the analytical measures are AIC = 186.12600, CAIC = 186.36130, BIC = 194.11630, HQIC = 189.36450, CM = 0.03276, AD = 0.20485, KS = 0.05085, and p-value = 0.94680.

From the above discussion, as well as the results in Tables 4 and 5, it is now obvious that the Z-FWE is the best choice for modeling the mortality rate data. Visual support for the numerical illustrations is presented in Figure 4. For the visual illustration of the Z-FWE distribution, the plots of the CDF, SF, PP, and QQ functions were obtained. These plots visually confirm the best fitting of the Z-FWE distribution.

7.   Forecasting models

In the previous section, we compared several distributions with the Z-FWE model using the mortality rate of patients due to the COVID-19 pandemic in Mexico. In this section, the study compares time series and machine learning techniques through forecasting on the same data set. Autoregressive Integrated Moving Average is a time series model, whereas Artificial Neural Network and group methods of data handling. For comparison purposes, we divide the data set into two parts: a training set and a testing set. The training set consists of 80 percent data and the testing set consists of 20 percent data, followed by ^[17]. The models are estimated on training data and then their performance is assessed on testing data. Below, the following methods are elaborated.

7.1.   Autoregressive integrated moving average Model

Numerous linear models have been devised in the literature on time series forecasting. Among them, the ARIMA model is the most well-known and frequently used over the previous few decades that has enjoyed fruitful applications in engineering, economics, and finance ^[18]. Basically, the ARIMA model has included the components, the autoregressive (AR) and moving average (MA). If the underlying time series is non-stationary, then it is transformed by using differencing approach ^[19].

ARIMA model specification is needed in several steps. In the initial step, we apply the autocorrelation function (ACF) and partial autocorrelation (PACF) to determine the order of AR and MA terms included in the model. In the second step, the estimated model is passed through diagnostic tests. If the estimated model is not suitable, a new tentative model should be diagnosed and will be followed by the same steps i.e., parameters estimation and model verification. This whole process is repeated until an adequate model is finally selected. In the last step, the final model can be utilized for forecasting aims ^[20]. Mathematically, the ARIMA model can be expressed as:

where $\beta_{t}$ and $M_{t}$ represent the error term and observed value at time t, respectively. $\vartheta (U) = 1-\sum_{u = 1}^{p}\vartheta_{u} U^{u}$ and $\varphi (U) = 1-\sum_{v = 1}^{q}\varphi_{v} U^{v}$ are polynomials in U of degree p and q, $\vartheta_{u} (1, 2, ..., p)$ and $\varphi_{v} (1, 2, ..., q)$ are the unknown parameters of the model. $\nabla = (1-U)$, where U is the backward shift operator, d denotes the length of difference, p and q indicate the length of AR and MA terms, respectively.

7.2.   Machine learning techniques

Qurban et al. ^[21] declared that nonlinear problems are commonly overcome by machine learning (ML) techniques. The most widely used ML techniques are Artificial neural networks (ANNs) and the group method of data handling (GMDH).

7.2.1.   Artificial neural network

The ANN models can model a wide range of non-linear problems. In comparison with other nonlinear models, ANN can approximate a huge group of nonlinear functions with good accuracy. This accuracy is achieved due to the parallel processing of the information from the data. The construction of the ANN model does not require any prior assumption regarding model form. Rather, the model is substantially ascertained by the characteristics of the data ^[22]. There are many types of ANN, but multilayer perceptron (MLP), is one of the most popular forms of ANN for forecasting in time series setup. The nexus between input and output layers has the following mathematical representation

where $\varphi_{uv}(u = 1, 2, ..., p, {\; \; } v = 1, 2, ..., q)$ and $\varphi_{v} (v = 0, 1, 2, ..., q)$ represent connection weights, p denotes input layers and q denotes the number of hidden layers. There is no hard and fast rule to select the input layers and hidden layers. The sole possible way is the error and trial approach to select the optimum number of p and q; see Khashei and Hajirahimi ^[23].

7.2.2.   Group method of data handling

The GMDH-type neural network can be illustrated as a collection of neurons where numerous pairs are related throughout every layer utilizing quadratic polynomial, which in turn induces new neurons all over the next layers. Such type of information can be utilized while connecting the input layers with the output layer. If the primary structure is tentative then it leads to a slight improvement in parameters estimation of a model. Primarily, the GMDH tries to achieve the hierarchic solution, by seeking a large set of simple models, keeping the best, and developing them iteratively to get a composition of functions like a statistical model ^[24]. The polynomial nodes or building blocks have often the quadratic form

M is the output of the model, $u_{i} (i = 1, 2)$ shows the input variables and the weight vector $\varphi.$

7.3.   Error metrics and Diebold Mariano test

The predictive power of the statistical models is assessed by using metrics computed from a holdout set (test data). From a statistical point of view, the forecast error is a more plausible criterion to judge the forecasting ability and choose the best method. The commonly used criteria such as root-mean-square error (RMSE), and mean absolute error (MAE) (Diebold and Mariano ^[25]) test are employed in this study.

7.3.1.   Root mean square error

The RMSE measures the magnitude of error in forecast comparison, and it is a quadratic scaled measure. The RMSE gives relatively more weight to large errors. A low value of RMSE for the model is preferred for a good forecast. It is given as

7.3.2.   Mean absolute error

The mean absolute error (MAE) is another widely implemented measure of forecasting accuracy. It is defined as

In the RMSE and MAPE formulas, $Y_{t}$ and $\hat{Y_{t}}$ indicate the actual and forecast values respectively, and H shows the forecast horizon.

7.3.3.   Diebold and Mariano test

When using different models for the prediction of a single variable, then it is better to check which model outperforms in terms of prediction. For this purpose, Diebold and Mariano (DM) test can be used. The DM statistic is a squared error difference mean of two competing predictive models, whose covariance matrix is obtained by accounting for the autocorrelation produced in the multistep forecasts. The DM test statistic tests the null hypotheses that two models have the same forecasting ability, i.e., to test forecast errors of two sets, say $e1i$ and $e2i$ having equal mean. A loss function is defined to evaluate the forecast error, mainly squared forecast error and absolute forecast error are used. The corresponding test statistics are used for hypothesis testing, that the two models have equal predictive performance, which is given by

where $\bar{D} = \frac{1}{p} \sum_{k = 1}^{p} D_{k}$, the $\bar{D}$ has an asymptotic variance is expressed as

where $\omega_{0}$ is the variance of $D_{k}$ and $\omega_{k}$ is the $k^{th}$ autocovariance of $D_{k}$ and p is the number of forecast steps, which is estimated by

7.4.   Empirical results

In this section, we elaborate on the findings of our prediction experiments along with a graphical representation. Figure 5 indicates the non-stationarity pattern in the series. In other words, the statistical properties such as mean, variance, and covariance of the original series are not constant over time. It represents a trend in the data, but we require a plain series having no trend. To achieve stationary series, we performed a differenced transformation. The ACF and PACF of the original and differenced series are given in Figure 6. It can be observed that the ACF associated with the original series is gradually dropping, but after transformation, the ACF declines at a swift pace. This shows that the COVID-19 series is the first difference stationary.

The output of the Box-Ljung test and Jarqua-Bera test are reported in Table 6. The corresponding p-values are greater than 5 percent, thus we cannot reject the null hypothesis, under which it is assumed that the residuals of the estimated model are random and normally distributed. In other words, it is stated that the residuals of the estimated model are normally distributed and uncorrelated. Thus, the model can be used for forecasting.

Another way to know the randomness and normality of the residuals of the fitted model is by considering the graphs of ACF, PACF, and QQ-plot of the residuals; see Figure 7. They reveal that the residuals are random and normally distributed.

The RMSE and MAE using the COVID-19 dataset are given in Table 7. We can observe that RMSE and MAE associated with machine learning (ML) algorithms such as ANN and GMDH are considerably lower than the time series model. Thus, it can be inferred that machine learning algorithms are superior to traditional time series models in terms of forecasting. Moreover, considering ML algorithms, the GMDH has produced a more accurate forecast than ANN. The forecast comparison is also portrayed in Figure 8. The plots in Figure 8 reveal that the ML algorithms, particularly GMDH, have captured the COVID-19 trend very well. Basically, these error metrics are loss functions, and in the field of statistics and econometrics, we are often interested in statistical differences. Therefore, apply the DM test for this purpose.

The values in Table 8 indicate the p-values for the DM test. Under the null hypothesis, it is assumed that models in row and column have the same forecasting accuracy. Under the rival hypothesis, the model in the column is more accurate than the row model. We found another piece of evidence in the form of a statistical test. From Table 8, it is quite clear that p-values are sufficiently small, which clearly demonstrates that the ML algorithms have high predictive power than ARIMA model.

8.   Concluding remarks

The COVID-19 pandemic has dramatically affected the economy, education, and health sectors, etc. Among them, the health sector is the most affected one. To have the best description and knowledge of the COVID-19 pandemic, numerous statistical studies have appeared. This paper offers a new statistical model for analyzing the mortality rate of the COVID-19 pandemic in Mexico. Some characterizations along with estimators of the new model were obtained. The new model was named the Z-FWE distribution and applied to the COVID-19 data in comparison with other statistical models. Based on the findings of this research, it is shown that the Z-FWE model was the best competitor for dealing with mortality rate data. Later, we compared the predictive power of ML techniques, i.e., ANN and GMDH, with the ARIMA. For comparison, we used RMSE, MAE, and DM tests. Based on all these criteria, it is inferred that ML techniques are superior in terms of forecasting than the ARIMA model. Furthermore, GMDH outperformed the ANN while getting small forecast errors. The study recommends that policymakers can utilize ML techniques, particularly GMDH, for forecasting the mortality rate of COVID-19 patients.

Acknowledgments

This research work was funded by Institutional Fund Projects under grant no. (IFPIP: 549-150-1443). The authors gratefully acknowledge technical and financial support provided by the Ministry of Education and King Abdulaziz University, DSR, Jeddah, Saudi Arabia.

Conflict of interest

All authors declare no conflicts of interest in this paper.