key: cord-0073234-uidevlfn
authors: Bothma, E.; Allison, J. S.; Visagie, I. J. H.
title: New classes of tests for the Weibull distribution using Stein’s method in the presence of random right censoring
date: 2022-01-09
journal: Comput Stat
DOI: 10.1007/s00180-021-01178-0
sha: e37bd86568f3eeb187b673409ec649020dc9b0dc
doc_id: 73234
cord_uid: uidevlfn

We develop two new classes of tests for the Weibull distribution based on Stein’s method. The proposed tests are applied in the full sample case as well as in the presence of random right censoring. We investigate the finite sample performance of the new tests using a comprehensive Monte Carlo study. In both the absence and presence of censoring, it is found that the newly proposed classes of tests outperform competing tests against the majority of the distributions considered. In the cases where censoring is present we consider various censoring distributions. Some remarks on the asymptotic properties of the proposed tests are included. We present another result of independent interest; a test initially proposed for use with full samples is amended to allow for testing for the Weibull distribution in the presence of censoring. The techniques developed in the paper are illustrated using two practical examples.

The Weibull distribution is often used in survival analysis as well as reliability theory, see e.g., Kalbfleisch and Prentice (2011) . This flexible distribution is a popular model which allows for constant, increasing and decreasing hazard rates. The Weibull distribution is also frequently applied in various engineering fields, including electrical and industrial engineering to represent, for example, manufacturing times, see Jiang and Murthy (2011) . As a result of its wide range of practical uses, a number of goodness-of-fit tests have been developed for the Weibull distribution; see e.g, Mann et al. (1973) , Tiku and Singh (1981) , Liao and Shimokawa (1999) , Cabaña and Quiroz (2005) as well as Krit (2014 The papers listed above deal with testing for the Weibull distribution in the full sample case; i.e., where all lifetimes are observed. However, random right censoring often occurs in the fields mentioned above. For example, we may study the duration that antibodies remain detectable in a patient's blood after receiving a specific type of Covid-19 vaccine, i.e. the duration of the protection that the vaccine affords the recipient. When gathering the relevant data, we will likely not be able to measure this duration in all of the patients. For example, some may leave the study by emigrating to a different country while still having detectable antibodies. In this case, the exact time of interest is not observed. This situation is referred to as random right censoring, see e.g., Cox and Oakes (1984) .

In the presence of censoring, testing the hypothesis that the distribution of the lifetimes is Weibull is complicated by the fact that an incomplete sample is observed. Balakrishnan et al. (2015) suggests a way to perform the required goodness-of-fit tests by transforming the censored sample to a complete sample. Another approach is to modify the test statistics used in the full sample case to account for the presence of censoring. Although fewer in number, tests for the Weibull distribution in the presence of random censoring are available in the literature. For example, Koziol and Green (1976) and Kim (2017) propose modified versions of the Cramér-von Mises test and the test proposed in Liao and Shimokawa (1999) , respectively, for use with censored data.

Throughout this paper we are primarily interested in the situation where censoring is present; the results relating to the full sample case are treated as special cases obtained when all lifetimes are observed. Before proceeding some notation is introduced. Let X 1 , . . . , X n be independent and identically distributed (i.i.d.) lifetime variables with continuous distribution function F and let C 1 , . . . , C n be i.i.d. censoring variables with distribution function H , independent of X 1 , . . . , X n . We assume non-informative censoring throughout. Let T j = min(X j , C j ) and

Note that in the full sample case T j = X j and δ j = 1 for j = 1, . . . , n.

Based on the observed pairs (T j , δ j ), j = 1, . . . n we wish to test the composite hypothesis

for some unknown λ > 0 and θ > 0. Here X ∼ W eibull(λ, θ ) refers to a Weibull distributed random variable with distribution function

This hypothesis is to be tested against general alternatives. We will make use of maximum likelihood estimation to estimate λ and θ . The log-likelihood of the Weibull distribution is

where d = n j=1 δ j . In the full sample case d = n. No closed form formulae for the maximum likelihood estimatesλ andθ exist, meaning that numerical optimisation techniques are required to arrive at parameter estimates.

Since the Weibull distribution has a shape parameter the first step for many goodness-of-fit tests for the Weibull distribution is to transform the data. If X ∼ W eibull(λ, θ ) then a frequently used transformation is log(X ), which results in a random variable that is type I extreme value distributed with parameters log(λ) and 1/θ . The resulting transformed random variable is part of a location scale family, which is a desirable result when performing goodness-of-fit testing. We therefore have that if X ∼ W eibull (λ, θ ) , then X (t) = θ(log(X ) − log(λ)) follows a standard type I extreme value distribution with distribution function

We denote a random variable with this distribution function by E V (0, 1). As a result the hypothesis in (1) holds, if, and only if, X (t) ∼ E V (0, 1). All of the test statistics considered make use of the transformed observed values

withλ andθ the maximum likelihood estimates of the Weibull distribution. Let

If X 1 , . . . , X n are realised from a Weibull(λ, θ ) distribution, then X (t) 1 , . . . , X (t) n will approximately follow an E V (0, 1) distribution see e.g. Kotz and Nadarajah (2000) . The resulting random variables are no longer independent. However, several properties of the classical testing procedures remain unaffected when performing this type of transformation; the interested reader is referred to Baringhaus and Henze (1991) as well as Gupta and Richards (1997) and the references therein for more details. The tests employed below are based on discrepancy measures between the calculated values of Y 1 , . . . , Y n and the standard type I extreme value distribution. The order statistics of Y 1 , . . . , Y n are denoted by Y (1) < · · · < Y (n) , while δ ( j) represents the indicator variable corresponding to Y ( j) .

The remainder of the paper is structured as follows. In Sect. 2, we propose two new classes of tests for the Weibull distribution for both the full sample and censored case. We also modify the test proposed in Krit (2014) to accommodate random right censoring. In the presence of censoring the null distribution of all the test statistics considered depend on the unknown censoring distribution, we therefore propose a parametric bootstrap procedure in Sect. 2 in order to compute critical values. Section 3 presents the results of a Monte Carlo study where the empirical powers of the newly proposed classes of tests as well as the newly modified test are compared to those of existing tests. The paper concludes in Sect. 4 with two practical applications; one concerning the survival times of patients diagnosed with a certain type of leukemia (no censoring is present in these data) and the other relates to observed leukemia remission times (in the presence of censoring). Some avenues for future research are also discussed.

Our newly proposed classes of tests are based on the following theorem, which characterises the standard type I extreme value distribution.

Theorem 1 Let W be a random variable with absolutely continuous density and assume that E e W < ∞. In this case W ∼ E V (0, 1), if, and only if, 

where the last step follows from the fact that f must be a density function and hence integrate to 1. This completes the proof.

Let w(t) be a non-negative, symmetric weight function. From Theorem 1, we have that

Note that the inclusion of the weight function, w, above is required to ensure that η is finite. Clearly η will be unknown because G is unknown. However, G can be estimated by the Kaplan-Meier estimator, G n , of the distribution function given by

More details about this estimator can be found in Kaplan and Meier (1958) , Efron (1967) as well as Breslow and Crowley (1974) . In the full sample case this estimator reduces to the standard empirical distribution function, G n (X ( j) ) = j/n. Let Δ j denote the size of the jump in G n (T ( j) );

Simple calculable expressions for the Δ j 's are

In the full sample case Δ j = 1/n, j = 1, . . . , n.

Estimating G by G n in (3), we propose the test statistic

where w a (t) is a weight function containing a user-defined tuning parameter a > 0. The null hypothesis in (1) is rejected for large values of S n,a .

Straightforward algebra shows that, if w a (t) = e −at 2 , then the test statistic simplifies to

and if w a (t) = e −a|t| the test statistic has the following easily calculable form

New goodness-of-fit tests containing a tuning parameter are often accompanied by a recommended value of this parameter; this choice is typically based on the finite sample power performance of the test. Another approach which may be used is to choose the value of the tuning parameter data-dependently; see e.g., Allison and Santana (2015) .

In this paper, we opt to use the values recommended in the literature for tests containing a tuning parameter. The weight functions specified above correspond to scaled Gaussian and Laplace kernels. These weight functions are popular choices found in the goodness-of-fit literature; see e.g., Meintanis and Iliopoulos (2003) , , Betsch and Ebner (2018) , Betsch and Ebner (2019) as well as Henze and Visagie (2020) . The popularity of these weight functions is, at least in part, due to the fact that their inclusion typically results in simple calculable forms for L 2 -type statistics which do not require numerical integration. As an alternative to the weight functions used, one may employ a symmetric uniform kernel as a weight function; see e.g., Fernández et al. (2008) . However, in the mentioned paper, the authors found that the computational time required for the test statistic obtained using the symmetric uniform kernel was substantial. This, coupled with the simple computational forms obtained using the weight functions defined above and the favourable power performance discussed in Sect. 3, motivated us to restrict our attention to the scaled Gaussian and Laplace kernels.

Although we do not derive the asymptotic results related to the proposed classes of test statistics we include some remarks in this regard. S n,a is a characteristic function based weighted L 2 -type statistic. The asymptotic properties of this class of statistics, in the complete sample case, are studied in detail in Feuerverger and Mureika (1977) , while more recent references include Baringhaus and Henze (1988) , Klar and Meintanis (2005) as well as Baringhaus et al. (2017) . A convenient setting for the derivation of the asymptotic properties of these tests is the separable Hilbert space of square integrable functions. Typically, the asymptotic null distribution of S n,a corresponds to that of ∞ −∞ |Z (t)| 2 w a (t)dt =: S a , where Z (·) is a zero-mean Gaussian process. The distribution of S a is the same as that of ∞ j=1 λ j U j , where U j are i.i.d. chi-squared random variables with parameter 1 and where λ j are eigenvalues of an integral operator (see e.g., Allison et al. (2021) ). These tests are consistent against a large class of fixed alternative distributions. Additionally, these tests are frequently consitent against contiguous alternatives converging to the null at a rate of n −1/2 .

In the case of random censoring, very little asymptotic results are available in the literature for test statistics of this type. Very recently, advances have been made in this regard; see e.g., Cuparić and Milošević (2021) , in which a test for exponentiality is considered based on so-called inverse probability censoring weights, where the authors derive the asymptotic properties of the given test. In addition Fernández and Rivera (2020) studied Kaplan-Meier U -and V -statistics in order to derive some asymptotic results relating to the lifetime distribution. Some of these results may be helpful in deriving the asymptotic properties of the tests proposed in this paper in future research.

The null distribution of each of the test statistics considered depends on the unknown censoring distribution, even in the case of a simple hypothesis, see D' Agostino and Stephens (1986) . Since we will not assume any known form of the censoring distribution, we propose the following parametric bootstrap algorithm to estimate the critical values of the tests.

1. Based on the pairs (T j , δ j ), j = 1, . . . , n estimate θ and λ byθ andλ, respectively, using maximum likelihood estimation. 2. Transform T j to Y j using the transformation in (2) for j = 1, . . . , n.

3. Calculate the test statistic, say W n := W (Y 1 , . . . , Y n ; δ 1 , . . . , δ n ). 4. Obtain a parametric bootstrap sample X * 1 , . . . , X * n by sampling from a Weibull distribution with parametersθ andλ. 5. Obtain a non-parametric bootstrap sample, C * 1 , . . . , C * n , by sampling from the Kaplan-Meier estimate of the distribution of C

The algorithm provided above is quite general and can easily be amended in order to test for any lifetime distribution in the presence of random censoring. In the absence of censoring there is no need to implement this algorithm; in this case, the critical values can be obtained via Monte Carlo simulation by sampling from any Weibull distribution and effecting the transformation discussed above.

In this section, we compare the power performances of the newly proposed tests to those of existing tests via a Monte Carlo simulation study. The existing tests used include the classical Kolmogorov-Smirnov (K S n ) and Cramér-von Mises (C M n ) tests. These tests have been modified for use with censored data, see Koziol and Green (1976) . The test introduced in Liao and Shimokawa (1999) is considered in the case of full samples. A modification making this test suitable for use with censored data is proposed in Kim (2017) ; we denote the test statistic by L S n in both the full sample and censored cases. The calculable forms of the test statistics mentioned above are

Krit (2014) proposes a test for the Weibull distribution in the full sample case. This test compares the empirical Laplace transform of the random variables resulting from the transformation in (2) to the Laplace transform of an E V (0, 1) random variable; ψ(t) = Γ (1 − t) for t < 1. Let ψ n be the empirical Laplace transform of the transformed observations, obtained using the Kaplan-Meier estimate of the distribution function;

The resulting test statistic is

where w a (t) = e at−e at is a weight function, a a user-specified tuning parameter and I some interval. Based on numerical considerations, Krit (2014) suggests that I = (−1, 0] should be used. The quantity in (6) can be approximated by a Riemann sum;

where m is the number of points at which the integrand is evaluated. In the numerical results shown below, we use a = −5 and m = 100 as recommended in Krit (2014) .

In the full sample case, G n , in (5), is taken to be the empirical distribution function. Upon setting Δ j = 1/n we obtain the test statistic in Krit (2014) . For each of the tests considered above, the null hypothesis in (1) is rejected for large values of the test statistics.

In the numerical results presented below, we use a nominal significance level of 10% throughout. Empirical powers are presented for sample sizes n = 50 and n = 100. The empirical powers for complete and censored samples are reported; censoring proportions of 10% and 20% are included. For each lifetime distribution considered, we report the powers obtained using three different censoring distributions. The first censoring distribution used is the exponential distribution, the parameter of which is chosen so as to obtain the specified level of censoring. The second censoring distribution used is the uniform distribution with support (0, m); again, m is chosen such that the required censoring level is achieved. The final censoring distribution used is the Koziol-Green model proposed in Koziol and Green (1976) . Denote the survival function of a given lifetime distribution by S. In this case, the Koziol-Green censoring distribution, indexed by β, has survival function S β (t). It can be shown that the censoring proportion is β/(1 + β). We chose the value of β so as to ensure that the required level of censoring is achieved. The alternative lifetime distributions considered are listed in Table 1 . Note that each of the alternatives considered have the same support as that of the Weibull distribution with the exception of the beta distribution which is restricted to the unit interval (0, 1) (Tables 10, 11, 12 and 13).

The obtained empirical powers are presented in Tables 2, 3 , 4, 5, 6, 7, 8 and 9. These tables report the percentages of 50 000 independent Monte Carlo samples that Weibull Lind(2) 10 10 11 11 11 11 11 10 10 11 Table 3 Heatmap of the estimated powers for the full sample case where n = 50 lead to the rejection of the null hypothesis, rounded to the nearest integer. For ease of comparison, the highest power in each line is printed in bold. Tables 2 and 4 contain the results relating to full samples. In order to ease visual comparison of the results obtained, we include so called "heatmaps", see Döring and Cramer (2019), of these results in Tables 3 and 5 . For each test considered, Tables 6, 7, 8 and 9 show three empirical powers against each lifetime distribution, corresponding to the three different censoring distributions used. In each case, the results for the exponential, uniform and Koziol-Green models are shown in the first, second and third lines, respectively. In order to reduce the computational cost associated with the numerical powers a warp-speed bootstrap procedure, see Giacomini et al. (2013) , is employed. This methodology has been employed by a number of authors in the literature to compare n,a , we include numerical powers in the cases where a is set to 1, 2 and 5 in the full sample case and set to 0.75, 1 and 2 in the presence of censoring. The difference between these two sets of choices is due to the fact that the newly proposed tests exhibits higher powers in the presence of censoring if slightly smaller values of a are used. All calculations are performed in R, see [44] . The LindleyR package is used to generate samples from censored distributions, see Mazucheli et al. (2016) . Parameter estimation is performed using the parmsurvfit package, see Jacobson et al. (2018) , while the tables are produced using the Stargazer package, see Hlavac (2018) . (2) n,2 10 10 10 8 9 9 9 10 10 9 10 10 10 9 8 8 10 9 9 8 W (2) 9 1 0 1 0 8 9 9 9 9 9 9 10 10 10 9 9 9 10 10 9 9 9 1 0 1 0 9 8 9 9 9 9 8 Γ (2 

First, we consider the results associated with the full sample case, given in Tables 2  and 4 , together with the heatmaps shown in Tables 3 and 5. All of the tests considered attain the nominal size for both sample sizes used. The tests associated with the highest powers are S (1) n,2 and S (1) n,5 , although L S n and S (2) n,5 also performs well. When analysing complete samples, we recommend using S (1) n,2 or S (1) n,5 . We now turn our attention to the powers achieved in the presence of censoring. The size of the tests are maintained closely for all sample sizes for censoring proportions of 10% and 20%, with the single exception of K R m,a in the case of small sample sizes. As expected, the powers generally increase with sample size and decrease marginally as the censoring proportion increases. Lind(0.5) 10 11 12 8 1 0 1 0 1 1 9 9 9 10 11 12 8 1 0 1 0 12 9 9 1 0 10 10 12 8 9 10 11 9 9 9

Lind(2) 10 10 10 7 9 9 10 9 9 8 10 10 10 6 9 9 1 0 9 9 9 10 10 11 7 9 9 9 9 9 8 W (1) 10 10 10 10 8 8 9 9 9 8 10 10 10 11 8 8 10 9 9 8 10 10 10 10 9 9 9 10 9 9

W (1.5) 10 10 10 9 9 9 9 10 10 9 10 10 10 9 9 9 9 9 9 8 10 10 10 11 8 8 9 9 9 8 W (2) 10 10 10 8 9 9 9 10 10 9 Lind(2) 10 10 11 10 9 9 10 9 9 9 10 10 11 12 9 9 1 0 9 9 8 10 10 11 10 9 9 10 9 9 9

Comparing the results associated with a sample size of 50, we see that S

(1) n,1 , K R n and L S n generally tend to provide the highest powers. However, it should be noted that K R n achieves very low power against certain alternatives; notably against the beta distributions considered when the censoring distribution is the Kozoil-Green model. When considering the empirical powers associated with samples of size 100, S

(1) n,2 exhibits the highest powers, followed by S (1) n,1 and L S n . Lind(0.5) 10 11 13 9 9 9 9 9 8 8 10 11 13 6 1 0 9 9 9 9 9 10 11 13 6 9 8 8 8 8 8 When compiling the numerical results, we also considered a wider range of values for the tuning parameter a than those reported in the table. Although some power variation is evident when varying a, the powers achieved by the newly proposed classes of tests are not particularly sensitive to the choice of the tuning parameter a. However, based on the observed numerical powers, we recommend using S (1) n,2 when testing the hypothesis in question in the presence of censoring.

In this section, we apply the tests used in Sect. 3 to test the hypothesis in (1) based on two real-world data sets. The first data set, reported in Table 10 , contains the survival , 47, 58, 74, 177, 232, 273, 285, 317, 429, 440, 445, 455, 468, 495, 497, 532, 571, 579, 581, 650, 702, 715, 779, 881, 900, 930, 968, 1077, 1109, 1314, 1334, 1367, 1534, 1712, 1784, 1877, 1886, 2045, 2056, 2260, 2429, 2509 , 5, 8, 8, 9, 10, 10, 10, 10, 10, 11, 12, 12, 12 * , 13, 14, 20, 20 * , 23, 23, 25, 25, 25, 28, 28, 28, 28, 29, 31, 31, 31, 32, 37, 40, 41, 41, 48, 48, 57, 62, 70, 74, 75, 89, 99, 100, 103, 124, 139, 143, 159 * , 161 * , 162, 169, 190 * , 195, 196 * , 197 * , 199 * , 205 * , 217 * , 219 * , 220, 245 * , 258 * , 269 * Kotze and Johnson (1983) as well as . This data set is not subject to censoring, i.e. all lifetimes are observed. The second data set contains the initial remission times of leukemia patients, in days; for more details see Lee and Wang (2003) , this data set can be found in Table 12 . These data contains censored observations, indicated using an asterisk. The original data were segmented into three treatment groups. However, Lee and Wang (2003) showed that the data do not display significant differences among the various treatments. As a result we treat the data as i.i.d. realisations from a single, censored, lifetime distribution. All reported p-values are estimated using one hunderd thousand bootstrap replications; these results are displayed in Tables 11 and 13, respectively. From the results of the practical example in Table 11 it is clear that none of the tests reject the null hypothesis that the survival times after a leukemia diagnosis are Weibull distributed at the 5% or 10% levels of significance. As a result, we conclude that the Weibull distribution is an appropriate model for these data.

The results associated with the initial remission times, in Table 13 , indicate that K S n , C M n and L S n reject the hypothesis in (1) at a 5% significance level. However, none of the remaining 7 tests considered result in a rejection of the null hypothesis at the 5% or 10% levels. We conclude that the Weibull distribution is likely to be an appropriate model for the observed times. The data set under consideration was also analysed in Bothma et al. (2020) , where the null hypothesis of exponentiality of the remission time was strongly rejected. The mentioned paper recommended that a more flexible distribution be used when modelling these data. The results above indicate that the additional flexibility of the Weibull (compared to the exponential) distribution indeed ensures that the Weibull distribution is a more appropriate model than the exponential for the initial remission times considered.

A number of interesting numerical phenomena are evident when considering the powers of the various tests. It is clear that the achieved powers and, therefore, the null distribution of the test statistic, is influenced by the shape of the censoring distribution. The effect of the censoring distribution on the critical values of the tests seem not to have been investigated in the literature to date. Some authors perform goodness-offit testing by enforcing a parametric assumption on the censoring distribution, see e.g., Kim (2017) . An additional consideration that seems to have been neglected in the literature is the effect on the null distribution of the test statistic, and hence the power of the test, of a specific assumption made in the Kaplan-Meier estimate of the distribution function. Some authors, in order to ensure that G n satisfies the requirements of a distribution, defines G n (t) = 1 for all t > X (n) regardless of whether or not the sample maximum is censored. We are currently investigating these open questions.

Funding Not applicable.

Availability of data and material All data used are in the public domain.

The authors declare that they have no conflict of interest.

Distribution-free goodness-of-fit tests for the pareto distribution based on a characterization

On a data-dependent choice of the tuning parameter appearing in certain goodness-of-fit tests

New weighted L 2 -type tests for the inverse Gaussian distribution

Testing the adequacy of semiparametric transformation models

An "apples-to-apples" comparison of various tests for exponentiality

An empirical analysis of some nonparametric goodness-of-fit tests for censored data

The limit distribution of weighted L2-goodness-of-fit statistics under fixed alternatives, with applications

A consistent test for multivariate normality based on the empirical characteristic function

A class of consistent tests for exponentiality based on the empirical laplace transform

Testing normality via a distributional fixed point property in the Stein characterization

A new characterization of the gamma distribution and associated goodness-of-fit tests

Kaplan-Meier based tests for exponentiality in the presence of censoring

A large sample study of the life table and product limit estimates under random censorship

Using the empirical moment generating function in testing for the Weibull and the type I extreme value distributions

New characterization-based exponentiality tests for randomly censored data

On the power of goodness-of-fit tests for the exponential distribution under progressive Type-II censoring

The two sample problem with censored data

A test for the two-sample problem based on empirical characteristic functions

The empirical characteristic function and its applications

A warp-speed method for conducting Monte Carlo experiments involving bootstrap estimators

Invariance properties of some classical tests for exponentiality

Testing for normality in any dimension based on a partial differential equation involving the moment generating function

stargazer: well-formatted regression and summary statistics tables

parmsurvfit: Parametric Models for Survival Data

A study of Weibull shape parameter: properties and significance

Nonparametric estimation from incomplete observations

Goodness-of-fit tests for randomly censored Weibull distributions with estimated parameters

Tests for normal mixtures based on the empirical characteristic function

New York Koziol JA, Green SB (1976) A Cramér-von Mises statistic for randomly censored data

Goodness-of-fit tests for the Weibull distribution based on the Laplace transform

New York Liao M, Shimokawa T (1999) A new goodness-of-fit test for type-i extreme-value and 2-parameter Weibull distributions with estimated parameters

A new goodness-of-fit test for the two-parameter Weibull or extreme-value distribution with unknown parameters

LindleyR: The Lindley Distribution and Its Modifications

Tests of fit for the Rayleigh distribution based on the empirical Laplace transform

Testing for serial independence in vector autoregressive models

An overview of goodness-of-fit tests for the Poisson distribution

R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing

Testing the two parameter Weibull distribution

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

W (1.5) 9 1 0 1 0 1 0 8 8 9 9 9 8 9 10 9 12 9 9 10 9 9 9 9 9 9 1 1 8 8 9 9 8 8 W (2) 10 10 10 8 9 9 9 9 9 8 10 10 10 10 8 8 9 9 9 8