Skewness of Maximum Likelihood Estimators in the Weibull Censored Data

: In this paper, we obtain a matrix formula of order n − 1/2 , where n is the sample size, for the skewness coefﬁcient of the distribution of the maximum likelihood estimators in the Weibull censored data. The present result is a nice approach to verify if the assumption of the normality of the regression parameter distribution is satisﬁed. Also, the expression derived is simple, as one only has to deﬁne a few matrices. We conduct an extensive Monte Carlo study to illustrate the behavior of the skewness coefﬁcient and we apply it in two real datasets.


Introduction
In its first appearance, the Weibull distribution [1] claimed its wide applicability. Survival analysis, reliability engineering, and extreme value theory are some of its applicability. To amplify the relevance of the Weibull, a regression structure is added to one of the parameters, i.e., the behavior of the distribution may be explained from covariates (explanatory variables) and unknown parameters to be estimated from observable data.
In statistical inference, it is often desirable to test if there are regression parameters statistically significant and the Wald test is commonly performed. Under standard regularity conditions, the null distribution of the Wald statistic is asymptotically chi-squared, a consequence of the maximum likelihood estimators (MLE) distribution. Therefore, the Wald test must be avoided if the sample size is not large enough, because the distribution of the MLE will be poorly approximated by the normal distribution.
Preventing the complexity of the statistical tests, the skewness coefficient (say γ) of the distribution of the MLE is an easy way to verify if the approximation to normality is adequate. A value of γ far from zero indicates a departure from the normal distribution. Pearson's standardized third cumulant defined by γ = κ 3 /κ 3/2 2 , where κ r is the rth cumulant of the distribution, is the most well-known measure of skewness. When γ > 0 (γ < 0) the distribution is positively (negatively) skewed and will have a longer (shorter) right tail and a shorter (longer) left tail. If the distribution is symmetrical, γ equals zero. However, there are in [2] (Exercise 3.26) asymmetrical distributions with as many zero-odd order central moments as desired, so, the value of γ must be interpreted with some caution.
In the statistical literature, there is not a closed-form for the skewness coefficient of γ of the MLE in several regression models. Ref. [3] obtained a general n −1/2 γ expression (say γ 1 ) for the distribution of the MLE, where n is the sample size. Following [3], several works have been developed in order to obtain the γ 1 coefficient. In the first, Ref. [4] determined its expression for the class of generalized linear models and, the last one, Ref. [5] defined the γ 1 for the varying dispersion beta regression model and showed that this coefficient for the distribution of the MLE of the precision parameter is relatively large in small to moderate sample sizes. This paper is the first focused on a censored model.
In this work, we derive the γ 1 coefficient of the distribution of the MLE of the linear parameters in the Weibull censored data, assuming σ known, as σ = 1/2 and 1, the Rayleigh and exponential models, respectively. We discuss the situation when σ is unknown, however, it can be replaced by a consistent estimator, and then we can turn back to the original situation. This type of procedure was performed, for instance, by [4].
The remainder of the paper is organized as follows. Section 2 defines the Weibull censored data. In Section 3, we obtain a simple matrix expression, of order n −1/2 , for the skewness coefficients of the distributions of the MLEs of the linear regression parameters. In Section 4, some Monte Carlo simulations are performed. Two applications are presented in Section 5. Concluding remarks are offered in Section 6.

