Flexible Log-Linear Birnbaum–Saunders Model

: Rieck and Nedelman (1991) introduced the sinh-normal distribution. This model was built as a transformation of a N(0,1) distribution. In this paper, a generalization based on a ﬂexible skew normal distribution is introduced. In this way, a more general model is obtained that can describe a range of asymmetric, unimodal and bimodal situations. The paper is divided into two parts. First, the properties of this new model, called ﬂexible sinh-normal distribution, are obtained. In the second part, the ﬂexible sinh-normal distribution is related to ﬂexible Birnbaum–Saunders, introduced by Martínez-Flórez et al. (2019), to propose a log-linear model for lifetime data. Applications to real datasets are included to illustrate our ﬁndings.

Rieck and Nedelman [1] developed the sinh-normal (SHN) distribution, which is given as the following transformation of a standard normal distribution γ > 0 is a shape parameter, ξ ∈ R is a location parameter and η > 0 is a scale parameter.

2.
The pdf of Y is unimodal for γ ≤ 2 and bimodal for γ > 2.

3.
The mean and variance of Y are E[Y] = ξ and Var[Y] = η 2 v(γ), where v(γ) is the variance when the scale parameter is equal to one (η = 1). There is no closed-expression for v(γ), but Rieck and Nedelman [1] provided asymptotic approximations for small and large values of γ.

4.
Let Y γ ∼ SHN(γ, ξ, η). Then, U γ = 2(Y γ −ξ) γη converges in law to a N(0, 1) distribution as The next lemma considers the particular case of a sinh-normal distribution in which the scale parameter is equal to 2, η = 2. This model is related to the Birnbaum-Saunders (BS) distribution. This result is the basis of the use of this distribution as a regression model, as shown in Section 5. Lemma 1.
Due to the result given in Lemma 1, the sinh-normal distribution is also named log-Birnbaum-Saunders distribution.
The following are particular cases of interest in the FSHN model: -If δ = 0, then the FSHN model reduces to the skew sinh-normal distribution studied by Leiva et al. [7]. -If λ = 0, then a symmetric model denoted by FSHN λ=0 (γ, ξ, η, δ) is obtained, which allows us to model symmetric bimodal data. -If δ = λ = 0, then the FSHN model reduces to the sinh-normal distribution introduced by Rieck and Nedelman [1].
Next, some results dealing with the inverse transformation to that introduced in (3) are given.

2.
If δ = 0, then we have the cdf of the skew-sinh-normal distribution.
The result in Corollary 1 for λ = 0 corresponds to the model studied by Olmos et al. [9], whereas the result for δ = 0 is a particular case of models studied by Vilca-Labra and Leiva-Sanchez [10].

