Modeling of daily confirmed Saudi COVID-19 cases using inverted exponential regression

1.   Introduction

Since the beginning of 2020, the world has been facing a coronavirus pandemic (COVID-19) at a rapid and alarming rate, representing a significant challenge for humanity and a serious threat to life. Although many countries have taken multiple and sometimes harsh measures to limit the spread of the virus and reduce its spread, the eyes of the public are now turning to scientists, doctors, and researchers of all scientific disciplines in the hope of finding a quick and successful treatment for this virus. Many authors have investigated the association of environmental and meteorological factors on the spread of COVID-19, see for example, Tello-Leal, et al. ^[1], Kodera, et al. ^[2], Meo, et al. ^[3], Casado-Aranda, et al. ^[4], Dogan, et al. ^[5], Fu, et al. ^[6] and Yuan, et al. ^[7].

McCullagh and Nelder ^[8] published a book on GLMs that led to their widespread use and appreciations. They extended the scoring method to maximum-likelihood estimation (MLE) in exponential families. Nelder and Pregibon ^[9] described methods of jointly estimating parameters of both link and variance functions. The iteratively reweighted least squares (IRLS) algorithm is amenable to statistics and measures that are common to all GLMs. Nelder and Wedderburn ^[10] used the Newton-Raphson process for regression coefficient estimates. Yuan and Bentler ^[11] reported that the convergence properties of the Fisher scoring algorithm are affected by many factors. One among the observed variables is multicollinearity. If the sample or model implied covariance matrix is close to singular, the Fisher scoring algorithm may have difficulty reaching a set of converged solutions. Nelder and Wedderburn ^[10] reported that the Newton-Raphson process with expected second derivatives is equivalent to Fisher's scoring technique. Additionally, de Jong and Heller ^[12] reported that the Newton-Raphson iteration equation leads to a sequence that often rapidly converges. This includes the D statistic along with specific residuals and influence measures. Liao ^[13] introduced a systematic way of interpreting commonly used probability models: logit, probit, and other GLMs. For recent works of using these models to the field of epidemiology and to the healthcare, see for example, Richardson and Hartman ^[14], Song, et al. ^[15], Mohamadou, et al. ^[16] and Trunfio, et al. ^[17].

The inverted exponential (IE) distribution that was introduced by Keller and Kamath ^[18]. The role of IE distributions is indispensable in many applications of reliability theory for its memoryless property and its constant failure rate. Dey ^[20] considered IE distribution as life distribution (see Abdel-Aty et al. ^[19], and Dey ^[20]). Singh et al. ^[21] obtained Bayes estimates for parameters of IE distribution by using informative and noninformative priors. They also compared the classical method with the Bayesian through the simulation study.

The Bayesian approach of the statistical modeling provides an alternative to standard GLMs. Posterior-mode estimation is an alternative to full posterior analysis or posterior mean estimation, which avoids numerical integrations or simulation methods. It was proposed by many authors (for more details, see Fahrmeir and Tutz ^[22] and Cepeda and Gamerman ^[23]). Dey et al. ^[24] described how to conceptualize, perform, and critique traditional GLMs from a Bayesian perspective, and how to use modern computational methods to summarize inferences using simulations. Olsson ^[25] gave an overview of GLMs and presented practical examples. The exponential family of distributions are discussed with maximum-likelihood estimation and ways of assessing the fit of the model. For the Bayesian estimation in this context, a useful asymmetric loss function known as the LINEX loss function was introduced by Varian (1975) and has been widely used by several authors. A highly used one is the zero-one loss function (for more details, see Sano et al. ^[26]). Robbins ^[27] has provided a more robust estimate and for estimating the parameters of prior distribution (hyperparameters), studies on the empirical Bayes (EB) method, Wei ^[28] has proposed the EB test of the regression coefficient, and working out the EB test decision rule by using kernel estimation of multivariate density function and its first-order partial derivatives, Singh ^[29] has proposed an EB approach in a multiple linear-regression model, Houston and Woodruff ^[30] have derived EB estimates using an m-group regression model to regress within-group estimates toward common values. More studies by many authors discussed EB test problems for parameters in a class of linear-regression model and other topics, e.g., Wind ^[31], Huang ^[32], Karunamuni ^[33], Yuan ^[34], Chen ^[35], Efron ^[36], Shao ^[37], Kim and Nembhard ^[38], and Jampachaisri et al. ^[39].

