An Asymmetric Bimodal Double Regression Model

Abstract

In this paper, we introduce an extension of the sinh Cauchy distribution including a double regression model for both the quantile and scale parameters. This model can assume different shapes: unimodal or bimodal, symmetric or asymmetric. We discuss some properties of the model and perform a simulation study in order to assess the performance of the maximum likelihood estimators in finite samples. A real data application is also presented.

1. Introduction

A wide range of phenomena can be defined more appropriately using probability distributions. This description is very useful because of the properties associated with a distribution: expectation, shape, range, etc. However, fitting data can be difficult when their distribution is bimodal, which occurs commonly in practice: in astrophysics, the metallicity of the globular cluster system in the Milky Way (see [1]); in ecology, the tree cover of moist savanna and tropical forest ecosystems (see [2]); in genetics, gene expression measurements (see [3]). Other practical examples of bimodality in data can be seen in [4,5,6]. In the literature, there are many proposals discussing bimodal distributions; e.g., the works of [7,8,9,10,11,12]. Bimodal data can be fitted by a mixture of two unimodal distributions. When the mixture is created from the same model, the main difficulty is the non-identifiability of the proposed mixture model. The traditional example is the mixture of normals. Alternatively, the most workable practical method is to use a distribution which already has bimodal properties. For the latter situation, we introduce the gamma–sinh Cauchy (GSC) distribution, proposed by [13]. We note that the initials GSC can be found in the literature as an acronym for Generalized Skew-Cauchy, and readers should be aware of this when reviewing the literature. This model has uni/bimodal properties. However, unlike the distributions discussed in the works mentioned above, in this model, one of the parameters can be interpreted as the q-th quantile under certain conditions. This is very convenient, because it allows covariates to be introduced into non-homogeneous populations directly. The probability density function (pdf) for the GSC distribution is given by

$$f(x;\mu ,\sigma ,\lambda ,\varphi )=\frac{\lambda cosh\left(\frac{x-\mu }{\sigma }\right)}{\sigma \pi \mathsf{\Gamma }\left(\varphi \right)\left\{1+{\left[\lambda sinh\left(\frac{x-\mu }{\sigma }\right)\right]}^2\right\}}{\left\{-log\left[\frac{1}{2}-\frac{1}{\pi }arctan\left\{\lambda sinh\left(\frac{x-\mu }{\sigma }\right)\right\}\right]\right\}}^{\varphi -1},$$

where $x,\mu \in \mathbb{R}$, $\sigma$, $\lambda$ and $\varphi >0$. The corresponding cumulative distribution function (cdf) is given by

$$F(x;\mu ,\sigma ,\lambda ,\varphi )=G\left(-log\left[\frac{1}{2}-\frac{1}{\pi }arctan\left\{\lambda sinh\left(\frac{x-\mu }{\sigma }\right)\right\}\right];\varphi \right),$$

where $G(·;\varphi )$ denotes the cdf of the gamma distribution with shape and scale parameters equal to $\varphi$ and 1 respectively. The GSC can be asymmetric or symmetric and, as the main advantage, its cdf has a closed-form expression which can be generated quickly in many different softwares. This is useful for generating random data, besides defining quantiles.

Regression models seek to describe the behavior of a variable of interest (or response) from covariables (explanatory variables). In general, a function called a link function links a characteristic of the response variable to the explanatory variables through parameters estimated from observed data. In our case, the response variable is bimodal and described by the GSC distribution, while the relationship between the response and explanatory variables is through the quantile. This type of relation is known as quantile regression. The literature on parametric models in the context of quantile regressions has increased considerably in recent years. For instance, for responses in the unit interval, see the works of [14,15]; for responses in the positive line see the work of [16]; and for responses in the real line see the works of [17,18].

The main advantage of quantile regression is that it is more robust against outliers. The advantage is that we can have a more informative approach to the response than simply modeling some specific measure of the population, such as the mean or median. The applicability of this method can be seen in ecology [19], econometrics [20], environmetrics [21], and medicine [22], for instance. In general terms, the distributions are not parametrized directly in terms of a general quantile $q\in (0,1)$ or any specific quantile, except for some particular cases. For instance, for the $N(\theta ,{\sigma }^2)$ model, $\theta$ represents the mean and the median of the distribution, but the q-th quantile is given by $\theta +z_q\sigma$, where $z_q$ is the q-th of the standard normal distribution. On the other hand, which quantile is of interest depends on the research. In some cases, small quantiles will be of interest, whereas in other contexts, large quantiles will be the focus. For instance, in a nutrition context, larger quantiles of weight are of interest because they enable the nutritionist to define which patients are at higher risk; on this basis, they can define special treatments for such patients.