Proof.
To obtain the maximums of f Y , we take the natural logarithm of f Y , except additive terms not depending of y. Thus, let us consider where s(y) = 2 γ sinh y−ξ η and c(y) = 2 γ cosh y−ξ η . Taking the first derivative of g(y) with respect to y, g (y) = 0, after straightforward calculations, (11) and (12) are obtained.
Then, the pdf of Y is symmetrical with respect to ξ and can be bimodal for δ < 0.
Proof. This result follows from the fact that, if y = ξ, then s(ξ) = 2 γ sinh(0) = 0 and the non-differentiability of the absolute value function at zero. Some illustrative plots for the pdfs in the FSHN model are given in Figures 1-3. In Figures 1 and 2, we mainly focus on the study of the shape parameters, δ and γ, so in every figure a positive or negative value of δ is fixed, and a range of possible values is considered for γ (less and greater than two), since our aim is to compare with the sinhnormal distribution, which is strongly unimodal for γ ≤ 2 and bimodal for γ > 2 [1]. In this way, the effect of δ and/or γ can be appreciated. Moreover, positive and negative values of λ are considered to assess the effect of the skewness parameter. In Figure 3, λ = 0 is fixed and therefore the symmetric submodel is plotted. Additional details are given next for every figure.  In Figure 1, without loss of generality for our purposes, the location and scale parameters are taken as ξ = 0 and η = 1. We intend to highlight the effect of δ on the bimodality of FSHN model. In Figure 1a, small values are considered for the skewness parameter, λ = −0.25 and 0.25, and a clearly negative value of δ = −2.75. All pdfs are bimodal. Note that the parameter γ ∈ {0.75, 1.5, 3.5, 5.5} and even for values of γ ≤ 2 the pdf is bimodal. These plots suggest that a negative value of δ induces bimodality and the effect of γ is attenuated. It is also possible to appreciate the effect of the sign of skewness parameter λ in Figure 1a; note that black and blue plots correspond to negative λ, whereas orange and red ones correspond to positive λ. On the other hand, Figure 1b, a positive value of δ = 1.5 is considered (λ = 1). We highlight that the distribution is unimodal even for values of γ > 2.
In Figure 2, a small value of δ is fixed (δ = 0.2 in Figure 2a and δ = −0.2 in Figure 2b); the plots suggest that, in this case, the bimodality depends on the value of γ. However, we must also be conscious that, in these plots λ = 0.4, for large values of λ unimodal distributions are nearly always obtained for any values of γ and δ. In Figure 3, the skewness parameter λ = 0 is fixed, therefore symmetric pdfs are always obtained. In Figure 3a, a negative value of δ = −0.75 is fixed and all plots are bimodal (even for values of the parameter γ < 2 in which the sinh-normal distribution is unimodal). In Figure 3b, δ = 0.75 is fixed and unimodal symmetric distributions are obtained.  Next, we show that the FSHN model is closed under location and scale changes.
and the result proposed is obtained.
Next, the limit behavior of FSHN model when γ approaches to 0 + is studied.
γη . By applying Proposition 6, Considering that lim Since φ(·) and Φ(·) are continuous differentiable functions, which is the pdf of a flexible skew-normal(δ, λ) distribution. That is, the result proposed.
In addition, note that this result is analogous to the property of the sinh-normal distribution given in Proposition 1.
Next, it is shown that the pth quantile of Y can be obtained from the pth quantile of Z ∼ FSN(δ, λ).
Then, the pth quantile of Y, y p , with 0 < p < 1 is given by where z p denotes the pth quantile of Z ∼ FSN(δ, λ).
Proof. It follows from (3) and the fact that the arcsinh function is monotonically increasing.
Submodel with λ = 0. For λ = 0, y p can be obtained from the quantiles of the N(0, 1) distribution.

Corollary 3.
If λ = 0, then As it can be seen in [1], recall that the log-Birnbaum-Saunders regression model is based on the relationship between the Birnbaum-Saunders and the sinh-normal distribution. In our proposal, a regression model is introduced next that is based on the relationship between the flexible Birnbaum-Saunders, studied in [2] and the FSHN distribution introduced in this paper. First, we recall the expression of the pdf for the flexible Birnbaum-Saunders.
Proof. This result is obtained by applying the transformation technique for random variables and taking into account (16).
In practice, it is usual to consider non-negative rv's following a FBS(α, β, δ, λ). Then, its neperian logarithm will follow a FSHN distribution, as established in next corollary.
If T depends on a known set of covariates then it is of interest to build a generalized linear model related to T [11]. This issue is studied in Section 5.

Moments
Let Y ∼ FSHN(γ, ξ, η, δ, λ). Recall that Thus, the moments of Y can be expressed in terms of the moments of arcsinh γ 2 Z with respect to Z ∼ FSN(δ, λ). To simplify the notation, let us introduce In particular, we have that