The Weibull Censored Data
We say that a continuous random variable T has Weibull distribution with scale parameter θ and shape parameter σ, or T ∼ WE(θ, σ), if its probability density function (pdf) is given by with t > 0, σ > 0 and θ > 0. From (1), we can observe two particular distributions: the exponential and the Rayleigh, where σ = 1 and σ = 1/2, respectively. In lifetime data, there is the censoring restriction, i.e, if T 1 , . . . , T n are a random sample from (1), instead of T i , we observe, under right censoring, t i = min(T i , L i ), where L i is the censoring time, independent of T i , i = 1, . . . , n. In this work, we consider an hybrid censoring scheme, where the study is finalized when a pre-fixed number, r ≤ n, out of n observations have failed, as well as when a prefixed time, say L 1 = . . . = L n = L, has been reached. The type I censoring is a particular case for r = n and the type II censoring appears when L 1 , . . . , L n = +∞. Additionally, we add the non-informative censoring assumption, i.e., the random variables L i does not depend on θ. Under this setup, the log-likelihood function has the form Usually, the regression modeling considers the distribution of Y i = log(T i ) instead of T i . The distribution of Y i is of the extreme value form with pdf given by where µ i = log θ i . The regression structure can be incorporated in (2) by making θ i = exp x i β , where β is a p-vector of unknown parameters and x i is a vector of regressors related to the ith observation. From this moment, we assume that σ is known, then, the log-likelihood function derived from (2) is given by The total score function and the total Fisher information matrix for β are, respectively, U β = σ −1 X W 1/2 v and K ββ = σ −2 X W X, where X = (x 1 , . . . , x n ) , the model matrix, assuming It can observed that the value of w i depends on the mechanism of censoring. That means denotes the rth order statistic from W 1 , . . . , W n and q = P W (r) ≤ log L i . Note that q = 1 and q = 0 for types I and II censoring, respectively. The proof is presented in the Appendix A. The MLE of β, β, is the solution of U β = 0. The β can not be expressed in closed-form. It is typically obtained by numerically maximizing the log-likelihood function using a Newton or quasi-Newton nonlinear optimization algorithm. Under mild regularity conditions and in large samples, approximately.

Skewness Coefficient
As discussed, the skewness coefficient is simple way to verify whether the approximation to normality is adequate. The model presented in (2) does not has a closed-form for this coefficient. The alternative is to apply the [3] result. These authors derived an approximation of order O(n −2 ) for the third cumulant of the MLE of the a-th regressor, i.e., . . , p, which can be expressed as where m Here, ∑ represents the summation over all combinations of parameters and over all the observations. From (3), after some algebra, we can express the third cumulant of the distribution ofβ for the Weibull censored data as where represents a direct product of matrices and 1 is a n-dimensional vector of ones. Finally, by (4) and the Fisher information matrix, the asymmetry coefficient of the distribution of β to order n −1/2 is given by for type II censoring, W = 0, then (5) reduces to γ 1 . More details about the involved expressions are presented in Appendix A. The study of asymptotic properties of the Weibull censored data was the goal of many papers. Refs. [6,7] derived the Bartlett and the Bartlett-type correction factors for likelihood ratio and score tests, respectively, for the exponential censored data. Ref. [8] generalized these previous for the Weibull censored data and also derived the Bartlett-type correction factors for the gradient test. Ref. [9] presented the asymptotic expansions up to order n −1/2 of the non null distribution functions of the likelihood ratio, Wald, Rao score and gradient statistics also for the censored exponential data. The result in expression (5) can be incorporated in this gallery.

Simulation Study
In this section, we compare the sample skewness coefficient (ρ) and the n −1/2 skewness coefficients evaluated in the true and estimated parameters ( γ 1 and γ 1 , respectively) of the distributions of the MLEs in the Weibull censored model. To draw the data, we consider three values for σ: 0.5, 1 and 3; five sample sizes: 20, 30, 40, 60 and 100; three values for the percent of censoring C: 10%, 25% and 50%; and two number of regressors p: 3 and 5, where we consider two vectors for β in each case: (−2, 0.5, 1) and (1, −0.75, 0.5) for p = 3 and (−2, 0.5, 1, −0.3, −0.5) and (1, −0.75, 0.5, −1, 0.8) for p = 5. For each combination of σ, β, % of censoring and sample size we considered 20,000 Monte Carlo replicates. Each vector of covariates x i considers an intercept term and the p − 1 remaining covariates were drawn independently from the standard normal distribution. Values from the Weibull model are drawn considering the inverse transformation method. Therefore, the greater n × C/100 values were censored at the observed (1 − C/100)-th quantile (a type II censoring scheme). For each sample, we considered the jackknife estimator for σ, say σ J . Therefore, the computation of γ 1 and γ 1 was performed considering (β, σ) and ( β, σ J ), the true and estimated parameters, respectively. Additionally, ρ is computed based on the 20,000 (marginal) skewness coefficient for the components of β. Table 1 summarizes the case β = (−2, 0.5, 1) (with p = 3 regressors) and C = 10%. The main conclusions are the following:  • The terms γ 1 and γ 1 are closer in all the considered combinations, suggesting that γ 1 approaches γ 1 in a reasonable way, even when the sample size is small.