In view of the importance of bimodal distributions, the chief object of this paper is to build on the quantile regression structure in the GSC distribution. Double regression has been studied quite extensively in the literature; for example, the authors of [23] consider a regression structure for both components based on a new parameterization indexed by mean and dispersion parameters; in [24], a regression model is proposed that is useful for situations where the variable of interest is continuous and restricted to the positive real line and is related to other variables through the mean and precision parameters; and in [25], a new parameterization of the gamma distribution is used that is indexed by mode and precision parameters.

The paper is organized as follows. In Section 2, we develop the GSC regression model with its properties. In Section 3, we perform a small-scale simulation and evaluate the point estimation. An application to real data, which illustrates the usefulness of the proposed model, is discussed in Section 4. Finally, conclusions are given in Section 5.

2. The GSC Regression Model

Gómez et al. [13] show that $F(\mu ;\mu ,\sigma ,\lambda ,\varphi )=G(log(2),\varphi )$ and then, for a fixed $q\in (0,1)$ such as

$$G(log(2),\varphi )=q,$$

(1)

the parameter $\mu$ represents the q-th quantile of the distribution. This equation can be solved numerically. Henceforth, we use the notation GSC${}_q(\mu ,\sigma ,\lambda )$ to refer to a random variable with GSC distribution where $\varphi$ is fixed as in (1), $\mu$ is the q-th quantile of the distribution, $\sigma$ is a scale parameter and $\lambda$ is a shape parameter. Figure 1 shows the relation between q and $\varphi$ and the regions where the model is unimodal or bimodal, depending on the modeled quantile q and the parameter $\lambda$. Note that, for $q\in (0,1)$, we have $({\lambda }_1,{\lambda }_2)\in {\mathbb{R}}_{+}^2$ such that the model is unimodal for ${\lambda }_1$ and bimodal for ${\lambda }_2$. In the next proposition, we enunciate a property related to the unimodality of the model.

Proposition 1.

The GSC${}_{0.5}(\mu ,\sigma ,\lambda )$model is unimodal for$\lambda \ge 1/\sqrt{2}$and bimodal for$0<\lambda <1/\sqrt{2}$.

Proof. 

Since $\mu$ and $\sigma$ are location and scale parameters, without a loss of generality, we can consider $\mu =0$ and $\sigma =1$. Deriving the pdf of the GSC${}_{0.5}(\mu =0,\sigma =1,\lambda )$ model in relation to x, we obtain that

$$f^{\prime }\left(x\right)=\frac{\partial f\left(x\right)}{\partial x}=\frac{sinh\left(x\right)(1-2{\lambda }^2-{\lambda }^2{sinh}^2\left(x\right))}{cosh\left(x\right)(1+{\lambda }^2{sinh}^2\left(x\right))}f\left(x\right).$$

Note that $f^{\prime }\left(x\right)=0$ if and only if $sinh\left(x\right)=0$ or $1-2{\lambda }^2-{\lambda }^2{sinh}^2\left(x\right)=0$. From the last equation, it follows that $sinh\left(x\right)=\pm \frac{\sqrt{1-2{\lambda }^2}}{\lambda }$, which has two solutions if $0<\lambda <1/\sqrt{2}$, one solution if $\lambda =1/\sqrt{2}$, and no solution if $\lambda >1/\sqrt{2}$. Therefore,

