Next Article in Journal
Insights into the 3D Slip Dynamics of Jeffrey Fluid Due to a Rotating Disk with Exponential Space-Dependent Heat Generation: A Case Involving a Non-Fourier Heat Flux Model
Next Article in Special Issue
Some Extensions of the Asymmetric Exponentiated Bimodal Normal Model for Modeling Data with Positive Support
Previous Article in Journal
Adaptive Backstepping Terminal Sliding Mode Control of Nonlinear System Using Fuzzy Neural Structure
Previous Article in Special Issue
The Slash Half-Normal Distribution Applied to a Cure Rate Model with Application to Bone Marrow Transplantation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

The Multivariate Skewed Log-Birnbaum–Saunders Distribution and Its Associated Regression Model

by
Guillermo Martínez-Flórez
1,*,†,
Sandra Vergara-Cardozo
2,
Roger Tovar-Falón
1,*,† and
Luisa Rodriguez-Quevedo
2
1
Departamento de Matemáticas y Estadística, Facultad de Ciencias Básicas, Universidad de Córdoba, Montería 230002, Colombia
2
Departamento de Estadística, Facultad de Ciencias, Universidad Nacional de Colombia, Bogotá 111321, Colombia
*
Authors to whom correspondence should be addressed.
These authors contributed equally to this work.
Mathematics 2023, 11(5), 1095; https://doi.org/10.3390/math11051095
Submission received: 19 October 2022 / Revised: 28 November 2022 / Accepted: 6 December 2022 / Published: 22 February 2023
(This article belongs to the Special Issue Probability, Statistics & Symmetry)

Abstract

:
In this article, a multivariate extension of the unit-sinh-normal (USHN) distribution is presented. The new distribution, which is obtained from the conditionally specified distributions methodology, is absolutely continuous, and its marginal distributions are univariate USHN. The properties of the multivariate USHN distribution are studied in detail, and statistical inference is carried out from a classical approach using the maximum likelihood method. The new multivariate USHN distribution is suitable for modeling bounded data, especially in the ( 0 , 1 ) p region. In addition, the proposed distribution is extended to the case of the regression model and, for the latter, the Fisher information matrix is derived. The numerical results of a small simulation study and two applications with real data sets allow us to conclude that the proposed distribution, as well as its extension to regression models, are potentially useful to analyze the data of proportions, rates, or indices when modeling them jointly considering different degrees of correlation that may exist in the study variables is of interest.

1. Introduction

Data whose response falls in the interval ( 0 , 1 ) such as indices, proportions, or rates appear very frequently in different fields of knowledge, mainly the areas of social sciences, engineering, economic sciences, and medicine. Some practical examples of these types of data are the proportion of patients who die from a certain disease or virus (SARS-CoV-2, Diabetes, HIV or Cancer) in a country or city; the Human Development Index or the illiteracy rate in a region or country; the proportion of deaths due to exposure to smoking or other exposure factors; the mortality rate from traffic accidents in a city; the percentage of items that do not meet the minimum requirements in an assembly line; and the portion of income that a family spends on entertainment.
For the analysis of data such as those described above, statistical methodologies developed from distributions with support in the interval ( 0 , 1 ) are required. In this sense, several probability distributions and regression models have been proposed; see Ferrari and Cribari-Neto [1], Kumaraswamy [2], Martínez-Flórez et al. [3,4], Kieschnick and Mccullough [5], Mazucheli et al. [6,7].
The univariate sinh-normal (SHN) distribution introduced by Rieck and Nedelman [8] has received special attention for modeling material-fatigue-related problems. The probability density function (pdf) of this distribution is given by:
φ ( y ) = 2 σ α cosh y γ σ ϕ 2 α sinh y γ σ , y R ,
where α > 0 and γ are shape and location parameters, respectively; σ > 0 is a scale parameter, and ϕ ( · ) is the pdf of the normal distribution. The pdf in (1) can also be written as:
φ ( y ) = b y ϕ ( b y )
where b y = 2 σ α cosh y γ σ is a derivative of b y = 2 α sinh y γ σ . The distribution function in (2) is denoted by SHN ( α , γ , σ ) . Several extensions of the SHN distribution have been studied by numerous authors, for example, Martínez-Flórez et al. [9] proposed the extended generalized SHN distribution, which has great applicability in the fit of datasets presenting high skewness and bimodality simultaneously; Lemonte [10] introduced an extension named skewed log-Birnbaum–Saunders (log-BS) regression model which is based on the asymmetric SHN distribution proposed by Leiva et al. [11]. The log-BS regression model is suitable for fitting data with high degrees of skewness. On the other hand, Moreno et al. [12] proposed a generalization of the Birnbaum–Saunders (BS) distribution [13] that affords flexibility for fitting data with greater skewness and kurtosis compared with other distributions.
Multivariate extensions of the SHN distribution have also been considered, for instance, by Martínez-Flórez et al. [14], Díaz-García and Domínguez-Molina [15], Lemonte [16], Marchant et al. [17] and, recently, by Martínez-Flórez et al. [18], among others.
For modeling material fatigue data, the most widely known distribution is the BS whose pdf is given by:
f T ( t ) = ϕ ( a t ) t 3 / 2 ( t + β ) 2 α β , t > 0 ,
where a t = 1 α t / β β / t , α > 0 is a shape parameter, and β > 0 is a scale parameter and the median distribution. The distribution in (3) is denoted by T BS ( α , β ) . Rieck and Nedelman [8] showed that, if Y SHN ( α , γ , σ = 2 ) , then T = exp ( Y ) follows a BS distribution with parameters α and β = exp ( γ ) . From this last relationship between the BS and SHN distributions, the SHN regression model can be formulated as follows: if x i = ( x i 1 , , x i p ) is a vector of covariates such that,
Y i = log ( T i ) = x i θ ,
for i = 1 , , n , and,
Y i SHN ( α , x i θ , σ )
then, the model in (4) is known as the log-linear BS regression model. More details about this regression model can be found in Rieck and Nedelman [8].
The BS distribution has great applicability to analyze data in several areas of knowledge, such as biology, medicine, engineering, etc.; however, so far, no extension of the BS distribution has been proposed to study the modeling of rates and proportions, i.e., of a random variable in the unit interval ( 0 , 1 ) from the BS model. In response to this special case, Mazucheli et al. [19] presented an extension of the BS distribution for fitting random variables in the unit interval ( 0 , 1 ) . The pdf of this model is given by:
f ( x ) = 1 2 x α β 2 π β log x 1 / 2 + β log x 3 / 2 exp 1 2 α 2 log x β + β log x + 2 ,
where x ( 0 , 1 ) , α > 0 is the shape parameter and β > 0 is a scale parameter. Based on the work of Mazucheli et al. [19], Martínez-Flórez et al. [3] studied the unit sinh-normal (USHN) distribution to deal with the problem of bounded observations on the interval ( 0 , 1 ) which has a pdf given by:
φ ( y ) = 1 ( 1 y ) log ( 1 y ) 1 2 σ α cosh log ( log ( 1 y ) ) γ σ × ϕ 2 α sinh log ( log ( 1 y ) ) γ σ ,
where y ( 0 , 1 ) , α > 0 is a shape parameter, γ is a location parameter, and σ > 0 is a scale parameter. The natural extension to the case of the model which considers covariates is the USHN linear regression (USHNR) model. The USHNR model is defined by considering a set of p explanatory variables that are denoted by x i = ( x i 1 , , x i p ) and, such that,
log log ( 1 Y i ) = x i θ + ε i , for i = 1 , , n
where θ = ( θ 1 , , θ p ) is a p dimensional vector of unknown parameters and ε i SHN ( α , 0 , σ ) . More details can be found in [3].
The main objective of this work is to introduce a new multivariate probability distribution capable of modeling data in the region ( 0 , 1 ) p . The new distribution is obtained from the extension of the univariate skewed unit-sinh-normal distribution and to do so, we rely on the conditionally specified distributions methodology introduced by Arnold. In addition, from the new distribution, we propose the multivariate unit-sinh-normal skewed regression model, which allows modeling data in the region ( 0 , 1 ) p through linear predictors. The new proposals are useful for the analysis of data on proportions, rates, or indices that arise in different fields of knowledge, such as those described at the beginning of this section. The results of a simulation presented in this work also show that these methodologies are viable alternatives to those existing in the current statistical literature.
This paper is organized as follows. In Section 2, the multivariate skew-normal distribution is revised, and its main properties are commented on. In Section 3, the new multivariate skewed unit sinh-normal distribution is introduced. Some properties are also derived, and the value of the coefficient correlation for the bivariate case is presented for some selected values of the parameter distribution. Section 4 presents the extension of the USHN to the case of the multivariate regression model and its respective statistical inference. Finally, two applications with real data to illustrate the applicability of the proposed methodologies and a small simulation study are presented in Section 5.

2. Multivariate Skew-Normal Distribution

