Log-epsilon-skew normal: A generalization of the log-normal distribution

1. Introduction

The epsilon-skew normal (ESN) distribution was introduced as a skew extension of the normal distribution in same spirit as Azzalini’s (1985) skew-normal (SN) distribution. In this paper we introduce and investigate properties of a family of non-negative random variables (r.v.’s) called the log-epsilon-skew normal (LESN) r.v.’s. A r.v. is said to have the LESN distribution if its logarithm has the ESN distribution. A key advantage of the LESN family is that it includes as a subfamily the log-normal distributions, which are well known to be in use in a wide spectrum of disciplines. Examples include astrophysics, Gandhi (2009), environmental sciences, Benning and Barnesa (2009), computer science, Doerr et al. (2013), economics, Cobb et al. (2013), biomedical, Feng et al. (2013), and radiology, Neti and Howell (2008). Eckhard et al. (2001) compared the use of the log-normal distribution across several different science disciplines.

The log-normal distribution has a long and rich history. Galton (1879) suggested the use of the log-normal distribution to analyze data for which the geometric mean is better than the arthimatic mean for estimating central tendency. McAlister (1879) derived the log-normal distribution at Galton’s suggestion. Finney (1941) examined the moments, moment estimation, and efficiency of the estimation. Heyde (1963), showed the log-normal distribution is not uniquely determined by its moments. It is known that the moment generating function for the log-normal distribution does not exist, see Romano and Siegel (1986). However, it is well-known that the moments of the log-normal distribution can be derived from the moment generating function of the normal distribution. For a review of the history of log-normal distribution and its applications as a generative law see Mitzenmacher (2004).

In terms of background, a LN(θ; σ) r.v. T is defined by its probability distribution function (p.d.f.).

where Φ(·) denotes the standard normal cumulative distribution function (c.d.f). Its c.d.f. is given by

In survival analysis, the log-normal distribution is used in accelerated failure time (AFT) models, see Collett (2003). The form of the log-normal hazard function is given by

The log-normal hazard function is unimodal starting at zero and increasing to a maximum, then decreases toward zero as t → ∞. Sweet (1990) gives mathematical properties of the hazard rate for the log-normal distribution.

One competing approach towards generalizing the log-normal distribution would be to utilize the skew normal distribution of Azzalini (1985, 1986) more recently extend by Domma et al. (2015). This approach has been utilized in practice by Azzalini and Kotz (2002) for modeling income data and Chai and Bailey (2008) for the purpose of modeling continuous data with a discrete component at zero. However, the mathematical properties have not been investigated.

This paper is organized as follows. In Section 2, we introduce the LESN distribution and discuss its properties. In Section 3, we examine the maximum likelihood estimates of the LESN parameters. In Section 4 we summarize the results of a simulation study to examine the behavior of the maximum likelihood estimates with and without the presence of censoring. We apply these methods to two real-life data sets in Section 5. In Section 6, we finish with some conclusions.

2. The log-epsilon-skew-normal distribution

For the further development in this section the ESN model is used as a starting point. Details of the epsilon-skew-normal (ESN) distribution, developed by Mudholkar and Hutson (2000), are given in the Appendix. Toward this end let T be a random variable such that log (T)~ESN(θ, σ, ϵ): It follows that T~LESN(θ, σ, ϵ). The pdf, cdf, and the quantile function of the LESN distribution are given as follows:

and,

respectively, where −1<ϵ<1 and Q₀(u) is given by Equation (A.3) in the Appendix. As can be seen at ϵ = 0, Equations (2.1) and (2.2) give Equations (1.1) and (1.2). Figure 1 shows some typical LESN probability density functions with θ = 0 and σ = 1. A horizontal line is drawn at the splice point x = e^θ= 1: We see that the LESN pdf is unimodal, with the mode before the splice point t = 1.

3. Maximum likelihood estimation

Let Y₁; Y₂, …, Y_n denote i.i.d. Y_i~LESN(θ, σ, ϵ): The loglikelihood for the i-th observation is given by

