The gamma-Maxwell regression for COVID-19 mortality rates of the 50 U.S. largest cities

Abstract

A new parametric regression model is developed based on the gamma-Maxwell distribution. Monte Carlo simulations show the accuracy of the maximum likelihood estimators. The proposed model explains COVID-19 mortality rates of the 50 U.S. largest cities.

Keywords

COVID-19 gamma-G family maximum likelihood Maxwell distribution regression model simulation

1. Introduction

In late 2019, a disease called SARS-CoV-2, known as COVID-19, spread and became a pandemic. A wide range of symptoms were reported, and people experienced difficulty breathing, fever, loss of taste or smell, sore throat, among others.1

After the outbreak, countless statistical applications were performed to pandemic data in the U.S. and other countries. A new quantile regression was constructed by Gallardo et al. (2022) for coronavirus mortality rates in the U.S. states. The article considered social, demographic and health variables that affected the first-wave mortality rate. Karmakar et al. (2021) explained COVID-19 incidence and mortality rates in U.S. in terms of county-level social and demographic variables. Moreover, Biazatti et al. (2022) presented applicability of the Weibull beta-prime distribution to three real COVID-19 data sets. In Brazil, Sousa et al. (2020) identified the risk factors associated with mortality and survival of COVID-19 cases in the Northeast region. In conclusion, the mortality is similar to countries with numerous of cases, with elderly people, and comorbidities presented a greater risk of death. For an extensive study, Duhon et al. (2021) estimated the initial growth rate of COVID-19 wide-world. They studied the association between the initial growth rate and non-pharmaceutical interventions as well as pre-existing country characteristics. Other related works can be found in Janke et al. (2021) and Stokes et al. (2021).

In this context, our main goal is a construction of a new regression model based on the gamma-Maxwell (GM) (Iriarte et al., 2017) distribution, which has advantages over other competitive models. We study the coronavirus mortality rates of the 50 U.S. largest cities to show the flexibility of the new regression. Presently, the USA has the highest number of COVID-19 deaths worldwide, more than one million.2 In addition to the geographical proportions, high rates of obesity and cardiac disease, the lack of a universal health care system hindered the American response to the pandemic. The application for regression modeling is done in which we investigate the influence of explanatory variables on the mortality rates from COVID-19. The paper aims to contribute to the literature with a new regression model as well as check the influence of disease in the USA. The behavior of the data set is illustrated in Fig. 1.

Figure 1.

Histogram and empirical density of COVID-19 mortality rates of the 50 U.S. largest cities.

Regression models are necessary when there are explanatory variables associated with the response variable. We present some global influence measures to assure the choice of the best fit regression model to the data set. Likewise, we establish residual analysis based on deviance residuals jointly with simulated envelopes.

The article is structured as follows. Section 2 reviews the GM distribution (Iriarte et al., 2017). Section 3 provides new mathematical properties. Section 4 determines the maximum likelihood estimates (MLEs), and provides a simulation study to verify their accuracy. Section 5 develops the GM regression model, and examines the consistency of the estimators. Section 6 defines deviance residuals for the fitted regression model, and investigates some diagnostic measures. Section 7 proves the advantage of the new regression model compared to other competitive models, and selects the best model for explaining the COVID-19 mortality rates of the 50 U.S. largest cities. The results indicate randomly residuals inside the simulated envelope, which is expected for a good fit. It also supports and presents some valuable findings. Finally, Section 8 ends with some remarks and suggests guidelines for future research.

2. Definition

The cumulative distribution function (cdf) $G(x)=G(x;\theta)$ of a Maxwell random variable $X$ has the form

$\displaystyle G(x)=\gamma_{1}\left(\frac{3}{2},\frac{x^{2}}{\theta^{2}}\right)% ,\,x>0,$ (1)

where $\theta>0$ is a scale parameter, $\gamma_{1}(a,z)=\Gamma(a)^{-1}\,\int_{0}^{z}u^{a-1}\,{\mathrm{e}}^{-u}\mathrm{% d}u$ is the incomplete gamma function ratio, and $\Gamma(a)=\int_{0}^{\infty}u^{a-1}\,{\mathrm{e}}^{-u}\mathrm{d}u$ is the gamma function.

The probability density function (pdf) of $X$ becomes

$\displaystyle g(x)=\frac{4}{\sqrt{\pi}}\,\frac{x^{2}}{\theta^{3}}\,{\rm e}^{-% \frac{x^{2}}{\theta^{2}}}.$ (2)

A random variable $Y\sim\text{GM}(\theta,\alpha)$ has the GM distribution (Castellares & Lemonte, 2015; Iriarte et al., 2017), where $\theta>0$ is a scale, and $\alpha>0$ is a shape parameter, if its cdf and pdf are (for $y>0$ )

$\displaystyle F(y)=\gamma_{1}\left\{\alpha,-\log\left[1-\gamma_{1}\left(\frac{% 3}{2},\frac{y^{2}}{\theta^{2}}\right)\right]\right\}$ (3)

and

$\displaystyle f(y)=\frac{4}{\sqrt{\pi}\,\Gamma(\alpha)}\,\frac{y^{2}}{\theta^{% 3}}\,{\rm e}^{-\frac{y^{2}}{\theta^{2}}}\left\{-\log\left[1-\gamma_{1}\left(% \frac{3}{2},\frac{y^{2}}{\theta^{2}}\right)\right]\right\}^{\alpha-1},$ (4)

respectively.

For $\alpha=1$ , the GM distribution is equal to the Maxwell. Figure 2 provides plots of the density, cumulative, hazard rate and survival functions of $Y$ for selected parameter values.

Figure 2.

GM distribution. (a) Pdf. (b) Cdf. (c) Hazard rate. (d) Survival function.

3. Properties