• For $0<\lambda <1/\sqrt{2}$, the equation $f^{\prime }\left(x\right)=0$ has three solutions. In this case, $f^{\prime }\left(x\right)>0$, $\forall x\in (-\infty ,\mathrm{asinh}\left(-\frac{\sqrt{1-2{\lambda }^2}}{\lambda }\right))\cup (0,\mathrm{asinh}\left(\frac{\sqrt{1-2{\lambda }^2}}{\lambda }\right))$ and $f^{\prime }\left(x\right)<0$, $\forall x\in (\mathrm{asinh}\left(-\frac{\sqrt{1-2{\lambda }^2}}{\lambda }\right),0)\cup (\mathrm{asinh}\left(\frac{\sqrt{1-2{\lambda }^2}}{\lambda }\right),+\infty )$. Then, $\mathrm{asinh}\left(\pm \frac{\sqrt{1-2{\lambda }^2}}{\lambda }\right)$ are the two modes of the distribution.
• For $\lambda \ge 1/\sqrt{2}$, the equation $f^{\prime }\left(x\right)=0$ has one solution. In this case, $f^{\prime }\left(x\right)>0$, $\forall x<0$ and $f^{\prime }\left(x\right)<0$, $\forall x>0$. Then, $x=0$ is the only mode of the distribution.

□

Suppose now that we are interested in modeling the q-th quantile of the distribution for a non-homogeneous population. We assume that, for a fixed $q\in (0,1)$, the q-th quantile of the distribution $\mu$ and the scale parameter $\sigma$ satisfy the following functional relations:

$${\mu }_i\left(q\right)={\mathbf{x}}_{1i}^{\top }{\mathit{\beta }}_1\left(q\right)\mathrm{and}log\left({\sigma }_i\left(q\right)\right)={\mathbf{x}}_{2i}^{\top }{\mathit{\beta }}_2\left(q\right),$$

(2)

where ${\mathit{\beta }}_1\left(q\right)={({\beta }_{11}\left(q\right),\dots ,{\beta }_{1p_1}\left(q\right))}^{\top }$ and ${\mathit{\beta }}_2\left(q\right)={({\beta }_{21}\left(q\right),\dots ,{\beta }_{2p_2}\left(q\right))}^{\top }$ are vectors of unknown regression coefficients such that $\left({\mathit{\beta }}_1{\left(q\right)}^{\top },{\mathit{\beta }}_2{\left(q\right)}^{\top }\right)\in {\mathbb{R}}^{p_1+p_2}$, with $p_1+p_2<n$; and ${\mathbf{x}}_{1i}^{\top }=(x_{11i},\dots ,x_{1p_{1}i})$ and ${\mathbf{x}}_{2i}^{\top }=(x_{21i},\dots ,x_{2p_{2}i})$ are the observations of the known regressors $p_1$ and $p_2$. Note that the vector ${\mathbf{x}}_{1i}$ is linked with the ${\mu }_i\left(q\right)$ parameter using the identity link; then, interpretations for the covariates in ${\mathbf{x}}_{1i}$ can be performed using the same idea as an ordinary linear regression. For instance, if the jth covariate $1\le j\le p_1$ is a quantitative variable, then, fixing the rest of the covariates, the q-th quantile of the distribution increases by ${\beta }_{1j}\left(q\right)$ units when $x_{1j}$ is increased to $x_{1j}+1$. Similarly, for the regression coefficients related to the scale parameter, after fixing the rest of the covariates, the scale of the distribution for the q-th quantile of the distribution is increased by $exp\left({\beta }_{2k}\left(q\right)\right)$ units when $x_{2k}$ is increased to $x_{2k}+1$, $1\le k\le p_2$. We highlight that in [13], a regression structure was assumed only for ${\mu }_i\left(q\right)$. However, the assumption that all the observations have the same scale parameter could be unrealistic, as each observation could have its own scale. For this reason, it seems reasonable to assume this double regression structure.

In this setting, the log-likelihood function for $\mathit{\psi }\left(q\right)=({\mathit{\beta }}_1^{\top }\left(q\right),{\mathit{\beta }}_2^{\top }\left(q\right),\lambda )$, up to a constant, is given by

$$\begin{array}{cc}\hfill \ell \left(\psi \left(q\right)\right)& =\sum _{i=1}^n\left[-log\left(1+{[\lambda sinh\left(z_i\left(q\right)\right)]}^2\right)\right.\hfill \\ \hfill & \left.+(\varphi -1)log\left\{-log\left[\frac{1}{2}-\frac{1}{\pi }arctan\{\lambda sinh\left(z_i\left(q\right)\right)\}\right]\right\}\right]+n\left[log\lambda -log\sigma \right],\hfill \end{array}$$

