Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
45,778 characters · 13 sections · 24 citation commands
Parametric quantile regression for income data
{ { Keywords.} {Income distributions $\cdot$ Quantile regression $\cdot$ Income data $\cdot$ Reparameterization.}}
\onehalfspacing
Income modeling plays an important role in determining workers' earnings, as well as being an important research topic in labor economics. In general, income data are modeled using mean-based regression models based on the normality assumption. Nevertheless, income is often unequally distributed, hence why this type of data usually has an asymmetric behavior and then the mean is not an appropriate central tendency measure. Therefore, quantile regression models are usually more useful in this context; see galarza2020, sanchez2021a and sauloetal:21.
Quantile regression models are robust alternatives to traditional mean-based models. That is because instead of focusing on the conditional mean, these models are based on the conditional quantile, such as median; see koenker:05. The quantile approach has the advantage of providing flexibility in modeling, as it allows considering the effects of explanatory variables throughout the spectrum of the dependent variable, thus also including the effect on the median, which is a measure of central tendency better than the mean in the asymmetric context.
Income modelling begins with pareto1897 propositions, establishing a law on how income distribution works. Later on, this suggested a distribution -- known as Pareto distribution -- and it has set a reference for other distributions, such as log-normal and gamma, to show their potential as for describing income distribution; see shirras1935, reed2003. Even though the Pareto, log-normal and gamma are the most frequently distributions applied to income data because of their abilities to describe this type of data, they have limitations. On the one hand, the Pareto model is appropriate to describe only the upper tail of the distribution. On the other hand, the log-normal and gamma distributions perform poor in describing both the upper and lower tails of the actual distributions. Income distributions such as Dagum and Singh-Maddala have outperformed the Pareto, log-normal and gamma distributions in terms of model fitting; see cramer1971, dagum2008.
Originally proposed by dagum1973,dagum1975, the Dagum distribution has flexibility to deal with distribution changes, nil and negative income, income range with non-predetermined positive minimum income start, and strictly decreasing and unimodal density functions. This distribution also shows good goodness of fit to income data and obeys a weak version of the Pareto law, i.e. it asymptotically approaches the Pareto distribution. The Dagum model accommodates both heavy tails and an interior mode, characteristics commonly found in income data, and not found singly in well-known distributions -- such as log-normal and Pareto; see kramer2002,dagum2008,kleiber2008.
The Singh-Maddala distribution was derived from the concept of hazard rate, an approach widely used in the reliability literature; see singh1976. This model also obeys the weak Pareto law, and one of its advantages is to be more flexible than other income distributions. The Dagum and Singh-Maddala distributions are special cases of the generalized beta distribution of the second kind (GB2); for more details on these models, one may refer to the works by kleiber1996,kleiber2008,kumar2017,hajargasht2012.
This work explores a parametric quantile regression approach for the Dagum and Singh-Maddala distributions. We first introduce reparameterizations of the Dagum and Singh-Maddala model by inserting quantile parameters, and then develop the new regression models. We then demonstrate that the proposed models outperform the recently proposed Birnbaum-Saunders quantile regression model sanchez2021b in terms of model fitting.
The rest of this paper proceeds as follows. In Section (ref), we describe the usual Dagum and Singh-Maddala distributions and propose reparameterizations of these distributions in terms of a quantile parameter. In this section, we also present some properties including mode, real and truncated moments. In Section (ref), we introduce the quantile regression models and also describe the parameter estimation by the maximum likelihood (ML) method. In Section (ref), we carry out a Monte Carlo simulation study to evaluate the performance of the estimators and generalized Cox-Snell (GCS) and random quantile (RQ) residuals. In Section (ref), we apply the Dagum and Singh-Maddala quantile regression models to a household income data set provided by the National Institute of Statistics of Chile, and finally in Section (ref), we provide some concluding remarks.
In this section, we describe the classical Singh-Maddala and Dagum distributions along with the proposed quantile-based reparameterizations of these distributions, which will be useful subsequently for developing the parametric quantile regression models. We also present some properties for each model, including mode, real and truncated moments.
If a random variable $Y$ follows a Singh-Maddala distribution with shape parameters $a,q>0$ and scale parameter $b>0$, denoted by $Y\sim \text{SM}(a, b, q)$, then the corresponding probability density function (PDF) and cumulative distribution function (CDF) are given by
and
respectively. The Singh-Maddala distribution includes as special cases the Lomax distribution when $a = 1$, and the log-logistic distribution when $q = 1$. If $Y$ follows a Singh-Maddala distribution, then $1/Y$ follows a Dagum distribution, and vice-versa.
The $\tau$-th quantile of $Y\sim \text{SM}(a, b, q)$ is obtained by inverting Equation (ref), which yields
From the quantile function (ref), we find that the most parsimonious way of conducting the reparametrization is using the scale parameter $b$, where we can then write
where $\gamma = q(\tau; a, b, q) > 0$. Then, the quantile-based Singh-Maddala PDF is given by
with notation $Y\sim {\rm QSM}(a,\gamma,q)$.
If $Y\sim {\rm QSM}(a,\gamma,q)$, then the following properties hold:
The PDF and CDF of a random variable $Y$ following a classical Dagum distribution with shape parameters $a,p>0$ and scale parameter $b>0$, denoted by $Y\sim \text{DA}(a, b, p)$, are given by
and
respectively. It is simple to observe that $ f_{\rm DA}(y; a, b, p) = (y/b)^{a(p-1)} f_{\rm SM}(y; a, b, p) $ and that when $p=1$ both densities coincide with the log-logistic distribution.
The $\tau$-th quantile of $Y\sim \text{DA}(a, b, p)$ is given by
By observing the three parameters of the classical Dagum distribution, the isolation of the scale according to the quantile would produce the simplest form of the new quantile-based Dagum distribution; it is represented as follows:
where $\gamma = q(\tau; a, b, p) > 0$. Then, the quantile-based Dagum PDF can be written as
If $Y\sim {\rm QDA}(a,\gamma,p)$, then the following properties hold:
Table (ref) presents the Singh-Maddala and Dagum distributions in their original and quantile-based versions. Figures (ref) and (ref) display different shapes of the quantile-based income distributions for different combinations of parameters, considering scenarios where $a$, $p$, $q$ and $\gamma$ are fixed. For Singh-Maddala, we can see that $a$ influences the kurtosis and skewness, while $q$ changes the kurtosis, as it decreases when $q$ increases. For Dagum, we see a similar pattern for $a$, changing both kurtosis and skewness, while $p$ affects the kurtosis.
Let $ Y_1, \ldots, Y_n $ be independent random variables such that each $ Y_i $, for $ i = 1,\ldots,n $, has PDF given by some reparameterized income distribution defined in Table (ref), for a fixed (known) probability $ \tau \in (0, 1)$ associated with the quantile of interest. Then, in the formulation of the Singh-Maddala and Dagum quantile regression models, the parameter $\gamma$ of $Y_i$ assumes the following functional relation:
where $\bm{\beta}(\tau) = (\beta_{0}(\tau), \ldots, \beta_{k}(\tau))^\top$ is the vector of the unknown regression coefficients, which are assumed to be functionally independent; $\bm{\beta}(\tau) \in \mathbb{R}^{(k+1)}$, with $k +1 < n$; and $\mathbf{x}_{i} = (x_{i1}, \ldots, x_{il})^\top$ is the observations of the $l$ known regressors, for $i = 1, \ldots, n$. In addition, we assume that the covariate matrices $\mathbf{X} = (\mathbf{x}_1, \ldots, \mathbf{x}_n)^\top$ has rank $l$. The link function $g: \mathbb{R}^+ \rightarrow \mathbb{R}$ in (ref) must be strictly monotone, positive, and at least twice differentiable, with $g^{-1}(\cdot)$ being the inverse function of $g(\cdot)$. Here, we chose to work with $\log$ as link since it is widely used and more flexible when it comes to simulation studies.
Consider a sample of size $n$, $Y_1,\ldots,Y_n$ say, such that $Y_i \sim \text{QSM}(a, \gamma_i, q)$. Then, the corresponding likelihood function for $\bm{\theta} = (\bm{\beta}(\tau)^{\top},a,q)^{\top}$, is
where $\gamma_i$ is as in (ref). By applying the logarithm in (ref), we obtain the log-likelihood function
Now, consider a sample of size $n$, $Y_1,\ldots,Y_n$ say, such that $Y_i \sim \text{QDA}(a, \gamma_i, p)$. Then, the corresponding likelihood function for $\bm{\theta} = (\bm{\beta}(\tau)^{\top},a,p)^{\top}$, is
where $\gamma_i$ is as in (ref). By applying the logarithm in (ref), we obtain the log-likelihood function
To obtain the ML estimate of $\bm{\theta}$, it is necessary to maximize the log-likelihood functions in (ref) and (ref). Therefore, we need to differentiate the log-likelihood functions to find the score vector $\dot{\ell}(\bm{\theta})$ and then equate it to zero, providing the likelihood equations. They are solved using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) quasi-Newton method, see mittelhammer2000. The method is implemented and applied using the R software. Under some regularity conditions Cox1979 and when $n$ is large, the asymptotic distribution of the ML estimator $\widehat{\bm{\theta}} = (\bm{\beta}(\tau)^{\top},a,q)^{\top}$ (QSM) or $\widehat{\bm{\theta}} = (\bm{\beta}(\tau)^{\top},a,p)^{\top}$ (QDA) follows asymptotically a multivariate normal distribution $$\widehat{\bm{\theta}} \,\dot{\sim}\, \textrm{N}_{k+3}(\bm{\theta}, {\bm{\Sigma}}^{-1}(\bm{\theta})),$$ where $\,\dot{\sim}\,$ means `approximately distributed' and ${\bm{\Sigma}}(\bm{\theta})$ is the expected Fisher information matrix, which is given by $$\bm{\Sigma}(\bm{\theta})= \mathbb{E}\left[- \ {\partial \ell \left(\bm{\theta}\right)\over \partial \bm{\theta} \; \partial \bm{\theta}^\top} \right].$$ A consistent estimator of $\bm{\Sigma}(\bm{\theta})$ is the estimated observed Fisher information matrix, given by $$\mathbf{K}(\widehat{\bm{\theta}})=- \ {\partial \ell \left(\bm{\theta}\right)\over \partial \bm{\theta} \;\partial \bm{\theta}^\top} \Big{|}_{\bm{\theta} = \widehat{\bm{\theta}}}.$$ Then, we can approximate $\bm{\Sigma}(\bm{\theta})$ by $\mathbf{K}(\widehat{\bm{\theta}})$.
Departures from regression models assumptions and goodness of fit are assessed by means of a residual analysis. Particularly, we use the generalized Cox-Snell (GCS) and randomized quantile (RQ) residuals:
where $F_Y$ is quantile-based Singh-Maddala or Dagum CDF, and $\widehat{\bm{\theta}}$ is the ML estimate of ${\bm{\theta}}$. If the model is correctly specified, the GCS is asymptotically standard exponential distributed, while the RQ is asymptotically standard normal distributed. With both residuals, graphical techquines, such as quantile-quantile (QQ) plots with simulated envelope, can be used to assess distributions assumptions.
In this section, we present Monte Carlo simulation studies for each reparameterized quantile model, considering different scenarios of parameters and sample sizes. The first part of the study consists in evaluating ML estimation performance, while the second evaluates the empirical distribution of the GCS and RQ residuals. Both studies consider simulated data generated from each one of the Singh-Maddala and Dagum quantile regression models according to
The Monte Carlo simulation experiments were performed using the R environment; see http://www.r-project.org.
The simulation scenario considers the following settings: sample sizes $n \in {50,100,150,250,600}$, vector of betas $\bm{\beta}(\tau) = (1, 0.5, 1.5)^{\top}$, quantiles $\tau \in \{0.10, 0.25, 0.50, 0.75, 0.90\}$, $(a,q)=(5, 1)$ (Singh-Maddala), and $(a,p)=(1, 0.5)$ (Dagum), with 500 Monte Carlo replications for each sample size. Covariate values $x_{1i},x_{2i}$ are obtained from a uniform distribution in the interval (0,1). To study the ML estimators, we use compute the relative bias (RB), root mean square error (RMSE) and the coverage probability (CP). We expect that, as sample size increases, the RB and RMSE reduces, and the CP approaches the 95% nominal level. The estimates of RB, RMSE and CP are computed from the Monte Carlo replicas as:
where $\theta$ and $\widehat{\theta}^{(i)}$ are the true parameter value and its respective $i$-th ML estimate, $m$ is the number of Monte Carlo replicas, $\mathcal{I}$ is an indicator function taking the value 1 if $\theta\in [L^{(i)}_{\widehat{\theta}},U^{(i)}_{\widehat{\theta}}]$, and 0 otherwise, where $L^{(i)}_{\widehat{\theta}}$ and $U^{(i)}_{\widehat{\theta}}$ are the $i$-th upper and lower limit estimates of the 95% confidence interval.
The results for Singh-Maddala models are shown in Figure (ref). It is possible to see that the simulations produced the expected outcomes. As the sample size increases, the RB and RMSE both decrease, and the CP tends to 95%. The results for the Dagum model are shown in Figure (ref). This figure presents results similar to those found for the Singh-Maddala model.
Here we show the performance of GCS and RQ residuals. We analyse the results with descriptive statistics (mean, median, standard deviation, coefficient of skewness and coefficient of kurtosis). The simulation scenario is exactly the same as in Subsection (ref). Figures (ref) and (ref) show the simulation results of the Singh-Maddala and Dagum models, respectively.
The reference values of mean, median, Sd, skewness and kurtosis are 1, 0.69, 1, 2 and 6, respectively, for GCS residual, and 0, 0, 1, 0 and 0, respectively, for RQ residual. From Figures (ref) and (ref), it is possible to verify that, as the sample size increases, the values tend to the expected results for each $\tau$. Therefore, we can use the both residuals to verify the fit of the proposed models.
In this section, we use the 2016 Chilean household income data set, provided by the National Institute of Statistics in Chile\footnote{Available at \url{https://www.ine.cl/estadisticas/sociales/ingresos-y-gastos/encuesta-suplementaria-de-ingresos}.} to illustrate the proposed parametric quantile regression models. This data set was also used by sanchez2021b, who introduced the Birnbaum-Saunders quantile regression model. While the Birnbaum-Saunders is not a distribution commonly used for income data, Singh-Maddala and Dagum are, so we assess if these models can produce better fits than the BS model.
The household income is the response variable ($Y$), whereas the covariates are the total income due to salaries ($X_1$), the total income due to independent work ($X_2$) and the total income due to retirements ($X_3$). The original dataset contains 107 variables, including the aforementioned, but these were selected based on economic and statistical criteria in relation to the response variable and descriptive analysis conducted by sanchez2021b. Moreover, all incomes are expressed in thousands of Chilean pesos\footnote{\ See \url{http://www.bancocentral.cl} for their equivalence in American dollars.}.
We report in Table (ref) descriptive statistics for the household income ($Y$). Figure (ref) shows the histogram along with usual and adjusted box plots rousseeuw2016. We observe that the household income data have a unimodal and right-skewed behavior, which i the precise needed scenario to uphold the usage of asymmetric distribution. Figure (ref) shows scatterplots (with correlation) between the household income ($Y$) and the covariates ($X_1$, $X_2$ and $X_3$). We observe that correlations are reasonable and significant, meanwhile the covariates have almost no linear correlation between each other.
We then analyze the household income data using the Singh-Maddala and Dagum quantile regression models, with regression structure expressed as\footnote{We use this specification in order to compare the results of the proposed models with those of the Birnbaum-Saunders quantile regression model.}
for $i,1,\ldots,100$. The proposed models are fitted using the function IncomeReg.fit, implemented in the R software R by the authors. The codes are available upon request.
Table (ref) presents the ML estimates, computed by the BFGS quasi-Newton method, standard errors (SEs) and Akaike (AIC) and Bayesian information (BIC) criteria values, for the Singh-Maddala and Dagum quantile regression models with $\tau=0.50$. As mentioned earlier, the results of the Birnbaum-Saunders quantile regression are presented as well. The results of Table (ref) reveal that the proposed Singh-Maddala and Dagum models provide better adjustments than the Birnbaum-Saunders model based on the values of log-likelihood, AIC and BIC. Particularly, the Singh-Maddala model has the lowest AIC and BIC values. The estimated parameters of the Birnbaum-Saunders, Dagum and Singh-Maddala models across $\tau$ are shown in Figure (ref). From this figure, we observe that the estimates associated with all the covariates tend to increase as $\tau$ increases, as expected.
The QQ plots with simulated envelope of the GCS and RQ residuals for the models considered in Table (ref) confirm the results presented in Table (ref); see Figure (ref). Similar results are obtained when considering $\tau = \{0.10,\ldots, 0.90\}$.
Figure (ref) shows 95% prediction intervals from the Birnbaum-Saunders, Dagum and Singh-Maddala quantile regression models for the household income data. The predictions were performed $20$-steps-ahead, namely, $20$-observations were not included in the estimation. From Figure (ref), we observe that 95%, 95% and 95% of the observations are within the limits of the prediction interval for the Birnbaum-Saunders, Dagum and Singh-Maddala models, respectively. Therefore, all the models provide values closer to the nominal 95% level.
In this paper, we have proposed parametric quantile regression models based on the Singh-Maddala and Dagum distributions. The proposed models are based on reparametrizations of the original distributions, by including the quantile as a parameter. The maximum likelihood method was used to estimate model parameters and Monte Carlo simulation studies were conducted in order to evaluate the performance of the estimators and the empirical distribution of the generalized Cox-Snell and random quantile residuals. Results showed that the estimates had good performance, and the residuals presented good agreement with their reference distributions. We applied the proposed models to a real data set, where we have modeled the household income as a function of the following covariates: total income due to salaries, total income due to independent work and total income due to retirements. The results were compared to those obtained by sanchez2021b, who proposed the Birnbaum-Saunders quantile regression model. We showed that both Singh-Maddala and Dagum models have better fit to data than Birnbaum-Saunders model, with Singh-Maddala also showing a slight superior performance than Dagum. Therefore, results were favorable to the usage of Singh-Maddala and Dagum quantile regression models. As part of future research, influence diagnostic tools can be investigated and also multivariate models can be studied.