In order to reduce the influence of outliers on the estimate, some robust measures were proposed in the literature. The common robust estimation method can be divided into several categories: M, MM, median, L1, Msplit, R, S, least-trimmed squares, and sign-constraint robust least squares estimation. Among these, Huber's M estimation has become one of the main robust estimation methods by virtue of its simple calculation and convenience to implement (see Li et al. ^[40]). The key aspect is the involvement of a loss function that is applied to data errors that was selected to less rapidly increase than the square loss function that is used in least-squares or maximum-likelihood procedures. There exist several well-known families of loss functions, such as Huber, Hampel, and Tukey's biweight (or bisquare) that can be used for the computation of M estimators (see Sinova and Aelst ^[41]).

A major contribution of this paper is to propose similarity to GLMs, except that the distribution of the response is not a member of the exponential family using the Bayesian approach. We interest in IE distribution of which the flexible distribution can describe different lifetimes from medicine, reliability, ecology, biological studies, and other areas. We propose the Bayesian and non Bayesian inverted exponential regression models to model and analyze Covid-19-related data with the aim of explaining the relationships between Covid-19 cases and environmental-related variables. The paper is organized as follows: In Section 2, we present an overview of GLMs and propose the IE-Reg model under a log link function. In Section 3, we perform IE-BReg and IE-EBReg under a gamma prior, log link, and two loss functions. We propose Huber's and Tukey's bisquare (biweight) function to improve Bayesian models. We also adopted the iteratively reweighted least squares (IRLS) algorithm to estimate the Bayesian regression coefficients. In Section 4, we apply IE-Reg, IE-BReg, and IE-EBReg models, including an estimation, and use criteria such as AIC, BIC, MSE, and D to the Saudi COVID-19 dataset collected from March 23 to June 21, 2020. Finally, Section 5 draws a succinct conclusion to the findings.

2.   Classical approach

Nelder and Wedderburn ^[10] introduced the class of GLMs, defined according to the assumption that $y_{1}, y_{2}, ...y_{n}$ are observations of the response variable, with density function $y_{i}$ as follows:

where $\psi(\cdot)$, $c(\cdot)$ are known functions, with $\theta_{i}$ being the canonical parameter. Link function $g(.)$, related to the regression coefficients, is given by

where $g(\mu_{i}) = \theta_{i}$, $\beta = (\beta_1, ..., \beta_p)^{\prime}$ is a vector of $p$ unknown regression parameters, $x^{\prime}_{i} = (x_{i1}, x_{i2}, ..., x_{ip})$ is a vector of explanatory variables, and $\eta_{i}$ is a linear predictor of vectors $x^{\prime}_{i}$ and $\beta$. Here, $g(.): (0, \infty)\longrightarrow R$ is a link function, which is a monotonic differentiable invertible function. The model given by (2.1) and (2.2) is called the GLM. The GLM class includes, as special cases, linear-regression and analysis-of-variance models, logit and probit models for quantal responses, log-linear models, and multinomial response models for counts (for more details, see McCaullagh and Nelder ^[42]).

Consider that the probability density function of IE distribution is as follows (see Abdel-Aty et al. ^[19]):

which has no mean and $\gamma$ is a scale parameter. The median value of the response variable is given by

Since mean does not exist, we use the median $\tilde{\mu}$ instead of it in the link function (see Das and Dey ^[43]). The cumulative function of IE distribution is given by

Let $y_i$ be a random sample from IE, and $\gamma_i = \tilde{\mu}_i\log(2)$, the log-likelihood function based on $y_i$, is given by

Regression coefficients are estimated using Fisher's scoring technique (for more details, see Nelder and Wedderburn ^[10], and McCaullagh and Nelder ^[42]). In order to develop the GLMs for our models, IE-Reg is similar to GLMs, except that the distribution of the response variable is not a member of the exponential family (Ferrari and Cribari-Neto ^[44]). We also suggest the logarithm (log) link functions of $g(.)$, in view of (2.2) as in the following lemma.

Lemma 2.1:

Let the response variable $Y$ have an IE distribution, $i = 1, 2, ..., n$, and let the link function of the form be

Thus, estimated coefficients $\hat{\beta}' = (\hat{\beta}_0, \hat{\beta}_1, \dots, \hat{\beta}_p)$ using Fisher's scoring technique at the $s$th iteration based on the IRLS process are given by

where $X$ is a covariates matrix, $\hat{\beta}^{(0)}_{j}$ is an initial vector, $Z' = (z_1, z_2, \dots, z_n)$, and

The procedure in (2.7) can be repeated until $|\hat{\beta}^{(s)}- \hat{\beta}^{(s-1)}|\leq \varepsilon$. IE-Reg model in this case is given by

Proof: See the appendices

3.   Bayesian approach

Diaconis and Ylvisker ^[43] introduced conjugate prior distribution for the exponential family, which, as in (2.1), can be shown as

where $k_{1}$ is a normalization constant, and $m, \mu_0$ are natural parameters. $\theta_i$ values are connected to the regression coefficients by link function $\eta_{i} = x^{\prime}_{i}\beta$ as