•
In general terms, γ 1 approaches well ρ for β 1 and β 2 . However, for β 0 the terms seem discrepant even for n = 100. This suggest that, as expected, when n increases the skewness coefficient of the MLE estimators for the components of β will be more symmetric.
Results suggest that, even with a moderate percentage of censored observations and small sample sizes, the distribution of the MLE for the components of β in the Weibull censored model are closer to the symmetry. The combinations of β, p and C not seem to affect the results. A simulation study showing this finding was omitted for the sake of brevity.

Applications
In this section we illustrate with two real dataset the application of the estimated skewness coefficient for the MLE estimators in the Weibull censored regression model. All the routine was performed in the statistical software R, [10]. Codes can be found in the personal website from the first author https://www.ufjf.br/tiago_magalhaes/downloads/.

Smokers Dataset
This dataset is related to a clinical trial on the effectiveness of triple-combination pharmacotherapy for tobacco dependence treatment conducted by the Cancer Institute of New Jersey and Robert Wood Johnson Foundation. The trial recruited 127 smokers 18 years or older with predefined medical illnesses from the local community. The outcome were the time (in days) to first relapse (return to smoking). The study lasted 182 days (26 weeks). Therefore, the times are subject to a censoring type I (32% of times were censored). We only considered the 113 patients where such observed time was positive (non-zero). Other measures were assigned randomly treatment group with levels combination or patch only (grp), age in years at time of randomization (age) and employment (full-time or non-full-time). We consider that time i ∼WE(θ i ; σ), where log θ i = X i β, β = (β intercept , β grp , β age , β employment ) and We estimated σ J = 1.617008 based on the jackknife method, which was used as known in all the computations. Table 2 shows the parameters estimates, their standard errors and the estimated skewness coefficients and Figure 1 shows the estimated density function based on 1000 bootstrap samples for the coefficients related to the covariates grp, age and employment. Note that the estimated skewness for all parameters were closer to zero, suggesting a symmetric distribution for the estimators which is corroborated by the estimated density based on the bootstrap.

Insulating Fluids Dataset
This dataset was presented in [11] on insulating fluids and it is related an accelerated test performed in order to determine the relationship between time (in minutes) to breakdown and voltage (in kilovolts). The authors assumed a regression structure based on the Weibull model and a common censoring time at L = 200 (type I censoring), i.e., time i ∼WEI(θ i , σ), where log θ i = X i β, i = 1, . . . , 76, X i = (β Intercept , β log-voltage ). We estimated σ J = 1.296704 based on the jackknife method, which was used as known. Table 3 shows the estimates, standard errors and estimated skewness coefficient for the MLE estimators and Figure 2 shows the estimated density function for β Intercept and β log-voltage . Newly, the estimated skewness for both parameters are closer to zero, suggesting a symmetric distribution for the estimators as also suggest the estimated density based on bootstrap.

Concluding Remarks
Since its beginning, the Weibull distribution and regression showed how it is important. In the frequentist context, this model depends strongly on the asymptotic properties of the MLE. Here, we presented an expression of the skewness that, in practical applications, can be used as an indicator of departure from the normal distribution of the MLE. Although the expression (3) entails a great deal of algebra, the final expression (5) of the skewness of the MLE distribution has a very nice form only involving simple operations on diagonal matrices and can be easily implemented into a statistical software, for instance, R, [10].
Author Contributions: All the authors contributed significantly to the present paper. Acknowledgments: We also acknowledge the referee's suggestions that helped us to improve this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
In this Section we provided some additional details related to the computation of the manuscript.