The Skewed-Elliptical Log-Linear Birnbaum–Saunders Alpha-Power Model

: In this paper, the skew-elliptical sinh-alpha-power distribution is developed as a natural follow-up to the skew-elliptical log-linear Birnbaum–Saunders alpha-power distribution, previously studied in the literature. Special cases include the ordinary log-linear Birnbaum–Saunders and skewed log-linear Birnbaum–Saunders distributions. As shown, it is able to surpass the ordinary sinh-normal models when ﬁtting data sets with high (above the expected with the sinh-normal) degrees of asymmetry. Maximum likelihood estimation is developed with the inverse of the observed information matrix used for standard error estimation. Large sample properties of the maximum likelihood estimators such as consistency and asymptotic normality are established. An application is reported for the data set previously analyzed in the literature, where performance of the new distribution is shown when compared with other proposed alternative models.


Introduction
When observed data does not follow a normal distribution, the use of the elliptical family of distributions is an important alternative. In this family, all distributions are of the symmetric type, including distributions with higher (or lower) kurtosis than the one for the normal distribution.
A random variable (RV) X is distributed according to the elliptical distribution with location parameter ξ ∈ R and scale parameter η > 0. We use the notation X ∼ EC(ξ, η; g).
Special cases of the distribution X ∼ EC(0, 1; g), are given by the Pearson type VII, Kotz, Student t, Cauchy, and the normal distributions, among others. Although the class of elliptic distributions is a good alternative for situations of departure from normality, it is not appropriate when observations follow an asymmetric distribution. These circumstances prompted the search for new distributions better suited to fit data with high asymmetry and kurtosis. It is then necessary to extend the family of elliptical distributions to a more general family of distributions typically called the skew-elliptical distribution.
The probability density function (pdf) for the standard scalar skew-elliptical model is given by f Y (y; λ) = 2 f (y)F(λy); y, λ ∈ R (1) where f (x) = c η g x−ξ η 2 , F is the cumulative distribution function (cdf) corresponding to f , and λ is an asymmetry parameter. This model is denoted by Y ∼ SE(0, 1; g, λ). The cdf of this model is given by Hence, for λ = 0, the standard elliptical model is obtained. A particular case of model (1) is the skew-normal distribution (see Azzalini [9]) with f = φ(·) and F = Φ(·) (where φ(·) is the pdf of standard normal distribution and Φ(·) is the respective cdf) with pdf given by, and cdf: where T(., .) is the well known Owen's function [10]. Another type of asymmetric distribution studied by Durrans [11], is the fractional order statistics model, with pdf where H is an absolutely continuous distribution function with pdf h and α > 0 is a parameter that controls asymmetry and kurtosis of the distribution. This model is denoted by AP(α). It has been used to fit data with high amounts of asymmetry and kurtosis. Moreover, Birnbaum and Saunders [12] introduced a new probability distribution for modeling the lifetime of certain structures under dynamic load called "fatigue life distributions", which is also asymmetric but only fits positive data. Its pdf is given by where γ > 0 is the shape parameter and τ > 0 is a scale parameter and the distribution median. We use the notation T ∼ BS(γ, τ). It is well known that the BS distribution becomes asymmetric as γ increases and symmetric around τ as γ gets close to zero. Another property includes kT ∼ BS(γ, kτ), for k > 0, and T −1 ∼ BS(γ, τ −1 ), as can be seen in Birnbaum and Saunders [12,13].
Díaz-García et al. [14] presented a generalization of this distribution for the case of symmetric-elliptic distributions and an extension to the case of asymmetric-elliptic distributions was studied by Vilca-Labra et al. [15], which is called doubly generalized Birnbaum-Saunders distribution.
Vilca-Labra and Leiva-Sanchez [15] investigated asymmetric extensions for the Birnbaum-Saunders model, with the normal distribution replaced by the elliptical asymmetric class of distributions [16]. Furthermore, Castillo et al. [17] considered the asymmetric epsilon-Birnbaum-Saunders model, Gómez et al. [18] considered an extension based on the slash-elliptical family of distributions, Martínez-Flórez et al. [19] studied the extension of this model to the case of exponentiated or alpha-power distributions, and Moreno-Arenas et al. [20] studied the extension to the case of hazard proportional family.
Rieck and Nedelman [21] developed the sinh-normal distribution, which is obtained as a transformation of the standard normal distribution after considering the RV Y = arcsinh(γZ/2)η + ξ, where Z ∼ N(0, 1), γ > 0 is a shape parameter, ξ ∈ R is a location parameter and η > 0 is a scale parameter. We use the notation Y ∼ SHN(γ, ξ, η). It is well known that if Y ∼ SHN(γ, ξ, 2) then the RV T = exp(Y) ∼ BS(γ, exp(ξ)). For this reason, the sinh-normal distribution is typically called the log-Birnbaum-Saunders distribution. A log-linear extension allowing incorporating covariates was also considered in Rieck and Nedelman [21] and can be written as where Y i is the log-survival time corresponding to the i-th experimental unit, β = (β 1 , β 2 , . . . , β p ) T an unknown parameter vector, x i = (x i1 , x i2 , . . . , x ip ) T a p-dimensional vector with the (known) values of the explanatory variables and i , are independent and identically distributed RV with i ∼ SHN(γ, 0, 2), i = 1, 2, . . . , n. Diagnostic analysis technics for the SHN model were developed by Galea et al. [22]. More recently, Barros et al. [23] extended this model for the case of error distributions with heavier tails emphasizing the use of the Student t distribution. They also conducted estimation and diagnostic studies for the model entertained.
Extensions for the sinh-normal model using an asymmetric setup were studied in Leiva et al. [24], where a skew-sinh-normal model is developed and used in a study of the air pollution in the city of Santiago de Chile. Some other asymmetric extensions of the sinh-normal models are reported in Lemonte [25], Santana et al. [26], and Martínez-Flórez et al. [19], who studied the extension of this model to the case of alpha-power distributions.
We consider in this paper an extension of the linear model (7) by replacing the sinhnormal distribution of the errors by the skewed-elliptical sinh alpha-power distribution.
The introduction of two extra shape parameters makes the error distribution more flexible, allowing it to incorporate additional kurtosis and asymmetry. Inference is conducted using the maximum likelihood approach while the observed information matrix is derived for computing standard errors. The new model is fitted to a real data set, revealing that it has great potential in applied scenarios.
The paper is organized as follows. In Section 2, the skewed-elliptical sinh-alphapower distribution is defined and some of its main properties are studied. In Section 3, the new model is defined is defined, maximum likelihood estimation is discussed, and the score function and observed information matrix for the model parameters are studied. In Section 4, we carry out a simulation study to assess the performance of the proposed estimators. In Section 5, we apply the proposed model to analyze two data sets. Concluding remarks are given in Section 6.