Due to the impossibility of obtaining closed-form mathematical properties of the GM distribution, we find a linear representation for its density function.

The exponentiated-Maxwell random variable $T_{a}\sim\exp-M(a)$ with power parameter $a>0$ has density $w_{a}(y)=a\,g(y)\,G(y)^{a-1}$ determined from Eqs (1) and (2) as

$\displaystyle w_{a}(y)=\frac{4a}{\sqrt{\pi}}\frac{y^{2}}{\theta^{3}}\,{\rm e}^% {-\frac{y^{2}}{\theta^{2}}}\gamma_{1}^{a-1}\left(\frac{3}{2},\frac{y^{2}}{% \theta^{2}}\right).$

Theorem 1. The GM density Eq. (4) admits the linear representation

$\displaystyle f(y)=\sum_{k=0}^{\infty}\frac{\varphi_{k}(a)}{(a+k)}\,w_{a+k}(y),$ (5)

where

$\displaystyle w_{a+k}(y)=\frac{4(a+k)}{\sqrt{\pi}}\,\frac{y^{2}}{\theta^{3}}\,% {\rm e}^{-\frac{y^{2}}{\theta^{2}}}\gamma_{1}^{a+k-1}\left(\frac{3}{2},\frac{y% ^{2}}{\theta^{2}}\right),$

is the density of $T_{a+k}$ , and $\varphi_{k}(a)$ is defined below.

Proof. The proof is straightforward. The pdf of $Y$ can be written as Eq. (5) from the main result in Castellares and Lemonte (2015), where $w_{a+k}(y)$ is the pdf of $T_{a+k}$ , $\varphi_{0}(a)=1/\Gamma(a)$ , $\varphi_{k}(a)=\frac{(a-1)}{\Gamma(a)}\psi_{k-1}(k+a-2)\,(k\geqslant 1)$ , and $\psi_{k}(x)$ are the Stirling polynomials: $\psi_{0}(x)=1$ , $\psi_{1}(x)=(x+1)/2$ , $\psi_{2}(x)=(3x+2)(x+1)/12$ , $\psi_{3}(x)=x(x+1)^{2}/8$ , etc. ∎

Hence, Eq. (5) is very useful to find structural properties of $Y$ from known properties of $T_{a}$ .

3.1 Moments

The $n$ th ordinary moments of $Y$ follows from Theorem 1

$\displaystyle\mu_{n}^{\prime}=E(Y^{n})=\sum_{k=0}^{\infty}\,\frac{\varphi_{k}(% a)}{(a+k)}\,E(T_{a+k}^{n}).$

Theorem 2. The $n$ th ordinary moment of $Y$ can be expressed as

$\displaystyle\mu_{n}^{\prime}=\sum_{k=0}^{\infty}\,\frac{\varphi_{k}(a)}{(a+k)% }\,\sum_{m,i=0}^{\infty}p_{m,i},$ (6)

where

$\displaystyle p_{m,i}=\frac{2(a+k)\theta^{3k+2i+n}}{\sqrt{\pi}}\,\Gamma(3(k+1)% /2+i+n/2)\,s_{m}(a+k-1)\,c_{m,i},$

and the quantities $s_{m}(a+k-1)$ and $c_{m,i}$ are defined below.

Proof. The $n$ th moment of $T_{a+k}$ has the form

$\displaystyle E(T_{a+k}^{n})=\frac{4(a+k)}{\sqrt{\pi}\theta^{3}}\,\int_{0}^{% \infty}\,x^{n+2}\,{\rm e}^{-\frac{y^{2}}{\theta^{2}}}\gamma_{1}^{a+k-1}\left(% \frac{3}{2},\frac{x^{2}}{\theta^{2}}\right)\,dx.$ (7)

The incomplete gamma function ratio admits the power series (Prataviera et al., 2020)

$\displaystyle\gamma_{1}^{a+k-1}\left(\frac{3}{2},\frac{x^{2}}{\theta^{2}}% \right)=\sum_{m=0}^{\infty}s_{m}(a+k-1)\,\gamma_{1}^{m}\left(\frac{3}{2},\frac% {x^{2}}{\theta^{2}}\right),$ (8)

where $s_{m}(a+k-1)=\sum_{j=m}^{\infty}(-1)^{j+m}\,\binom{a+k-1}{j}\,\binom{j}{m}$ .

Likewise, from Prataviera et al. (2020), the power series holds

$\displaystyle\gamma_{1}\left(\frac{3}{2},\frac{x^{2}}{\theta^{2}}\right)=\sum_% {i=0}^{\infty}a_{i}\,x^{2i+3},$ (9)

where $a_{i}=a_{i}(\theta)=\frac{(-1)^{i}}{(3/2+i)\theta^{2i+3}\Gamma(3/2)i!}$ (for $i\geqslant 0$ ).

Applying Eq. 0.314 from Gradshteyn and Ryzhik (2014) leads to

$\displaystyle\left(\sum_{i=0}^{\infty}a_{i}\,z^{i}\right)^{m}=\sum_{i=0}^{% \infty}c_{m,i}\,z^{i},$ (10)

where $c_{m,i}=c_{m,i}(\theta)$ can be determined from the recurrence relation (for $i\geqslant 1$ )

$\displaystyle c_{m,i}=(i\,a_{0})^{-1}\sum_{r=1}^{i}[(m+1)r-i]\,a_{r}\,c_{m,i-r},$

and $c_{m,0}=a_{0}^{m}$ .

From Eq. (10),

$\displaystyle\gamma_{1}^{m}\left(\frac{3}{2},\frac{x^{2}}{\theta^{2}}\right)=% \sum_{i=0}^{\infty}c_{m,i}\,x^{2i+3m}.$ (11)

Combining Eqs (11) and (8) leads to