Posterior distribution of $\theta_i$ is given by

Das and Dey ^[43] suggested a Jacobian of transformation and rewrote (3.3) with term $\eta_i$, as

where $k_{2}$ is a normalization constant, and $\frac{\partial g^{\ast}(\eta_i)}{\partial \eta_i}\neq 0$. They used a zero-one loss function to attain the posterior mode of (3.4) as $\hat{\eta_i} = h(y_{i})$; hence, estimated coefficients $\hat{\beta}^{*} = (\hat{\beta}^{*}_0, \hat{\beta}^{*}_1, ..., \hat{\beta}^{*}_p)^{\prime}$ are given by

where $\hat{\beta}^{*}$ is the least-squares estimates, and $\hat{\eta}^{\prime} = (\hat{\eta}_1, \hat{\eta}_2, ..., \hat{\eta}_n)$ (for more details, see Das and Dey ^[43], and Das and Dey ^[45]).

3.1.   IE-BReg and IE-EBReg models

In order to develop a Bayesian approach, we suggest Bayesian and empirical Bayesian regression models (IE-BReg and IE-EBReg) that are similar to Bayesian GLMs, except that the distribution of the response variable is not a member of the exponential family. We used the general form of the posterior in (3.4), and since $g(\cdot)$ is a monotonic differentiable function, we then attain posterior Bayes estimates. In addition, we use log link function with zero-one and LINEX loss functions to be appropriate of Bayes estimates. IE-BReg and IE-EBReg estimators correspond to the log link function using different loss functions, as in the following lemmas.

Lemma 3.1: $\mbox{(IE-BReg and IE-EBReg models based on zero-one loss function)}$

Let the response variable $Y$ have an IE distribution, and the link function of the form be as in (2.6). Consider that $\tilde{\mu}$ has a gamma prior $G\left(\alpha, \lambda\right)$ with the following density function

Thus, the posterior mode of $\eta_i$ by using zero-one loss function can be derived by solving the following equation:

Estimated coefficients $\hat{\beta}^{*}$ are given as in (3.5). In this case, IE-BReg and IE-EBReg models are given by

In the case of IE-EBReg, empirical Bayes estimates for unknown prior distribution parameter $\lambda$ are given by MLE from the data. Therefore, IE-EBReg estimates are found by placing these estimated prior distribution parameter into Equation (3.7) by $\hat{\lambda}$.

Proof: See the appendices

Lemma 3.2

Let the response variable $Y$ have an IE distribution, and the link function of the form be as in (2.6). Consider that $\tilde{\mu}$ has a gamma prior $G\left(\alpha, \lambda\right)$ with the density function as in (3.6). Using Jacobian transformation from $\tilde{\mu}_i$ to $\eta_i$, the posterior function of $\eta_i$ can be written as in the following equation

Thus, the scale parameter of the prior $G\left(\alpha, \lambda\right)$ is less than or equal one ($\lambda \leq 1$) and its variance ($\sigma^2$) is greater than or equal the mean ($\mu$).

Proof: See the appendices

Lemma 3.3: $\mbox{(IE-BReg and IE-EBReg models based on LINEX loss function)}$

Let the response variable $Y$ have an IE distribution, and let the link function of the form be as given in (2.6). Consider that $\tilde{\mu}$ has a gamma prior with density function as given in (3.6). As a result, the posterior Bayes estimates of $\eta_i$, by using the LINEX loss function, can be derived as follows

The estimated coefficients $\hat{\beta}^{*}$ and the IE-BReg and IE-EBReg models and their estimates in this case are given as in (3.5) and (3.8), respectively. In the caes of IE-EBReg model, the regression coefficients estimates is found by placing these estimated prior distribution parameter into Equation (3.9) by $\hat{\lambda}$.

Proof: See the appendices

Lemma 3.4:

Let the response variable $Y$ have an IE distribution, and the link function of the form be as in (2.6). Consider that $\tilde{\mu}$ has a gamma prior $G\left(\alpha, \lambda\right)$ with the density function as in (3.6). Using Jacobian transformation from $\tilde{\mu}_i$ to $\eta_i$, the posterior function of $\eta_i$ can be written as in the following equation

Thus, the scale parameter of the prior $G\left(\alpha, \lambda\right)$ is greater than or equal one ($\lambda \geq 1$) and its variance ($\sigma^2$) is less than or equal the mean ($\mu$).

Proof: See the appendices

3.2.   Robust IE-BReg and IE-EBReg models

M-estimation is considered to be the most common method of robust regression. It was proposed by Huber ^[46] in the presence of outliers, and it is more efficient than ordinary least squares (OLS) (Rousseeuw and Leroy ^[47], and Chang, et al. ^[48]). The Huber's function takes the following form (Huber ^[46,49]):