Skew-Elliptical Alpha-Power Model
We now extend the model AP(α). This new model has pdf given by where H is the distribution function corresponding to h. We call it the skew-elliptical alpha-power distribution and we will denote it by Z ∼ APSE(0, 1; g, λ, α). Moments of the RV Z have no closed form, but under a variable change, the r-th moment of the RV Z can be written as where H −1 is the inverse of the function H. If h is of the form (3); then, we have the skewnormal alpha-power model of parameter λ and α. This model is denoted PESN(λ, α). Special cases of model PESN occur with α = 1 so that the skew-normal model φ SN (x) follows. On the other hand, with λ = 0, the model with pdf ϕ Φ (x), that is, Durrans's generalized normal model, follows. The ordinary standard normal model is also a special case that follows by taking α = 1 and λ = 0, that is, ϕ SE (x; 0, 1) = φ(x). Notice from Figure 1a,b below that α and λ affect both distribution asymmetry and kurtosis; hence, the model proposed seems more flexible than the models by Azzalini [9] and Durrans [11].

Skew-Elliptical Sinh-Alpha-Power Model
In this subsection, we present the pdf for the skew-elliptical sinh-alpha-power distribution. Expressions are derived for the model mean and variance and consequences of changing parameter values in the density are reported.
The skew-elliptic sinh-alpha-power distribution is defined as a transformation of the skew-elliptical alpha-power distribution, that is, where γ, α ∈ R + and λ ∈ R are shape parameters, ξ ∈ R is a location parameter and η > 0 is a scale parameter. The pdf of Y is given by which we denote by Y ∼ SAPSE(γ, ξ, η; g, λ, α). When h = φ SN and H = Φ SN , the skew-sinh-normal power distribution follows, with pdf given by which we denote by Y ∼ SPESN(γ, ξ, η, λ, α). Particular case: • sinh-normal model, when λ = 0 and α = 1, • skew sinh-normal distribution, when α = 1 • sinh-normal power-normal distribution, when λ = 0 The cdf of Y is given by with Φ SN (·) as above. Accordingly to (12), the inversion method can be used for generating from an RV with distribution SPESN(γ, ξ, η, λ, α). That is, if U ∼ U(0, 1), the uniform distribution, then the RV Y = ξ + η arcsinh γ 2 Φ −1 SN (U 1/α ) is distributed according to the SPESN distribution with parameters γ, ξ, η, λ and α where Φ −1 SN is the inverse of the skew-normal distribution. Figure 2 reveals the behavior of the pdf in (11) for ξ = 0 and η = 2 and different values for γ, λ, and α. Let Y ∼ SPESN(γ, ξ, η, λ, α). Then, the mode (or modes) of the RV Y is the solution of the nonlinear equation Figure 2, the RV Y is asymmetric around location parameter ξ and has a single mode. Figure 3 shows the bimodal case. Moreover, expectation and variance of Y are given by where c j and V are functions of γ, λ, and α, given by and V(γ, λ, α) is the variance of the RV arcsinh(γZ/2). There are no closed form expressions for c j and V.
The Y moment-generating function can be written in the form where E(T r ) can be obtain using the expression given in Theorem (2.6) de Martínez-Flórez et al. (2020) [27]. Let Then, we can write as γ → 0. Therefore, U converges in distribution to a standard PESN distribution. The survivor function, cumulative risk function, inverted risk function, and risk function for model SPESN are given, respectively, by: where r SESN (t) and R SESN (t) are the indices for the risk and inverted risk for the skewed sinh-normal model, that is, the inverse risk rate is proportional to the risk rate for the skewed sinh-normal distribution. Therefore, the intervals where R(t) is decreasing or increasing are the same intervals where R SSN (t) is decreasing or increasing.   The density generator of the generalized Student t, type I logistic, type II logistic, and power exponential are, respectively, given by As a natural step, we now introduce the skew-alpha-power BS model from the skewsinh-power-normal model. This result is stated in the following theorem. ξ, 2, λ, α), then RV T = exp(Y) follows the skewed Birnbaum-Saunders alpha-power distribution, with shape parameters γ, λ, α and scale parameter exp(ξ). We denote ∼ PESNBS(γ, exp(ξ), λ, α).
As in the classical BS model and the asymmetric BS model of Vilca-Labra et al. [15], this new model is an alternative for modeling the lifetime of certain structures under dynamic load on joint observations with asymmetry and/or kurtosis outside the range permitted by these distributions. Properties of the PESNBS model can be obtained from the properties derived for the BS and skew-BS models. τ, λ, α). Then, 1.

