The Dual-Dagum family of distributions: Properties, regression and applications to COVID-19 data

1.Introduction

Over the last few decades, many generators have been studied in the distribution theory literature. Two generators that stand out are the Beta-G (B-G) (Eugene et al., 2002) and Kumaraswamy-G (Kw-G) (Cordeiro & de Castro, 2011) classes.

Regarding the B-G family, we can say that, although it contains the incomplete and complete beta functions, its flexibility in terms of adjustment to real data is widespread. Several authors introduced new distributions in this family in different contexts: cancer recurrence (Paranaíba et al., 2011), waiting times before service of 100 bank customers (Abd El-Bar & Ragab, 2015), test on the endurance of deep groove ball bearings (Abu-Zinadah & Bakoban, 2017), survival times of 33 patients suffering from acute Myelogeneous Leukaemia (Mead et al., 2017), among others. More than one-hundred different published distributions in this class can be found to date.

The second family stands out because of the simplicity of its density function, which does not include complicated functions. Further, its suitability for the most diverse types of data sets is widely discussed in the literature. We can cite, for example, the work that originated this family and used data from adult numbers of T. confusum cultured at 29∘C (Eugene et al., 2002). Since then, many sub-models in this family have been proposed for various types of data sets: fatigue (Cordeiro et al., 2010), body mass index (Mameli, 2015), component lifetime (Cordeiro et al., 2016), among others.

We know through an analysis of these works that the fits of both classes to real data have a better performance compared to other known classes. We can note that the data sets studied in the aforementioned works are of different types. However, many authors end up repeating the same data sets used in previous works by other authors.

In this sense, we define a new class from the Dagum distribution (Dagum, 1975) and use data bases never published before. The data bases in question concern a very current topic: COVID-19. We understand the importance of studies on this pandemic that impacted the world, and then use COVID-19 data from two cities in Brazil.

The remainder of the paper is organized as follows. Section 2 defines the new family. In Section 3, we present some of its generated distributions. The main properties of the new family are reported in Section 4. Estimation including the case of censoring is addressed in Section 5. A simulation study is done in Section 6. In Section 7, we construct the Log-Dual-Dagum-Weibull regression, and estimate the parameters. Two applications to real data are reported in Section 8, including a regression application and a sensitivity analysis. Conclusions end the paper in Section 9.

2.The new family

The new generator is defined based on the survival function of the Dagum distribution (Dagum, 1977). Kleiber and Kotz (2003) and Kleiber (2008) analyzed characteristics and properties of this distribution. The Dagum distribution presents forms of the increasing, decreasing, bathtub and inverted bathtub risk function (Domma, 2002). This behavior has aroused the interest of several authors to study it in survival analysis (Domma et al., 2011a, b). In this sense, we propose the Dual-Dagum-G (DDa-G) family.

Let W be a Dagum random variable with two positive shape parameters a and b. The cumulative distribution function (cdf) of the DDa-G family (for x∈ℝ) is

where G⁢(x)=G⁢(x;𝜽) is the baseline cdf, and 𝜽 is its parameter vector.

Henceforth, Eq. (1) refers to the random variable X∼DDa-G⁢(a,b).

The probability density function (pdf) of X has the form

where g⁢(x)=∂⁡G⁢(x)/∂⁡x.

Equations (1) and (2) do not involve complicated mathematical functions, which is an advantage of this family when compared, for example, with the Beta generator.

The hazard rate function (hrf) of X is

3.Special DDa-G distributions

3.2Dual-Dagum-log-logistic (DDa-LL)

The cdf of the log-logistic (LL) distribution is (for x,α,β>0)

G⁢(x)=1-[1+(xα)β]-1.

Inserting this expression and its derivative in Eq. (2) leads to the DDa-LL density (for x>0)

fDDa-LL⁢(x)=a⁢b⁢β⁢xβ-1⁢[1+(xα)]-2⁢{1+[-log⁡(1-[1+(xα)β]-1)]-a}-b-1αβ⁢{1-[1+(xα)β]-1}⁢[-log⁡(1-[1+(xα)β]-1)]a+1.

Figure 1.