The Flexible Log-Birnbaum-Saunders Regression Model
In this section the flexible log-Birnbaum-Saunders regression model is introduced. Let T 1 , . . . , T n be independent positive continuous random variables such as Let us assume that the distribution of T i , proposed in (27), depends on a set of known explanatory variables, X i,1 , . . . , X i,p−1 , satisfying for i = 1, . . . , n 1.
The shape and skewness parameters in (27) Let us consider Y i = log(T i ). By applying Corollary 4, Y i ∼ FSHN(γ, x T i β, 2, δ, λ). Thus, Y i can be written as a linear model where the error term i ∼ FSHN(γ, 0, 2, δ, λ) (Proposition 6 is applied) and i s are independent. As for the expectation, variance of i and covariance of i and j (i = j), taking into account (19), we have for i = 1, . . . , n Taking where the initial vector of unknown parameters has been partitioned into β = (β 0 , β 1 T ) T with β 1 = (β 1 . . . , β p−1 ) T .
Particular cases of interest. 1.
If δ = λ = 0, then we have the log-Birnbaum-Saunders regression model proposed by Rieck and Nedelman [1], which is based on the sinh-normal distribution [1].
6. Applications 6.1. Illustration 1 Dataset 1 consists of n = 116 observations of daily ozone mean concentration in the atmosphere (in ppb= ppm × 1000) in New York taken from May to September 1973 [14]. The average concentration of pollutants in the air is of interest in epidemiological studies due to its serious adverse effects on the human health. It is usually assumed that these data are independent, and therefore they do not require a cyclic trend analysis, see for instance the works of Gokhale and Khare [15], Nadarajah [16] and Leiva et al. [7], where this dataset is also studied.
Descriptive summaries. The sample mean, variance, skewness and kurtosis coefficients are:ȳ = 42.1293, s 2 = 1088.2010, √ b 1 = 1.2098 and b 2 = 1.1122. From these summaries and Figure 4, we can conclude that this dataset is positively (or right) skewed. Thus, sinh-normal, skew sinh-normal, Birnbaum-Saunders, extended Birnbaum-Saunders [7] and FSHN distributions can be considered as models for this dataset. Their performance is compared throughout the Akaike information criterion, AIC = −2 l(·) + 2k and Bayesian information criterion, BIC = −2 l(·) + log(n)k, where l(·) denotes the log-likelihood function evaluated at the MLEs of parameters and k is the number of parameters in the model. From the results in Table 1, the FSHN distribution provides the best fit to these data since its AIC and BIC are the smallest [17,18]. On the other hand, since the sinh-normal, skew sinh-normal and FSHN are nested models, they can be compared by using the likelihood ratio test, [19].
First, let us consider the test H 0 : (δ, λ) = (0, 0) vs H 1 : (δ, λ) = (0, 0), which is equivalent to testing the sinh-normal (SHN) versus a FSHN distribution. The associated likelihood ratio statistic is ξ, η, δ, λ) , which is asymptotically distributed as a chi-square variable with two degrees of freedom (df). After substituting the estimated values of the parameters, we obtain −2 log(Λ 1 ) = 60.261, which is greater than the 5% chi-square critical value with 2 df, which is 5.99. Therefore, the FSHN model is preferred to the SHN model for this dataset. Second, we compare the FSHN model to the skew sinh-normal (skew-SHN) model. Let us now consider the following likelihood ratio statistic which is asymptotically distributed as a chi-square variable with d f = 1, χ 2 1 . After substituting the estimated values, we obtain −2 log(Λ 2 ) = 7.182, which is greater than the 5% chi-square critical value with 1 df, which is 3.84. Therefore, FSHN is preferred to skew-SHN model for this dataset. Figure 4 presents the histogram and the fitted pdfs. Note that FSHN and extended BS provide a good fit to this dataset.

Illustration 2
Dataset 2 consists of n = 500 observations of fetus weight (in grams) before birth obtained by ultrasound technique. Descriptive summaries are given in Table 2. The data are available at http://www.mat.uda.cl/hsalinas/cursos/2011/R/weight.rar (accessed on 17 May 2021). The variables b.weigt (fetal weight in grams) and a.weigt (birth weight in grams) can be found there. We study the b.weight variable. In Table 2, note that the mean and median are similar and the sample asymmetry coefficient √ b 1 is close to zero. Therefore, the distribution is fairly symmetrical. This fact can be seen in the histogram plotted in Figure 5, where we can also appreciate the bimodality of this dataset. Bimodality is due to the difference in weight, at the gestation stage, between males and females. Therefore, a symmetric bimodal model must be considered. In the literature, there are few models which exhibit both properties [20]. We propose the bimodal normal [21], the two-piece skew normal [22], the sinh-normal with γ > 2 and the submodel FSHN with λ = 0. In addition, we consider a mixture of two normal distributions whose pdf is with 0 < α < 1, µ i ∈ R and σ i > 0, i = 1, 2.  The estimated parameters in these models, along with their standard errors in parentheses, are given in Table 3. In this table, the following abbreviations are used: SHN (sinh-normal), TN (two-piece skew normal [22]), BN (bimodal normal [21]), MN (mixture of two normals given in (35)), FSHN λ=0 (FSHN with λ = 0). From AIC and BIC in Table 3, we can conclude that the mixture of two normals and the FSHN λ=0 model provide the best fit to this dataset. However, it is of interest to point out that the use of mixtures is a quite controversial issue in statistics, mainly due to non-identifiability problems (see, e.g., [23]).
It is obtained that −2 log(Λ 1 ) = 15.868, which is greater than the χ 2 1,0.05 = 3.84. Thus, H 0 is rejected, leading to the conclusion that the proposed FSHN λ=0 model fits better than the SHN considered in this case.
In Figure 5, the histogram of bebe weights is plotted along with the proposed bimodal models. There, the good fit provided by our proposal can be checked.

