A simulation-based goodness-of-fit test for survival data

Statistics & Probability Letters 47 (2000) 403 – 410

A simulation-based goodness-of- t test for survival data Gang Lia;∗;1 , Yanqing Sunb a Department

of Biostatistics, School of Public Health, University of California, Los Angeles, CA 90095-1772, USA b University of North Carolina at Charlotte, USA Received June 1999

Abstract To check the validity of a parametric model for survival data, a number of supremum-type tests have been proposed in the literature using Khmaladze’s (1993, Ann. Statist. 18, 582– 602) transformation of a test process. However, such a transformation is usually very complicated and lacks a clear interpretation. Information could also be lost through transformation. In this note, we propose a simulation-based supremum-type test directly from the original test process using an idea originally introduced by Lin et al. (1993, Biometrika 80, 557–572). The test is developed under the framework of Aalen’s (1978, Ann. Statist. 6, 701–726) multiplicative intensity counting process model, and therefore applies to a number of survival models including those with very general forms of censoring and truncation. By comparing the observed test process with a set of simulated realizations of an approximating process, our method can be used as a graphical tool as well as a formal test for checking the adequacy of the assumed parametric model. We establish consistency of the resulting test under any xed alternative. Its performance is investigated in a simulation study. Illustrations are given using c 2000 Elsevier Science B.V. All rights reserved some real data sets. Keywords: Censoring; Kolmogorov–Smirnov test; Monte Carlo method; Product-limit estimator; Truncation; Weak convergence

1. Introduction Let (N1 ; : : : ; Nn ) be a multivariate counting process with intensity process (1 ; : : : ; n ). Assume the Aalen (1978) multiplicative intensive model i (t) = Yi (t)h(t);

i = 1; : : : ; n;

(1)

where h(t) is a hazard function and Yi (t) is a predictable {0; 1}-valued process, indicating (by the value 1) that the ith subject is at risk at time t. Model (1) is known to encompass a wide range of models in survival analysis. An important example is the random censorship model in which one observes i.i.d. copies {X˜ i ; i }; i = 1; : : : ; n, of X˜ = min(X; C);  = I (X 6C), where X is the survival time with distribution ∗

Corresponding author. Tel.: +1-310-206-5865; fax: +1-310-267-2113. E-mail address: [email protected] (G. Li)

1

The research of Gang Li was partially supported by a National Institute of Health grant.

c 2000 Elsevier Science B.V. All rights reserved 0167-7152/00/$ - see front matter PII: S 0 1 6 7 - 7 1 5 2 ( 9 9 ) 0 0 1 8 6 - 8

404

G. Li, Y. Sun / Statistics & Probability Letters 47 (2000) 403 – 410

F(t); C is the censoring time with distribution G(t), and X and C are independent. Let Ni (t) = I (X˜ i 6t; i = 1) and Yi (t) = I (X˜ i ¿t); i = 1; : : : ; n. Then (N1 ; : : : ; Nn ) is a multivariate counting process with Ni having intensity process i (t) = Yi (t)h(t), where h(t) is the hazard function of X ; see Aalen (1978) and Gill (1980). Another useful example is the random truncation model. A randomly left-truncated data consists of n i.i.d. draws (X1∗ ; T1∗ ); : : : ; (Xn∗ ; Tn∗ ), from the conditional distribution of (X; T ), given that X ¿ T . Here X is the survival time with hazard function h(t) and T is a truncation variable. Let Ni (t) = I {Ti∗ ¡ Xi∗ 6t} and Yi (t) = I {Ti∗ ¡ t6Xi∗ }. Then (N1 ; : : : ; Nn ) is a multivariate counting process with Ni having intensity process i (t) = Yi (t)h(t); see Keiding R t and Gill (1990) and Aalen and Johansen (1978). Let F(t) = 1 − exp{− 0 h(s) ds} be the distribution function. The aim of this paper is to develop a supremum-type goodness-of- t procedure for testing the null hypothesis H0 that F belongs to a parametric family F = {F(·; );  ∈ }, where  is an open set in the Euclidean space Rq ; q¿1. There is an extensive literature on the general problem of testing t of a parametric family. Chi-square-type tests, wherein the time space is divided into cells and in each cell the number of observed events is compared with what is expected under the null hypotheses, have been developed by various of authors under dierent censoring schemes. Here we mention Moore and Spruill (1975), Habib and Thomas (1986), Akritas (1988), Kim (1993), Mihalko and Moore (1980), Hjort (1990) and Li and Doss (1993) among others. The chi-square-type tests are simple and easy to use, but involves subjective choice of a partition. Grouping data could also lead to loss of information. The Kolmogorov–Smirnov supremum-type tests, originated from Kolmogorov (1933) and Doob (1949), have also been studied by a number of authors including Khmaladze (1993), Hjort (1990), Andersen et al. (1992), Nikabadze and Sun (1997). A supremum-type test is typically based on certain test √ and Stute (1997) ˆ where Fn is a nonparametric estimate of F and ˆ is an estimate of . process such as n{Fn (t) − F(t; )}, However, the limiting distribution of such a test process is usually intractable. For this reason, supremum-type tests developed in the past were based on some transformation of the test process using Khmaladze’s (1993) idea. Although the transformation approach produces asymptotically valid goodness-of- t tests, it lacks a clear statistical interpretation. It also does not provide information on the pattern of departure from the assumed model since such features are lost through transformation. In this note, we propose an alternative approach for deriving supreme-type tests.√Instead of constructing a ˆ Using transformation, we develop a supremum test directly from the original test process n{Fn (t) − F(t; )}. an idea from Lin et al. (1993) for checking the Cox model, we approximate the null distribution of the test process by generating a large number of sample paths from a zero-mean Gaussian process. The observed pattern of the test process is then compared with the simulated sample paths of the approximating process. This leads to a formal goodness-of- t test as well as a graphical tool for checking the adequacy of the assumed parametric model. In Section 2 we describe the simulation-based goodness-of- t test. In Section 3, we show that the proposed test is consistent. In Section 4 we illustrate our method using some real data sets. We also present a simulation study to investigate its performance for nite samples. 2. The test Pn Pn De ne N (t) = i=1 Ni (t) and Y (t) = i=1 Yi (t). We assume hereafter that Y (t)=n converges uniformly to some limit y(t) with probability 1. For example, y(t) = {1 − F(t−)}{1 − G(t−)} for the random censorship model and y(t) = P(T ¡ t6X ) for the random left-truncation model. Let  Y  dN (s) 1− Fn (t) = 1 − Y (s) 06s6t

G. Li, Y. Sun / Statistics & Probability Letters 47 (2000) 403 – 410

405

be the product-limit estimator of F (Gill and Johansen, 1990). Let ˆ be the maximum likelihood estimator of  obtained by maximizing the parametric log-likelihood Z ln () =

∞

0

{log h(t; ) dN (t) − h(t; )Y (t) dt};

where h(t; )=f(t; )={1−F(t; )} and f(t; ) are the hazard function and density function of F(t; ), respectively. De ne the following test process ˆ Zn (t) = n1=2 {Fn (t) − F(t; )};

t¿0:

Clearly large deviations from 0 would indicate evidence against H0 . The null limiting distribution of Zn has been well understood in the literature. For example, Sun (1997, Theorem 3.1) showed that under H0 ; Zn converges weakly to a mean zero Gaussian process Z in the Skorohod space D[0; ] for some  ¿ 0, provided that certain regularity conditions hold. The covariance function of Z is given by Z cov{Z(t1 ); Z(t2 )} = S(t1 )S(t2 )

0

min{t1 ; t2 }

h(s) ds − {F0 (t1 )}T I −1 ()F0 (t2 ); y (s)

where S = 1 − F; y (t) denotes the limit of R ∞Y (t)=n under H0 ; I () = (Iij ()) is the information matrix for  with the (i; j)th element given by Iij () = 0 h0i (s)h0j (s)y (s) ds=h(s; ), and F0 (t) and h0 (t) are the vectors of partial derivatives of F(t; ) and h(t; ) with respect to . However, it is dicult to derive a Kolmogorv– Smirnov-type test directly from the weak convergence result since the distribution of supt∈[0; ] |Z(t)| is intractable. To approximate the distribution of Zn , we rst note that, under H0 ; Zn (t) has the following martingale representation: Z Zn (t) = S(t)

t

n1=2

Pn

0

dMi (s) − n−1=2 {F0 (t)}T I −1 Y (s) i=1

Z

∞



0

h0 (s) h(s; )

X n

dMi (s) + op (1);

(2)

i=1

Rt where Mi (t)=Ni (t)− 0 Yi (s)h(s; ) ds; i =1; : : : ; n, are orthogonal local square integrable martingales. This can ˆ and an application of Taylor series expansion; be obtained from the martingale representations of Fn and , see, e.g., Sun (1997). Based on (2), we de ne the following process: Zn∗ (t)

ˆ = S(t; )

Z

t

n1=2

0

Pn

i dNi (s) −1 − n−1=2 {Fˆ0 (t)}T Iˆ Y (s)

i=1

Z 0

∞

(

h0ˆ(s) ˆ h(s; )

)

n X

i dNi (s);

(3)

i=1

where 1 ; : : : ; n are i.i.d. standard normal random variables independent of the data, Fˆ0 and h0ˆ , are obtained by replacing  by ˆ in F 0 and h0 , respectively, and Iˆ = (Iˆij ) is given by 

2 1 l () @ n : Iˆij = − n @i @j =ˆ



406

G. Li, Y. Sun / Statistics & Probability Letters 47 (2000) 403 – 410

Note that Zn∗ (t) is obtained replacing S; I and  in (2) by their sample estimates, and by replacing Mi (t) with i Ni (t), an idea originally used by Lin et al. (1993) to approximate cumulative sums of martingale-based residuals under Cox’s model. The motivation of approximating Mi (t) by i Ni (t) is that Mi (t) has mean 0 and variance E{Ni (t)}. It can be shown that, conditional on the data, Zn∗ converges weakly to Z under H0 , where Z is limiting process of Zn described earlier. This can be proved by noting that Zn∗ is a Gaussian process whose covariance converges to that of Z and by verifying a tightness criterion of Billingsley (1968, p. 128). Thus, under H0 , the distribution of Zn can be approximated by the conditional distribution of Zn∗ given the data. To obtain an -level Kolmogorov–Smirnov test for H0 , we rst generate a large number, say L = 2000, of independent samples of size n from the standard normal distribution. Then, let cn () be the 100(1 − )th percentile of ∗ ∗ (t)|; : : : ; sup |ZnL (t)|; sup |Zn1

06t6

where

∗ Znl

06t6

is calculated using (3) based on the lth normal sample, l = 1; : : : ; L.

Reject H0 if

sup |Zn (t)| ¿ cn ():

(4)

06t6

By plotting Zn (t) against time t together with a reference set of sample paths of Zn∗ (t), one can examine how unusual the observed pattern of Zn (t) is, and conclude whether or not the parametric model F provides good t to the data. In Section 4, we give some illustrations using both simulated and real data. 3. Consistency We show that the test based on (4) is consistent in the sense that its rejection power tends to 1 as n → ∞ under any xed alternative Ha : F = F0 6∈ F. Test (4) is based on the fact that, under H0 , the process Zn∗ (t) has the same limiting process Z(t) as Zn (t). When H0 does not hold, Zn∗ (t) will generally not converge to Z(t). However, the following theorem shows that Zn∗ (t) still converges weakly to some Gaussian process under any xed alternative Ha : F = F0 6∈ F. Theorem 3.1. Let Mn () = ln ()=n and let Z ∞ {h0 (s) log h (s) − h (s)}y0 (s) ds; M0 () = EF0 {Mn ()} = 0

where h0 (s) is the hazard function of F0 (s) and y0 (s) is the limit of Y (s)=n under F = F0 . Let ˜ be a value that maximizes M0 (). Suppose that a:s:

sup |Mn () − M0 ()| → 0

∈

and

sup

˜ :k−k¿

˜ for any  ¿ 0: M0 () ¡ M0 ()

(5)

If F = F0 6∈ F; then; conditional on the data; d Zn∗ → Z˜

as n → ∞

for almost all sequences (N1 ; Y1 ); (N2 ; Y2 ); : : : , where Z˜ is a mean zero Gaussian process with covariance function given by Z min{t1 ; t2 } h0 (s) ˜ ˜ ˜ 2 )} = S(t1 ; )S(t ˜ 1 ); Z(t ds ; ) cov{Z(t 2 y0 (s) 0

G. Li, Y. Sun / Statistics & Probability Letters 47 (2000) 403 – 410

Z

˜ + {F0˜(t1 )}T I −1 () Z −

t1 0

(

h0˜(s)

h˜(s)

)T

0

∞

(

h0˜(s)

h˜(s) 

h0 (s) ds I

−1

)(

h0˜(s)

)T

h˜(s)

˜ 0˜(t2 ) − ()F 

407

 ˜ 0˜(t2 ) h0 (s)y0 (s) ds × I −1 ()F 

Z 0

t2

(

h0˜(s)

h˜(s)

)T



˜ 0˜(t1 ): h0 (s) ds I −1 ()F 

Proof. Similar to Van der Vaart (1998, Theorem 5:7), it can be shown that under the conditions (5), ˆ converges to ˜ with probability 1. Using this fact, one can verify that the covariance function of Zn∗ (t) ˜ converges to that of Z(t) with probability 1. Weak convergence of Zn∗ (t) can then be established by noting ∗ that Zn (t) is a Gaussian process and by verifying a tightness criterion of Billingsley (1968, p. 128). A key to the above proof is that, even when H0 does not hold, ˆ still converges to some value ˜ ∈ . ˜ can be regarded as the projection of F0 in F. Here F(t; ) Corollary 3.1. Under Ha : F = F0 6∈ F;   P sup |Zn (t)| ¿ cn () → 1 06t6

as n → ∞. Proof. It follows from Theorem 3.1 that, if F = F0 6∈ F, then cn () converges to a nite limit that is equal ˜ On the other hand, to the 100(1 − )th percentile of sup06t6 |Z(t)|. ˆ Zn (t) = n1=2 {Fn (t) − F0 (t)} + n1=2 {F0 (t) − F(t; )}; where the rst term is known to converge weakly to a Gaussian process and the second term goes to ∞ with a:s: ˜ 6= F0 (t). Therefore, ˆ → F(t; ) probability 1 since F(t; )   P sup |Zn (t)| ¿ cn () → 1: 06t6

4. Numerical studies We rst illustrate the simulation-based test (4) with two real data sets. Then we present results from a small simulation study on its nite-sample performance in terms of signi cance level and rejection power. The rst data set came from a prospective clinical study in which 205 individuals with malignant melanoma had radical skin surgery at Odense University Hospital during the period 1962–1977 and were followed until the end of 1977. Survival time was censored by death from other causes and by the termination of the follow-up period. Of the 126 female patients, 28 were observed to die from the disease, and of the 79 male patients, 29 were observed to die from the disease. The data and more detailed description can be found in Andersen et al. (1992, p. 11). These authors found that sex was a signi cant risk factor, with men having higher risk than women; thus it makes sense to stratify by sex. Their analysis also suggested that the hazard rate for men is a constant or the survival time follows an exponential distribution. We applied our procedure to test this hypothesis for men. By using L = 5000 realizations to evaluate the distribution of the supremum test statistic under the null hypothesis, we obtained a p-value of 0.54. We also plotted the observed test process Zn (t) for men, along with 10 simulated realizations of the process Zn∗ (t) in Fig. 1. Clearly, Zn (t) falls within the range of the 10 simulated reference paths. Therefore, both the p-value and the plot showed no evidence against exponentiality.

408

G. Li, Y. Sun / Statistics & Probability Letters 47 (2000) 403 – 410

Fig. 1. The solid line is the realization of the test process calculated from the melaloma data for men; dashed lines are 10 realizations of the approximate limiting test process under the null hypothesis.

Table 1 The observed levels=powers (%) of the supremum test (4) at the 0.05 nominal level for testing exponentiality True distributions Censoring (%)

Size n

UE

WB (0.6)

50

50 60 75 100 125 150 200

5.9 7.2 5.7 7.4 6.2 5.2 6.3

52.8 63.8 73.8 84.2 90.4 95.4 97.7

20

50 60 75 100 125 150 200

4.3 5.3 5.9 4.9 6.4 5.2 5.6

71.6 83.6 91.6 97.9 99.7 100 100

WB (1.6)

GM (0.6)

GM (1.8)

LN (1.5)

22.3 27.0 31.5 43.4 58.1 73.7 92.9

31.3 36.7 44.2 61.1 67.1 75.3 87.9

17.2 18.7 23.6 29.4 37.9 48.1 69.7

21.6 29.6 36.0 46.8 56.6 62.0 76.7

90.0 93.1 98.4 99.7 99.9 100 100

32.1 44.2 54.7 74.7 84.2 91.0 97.7

67.6 77.2 83.1 94.7 96.1 99.4 99.8

49.9 61.5 73.0 89.6 95.0 96.8 99.8

To illustrate how our method would behave when H0 does not hold, we tested for exponentiality with a data that is known to be from a dierent distribution. The data are given in Table 1 of Hollander et al. (1997) concerning the number of days between a manuscript’s submission to the Theory and Methods Section of JASA and its rst review or the cut-o date, along with a censoring indicator (1 if a paper received its rst review by the cut-o date; 0 otherwise). Among the 432 submissions during 1994, 275 are uncensored and 157 are censored. Our test yielded a p-value of 0.004, which shows strong evidence that the exponential family

G. Li, Y. Sun / Statistics & Probability Letters 47 (2000) 403 – 410

409

Fig. 2. The solid line is the realization of the test process calculated from the time-to- rst-review data; dashed lines are 10 realizatins of the approximate limiting test process under the null hypothesis.

does not t well. We also plotted the observed test process Zn (t) along with 10 simulated realizations of Zn∗ (t); see Fig. 2. The plot shows that Zn (t) deviates severely from the 10 reference paths. In fact, the con dence bands of the survival function obtained by Hollander et al. (1997) suggested that the time-to- rst-review is uniformly distributed. A simulation study was conducted to investigate the level and rejection power of test (4). The null hypothesis was H0 : F ∈ F = {F(t; ) = 1 − exp(−t);  ¿ 0}. We generated failure times from the unit exponential (UE) with  = 1, the Weibull distribution (WB( )) with survival function S(t) = exp(−t ) for = 0:6 and 1.6, the gamma distribution (GM( )) with density function f(t) = t −1 exp(−t)= ( ) for = 0:6 and 1.8 and the log-normal distribution (LN( )) with survival function S(t) = 1 − ( −1 log t) for = 1:5, where  is the standard normal distribution function. The censoring times were generated from the exponential distribution with parameter adjusted to give prescribed censoring rate. The observed levels and powers of test (4) are given in Table 1. Each entry was computed using 1000 samples. For each sample, the asymptotic distribution of sup06t60 |Z˜ n (t)| was evaluated based on L=2000 realizations of Zn∗ (t). We took  to be the 80th percentile of the unit exponential distribution. It is seen from Table 1 that the observed level (corresponding to UE) of test (4) is close to the nominal level 0.05 in most cases. The rejection power is also quite satisfactory except for small samples (n660) with heavy (50%) censoring. References Aalen, O.O., 1978. Nonparametric inference for a family of counting processes. Ann. Statist. 6, 701–726. Aalen, O.O., Johansen, S., 1978. An empirical transition matrix for nonhomogeneous Markov chains based on censored observations. Scand. J. Statist. 5, 141–150. Akritas, M.G., 1988. Pearson-type goodness-of- t tests: the univariate case. J. Amer. Statist. Assoc. 83, 222–230. Andersen, P.K., Borgan, I., Gill, R.D., Keiding, N., 1992. Statistical Models Based on Counting Processes. Springer, New York. Billingsley, P., 1968. Convergence of Probability Measures. Wiley, New York.

410

G. Li, Y. Sun / Statistics & Probability Letters 47 (2000) 403 – 410

Doob, J., 1949. Heuristic approach to the Kolmogorov–Smirnov theorems. Ann. Math. Statist. 20, 393–403. Gill, R.D., 1980. Censoring and Stochastic Integrals. Math. Centre Tracts, Vol. 124. Mathematical Centre, Amsterdam. Gill, R.D., Johansen, S., 1990. A survey of product-integration with a view towards application in survival analysis. Ann. Statist. 18, 1501–1555. Habib, M.G., Thomas, D.R., 1986. Chi-square goodness-of- t tests for randomly censored data. Ann. Statist. 14, 759–765. Hjort, N.L., 1990. Goodness of t tests in Models for life history data based on cumulative hazard rates. Ann. Statist. 18, 1221–1258. Hollander, M., McKeague, I.W., Yang, J., 1997. Likelihood ratio-based con dence bands for survival functions. J. Amer. Statist. Assoc. 92, 215–226. Keiding, N., Gill, D., 1990. Random truncation models and Markov processes. Ann. Statist. 18, 582–602. Khmaladze, E.V., 1993. Goodness of t problem and scanning innovation martingales. Ann. Statist. 21, 798–829. Kim, J.H., 1993. Chi-square goodness-of- t tests for randomly censored data. Ann. Statist. 21, 1621–1639. Kolmogorov, A.N., 1933. Sulla determinazione empirica di una legge di distribuzione. Giorn. Instit. Ital. Attuari. 4, 83–91. Li, G., Doss, H., 1993. Generalized Pearson–Fisher chi-square goodness-of- t tests, with applications to models with life history data. Ann. Statist. 21, 772–797. Lin, D.Y., Wei, L.J., Ying, Z., 1993. Checking the Cox model with cumulative sums of martingale-based residuals. Biometrika 80, 557–572. Mihalko, D., Moore, D., 1980. Chi-square tests for Type II censored data. Ann. Statist. 8, 625–644. Moore, D.S., Spruill, M.C., 1975. Uni ed large-sample theory of general chi-square test of t. Ann. Statist. 3, 599–616. Nikabadze, A., Stute, W., 1997. Model checks under random censorship. Statist. Probab. Lett. 32, 249–259. Sun, Y., 1997. Weak convergence of the generalized parametric empirical processes and goodness-of- t tests for Parametric Models. Comm. Statist. A – Theory Methods 26 (10), 2393–2413. Van der Vaart, A.W., 1998. Asymptotic Statistics. Cambridge University Press, New York.