$\displaystyle E(T_{a+k}^{n})=\frac{4}{\sqrt{\pi}\theta^{3}}\,\sum_{m,i=0}^{% \infty}(a+k)\,s_{m}(a+k-1)\,c_{m,i}\,\int_{0}^{\infty}x^{n+2+2i+3m}\,{\rm e}^{% -\frac{y^{2}}{\theta^{2}}}\,dx.$

By solving the integral, Eq. (7) becomes

$\displaystyle E(T_{a+k}^{n})=\sum_{m,i=0}^{\infty}p_{m,i},$

where

$\displaystyle p_{m,i}=\frac{2(a+k)\theta^{3m+2i+n}}{\sqrt{\pi}}\,\Gamma(3(m+1)% /2+i+n/2)\,s_{m}(a+k-1)\,c_{m,i}.$

Finally,

$\displaystyle\mu_{n}^{\prime}=E(Y^{n})=\sum_{k=0}^{\infty}\,\frac{\varphi_{k}(% a)}{(a+k)}\,\sum_{m,i=0}^{\infty}p_{m,i}.$ (12)

∎

3.2 Generating function

The generating function (gf) of $Y$ , say $M_{Y}(\cdot)$ , can be determined from the gf of the exp-Maxwell $(a+k)$ random variable. It follows from Theorem 1

$\displaystyle M_{Y}(-t)=E\left(\mathrm{e}^{-tY}\right)=\sum_{k=0}^{\infty}\,% \frac{\varphi_{k}(a)}{(a+k)}\,M_{T_{a+k}}(-t).$

Theorem 3. The gf of the GM distribution has the form

$\displaystyle M_{Y}(-t)=\sum_{k=0}^{\infty}\frac{\varphi_{k}(a)}{(a+k)}\,\sum_% {m,i=0}^{\infty}t_{m,i}\,J_{s}(t),$ (13)

where

$\displaystyle t_{m,i}=\frac{2(a+k)}{\sqrt{\pi}\theta^{3}}\sum_{m,i=0}^{\infty}% s_{m}(a+k-1)\,c_{m,i},$ $\displaystyle J_{s}(t)=(-1)^{s}\sqrt{\pi}\frac{\partial^{s}}{\partial t^{s}}% \left[{\rm e}^{\frac{t^{2}\theta^{2}}{4}}\text{erfc}\left(\frac{t}{2\sqrt{% \theta^{-2}}}\right)\right],\quad s=2(i+1)+3m,$

and all other quantities are defined below.

Proof. The gf $M_{Y}(-t)$ can be expressed as

$\displaystyle M_{Y}(-t)=\frac{4(a+k)}{\sqrt{\pi}\theta^{3}}\int_{0}^{\infty}x^% {2}\,{\rm e}^{-tx-\frac{x^{2}}{\theta^{2}}}\,\gamma_{1}^{a+k-1}\left(\frac{3}{% 2},\frac{x^{2}}{\theta^{2}}\right)\,dx.$

Inserting Eqs (11) and (8) in the last expression

$\displaystyle M_{Y}(-t)=\frac{4(a+k)}{\sqrt{\pi}\theta^{3}}\,\sum_{m,i=0}^{% \infty}s_{m}(a+k-1)\,c_{m,i}\,\int_{0}^{\infty}x^{2+2i+3m}\,{\rm e}^{-tx-\frac% {x^{2}}{\theta^{2}}}\,dx.$

The above integral can be written, based on Eq. (2.3.15.7) in Prudnikov and Marichev (1986), as

$\displaystyle J_{s}(t)=\int_{0}^{\infty}x^{2(i+1)+3m}\,{\rm e}^{-tx-\frac{x^{2% }}{\theta^{2}}}\,dx=\frac{(-1)^{s}}{2}\,\sqrt{\pi\theta^{2}}\,\frac{\partial^{% s}}{\partial t^{s}}\left[{\rm e}^{\frac{t^{2}\theta^{2}}{4}}\,\text{erfc}\left% (\frac{t}{2\sqrt{\theta^{-2}}}\right)\right],$

where $s=2(i+1)+3m$ , and $\text{erfc}(w)$ is the complementary error function.

Hence,

$\displaystyle M_{Y}(-t)=\sum_{m,i=0}^{\infty}t_{m,i}\,J(t,s),$

where

$\displaystyle t_{m,i}=\frac{2(a+k)}{\sqrt{\pi}\theta^{3}}\sum_{m,i=0}^{\infty}% s_{m}(a+k-1)\,c_{m,i}$

and

$\displaystyle J_{s}(t)=(-1)^{s}\sqrt{\pi}\frac{\partial^{s}}{\partial t^{s}}% \left[{\rm e}^{\frac{t^{2}\theta^{2}}{4}}\text{erfc}\left(\frac{t}{2\sqrt{% \theta^{-2}}}\right)\right].$

Finally the gf of $Y$ reduces to

$\displaystyle M_{Y}(-t)=\sum_{k=0}^{\infty}\frac{\varphi_{k}(a)}{(a+k)}\,\sum_% {m,i=0}^{\infty}t_{m,i}\,J_{s}(t).$ (14)

∎

4. Estimation

The log-likelihood function for $\bm{\eta}=(\theta,\alpha)^{\text{T}}$ given a random sample $y_{1},\ldots,y_{n}$ from $Y\sim\text{GM}(\theta,\alpha)$ reduces to

$\displaystyle l_{n}(\bm{\eta})=n\log\left(\frac{4}{\sqrt{\pi}}\right)-n\log(% \Gamma(n))+\frac{1}{\theta^{3}}\sum_{i=1}^{n}y_{i}^{2}-\frac{1}{\theta^{2}}% \sum_{i=1}^{n}y_{i}^{2}-(\alpha-1)\sum_{i=1}^{n}\log\left[1-\gamma_{1}\left(% \frac{3}{2},\frac{y_{i}^{2}}{\theta^{2}}\right)\right].$ (15)

The MLEs can be found by solving the equations given in Iriarte et al. (2017), or $l_{n}(\bm{\eta})$ can be maximized numerically using the optim routine in R Core Team (2013).

We study the accuracy of these estimators by generating $N=$ 1,000 samples of sizes $n=$ 50, 100, 250, 500, 750, 1000 for the true parameters (2.75, 2.00) (scheme 1), and (1.85, 0.95) (scheme 2).

Table 1

Simulations from schemes 1 and 2

Scheme 1
Par	$n=$ 50			$n=$ 100			$n=$ 250
	AE	Bias	MSE	AE	Bias	MSE	AE	Bias	MSE
$\theta$	2.68	$-$ 0.07	0.07	2.70	$-$ 0.05	0.03	2.72	$-$ 0.03	0.01
$\alpha$	2.17	0.17	0.25	2.11	0.11	0.10	2.06	0.06	0.03
Par	$n=$ 500			$n=$ 750			$n=$ 1000
$\theta$	2.72	$-$ 0.03	0.01	2.73	$-$ 0.02	0.00	2.73	$-$ 0.02	0.00
$\alpha$	2.05	0.05	0.02	2.03	0.03	0.01	2.03	0.03	0.01
Scheme 2
Par	$n=$ 50			$n=$ 100			$n=$ 250
	AE	Bias	MSE	AE	Bias	MSE	AE	Bias	MSE
$\theta$	1.66	$-$ 0.19	0.06	1.68	$-$ 0.17	0.04	1.69	$-$ 0.16	0.03
$\alpha$	1.12	0.17	0.08	1.09	0.14	0.04	1.06	0.11	0.02
Par	$n=$ 500			$n=$ 750			$n=$ 1000
$\theta$	1.69	$-$ 0.16	0.03	1.69	$-$ 0.16	0.03	1.70	$-$ 0.15	0.03
$\alpha$	1.06	0.11	0.02	1.05	0.10	0.01	1.05	0.10	0.01

Table 1 reports the average estimates (AEs), biases and mean square errors (MSEs) for each scheme. The AEs converge to the true parameters, and the biases and MSEs tend to zero when $n$ increases, thus the consistency criterion holds.

5. The GM regression model

5.1 Definition

The systematic component of the GM regression model considers that the parameter $\theta$ in Eq. (4) varies across observations (for $i=1,\ldots,n$ ) as

$\displaystyle g(\theta_{i})=\exp(\bm{x}_{i}^{\text{T}}\bm{\beta}),$ (16)

where $g(\cdot)$ is a twice continuously differentiable log-linear link function, and $\bm{\beta}=(\beta,\ldots,\beta_{p})^{\text{T}}$ is the parameter vector of dimension $p$ associated with the explanatory variables $\bm{x}_{i}^{\text{T}}=(x_{i1},\ldots,x_{ip})$ . It is assumed that the components of $\bm{\beta}$ are independent.

Consider a sample $(y_{1},\bm{x_{1}}),\ldots,(y_{n},\bm{x_{n}})$ of $n$ independent observations. The log-likelihood function for the vector of parameters $\bm{\eta}=(\beta^{\text{T}},\alpha)^{\text{T}}$ from regression model Eq. (16) has the form

$\displaystyle l(\bm{\eta})=n\log\left(\frac{4}{\sqrt{\pi}}\right)-n\log(\Gamma% (n))+\frac{1}{\theta_{i}^{3}}\sum_{i=1}^{n}y_{i}^{2}-\frac{1}{\theta_{i}^{2}}% \sum_{i=1}^{n}y_{i}^{2}-(\alpha-1)\sum_{i=1}^{n}\log\left[1-\gamma_{1}\left(% \frac{3}{2},\frac{y_{i}^{2}}{\theta_{i}^{2}}\right)\right].$ (17)

The parameter vector $\bm{\eta}=(\beta^{\text{T}},\alpha)^{\text{T}}$ in model Eq. (16) can be estimated by numerical maximization via the optim routine in R Core Team (2013). Additionally, to compare nested models of the GM regression model (when $\alpha=1$ , Maxwell regression model), we adopt the likelihood ratio (LR) statistic.

5.2 Simulation study for the regression model

The accuracy of the MLEs in the GM regression model is evaluated based on the measures (Cordeiro et al., 2017): bias, mean square error (MSE), estimated average length (AL), and coverage probability (CP). We generate $N=$ 1,000 samples of sizes $n=$ 50, 55, …, 1000 from the GM quantile in Iriarte et al. (2017) for $\alpha=$ 2.00, $\beta_{0}=$ 1.50 and $\beta_{1}=$ 0.80.

Figure 3.

Plots of the measure values for the parameters. (a–c) Biases. (d–f) MSEs. (g–i) ALs.

Figure 4.

Plots of the CPs for the parameters.

Figures 3 and 4 show these measure values versus the sample size. The biases, MSEs and ALs decrease to zero when $n$ increases. Further, the CPs are very close to 0.95 when $n$ increases. These results confirm the consistency of the MLEs.

6. Global influence and residual analysis

We present diagnostic measures (Cook & Weisberg, 1982) and residual analysis (Venables & Ripley, 2013) to investigate if the model provides a good representation of the data, i.e., if the sample contains outliers or influential observations.

We employ some measures based on case deletion in the systematic component Eq. (16) to find influential observations in the regression model

$\displaystyle g(\eta_{s})=\exp(\bm{x_{s}}^{\text{T}}\beta),\quad s=1,\ldots,n,% \quad s\neq i.$ (18)

Here, the changes in the parameter estimates coming from exclusion of the $i$ th original observation, where $(i)$ is written as subscript. Then, $l_{(i)}(\eta)$ denotes the log-likelihood function for $\eta$ from model Eq. (18), where the MLE is given by $\hat{\eta}_{(i)}=(\hat{\beta}_{1(i)}^{\text{T}},\ldots,\hat{\beta}_{p(i)}^{% \text{T}})$ .

By comparing the difference between $\hat{\eta}_{(i)}$ and $\hat{\eta}$ , we can measure if, when deleting the $i$ th observation, it influences the parameter estimates. So, this observation can be considered as influential.

The first global influence measure is the Generalized Cook distance (GCD), namely

$\displaystyle\textit{GCD}_{i}=(\hat{\eta}_{(i)}-\hat{\eta})^{\text{T}}[\bm{% \ddot{J}(\hat{\eta})}](\hat{\eta}_{(i)}-\hat{\eta}),$ (19)

where $\bm{J(\hat{\eta})}$ is the observed information matrix.

The second influence measure is the likelihood distance

$\displaystyle\textit{LD}_{i}=2\left\{l(\hat{\eta})-l\left(\hat{\eta}_{(i)}% \right)\right\}.$ (20)

Moreover, the analysis of residuals is an efficient way to check the model adequacy and verify if there are incompatibilities to the response distribution. The deviance residuals for the GM regression are defined by

$\displaystyle r_{D_{i}}=\text{sgn}\left(\hat{r}_{M_{i}}\right)\sqrt{-2[\hat{r}% _{M_{i}}+\log(1-\hat{r}_{M_{i}})]},$ (21)

where

$\displaystyle\hat{r}_{M_{i}}=1+\log\left\{1-\gamma_{1}\left[\alpha,-\log\left(% 1-\gamma_{1}\left(\frac{3}{2},\frac{y_{i}^{2}}{\theta_{i}^{2}}\right)\right)% \right]\right\},$ (22)

are the martingale residuals, $\text{sgn}(\cdot)$ takes value $+$ 1, when the argument is positive, and $-$ 1, otherwise.

The normal probability plot serves to assess the normality assumption of the residuals. The empirical distribution of the deviance residuals agrees with the standard normal distribution when the sample size increases.

We construct envelopes to interpret the probability normal plot of the deviance residuals by means of simulated confidence bands (Atkinson, 1985). When most of the points lie randomly distributed into these bands, the regression model can be considered well-fitted to the data.

7. Application: COVID-19 mortality rates of the 50 largest cities

For showing the utility of the new regression model compared to other competitive models, we provide an application for COVID-19 mortality rates in the USA. These rates are calculated for the 50 biggest cities from the Center for Disease Control and Prevention (CDC) (CDC, 2022).

The dependent variable and explanatory variables obtained from County Health Rankings Model (2022a), U.S. Department of Agriculture (2022) and County Health Rankings Model (2022b) are described below:

1.
MR: Mortality rate (dependent variable).
2.
DB: Diabetes prevalence (percentage of adults with diabetes aged 20) (data from 2019).
3.
UN: Unemployment rate (percentage of average annual unemployment) (data from 2021).
4.
LE: Life expectancy (data from 2018–2020).

Table 2
Summary of statistics

Statistics

Variable Mean Median Skewness Kurtosis Min. Max CV (%)

MR 0.01 0.01 0.49 2.86 0.00 0.02 34.09

DB 0.09 0.09 $-$ 0.16 2.37 0.05 0.13 20.69

UN 0.04 0.03 1.74 6.72 0.02 0.08 30.15

LE 79.80 79.55 0.12 2.48 75.50 85.90 3.00

Table 2 provides some descriptive statistics from this dataset.

We compare the GM model with other alternative generators by taking the Maxwell and Weibull baselines: Maxwell, alpha-power-Maxwell (Erdogan et al., 2021), transmuted-Maxwell (Iriarte & Astorga, 2014), odd log-logistic Maxwell (Prataviera et al., 2020), Weibull, beta-Maxwell (Amusan, 2010) and Maxwell-Weibull (Ishaq & Abiodun, 2020).

The alpha-power-Maxwell (APM) and odd log-logistic Maxwell (OLLM) densities are (for $x>0$ )

$\displaystyle f_{\text{APM}}(x;\theta,\alpha)=\begin{cases}\frac{\log\alpha}{% \alpha-1}\,\frac{4}{\Gamma(1/2)}\,\left(\frac{x}{\theta}\right)^{2}\,{\rm e}^{% -\frac{x^{2}}{\theta^{2}}}\,\alpha^{\frac{1}{\Gamma(3/2)}\gamma\left[\left(% \frac{x}{\theta}\right)^{2},\frac{3}{2}\right]}&\alpha\neq 1\\ f_{\text{Maxwell}}(x;\theta)&\alpha=1\end{cases}$ (23)

and

$\displaystyle f_{\text{OLLM}}(x;\theta,\alpha)=\frac{4\alpha}{\sqrt{\pi}\theta% ^{3}}\,x^{2}\,{\rm e}^{-\frac{x^{2}}{\theta^{2}}}\,\frac{\gamma_{1}^{\alpha-1}% (3/2,x^{2}/\theta^{2})[1-\gamma_{1}(3/2,x^{2}/\theta^{2})]^{\alpha-1}}{\{% \gamma_{1}^{\alpha}(3/2,x^{2}/\theta^{2})+[1-\gamma_{1}(3/2,x^{2}/\theta^{2})]% ^{\alpha}\}^{2}},$ (24)

respectively, where all parameters are positive.

The transmuted-Maxwell (TM) density is (for $x>0$ )

$\displaystyle f_{\text{TM}}(x;\theta,\lambda)=\frac{4}{\sqrt{\pi}}\,\frac{x^{2% }}{\theta^{3}}\,{\rm e}^{-\frac{x^{2}}{\theta^{2}}}\left[1+\lambda-2\lambda% \gamma_{1}(3/2,x^{2}/\theta^{2})\right],$ (25)

for $\alpha>0$ and $|\lambda|\leqslant 1$ .

Finally, the densities of the beta-Maxwell (BM) and Maxwell-Weibull (MW) are (for $x>0$ )

$\displaystyle f_{\text{BM}}(x;\theta,a,b)=\frac{1}{B(a,b)}\,\frac{4}{\sqrt{\pi% }}\,\frac{x^{2}}{\theta^{3}}\,\left[\gamma_{1}(3/2,x^{2}/\theta^{2})\right]^{a% -1}\left[1-\gamma_{1}(3/2,x^{2}/\theta^{2})\right]^{b-1}$ (26)

and

$\displaystyle f_{\text{MW}}(x;\theta,a,b)=\frac{1}{\theta^{3}}\,\left[1-\gamma% _{1}(3/2,x^{2}/\theta^{2})\right]^{-2}\,\frac{4}{\sqrt{\pi}}\,\frac{x^{2}}{% \theta^{3}}\,\left[\frac{\gamma_{1}(3/2,x^{2}/\theta^{2})}{1-\gamma_{1}(3/2,x^% {2}/\theta^{2})}\right]^{2}\exp\left[-\frac{1}{\theta^{2}}\left(\frac{\gamma_{% 1}(3/2,x^{2}/\theta^{2})}{1-\gamma_{1}(3/2,x^{2}/\theta^{2})}\right)^{2}\right]$ (27)

respectively, and all parameters are positive.

The MLEs for all models are calculated using the goodness.fit function from AdequacyModel package (Marinho et al., 2019) available in R Core Team, (2013) with the BFGS method. The best fitted model is chosen based on well-known measures (some with only initials): AIC, CAIC, BIC, HQIC, Cramér-von Mises ( $\text{W}^{}$ ), Anderson-Darling ( $\text{A}^{}$ ) and Kolmogorov-Smirnov (KS) (with its $p$ -value in parentheses).

Table 3
Findings from the fitted models to COVID-19 mortality rates

Model Parameters AIC CAIC BIC HQIC $\text{W}^{}$ $\text{A}^{}$ KS

Maxwell( $\theta$ ) 0.01 $-$ 416.34 $-$ 416.26 $-$ 414.43 $-$ 415.61 0.06 0.35 0.14

(0.00) (0.25)

GM( $\bm{\theta,\alpha}$ ) 0.01 1.64 ${-}$ 419.56 ${-}$ 419.31 ${-}$ 415.74 ${-}$ 418.11 0.06 0.34 0.10

(0.00) (0.27) (0.74)

APM( $\theta,\alpha$ ) 0.01 7.30 $-$ 417.80 $-$ 417.54 $-$ 413.97 $-$ 416.34 0.07 0.45 0.10

(0.00) (6.02) (0.71)

TM( $\theta,\lambda$ ) 0.01 $-$ 0.91 $-$ 419.12 $-$ 418.86 $-$ 415.29 $-$ 417.66 0.06 0.35 0.11

(0.00) (0.26) (0.59)

OLLM( $\theta,\alpha$ ) 0.01 1.33 $-$ 419.25 $-$ 418.99 $-$ 415.42 $-$ 417.79 0.06 0.35 0.10

(0.00) (0.16) (0.74)

Weibull( $a, b$ ) 3.20 0.01 $-$ 417.66 $-$ 417.40 $-$ 413.84 $-$ 416.20 0.08 0.49 0.11

(0.34) (0.00) (0.63)

BM( $\theta,a,b$ ) 0.012 1.47 2.86 $-$ 417.30 $-$ 416.78 $-$ 411.57 $-$ 414.52 0.22 1.40 0.83

(0.00) (0.28) (1.86) (0.00)

MW( $\theta,a,b$ ) 0.00 1.26 0.47 $-$ 416.79 $-$ 416.28 $-$ 411.05 $-$ 414.61 0.07 0.40 0.10

(0.00) (0.13) (0.18) (0.73)

The results for all statistics are reported in Table 3 with the MLEs and standard errors (SEs) in parentheses. They reveal that the GM distribution is the best fitted model. Further, plots of the estimated pdfs and cdf for the three best models in Fig. 5 confirm that the GM distribution is the most suitable model for the current COVID-19 mortality rate data.

Table 4
Fitted GM regression to COVID-19 in the 50 biggest cities in USA

Parameters Estimates SEs $p$ -value

Intercept $-$ 12.19 1.89 0.00

DB 10.92 2.57 0.00

UN 16.31 4.72 0.00

LE 0.07 0.02 0.00

$\alpha$ 2.41 0.51 0.00

Figure 5.
(a) Histogram and estimated pdfs. (b) Empirical and estimated cdfs.

The LR statistic for testing $H_{0}:\alpha=$ 1 against $H_{1}:H_{0}$ is false, i.e., the GM versus Maxwell distribution is $w=$ 5.22 with $p$ -value 0.0223. It then provides clear evidence of the significance of the extra parameter to model the current data.

Figure 6.
Profile log-likelihood functions from the fitted GM regression model to COVID-19 data. Parameters: (a) $\hat{\alpha}$ . (b) $\hat{\beta_{0}}$ . (c) $\hat{\beta_{1}}$ . (d) $\hat{\beta_{2}}$ . (e) $\hat{\beta_{3}}$ .

Figure 7.
(a) GCD for the fitted GM regression model. (b) LD for the fitted GM regression model.

Nest, we consider the systematic component (for $i=1,\ldots,50$ ):

$\displaystyle\theta_{i}=\exp(\text{Intercept}+x_{i1}\,\text{DB}+x_{i2}\,\text{% UN}+x_{i3}\,\text{LE}).$ (28)