Observed Information Matrix
The observed information matrix follows from minus the second derivatives of the log-likelihood function, J(θ) = − θθ , where θθ is the Hessian matrix with parameters replaced by their MLEs. The elements of the observed information matrix are denoted by j β j β k , j γβ j , . . . , j αα , and are given by The expected (or Fisher) information matrix is n −1 times the expected value of the elements of the observed information matrix.
For λ = 0 and α = 1 and using the approximation in Chaibub Neto and Branco [28], we can write the expected Fisher information matrix as x exp(−t 2 )dt being the error function (see Prudnikov et al. [29]) and It can be verified that the columns (lines) of the matrix I(θ) are linearly independent; hence, it is invertible. Therefore, for large samples, the MLEθ of θ is asymptotically normal, that is,θ A → N p+4 (θ, I F (θ) −1 ), resulting that the asymptotic variance of the MLEθ is the inverse of I F (θ), which we denote by Σ θ = I F (θ) −1 , which is estimated replacing parameters by their MLEs.

Simulation Study
In this section, a simulation study was performed to evaluate the performance of the ML estimators for the model parameter.
To check the good behavior of ML estimators in the SPESN model, a simulation study is carried out. In this study, 1000 samples of sample size n = 100, 200, and 300 have been generated from the SPESN model. The aim of this simulation is to study the performance of the ML estimators for the parameters using the proposed procedure.
For each generated sample, the ML estimates for each parameter were computed numerically. The bias, standard error (se) and mean squared error (MSE) of the estimates for the parameters (β 0 , β 1 , γ, λ, α) are given as summaries in Table 1. It can be seen there that the ML estimates are quite stable and close to the real values for the sample sizes under consideration. As expected, the bias and standard deviation of our estimates decrease as the sample size n increases.  The data sets are of times to failure (T) in rolling contact fatigue of ten hardened steel specimens tested at each of four values of four contact stress points (x). The data were obtained using a four-ball rolling contact test rig at the Princeton Laboratories of Mobil Research and Development Co. The data set was given initially in McCool [30] and reported in Chan et al. [31]. These data sets were also analyzed in Lemonte [25], who, like Chan et al. [31], considered the regression model For the regression model described above, we propose using the sinh-normal (SHN) distribution, skewed sinh-normal (SESN), and skewed sinh power-normal (SPESN) distributions for i , i = 1, 2, . . . , n. To compare the fitness of these models, we use the Akaike [32] criterion, namely where k is the number of parameters in the model. According to this criterion, the model that best fits the data is the one with the lowest AIC. The SESN model can be obtained as a special case of the SPESN model when α = 1. The estimates for the parameters of these models are presented in Table 2. According to the AIC criterion, the SPESN model presents the best fit to the data set. Making the correction in the intercept of the SPESN model, we find thatβ * 0 = 0.289, so we estimate the vectorθ * = (0.289, −13.602) . We can conclude that the regression model with SPESN error distribution provides a better fit than the regression model with SESN error distribution.
Note that the skewed log-linear Birnbaum-Saunders alpha-power model provides better fit to the data set the corresponding ones for the remaining models.
To confirm the good fit of the distributions used for the error term, we plotted the transformed standardized residual scale Z i = (2/γ) sinh(Y i − x β)/2 for the distribution of the estimated errors. Under this scale, the distribution of Z is normal for the SHN model, while for the SESN model, Z is the ordinary SN distribution and for the SPESN model, Z is PESN distribution. Figure 5 depicts the distribution for the scaled residuals Z for the set of models with the corresponding theoretical distributions.  Figure 6 shows the qq-plots for the scaled residuals Z for the set of models SSN and SESN with the indicated theoretical distributions. One can see the good fit of the model with errors SPESN. Thus, this model is presented as a viable alternative to study data when the distribution of the response variable is asymmetric.