Denote the loglikelihood for the i-th observation for the ESN(θ, σ, ϵ) as l_ESN (θ, σ², ϵ|x_i), then l_ESN (θ, σ², ϵ|y_i) = log y_{i +} l_ESN (θ, σ², ϵ|log y_i). We denote the maximum likelihood estimates of (θ, σ², ϵ) as $\left(\widehat{\theta },{\widehat{\sigma }}^2,\widehat{ϵ}\right)$.

Let y₍₁₎ ≤ y₍₂₎ ≤ … ≤ y_(n) denote the order statistic of a sample from the LESN(θ, σ, ϵ) population. Also denote y₀ = 0 and y_n+1 = ∞. We apply some of the results from Mudholkar and Hutson (2000) to the LESN(θ, σ, ϵ) distribution.

Lemma 3.1. There exists and integer k = k(y₍₁₎, …, y_(n)), 0 ≤ k ≤ n, such that the loglikelihood l(θ, σ², ϵ) can be expressed as

Proof. The existence of k is obvious. The form of the loglikelihood follows from the p.d.f. of LESN(θ, σ, ϵ) given by Equation (2.1). Note that the case k = 0 corresponds to the case of the p.d.f. when ϵ = −1 and the case k = n corresponds to the case of the p.d.f. when ϵ = 1. In other words for k = 0 and k = n the loglikelihood comes directly from the half-normal distributions.

The estimation of k, i.e. identification of the interval containing the maximum likelihood estimate of the mode, is helped by the following two lemmas.

Lemma 3.2. If k = 0 or k = n the maximum likelihood estimate $\left(\widehat{\theta },{\widehat{\sigma }}^2,\widehat{ϵ}\right)$ is given by

where $s_0^2={\sum }_{i=2}^n{\left(\log y_{\left(i\right)}-\log y_{\left(1\right)}\right)}^2/\left(4n\right)$ and $s_n^2={\sum }_{i=2}^{n-1}{\left(\log y_{\left(i\right)}-\log y_{\left(n\right)}\right)}^2/\left(4n\right)$.

Proof. If k = 0, then ϵ = −1; and in the loglikelihood, x_(i) ≥ θ, i = 1, 2, …, n. A maximization of the loglikelihood then yields $\widehat{\theta }=x_{\left(1\right)}$, the sample minimum, and ${\widehat{\sigma }}_0^2=s_0^2$.

The proof of Lemma 3.2 at k = n follows similarly.

Lemma 3.3. If 1 ≤ k<n then the local minima of

j = 1, 2, …, n−1, if any, determines the local maxima of the loglikelihood.

Proof. If 1 ≤ k<n then, for a fixed σ²; the loglikelihood l(θ, σ², ϵ) at Equation (3.1) is maximized by examing the local minima of

where −1<ϵ<1, and j = 1; 2, … n−1. Now g_j(y, θ; ϵ) may be minimized sequentially, first with respect to ϵ for a fixed θ (and consequently for fixed k), and then with respect to θ. Thus, equating the derivative of g_j(y, θ, ϵ) with respect to ϵ to zero we have:

Eliminating ϵ between Equations (3.5) and (3.6), and simplifying, we get the minimum of g_j(y, θ, ϵ) with respect to ϵ as

Thus the problem is reduced to minimizing g_j(y, θ), or equivalently h_j(y, θ)= (4g_j(y, θ))^1/3, with respect to θ subject to log (y_(j)) ≤ θ log (y_(j+1)). Clearly, h_j(y, θ) is a differentiable function of θ with derivative