Table 4 lists the results from the fit of the previous GM regression model. For a significance level of 5%, all explanatory variables are significant.

Profile log-likelihood plots versus some parameter values (with all other parameter estimates fixed) are reported in Fig. 6. These plots can provide confidence intervals for the parameters.

Figure 8.
Normal probability plot of $r_{D}$ ’s with simulated envelope.

The GCD and the LD measures identify possible influential observations, nevertheless these do not affect the model, as illustrated in Fig. 7. The plot of the deviance residual ( $r_{D}$ ’s) with simulated envelope in Fig. 8 reveal that they lie inside the bands randomly, thus indicating a good fit of the GM regression model.

Finally, from the parameter estimates reported in Table 4, the GM regression model becomes

$\displaystyle\hat{\theta}_{i}=\exp(-12.19+10.92\,\text{DB}+16.31\,\text{UN}+0.% 07\,\text{LE}).$ (29)

Some conclusions follow from the above equation:

•
The DB prevalence is extremely significant ( $p$ -value 0.00). Corona et al. (2021) provided an overview and a meta-analysis of main predictors, and, as expected of another research, diabetes was the best predictor of mortality rate. According to study (Stokes et al., 2020), counties with high levels of diabetes had excess mortality not assigned for COVID-19.
•
The UN rate gives a $p$ -value of 0.00, and its estimate is positive. So, the MR increases in cities with a larger unemployment rate, which corroborates with the research in Paul et al., (2021) and Mirahmadizadeh et al. (2022). The first showed increased mortality in high unemployment urban cities, while the second, with a cohort study, examined unemployment as a risk factor for COVID-19 mortality rate in Iran.
•
The LE is significant at the 1% level, and its estimate is positive, which indicates that the MR is slightly higher in cities with high life expectancy. Similarly, Notari and Torrieri (2022) showed that in the European region, the life expectancy is an increase factor for COVID-19 mortality. Further, Wang et al. (2020) also determined a positive correlated between life expectancy and the initial transmission growth rate of COVID-19.

8. Conclusions

Model	Parameters	AIC	CAIC	BIC	HQIC	$\text{W}^{*}$	$\text{A}^{*}$	KS
Maxwell( $\theta$ )	0.01			$-$ 416.34	$-$ 416.26	$-$ 414.43	$-$ 415.61	0.06	0.35	0.14
	(0.00)									(0.25)
GM( $\bm{\theta,\alpha}$ )	0.01	1.64		${-}$ 419.56	${-}$ 419.31	${-}$ 415.74	${-}$ 418.11	0.06	0.34	0.10
	(0.00)	(0.27)								(0.74)
APM( $\theta,\alpha$ )	0.01	7.30		$-$ 417.80	$-$ 417.54	$-$ 413.97	$-$ 416.34	0.07	0.45	0.10
	(0.00)	(6.02)								(0.71)
TM( $\theta,\lambda$ )	0.01	$-$ 0.91		$-$ 419.12	$-$ 418.86	$-$ 415.29	$-$ 417.66	0.06	0.35	0.11
	(0.00)	(0.26)								(0.59)
OLLM( $\theta,\alpha$ )	0.01	1.33		$-$ 419.25	$-$ 418.99	$-$ 415.42	$-$ 417.79	0.06	0.35	0.10
	(0.00)	(0.16)								(0.74)
Weibull( $a, b$ )	3.20	0.01		$-$ 417.66	$-$ 417.40	$-$ 413.84	$-$ 416.20	0.08	0.49	0.11
	(0.34)	(0.00)								(0.63)
BM( $\theta,a,b$ )	0.012	1.47	2.86	$-$ 417.30	$-$ 416.78	$-$ 411.57	$-$ 414.52	0.22	1.40	0.83
	(0.00)	(0.28)	(1.86)							(0.00)
MW( $\theta,a,b$ )	0.00	1.26	0.47	$-$ 416.79	$-$ 416.28	$-$ 411.05	$-$ 414.61	0.07	0.40	0.10
	(0.00)	(0.13)	(0.18)							(0.73)

Parameters	Estimates	SEs
Intercept	$-$ 12.19	1.89
DB	10.92	2.57
UN	16.31	4.72
LE	0.07	0.02
$\alpha$	2.41	0.51

We proposed a new gamma-Maxwell regression model, and obtained some new mathematical properties of the response distribution. The model parameters were estimated by the maximum likelihood method, and the consistency of the estimators was proved using Monte Carlo simulations. Diagnostic analysis and deviance residuals were provided. We proved that the new model is more flexible than some others by analyzing COVID-19 mortality rates in the fifty largest U.S. cities.

Footnotes

https://www.cdc.gov/coronavirus/2019-ncov/symptoms-testing/symptoms.html.

https://covid.cdc.gov/covid-data-tracker/#datatracker-home.

References

Amusan

G. E.

(2010). The beta Maxwell distribution. Theses, Dissertations and Capstones, 5.

Atkinson

A. C.

(1985). Plots, transformations, and regression: An introduction to graphical methods of diagnostic regression analysis.

Biazatti

E. C.

Cordeiro

G. M.

Rodrigues

G. M.

Ortega

E. M.

& de Santana

L. H.

(2022). A Weibull-beta prime distribution to model COVID-19 data with the presence of covariates and censored data. Stats, 5, 1159-1173.

Castellares

& Lemonte

A. J.

(2015). A new generalized Weibull distribution generated by gamma random variables. Journal of the Egyptian Mathematical Society, 23, 382-390.

Center for Diseases Control and Prevention (CDC) (2022). Retrieved April 11, 2022, from https://COVID.cdc.gov/COVID-data-tracker/#county-view?list_select_state=all_states&list_select_county=all_counties.

Cook

R. D.

& Weisberg

(1982). Residuals and influence in regression.