The multivariate skew-normal (SN) distribution was studied by Arnold et al. [20] by using the theory of conditionally specified distributions; see [21]. The construction of the multivariate SN distribution is as follows: for each j = 1 , 2 , , p , define the vector Z ( j ) to be the ( p 1 ) dimensional random vector obtained from Z by deleting Z j . In parallel, for a real vector z = ( z 1 , z 2 , , z p ) , z ( j ) is obtained by deleting the jth coordinate z j of z . Now, suppose that, for each j = 1 , 2 , , p , the conditional distribution of Z j given Z ( j ) = z ( j ) is a SN distribution with a parameter which is a function of z ( j ) . Thus, it is assumed for each j that
Z j Z ( j ) = z ( j ) SN λ j j z j .
The joint pdf of Z = ( Z 1 , Z 2 , , Z p ) is given by
f Z ( z ) = 2 j = 1 p ϕ ( z j ) Φ λ j = 1 p z j .
In the distribution (8), the marginal densities follow a standard normal distribution, that is, for j = 1 , 2 , , p , and
Z j N ( 0 , 1 ) ,
the conditional distribution follows a SN distribution, (see Azzalini [22]) of parameter λ j j z j , with the pdf given by:
f ( Z j Z ( j ) = z ( j ) ) = 2 ϕ ( z j ) Φ λ j j z j
From this distribution, Lemonte et al. [23] presented the multivariate Birnbaum–Saunders distribution, whose joint pdf is given by:
f T 1 , , T p ( t 1 , , t p ) = 2 j = 1 p ϕ ( a t j ) Φ λ j = 1 p a t j j = 1 p t j 3 / 2 ( t j + η j ) 2 α j η j , t 1 , , t p > 0 ,
where
a t j ( α j , η j ) = a t j = 1 α j t j η j η j t j ,
for j = 1 , 2 , , p , with α j > 0 and η j > 0 being the shape and scale parameters, respectively. The distribution (9) is denoted by MVBS ( α , η , λ ) .
Another extension based on the multivariate SN model of Arnold et al. [20] is the multivariate asymmetric SHN distribution, studied by Martínez-Flórez et al. [24], whose joint pdf is given by:
f Y 1 , , Y p ( y 1 , , y p ) = 2 j = 1 p b j j = 1 p ϕ ( b j ) Φ λ j = 1 p b j , y 1 , , y p R ,
where b j = 2 α j sinh Y j γ j σ j , and b j = 2 α j σ j cosh Y j γ j σ j for j = 1 , 2 , , p is a derivative of b j with respect to Y j ; α j > 0 and σ j > 0 are shape and scale parameters, respectively, and γ j , λ R are location and asymmetry parameters, respectively. The distribution in (11) is denoted by MVSHN ( α , γ , σ , λ ) , with α = ( α 1 , , α p ) , γ = ( γ 1 , , γ p ) and σ = ( σ 1 , , σ p ) .
Although MVBS and MVSHN distributions, which are defined in R 2 p + 1 and R 3 p + 1 , respectively, can be used to fit sets of random variables whose domain is the unit interval ( 0 , 1 ) , these are not appropriate given the support of these distributions and the support of a bounded random variable vector. In the statistical literature, there are few distributions studied to fit sets of variables in the unit interval, that is, whose domain of definition is ( 0 , 1 ) p , which can be useful to fit rates and proportions. The interest for these type of distributions has been very little; we highlight the works of Cepeda et al. [25], Souza and Moura [26] and Lemonte and Moreno-Arenas [27], among others.

3. Multivariate Skewed Unit-Sinh-Normal Distribution