where $k$ is the tuning constant, $r$ is the residual corresponding to the observation in OLS, and $\rho(\cdot)$ is the objective function that satisfies certain properties. Often, $\rho(\cdot)$ can be formed by using a linear combination of the residuals. Defining function $\frac{\partial}{\partial r}\rho(r)$ and the corresponding weight function in this case is as follows:

Another M-estimation function is the Tukey bisquare's (biweight) function. This is based on Tukey's function, taking the form of that in Sinova and Van Aelst ^[41]

where $k$ is the tuning constant and $r$ is the residual corresponding to the observation in OLS. Defining function $\frac{\partial}{\partial r}\rho(r) = \psi(r)$ and the corresponding weight function in this case is given as follows:

To make the IE-BReg and IE-EBReg models are robust, we suggest the Huber's and biweight functions for these models based on an adopted IRLS algorithm. There are also many other versions of the M-estimation function that could be used here.

Lemma 3.5: $\mbox{(IE-BReg and IE-EBReg models based on M-estimation functions)}$

Let the response variable $Y$ have an IE distribution, and let the link function of the form be as given in (2.6). Consider that $\tilde{\mu}$ has a gamma prior with density function is given as in (3.6). Using the Jacobian transformation from $\tilde{\mu}_i$ to $\eta_i$ and using the log link function, we have the posterior distribution of $\eta_i$ is given as in (3.10). Thus, the estimated coefficients $\hat{\beta}^{*} = (\hat{\beta}^{*}_0, \hat{\beta}^{*}_1, ..., \hat{\beta}^{*}_p)^{\prime}$ are given as

where $\hat{\eta}_i = h(y)$ and $\hat{\eta}^{\prime} = (\hat{\eta}_1, \hat{\eta}_2, ..., \hat{\eta}_n)$ are the posterior Bayes estimates of $\eta_i$ using the zero-one or LINEX loss functions, and $W = diag(w_{1}, w_{2}, \dots, w_{n})$, $w_{i}$ are the selected weights depending on M-estimation functions. In this case, coefficients are estimated using the adopted IRLS Algorithm.

An adopted algorithm based on IRLS and M-estimation is employed as follows:

Equation (3.15) is solved using an adopted algorithm on the basis of the standard IRLS algorithm (for more details, see Maronna et al. ^[50], Wen and Liu ^[51], and Kikuchi et al. ^[52]). This algorithm is proposed for the solution of the IE-BReg and IE-EBReg estimates, that is employed in the following steps:

(i) Setting the iteration counter at $q = 0$, finding an initial estimates of regression coefficients $\hat{\beta}_{j}^{*(q)}, \ j = 0, 1, 2, ..., p-1$ using IE-Reg estimates.