Shapes of the (a) DDa-W⁢(a,b,k,λ) and (b) DDa-LL⁢(a,b,β,α) densities.

Figure 2.

Shapes of the (a) DDa-W⁢(a,b,k,λ) and (b) DDa-LL⁢(a,b,β,α) hazard rates.

Figures 1 and 2 display shapes of the pdf and hrf of the previous generated models, which show their flexibility in fitting data with different shapes. For example, the Weibull pdf presents only decreasing and unimodal shapes, whereas the DDaW pdf has an extra shape: decreasing-increasing-decreasing.

4.Some properties

4.2Quantile function

Let QG⁢(u)=G-1⁢(u) be the quantile function (qf) of G (for 0<u<1). Inverting F⁢(x)=u in Eq. (1), the qf of X becomes

(7)

x=Q⁢(u)=F-1⁢(u)=QG⁢(e-[-1+(1-u)-1/b]-1/a).

Equation (7) reveals that the qf of the proposed family is a function of the baseline qf.

The skewness and kurtosis of X can be calculated from the first four ordinary moments (see Section 4.3). Alternatively, approximations for the skewness (S) and kurtosis (K) can be obtained from Eq. (7) and the formulae

S=Q⁢(3/4)-2⁢Q⁢(1/2)+Q⁢(1/4)Q⁢(3/4)-Q⁢(1/4)

and

K=Q⁢(7/8)-Q⁢(5/8)+Q⁢(3/8)-Q⁢(1/8)Q⁢(6/8)-Q⁢(2/8),

reported by Kenney and Keeping (1961) and Moors (1988), respectively.

Figures 3 and 4 display the skewness and kurtosis of the DDa-W distribution as functions of both a and b.

Figure 3.

Bowley’s skewness of the DDa-W distribution. (a) as function of a (b=λ=2) and (b) as function of b (a=λ=2).

Figure 4.

Moors’ kurtosis of the DDa-W distribution. (a) as function of a (b=λ=2) and (b) as function of b (a=λ=2).

5.Estimation

The estimation of the unknown parameters of the DDa-G distribution is performed by the maximum likelihood method. Let x1,…,xn be a sample from Eq. (2). Let 𝜽 be the q×1 parameter vector in G⁢(⋅). The log-likelihood function log⁡L=log⁡L⁢(a,b,𝜽) has the form

The R software has the AdequacyModel computational library (Marinho et al., 2019) as a good alternative for maximizing log⁡L⁢(a,b,𝜽), obtain the maximum likelihood estimates (MLEs), their standard errors, and some statistical measures to evaluate the adequacy of a fitted distribution.

6.Simulation study

We adopt the exponential (E) baseline (with the expected value λ) for the simulations to assess the accuracy of the MLEs in the DDa-E(a,b,λ) model. In order to generate observations from the DDa-E distribution, we adopt the inversion method following the steps:

We consider 2,000 Monte Carlo replications and the BFGS algorithm in the R software for maximizing the log-likelihood, obtain the MLEs and their averages, biases and mean square errors (MSEs). The simulation process is carried out as below:

Table 1 reports these findings. The average estimates converge to the true parameter values and the biases decrease when n increases (for all scenarios).

7.Log-Dual-Dagum-Weibull (LDDa-W) regression

If X follows density Eq. (3), then Y=log⁡(X) has the LDDa-W distribution, which reparameterized in terms of k=σ-1 and λ=eμ, has the form (for y∈ℝ)

where μ∈ℝ, σ>0, a>0 and b>0.

The survival function corresponding to Eq. (8) is

The pdf of the standardized random variable Z=(Y-μ)/σ (for z∈ℝ) is

The lifetimes xi are affected by known explanatory variables 𝐯i=(vi⁢1,…,vi⁢p)⊤ in many applications. Consider a sample (y1,𝐯1),…,(yn,𝐯n) of independent observations, where yi=min⁡{log⁡(xi),log⁡(ci)}, and the log-lifetime log⁡(xi) and the log-censoring log⁡(ci) are assumed independent (under non-informative censoring).

The LDDa-W regression for the response variable yi is defined by