where $z_i\left(q\right)=(x_i-\mu \left(q\right))/{\sigma }_i\left(q\right)$. The maximum likelihood (ML) estimator of $\mathit{\psi }\left(q\right)$, say $\widehat{\mathit{\psi }}\left(q\right)$, is obtained by maximizing $\ell \left(\mathit{\psi }\left(q\right)\right)$ in relation to $\mathit{\psi }\left(q\right)$. For this model, such a maximization procedure does not provide a closer form, meaning that numerical procedures need to be implemented. Specifically, we use the Broyden–Fletcher–Goldfarb–Shanno (BFGS) quasi-Newton method; see [26] (p. 199). This procedure is implemented in the R software [27]. The programs are available on request. Finally, under regularity conditions, $\widehat{\mathit{\psi }}\left(q\right)$ satisfies that

$$\sqrt{n}{\left[\widehat{H}\left(\widehat{\mathit{\psi }}\right)\right]}^{-1}\left(\widehat{\mathit{\psi }}-\mathit{\psi }\right)\to N_{p_1+p_2}(\mathbf{0},{\mathit{I}}_{p_1+p-2}),asn\to +\infty ,$$

(3)

where $N_p(\mathbf{0},{\mathit{I}}_p)$ denotes the standard multivariate distribution and $\widehat{H}\left(\widehat{\mathit{\psi }}\right)$ denotes the estimated Hessian matrix of the log-likelihood function in relation to $\mathit{\psi }$.

3. Simulation Study

In this section, we present a simulation study to evaluate the performance of ML estimates in finite samples. The computational procedure is implemented using R software [27]. Values of the GSC distribution were drawn using inverse transform sampling. We considered a scheme with two covariates, both simulated from the uniform distribution between −2 and 2. We considered combinations of values for the quantile $q:0.10,0.50$ and $0.75$; vectors for ${\mathit{\beta }}_1:(1,-1,0.5),(-0.5,-2,1),(1,-1,0.5)$ and $(-0.5,-2,1)$; vectors for ${\mathit{\beta }}_2:(-1,1.6,-0.5),(-1,1.6,-0.5)$ and $(1,0.7,-0.3)$; and the parameter $log\left(\lambda \right):-1.39$ and $1.61$. We also considered three sample sizes: 100, 200 and 500. Based on 5000 replicates, we compute the mean of the estimated bias for each estimator (bias), the mean of the estimated standard errors (SE), the root of the estimated mean squared error (RMSE), and the 95% coverage probabilities (CP).

Table 1. Estimated bias, SE, RMSE, and 95%CP for the ML estimators for the GSC double regression model under different scenarios based on 5000 Monte Carlo replicates.

					$\mathit{n}=100$				$\mathit{n}=200$				$\mathit{n}=500$
q	${\mathit{\beta }}_1$	${\mathit{\beta }}_2$	$log\left(\mathit{\lambda }\right)$	Parameter	Bias	SE	RMSE	CP	Bias	SE	RMSE	CP	Bias	SE	RMSE	CP
0.10	(1, −1, 0.5)	(−1, 1.6, −0.5)	−1.39	${\beta }_{10}$	0.0146	0.0514	0.0648	0.8694	0.0066	0.0329	0.0363	0.9198	0.0025	0.0200	0.0208	0.9394
				${\beta }_{11}$	0.0060	0.0275	0.0340	0.8786	0.0029	0.0177	0.0191	0.9246	0.0010	0.0106	0.0110	0.9398
				${\beta }_{12}$	−0.0011	0.0078	0.0099	0.8676	−0.0003	0.0058	0.0066	0.9086	−0.0002	0.0038	0.0040	0.9410
				${\beta }_{20}$	−0.0544	0.1004	0.1184	0.9006	−0.0257	0.0699	0.0750	0.9314	−0.0085	0.0438	0.0449	0.9424
				${\beta }_{21}$	0.0194	0.0505	0.0598	0.9040	0.0083	0.0357	0.0385	0.9314	0.0033	0.0224	0.0236	0.9390
				${\beta }_{22}$	−0.0089	0.0457	0.0503	0.9222	−0.0051	0.0332	0.0359	0.9268	−0.0009	0.0206	0.0210	0.9462
				$log\left(\lambda \right)$	−0.0983	0.3151	0.3520	0.9332	−0.0498	0.2163	0.2273	0.9460	−0.0162	0.1340	0.1360	0.9498
			1.61	${\beta }_{10}$	0.0003	0.0449	0.0536	0.8920	−0.0002	0.0287	0.0310	0.9310	0.0002	0.0175	0.0179	0.9428
				${\beta }_{11}$	0.0003	0.0255	0.0311	0.8866	−0.0001	0.0164	0.0181	0.9222	0.0001	0.0099	0.0102	0.9410
				${\beta }_{12}$	−0.0001	0.0069	0.0098	0.8308	0.0001	0.0053	0.0064	0.8950	0.0000	0.0036	0.0039	0.9334
				${\beta }_{20}$	−0.0627	0.1116	0.1311	0.8904	−0.0315	0.0791	0.0863	0.9162	−0.0113	0.0501	0.0514	0.9398
				${\beta }_{21}$	0.0113	0.0460	0.0504	0.9264	0.0047	0.0337	0.0352	0.9386	0.0014	0.0215	0.0216	0.9472
				${\beta }_{22}$	−0.0058	0.0467	0.0503	0.9280	−0.0023	0.0336	0.0350	0.9388	−0.0012	0.0211	0.0215	0.9464
				$log\left(\lambda \right)$	−0.1492	0.3012	0.3461	0.9266	−0.0777	0.2073	0.2277	0.9312	−0.0274	0.1288	0.1333	0.9446