(ii) The initial residuals $r_{(i)}^{*(q)} = Y_{i}-e^{(X'_{i}\hat{\beta}_{j}^{*(q)})}$ are based on the log link function that is given as in (2.6), and calculate an initial scale estimate $s^{*(q)} = 1.4826(median|r_{i}^{*(q)}|)$.

(iii) An initial standardized residuals ${u_{i}^{*(q)}} = \frac{r_{i}^{*(q)}}{s^{*(q)}}$ are calculated and used to calculate initial estimates for the weight function. Preliminary weights are $w_{i}^{*(q)} = w(u_{i}^{*(q)})$.

(iv) Calculate $\hat{\lambda}$ as the prior distribution $\pi(\tilde{\mu})$ is $G\left(\alpha, \lambda\right)$ in the case of IE-EBReg using MLE estimates, or the scale parameter $\lambda$ is known in the case of IE-BReg model.

(v) Calculate Bayes estimates $\hat{\eta}_i = h(y); i = 1, ..., n$ using the prior $G\left(\alpha, \lambda\right)$ and zero-one or LINEX loss functions.

(vi) Using weights from Steps i-iii and Steps iv and v to find estimators in (3.15).

(vii) Set $q = q+1$; then, go to Step ii. Steps (ii) to (vii) are repeated until the estimate of $\hat{\beta}^{*(q)}$ is stabilized from the previous iteration, which means: $|\hat{\beta}^{*(q+1)}- \hat{\beta}^{*(q)}|\leq \varepsilon$.

Under regularity conditions, estimator $\hat{\beta}^{*}$ has asymptotically normal distribution $\hat{\beta}^{*} \equiv N(\beta^{*}, \left(X^{\prime}WX\right)^{-1})$ (see Houston and Woodruff ^[30], Tellinghuisen ^[53]).

4.   Data analysis

In this section, we apply the IE-Reg, IE-BReg, and IE-EBReg models for the daily confirmed COVID-19 cases in Saudi Arabia. The relevant dataset in this application is COVID-19 data from Saudi Arabia in 2020. These data contain 91 observations from March 23-June 21, 2020 in which the response variable $Y$ is the daily confirmed COVID-19 cases in Saudi Arabia. Explanatory variables are: $X_1$; daily recovered COVID-19 cases, $X_2$, daily critical COVID-19 cases; $X_3$, daily active COVID-19 cases; $X_4$, tests per million (PCR tests); $X_5$, curfew hours per day; $X_6$, maximal temperature in Celsius per day; $X_{7}$, maximal relative humidity ($\%$); $X_{8}$, maximal wind speed in miles per hour (mph); and $X_{9}$, maximal pressure in hectopascal (hPa). In the case of variables $X_i, i = 6, ..., 9$ are the average for the cities of Riyadh, Jeddah, and Dammam that have had the highest number of confirmed and death cases. This dataset was taken from the Ministry of Health of Saudi Arabia (COVID-19 dashboard) and the Ministry of Health's Twitter account.

Lemmas in Sections 2 and 3 were applied to these data. IE-Reg, IE-BReg, and IE-EBReg models based on log link and loss functions were used. Bayes coefficients were obtained using a gamma prior $G\left(\alpha, \lambda\right)$ to a known shape $\alpha$ and an unknown scale $\lambda$ parameter. We used generated data from a gamma prior distribution to estimate $\lambda$ in the case of IE-BReg model and use our data in the case of IE-EBReg model. We also compared the performance of all these models. Different plots, such as the quantile-quantile (Q-Q) plot, the empirical cumulative distribution function (ECDF), and box plot, were proposed to aid in distributional assessment and identify outliers. In addition, Huber's and biweight functions are suggested in the case of Lemma (3.5) to avoid such distortions due to outliers (for more details, see Sinova and Van Aelst ^[41]). The backward-selection method is used in the IE-EBReg model to remove the input variable (see Table 3). Modeling performance is measured in terms of some criteria, such as AIC, BIC, D, D/df (divided by its degrees of freedom), and mean square error (MSE) (de Jong and Heller ^[12]). We also used Thiel's inequality coefficient (TIC) to compare the prediction accuracy of the selected models (Leuthold ^[54], and Niu et al. ^[55]). To check the adequacy for the selected models, we consider the deviance residuals (see McCaullagh and Nelder ^[42]). The predictive results of these models and other numerical results are shown in Tables 1-8. R software was used to carry out calculations. In order to compare with known distributions, the glm() function in "stats" was used to fit the GLMs (Faraway ^[56]). Functions qqPlot(), ecdf(), boxplot, and ks.test() in R package "stats" were used for the assessment distributions (Fox and Weisberg ^[57]).

K-S test is used to recommend the IE distribution as a good fit for the daily confirmed COVID-19 cases in Saudi Arabia comparing with some other distributions as given below:

Based on the results obtained from the above table, the large $p$-value for the test indicated that IE distribution fit the response variable in the given data quite well. Figure 1 provides the Q-Q plot, ECDF, and fitted functions of the selected models, and it is clear that IE distribution fits these data well.

Figure 2 presents box plots corresponding to each of the Saudi Arabia COVID-19 dataset variables, and the chart also maps outliers that exceeded the values of fences. The plot also displays the maximum, minimum, and median of the data, along with the first and third quantile. The outliers could also be identified as in Figure 2 (shown as unfilled circles) in explanatory variable $X_2$, and we can observe a big difference between the maximal value and the rest of the observation, which is beyond the outer fence. This case could be considered a very extreme outlier. The chart also shows more outliers on variables $X_3$, $X_5$, and $X_8$. For the $X_2$ and $X_3$ box plots, variables having outliers were much higher than the third quartile or much lower than the first quartile of the box was. Possible reasons for this deviation may be when daily active or critical cases exceeded the recovered cases, or daily recovered cases exceeded active or critical cases. There is no outlier in the daily-confirmed-cases variable in the context of box-plot analysis, and this variable was slightly asymmetrically distributed with a relatively heavy tail.

Table 1 presents that the MSE value of the IE-Reg, gamma models was large, while the MSE values of the inverse Gaussian was very large. Table 2 shows that the MSE value of the IE-Reg model was large, while the MSE values of the IE-BReg and IE-EBReg models were small. In the case of IE-Reg model, $\hat{\beta}^{(s)}$ was stabilized when Fisher's scoring procedure converged at $s = 15$ because of $|\hat{\beta}^{(15)}- \hat{\beta}^{(14)}| < 0$. In the case of the IE-BReg and IE-EBReg models, the fitting results based on Huber's function had the largest MSE values compare to the MSE values based on biweight function. The results in this table also show that the MSE values based on biweight function and without use any weight function were almost similar, but fitting results based on biweight function were best in terms of AIC, and D statistics. Comparing between Bayesian fitting results, we observed that the results based on LINEX loss function were better than those of the zero-one loss function. Gamma prior variance was larger than the mean in the case of used the zero-one loss function, while gamma prior variance was less than the mean in the case of used the LINEX loss function.

Table 2 also shows that the IE-EBReg model based on gamma prior $G(2.610, 1.420)$ and biweight function was the best on our dataset, since D and D/df $(D/df = 0.380), df = 83$ were acceptable with a low MSE compared to that of the other models. Tables 1 and 2 show that the $D/df$ of this model was very close to 1, indicating that the fitting degree was very good. The value of z statistic was large, so there was also a significant relationship among the variables.

Through this comparison, we can conclude that the predictive model was IE-EBReg based on biweight and LINEX functions, and it is given as follows:

In the side of the adopted IRLS procedure based on M-estimation, $\hat{\beta}^{*(s)}$ is stabilized and it converged at $s = 15$ because of $|\hat{\beta}^{*(15)}- \hat{\beta}^{*(14)}| < 4.2e-06$. The performance of the selected model was indicative that the adopted algorithm works well. The actual and predicted for the IE-EBReg model compared to IE-Reg are represented graphically in Figure 3. Figure 4 show that the deviance residuals against the indices of the observations suggest that the residuals for the IE-EBReg model are randomly scattered around zero. From Table 2 and 3, we can conclude that, for the selected model among possible IE-EBReg models based on the smallest AIC = 54.151 and lowest MSE = 0.230 with eight variables, this fitting result was the best when using level of significance $\alpha = 0.05$.

Table 4 presents the fitting results for the IE-EBReg model based on biweight function that is the best based on relative errors. For the prediction results, Table 5 shows that the IE-EBReg model based on Huber's function compared to that based on biweight function had the worst prediction accuracy because the TIC value was closer to 1, while the TIC value for the IE-EBReg model based on biweight function was closer to 0.

According to the residuals results of the IE-BReg and IE-EBReg models, it can be found that the errors of the first data are too large, which seriously affects of results, see Table 6 (Case 1). To make a precision comparison, we removed one observation form data $i = 1$ that have large error and we reestimated the model after removing from the data (Case 2). Furthermore, we replaced $y_{1}$ with the mean of observation from $i = 2$ to $i = 11$ and we reestimated the model after modifying the data (Case 3). We also replaced $y_{i}, i = 1, 4, 5, 6, 7, 9$ with the mean of observation from $i = 2$ to $i = 11$ and we reestimated the model after modifying the data (Case 4). However, we removed the observations $i = 1, 4, 5, 6, 7, 9$ and we reestimated the model after removing observations (Case 5). The relative changes in the parameter estimates are presented in Table 6. The the deviance residuals against the indices of the observations suggest that the residuals are randomly scattered around zero at level of significance $\alpha = 0.001$ (Cox Stuart test) for the IE-BReg and IE-EBReg models based on the priors $G(2.3537, 1.410)$ and $G(2.3537, 1.493)$.

Table 6 and 7 show that, the IE-EBReg model (Case 4) is the best for our data, and it is given as follows:

In the side of the adopted IRLS procedure based on M-estimation, $\hat{\beta}^{*(s)}$ is stabilized and it converged at $s = 15$ because of $|\hat{\beta}^{*(15)}- \hat{\beta}^{*(14)}| < 0$. The actual and predicted for the IE-EBReg model compared to IE-Reg model is represented graphically in Figure 5. From Table 8, we can conclude that, this model has smallest AIC = 50.6634 and a low MSE = 0.2014, and there was also a significant relationship among the variables when using level of significance $\alpha = 0.05$. Figure 6 show that the deviance residuals against the indices of the observations suggest that the residuals for the IE-EBReg model (Case 4) are randomly scattered around zero. In comparison between IE-EBReg models that were shown in (4.1) and (4.2), and based on the results on Table 2 and 6, we can conclude that the IE-EBReg model in (4.2) is the best for our data.

Figures 3 shows that the IE-EBReg model based on biweight and LINEX loss function generally fits the dataset better than other models. This also clear from Figure 4 as the plot of the deviance residuals against the indices of the observations suggest that the residuals are randomly scattered around zero at level of significant 0.001.

Figure 5, shows that the IE-EBReg model (Case 4) which based on biweight and LINEX function fits the dataset better. Also, Figure 6 shows the deviance residuals against the indices of the observations suggest that the residuals are randomly scattered around zero at level of significant 0.001.

5.   Conclusions

In this paper the regression models IE-Reg, IE-BReg and IE-EBReg for modeling the daily confirmed Saudi COVID-19 cases with some envionmental-related variables (covariates) are considered. Zero-one and LINEX loss functions were used to attain the Bayesian and empirical Bayesian estimates based on the log-link function. In a non-Bayesian approach, parameter estimation is done by the Fisher scoring technique, and closed-form expressions are provided for the score function, and for Fisher's information matrix and its inverse.In the Bayesian approach, parameter estimation is performed using a gamma prior distribution, Jacobian transformation, and least-squares estimates. The IE-Reg, IE-BReg, and IE-EBReg models were compared to find which model predicted better. To deal with outlier problems, IE-BReg based on Huber's and biweight functions, and the adopted algorithm based on IRLS to find the estimates, were proposed. For distributional assessment, Q-Q, ECDF, box plots, and the KS test were applied. Some criteria, namely, AIC, BIC, D, D/df, and MSE, were also computed for all regression models.

According to the results of the application, it was concluded that IE-BReg and IE-EBReg with a log link function performed the best in terms of the AIC, BIC, D, D/df, and MSE statistics, so they are recommended for these data. In contrast, IE-Reg showed poor results compared with those of the other models. Results indicated that the IE-EBReg model is highly capable of improving the performance of regression models to a greater extent in the prediction of daily confirmed COVID-19 cases in Saudi Arabia. Finally, it is found the following regressors are significant for the model: Explanatory variables are: X1, daily recovered COVID-19 cases; X3, daily active COVID-19 cases; X4, tests per million (PCR tests); X5, curfew hours per day; X6, maximal temperature in Celsius per day; X7, maximal relative humidity ($\%$); and X9, maximal pressure in hectopascal (hPa).

Acknowledgments

The authors would like to thank the editor and referees for their helpful comments, which improved the presentation of the paper. Also, the authors would like to extend their sincere appreciation to the Deanship of Scientific Research at King Saud University for its funding this Research Group (RG -1435-056).

Conflict of interest

The authors declares that they have no conflicts of interest.

Data availability

The data is freely available at: https://covid15.moh.gov.sa/ https://twitter.com/SaudiMOH?s=20/ https://www.worldometers.info/coronavirus/ https://www.wunderground.com/.

Appendices

5.1.   Proof of Lemma 2.1:

Suppose that, in $y_i\sim f(y_i; \gamma)$, as in (2.3), the log-likelihood function based on $y_i, \ \ i = 1, 2, ..., n$ is given as in (2.5). Link function $\log(\tilde{\mu}_i)$ connecting the $\tilde{\mu}_i$ with linear model $x'_i\beta$ in this case is given as in (2.6). Score function $U_r$ for log likelihood is written from one observation as

From (2.5) and (2.6), we have

which can be written in matrix notation as

Taking the second derivatives of $l(\beta)$, we have

hence,

Since $E\left(\frac{1}{y_i}\right) = \frac{1}{\tilde{\mu}_{i}\log(2)}$ and $\frac{\partial \tilde{\mu}_i} {\partial \beta_j} = x_{ij}\tilde{\mu}_i$, then $I_{rj} = -E\left(\frac{\partial U_r(\beta)}{\partial \beta_j}\right) = \sum^n_i x_{ij} x_{ir},$ and Fisher's information matrix is given as

where $W\left(\hat{\beta}\right)$ is the unit matrix. Fisher's scoring process can be applied to obtain

From (5.3) and (5.4), we have

and

Thus, the estimated coefficients $\hat{\beta}$ are given by

as given in (2.7). To derive the MLS of $\beta$, IRLS algorithm is used. Under regularity conditions on the likelihood function, the MLE $\hat{\beta}^{(s)}$ is asymptotically normal, unbiased, and efficient, with covariance matrix equal to the inverse of Fisher's information matrix (Houston and Woodruff ^[30]). Therefore, asymptotically,

where $\left(X^{\prime}WX\right)^{-1}$ is the inverse of Fisher's information matrix.

5.2.   Proof of Lemma 3.1:

Suppose that $y_i\sim f(y_i; \gamma)$ is as in (2.3), $\gamma = \tilde{\mu}\log(2)$; then, the density function of $y_{i}$ is given by

Consider a gamma prior for $\tilde{\mu}_i$, which can be written as in (3.6). Posterior distribution of $\tilde{\mu}_i$ is given by

Using Jacobian transformation from $\tilde{\mu}_i$ to $\eta_i$, and using the log link function that is given as in (2.6), we have

where $g^{\ast}(\eta_i) = e^{\eta_i} = \tilde{\mu}_i$ and $\frac{\partial g^{\ast}(\eta_i)}{\partial\eta_i} = e^{\eta_i}$. Then,

Taking the derivative of the log posterior, we have

hence, the posterior mode of $\eta_i$ is given as in (3.7). Thus, the estimated coefficients $\hat{\beta}^{*}$, IE-BReg, and IE-EBReg models in this case, are given as in (3.5) and (3.8), respectively.

In the case that prior distribution parameter $\lambda$ is unknown, for this estimation task, $\hat{\lambda}$ estimates are obtained via numerical maximization of the following marginal likelihood (Shao ^[37], Dikici et al. ^[58], and Coluccia et al. ^[59])

Thus, the IE-EBReg estimate is found by placing these estimated prior distribution parameter into Equation (3.7) by $\hat{\lambda}$.

5.3.   Proof of Lemma 3.2:

Suppose that $y_i\sim f(y_i; \gamma)$ is as in (2.3), $\gamma = \tilde{\mu}\log(2)$; then, the density function of $y_{i}$ is given as in (5.6). Consider a gamma prior $G(\alpha, \lambda)$ for $\tilde{\mu}_i$, which can be written as in (3.6). Posterior distribution of $\tilde{\mu}_i$ is given as in (5.7). Using Jacobian transformation from $\tilde{\mu}_i$ to $\eta_i$, zero-one loss function and using the log link function as in (2.6), we have

Because of $0 < \alpha < \infty$, thus,

hence, we get

Since $0 < y_{i} < \infty$ for every $i = 1, ..., n$, thus

From (5.14) and (5.15), we find $\lambda \leq \frac{1}{e^{\hat{\eta}_i}}$, but $e^{\hat{\eta}_i}\leq 1$, because of equation (5.13) and the fact $e^{\hat{\eta}_i} < e^{\hat{\eta}_i}\left(\lambda+\frac{\log(2)}{y_i}\right)$. Now, by using $\frac{1}{e^{\hat{\eta}_i}}\geq 1$ and $\lambda \leq \frac{1}{e^{\hat{\eta}_i}}$, then we obtain $\lambda \leq 1$. Hence, $\frac{1}{\lambda}\geq 1$ and $\frac{\alpha}{\lambda^2}\geq \frac{\alpha}{\lambda}$, thus the variance of $\tilde{\mu}$ is greater than or equal the mean.

5.4.   Proof of Lemma 3.3:

Suppose that $y_i\sim f(y_i; \gamma)$ is as it is in (2.3), $\gamma = \tilde{\mu}\log(2)$, and the density function of $y_{i}$ is given as in (2.3). Consider $\tilde{\mu}_i$ has a gamma prior with density function, which can be written as in (3.6). Using the posterior distribution of $\eta_i$ that is given in (5.9), we have

and

Using the LINEX loss function, we have

Thus, the posterior Bayes estimates of $\eta_i$ by using the LINEX loss function are given as in (3.9).

In the case that prior distribution parameter $\lambda$ is unknown for this estimation task, $\hat{\lambda}$ estimates are obtained via numerical maximization of the marginal likelihood of Equation (5.11). As a result, the IE-EBReg estimates are found by placing these estimated prior distribution parameters into Equation (3.9) by $\hat{\lambda}$.

5.5.   Proof of Lemma 3.4:

Suppose that $y_i\sim f(y_i; \gamma)$ is as in (2.3), $\gamma = \tilde{\mu}\log(2)$; then, the density function of $y_{i}$ is given as in (5.6). Consider a gamma prior $G(\alpha, \lambda)$ for $\tilde{\mu}_i$, which can be written as in (3.6). Posterior distribution of $\tilde{\mu}_i$ is given as in (5.7). Using Jacobian transformation from $\tilde{\mu}_i$ to $\eta_i$, LINEX loss function and using the log link function as in (2.6), we have

Because of $0 < \alpha < \infty$, thus $0 < e^{-\alpha\hat{\eta}_i} < 1$, and

Using (5.20) and (5.21), we get

hence

and,

Thus, $\lambda < \frac{\log(2)}{y_i}$ and $\left(\lambda^{\alpha-1}-1\right) < \frac{\log(2)}{y_i}$. Hence, $\frac{1}{\lambda} < 1+\frac{\log(2)}{y_i}$. Since $0 < y_{i} < \infty$ for every $i = 1, ..., n$, thus $1 < 1+ \frac{\log(2)}{y_i} < \infty$. Now, by using $\frac{1}{\lambda} < 1+\frac{\log(2)}{y_i}$ and $1 < 1+\frac{\log(2)}{y_i}$, then we obtain $\frac{1}{\lambda} \leq 1$, $\lambda \geq 1$ and $\frac{\alpha}{\lambda^2}\leq \frac{\alpha}{\lambda}$. Thus, the variance of $\tilde{\mu}$ is smaller than or equal the mean.