Cordeiro

G. M.

Alizadeh

Ozel

Hosseini

Ortega

E. M. M.

& Altun

(2017). The generalized odd log-logistic family of distributions: Properties, regression models and applications. Journal of Statistical Computation and Simulation, 87, 908-932.

Corona

Pizzocaro

Vena

Rastrelli

Semeraro

Isidori

A. M.

, …& Maggi

(2021). Diabetes is most important cause for mortality in COVID-19 hospitalized patients: Systematic review and meta-analysis. Reviews in Endocrine and Metabolic Disorders, 22, 275-296.

County Health Rankings Model (2022). Health Outcomes – Diabetes Prevalence. Retrieved April 13, 2022, from https://www.countyhealthrankings.org/explore-health-rankings/measures-data-sources/county-health-rankings-model/health-outcomes/quality-of-life/diabetes.

10.

County Health Rankings Model (2022). Health Outcomes – Life Expectancy. Retrieved April 13, 2022, from https://www.countyhealthrankings.org/explore-health-rankings/measures-data-sources/county-health-rankings-model/health-outcomes/length-of-life/life-expectancy.

11.

Duhon

Bragazzi

& Kong

J. D.

(2021). The impact of non-pharmaceutical interventions, demographic, social, and climatic factors on the initial growth rate of COVID-19: A cross-country study. Science of The Total Environment, 760, 144325.

12.

Erdogan

Bagci

Arslan

& Celik

H. E.

(2021). Alpha power Maxwell distribution: Properties and application. Journal of Mathematical Modeling, 9, 585-598.

13.

Gallardo

D. I.

Bourguignon

Gómez

Y. M.

Caamaño-Carrillo

& Venegas

(2022). Parametric quantile regression models for fitting double bounded response with application to COVID-19 mortality rate data. Mathematics, 10, 2249.

14.

Gradshteyn

I. S.

& Ryzhik

I. M.

(2014). Table of integrals, series, and products. Academic press.

15.

Iriarte

Y. A.

& Astorga

J. M.

(2014). Transmuted Maxwell probability distribution. Revista IntegraciÃn, 32, 211-221.

16.

Iriarte

Y. A.

Astorga

J. M.

Bolfarine

& Gómez

H. W.

(2017). Gamma-Maxwell distribution. Communications in Statistics-Theory and Methods, 46, 4264-4274.

17.

Ishaq

A. I.

& Abiodun

A. A.

(2020). The Maxwellâ€“Weibull distribution in modeling lifetime datasets. Annals of Data Science, 7, 639-662.

18.

Janke

A. T.

Mei

Rothenberg

Becher

R. D.

Lin

& Venkatesh

A. K.

(2021). Analysis of hospital resource availability and COVIDâ€19 mortality across the United States. Journal of Hospital Medicine, 16, 211-214.

19.

Karmakar

Lantz

P. M.

& Tipirneni

(2021). Association of social and demographic factors with COVID-19 incidence and death rates in the US. JAMA Network Open, 4, e2036462-e2036462.

20.

Marinho

P. R. D.

Silva

R. B.

Bourguignon

Cordeiro

G. M.

& Nadarajah

(2019). AdequacyModel: An R package for probability distributions and general purpose optimization. PloS One, 14, e0221487.

21.

Mirahmadizadeh

Shamooshaki

M. T. B.

Dadvar

Moradian

M. J.

& Aryaie

(2022). Unemployment and COVID-19-related mortality: A historical cohort study of 50,000 COVID-19 patients in Fars, Iran. Epidemiology and Health, 44.

22.

Notari

& Torrieri

(2022). COVID-19 transmission risk factors. Pathogens and Global Health, 116, 146-177.

23.

Paul

Arif

Pokhrel

& Ghosh

(2021). The association of social determinants of health with COVIDâ€19 mortality in rural and urban counties. The Journal of Rural Health, 37, 278-286.

24.

Prataviera

Ortega

E. M.

& Cordeiro

G. M.

(2020). A new bimodal Maxwell regression model with engineering applications. Applied Mathematics & Information Sciences, 14, 817-831.

25.

Prudnikov

A. P.

& Marichev

O. I.

(1986). Integrals and series: special functions (Vol. 2). CRC press.

26.

R Core Team (2013). R: A language and environment for statistical computing.

27.

Sousa

G. J. B.

Garces

T. S.

Cestari

V. R. F.

Florêncio

R. S.

Moreira

T. M. M.

& Pereira

M. L. D.

(2020). Mortality and survival of COVID-19. Epidemiology & Infection, 148.

28.

Stokes

Lundberg

D. J.

Hempstead

Elo

I. T.

& Preston

S. H.

(2020). Assessing the impact of the COVID-19 pandemic on US mortality: A county-level analysis. MedRxiv.

29.

U.S. Department of Agriculture (2022). County-level data sets. Retrieved April 15, 2020, from https://www.ers.usda.gov/data-products/county-level-data-sets/unemployment/.

30.

Venables

W. N.

& Ripley

B. D.

(2013). Modern applied statistics with S-PLUS. Springer Science & Business Media.

31.

Wang

X. Q.

Song

Yang

Chen

R. J.

Zheng

Y. L.

H. Y.

, …& Chen

P. J.

(2020). Association between ageing population, median age, life expectancy and mortality in coronavirus disease (COVID-19). Aging (Albany NY), 12, 24570.

	Statistics
Variable	Mean	Median	Skewness	Kurtosis	Min.	Max	CV (%)
MR	0.01	0.01	0.49	2.86	0.00	0.02	34.09
DB	0.09	0.09	$-$ 0.16	2.37	0.05	0.13	20.69
UN	0.04	0.03	1.74	6.72	0.02	0.08	30.15
LE	79.80	79.55	0.12	2.48	75.50	85.90	3.00