Application 2: Patients with Advanced Lung Cancer Status
The second application reported is related to the censored data set presented in Kalbfleisch and Prentice [33]. The uncensored data set under study corresponds to 128 male patients with advanced lung cancer status. The study aims at explaining the survival time T (in days) using a regression model. The explanatory variables that we consider were the following: a randomized measure of patients condition (Karnofsky). The data set was also analyzed in Martínez-Flores et al. [19]. Consider the regression model where (x 1 ) is 10-30 fully randomized, 40-60 partially, and 70-90 patient's ability for selfcare and x 4 histologic tumor squamous type: dichotomized as 1 for "yes" and 0 for "no" The estimates for the parameters of these models are presented in Table 3. According to the AIC criterion, the SPESN model presents the best fit to the data set. We can conclude that the regression model with SPESN error distribution provides a better fit than the regression model with SESN error distribution.
We consider now the problem of testing the null hypothesis of no difference between the SPESN model and the traditional SHN model, using the likelihood ratio statistics (models are nested) Λ 1 = L SHN (θ) . After numerical evaluations, we obtain −2 log(Λ 1 ) = 169.422, which is greater than the 5% critical value 5.991. Therefore, the null hypothesis is rejected and we conclude that the SPESN model (which involves two extra parameters, making it more flexible in terms of asymmetry and kurtosis) fits the data better (in fact, much better) than the SHN model. The SPESN model is also compared to the SESN model by testing the hypothesis, using the likelihood ratio statistics Λ 2 = SESN (θ) After numerical evaluations, we obtain −2 log(Λ 2 ) = 39.495, which is greater than critical value χ 2 1,95% = 3.841, so that H 02 is rejected at the 5% critical value. Note that skewed log-linear Birnbaum-Saunders alpha-power model provides better fit to the data set the corresponding ones for the remaining models.
To confirm the good fit of the distributions used for the error term, we plotted the transformed standardized residual scale Z i = (2/γ) sinh(Y i − x β)/2 for the distribution of the estimated errors. Under this scale, the distribution of Z is normal for the SHN model, while for the SESN model, Z is the ordinary SN distribution and for the SPESN model, Z is PESN distribution. Figure 7 depicts the distribution for the scaled residuals Z for the set of models with the corresponding theoretical distributions.  Figure 8 shows the qq-plots for the scaled residuals Z for the set of models SSN and SESN with the indicated theoretical distributions. One can see the good fit of the model with errors SPESN. Thus, this model is presented as a viable alternative to study data when the distribution of the response variable is asymmetric.

Conclusions
This paper introduced an extension of the log-linear Birnbaum-Saunders model by replacing the sin-normal distribution for the error term by the skew-sinh-power-normal distribution. Some features of the new model are: • The introduction of two extra shape parameters makes the error distribution more flexible, allowing it to incorporate additional kurtosis and asymmetry. • Maximum likelihood properties of large samples such as consistency and asymptotic normality were established. • A simulation study was performed to evaluate of the ML estimations. As expected, the bias and standard deviation of our estimation decrease as the sample size n increases. • In the applications, AIC criteria statistics are used. These criteria indicate that the model that best fit the data is SPESN model. • A Bayesian approach will be worked on in a future work.