It follows that the solutions of $h_j^{\text{'}}\left(\mathbf{y},\theta \right)=0$, j = 1,2, …, n−1, if any, correspond to the local stationary points of the loglikelihood.

It is observed that in general the function g_j(y, θ) is monotone in most intervals. In other words, most of the half-open intervals [log(y_(j)), log (y_(j+1))), do not contain the local minima of g_j(y, θ). The problem of identifying the intervals which have the local minima can be efficiently handled by graphing the continuously differentiable (except possibly at θ = log (y₍₁₎) and log (y_(n))) function

where I_j(θ) is the indicator function for the half-open interval [log(y_(j)), log (y_(j+1))), and g_j(y, θ) is given by Equation (3.7). However, the isolation of the global maximum of the likelihood from the set of the local minima identified as in Lemma 3.3 requires the values of the likelihood at these points. Towards this end we have the following:

Lemma 3.4. Let θ₀ denote a solution of

j = 1, 2, …, n−1. Then the corresponding values ${\sigma }^2={\sigma }_0^2$ and ϵ = ϵ₀ which maximize the loglikelihood are given by

Proof. Follows directly from Equations (3.1) and (3.6).

The value of the integer j, i.e. the interval containing the maximum likelihood estimator of the mode, is determined by the value of the loglikelihood over the set of all $\left({\theta }_0,{\sigma }_0^2,{ϵ}_0\right)$ given by Lemmas lem:lem2 and lem:lem3.

Lemma 3.5. The global maximum of the loglikelihood is obtained by examining the following local maxima for the set k = 0, k = n, and the set $\left({\theta }_0,{\sigma }_0^2,{ϵ}_0\right)$ given by Lemmas 3.2 and 3.4:

where $s_0^2$ and $s_1^2$ are given by Lemma 3.2 and ${\sigma }_0^2$ is given by Lemma 3.4.

For ϵ = −1 θ becomes a threshold parameter with log (y₍₁₎) as a consistent and reasonable estimator. The scale parameter may then be estimated by substituting log (y₍₁₎) for θ and maximizing the likelihood. The asymptotic theory of estimation in the presence of a threshold parameter, e.g. as discussed by R.L. Smith (1985), is then operative. The analogous observation applies when ϵ = 1.

Theorem 3.6. As n → ∞, the maximum likelihood estimator $\sqrt{n}\left({\widehat{\theta }}_{ml}-\theta ,{\widehat{\sigma }}_{ml}2-{\sigma }^2,{\widehat{ϵ}}_{ml}-ϵ\right)$ is a centered trivariate normal distribution with variance-covariance matrix

Proof. The regularity conditions for the validity of the large sample theory for maximum likelihood estimation (e.g. Rao 1973, Section 5g) are obviously satisfied. Also, ∑_ml is the inverse of the expectation of the total sample information matrix given by the density in Equation (2.1) and may be computed readily. The proof is completed by taking expectations using the score equations. Theorem 3.7. If ϵ = ϵ₀ then as n → ∞, $\sqrt{n}\left({\widehat{\theta }}_{ml}\left({ϵ}_0\right)-\theta ,{\widehat{\sigma }}_{ml}2\left({ϵ}_0\right)-{\sigma }^2\right)$ is a centered bivariate normal distribution with variance-covariance matrix

In the context of LESN modeling the parameter ϵ may be interpreted as a measure of distance of the distribution from log-normality with respect to LESN alternatives. Hence, the test

may be used as a diagnostic tool with respect to the appropriateness of testing the appropriateness of an underlying log-normal model versus a broad class of LESN alternatives with the large sample p-value given via the results of Theorem 3.7. A parametric bootstrap approach given by Ristic et al. (2018) may be used to test the validity of the model given real datasets.

4. Monte Carlo experiments

In this section we outline two Monte Carlo experiments performed in order to study the behavior of the maximum likelihood based estimates of the LESN parameters and summarize its results. The first study examines the uncensored case and the second considers the randomly censored data. Starting values were given by considering the maximum likelihood estimates for θ and σ given by a log-normal distribution and setting ϵ = 0 with parameter estimation carried out using classic Newton-Raphson methodologies.

In the uncensored case, we let T~LESN(0, 1, ϵ) with sample sizes n = 25, 50, 100, 200, and 500. We utilized 10000 simulations under each scenario. We let ϵ range from 0 to 0.75 by 0.25, which includes the special case of ϵ = 0 corresponding to the log-normal distribution.

Tables 1 and 4 summarize the means ± standard deviations of the MLE’s for each of the parameter configurations estimated on the basis of 10000 replications in the simulation study. As expected, the finite bias approaches zero and the variance of the parameter estimates decrease as the sample size increases. Hence, via simulation we have that the maximum likelihood estimates of the LESN parameters are well-behaved, thus validating the results at Equation (3.14).

For the censored case, we needed to simulate T_i = min(X_i, C_i); where X_i denotes the failure time and C_i denotes the right censoring time for subject i = 1, 2, …, n. We assume that the censoring times, C_i, are identically independently exponential random variables with mean 1/λ. That is C_i~Expo(λ). We fix μ = 0 and σ = 1. Let ϵ = −0:75 to 0.75 by 0.25.

Values of λ were chosen to obtain samples with expected censoring proportions of 10%, 20%, 30%, 40%, and 50%. Finding the values of λ analytically would require solving

for λ, where c is the expected censoring proportion. This problem is equivalent to choosing a value of λ so that E(δ_i)= 1−c; where δ_i = 1 if X_i ≤ C_i and 0 otherwise. This was done via simulations.

Tables 1–4 summarize the maximum likelihood parameter estimation of the LESN model parameters θ, σ, and ϵ. We see there is finite bias in the parameter estimates. As expected, the bias approaches zero and the variance of the parameter estimates decrease as the sample size increases. This bias does increase as the censoring proportions increase as would be expected.

5. Two examples

For the purpose of illustration we first look at the Mayo Clinic primary biliary cirrhosis data that can be downloaded from http://www.umass.edu/statdata/statdata/data/, Fleming and Harrington (1991). In particular, we examine 418 serum bilirubin (mg/dl) measurements. Some basic statistics are reported in Table 5. In Table 6 maximum likelihood estimators are presented for log-normal and LESN. It is interesting to note the p-value for testing the hypothesis Equation (3.16) is p < 0.0001. This suggests that the log-normal model may not be appropriate. Starting values for both examples were given by considering the maximum likelihood estimates for θ and σ given by a log-normal distribution and setting ϵ = 0 with parameter estimation carried out using classic Newton-Raphson methodologies. The p-value calculations follow directly from Theorem 3.7.

Figure 3a depicts histogram of the serum bilirubin measurements and model fitting for log-normal and LESN models. We see that the LESN model provides a better fit. This is supported by Figure 3b which depicts a histogram of the natural logarithm of the serum bilirubin measurements and model fitting for normal and ESN models. Even after the log transformation, the data are clearly right skewed.

For our second example, we look at data that was part of a large clinical trial carried out by the Radiation Therapy Oncology Group in the United States, Kalbfleisch and Prentice (1980). The data may also be downloaded from http://www.umass.edu/statdata/statdata/data/. In particular, we examine the survival time in days from date of diagnosis in 195 patients with squamous carcinoma in the oropharynx. There were 53 patients (27%) whose survival times were right censored. In Table 7 maximum likelihood estimators are presented for log-normal and LESN. It is interesting to note that the p-value for testing Equation (3.16) is p = 0.0174. As in the previous example, this suggests that the log-normal model may not be appropriate. Figure 4 shows the estimated hazard function for log-normal and LESN models.

6. Conclusions and discussion

In this paper we have proposed a generalization of the log-normal distribution called the log-epsilon-skew-normal distribution, obtained some of its basic properties, and developed its adaptation for the analysis of the randomly censored survival data including the maximum likelihood estimation, including the associated information matrices, of its parameters in uncensored and random censored cases. The MLE’s in the two cases are also studied empirically for the moderate sample-size situation and the application are illustrated by two examples. It is noted that the entire development is based on the ESN distribution. We can adapt methodology developed by Hutson (2004) to extend the use of LESN in log-transformed response variables in common regression problems. It is to be noted that if the sd parameter in the ESN distribution is replaced by an independently distributed $\sqrt{{\chi }_v^2/v}$ variable or, alternatively, treating student’s t-distribution in place of the N(μ, σ) distribution as the normal distribution in the ESN definition, we get an epsilon skew-t distribution (Gómez et al. 2007). This model, which includes ES-Cauchy as a particular case, is more flexible in terms of both skewness and kurtosis. This model together its log-modification is currently under investigation. Similar developments involving the ES-logistic and LES-Logistic which may be more appte for the quantal response analysis are also being pursued.