0.50	(−0.5, −2, 1)	(−1, 1.6, −0.5)	−1.39	${\beta }_{10}$	0.0126	0.0512	0.0636	0.8744	0.0073	0.0331	0.0368	0.9182	0.0023	0.0200	0.0204	0.9398
				${\beta }_{11}$	0.0050	0.0275	0.0338	0.8812	0.0032	0.0178	0.0195	0.9230	0.0009	0.0106	0.0108	0.9422
				${\beta }_{12}$	−0.0011	0.0077	0.0098	0.8628	−0.0005	0.0058	0.0067	0.9118	−0.0003	0.0039	0.0041	0.9350
				${\beta }_{20}$	−0.0576	0.1001	0.1197	0.8924	−0.0239	0.0700	0.0748	0.9288	−0.0087	0.0438	0.0446	0.9446
				${\beta }_{21}$	0.0186	0.0499	0.0586	0.9074	0.0080	0.0358	0.0381	0.9342	0.0036	0.0224	0.0236	0.9376
				${\beta }_{22}$	−0.0084	0.0455	0.0489	0.9274	−0.0044	0.0333	0.0352	0.9346	−0.0011	0.0207	0.0207	0.9492
				$log\left(\lambda \right)$	−0.1142	0.3149	0.3507	0.9382	−0.0432	0.2162	0.2242	0.9448	−0.0154	0.1340	0.1391	0.9422
		(1, 0.7, −0.3)	−1.39	${\beta }_{10}$	−0.0068	0.0529	0.0668	0.8808	−0.0020	0.0346	0.0384	0.9198	−0.0008	0.0210	0.0220	0.9428
				${\beta }_{11}$	−0.0025	0.0297	0.0382	0.8720	−0.0005	0.0196	0.0223	0.9164	−0.0003	0.0118	0.0124	0.9428
				${\beta }_{12}$	0.0002	0.0082	0.0121	0.8276	0.0001	0.0063	0.0079	0.8782	0.0000	0.0043	0.0047	0.9342
				${\beta }_{20}$	−0.0667	0.1081	0.1297	0.8862	−0.0275	0.0765	0.0813	0.9306	−0.0102	0.0482	0.0497	0.9402
				${\beta }_{21}$	0.0154	0.0511	0.0574	0.9182	0.0048	0.0375	0.0393	0.9352	0.0023	0.0238	0.0245	0.9430
				${\beta }_{22}$	−0.0070	0.0512	0.0560	0.9208	−0.0047	0.0371	0.0386	0.9426	−0.0010	0.0232	0.0234	0.9510
				$log\left(\lambda \right)$	−0.1579	0.3161	0.3625	0.9240	−0.0694	0.2164	0.2292	0.9398	−0.0227	0.1342	0.1375	0.9440
0.75	(1, −1, 0.5)	(1, 0.7, −0.3)	1.61	${\beta }_{10}$	−0.0061	0.5813	0.6422	0.9130	−0.0004	0.4112	0.4314	0.9352	−0.0020	0.2552	0.2604	0.9442
				${\beta }_{11}$	0.0032	0.3709	0.4267	0.9004	0.0015	0.2696	0.2875	0.9296	0.0000	0.1632	0.1668	0.9426
				${\beta }_{12}$	0.0059	0.2288	0.2930	0.8690	0.0061	0.1650	0.1863	0.9152	−0.0001	0.1074	0.1115	0.9414
				${\beta }_{20}$	−0.0629	0.1118	0.1297	0.8912	−0.0292	0.0792	0.0866	0.9180	−0.0127	0.0500	0.0510	0.9422
				${\beta }_{21}$	0.0079	0.0458	0.0497	0.9274	0.0033	0.0338	0.0351	0.9400	0.0019	0.0214	0.0219	0.9490
				${\beta }_{22}$	−0.0037	0.0467	0.0503	0.9292	−0.0012	0.0337	0.0354	0.9382	−0.0007	0.0211	0.0213	0.9462
				$log\left(\lambda \right)$	−0.1481	0.3016	0.3408	0.9286	−0.0727	0.2073	0.2253	0.9362	−0.0313	0.1289	0.1305	0.9508
	(−0.5, −2, 1)	(1, 0.7, −0.3)	1.61	${\beta }_{10}$	−0.0043	0.1204	0.1302	0.9226	−0.0013	0.0834	0.0886	0.9322	−0.0007	0.0516	0.0519	0.9474
				${\beta }_{11}$	0.0026	0.0772	0.0880	0.9102	0.0009	0.0547	0.0594	0.9242	0.0005	0.0330	0.0333	0.9462
				${\beta }_{12}$	−0.0013	0.0494	0.0627	0.8624	−0.0006	0.0338	0.0366	0.9208	0.0001	0.0217	0.0224	0.9370
				${\beta }_{20}$	−0.0547	0.2266	0.2380	0.9246	−0.0307	0.1588	0.1659	0.9306	−0.0099	0.1000	0.0997	0.9500
				${\beta }_{21}$	0.0264	0.1123	0.1207	0.9282	0.0098	0.0805	0.0833	0.9414	0.0045	0.0504	0.0517	0.9428
				${\beta }_{22}$	−0.0142	0.1147	0.1221	0.9342	−0.0080	0.0803	0.0843	0.9382	−0.0014	0.0496	0.0508	0.9470
				$log\left(\lambda \right)$	−0.0064	0.2947	0.2999	0.9452	−0.0107	0.2041	0.2084	0.9432	−0.0028	0.1282	0.1277	0.9488

summarizes the results. From Table 1, it can be observed that the bias, SE, and RMSE for all the parameters tend to approach zero when the sample size is increased, showing that the ML estimates obtained are asymptotically consistent. On the other hand, the CP values are closer to the nominal values used in their construction (95%), suggesting that the asymptotic distribution in Equation (3) is reasonable, even in finite samples.

4. Application

To illustrate the GSC double regression model, we consider the Australian data set available in the package sn in R [28], which includes data on 202 athletes collected at the Australian Institute of Sport. Codes were performed in [27] and are available upon request. Our main aim is to explain the body fat percentage (Bfat) in terms of the body mass index (bmi) and the lean body mass (lbm). Particularly, we consider Bfat${}_i\left(q\right)\sim GSC({\mu }_i\left(q\right),{\sigma }_i\left(q\right),\lambda ,\varphi \left(q\right))$, where $\varphi \left(q\right)$ satisfies (1), $q\in (0,1)$ and for $i=1,\dots ,202$, we have that

$${\mu }_i\left(q\right)={\beta }_{11}\left(q\right)+{\beta }_{12}\left(q\right){\mathtt{bmi}}_i+{\beta }_{13}\left(q\right){\mathtt{lbm}}_i \mathrm{and} {\sigma }_i\left(q\right)={\beta }_{21}\left(q\right)+{\beta }_{22}\left(q\right){\mathtt{bmi}}_i+{\beta }_{23}\left(q\right){\mathtt{lbm}}_i.$$

In other words, the bmi and lbm explain both the q-th quantile of Bfat and the scale of the distribution. The same structure of covariates was considered in [13], but without modeling the scale parameter; i.e., considering ${\beta }_{22}\left(q\right)={\beta }_{23}\left(q\right)=0$. We refer to those models as GSC and GSC${}_0$ for the cases where $\sigma$ is modeled and not modeled, respectively. We considered $q\in \{0.1,0.25,0.5,0.75,0.9\}$. Our approach is compared with the skewed Laplace (SKL) model in [29].