Illustration 3
Dataset 3 consists of 40 independent observations, which correspond to the failure time, T, for hardened steel specimens in a rolling contact fatigue test. The observations were taken at each of four values of contact stress, which is the covariate x in our proposal. This dataset can be seen in the work of Chan et al. [24]. Let Y i = log T i and consider the regression model We fit the log-Birnbaum-Saunders, the log-skew-Birnbaum-Saunders proposed by Lemonte et al. [13] and the log-flexible Birnbaum-Saunders introduced in this paper to this dataset. The estimated parameters along with the AIC and BIC are given in Table 4. According to AIC and BIC, log-flexible Birnbaum-Saunders provides a better fit than the other ones, so this model would be preferred.

Simulation
Next, a simulation study is presented to illustrate the performance of our results. Two features are studied: (1) the global performance of MLEs of parameters if the sample size increases; and (2) the effect of varying every shape parameter, γ, δ and λ, on the performance of the rest of estimators under consideration. Both issues are studied in the use of the flexible sinh-normal distribution as a regression model.
As sample sizes, we consider n ∈ {25, 50, 75, 100, 500}, m = 5000 simulations and values of the parameters (γ, β 0 , β 1 , δ, λ) = (1, 1.75, 1.25, −1.5, 2). As x, covariate values of a U(0, 1) distribution are considered. Values of the FSHN were generated from values of the FSN distribution, which were obtained by using the stochastic representation of FSN model given in [5]. The results are listed in Tables 5-7. As for measures of performance, the mean of estimates, their standard errors (sd) in parentheses, the relative bias and the square root of the mean squared error (MSE) are given. From the results in Tables 5-7, we can conclude that the MLEs are biased; in general, the standard error, relative bias and the squared root of MSE decrease if the sample size n increases. We highlight that, in this simulation, β 1 exhibits a good behavior, whereas the greater variability corresponds to λ.
Simulations were carried out by using the optim function of software R [25] by applying the Nelder-Mead method.

Simulation II
Next, a sensitivity study about the effect of varying the shape parameters (and the sample size) on the estimates of other parameters is presented.
Summaries varying the γ parameter are given in Table 8. Similarly, results varying the parameters δ and λ are given in Tables 9 and 10, respectively. For β 0 and β 1 , it can be seen that the relative biases are small. Moreover, we can appreciate that the relative bias, standard error and √ MSE decrease when the sample size increases, especially for β 1 . δ and λ exhibit greater variability, but their relative bias, sd and √ MSE also decrease when the sample size increases. As for γ, we can point out that this estimator is also well behaved when n increases. Its relative bias for small sample sizes, n = 25, is small and the √ MSE is a moderate value.
As final conclusion, we point out that these simulation studies cover a variety of situations as for the shapes of the FSHN distribution (unimodal and bimodal) and suggest that the estimators of the parameters are consistent when the sample size increases [26].

Discussion
The BS distribution is an asymmetric model used for survival time data and material lifetime subject to stress as it can be seen in [27]. Due to its practical and theoretical interest, a number of generalizations can be found in literature. We can cite the extensions provided in [28] to the family of elliptical distributions; in [10] based on the elliptical asymmetric distributions; and the extended Birnbaum-Saunders (EBS) distribution introduced in [7].
Other generalizations intend to solve specific deficiencies observed when this model is fitted to a dataset. In this sense, we can cite the epsilon-Birnbaum-Saunders model introduced in [8] to accommodate outliers; the extension based on the slash-elliptical family of distributions given in [3]; and the generalized modified slash Birnbaum-Saunders, which is based on the work in [29] and proposed in [30]. All these extensions are appropriate to fit data with greater or smaller asymmetry (or kurtosis) than that of the usual BS model, but they are not appropriate for fitting bimodal data. This issue is of great interest, as discussed by Olmos et al. [9], Bolfarine et al. [31], Martínez-Flórez et al. [2] and Elal-Olivero et al. [32]. The flexible BS distribution was proposed to model skewness and to fit data with and without bimodality. In this paper, this model is spread out to be used as a regression model in those situations in which regression models based on other generalizations of BS, such as those proposed in [7,13,[33][34][35][36][37], do not provide a satisfactory fit.