where zi is the random error with density Eq. (9), γ=(γ1,…,γp)⊤, σ>0 is a scale, and a>0 and b>0 are shape parameters. The location vector 𝝁=(μ1,…,μn)⊤, where μi=𝐯i⊤⁢γ, is 𝝁=𝐕⁢γ, where 𝐕=(𝐯1,…,𝐯n)⊤ is a known model matrix. The log-DDa-exponential (LDDa-E) regression is given by Eq. (10) with σ=1.

Let F and C be the sets representing the observed lifetimes and censoring times, respectively. The total log-likelihood function for 𝜽=(a,b,σ,γ⊤)⊤ can be determined from Eqs (9) and (10) as

where q is the observed number of failures. The MLE 𝜽^ of 𝜽 can be found by maximizing Eq. (11).

8.Applications

The Weibull and Birnbaum-Saunders distributions are taken as baselines to prove the flexibility of the new family. The data sets were obtained from the open data portal of the Federal Government linked to the Ministry of Health and comprise events from 2020–2021 (accessed on August 23, 2021). The data portal is available at https://dados.gov.br/dataset/bd-srag-2020.

All computations are done in R using

GenSA

MASS

AdequacyModel

goodness.fit()

SANN

MASS

The new distributions are compared with well-known models belonging to the Kw-G and B-G classes using the statistics: Cramér-von Mises (W*), Anderson-Darling (A*), Kolmogorov- Smirnov (K⁢S), and p-values of the K⁢S test. We present below the alternative densities for the applications.

In the following, we calculate descriptive statistics, MLEs, their standard errors (SEs) and adequacy statistics to compare the fitted distributions to the data sets.

8.2Regression modeling applied to COVID-19 data in Curitiba (Brazil)

The study comprises the time (in days) elapsed from the date of hospitalization until death by the coronavirus, of 1,951 patients in Curitiba-PR, with all observations failing, that is, censored times were not considered in the study, with occurrences of death in 2020 and 2021.

The explanatory variables are (for i=1,…,1951):

• • vi⁢1: sex (1 = masculine, 2 = feminine);

• • vi⁢2: age (in years).

The computational part is developed in R using

survival

library, with

optim()

function, and the “

SANN

” method.

We adopt the Akaike information criterion (AIC), corrected Akaike information criterion (CAIC), and Bayesian information criterion (BIC) to choose the appropriate model. We compare the fits of the LDDa-W Eq. (8) with the log-Kumaraswamy-Weibull or Kumaraswamy Gumbel (Kw-Gu) (Cordeiro et al., 2012), log-beta Weibull (LBW) and log-Weibull (LW) models. The densities for the alternative regressions are reported below:

Table 3

Estimation results from some fitted regressions to the COVID-19 data in Curitiba, and the adequacy measures

	LDDa-W		Kw-Gu		LBW		LW
Parameter	Estimate	SE	Estimate	SE	Estimate	SE	Estimate	SE
γ0	5.635*	0.419	3.496*	0.311	3.378*	0.508	4.090*	0.100
γ1	-0.122*	0.035	-0.147*	0.037	-0.068**	0.039	-0.193*	0.035
γ2	-0.009*	0.001	-0.013*	0.001	-0.012*	0.001	-0.013*	0.001
σ	2.355	0.349	1.419	0.186	1.973	0.608	0.758	0.012
a	5.876	0.709	3.076	0.795	6.065	3.688	–	–
b	1.677	0.348	1.813	0.234	3.492	1.124	–	–
AIC	4880.25		4882.50		4895.15		4990.87
CAIC	4880.29		4882.54		4895.19		4990.89
BIC	4913.70		4915.95		4928.61		5013.18

p*-value < 0.0001; p**-value < 0.1.

Figure 5.

Estimated DDa-W, DDa-BS, Kw-W and Beta-W densities.

• • The LW (or Gumbel) density function

  f⁢(y;μ,σ)=1σ⁢exp⁡{(y-μσ)-exp⁡(y-μσ)},

  where y∈ℝ, σ>0 and μ∈ℝ; see Johnson et al. (1995).