Following the idea of Arnold et al. [20], Lemonte et al. [23] and Martínez-Flórez et al. [24], in this section, a multivariate extension of the SHN distribution to fit vectors of rates and proportions is proposed, which is named multivariate skewed USHN distribution (MVSUSHN). The construction of the MVSUSHN is as follows: for j = 1 , 2 , , p , let
Y j = 1 exp exp γ j + σ j sinh 1 α j Z j 2 ,
where Z j N ( 0 , 1 ) for j = 1 , 2 , , p . Then, the joint pdf of the vector with MVSUSHN distribution is given by,
f Y 1 , , Y p ( y 1 , , y p ) = 2 j = 1 p b j j = 1 p ϕ ( b j ) Φ λ j = 1 p b j , y 1 , , y p ( 0 , 1 ) ,
where
b j = 2 α j sinh log ( log ( 1 y j ) ) γ j σ j
and
b j = 2 α j σ j ( 1 y j ) ( log ( 1 y j ) ) cosh log ( log ( 1 y j ) ) γ j σ j
for j = 1 , 2 , , p is the derivative of b j with respect to y j ; α j > 0 and σ j > 0 are shape and scale parameters, and γ j and λ R are location and asymmetry parameters, respectively. The MVSUSHN is denoted as MVSUSHN ( α , γ , σ , λ ) , with α = ( α 1 , , α p ) , γ = ( γ 1 , , γ p ) and σ = ( σ 1 , , σ p ) . For λ = 0 , the case of independence is obtained, that is, the product of the pdf of USHN random variables studied by Martínez-Flórez et al. [3]. It follows that the parameter λ is directly associated with the correlation parameter. The Figure 1 shows the contours of the bivariate skewed USHN (BVSUSHN) distribution for some selected values of the parameters, while the Figure 2 presents the shape of the density function for particular values of the parameter of the the BVSUSHN distribution.
The following theorem provides the marginal and conditional distributions of the MVSUSHN distribution.
Theorem 1.
If ( Y 1 , Y 2 , , Y n ) M V S U S H N ( α , γ , σ , λ ) then,
(1)
Y j U S H N ( α j , γ j , σ j ) for j = 1 , 2 , , n .
(2)
The conditional pdf of Y j Y ( j ) = y ( j ) is given by
f Y j Y ( j ) ( y j Y ( j ) = y ( j ) ) = 2 b j ϕ ( b j ) Φ λ j = 1 n b j .
(3)
The cumulative distribution function (cdf) of Y j Y ( j ) = y ( j ) is given by
P ( Y j y j Y ( j ) = y ( j ) ) = Φ ( b j ) 2 T b j , λ j j b j ,
where T ( · ) is the Owen function; see [28].
Proof. 
(1)
For k = 1 , 2 , , p and applying the integral over all subindex k (given by j ), other than j, we obtain,
f Y j ( y j ) = ( 0 , 1 ) j j ( 0 , 1 ) 2 k = 1 p b k ϕ ( b k Φ λ k = 1 p b k j j d y j = b j ϕ ( b j ) ( 0 , 1 ) j j ( 0 , 1 ) 2 j j b j ϕ ( b j ) Φ ( λ b j ) j j b j j j d y j .
Now, using the transformation Z j = 2 α j sinh log ( log ( 1 Y j ) ) γ j σ j for all j j
f Y j ( y j ) = b j ϕ ( b j ) R j j R 2 j j ϕ ( z j ) Φ ( λ z j ) j j z j j j d z j = b j ϕ ( b j ) ( 1 ) = b j ϕ ( b j ) .
where the second last result follows from Arnold et al. [20].
(2)
Let
f Z j Z ( j ) ( Z j Z ( j ) = z ( j ) ) = 2 ϕ ( z j ) Φ λ j j z j z j ,
then, with the transformation Y j = 1 exp exp γ j + σ j sinh 1 α j Z j σ j , it is found that Z j = 2 α j sinh log ( log ( 1 Y j ) ) γ j σ j = b j and d Z j d Y j = b j and, by the Transformation Theorem, it follows:
f Y j Y ( j ) ( Y j Y ( j ) = y ( j ) ) = 2 b j ϕ ( b j ) Φ λ j = 1 p b j .
(3)
It has that
P ( Y j y j Y ( j ) = y ( j ) ) = y j f Y j ( t j Y ( j ) = y ( j ) ) d t j
through the transformation Z j = 2 α j sinh log ( log ( 1 T j ) ) γ j σ j = b j , it follows that
P ( Y j y j Y ( j ) = y ( j ) ) = b j 2 ϕ ( z j ) Φ λ j j b j z j d z j = Φ ( b j ) 2 T ( b j , λ j j b j )
where the last equality follows the properties of the cdf of the SN distribution, which is widely known in the literature.
 □
Martínez-Flórez et al. [3] showed that, when α 0 , the random variable
log ( log ( 1 Y ) ) γ α σ / 2 ,
converges to a standard normal distribution. From here, if α j 0 with j = 1 , 2 , , p , then
Y MVSN ( γ , σ , λ ) .
Furthermore, if X j = log ( log ( 1 Y j ) ) for j = 1 , 2 , , p , then
X MVSHN ( α , γ , σ , λ ) .
If Y MVSUSHN ( α , γ , σ , λ ) then, Z j = 2 α j sinh log ( log ( 1 Y j ) ) γ j σ j N ( 0 , 1 ) , for all j = 1 , 2 , , p , then,
Z = ( Z 1 , Z 2 , , Z p ) MVSN ( 0 , I , λ ) .
(see [3]). It can be seen that, if σ 1 = σ 2 = = σ q = 2 , the multivariate standard USHN distribution follows, which is denoted by MVSUSHN ( α , γ , 2 1 q , λ ) . In this case, the marginals are standard USHN, and the variables X j = log ( log ( 1 Y j ) ) for j = 1 , 2 , , p follow the SHN distribution of Rieck and Nedelman [8].
In order to study the unimodal distribution, let p = 2 , and suppose that α 1 = α 2 , γ 1 = γ 2 = 0 , and σ 1 = σ 2 = 1 . By differentiating the logarithm of the conditional distributions and equaling to zero, it has
log f ( y 2 y 1 ) y 2 = ( 1 + log ( 1 y 2 ) ) + b y 2 b y 2 b y 2 b y 2 + λ b y 2 b y 1 ϕ ( λ b y 1 b y 2 ) Φ ( λ b y 1 b y 2 ) , log f ( y 1 y 2 ) y 1 = ( 1 + log ( 1 y 1 ) ) + b y 1 b y 1 b y 1 b y 1 + λ b y 1 b y 2 ϕ ( λ b y 1 b y 2 ) Φ ( λ b y 1 b y 2 ) ,
the following equations are obtained:
b y 2 b y 2 b y 2 b y 2 + λ b y 2 b y 1 ϕ ( λ b y 1 b y 2 ) Φ ( λ b y 1 b y 2 ) = ( 1 + log ( 1 y 2 ) )
b y 1 b y 1 b y 1 b y 1 + λ b y 1 b y 2 ϕ ( λ b y 1 b y 2 ) Φ ( λ b y 1 b y 2 ) = ( 1 + log ( 1 y 1 ) )
Multiplying the equation in (15) by b y 2 b y 1 2 and the equation in (16) by b y 1 b y 2 2 , and subtracting these two results, it follows that
b y 2 2 b y 1 2 b y 2 2 b y 2 2 b y 1 2 b y 2 2 b y 1 2 + b y 2 2 b y 1 2 b y 1 2 = ( 1 + log ( 1 y 1 ) ) b y 2 2 b y 1 b y 1 ( 1 + log ( 1 y 2 ) ) b y 2 b y 2 b y 1 2 .
Note that by letting y 2 = y 1 and substituting in (17), it follows:
b y 1 2 b y 1 2 b y 1 2 b y 1 4 b y 1 2 b y 1 2 + b y 1 4 b y 1 2 = ( 1 + log ( 1 y 1 ) ) b y 1 3 b y 1 ( 1 + log ( 1 y 1 ) ) b y 1 b y 1 3 0 = 0
Therefore, y 1 = y 2 is a trivial solution of the Equation (17). Then, by replacing y 1 = y 2 in (15), it has
b y 2 ( 1 b y 2 ) Φ ( λ b y 2 ) + ( 1 + log ( 1 y 2 ) ) b y 2 Φ ( λ b y 2 ) + λ b y 2 b y 2 ϕ ( λ b y 2 ) = 0
from which results the function
g ( y 2 ; λ ) = b y 2 ( 1 b y 2 ) Φ ( λ b y 2 ) + ( 1 + log ( 1 y 2 ) ) b y 2 Φ ( λ b y 2 ) + λ b y 2 b y 2 ϕ ( λ b y 2 ) .
Then, by applying the transformation Y j = 1 exp ( exp ( σ j arcsinh ( α j Z j / 2 ) + γ j ) ) , the bivariate distribution with conditional SN distributions is obtained.
The bivariate distribution with conditional asymmetric USHN is a one-to-one transformation. The λ values for which the SBVUSHN distribution is unimodal are the same as for which the bivariate SN distribution is unimodal and, according to Arnold et al. [20], the bivariate SN distribution is unimodal for λ π / 2 . One can note that the equation g ( y 2 ; λ ) is similar to the equation found by Arnold et al. [20] for the BVSN distribution. The modes of the BVSUSHN distribution can be obtained by solving the equation g ( y 2 ; λ ) = 0 and Y 1 Y 2 = 0 .

Moments and Correlation

The covariance for the random variables Y j and Y j is given by:
cov ( Y j , Y j ; λ ) = E ( Y j Y j ) E ( Y j ) E ( Y j )
where the moment product E ( Y j Y j ) for two random variables is given by
E ( Y j Y j ) = 2 ( 0 , 1 ) ( 0 , 1 ) y j y j b y j b y j ϕ ( b y j ) ϕ ( b y j ) Φ λ b y j b y j d y j d y j
and
E ( Y r ) = j = 0 r l = 0 r j ( 1 ) j + l ( j e γ ) l l ! k a 1 ( α 2 ) + k b 1 ( α 2 ) k 1 / 2 ( α 2 )
with a = r σ + 1 2 , b = r σ 1 2 , and k λ ( · ) being the third-order function of Bessel defined by
k λ ( v ) = 1 2 v 2 λ 0 u λ 1 e u v 2 4 u d u .
A proof of the result in (18) can be seen in the Appendix A. Using (18) the variances: var ( Y i ; λ ) = E ( Y i 2 ) E 2 ( Y i ) are obtained for i = j , j ; thus, the correlation coefficient is obtained from:
cor ( Y j , Y j ; λ ) = cov ( Y j , Y j ; λ ) var ( Y j ; λ ) var ( Y j ; λ ) .
To compute cor ( Y j , Y j ; λ ) , it is necessary to use numerical methods to determine the simple moments and the product moments. It can be shown that for this distribution cov ( Y j , Y j ; λ ) = cov ( Y j , Y j ; λ ) , whereby cor ( Y j , Y j ; λ ) = cor ( Y j , Y j ; λ ) . To study the range of the correlation coefficient, cor ( Y j , Y j ; λ ) was evaluated for some parameter values. The values taken for the parameters: σ 1 = σ 2 = 1 and γ 1 = γ 2 = 0 , α 1 , α 2 , and λ varying. The Table 1 shows the values for the parameters α 1 , α 2 , λ , and the values of cor ( Y j , Y j ; λ ) for a pair of variables Y j and Y j . The values were obtained for the case of λ > 0 and the case of λ < 0 results for symmetry, given that cor ( Y j , Y j ; λ ) = cor ( Y j , Y j ; λ ) .
As can be seen, the case of independence occurs for λ = 0 , i.e., cor ( Y j , Y j ; 0 ) = 0 , then according to the Table 1, for the values | cor ( Y j , Y j ; λ ) | 0.9938 , which leads us to the conclusion that for this model | cor ( Y j , Y j ; λ ) | 1.0 .

4. Multivariate Skewed USHN Regression Model

This section presents an extension of the USHN regression model for the case of multiple bounded response variables (rates and proportions). Suppose that we have q variables measuring rates or proportions in a sample of size n , i.e., for i = 1 , 2 , , n , we have the vector of dimension q × 1
y i = ( y i 1 , y i 2 , , y i q ) .
Assume also that there are p explanatory variables X 1 , X 2 , , X p where, for i = 1 , 2 , , n ,
X i = ( x i 1 , x i 2 , , x i q ) ,
and there is a matrix q × p associated to the ith observed response y i with
x i j = ( x i j 1 , x i j 2 , , x i j p ) ,
for j = 1 , 2 , q , a p dimensional vector of values of the explanatory variables. For the vector of response variables, we use the operator vec ( · ) , which transforms matrices into a column vector from the columns of the matrix. Therefore,
y = vec ( y 1 , y 2 , , y n ) .
Thus, the MVSUSHN regression model is given by
Z i = X i β + ε i , i = 1 , 2 , , n ,
where z i j = log ( log ( 1 y i j ) ) ) , being β is a vector of unknown parameters of dimension p, and the vectors ε i for i = 1 , 2 , , n are vectors of independent and identically distributed random variables such that ε i MVSHN ( α , 0 q , Σ , λ ) where Σ = diag ( σ 1 , σ 2 , , σ q ) , and 0 q , a vector of zeros of dimension q . It follows that, z i SMVSHN ( α , X i β , Σ , λ ) .
From this result, we have by the theorem shown above that y i j USHN ( α j , x i j β j , σ j ) for j = 1 , 2 , , q i.e., each marginal follows a USHN regression model. Thus, defining X = diag ( X 1 , X 2 , , X q ) , where X j for j = 1 , 2 , , q is a matrix of size n × p j and X of dimension n q × ( p 1 + p 2 + + p q ) , then the MVSUSHN can be represented by
Z = X β + ε ,
where Z i j = X i j β j + ε i j , i = 1 , 2 , , n , j = 1 , 2 , , q taking values z i j = log ( log ( 1 y i j ) ) ) , with β = vec ( β 1 , β 2 , , β q ) with β j a vector of dimension p j × 1 , i.e., β is a vector of dimension p × 1 being p = p 1 + p 2 + + p q , ε = vec ( ε 1 , ε 2 , , ε q ) is an error vector with ε j MVSHN ( α , 0 n q , Σ , λ ) where ε j = ( ε 1 j , ε 2 j , , ε n j ) with ε i j SHN ( α j , 0 , σ j ) it follows that z i j SHN ( α j , X i j β j , σ j ) for j = 1 , 2 , , q .

Statistical Inference