Table 2. AIC for GSC, GSC${}_0$, and SKL models in the athlete data set for different quantiles. We also present the statistical p-value.

	AIC			log-Likelihood		LRT		KS
$\mathit{\tau }$	GSC${}_0$	GSC	SKL	GSC${}_0$	GSC	Statistical	p-Value	p-Value
0.10	1168.54	1154.64	1194.28	−574.27	−563.32	21.90	<0.0001	0.988
0.25	1172.72	1164.55	1172.70	−576.36	−568.27	16.17	0.0003	0.646
0.50	1174.74	1171.83	1182.66	−577.37	−571.91	10.91	0.0043	0.180
0.75	1171.50	1176.51	1221.65	−575.75	−574.25	3.00	0.2235	0.839
0.90	1223.71	1229.75	1280.45	−601.85	−600.88	1.95	0.3768	0.004

shows the Akaike information criterion (AIC [30]) for the three models. We also present the statistic for the likelihood ratio test (LRT) to test $H_0:{\beta }_{22}\left(q\right)={\beta }_{23}\left(q\right)=0$ versus $H_0:{\beta }_{22}\left(q\right)\ne 0$ or ${\beta }_{23}\left(q\right)\ne 0$, for all the quantiles considered. In addition, we also compute the quantile residuals [31] for the GSC model. If the model is correctly specified, the residual should be a random sample from the standard normal distribution. We checked this assumption with the traditional Kolmogorov–Smirnov test. Note that the GSC presents the lowest AIC for the quantiles up to the median, and GSC${}_0$ presents the lowest AIC for the rest of the quantiles. This is explained because, according to the LRT, the coefficients related to the bmi and lbm variables are not significant (under any common level of significance) for modeling the scale parameter for $q=0.75$ and $q=0.9$, while they are significant for the rest of the quantiles. Finally, based on the quantile residuals, the GSC double regression model seems to be appropriate for modeling all quantiles, except the largest. Figure 2 also shows the regression coefficients in terms of the quantiles and their respective 95% confidence intervals. Note that ${\beta }_{12}\left(q\right)$ and ${\beta }_{13}\left(q\right)$ are significant (based on 5% significance) for all the quantiles considered; i.e., the bmi and lbm variables are relevant for explaining the different quantiles of Bfat. Specifically, we can obtain the following interpretations for ${\beta }_{12}\left(q\right)$ and ${\beta }_{13}\left(q\right)$:

• (Interpreting ${\beta }_{12}\left(q\right)$) For a fixed lbm, for athletes in the lowest 10% of Bfat, the Bfat is increased by 0.8572 units (95% confidence interval 0.4191; 1.2954) for each unit increase in bmi, and for athletes in the highest 90% of Bfat, the Bfat is increased by 2.5834 units (95% confidence interval 2.3227; 2.8442) for each unit increase in bmi.
• (Interpreting ${\beta }_{13}\left(q\right)$) For a fixed bmi, for athletes in the lowest 10% of Bfat, the Bfat is decreased by 0.3039 units (95% confidence interval −0.4093; −0.1984) for each unit increase in lbm, and for athletes in the highest 10% of Bfat, the Bfat is decreased by 0.5511 units (95% confidence interval −0.6077; −0.4945) for each unit increase in lbm.

We highlight the large difference between the interpretations for athletes in the lowest 10% and highest 10% of Bfat.

Finally, Figure 3 shows the different estimated pdf values for the Bfat under different scenarios for bmi and lbm. We note that these estimated pdf values assumed different shapes—unimodal, bimodal symmetric, and bimodal asymmetric—justifying the use of the double regression GSC model in this example. Finally, Figure 4 shows the pairs $(\lambda ,q)$ for the five quantiles modeled, identifying the unimodal and bimodal cases.

5. Conclusions

In this paper, we present a new extension of the GSC model, introducing a double regression structure to model both the quantile and the scale of the distribution. This structure produces a competitive model for modeling heterogeneous populations with different shapes: unimodal symmetric, unimodal asymmetric, bimodal symmetric, and bimodal asymmetric. The illustration with a real data set shows that the model provides better performance than other proposals in the literature. A limitation of the model is that the shape parameter is common for all the observations. Further extensions should include covariates in this parameter also.