• • The LBW density function

  f⁢(y;a,b,μ,σ)=1σ⁢B⁢(a,b)⁢exp⁡{(y-μσ)-b⁢exp⁡(y-μσ)}⁢{1-exp⁡[-exp⁡(y-μσ)]}a-1,

  where y∈ℝ, σ>0 and μ∈ℝ. For more details, see Ortega et al. (2013).

• • The Kw-Gu density function

  f⁢(y;a,b,μ,σ)=a⁢bσ⁢exp⁡{(y-μσ)-exp⁡(y-μσ)}⁢{1-exp⁡[-exp⁡(y-μσ)]}a-1×[1-(1-exp⁡[-exp⁡(y-μσ)])a]b-1,

  where y∈ℝ, σ>0 and μ∈ℝ.

The failure rate function is useful to aid in model identification more suitable for the variable time. In this context, the TTT plot (not shown here) for the data under study shows an increasing appearance for the most part, but due to its final behavior, it indicates an inverting bathtub risk function. The descriptive statistics for the time until death for COVID-19 data in Curitiba include: mean = 16.86, SD = 15.54, skewness = 3.75, kurtosis = 29.77, and minimum and maximum values of 1 and 205, respectively.

Next, we provide results from the fit of the regression

(12)

yi=γ0+γ1⁢vi⁢1+γ2⁢vi⁢2+σ⁢zi,

where z1,…,z1951 are independent random variables with density function Eq. (9).

Table 3 provides some findings from the regressions fitted to the current data. They indicate that LDDa-W model provides the best fit to the data. Further, all covariates (vi⁢1: sex and vi⁢2: age) are significant at 1%, and then the categories of these explanatory variables are statistically different.

Thus, the time to death decreases when the age increases. Regarding the patient’s gender, male patients present smaller time until death than female patients, since the estimate of its coefficient is negative.

After the LDDa-W regression estimation, the plots of the empirical and estimated survival functions support the model adequacy to these data.

Also, as part of the analysis, it is important to verify if there are observations influencing the model’s adjustment. A sensitivity analysis was carried out to investigate this fact using the Cook’s distance and will be presented below.

8.2.1Sensitivity analysis and influential observations

Under the Generalized Cook Distance (Cook, 1977), the observations #349, #826 and #897 are the ones that stand out the most, thus indicating that they can be possible influential observations.

The observation #349 refers to a female individual, aged 84 years and with a time of hospitalization until death of 172 days. The observation #826 is identified as a 78-year-old male and has a time to death of 6 days. And the observation #897 refers to the male individual aged 61 years old, whose hospitalization time until death was only 3 days. The observations #349, #826 and #897 represent individuals with peculiar behaviors, but do not show signs of error in data collection or transcription, and therefore must be kept in the database. The final model is given in Eq. (12).

The impact of possible influential observations detected should be analyzed in order to assess the estimates and sensitivity of the model. This analysis considers new estimates for the model parameters from sub-samples referring to the withdrawal of these observations individually and in groups.

It is considered that the changes in the estimated values for the parameters are not very expressive and there is significance of the explanatory variables when considering the level of 10%. In addition, there was no change in the sign of the coefficient of the explanatory variables, so the inclusion or exclusion of the identified observations does not presuppose changes in the interpretation of the results.

9.Conclusions

One of the main objectives of distribution theory is to define a family of models to better explain lifetime phenomenon in several areas of knowledge. We proposed the Dual-Dagum-G (“DDa-G”) family, which can generalize all classical continuous distributions. Its parameters are estimated by maximum likelihood, and a simulation study showed the consistency of the estimators. We showed the flexibility of the new family by means of two real COVID-19 data sets. We proved that the new Log-Dual-Dagum-Weibull (LDDaW) regression outperformed regressions based on well-known Kumaraswamy-G and Beta-G generators. After verifying the good fit of the new regression, a sensitivity analysis was performed, where it was possible to verify the occurrence of influential observations. As future work, it could be interesting to investigate other methods of sensitivity analysis, such as the local influence, if the results obtained through the Cook Distance prevail and still carry out a residual analysis for the new regression model.

Note: Computer codes in

R