For i = 1 , 2 , , n and j = 1 , 2 , , q , define: δ 1 = ( δ 11 , δ 21 , , δ n 1 ) ,   δ 2 = ( δ 12 , δ 22 , , δ n 2 ) and δ 3 = ( δ 13 , δ 23 , , δ n 3 ) ; with δ i 1 = ( δ i 11 , δ i 12 , , δ i 1 q ) ,   δ i 2 = ( δ i 21 , δ i 22 , , δ i 2 q ) and δ i 3 = ( δ i 31 , δ i 32 , , δ i 3 q ) where
δ i 1 j = 2 α j σ j ( 1 y i j ) ( log ( 1 y i j ) ) cosh z i j x i j β σ j , δ i 2 j = 2 α j sinh z i j x i j β σ j
and
δ i 3 j = Φ λ j = 1 p δ i 2 j .
for j = 1 , 2 , , p . The log-likelihood function for the parameter vector θ = ( α , β , Σ , λ ) is
( θ ) = i = 1 n i ( θ ) ,
where
i ( θ ) = q 2 log ( 2 π ) + i = 1 q σ j + i = 1 q log ( ( 1 y i j ) ( log ( 1 y i j ) ) ) + i = 1 q log ( δ i 1 j ) 1 2 i = 1 q δ i 2 j 2 + log ( δ i 3 j ) .
To obtain the score function, denoted U ( θ ) , we took the derivative of the log-likelihood function with respect to each of the parameters, so the elements of the score function are given by
U ( β j k ) = 1 σ j i = 1 n x i j k δ i 1 j δ i 2 j δ i 2 j δ i 1 j + λ σ j i = 1 n x i j k δ i 1 j ω i j j δ i 2 j , k = 1 , , p j ,
U ( α j ) = n α j + 1 α j i = 1 n δ i 2 j 2 λ α j i = 1 n ω i j = 1 q δ i 2 j ,
U ( σ j ) = n σ j 1 σ j i = 1 n v i j tanh ( v i j ) + 1 σ j i = 1 n v i j δ i 1 j δ i 2 j λ σ j i = 1 n δ i 1 j ω i j j δ i 2 j ,
U ( λ ) = i = 1 n ω i j = 1 q δ i 2 j ,
where v i j = ( log ( log ( 1 y i j ) x i j β j ) / σ j , for i = 1 , , n , and j = 1 , 2 , p and; w i = ϕ λ j = 1 q δ i 2 j / Φ λ j = 1 q δ i 2 j .
The maximum likelihood estimator (MLE) for β j 1 , , β j p j , α j and σ j are the solutions to the equations U ( β j k ) = 0 , U ( α j ) = 0 and U ( σ j ) = 0 for j = 1 , 2 , , q , k = 1 , 2 , , p j , which require numerical procedures. To start the iterative process, the least squares estimates can be used for the β j , i.e., β ˜ j 0 = ( X j X j ) 1 X j z j , from where σ ^ j 0 2 = 1 n p j i = 1 n ( z i j x j β ˜ j 0 ) 2 , while for the α j , the initial values could be implemented α ^ j 0 = α ˜ j , where α ˜ j = 4 n i = 1 n sinh z i j x j β ˜ j 0 σ ¯ j 0 . . With these initial values, and assuming them to be the true values of the parameters, a one-dimensional function for the parameter λ can be obtained, which can be estimated by some numerical method such as uniroot from R Development Core Team [29].
The elements of the observed information matrix that are calculated as minus the second derivative of the log-likelihood function with respect to the parameters, denoted by I θ j θ k , are given by
I θ j θ k = ( θ ) θ j θ k
The explicit expressions of these elements are presented in the Appendix B. The elements of the Fisher information matrix (denoted κ θ j θ k ) are given by the expectation of the elements of the observed information matrix, that is κ θ j θ k = E ( I θ j θ k ) , these are calculated numerically, therefore, the information matrix is expressed by
κ θ θ = E ( θ ) θ j θ k = κ α α κ α β κ α Σ κ α λ κ α β T κ β β κ β Σ κ β λ κ α Σ T κ β Σ T κ Σ Σ κ Σ λ κ α λ T κ β λ T κ Σ λ T κ λ λ
When λ = 0 , the case of the independence of univariate USHN distribution is obtained; thus, it follows that: κ β β = bloq . diag ( c ( α 1 ) x 1 x 1 / 4 , c ( α 2 ) x 2 x 2 / 4 , , c ( α q ) x q x q / 4 ) , is a diagonal block matrix where c ( α j ) = 1 + 4 α j 2 2 π / α j 2 1 erf [ ( 2 / α j 2 ) 1 / 2 ] exp ( 2 / α j 2 ) and erf ( x ) is the error function given by: erf ( x ) = 2 π 0 x e t 2 d t , see Rieck and Nedelman [8]; κ α α = diag ( 2 / α 1 2 , 2 / α 2 2 , , 2 / α q 2 ) ,   κ Σ Σ = diag ( κ σ 1 σ 1 , κ σ 2 σ 2 , , κ σ q σ q ) , where κ σ j σ j = a 2 ( α j , σ j ) σ j 2 + 2 b ( α j , σ j ) d ( α j , σ j ) σ j 2 with a l ( α j , σ j ) = E v j l 2 δ 2 j 2 + 4 α j 2 1 + δ 2 j 2 δ 2 j 2 + 4 / α j 2 , b ( α j , σ j ) = E ( v j δ 1 j δ 2 j ) and d ( α j , σ j ) = E v j δ 1 j δ 2 j . Expectations in the above expressions must be calculated numerically, κ α β = 0 is a matrix of zeros, κ α Σ = diag ( κ α 1 σ 1 , κ α 2 σ 2 , , κ α q σ q ) with κ α j σ j = 2 α j σ j b ( α j , σ j ) , κ β Σ = diag ( κ β 1 σ 1 , κ β 2 σ 2 , , κ β q σ q ) with κ β j σ j = a 1 ( α j , σ j ) 2 σ j x j , κ α λ = 0 is a vector of size q ,   κ β λ = 0 is a vector of size p 1 + p 2 + + p q   κ Σ λ = 0 is a vector of size q and κ λ λ = 2 π .
The rows (or columns) of the matrix κ θ θ are linearly independent, so the determinant is different than zero; this guarantees the existence of the inverse of κ θ θ . Hence, for large samples, the MLE θ ^ of θ is asymptotically normal, that is,
θ ^ D N p + 2 q + 1 ( θ , κ θ θ 1 ) ,
resulting that the asymptotic variance of the MLE θ ^ is the inverse of I ( θ ^ ) . The approximation to the N p + 2 q + 1 ( θ , κ θ θ 1 ) can be used to construct confidence intervals for α j , β j k j , σ j y λ , these are given by α j ^ z 1 ξ / 2 κ ^ ( α ^ ) , β ^ j k j z 1 ξ / 2 κ ^ ( β ^ j k j ) and λ ^ z 1 ξ / 2 κ ^ ( λ ^ ) , where κ ^ ( · ) is on the diagonal of the matrix κ θ θ 1 for each parameter and z 1 ξ / 2 is the quantil 100 ( δ / 2 ) % of the standard normal distribution.
The hypothesis of interest
H 0 : λ = 0 versus λ 0
can be tested by the statistic test (for large n)
2 [ ^ ( θ ^ * , λ ^ ) ^ ( θ ^ * , 0 ) ] χ 1 2
where θ * is the vector θ without the parameter λ , and θ ^ * is the MLE of θ restricted on H 0 .

5. Numerical Results

In this section, the results of the simulation study and two applications to illustrate the applicability of the proposed models are presented.

5.1. Simulation Study

In order to study the behavior of the MLE of the parameter vector of the MVSUSHN regression model, a small Monte Carlo simulation study with covariates in the model was carried out. We considered X 1 i N ( 0 , 1 ) and X 2 i U ( 0 , 1 ) for i = 1 , 2 , , n . For λ = 0 , the case of independence; the results for each model can be seen in the studies conducted by Martínez-Flórez et al. [3].
We used the bivariate MVSUSHN ( ( α 1 , α 2 ) , ( β 10 , β 11 , β 20 , β 21 ) , ( σ 1 , σ 2 ) , λ ) regression model, and we took the values: α 1 = 1.5 , α 2 = 0.75 , β 10 = 0.75 , β 11 = 0.50 , β 20 = 0.50 , and β 21 = 1.5 and σ 1 = σ 2 = 2 , while λ = 1.75 , 3.5 , and 5.25 . The sample size was n = 40 , 80 , 120 , and 200, and the number of iterations was 5000. We studied the absolute value of the relative bias (RB), root of mean square error (RMSE), length of confidence interval (LCI), and coverage probability (CP). To generate the samples, we performed the following algorithm:
(1)
Generate a uniform random U 1 U ( 0 , 1 ) and a random number x 1 with distribution N ( 0 , 1 ) .
(2)
Generate ε 1 = 2 arcsinh ( α 1 Φ 1 ( U 1 ) / 2 ) with Φ 1 ( · ) the inverse of the standard normal function.
(3)
Let y 1 = 1 exp ( exp ( β 10 + β 11 x 1 + ε 1 ) ) .
(4)
Compute b = ( 2 / α 1 ) sinh ( ( log ( log ( 1 y 1 ) ) ( β 10 + β 11 x 1 ) ) / 2 ) .
(5)
Generate another uniform random number (independent of U 1 ) U 2 U ( 0 , 1 ) and x 2 also with distribution U ( 0 , 1 )
(6)
Compute the error ε 2 such that ε 2 = 2 arcsinh ( α 2 Φ SN 1 ( U 2 , 0 , 1 , λ b 1 ) / 2 ) , where Φ SN 1 ( · , 0 , 1 , · ) is the inverse function of the standard skew-normal and arcsinh ( · ) is the inverse of the hyperbolic sine function.
(7)
Let y 2 = 1 exp ( exp ( β 20 + β 21 x 2 + ε 2 ) ) . This algorithm is generated n times, finally obtaining the USHN bivariate random sample.
The results of the simulations can be seen in the Table 2. In general, it can be seen that the RB, RSME, and LCI of the model parameters decrease as n increases; this decrease is slower for some parameters. It can also be seen that the CP increases as n increases.

5.2. Illustration 1

To show the relevance of the MVSUSHN distributio, a real data set from a study conducted by Freeman [30] on drunk driving legislation and traffic fatalities in 48 states in the United States of America (USA) during the period from 1980 to 2004 is considered. The database is available in the wooldridge library by Shea and Brown [31] of the software R Development Core Team [29] under the name of driving, and it contains information associated with current legislation, accident records, and some demographic characteristics. For this illustration, the unemployment rate variables ( y 1 ) and the percent of the population aged 14 through 24 ( y 2 ) were used. The bivariate beta (BVBeta) model of Cepeda et al. [25], the bivariate Johnson SB (BVJSB) model of Lemonte and Moreno-Arenas [27], and the BVSUSHN model were fitted. The BVBeta distribution of Cepeda et al. [25] is based on the e Farlie–Gumbel–Morgentern copula; see Nelsen [32], and their joint pdf is given by
f X 1 X 2 ( x 1 , x 2 ) = f X 1 ( x 1 ) f X 2 ( x 2 ) ( 1 + θ ( 1 2 F X 1 ( x 1 ) ) ( 1 2 F X 2 ( x 2 ) ) ,
where f X j ( x j ) and F X j ( x j ) correspond, respectively, to the pdf and cdf of the beta distributions of parameters α j and β j for j = 1 , 2 and θ ( 1 , 1 ) . To compare models, we used the Akaike information criterion (AIC) of [33] and the Bayesian information criterion (BIC) of [34], defined, respectively, by
A I C = 2 ( θ ^ ) + 2 p and B I C = 2 ( θ ^ ) + n log ( n ) ,
where p is the number of parameters and ( · ) is the log-likelihood function evaluated at the MLEs of parameters. The best model is the one with the smallest AIC or BIC.
For fitting the bivariate model, we used the optim function of R Development Core Team [29]. The parameter estimates of these models, accompanied by their standard errors in parentheses, obtained using the maximum likelihood method, are given in the Table 3. According to the AIC and BIC criteria, the BVSUSHN model presents the best fit.
The graphs in the Figure 3 show the contours of the fitted models. For the BVSUSHN distribution, we have that X j = log ( log ( 1 Y j ) ) USHN ( α j , γ j , σ j ) for j = 1 , 2 and ( X 1 , X 2 ) BVSUSHN ( α 1 , α 2 , γ 1 , γ 2 , σ 1 , σ 2 , λ ) , then it follows that
W j = 2 α j sinh X j γ j σ j N ( 0 , 1 ) , j = 1 , 2
and
( W 1 , W 2 ) SBSN ( 0 2 , I 2 , λ )
where 0 2 and I 2 are a column vector of size 2 and an identity matrix of size 2 × 2 , respectively. The statistic of the Kolmogorov–Smirnov (KS) goodness-of-fit test joint with the respective p -values for the marginal distributions for the BVSJB, BVBeta, and BVSUSHN distribution are presented in Table 3. From here, it can be see that the BVSUSHN distribution shows a good fit compared with the BVSJB and BVBeta models.

5.3. Illustration 2

In the second illustration, the USHN bivariate regression model is fitted. The real data were taken from http://www.pe.undp.org (accessed on 31 August 2022), and they correspond to measurements made on 195 districts in Peru. In this illustration, the interest is to model the human development index ( Y 1 ) and illiteracy rate ( Y 2 ) as functions of the proportion of people with high poverty level (HPL) by using the MVSUSHN regression model. As far as poverty is concerned, poor are identified as those unable to obtain minimum required calorie per day to keep body and soul together. We have that log ( log ( 1 Y 1 i ) ) = β 10 + β 11 H P L + ε 1 i and log ( log ( 1 Y 2 i ) ) = β 20 + β 21 H P L + ε 2 i with ( ε 1 , ϵ 2 ) BVSHN ( ( α 1 , α 2 ) , ( 0 , 0 ) , diag ( σ 1 , σ 2 ) , λ ) .
The MLEs of the vector of model parameters, with standard errors in parenthesis are given by α ^ 1 = 0.3938 ( 0.1539 ) ,   β ^ 10 = 0.0269 ( 0.0075 ) ,   β ^ 11 = 0.5734 ( 0.0277 ) ,   σ ^ 1 = 0.3765 ( 0.1384 ) ,   α ^ 2 = 0.1303 ( 0.0654 ) ,   β ^ 20 = 2.8999 ( 0.0500 ) ,   β ^ 21 = 3.5010 ( 0.1832 ) ,   σ ^ 2 = 7.8736 ( 3.9272 ) , and λ ^ = 3.1520 ( 0.6167 ) .
W j = 2 α j sinh log ( log ( 1 Y j i ) ) β j 0 β j 1 H P L σ j N ( 0 , 1 ) , j = 1 , 2 ,
and
( W 1 , W 2 ) SBSN ( 0 2 , I 2 , λ ) ,
where 0 2 and I 2 are a column vector of size 2 and an identity matrix of size 2 × 2 . To study the model fit, we perform the Kolmogorov–Smirnov test for the bivariate vector ( ε 1 , ε 2 ) .
For the multivariate Kolmogorov–Smirnov test of goodness of fit proposed by Justel [35], special for the case of a bivariate distribution, which we denote by BKS (bivariate Kolmogorov–Smirnov), the statistic is given by
d n = sup ( x 1 , x 2 ) R 2 F n ( x 1 , x 2 ) F ( x 1 , x 2 )
where F n is the empirical distribution function of the sample, and F is some specified distribution function. When the F distribution is unknown, the Kolmogorov–Smirnov statistic is defined by
d n ( F ) = max D 1 , D 2 ,
where
D 1 = sup ( x 1 , x 2 ) R 2 G n ( y 1 , y 2 ) y 1 × y 2
by using the transformations y 1 = F X 1 ( x 1 ) , y 2 = F X 2 X 1 ( x 2 x 1 ) , and
D 2 = sup ( x 1 , x 2 ) R 2 G n ( y 2 , y 1 ) y 2 × y 1
by using the transformations y 2 = F X 2 ( x 2 ) and y 1 = F X 1 X 2 ( x 1 x 2 ) , where G n is the empirical distribution function of the sample. For the special case of the BVSUSHN model,
d n ( BPSN ) = max 0.08622004 , 0.0750111 = 0.08622004 ,
which is less than 0.1265 (for n = 200 ), which is the critical value given by Justel [35] at level of 5%. Therefore, it is concluded that the MVSUSHN model fits the data set well.
We also performed univariate Kolmogorov–Smirnov goodness-of-fit tests for W j with j = 1 , 2 yielding the test statistics of D 1 = 0.5282 with p value = 0.949 and D 2 = 0.5076 with p value = 0.9898 , indicating that the marginals show a good fit.
The Figure 4 shows the envelope plots for the marginal and contour distributions for the residuals of the fitted model. For the envelope plot, we used the martingale residual transformation, r M T i , proposed by Barros et al. [36]. These residuals are defined by
r M T i = sgn ( r M i ) 2 [ r M i + α i log ( δ i r M i ) ] ; i = 1 , 2 , , n
where r M i = δ i + log ( S ( e i , θ ^ ) ) is the martingale residual proposed by Ortega et al. [37], where δ i = 0 , 1 indicates whether the i-th observation is censored or not, respectively, sgn ( r M i ) denotes the sign of r M i , and S ( e i ; θ ^ ) represents the survival function evaluated at e i , where θ ^ are the MLE for θ .

6. Concluding Remarks

Diverse distributions to deal with the problem of bounded data on the interval ( 0 , 1 ) were proposed, with great applicability in all fields of knowledge, especially in the social sciences, humanities, medicine, and engineering, among others. However, few proposals have been developed in the statistical literature to jointly model two or more variables, such as those described above, and especially that incorporate covariates to explain the variability of the variables of interest.
In this paper, a new multivariate distribution was introduced from the conditionally specified distributions methodology useful for modeling responses in the ( 0 , 1 ) p region jointly. The new distribution, which is absolutely continuous, is called the skewed log-Birnbaum–Saunders distribution and is also extended to the case of regression models. For the multivariate distribution, the marginal densities and conditional distributions were presented, and the Fisher information matrix of the multivariate regression model was also presented. For the estimation of the parameters in the models, a classical approach was used together with the maximum likelihood method. A small Monte Carlo simulation study was carried out to study the benefits and limitations of the new methodologies, which allows us to conclude that the parameter estimators behave asymptotically well. Two applications with real data to illustrate the usefulness of the introduced methodologies showed great flexibility to model data in ( 0 , 1 ) p for the particular case p = 2 , which makes them excellent alternatives to existing methodologies in the statistical literature.

Author Contributions

Conceptualization, G.M.-F., S.V.-C., and R.T.-F.; data curation, G.M.-F., S.V.-C., and R.T.-F.; formal analysis, G.M.-F., S.V.-C., R.T.-F., and L.R.-Q.; funding acquisition, G.M.-F., and R.T.-F.; investigation, G.M.-F., S.V.-C., R.T.-F., and L.R.-Q.; resources, G.M.-F., and R.T.-F.; software, G.M.-F., and R.T.-F.; supervision, G.M.-F., and R.T.-F.; validation, G.M.-F., S.V.-C., R.T.-F., and L.R.-Q.; visualization, G.M.-F., S.V.-C., R.T.-F., and L.R.-Q.; writing—original draft, G.M.-F., S.V.-C., and R.T.-F.; writing—review and editing, G.M.-F., S.V.-C., and R.T.-F. All authors have read and agreed to the published version of the manuscript.

Funding

The research of G. Martínez-Flórez and R. Tovar-Falón was supported by project: Resolución de Problemas de Situaciones Reales Usando Análisis Estadístico a través del Modelamiento Multidimensional de Tasas y Proporciones; Esquemas de Monitoreamiento para Datos Asimétricos no Normales y una Estrategia Didáctica para el Desarrollo del Pensamiento Lógico-Matemático. Universidad de Córdoba, Colombia, Acta de Compromiso Número FCB-05-19.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Details about available data are given in Section 5.

Acknowledgments

G. Martínez-Flórez and R. Tovar-Falón acknowledge the support given by Universidad de Córdoba, Montería, Colombia. S. Vergara-Cardozo recognizes the support given by Universidad Nacional de Colombia, Sede Bogotá.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A. Expected Value of the LSHN Distribution

This Appendix presents the derivation of the expected value of a random variable with LSHN distribution, which is used to obtain the Equation (18).
Let Z SHN ( α , γ σ ) , then X = exp ( Z ) LSHN ( α , γ σ ) , where LSHN denotes the pdf of a non-negative SHN distribution; see Martínez-Flórez et al. [3].
Thus, if X LSHN , then Z = log ( X ) SHN . Following Rieck and Nedelman [8], we have that
X = exp ( Z ) X r = exp ( r Z ) E ( X r ) = E ( e r Z ) E ( X r ) = k = 0 ( r ) k k ! E ( Z k ) E ( X r ) = k = 0 ( r ) k k ! e k r k a ( α 2 ) + k b ( α 2 ) k 1 / 2 ( α 2 )
where a = ( r σ + 1 ) / 2 , b = ( r σ 1 ) / 2 and k Λ ( · ) is the third-order function of Bessel.
Now, if X LSHN ( α , γ σ ) , then Y = 1 exp ( X ) USHN ( α , γ , σ ) .
Hence,
E ( Y n ) = E [ ( 1 e X ) n ] = j = 0 n n j ( 1 ) j E ( e j X ) = j = 0 n l = 0 n j ( 1 ) j + 1 ( j e γ ) l l ! k a 1 ( α 2 ) + k b 1 ( α 2 ) k 1 / 2 ( α 2 )
The last term is obtained by using Taylor expansion for e j X for the jth moment of the LSHN ( α , γ , σ ) .

Appendix B. Elements of the Observed Information for the SMVSHN Regression Model

This Appendix presents the elements of the observed information, which are calculated from Equation (26).
I α j α j = n α j 2 + 3 α j 2 i = 1 n δ i 2 j 2 + λ α j 2 i = 1 n j = 1 q δ i 2 j ω i 2 α j + λ j = 1 q δ i 2 j λ j = 1 q δ i 2 j + ω i ,
I α j α j = λ α j α j i = 1 n j = 1 q δ i 2 j ω i 1 + λ j = 1 q δ i 2 j λ j = 1 q δ i 2 j + ω i ,
I α j β j k = 2 α j σ j i = 1 n x i j k δ i 1 j δ i 2 j + λ α j σ j i = 1 n x i j k δ i 1 j j j δ i 2 j ω i 1 + λ j = 1 q δ i 2 j λ δ i 2 j j j δ i 2 j + ω i ,
I α j β j k = λ α j σ j i = 1 n x i j k δ i 1 j j j δ i 2 j ω i 1 + λ j = 1 q δ i 2 j λ δ i 2 j j j δ i 2 j + ω i ,
I α j σ j = 2 α j σ j i = 1 n v i j δ i 1 j δ i 2 j + λ α j σ j i = 1 n v i j δ i 1 j j j δ i 2 j ω i 1 + λ λ δ i 2 j j j δ i 2 j + ω i ,
I α j σ j = λ α j σ j i = 1 n v i j δ i 1 j j j δ i 2 j ω i 1 + λ λ δ i 2 j j j δ i 2 j + ω i ,
I α j λ = 1 α j i = 1 n j = 1 q δ i 2 j ω i 1 λ j = 1 q δ i 2 j λ j = 1 q δ i 2 j + ω i ,
I β j k β j k = 1 σ j 2 i = 1 n x i j k x i j k 2 δ i 2 j 2 + 4 α j 2 1 + δ i 2 j 2 δ i 2 j 2 + 4 / α j 2 + λ σ j 2 i = 1 n x i j k x i j k j j δ i 2 j ω i δ i 2 j + λ δ i 1 j 2 j j δ i 2 j λ δ i 2 j j j δ i 2 j + ω i ,
I β j k β j k = λ σ j σ j i = 1 n x i j k x i j k δ i 1 j δ i 1 j ω i l j , j δ i 2 l + λ j j δ i 2 j j j δ i 2 j λ δ i 2 j j j δ i 2 j + ω i
I β j k σ j = 1 2 σ j i = 1 n x i j k v i j δ i 1 j 2 + δ i 2 j 2 s e c h 2 v i j λ σ j 2 i = 1 n x i j k ω i δ i 1 j + δ i 2 j v i j j j δ i 2 j + λ 2 σ j 2 i = 1 n x i j k v i j δ i 2 j 2 ω i j j δ i 2 j 2 λ δ i 2 j j j δ i 2 j + ω i
I β j k σ j = λ σ j σ j i = 1 n x i j k v i j δ i 1 j δ i 1 j ω i l j , j δ i 2 l + λ j j δ i 2 j j j δ i 2 j λ δ i 2 j j j δ i 2 j + ω i
I β j k λ = 1 σ j i = 1 n x i j k δ i 1 j j j δ i 2 j ω i 1 λ j = 1 q δ i 2 j λ j = 1 q δ i 2 j + ω i ,
I σ j σ j = 2 σ j 2 i = 1 n v i j δ i 1 j δ i 2 j δ i 2 j δ i 1 j + 1 σ j 2 i = 1 n v i j 2 2 δ i 2 j 2 + 4 α j 2 1 + δ i 2 j 2 δ i 2 j 2 + 4 / α j 2 + λ σ j 2 i = 1 n j j δ i 2 j ω i 2 v i j δ i 1 j v i j 2 δ i 2 j + λ v i j 2 δ i 1 j 2 j j δ i 2 j λ δ i 2 j j j δ i 2 j + ω i
I σ j σ j = λ σ j σ j i = 1 n v i j v i j δ i 1 j δ i 1 j ω i l j , j δ i 2 l + λ j j δ i 2 j j j δ i 2 j λ δ i 2 j j j δ i 2 j + ω i
I σ j λ = 1 σ j i = 1 n v i j δ i 1 j j j δ i 2 j ω i 1 λ j = 1 q δ i 2 j λ j = 1 q δ i 2 j + ω i ,
I λ λ = i = 1 n j = 1 q δ i 2 j 2 ω i λ j = 1 q δ i 2 j + ω i ,

References

  1. Ferrari, S.; Cribari-Neto, F. Beta regression for modelling rates and proportions. J. Appl. Stat. 2004, 31, 799–815. [Google Scholar] [CrossRef]
  2. Kumaraswamy, P. A generalized probability density function for double-bounded random processes. J. Hydrol. 1980, 46, 79–88. [Google Scholar] [CrossRef]
  3. Martínez-Flórez, G.; Tovar-Falón, R. Regression Models Based on the Unit Sinh-Normal Distribution. Mathematics 2021, 9, 1231. [Google Scholar] [CrossRef]
  4. Martínez-Flórez, G.; Azevedo-Farias, R.B.; Tovar-Falón, R. New Class of Unit-Power-Skew-Normal Distribution and Its Associated Regression Model for Bounded Responses. Mathematics 2022, 10, 3035. [Google Scholar] [CrossRef]
  5. Kieschnick, R.; Mccullough, B.D. Regression analysis of variates observed on (0,1). Stat. Model. 2003, 3, 193–213. [Google Scholar] [CrossRef] [Green Version]
  6. Mazucheli, J.; Menezes, A.F.B.; Ghitany, M.E. The unit-Weibull distribution and associated inference. J. Appl. Probab. Stat. 2018, 13, 1–22. [Google Scholar]
  7. Menezes, A.F.B.; Mazucheli, J.; Dey, S. The unit-logistic distribution: Different methods of estimation. Pesqui. Oper. 2018, 38, 555–578. [Google Scholar] [CrossRef]
  8. Rieck, J.R.; Nedelman, J.R. A log-linear model for the Birnbaum-Saunders distribution. Technometrics 1991, 33, 51–60. [Google Scholar]
  9. Martínez-Flórez, G.; Elal-Olivero, D.; Barrera-Causil, C. Extended Generalized Sinh-Normal Distribution. Mathematics 2021, 9, 2793. [Google Scholar] [CrossRef]
  10. Lemonte, A.J. A log-Birnbaum-Saunders regression model with asymmetric errors. J. Stat. Comput. Simul. 2011, 82, 1775–1787. [Google Scholar] [CrossRef] [Green Version]
  11. Leiva, V.; Vilca-Labra, F.; Balakrishnan, N.; Sanhueza, A. A skewed sinh-normal distribution and its properties and application to air pollution. Commun. Stat. Theory Methods 2010, 39, 426–443. [Google Scholar] [CrossRef]
  12. Moreno-Arenas, G.; Martínez-Flórez, G.; Barrera-Causil, C. Proportional Hazard Birnbaum-Saunders Distribution with Application to the Survival Data Analysis. Rev. Colomb. Estad. 2016, 39, 129–147. [Google Scholar] [CrossRef]
  13. Birnbaum, Z.W.; Saunders, S.C. A new family of life distributions. J. Appl. Probab. 1969, 6, 319–327. [Google Scholar] [CrossRef]
  14. Martínez-Flórez, G.; Bolfarine, H.; Gómez, H.W. The Log-Linear Birnbaum-Saunders Power Model. Methodol. Comput. Appl. Probab. 2017, 19, 913–933. [Google Scholar] [CrossRef]
  15. Díaz–García, J.A.; Domínguez–Molina, J.R. Some generalisations of Birnbaum–Saunders and sinh–normal distributions. Int. Math. Forum 2014, 1, 1709–1727. [Google Scholar]
  16. Lemonte, A. A Multivariate Birnbaum-Saunders regression model. J. Stat. Comput. Simul. 2013, 12, 2244–2257. [Google Scholar] [CrossRef]
  17. Marchant, C.; Leiva, V.; Cysneiros, F. A multivariate log-linear model for Birnbaum-Saunders distributions. IEEE Trans. Reliab. 2016, 65, 816–864. [Google Scholar] [CrossRef]
  18. Martínez-Flórez, G.; Azevedo-Farias, R.B.; Tovar-Falón, R. An exponentiated multivariate extension for the Birnbaum-Saunders log-linear model. Mathematics 2022, 10, 1299. [Google Scholar] [CrossRef]
  19. Mazucheli, J.; Menezes, A.; Dey, S. The unit-Birnbaum-Saunders distribution with applications. Chil. J. Stat. 2018, 9, 47–57. [Google Scholar]
  20. Arnold, B.C.; Castillo, E.; Sarabia, J.M. Conditionally specified multivariate skewed distributions. Sankhya A 2002, 64, 206–226. [Google Scholar]
  21. Arnold, B.C.; Castillo, E.; Sarabia, J.M. Conditionally specified distributions. In Lecture Notes in Statistics; Berger, J., Fienberg, J., Gani, J., Krickeberg, I., Singer, B., Eds.; Springer: New York, NY, USA, 1992; Volume 73. [Google Scholar]
  22. Azzalini, A. A class of distributions which includes the normal ones. Scand. J. Stat. 1985, 12, 171–178. [Google Scholar]
  23. Lemonte, A.J.; Martínez-Flórez, G.; Moreno-Arenas, G. Multivariate Birnbaum-Saunders distribution: Properties and associated inference. J. Stat. Comput. Simul. 2015, 85, 374–392. [Google Scholar] [CrossRef]
  24. Martínez-Flórez, G.; Azevedo-Farias, R.; Moreno-Arenas, G. Multivariate log-Birnbaum-Saunders regression models. Commun. Stat. Theory Methods 2017, 46, 10166–10178. [Google Scholar] [CrossRef]
  25. Cepeda-Cuervo, E.; Achcar, J.A.; Garrido-Lopera, L. Bivariate beta regression models: Joint modeling of the mean, dispersion and association parameters. J. Appl. Stat. 2014, 41, 677–687. [Google Scholar] [CrossRef]
  26. Souza, D.F.; Moura, F.A.S. Multivariate Beta Regression with Application in Small Area Estimation. J. Off. Stat. 2016, 32, 747–768. [Google Scholar] [CrossRef] [Green Version]
  27. Lemonte, A.J.; Moreno-Arenas, G. On a multivariate regression model for rates and proportions. J. Appl. Stat. 2019, 46, 1084–1106. [Google Scholar] [CrossRef]
  28. Tables for computing bi-variate normal probabilities. Ann. Math. Stat. 1976, 27, 1075–1090.
  29. R Development Core Team. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2021; Available online: http://www.R-project.org (accessed on 31 August 2022).
  30. Freeman, D.G. Drunk Driving Legislation and Traffic Fatalities: New Evidence on BAC 08 Laws. Contemp. Econ. Policy 2007, 25, 293–308. [Google Scholar] [CrossRef]
  31. Shea, J.M.; Brown, K.J. Wooldridge: 115 Data Sets. In Introductory Econometrics: A Modern Approach, 7e; Wooldridge, J.M., Ed.; R package version 1.4-2; Cengage Learning: Boston, MA, USA, 2022. [Google Scholar]
  32. Nelsen, R.B. An Introduction to Copulas; Springer Series in Statistics; Springer: New York, NY, USA, 2006; p. 272. [Google Scholar]
  33. Akaike, H. A new look at statistical model identification. IEEE Trans. Autom. Control. 1974, AU-19, 716–722. [Google Scholar] [CrossRef]
  34. Schwarz, G. Estimating the dimension of a model. Ann. Stat. 1978, 6, 461–464. [Google Scholar] [CrossRef]
  35. Justel, A.; Peña, D.; Zamar, Z. A multivariate Kolmogorov-Smirnov test of goodness of fit. Stat. Probab. Lett. 1997, 35, 251–259. [Google Scholar] [CrossRef] [Green Version]
  36. Barros, M.; Galea, M.; Gonzalez, M.; Leiva, V. Influence diagnostics in the tobit censored response model. Stat. Methods Appl. 2010, 19, 379–397. [Google Scholar] [CrossRef]
  37. Ortega, E.M.; Bolfarine, H.; Paula, G.A. Influence diagnostics in generalized log-gamma regression models. Comput. Stat. Data Anal. 2003, 42, 165–186. [Google Scholar] [CrossRef]
Figure 1. Graphs of contours for the BVSUSHN distribution for: (a) BVSUSHN ( 0.5 , 1.5 , 0 , 0 , 1 , 1 , 1.5 ) , (b) BVSUSHN ( 0.25 , 0.75 , 0 , 0 , 1 , 1 , 0.75 ) , and (c) BVSUSHN ( 1.25 , 1.75 , 0 , 0 , 1 , 1 , 0.75 ) .
Figure 1. Graphs of contours for the BVSUSHN distribution for: (a) BVSUSHN ( 0.5 , 1.5 , 0 , 0 , 1 , 1 , 1.5 ) , (b) BVSUSHN ( 0.25 , 0.75 , 0 , 0 , 1 , 1 , 0.75 ) , and (c) BVSUSHN ( 1.25 , 1.75 , 0 , 0 , 1 , 1 , 0.75 ) .
Mathematics 11 01095 g001
Figure 2. Graphs of density function of the BVSUSHN distribution for: (a) BVSUSHN ( 0.5 , 1.5 , 0.0 , 0.0 , 1.0 , 1.0 , 1.5 ) , (b) BVSUSHN ( 0.25 , 0.75 , 0.0 , 0.0 , 1.5 , 1.5 , 0.75 ) , (c) BVSUSHN ( 1.25 , 1.75 , 0.0 , 0.0 , 1.0 , 1.0 , 0.75 ) , and (d) BVSUSHN ( 2.5 , 2.5 , 0.0 , 0.0 , 1.0 , 1.0 , 2.5 ) .
Figure 2. Graphs of density function of the BVSUSHN distribution for: (a) BVSUSHN ( 0.5 , 1.5 , 0.0 , 0.0 , 1.0 , 1.0 , 1.5 ) , (b) BVSUSHN ( 0.25 , 0.75 , 0.0 , 0.0 , 1.5 , 1.5 , 0.75 ) , (c) BVSUSHN ( 1.25 , 1.75 , 0.0 , 0.0 , 1.0 , 1.0 , 0.75 ) , and (d) BVSUSHN ( 2.5 , 2.5 , 0.0 , 0.0 , 1.0 , 1.0 , 2.5 ) .
Mathematics 11 01095 g002
Figure 3. Contour plots for the fitted models. (a) VBJSB, (b) VBBeta, and (c) SMVUSHN.
Figure 3. Contour plots for the fitted models. (a) VBJSB, (b) VBBeta, and (c) SMVUSHN.
Mathematics 11 01095 g003
Figure 4. Envelope plots for the marginals and contour for the fitted model. (a) W 1 , (b) W 2 , and (c) MVSUSHN.
Figure 4. Envelope plots for the marginals and contour for the fitted model. (a) W 1 , (b) W 2 , and (c) MVSUSHN.
Mathematics 11 01095 g004
Table 1. Correlation coefficient for BVUSHN distribution.
Table 1. Correlation coefficient for BVUSHN distribution.
λ α 1 / α 2 0.0850.50.751.52.253.05.07.510.0
0.0950.51240.34010.34100.34100.33880.33470.30770.28460.2759
0.450.35130.31600.31690.31730.31550.31170.27160.22390.1966
0.750.35240.31710.31810.31830.31650.31260.27220.22430.1969
0.51.00.35280.31770.31860.31880.31680.31280.27220.22430.1969
2.00.35070.31630.31710.31690.31480.31060.26970.22200.1949
3.00.34560.31190.31260.31220.30990.30570.26500.21800.1915
5.00.32610.27180.27220.27130.26900.26500.23180.19640.1775
7.50.31700.22400.22430.22340.22130.21800.19640.17740.1693
10.00.31880.19670.19690.19620.19430.19150.17750.16930.1684
0.0950.82580.59500.60030.61230.61770.61830.58090.52520.4904
0.450.61920.52690.53210.54420.55040.55170.50700.43430.3845
0.750.62600.53310.53830.55040.55650.55760.51200.43810.3877
1.00.63110.53790.54310.55510.56110.56210.51560.44080.3900
1.52.00.64250.54970.55480.56620.57170.57220.52320.44600.3941
3.00.64460.55270.55760.56850.57350.57360.52320.44510.3932
5.00.61400.50790.51200.52070.52410.52320.47740.40960.3661
7.50.56950.43490.43810.44440.44640.44510.40960.36100.3315
10.00.54340.38500.38770.39300.39450.39320.36610.33150.3117
0.0950.93170.67330.68120.70130.71310.71840.68780.62820.5857
0.450.69750.58750.59450.61280.62420.62940.59000.51630.4623
0.750.70730.59580.60300.62140.63270.63790.59760.52240.4677
1.00.71540.60280.60990.62840.63970.64480.60360.52730.4718
2.52.00.73750.62220.62940.64780.65880.66350.61940.53960.4820
3.00.74690.63070.63790.65590.66650.67080.62500.54320.4847
5.00.72430.59120.59760.61320.62190.62500.58230.50920.4579
7.50.67720.51720.52240.53480.54140.54320.50920.45340.4157
10.00.64300.46320.46770.47810.48340.48470.45790.41570.3882
0.0950.97650.70380.71280.73740.75310.76160.73670.68010.6363
0.450.72670.61350.62140.64280.65700.66450.62900.55670.5025
0.750.73790.62270.63070.65230.66660.67420.63790.56420.5090
1.00.74740.63070.63870.66050.67490.68240.64530.57050.5144
3.52.00.77550.65420.66240.68450.69880.70610.66620.58740.5288
3.00.78950.66590.67420.69620.71020.71720.67550.59450.5345
5.00.77300.63030.63790.65760.66980.67550.63670.56330.5098
7.50.72970.55790.56420.58050.59020.59450.56330.50640.4656
10.00.69510.50350.50900.52300.53110.53450.50980.46560.4348
0.0950.99380.72280.73240.76010.77900.79030.77070.71820.6763
0.450.74460.63230.64050.66430.68080.69030.65830.58810.5346
0.750.75630.64180.65020.67430.69100.70060.66810.59670.5421
1.00.76690.65050.65900.68340.70020.70990.67670.60400.5486
5.02.00.79920.67720.68600.71100.72810.73770.70210.62540.5672
3.00.81700.69170.70060.72580.74280.75230.71490.63580.5760
5.00.80610.65970.66810.69130.70660.71490.68030.60840.5544
7.50.76720.58940.59670.61640.62910.63580.60840.55210.5102
10.00.73500.53570.54210.55940.57030.57600.55440.51020.4776
Table 2. Relative bias, root of mean square error, length of confidence interval, and coverage probability for the MVSUSHN regression model.
Table 2. Relative bias, root of mean square error, length of confidence interval, and coverage probability for the MVSUSHN regression model.
λ = 1.75 λ = 3.50 λ = 5.25
n θ i RBRMSELCICPRBRMSELCICPRBRMSELCICP
40 α 1 0.22790.65294.39440.99920.18580.56324.23471.00000.15880.48974.16031.0000
α 2 0.57970.38823.95571.00000.41860.24543.79131.00000.27290.17523.77301.0000
β 10 0.16480.04460.64710.72580.10590.03010.55430.77230.08440.02260.50620.8373
β 11 0.06310.02950.62830.92230.04140.02180.54720.92720.03140.01890.50460.9346
β 20 0.02250.03130.60370.89490.02100.02280.52510.91030.01620.01660.48330.9279
β 21 0.17250.16811.15390.68160.10470.09560.96310.77420.08380.06850.88420.7535
σ 1 0.05460.45984.50550.89340.03080.44184.61310.92720.02840.39354.49960.9519
σ 2 0.33590.69079.37600.70980.26931.00026.52280.80360.34441.58706.21030.9096
λ 0.17050.75523.47340.84460.20844.61588.52850.75510.20876.244010.38600.7217
80 α 1 0.04220.18412.66641.00000.04240.17322.64721.00000.03960.11092.72761.0000
α 2 0.12960.07982.68921.00000.22380.06572.80551.00000.29740.07272.78641.0000
β 10 0.15800.02700.45980.82390.10480.01670.39870.86830.07840.00990.35910.8871
β 11 0.05490.01410.43740.92820.04030.00910.37400.94240.03060.00790.34400.9392
β 20 0.01730.01180.42930.92760.01520.00920.36830.93800.00510.00680.32690.9286
β 21 0.16980.10330.79480.76580.10400.05410.66370.83200.07970.03360.57860.8199
σ 1 0.03010.22653.03650.99730.02890.18682.94270.99930.02230.11922.75081.0000
σ 2 0.10440.48677.97240.96060.23630.84635.80321.00000.22651.27975.88601.0000
λ 0.10380.19311.88140.90430.17820.74324.01990.91700.16741.49696.68110.8465
120 α 1 0.03630.09032.29341.00000.03130.07932.19801.00000.02080.10482.07281.0000
α 2 0.11720.04132.68551.00000.17600.05942.31821.00000.18530.06132.40981.0000
β 10 0.15220.01960.38020.82480.10370.01120.32170.87550.07760.00900.28990.9029
β 11 0.05370.00780.35330.92940.03750.00620.30080.94310.03060.00510.27620.9444
β 20 0.00430.00920.36350.94080.01160.00560.29640.94480.00200.00520.26720.9410
β 21 0.16750.07900.64670.86060.10420.03570.52600.84990.07940.03040.47650.8885
σ 1 0.01030.07932.27901.00000.00340.06902.14091.00000.01040.10542.27881.0000
σ 2 0.05990.43956.86281.00000.09900.57983.12581.00000.20550.52804.78021.0000
λ 0.09090.13701.56470.91510.08930.55413.25420.92740.11991.35914.58930.9413
200 α 1 0.00550.06111.57851.00000.00730.05731.56891.00000.00120.06031.55901.0000
α 2 0.10220.03501.70301.00000.03160.03161.82311.00000.14460.02801.92421.0000
β 10 0.12990.01860.28720.87660.09050.01000.24900.87950.06560.00640.22200.9357
β 11 0.03600.00550.26930.95520.03140.00370.23230.94810.02600.00290.21030.9487
β 20 0.00140.00530.26910.94640.00510.00350.22800.94380.00110.00280.20410.9455
β 21 0.14070.07010.49150.91870.08590.03440.40860.90200.07240.02390.36280.8919
σ 1 0.00630.06051.61671.00000.00280.06451.65781.00000.00950.06731.69351.0000
σ 2 0.05270.25035.34911.00000.05080.19612.21761.00000.06720.17952.56551.0000
λ 0.08640.12101.05790.94970.08220.41122.15100.95800.07340.87403.26510.9598
Table 3. MLE (SE) for BVBeta, BVSJB, and BVSUSHN models.
Table 3. MLE (SE) for BVBeta, BVSJB, and BVSUSHN models.
ParametersBVSJBBVBetaBVSUSHN
α 1 0.1554 ( 0.0010 ) 57.0023 ( 2.2745 ) 1.5520 ( 0.1731 )
α 2 0.0624 ( 0.0018 ) 312.5664 ( 12.5410 ) 1.7849 ( 0.0035 )
β 1 5.0188 ( 0.9010 ) 8.2716 ( 0.3220 ) 0.2001 ( 0.0163 )
β 2 2.7763 ( 0.1020 ) 128.3064 ( 5.1770 ) 0.1729 ( 0.1590 )
σ 1 1.7204 ( 0.6972 ) 2.8302 ( 0.0097 )
σ 2 8.0839 ( 0.2439 ) 4.0036 ( 3.6538 )
λ 3.2156 ( 0.5746 ) 0.9995 ( 0.0001 ) 0.6239 ( 0.0547 )
KS test ( p -value) D 1 = 0.14 ( 6 × 10 10 ) D 1 = 0.12 ( 5 × 10 8 ) D 1 = 0.5725 ( 0.8987 )
D 2 = 0.07 ( 0.00234 ) D 2 = 0.08 ( 0.00234 ) D 2 = 0.5175 ( 0.9518 )
AIC 10 , 980.07 12 , 446.15 12 , 559.71
BIC 10 , 944.44 12 , 440.69 12 , 524.08
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Martínez-Flórez, G.; Vergara-Cardozo, S.; Tovar-Falón, R.; Rodriguez-Quevedo, L. The Multivariate Skewed Log-Birnbaum–Saunders Distribution and Its Associated Regression Model. Mathematics 2023, 11, 1095. https://doi.org/10.3390/math11051095

AMA Style

Martínez-Flórez G, Vergara-Cardozo S, Tovar-Falón R, Rodriguez-Quevedo L. The Multivariate Skewed Log-Birnbaum–Saunders Distribution and Its Associated Regression Model. Mathematics. 2023; 11(5):1095. https://doi.org/10.3390/math11051095

Chicago/Turabian Style

Martínez-Flórez, Guillermo, Sandra Vergara-Cardozo, Roger Tovar-Falón, and Luisa Rodriguez-Quevedo. 2023. "The Multivariate Skewed Log-Birnbaum–Saunders Distribution and Its Associated Regression Model" Mathematics 11, no. 5: 1095. https://doi.org/10.3390/math11051095

APA Style

Martínez-Flórez, G., Vergara-Cardozo, S., Tovar-Falón, R., & Rodriguez-Quevedo, L. (2023). The Multivariate Skewed Log-Birnbaum–Saunders Distribution and Its Associated Regression Model. Mathematics, 11(5), 1095. https://doi.org/10.3390/math11051095

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop