- Email: [email protected]
An optimal completion of the product limit estimator
ARTICLE IN PRESS
Statistics & Probability Letters 76 (2006) 913–919 www.elsevier.com/locate/stapro
An optimal completion of the product limit estimator Zhiqiang Chen,1, Eswar Phadia Department of Mathematics, The William Paterson University of New Jersey, Wayne, NJ 07470, USA Received 4 June 2003; accepted 7 September 2005 Available online 14 November 2005
Abstract It is well known that the product limit estimator is undefined beyond the largest observation if it is censored. Some completion methods are suggested in the literature (see e.g. [Efron, 1967. The two sample problem with censored data. Proceedings of the 5th Berkeley Symposium] and [Gill, 1980. Censoring and stochastic integrals. Mathematical Centre Tract No. 124, Mathematisch Centrum, Amsterdam]). In this note, we propose a completion method that is optimal in the sense that the expected value of the integrated squared error loss function is minimized. This method yields an estimator that falls between the above two extremes and possesses the same large sample properties. New bounds for the biases are also derived for the above-mentioned cases. r 2005 Elsevier B.V. All rights reserved. MSC: primary 62C15; 62D05 Keywords: Bias bound; Censored data; Kaplan–Meier estimator; Loss function; Proportional hazard model
1. Introduction The product limit (PL) estimator introduced by Kaplan and Meier (1958) has been used extensively in practice even though it is known that the estimator is undefined beyond the last observation (in the sense of ordered observations) if it happens to be a right censored observation. To overcome this shortfall, Efron (1967) suggested that it be defined as zero beyond the last observation, whether it is censored or not. He termed such an estimator as self-consistent. Gill (1980) on the other hand recommended that it be defined as the value of the estimator at the last observation. Some other completion methods have been suggested in the literature without adequate justification. In this note, we propose a completion method that is optimal in the sense that the expected value of the integrated squared error loss function is minimum. In our treatment, we use the exact formula of the kth moment derived by Phadia and Shao (1999). This approach yields an estimator that falls between the above two extremes but has the same large sample properties. New bounds for biases are also derived. Although our focus is to concentrate on the proportional hazard model, a general approach is adopted. Corresponding author. Tel.:+1 973 720 3382; fax: +1 973 720 3622. 1
E-mail address: [email protected] (Z. Chen). Partially supported by CFR, College of Science and Health, William Paterson University.
0167-7152/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.spl.2005.10.023
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
914
The organization of the paper is as follows. In Section 2 we introduce some notations and preliminaries. The main result is presented in Section 3. Bias bounds are then presented in Section 4. 2. Notations and preliminaries Let X i ; i ¼ 1; 2; . . . ; n, be i.i.d. sample from F , and Y i ; i ¼ 1; 2; . . . ; n i.i.d. (G), be the censoring variables, independent of X i , and we observe, as in the context of survival analysis, Z i ¼ minfX i ; Y i g and di ¼ IðX i pY i Þ; i ¼ 1; 2; . . . ; n. We assume that both F and G are continuous and defined on the interval ½0; 1Þ. If we further assume that di and Z i are independent, we have the proportional hazard model. Kaplan and Meier (1958) introduced the PL estimator for the survival function SðtÞ ¼ 1 F ðtÞ, defined as IðZi:n ptÞ di:n n1 b SðtÞ ¼ P1 1 , niþ1 where Z i:n denotes the ith ordered observation among Z 0i s and di:n corresponds to Z i:n . However, for t4Z n:n , and dn:n ¼ 0, the estimator was undefined. Efron (1967) introduced the notion of self-consistency and suggested that it be defined as 0, while Gill (1980) defined it as IðZi:n ptÞ di:n ð1 dn:n ÞP 1 ; for t4Zn:n . niþ1 cE ðtÞpSðtÞpE ScG ðtÞ, where S cE and ScG stands for the Efron’s and Gill’s versions, It is well known that E S respectively. Since the survival function is monotone, any sensible modification of the last term of the PL estimator will be IðZi:n ptÞ di:n 1 ; for t4Zn:n , Sbc ðtÞ ¼ cð1 dn:n ÞPn1 1 niþ1 where c will be a quantity between 0 and 1. Clearly, the extreme values of c yield Efron’s and Gill’s versions, respectively. The goal of this article is to suggest a completion which minimizes the mean squared error loss R1 EðSb SÞ2 ¼ EðFb F Þ2 . But instead of point-wise consideration, we will consider the loss function 0 ðFbðtÞ R1 R1 F ðtÞÞ2 dF ðtÞ and try to determine a constant c which minimizes E ðFbðtÞ F ðtÞÞ2 dF ðtÞ ¼ EðSb SÞ2 dS. 0
0
k
The results of Phadia and Shao (1999) on the exact formula for E Sb ðtÞ will be used. Rt Rt For this purpose, let H 1 ðtÞ ¼ PðZpt; X 4Y Þ ¼ 0 S dG, H 2 ðtÞ ¼ PðZpt; X pY Þ ¼ 0 G dF , and HðtÞ ¼ k PðZptÞ where G ¼ 1 G. For positive integer k, define fðkÞ j ðtÞ ¼ H 1 ðtÞ þ ððn jÞ=ðn j þ 1ÞÞ H 2 ðtÞ. For i ¼ 1; 2; . . . ; n, denote Z i:n ¼ Z i , such that Z 1 pZ 2 p pZ n are the ordered observations and di corresponds to Z i . Then Phadia and Shao (1999) showed that Z n k ni E Sb ðtÞ ¼ S H i! Iðzi otÞPij¼1 dðfðkÞ j ðzj ÞÞ, i
where the summation extends to n 1 for Efron’s version and to n for Gill’s version. They also used linear approximation for the integral and obtained Z i! Iðzi otÞPij¼1 dðfðkÞ e ij¼1 fðkÞ (2.1) j ðzj Þ¼P j ðtÞ k
yielding E Sb ðtÞ ¼ SðniÞH
ni
Pij¼1 fðkÞ j ðtÞ.
3. Optimal completion results We are now ready to derive the main result.
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
Theorem 1. The critical value of c such that the corresponding completion minimizes E R R1 ð1Þ n 0 SðtÞ½ Iðzn otÞPj¼1 dðfj ðzj ÞÞ dSðtÞ c ¼ R1 R , ð2Þ n 0 ½ Iðzn otÞPj¼1 dðfj ðzj ÞÞ dSðtÞ
915
R1 0
ðFbðtÞ F ðtÞÞ2 dF ðtÞ is (3.1)
which is approximately R1 ð1Þ n 0 SðtÞPj¼1 fj ðtÞ dSðtÞ . R 1 n ð2Þ 0 Pj¼1 fj ðtÞ dSðtÞ Proof. Differentiating the loss function with respect to c, we get Z 1 Z 1 Z 1 2 d d d 2 d dSðtÞ b fEðSb ðtÞÞg dSðtÞ. EðSðtÞ SðtÞÞ dSðtÞ ¼ 2SðtÞ fEðSðtÞÞg dc dc dc 0 0 0 R1 R1 R R1 b b2 Since 0 2SðtÞ d=dcfEðSðtÞÞg dSðtÞ ¼ 2 0 SðtÞn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ dSðtÞ and 0 d=dc fEðS ðtÞÞg dSðtÞ ¼ R1 R 2c 0 n! Iðzn otÞPnj¼1 dðfð2Þ j ðzj ÞÞ dSðtÞ, we get the critical value R R1 ð1Þ n 0 SðtÞ½ Iðzn otÞPj¼1 dðfj ðzj ÞÞ dSðtÞ c ¼ R1 R 40. ð2Þ n 0 ½ Iðzn otÞPj¼1 dðfj ðzj ÞÞ dSðtÞ The approximate value is obtained trivially by making the substitution given in (2.1). & R R1 cc F ðtÞÞ2 dF ðtÞ ¼ 1 EðSbc SÞ2 dS is a quadratic function of c, we see from the above Since E 0 ðF 0 proof that, if c 2 ð0; 1, then Sbc is the optimal completion. However, if c41, Gill’s version is the optimal because the constant c has to be in the interval of 0 and 1 in order for the completion to make sense. It is believed that the above critical value c is between 0 and 1, but it is rather difficult to prove it in general. However, it is true for the proportional hazard model as shown in Corollary 3. In any case, we define the completion to be optimal for values of c in the interval ð0; 1. Remark 1. Our optimal completion of the PL estimator and the two versions mentioned above differ only in the value beyond the last observation. In fact they are all very close and the differences tend to zero exponentially fast. This can be seen for instance as follows (see Phadia and Van Ryzin, 1980). With E i ¼ fZ i puoZ iþ1 g and Z nþ1 ¼ 1, we have cE ðtÞ Sbc ðtÞj2 jE n g PðE n Þ. cE ðtÞ Sbc ðtÞj2 g ¼ EfjS EfjS The first factor on the right-hand side is finite and PðE n Þ ¼ PðZ n puÞ ¼ ð1 HðuÞÞn ¼ exp½n lnðð1 HðuÞÞ1 Þ!0
if HðuÞ40.
Thus all of their asymptotic properties should be the same. Similar results hold for SbG . Under the proportional hazard model, Z and d are independent and therefore H 2 ðtÞ ¼ PðZpt; d ¼ 1Þ ¼ Pðd ¼ 1ÞPðZptÞ ¼ gHðtÞ, H 1 ðtÞ ¼ ð1 gÞHðtÞ, where g ¼ Pðd ¼ 1Þ ¼ PðX pY Þ. The next proposition gives a neat expression for c. Theorem 2. Under the proportional hazard model, the critical value of c in (3.1) reduces to 1 iði gÞði þ gÞ c ¼ Pn1 . 2 ði þ 2gÞ½ði gÞ2 þ gð1 gÞ
(3.2)
Proof. By substituting the values of H 1 ðtÞ and H 2 ðtÞ, it is easy to see that the constant c derived in Theorem 1 reduces to R1 SðtÞH n ðtÞ dSðtÞ n ðð1 gÞ þ ðn jÞ=ðn j þ 1ÞgÞ c ¼ 0R 1 n P1 . ðð1 gÞ þ ððn jÞ=ðn j þ 1ÞÞ2 gÞ 0 H ðtÞ dSðtÞ
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
916
R1 R1 To evaluate 0 H n ðtÞ dSðtÞ and 0 SðtÞH n ðtÞ dSðtÞ, note that SG ¼ 1 H, and dH ¼ 1=gG dS, where H ¼ 1 H. Now using the technique of integration by parts, we have Z 1 Z 1 n H ðtÞ dSðtÞ ¼ nSðtÞH n1 ðtÞ dHðtÞ 0 0 Z n 1 n1 ¼ H ðtÞSðtÞGðtÞ dSðtÞ g 0 Z 1 Z 1 n ¼ H n1 ðtÞ dSðtÞ H n ðtÞ dSðtÞ . g 0 0 Therefore by induction, Z 1 H n ðtÞ dSðtÞ ¼
Z 1 n H n1 ðtÞ dSðtÞ nþg 0 Z 1 n n1 ¼ H n2 ðtÞ dSðtÞ n þ g ðn 1Þ þ g 0 Z 1 n! ¼ n dSðtÞ P1 ði þ gÞ 0 n! ¼ n . P1 ði þ gÞ
0
Similarly, Z 1
Z 1 1 2 S ðtÞnH n1 ðtÞ dHðtÞ 2 0 Z n 1 ¼ SðtÞH n1 ðtÞSðtÞGðtÞ dSðtÞ 2g 0 Z n 1 ¼ SðtÞ½H n1 ðtÞ H n ðtÞ dSðtÞ, 2g 0
SðtÞH n ðtÞ dSðtÞ ¼
0
and by induction, we obtain Z 1 SðtÞH n ðtÞ dSðtÞ ¼
Z 1 n SH n1 dS n þ 2g 0 Z 1 n! ¼ n S dS P1 ði þ 2gÞ 0 n! ¼ . 2Pn1 ði þ 2gÞ
0
Hence R1 n 1 ði þ gÞ 0R SðtÞH ðtÞ dSðtÞ . ¼ Pn1 1 n 2 ði þ 2gÞ H ðtÞ dSðtÞ 0 Now we will simplify Pn1
ðð1 gÞ þ ðn jÞ=ðn j þ 1ÞgÞ , ðð1 gÞ þ ððn jÞ=ðn j þ 1ÞÞ2 gÞ
Pn1
2 ðð1 gÞ þ ðn jÞ=ðn j þ 1ÞgÞ n ½ðn j þ 1Þ ð1 gÞ þ ðn j þ 1Þðn jÞg ¼ P 1 ðð1 gÞ þ ððn jÞ=ðn j þ 1ÞÞ2 gÞ ½ðn j þ 1Þ2 ð1 gÞ þ ðn jÞ2 g
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
917
½i2 ð1 gÞ þ iði 1Þg ½i2 ð1 gÞ þ ði 1Þ2 g iði gÞ , ¼ Pn1 ½ði gÞ2 þ gð1 gÞ
¼ Pn1
where in the second equality above, we substituted n j þ 1 ¼ i. Putting together these terms we get (3.2).
&
Thus, the optimality constant c depends on the proportional hazard rate g ¼ Pðd ¼ 1Þ. Since d is observable, g and hence, c and Sbc can easily be estimated from the given data. Our next corollary shows that c is within the (0; 1) range as we expected, and therefore, Sbc is optimal. Corollary 3. Under the proportional hazard model, the optimality constant c, a function of g, satisfies 1 4 g2 4 g2 oco , 4 4 g2 þ 5gð1 gÞ 4 g2 þ 5gð1 gÞ and hence c 2 ð0; 1Þ. Proof. Since cðgÞ ¼
ð1 þ gÞ 2ð2 gÞð2 þ gÞ iði gÞði þ gÞ Pn3 2ð1 þ 2gÞ ð2 þ 2gÞð4 3gÞ ði þ 2gÞ½ði gÞ2 þ gð1 gÞ
¼
4 g2 1 iði gÞði þ gÞ Pn 4 g2 þ 5gð1 gÞ 2 3 ði þ 2gÞ½ði gÞ2 þ gð1 gÞ
¼
4 g2 f ðgÞ; 4 g2 þ 5gð1 gÞ
say.
Denote iði 2gÞ D 1 LðgÞ ¼ Pn3 2 ði gÞ2 þ gð1 gÞ
and
iði þ gÞ D 1 . UðgÞ ¼ Pn3 2 ði þ 2gÞði gÞ
Then we have 1 iði gÞði þ gÞ LðgÞpf ðgÞ ¼ Pn3 pUðgÞ. 2 ði þ 2gÞ½ði gÞ2 þ gð1 gÞ Simple calculation shows that 1 2 1 2gð2i þ gÞ 0 n þ 40, ¼ Sn3 ½ln 2UðgÞ ¼ S3 ði þ gÞ ði þ 2gÞ ði gÞ ði gÞði þ gÞði þ 2gÞ so UðgÞ is an increasing function of g when g 2 ð0; 1Þ, therefore 1 iði þ 1Þ 1 n 4 n ¼ ¼ o1. f ðgÞpUðgÞpUð1Þ ¼ Pn3 2 ði þ 2Þði 1Þ 2 2 n þ 2 n þ 2 Similarly, ½ln 2LðgÞ0 o0, so f ðgÞXLðgÞXLð1Þ ¼ n=4ðn 1Þ414. & 4. Bias bounds cE ðtÞ and ScG ðtÞ be as defined earlier, the Efron’s and Gill’s completion, respectively, and Sbc ðtÞ be the Let S b SðtÞ be the bias with suffix indicating optimal completion considered in this article. Further, let BðtÞ ¼ E SðtÞ cE ðtÞpSðtÞ, and Klein (1991) proved that the bias of a particular completion. Efron (1967) showed that E S Rt n SðtÞpE ScG ðtÞ. A special case of Zhou (1988) shows that jBG jp 0 H ðtÞ dF ðtÞ. We now give a different bias bound.
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
918
Proposition 4. For both Efron’s and Gill’s completion, the bound on bias B is Z 1 ð1Þ n n jBjon! Iðzn otÞPj¼1 dðfj ðzj ÞÞ Pj¼1 HðtÞ H 2 ðtÞ , j and for the completion discussed in this article, the bound is Z 1 ð1Þ n n jBc jp maxfc; 1 cgn! Iðzn otÞPj¼1 dðfj ðzj ÞÞ maxfc; 1 cgPj¼1 HðtÞ H 2 ðtÞ . j Proof. We have Z c 0oBG ¼ E S E ðtÞ SðtÞ þ n! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ Z on! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ, and cE ðtÞ SðtÞ4 n! 04BE ¼ E S
Z
Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ,
Therefore, R jBjon! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ for both Efron’s and Gill’s versions. Now Z c Bc ¼ E S G ðtÞ SðtÞ þ ðc 1Þn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ Z X ð1 cÞn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ, and Z c Bc ¼ E S E ðtÞ SðtÞ þ cn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj Þ Z pcn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ. R So, jBc jp maxfc; 1 cgn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ.
&
In the case of proportional hazard model, Z 1 n n g , n! Iðzn otÞPnj¼1 dðfð1Þ ðz Þ ¼ H ðtÞP 1 j j 1 i (see Chen et al., 1982), therefore we have the following. Corollary 5. In the case of proportional hazard model, the bounds on biases reduces to 1 1 jBjoH n ðtÞPn1 1 g and jBc jp maxfc; 1 cgH n ðtÞPn1 1 g . i i The above bias bounds are easy to obtain. However, they are new to the best of our knowledge. To see how good these bounds are, we compareRthem with Zhou’s (1988) result. In the case of propositional hazard model, t earlier computation shows that 0 H n ðxÞ dF ðxÞ ! n!=Pn1 ði þ gÞ as t ! 1. Since H n ðtÞPn1 ð1 ð1=iÞgÞ ! ð1=n!ÞPn1 ði gÞ, as t ! 1, and clearly ð1=n!ÞPn1 ði gÞon!=Pn1 ði þ gÞ when g 2 ð0; 1Þ, the new bias bound given above The question whether the new bound is always better (that is, whether R t n is sharper Rfor t ! 1. ð1Þ n 0 H ðtÞ dF ðtÞXn! Iðzn otÞPj¼1 dðfj ðzj ÞÞ) for all t40, is still open.
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
919
References Chen, Y.Y., Hollander, M., Langberg, N.A., 1982. Small-sample results for the Kaplan–Meier estimator. J. Amer. Statist. Assoc. 77, 141–144. Efron, B., 1967. The two sample problem with censored data. Proceedings of the 5th Berkeley Symposium vol. 4, pp. 831–852. Gill, R.D., 1980. Censoring and stochastic integrals. Mathematical Centre Tract No. 124. Mathematisch Centrum, Amsterdam. Kaplan, E.L., Meier, P., 1958. Nonparametric estimation from incomplete observations. J. Amer. Statist. Assoc. 53, 457–481. Klein, J.P., 1991. Small sample moments of some estimators of the variance of the Kaplan–Meier and Nelson–Aalen estimators. Scand. J. Statist. 18, 333–340. Phadia, E., Shao, Y., 1999. Exact moments of the product limit estimator. Statist. Probab. Lett. 41, 277–286. Phadia, E.G., Van Ryzin, J., 1980. A note on convergence rates for the product limit estimator. Ann. Statist. 8, 673–678. Zhou, M., 1988. Two sided bias bounds of the Kaplan–Meier estimator. Probab. Theory Related Fields 79, 165–173.
Statistics & Probability Letters 76 (2006) 913–919 www.elsevier.com/locate/stapro
An optimal completion of the product limit estimator Zhiqiang Chen,1, Eswar Phadia Department of Mathematics, The William Paterson University of New Jersey, Wayne, NJ 07470, USA Received 4 June 2003; accepted 7 September 2005 Available online 14 November 2005
Abstract It is well known that the product limit estimator is undefined beyond the largest observation if it is censored. Some completion methods are suggested in the literature (see e.g. [Efron, 1967. The two sample problem with censored data. Proceedings of the 5th Berkeley Symposium] and [Gill, 1980. Censoring and stochastic integrals. Mathematical Centre Tract No. 124, Mathematisch Centrum, Amsterdam]). In this note, we propose a completion method that is optimal in the sense that the expected value of the integrated squared error loss function is minimized. This method yields an estimator that falls between the above two extremes and possesses the same large sample properties. New bounds for the biases are also derived for the above-mentioned cases. r 2005 Elsevier B.V. All rights reserved. MSC: primary 62C15; 62D05 Keywords: Bias bound; Censored data; Kaplan–Meier estimator; Loss function; Proportional hazard model
1. Introduction The product limit (PL) estimator introduced by Kaplan and Meier (1958) has been used extensively in practice even though it is known that the estimator is undefined beyond the last observation (in the sense of ordered observations) if it happens to be a right censored observation. To overcome this shortfall, Efron (1967) suggested that it be defined as zero beyond the last observation, whether it is censored or not. He termed such an estimator as self-consistent. Gill (1980) on the other hand recommended that it be defined as the value of the estimator at the last observation. Some other completion methods have been suggested in the literature without adequate justification. In this note, we propose a completion method that is optimal in the sense that the expected value of the integrated squared error loss function is minimum. In our treatment, we use the exact formula of the kth moment derived by Phadia and Shao (1999). This approach yields an estimator that falls between the above two extremes but has the same large sample properties. New bounds for biases are also derived. Although our focus is to concentrate on the proportional hazard model, a general approach is adopted. Corresponding author. Tel.:+1 973 720 3382; fax: +1 973 720 3622. 1
E-mail address: [email protected] (Z. Chen). Partially supported by CFR, College of Science and Health, William Paterson University.
0167-7152/$ - see front matter r 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.spl.2005.10.023
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
914
The organization of the paper is as follows. In Section 2 we introduce some notations and preliminaries. The main result is presented in Section 3. Bias bounds are then presented in Section 4. 2. Notations and preliminaries Let X i ; i ¼ 1; 2; . . . ; n, be i.i.d. sample from F , and Y i ; i ¼ 1; 2; . . . ; n i.i.d. (G), be the censoring variables, independent of X i , and we observe, as in the context of survival analysis, Z i ¼ minfX i ; Y i g and di ¼ IðX i pY i Þ; i ¼ 1; 2; . . . ; n. We assume that both F and G are continuous and defined on the interval ½0; 1Þ. If we further assume that di and Z i are independent, we have the proportional hazard model. Kaplan and Meier (1958) introduced the PL estimator for the survival function SðtÞ ¼ 1 F ðtÞ, defined as IðZi:n ptÞ di:n n1 b SðtÞ ¼ P1 1 , niþ1 where Z i:n denotes the ith ordered observation among Z 0i s and di:n corresponds to Z i:n . However, for t4Z n:n , and dn:n ¼ 0, the estimator was undefined. Efron (1967) introduced the notion of self-consistency and suggested that it be defined as 0, while Gill (1980) defined it as IðZi:n ptÞ di:n ð1 dn:n ÞP 1 ; for t4Zn:n . niþ1 cE ðtÞpSðtÞpE ScG ðtÞ, where S cE and ScG stands for the Efron’s and Gill’s versions, It is well known that E S respectively. Since the survival function is monotone, any sensible modification of the last term of the PL estimator will be IðZi:n ptÞ di:n 1 ; for t4Zn:n , Sbc ðtÞ ¼ cð1 dn:n ÞPn1 1 niþ1 where c will be a quantity between 0 and 1. Clearly, the extreme values of c yield Efron’s and Gill’s versions, respectively. The goal of this article is to suggest a completion which minimizes the mean squared error loss R1 EðSb SÞ2 ¼ EðFb F Þ2 . But instead of point-wise consideration, we will consider the loss function 0 ðFbðtÞ R1 R1 F ðtÞÞ2 dF ðtÞ and try to determine a constant c which minimizes E ðFbðtÞ F ðtÞÞ2 dF ðtÞ ¼ EðSb SÞ2 dS. 0
0
k
The results of Phadia and Shao (1999) on the exact formula for E Sb ðtÞ will be used. Rt Rt For this purpose, let H 1 ðtÞ ¼ PðZpt; X 4Y Þ ¼ 0 S dG, H 2 ðtÞ ¼ PðZpt; X pY Þ ¼ 0 G dF , and HðtÞ ¼ k PðZptÞ where G ¼ 1 G. For positive integer k, define fðkÞ j ðtÞ ¼ H 1 ðtÞ þ ððn jÞ=ðn j þ 1ÞÞ H 2 ðtÞ. For i ¼ 1; 2; . . . ; n, denote Z i:n ¼ Z i , such that Z 1 pZ 2 p pZ n are the ordered observations and di corresponds to Z i . Then Phadia and Shao (1999) showed that Z n k ni E Sb ðtÞ ¼ S H i! Iðzi otÞPij¼1 dðfðkÞ j ðzj ÞÞ, i
where the summation extends to n 1 for Efron’s version and to n for Gill’s version. They also used linear approximation for the integral and obtained Z i! Iðzi otÞPij¼1 dðfðkÞ e ij¼1 fðkÞ (2.1) j ðzj Þ¼P j ðtÞ k
yielding E Sb ðtÞ ¼ SðniÞH
ni
Pij¼1 fðkÞ j ðtÞ.
3. Optimal completion results We are now ready to derive the main result.
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
Theorem 1. The critical value of c such that the corresponding completion minimizes E R R1 ð1Þ n 0 SðtÞ½ Iðzn otÞPj¼1 dðfj ðzj ÞÞ dSðtÞ c ¼ R1 R , ð2Þ n 0 ½ Iðzn otÞPj¼1 dðfj ðzj ÞÞ dSðtÞ
915
R1 0
ðFbðtÞ F ðtÞÞ2 dF ðtÞ is (3.1)
which is approximately R1 ð1Þ n 0 SðtÞPj¼1 fj ðtÞ dSðtÞ . R 1 n ð2Þ 0 Pj¼1 fj ðtÞ dSðtÞ Proof. Differentiating the loss function with respect to c, we get Z 1 Z 1 Z 1 2 d d d 2 d dSðtÞ b fEðSb ðtÞÞg dSðtÞ. EðSðtÞ SðtÞÞ dSðtÞ ¼ 2SðtÞ fEðSðtÞÞg dc dc dc 0 0 0 R1 R1 R R1 b b2 Since 0 2SðtÞ d=dcfEðSðtÞÞg dSðtÞ ¼ 2 0 SðtÞn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ dSðtÞ and 0 d=dc fEðS ðtÞÞg dSðtÞ ¼ R1 R 2c 0 n! Iðzn otÞPnj¼1 dðfð2Þ j ðzj ÞÞ dSðtÞ, we get the critical value R R1 ð1Þ n 0 SðtÞ½ Iðzn otÞPj¼1 dðfj ðzj ÞÞ dSðtÞ c ¼ R1 R 40. ð2Þ n 0 ½ Iðzn otÞPj¼1 dðfj ðzj ÞÞ dSðtÞ The approximate value is obtained trivially by making the substitution given in (2.1). & R R1 cc F ðtÞÞ2 dF ðtÞ ¼ 1 EðSbc SÞ2 dS is a quadratic function of c, we see from the above Since E 0 ðF 0 proof that, if c 2 ð0; 1, then Sbc is the optimal completion. However, if c41, Gill’s version is the optimal because the constant c has to be in the interval of 0 and 1 in order for the completion to make sense. It is believed that the above critical value c is between 0 and 1, but it is rather difficult to prove it in general. However, it is true for the proportional hazard model as shown in Corollary 3. In any case, we define the completion to be optimal for values of c in the interval ð0; 1. Remark 1. Our optimal completion of the PL estimator and the two versions mentioned above differ only in the value beyond the last observation. In fact they are all very close and the differences tend to zero exponentially fast. This can be seen for instance as follows (see Phadia and Van Ryzin, 1980). With E i ¼ fZ i puoZ iþ1 g and Z nþ1 ¼ 1, we have cE ðtÞ Sbc ðtÞj2 jE n g PðE n Þ. cE ðtÞ Sbc ðtÞj2 g ¼ EfjS EfjS The first factor on the right-hand side is finite and PðE n Þ ¼ PðZ n puÞ ¼ ð1 HðuÞÞn ¼ exp½n lnðð1 HðuÞÞ1 Þ!0
if HðuÞ40.
Thus all of their asymptotic properties should be the same. Similar results hold for SbG . Under the proportional hazard model, Z and d are independent and therefore H 2 ðtÞ ¼ PðZpt; d ¼ 1Þ ¼ Pðd ¼ 1ÞPðZptÞ ¼ gHðtÞ, H 1 ðtÞ ¼ ð1 gÞHðtÞ, where g ¼ Pðd ¼ 1Þ ¼ PðX pY Þ. The next proposition gives a neat expression for c. Theorem 2. Under the proportional hazard model, the critical value of c in (3.1) reduces to 1 iði gÞði þ gÞ c ¼ Pn1 . 2 ði þ 2gÞ½ði gÞ2 þ gð1 gÞ
(3.2)
Proof. By substituting the values of H 1 ðtÞ and H 2 ðtÞ, it is easy to see that the constant c derived in Theorem 1 reduces to R1 SðtÞH n ðtÞ dSðtÞ n ðð1 gÞ þ ðn jÞ=ðn j þ 1ÞgÞ c ¼ 0R 1 n P1 . ðð1 gÞ þ ððn jÞ=ðn j þ 1ÞÞ2 gÞ 0 H ðtÞ dSðtÞ
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
916
R1 R1 To evaluate 0 H n ðtÞ dSðtÞ and 0 SðtÞH n ðtÞ dSðtÞ, note that SG ¼ 1 H, and dH ¼ 1=gG dS, where H ¼ 1 H. Now using the technique of integration by parts, we have Z 1 Z 1 n H ðtÞ dSðtÞ ¼ nSðtÞH n1 ðtÞ dHðtÞ 0 0 Z n 1 n1 ¼ H ðtÞSðtÞGðtÞ dSðtÞ g 0 Z 1 Z 1 n ¼ H n1 ðtÞ dSðtÞ H n ðtÞ dSðtÞ . g 0 0 Therefore by induction, Z 1 H n ðtÞ dSðtÞ ¼
Z 1 n H n1 ðtÞ dSðtÞ nþg 0 Z 1 n n1 ¼ H n2 ðtÞ dSðtÞ n þ g ðn 1Þ þ g 0 Z 1 n! ¼ n dSðtÞ P1 ði þ gÞ 0 n! ¼ n . P1 ði þ gÞ
0
Similarly, Z 1
Z 1 1 2 S ðtÞnH n1 ðtÞ dHðtÞ 2 0 Z n 1 ¼ SðtÞH n1 ðtÞSðtÞGðtÞ dSðtÞ 2g 0 Z n 1 ¼ SðtÞ½H n1 ðtÞ H n ðtÞ dSðtÞ, 2g 0
SðtÞH n ðtÞ dSðtÞ ¼
0
and by induction, we obtain Z 1 SðtÞH n ðtÞ dSðtÞ ¼
Z 1 n SH n1 dS n þ 2g 0 Z 1 n! ¼ n S dS P1 ði þ 2gÞ 0 n! ¼ . 2Pn1 ði þ 2gÞ
0
Hence R1 n 1 ði þ gÞ 0R SðtÞH ðtÞ dSðtÞ . ¼ Pn1 1 n 2 ði þ 2gÞ H ðtÞ dSðtÞ 0 Now we will simplify Pn1
ðð1 gÞ þ ðn jÞ=ðn j þ 1ÞgÞ , ðð1 gÞ þ ððn jÞ=ðn j þ 1ÞÞ2 gÞ
Pn1
2 ðð1 gÞ þ ðn jÞ=ðn j þ 1ÞgÞ n ½ðn j þ 1Þ ð1 gÞ þ ðn j þ 1Þðn jÞg ¼ P 1 ðð1 gÞ þ ððn jÞ=ðn j þ 1ÞÞ2 gÞ ½ðn j þ 1Þ2 ð1 gÞ þ ðn jÞ2 g
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
917
½i2 ð1 gÞ þ iði 1Þg ½i2 ð1 gÞ þ ði 1Þ2 g iði gÞ , ¼ Pn1 ½ði gÞ2 þ gð1 gÞ
¼ Pn1
where in the second equality above, we substituted n j þ 1 ¼ i. Putting together these terms we get (3.2).
&
Thus, the optimality constant c depends on the proportional hazard rate g ¼ Pðd ¼ 1Þ. Since d is observable, g and hence, c and Sbc can easily be estimated from the given data. Our next corollary shows that c is within the (0; 1) range as we expected, and therefore, Sbc is optimal. Corollary 3. Under the proportional hazard model, the optimality constant c, a function of g, satisfies 1 4 g2 4 g2 oco , 4 4 g2 þ 5gð1 gÞ 4 g2 þ 5gð1 gÞ and hence c 2 ð0; 1Þ. Proof. Since cðgÞ ¼
ð1 þ gÞ 2ð2 gÞð2 þ gÞ iði gÞði þ gÞ Pn3 2ð1 þ 2gÞ ð2 þ 2gÞð4 3gÞ ði þ 2gÞ½ði gÞ2 þ gð1 gÞ
¼
4 g2 1 iði gÞði þ gÞ Pn 4 g2 þ 5gð1 gÞ 2 3 ði þ 2gÞ½ði gÞ2 þ gð1 gÞ
¼
4 g2 f ðgÞ; 4 g2 þ 5gð1 gÞ
say.
Denote iði 2gÞ D 1 LðgÞ ¼ Pn3 2 ði gÞ2 þ gð1 gÞ
and
iði þ gÞ D 1 . UðgÞ ¼ Pn3 2 ði þ 2gÞði gÞ
Then we have 1 iði gÞði þ gÞ LðgÞpf ðgÞ ¼ Pn3 pUðgÞ. 2 ði þ 2gÞ½ði gÞ2 þ gð1 gÞ Simple calculation shows that 1 2 1 2gð2i þ gÞ 0 n þ 40, ¼ Sn3 ½ln 2UðgÞ ¼ S3 ði þ gÞ ði þ 2gÞ ði gÞ ði gÞði þ gÞði þ 2gÞ so UðgÞ is an increasing function of g when g 2 ð0; 1Þ, therefore 1 iði þ 1Þ 1 n 4 n ¼ ¼ o1. f ðgÞpUðgÞpUð1Þ ¼ Pn3 2 ði þ 2Þði 1Þ 2 2 n þ 2 n þ 2 Similarly, ½ln 2LðgÞ0 o0, so f ðgÞXLðgÞXLð1Þ ¼ n=4ðn 1Þ414. & 4. Bias bounds cE ðtÞ and ScG ðtÞ be as defined earlier, the Efron’s and Gill’s completion, respectively, and Sbc ðtÞ be the Let S b SðtÞ be the bias with suffix indicating optimal completion considered in this article. Further, let BðtÞ ¼ E SðtÞ cE ðtÞpSðtÞ, and Klein (1991) proved that the bias of a particular completion. Efron (1967) showed that E S Rt n SðtÞpE ScG ðtÞ. A special case of Zhou (1988) shows that jBG jp 0 H ðtÞ dF ðtÞ. We now give a different bias bound.
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
918
Proposition 4. For both Efron’s and Gill’s completion, the bound on bias B is Z 1 ð1Þ n n jBjon! Iðzn otÞPj¼1 dðfj ðzj ÞÞ Pj¼1 HðtÞ H 2 ðtÞ , j and for the completion discussed in this article, the bound is Z 1 ð1Þ n n jBc jp maxfc; 1 cgn! Iðzn otÞPj¼1 dðfj ðzj ÞÞ maxfc; 1 cgPj¼1 HðtÞ H 2 ðtÞ . j Proof. We have Z c 0oBG ¼ E S E ðtÞ SðtÞ þ n! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ Z on! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ, and cE ðtÞ SðtÞ4 n! 04BE ¼ E S
Z
Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ,
Therefore, R jBjon! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ for both Efron’s and Gill’s versions. Now Z c Bc ¼ E S G ðtÞ SðtÞ þ ðc 1Þn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ Z X ð1 cÞn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ, and Z c Bc ¼ E S E ðtÞ SðtÞ þ cn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj Þ Z pcn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ. R So, jBc jp maxfc; 1 cgn! Iðzn otÞPnj¼1 dðfð1Þ j ðzj ÞÞ.
&
In the case of proportional hazard model, Z 1 n n g , n! Iðzn otÞPnj¼1 dðfð1Þ ðz Þ ¼ H ðtÞP 1 j j 1 i (see Chen et al., 1982), therefore we have the following. Corollary 5. In the case of proportional hazard model, the bounds on biases reduces to 1 1 jBjoH n ðtÞPn1 1 g and jBc jp maxfc; 1 cgH n ðtÞPn1 1 g . i i The above bias bounds are easy to obtain. However, they are new to the best of our knowledge. To see how good these bounds are, we compareRthem with Zhou’s (1988) result. In the case of propositional hazard model, t earlier computation shows that 0 H n ðxÞ dF ðxÞ ! n!=Pn1 ði þ gÞ as t ! 1. Since H n ðtÞPn1 ð1 ð1=iÞgÞ ! ð1=n!ÞPn1 ði gÞ, as t ! 1, and clearly ð1=n!ÞPn1 ði gÞon!=Pn1 ði þ gÞ when g 2 ð0; 1Þ, the new bias bound given above The question whether the new bound is always better (that is, whether R t n is sharper Rfor t ! 1. ð1Þ n 0 H ðtÞ dF ðtÞXn! Iðzn otÞPj¼1 dðfj ðzj ÞÞ) for all t40, is still open.
ARTICLE IN PRESS Z. Chen, E. Phadia / Statistics & Probability Letters 76 (2006) 913–919
919
References Chen, Y.Y., Hollander, M., Langberg, N.A., 1982. Small-sample results for the Kaplan–Meier estimator. J. Amer. Statist. Assoc. 77, 141–144. Efron, B., 1967. The two sample problem with censored data. Proceedings of the 5th Berkeley Symposium vol. 4, pp. 831–852. Gill, R.D., 1980. Censoring and stochastic integrals. Mathematical Centre Tract No. 124. Mathematisch Centrum, Amsterdam. Kaplan, E.L., Meier, P., 1958. Nonparametric estimation from incomplete observations. J. Amer. Statist. Assoc. 53, 457–481. Klein, J.P., 1991. Small sample moments of some estimators of the variance of the Kaplan–Meier and Nelson–Aalen estimators. Scand. J. Statist. 18, 333–340. Phadia, E., Shao, Y., 1999. Exact moments of the product limit estimator. Statist. Probab. Lett. 41, 277–286. Phadia, E.G., Van Ryzin, J., 1980. A note on convergence rates for the product limit estimator. Ann. Statist. 8, 673–678. Zhou, M., 1988. Two sided bias bounds of the Kaplan–Meier estimator. Probab. Theory Related Fields 79, 165–173.