Lifetime Data Analysis manuscript No. (will be inserted by the editor)

A cohort-weighted Kaplan-Meier statistic for addressing random non-homogeneity in survival comparisons

January 11, 2017

Keywords Survival analysis · Kaplan-Meier · Heterogeneous distribution · Nonparametric · Hypothesis test · Asymptotic analysis

1 Introduction Survival analysis addresses the classical statistical problem of determining characteristics of the waiting time unJ.C. Chang Epidemiology and Biostatistics Section Rehabilitation Medicine Department The National Institutes of Health, Clinical Center Bethesda, Maryland 20892 U.S.A. E-mail: [email protected] A. Heuser and M. Huynh Impaq International LLC Columbia, MD 21044 U.S.A. E-mail: [email protected] E-mail: [email protected]

Population 1 Population 2 1

1

rt ho

survival probability

Abstract The Kaplan-Meier product-limit estimator is a simple and powerful tool in time to event analysis. However, by design, it is agnostic to the influence of covariates. Hence it is not suited for resolving issues of heterogeneity and differential censoring that may feature in real applications, except through extensions or modifications. A specific example of these issues occurs in longitudinal studies comparing populations where the underlying survival functions are non-stationary. Motivated by this problem, we introduce a weighted product limit estimator for the population survival function under random representation of constituent cohorts. Based on this estimator we provide a test statistic for comparing populations. We derive the asymptotic behavior of the statistic based on an empirical process of productlimit estimators.

survival probability

arXiv:1701.02424v1 [stat.ME] 10 Jan 2017

Aaron Heuser · Minh Huynh · Joshua C. Chang

(z )

try en

co

(a)

surv ival time

2

surv

1

t

rt ho co

(b)

ival time

1

2

3

4

t

Fig. 1 Inhomogeneity of survival within populations can result due to at least two reasons. In (a), inhomogeneity results from a categorical covariate that influences survival statistics. In (b), inhomogeneity results from non-stationarity, where cohorts of individuals are sampled at different times. In this case, the problem of progressive censoring is apparent because later cohorts have not been observed as long.

til an event, canonically death, from observations of their occurrence sampled from within a population. This problem is not trivial as the expected waiting time is typically dependent on the time-already-waited. For instance, a hundred-year-old can be more certain of surviving to his or her one hundred and-first birthday than a newborn might reasonably be. However, the comparison may shift in the newborn’s favor for the living to one-hundred and twenty-one, particularly in light of medical advances that make survival probabilities non-stationary. Parametric approaches for assembling survival curves are usually not flexible enough to capture this complexity. One simple approach to this problem was pioneered by the work of Kaplan and Meier [8]. Their productlimit estimator [5, 6, 16, 17] is a non-parametric statistic that is used for inferring the survival function for members of a population from observed lifetimes. This method is particularly useful in that it naturally handles the presence of right censoring, where some event-times

2

are only partially observed because they fall outside the observation window. It was not, however, designed to account for varying subpopulations that may yield nonhomogeneity in overall population survival (Fig. 1). For instance, in the example given above, subpopulations for survival characteristics may be defined by birth year or entry cohort of a subject in a particular study (Fig. 1). Several existing statistical methods address variants of this limitation. A natural approach is to consider the varying subpopulations as defining underlying covariates, thus laying the framework for a proportional hazards model. The assumption of proportional hazards is quite strong. When considering time-dependent statistics (as in the motivational example), it is violated in all but a few specific cases. Likewise, frailty models, first developed by Hougaard (cf. [7]), and extended by Aalen (cf. [1]), assume multivariate event distributions, but also make assumptions on the underlying event distributions and assume proportional hazards. Other existing methods, such as bivariate survival analysis (cf. [11]), consider the time to observation and the time to event as conditionally independent random times. Underlying these methods is the assumption that upon the time of observation, all individuals will then have a similar event time distribution, thus failing to acknowledge the temporal changes. Lastly, in the work of Pepe and Fleming (cf. [12, 13]), a class of weighted Kaplan-Meier statistics is introduced. Though these statistics exhibit the same limitations as in the standard Kaplan-Meier case, it should be noted that [13] introduces the stratified weighted Kaplan-Meier statistic. The statistic presented here is a priori quite similar, but instead of a weighting function, includes the empirical prevalence. In doing so, the weight is no longer independent of the event time estimate, and thus requires much different methods of proof. We thus consider the overall survival distribution for a population of individuals with sub-populations that exhibit non-homogeneous survival distributions. Through this consideration, a new test statistic, based upon the empirical process of product-limit estimators is developed. Through constructive methods, this test-statistic compares survival distributions among the distinct subpopulations, and weights according to distribution of the identified subgroups.

Aaron Heuser et al.

assumption that survival is conditional on cohort z γ and population. One representation of the marginal survival probabil-  (i) ity for members of population i, θt = P T γ > t | γ ∈ Γ (i) , is found by conditioning on cohort, (i)

θt =

d o n X P T γ > t | z γ = z, γ ∈ Γ (i) {z } z=1 | (i)

Sz,t

n o × P z γ = z | γ ∈ Γ (i) , | {z }

(2.1)

(i)

qz

(i)

where Sz,t represents the survival function for individuals of cohort z in population i. We use this representation of the survival probability as motivation to formulate an estimator for the population-average survival functions (i) θˆt =

d X

(i)

qˆz(i) Sˆz,t ,

(2.2)

z=1 (i)

(i)

where qˆz and Sˆz,t are estimators of the cohort prevalence and cohort-wise survival respectively. Our concern is the general situation where random samples of size n(i) are chosen from each of the respective populations. Within these samples, the number of indi(i) viduals within each cohort, nz , is counted, from which an estimator of the cohort distribution is obtained, (i)

qˆz(i) =

nz . n(i)

(2.3)

2 Statistical method

In turn, we assume that the cohort-level survival func(i) tions Sˆz,t are estimated independently using the productlimit estimator. Note that since the product limit esti(i) mator is not a linear functional of sampled lifetimes, θˆt is distinct from the estimator obtained by applying the product limit estimator directly on all n(i) samples of population i. To prevent confusion, we denote all direct applications of the product-limit estimator using Sˆ and all instances of weighted sums of product limit estimaˆ tors using the Greek letter θ. With these elements in place, we define our test statistic s Z τ   n(1) n(2) ˆ(1) − θˆ(2) , ˆ Θ= dt θ (2.4) t t n(1) + n(2) 0

Suppose Γ (1) and Γ (2) are disjoint populations of individuals where each individual belongs to exactly one of d distinct cohorts labeled z ∈ Zd . For randomly selected individuals γ ∈ Γ (i) within population i, we desire to understand the statistics of the event time T γ under the

where τ = inf {τz : z ∈ Zd }, and τz denotes the time at which cohort z is censored in observations. We state here the main result of the paper – the asymptotic behavior of this statistic within a null-hypothesis significance testing framework.

A Survival Estimator and Statistic for Nonparametric Heterogeneous Survival Analysis (i)

Theorem 1 Let Cz,t denote the probability that a z-type individual has not yet been censored at time t ≥ 0, and (i) qz denote the probability that an individual in population i is of cohort z, and let p(i) = n(i) /(n(1) +n(2) ). Supd (1) (2) ˆ− pose that θt = θt . Then Θ → N (0, σ 2 ), as n(i) → ∞, with

σ2 =

2 X i=1

+

 d X (1 − p(i) )  qz(i) φ2z − z=1

d Z τz X z=1

 dSz,t Wz,t

0

φz,t Sz,t

d X

!2  qz(i) φz



z=1

2 ,

where for 0 ≤ t∧τz , where τz is theRtime at which samples τ of cohort z are censored, φz,t = t z ds Sz,s , φz ≡ φz,0 , Sz,t is the survival function for the pooled data of cohort z, and ! (1) (2) (2) (1) p(1) Cz,t− qz + p(2) Cz,t− qz Wz,t = . (2) (1) Cz,t− Cz,t− The variance σ 2 may be consistently estimated by  !2  2 d d X X X σ ˆ2 = (1 − p(i) )  qˆz(i) φˆ2z − qˆz(i) φˆz  z=1

i=1

+

d Z τz X z=1

0

ˆ z,t dSˆz,t W

z=1

φˆz,t Sˆz,t

!2 ,

(2.5)

where for 0 ≤ t ∧ τz , Sˆz,t is the product-limit estimate of the pooled data for cohort z, Z τz φˆz,t = ds Sˆz,s , (2.6) t (i) Cˆz,t is the product-limit estimate associated to the event of censoring for cohort z within population i, φˆz ≡ φˆz,0 , and ! (2) (1) (1) (2) p(1) Cˆz,t− qˆz + p(2) Cˆz,t− qˆz ˆ z,t = W . (1) (2) Cˆz,t− Cˆz,t−

3

At any given time t ≥ 0, each particle will have exactly one associated state x in the set Z4 , referring respectively to states of (1) dormancy, (2) activity, (3) inactivity, (4) censored. Assume that the path of any particle is statistically dependent upon its particular subsystem, and that given the respective subsystems of any two particles, their resulting paths are statistically independent. Assume further that at a reference time t = 0, all particles enter into the active state (x = 2). Let d ∈ N and τ ∈ (0, ∞) be fixed. We will assume the existence of a collection of individuals Γ , assumed to be infinite in size, where each individual γ ∈ Γ exhibits a c`adl`ag path-valued state xγt , for t ≥ 0. For each γ ∈ Γ , xγt is determined by the individuals particle type z γ and a random jump time ξ γ . The particle type z γ is distributed in the population through the probability γ mass P(zP = z) = qz , where q = (q1 , . . . , qd ) ∈ (0, 1)d d satisfies z=1 qz = 1. Let St = (S1,t , . . . , Sd,t ) be the survival vector Sz,t = P {Tz > t}, which is assumed continuous for t ≥ 0. Suppose that it is desired to understand the event probabilities for randomly selected γ ∈ Γ . Given a random sample γ1 , . . . , γn , n ∈ N of individuals, let n = (n1 , . . . , nd ) where nz is the number of drawn of cohort z. In considering the event time probabilities of each subgroup, the random number of initial particles excludes the use of many well established results in survival analysis. Therefore, we begin with a somewhat restricted framework, and assume a known number of initial individuals of each type. Assume the sample contains a known number nz of individuals of cohort z, and let µjz,t be the number of the initial cohort z individuals who are in state j ∈ Z4 at time t. Denote the z-type cumulative hazard by Λz,t and respectively define the z-type cumulative hazard and survival estimates by Z t dµ3z,s Λˆz,t = (2.7) 2 0 µz,s−  Y Sˆz,t = 1 − dΛˆz,s . (2.8) s≤t

2.1 Empirical process framework To prove Theorem 1, we turn to a modeling framework that will provide us the asymptotic statistics of the product limit estimator. In the appendix, we build on this framework in a series of Lemmata in order to prove the main result. Consider a closed particle-system, such that according to a predefined set of characteristics, the system can be subdivided into mutually exclusive subsystems. Note that we will restrict this discussion to only a single population or particles. These arguments will extend to multiple populations as mentioned in this manuscript by treating separate populations as independent.

Define further Bz,t =

√

nz

Sˆz,t − Sz,t Sz,t

ˆz,t = B ˆz,τ for all t ≥ τz . and note that Sˆz,t = Sˆz,τz and B z From [5], it follows that {Bz,t : t ≥ 0} is a mean-zero square-integrable martingale with Meyer bracket process !2 Z t∧τz 1{µ2 >0} Sˆz,s− z,s− dΛz,s , hBz , Bw it = δzw nz 2 S µ z,s 0 z,s− where t ∧ τz = min{t, τz }, and δ(·,·) is the Kroenicker delta function.

4

Aaron Heuser et al.

q2 = 0

q2 = 0.25

q2 = 0.30 1 survival (St )

density (πt )

0.8 0.6 0.4 0.2 0

q2 = 0.35

0

5

10

15

time

0.5

0

0

5

10

15

time

Fig. 2 Admixture test distributions used in simulated investigations of our estimator. Populations formed using q2 ∈ [0, 1) admixtures of (1 − q2 )exponential(λ = 5−1 ) and q2 Weibull(k = 5, λ = 1) event time distributions. Event time density functions πt and corresponding survival functions St are shown for various values of q2 .

Proof (of Theorem 1) To prove the main theorem, we build upon the modeling framework previously mentioned, starting with a deterministic sample size N ∈ N and replacing it at the cohort-level with sample sizes given by a random vector N ∈ Zd , where Nz = N az for a = (a1 , . . . , ad ) ∈ ∆d , the d−dimensional unit simplex, chosen in a sufficiently small neighborhood V of q. We define Sˆz,t (az ), Λˆz,t (az ), and Bz,t (az ), as above, but under the assumption that the initial number of ztype individuals is Nz . Note that replacing az with qz will describe the case of random subsystems in the main theorem. Therefore, work will first be done in the restricted case. Then through an application of the MannWald Theorem (cf. [3]), results for the case of random subsystems will be derived. Considering the above model for the empirical survival function, convergence of the statistic in Eq. 2.4 follows immediately from Corollary 10 in the appendix. The consistency of σn2 follows from theorem 4.2.2 of [5] and the Glivenko-Cantelli theorem. t u 3 Simulated investigation Our overall statistical method consists of two parts. The first part is an alternative method of piecing together cohort-level survival curves to assemble an overall survival function for a population. The second part is a method for comparing survival functions created in this way based on the area between the curves. To examine our overall method, we turned to simulations. As test populations, we examined admixtures of exponential and Weibull distributions for the event time, and compared survival in these mixture populations to survival of a population of purely exponential event times (Fig. 2). Population 1 consists of individuals having an exponentially distributed lifetime with a mean of λ−1 = 4 years. Population 2 consists of two types of individuals: those who have an exponentially distributed lifetime

with a mean of 5 years (type z = 1), and those of type z = 2 who have a Weibull distributed lifetime with shape parameter k = 5 and scale parameter λ = 1. Since Population 1 is homogeneous, we only track subpopulations of Population 2 - we drop the superscript and denote the proportion of Population 2’s members of type 2 by q2 . It is most instructive to examine our method in the neighborhood where both populations have approximately the same expectation value for the event time, which occurs for q2 ≈ 0.245. For this reason, we chose values near 0.25 for our simulations. To compare our new reweighted Kaplan-Meier method (Eq. 2.2) to the standard Kaplan-Meier estimator, we estimated survival for the admixed population for q2 = 0.25, using various sample sizes. In Fig. 3, we present example reconstructions using these two methods. The estimator variance was approximated using 10, 000 resamplings of sample size n of the admixed population, for each value of n. The estimation error, as defined by mean-squared difference between the reconstruction and the true survival function, was approximated in the same manner. To better-understand the performance of our test statistic (Eq. 2.4), we evaluated its statistical power against that of other test statistics in distinguishing between Population 1 and Population 2 for various values of q2 . For samples of size n(i) ∈ {30, 50, 100, 200, 1000} taken from each population, we performed 1000 null hypothesis statistical tests using our method, the log-rank method [2], and the standard Kaplan-Meier Wilcoxon signed-rank difference-of-mean methods [15, 18]. The power of the test, or the proportion of times that the null hypothesis was correctly rejected, is shown in Fig. 4.

4 Discussion and Conclusion In this manuscript we have proposed a test statistic that uses a cohort-averaged survival function estimator in order to make cross-population comparisons of survival within a null hypothesis significance testing framework. The proposed survival estimator was an empirically-weighted average of cohort-level product-limit estimates. The test statistic involved computation of the area between estimated survival functions for two populations. By invoking an empirical stochastic process, we proved asymptotic normality of this test statistic. Using simulations, we contrasted the weighted survival estimator against the pure Kaplan-Meier estimator. It is seen, in Fig. 3, that this new estimator has comparable performance to the pure Kaplan-Meier estimator at large sample sizes. Asymptotically, both estimators converge to the true survival function. At small sample sizes, there are differences. This new estimator appears to have smaller variance at the cost of larger error

A Survival Estimator and Statistic for Nonparametric Heterogeneous Survival Analysis

Our method (θˆt )

estimated survival

n = 10

5

Pure Kaplan-Meier (Sˆt )

n = 20

n = 40

n = 80

1

1

1

1

0.8

0.8

0.8

0.8

0.6

0.6

0.6

0.6

0.4

0.4

0.4

0.4

0.2

0.2

0.2

0.2

0

0

0

0

5

10

15

0

estimator variance

·10−2

10

15

0

·10−2

5

10

15

0

0

5

10

15

5

10

15

5

10

15

·10−3

·10−3 6

3

4

2

2

1

1

2

0.5

1

0

0

0 0

5

10

0

15

·10−2

3 estimator error

5

1

1

0.5

10

0

5

10

15

5

10

15

·10−3

6

3

4

2

2

1 0

0 0

0

·10−3

0

0

0 0

15

·10−2

1.5

2

5

5

10

15 time

0

5

10

15

0

Fig. 3 Comparing estimators of survival. The survival estimation method of Eq. 2.2 compared to pure Kaplan-Meier for a population containing an admixture of (1 − q2 )Exponential(1/5) and q2 Weibull(1, 5) individuals, where q2 = 0.25. At a given sample size n, the survival estimates are obtained (top row: examples shown and contrasted). The estimator variance and mean square error were approximated using 10, 000 resamplings for each of the sample sizes.

Our method

Kaplan-Meier mean-difference q2 = 0.30

q2 = 0.25

power

0.8 0.6 0.4 0.2 0

500

Log-rank

1,000

q2 = 0.35

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2 0

500

1,000

0

500

1,000

sample size Fig. 4 Simulated power computation comparing exponentially distributed lifetimes against a mixture of q2 Weibull and (1 − q2 ) exponential distributions, where q2 determines the amount of mixing. A larger value of q2 implies more real difference between the survival functions of the two populations. The power of our method (black) is compared to the power of Kaplan-Meier difference-of-mean (blue) and Log-rank (red) methods. (More power is better).

6

at earlier times. This error at earlier times is mitigated by decreased error in tail estimation. Hence, dependent on costs, for small samples, this new estimator may be preferable to the pure Kaplan-Meier estimator. In simulations of the test statistic derived from the new survival estimator, we saw superior performance compared to existing methods. In Fig. 4, it is seen that in ˆ was better at distinguishing all cases, the test statistic Θ between the two populations than either the Wilcoxon signed-rank test or the log-rank test. The relatively-high statistical power of this new statistic is due to tighter variation in the test-statistic. In nearly all cases (> 99.5%), the estimator variance for our method was less than that of the other two tests (not shown). A variant of this method was used in Rasch et al [14] in order to classify physical disorders based on severity for the sake of prioritization of processing for disability claims. Since the underlying survival surface is non-stationary, and the fixed observation windows create progressive censoring, that paper illustrates the utility of this statistical method. In that manuscript, the cohorts were defined based on binned application times and a heuristic “survival surface” was estimated. Although the most direct and natural applications of the method that we have presented here involve discretely-indexed covariates, it is possible to use this method for continuously-indexed covariates such as time by employing the binning strategy used in Rasch et al [14]. This approach is particularly fruitful if the sampling windows are coarse and there is clear separation between cohorts to maintain statistical independence. In this situation, it may be unreasonable to expect to reconstruct a full continuous surface for survival. Nonetheless, a possible future extension of this method might involve replacing the sum of Eq. 2.1 with an integral and using statistical regularization tools [4] in order to infer true continuously-indexed surfaces.

5 Acknowledgements The authors would like to thank the United States Social Security Administration for providing opportunity and funding for this research, Dr. Leighton Chan and Dr. Elizabeth Rasch for insightful discussions, guidance and support, Dr. Pei-Shu Ho for help obtaining data, and the Epidemiology and Biostatistics section at NIH Clinical Center, Rehabilitation Medicine Department for supporting our work. This research was supported in part by the Intramural Research Program of the NIH, Clinical Center.

Aaron Heuser et al.

References 1. Aalen OO (1994) Effects of frailty in survival analysis. Statistical Methods in Medical Research 3(3):227–243 2. Berty HP, Shi H, Lyons-Weiler J (2010) Determining the statistical significance of survivorship prediction models. Journal of evaluation in clinical practice 16(1):155–165 3. Billingsley P (1968) Convergence of Probability Measures. John Wiley and Sons, New York 4. Chang JC, Savage VM, Chou T (2014) A pathintegral approach to bayesian inference for inverse problems using the semiclassical approximation. Journal of Statistical Physics 157(3):582–602 5. Gill R (1980) Censoring and Stochastic Integrals. Mathematisch Centrum, Amsterdam 6. Gill R (1981) Large Sample Behaviour of the Product-limit Estimator on the Whole Line. MC. SW, Stichting Math. Centrum 7. Hougaard P (1984) Life table methods for heterogeneous populations: distributions describing the heterogeneity. Biometrika 71(1):75–83 8. Kaplan E, Meier P (1958) Nonparametric estimator from incomplete observations. Journal of the American Statistical Association 53(282):457–481 9. Karatzas I, Shreve S (2000) Brownian Motion and Stochastic Calculus. Springer-Verlag, New York 10. Lenglart E (1977) Relation de domination entre deux processus. Annales de l’Institut Henri Poincar´e 13:171–179 11. Lin D, Ying Z (1993) A simple nonparametric estimator of the bivariate survival function under univariate censoring. Biometrika 80(3):573–581 12. Pepe M, Fleming T (1989) Weighted kaplan-meier statistics: A class of distance tests for censored survival data. Biometrics 45(2):497–507 13. Pepe MS, Fleming TR (1991) Weighted kaplan-meier statistics: Large sample and optimality considerations. Journal of the Royal Statistical Society Series B (Methodological) 52:341–352 14. Rasch EK, Huynh M, Ho PS, Heuser A, Houtenville A, Chan L (2014) First in line: prioritizing receipt of social security disability benefits based on likelihood of death during adjudication. Medical care 52(11):944–950 15. Schoenfeld D (1981) The asymptotic properties of nonparametric tests for comparing survival distributions. Biometrika 68(1):316–319 16. Shorack G, Wellner J (1986) Empirical processes with applications to statistics. Wiley, New York 17. van der Vaart AW (1996) New donsker classes. The Annals of Probability 24(4):2128–2140 18. Wilcoxon F (1945) Individual comparisons by ranking methods 1(6):80–83

A Survival Estimator and Statistic for Nonparametric Heterogeneous Survival Analysis

7

A Preliminaries and Notation

Since any survival function is monotone, an immediate result that follows from the above is assumption is

Given any pair of random elements X, Y , we denote equality in a distributional sense by X ≈ Y . Let P be a probability measure on the measurable space (X, A). The empirical measure generated by the sample of random elements x1 , . . . , xn , n ∈ N is given by

c < Sz,τz ≤ Sz,t ≤ 1,

Pn = n−1

n X

δxi ,

(A.1)

where for any x ∈ X, and any A ∈ A,  1, x ∈ A, δx (A) = 0, x ∈ / A.

(A.2)

Note that alternatively, when needed, one may write δx (A) as the indicator function 1A (x) on the set A. Furthermore, in the case that A = {k}, k ∈ Z, and x ∈ Z, we write δx (A) ≡ δx,k . Given H, a class of measurable functions h : X → R, the empirical measure generates the map H → R given by h 7→ Pn h, where for any signed measure Q and measurable function R h, we use the notation Qh = dQ h. Furthermore, define the H-indexed empirical process Gn by Gn h =

√

n 1 X n (Pn − P) h = √ (h(Ai ) − Ph) , n i=1

(A.3)

and with thePempirical process, identify the signed measure Gn = n−1/2 n i=1 (δAi − P). Note that for a measurable function h, from the law of large numbers and the central limit theorem, it follows that
d → N 0, P (h − Ph)2 , provided Ph Pn h →a.s Ph, and Gn h − d

exists and Ph2 < ∞, and where “− →” denotes convergence in distribution. In addition to the preceding notation, given the elements f , and fn , n ∈ N, we also denote respectively, convergence in probability and in distribution, of fn to f , by p fn − → f. For any map x : H → Rk , k ∈ N, define the uniform norm kxkH by kxkH = sup {|x(h)| : h ∈ H},

(A.4)

 `∞ (H) = x : H → Rk : kxkH < ∞ . d

If for some tight Borel measurable element G ∈ `∞ (H), Gn − → G, in `∞ (H), we say that H is a P-Donsker class.

B Convergence Theorems In order to guarantee convergence of the estimator, we make the following assumptions (based upon an initially known sample size distribution). Assumption 2 We assume that the initial sample is chosen large enough to ensure that individuals of cohort z, at state 1, exist at all points t ∈ [0, τz ], z ∈ {1, . . . , d}. That is, inf µ1,z,τz − > 0,

Assumption 3 It is assumed that as n becomes large, the sample size for each individual type will grow to infinity. That is, lim

inf

a.s.

µ1,z,τz − = ∞,

a.s.

Assumption 4 For each z ∈ {1, . . . , d} there exists a nonincreasing continuous function mz : [0, ∞) → (0, 1] such that µ1,z,t lim sup − mz,t = 0 a.s. n→∞ t≥0 naz Note that in the case of fixed censoring, that is, in the case that censoring exists only at time τ , the above is satisfied by mz,t = Sz,t . In the general case, mz,t can be seen as the probability that an individual of cohort z has not yet left state 1. That is, mz,t is the probability that an individual has not left due to censoring or death by time t, and so mz,t = Sz,t Cz,t− , where Cz,t is the probability that censoring has not occured by time t. To prove the main theorem, we now present a series of lemmata. Lemma 5 If qˆ and Sˆz,s− (ˆ qz ) are defined as in the previous section, then Z t∧τz d   p √ X n (ˆ qz − qz ) ds Sˆz,s− (ˆ qz ) − Sz,s − → 0, z=1

0

as n → ∞, uniformly in t ≥ 0. Proof It is claimed that to prove the statement of the lemma, it suffices to show that !2 Sˆz,t− (ˆ qz ) − Sz,t p sup − → 0, (B.2) Sz,t t≥0 uniformly in t ≥ 0, for each z = 1, . . . , d. Indeed, for if the above holds, then Z t∧τz   p qz ) − Sz,s − → 0, ds Sˆz,s− (ˆ 0

and in the case that H ⊂ R, write k·kH ≡ k·k∞ . A class H for which kPn − PkH → 0 is called a P-Glivenko-Cantelli class. Denote by `∞ (H) the class of uniformly bounded functions on H. That is, for a general k ∈ N,

z∈Nd

(B.1)

for some constant c > 0.

n→∞ z∈Nd ,a∈V

i=1

t ≥ 0,

uniformly in t ≥ 0. Since the central limit theorem implies that √ d n(ˆ qz − qz ) − → N (0, qz (1 − qz )), each term in the sum would converge in probability to 0, uniformly in t ≥ 0. And so, if EN denotes the expectation given N , we have that !2 Sˆz,t− (ˆ qz ) − Sz,t 1 E = E EN (Bz,t (N ))2 Sz,t N ! Z t∧τz 1 dΛz,s Sˆz,s− = E EN N N µ1z,s− Sz,s 0 ! Z t∧τz dΛz,s Sˆz,s− =E 1 µ Sz,s 0 z,s− ≤ CE (µ1,z,τz )−1 , for some arbitrary constant C. From Lenglart’s inequality (cf. [10]), ( !2 )   Sˆz,t− (ˆ qz ) − Sz,t η C >  ≤ + P µ1,z,τz − < , P sup Sz,t  η t for any arbitrary η,  > 0. Therefore, from Assumption 3, since Nz → ∞ a.s., the desired result follows.

8

Aaron Heuser et al. d d Then ζˆt1 − → ζt1 and ζˆt2 (a) − → ζt2 (a), in the space of compactly supported functions DRd [0, ∞) as n → ∞, for each a ∈ V , 1 1 where ζt1 = (ζ1,t , . . . , ζd,t ) is the mean-zero square-integrable Gaussian process defined by

For any t ≥ 0, r Z t∧τ √ 1,β ˆ t = n2 Θ ds n1 (θˆs− − θs1,β ) n 0 r Z n1 t∧τ √ 2,β − θs2,β ). ds n2 (θˆs− − n 0



1 ζz1 , ζw



(B.8)

t

Z

ˆ For a general survival function θ, with respective estimate θ, define Yˆt by Z t∧τ  √  β Yˆt = ds n θˆs− (B.3) − θsβ , t ≥ 0. 0

t∧τz

= −qz qw

t∧τw

 Z

ds Sw,s

ds Sz,s 0

0

Z



2

t∧τz

+ δz,w qz

ds Sz,s

,

0 2 2 (a), . . . , ζd,t ) is given by and ζt2 (a) = (ζ1,t

If the p process Yˆ converges in distribution to some Y ∼ N (0, σ ), since ni /n converges to p(i) , i = 1, 2, it follows that 2

d

ˆt − Θ →

p

p 2 2 p(2) Yt1 − p(1) Yt2 ≈ N (0, p(2) σ1,t + p(1) σ2,t ).

√

Z

t∧τz

  ds qˆz Sˆz,s− (ˆ qz ) − qz Sz,s 0 Z t∧τz   √ = n(ˆ qz − qz ) ds Sˆz,s− (ˆ qz ) − Sz,s 0 Z t∧τz √ + n(ˆ ds Sz,s qz − qz ) 0 Z t∧τz √ + nqz ds (Sˆz,s− (ˆ qz ) − Sz,s ) n

Z

t∧τz

˜z,s φz,s − φz,t∧τ B ˜z,t∧τ dB z z

 .

(B.9)

0

The processes ζˆ1 and ζˆ2 (a) are independent, and there exist a Skorohod representations such that

P ˆ Note that Yˆt = d z=1 Zz,t , where ˆz,t = Z

qz 2 ζz,t (a) = √ az

(B.4)

1 1 sup ζˆz,t − ζz,t → 0,

t≥0

and sup t≥0,a∈V

ˆ2 2 (a) → 0, ζz,t (a) − ζz,t

almost surely as n → ∞.

0

Therefore, if it can be shown that Z d √ X n (ˆ qz − qz )

t∧τz

  p ds Sˆz,s− (ˆ qz ) − Sz,s − → 0,

0

z=1

uniformly in t, then convergence of (Yˆt : t ≥ 0) is dependent only upon the convergence of the d-dimensional vector-valued ˆ q ) given by process ζ(ˆ ζˆz,t (a) =

√

Z

sup a∈V,t≥0

t∧τz

n(ˆ qz − qz ) ds Sz,s 0 Z t∧τz √ ds (Sˆz,s− (az ) − Sz,s ), + nqz

Proof To begin note that independence follows immediately from the independence of the respective limiting processes. Since N is a multinomial random variable, (B.8) follows from the central limit theorem. In the case of ζˆt2 (a), we first consider Bt (a). An application of Lenglart’s inequality, very similar to that in the proof of Lemma 5, along with Assumption 3, shows that

(B.5) sup a∈V,t≥0

with a = (a1 , . . . , ad ) ∈ (0, 1)d chosen in a sufficiently small neighborhood V of q. This decomposition will thus lead to the main theorem. To show the desired convergence of ζˆt (ˆ q ), we first focus on convergence of ζˆt (a). R Let φz,t = tτz ds Sz,s and write ζˆt (a) = ζˆt1 + ζˆt2 (a), where √

(−dφz,s ),

−

1,z,t−

1 p − → 0, mz,t

naz µ1,z,t−

Sˆz,t− (a) Sz,t

!2

(B.6) p

and t∧τz

(−dφz,s )Bz,s (az ),

(B.7)

0

Lemma 6 Suppose that

1 , mz,t

p

− →

hBz (a), Bw (a)it − → δz,w Z

as n → ∞.

It follows that

0

qz 2 ζˆz,t (a) = √ az

naz µ

uniformly in t ≥ 0, and since mz,t = Sz,t Cz,t− ,

t∧τz

Z n (ˆ qz − qz )

as n → ∞.

Moreover, from Assumption 4,

0

1 ζˆz,t =

ˆ p → 0, Sz,t− (a) − Sz,t −

n o n o ζˆt1 (a) : t ≥ 0 and ζˆt2 (a) : t ≥ 0

are the processes respectively defined by equations (B.6) and ˜ is the d-dimensional mean-zero Gaussian (B.7), and that B process defined by Z t∧τz D E dΛz,s ˜z , B ˜w = δz,w B . t S z,s Cz,s− 0

Z

t∧τz

0

dΛz,s . Sz,s Cz,s−

d ˜ Therefore, from theorem 4.2.1 of [5], B(a) − → B, and there exists a Skorohod representation of B(a) such that

sup t≥0,a∈V

˜z,t → 0, Bz,t (a) − B

almost surely as n → ∞. Since almost sure convergence of B(a) implies almost sure convergence of bounded functionals of B(a), the desired convergence of ζˆ2 (a) follows from Theorem 2.1 of [6].

A Survival Estimator and Statistic for Nonparametric Heterogeneous Survival Analysis n o ˆ Corollary 7 If the process ζ(a) = ζˆt (a) is defined by equation (B.5), then d X

d X

d ζˆz (a) − →

ζz,t (a)

=

d X

Since Sτz > 0 for all z, from Doob’s martingale inequality (cf. [9]), d X

(B.10) E sup

z=1

z=1

t≥0 1 ζz,t

+

1 Proof From the previous theorem we may assume that ζˆz,t → 1 2 2 ˆ ζz,t and ζz,t (a) → ζz,t (a) almost surely, uniformly for a ∈ V and t ≥ 0. Therefore

1 1 −√ √ az bz

2

Define the map g : V × ` g(a, f ) = f (a, ·), then ζz,t

z=1

∞

d N X , ζz n z=1

=g

(V × [0, ∞)) → `

∞

p 1 √ ( az − bz )2 az bz √ 2 √ p √ az + bz √ ≤ δ −2 ( az − bz )2 √ az + bz 1 ≤ 3 (az − bz )2 . 4δ

Combining the above two results gives

([0, ∞)) by d X

E sup

!

t≥0

.

Furthermore, if for any (a1 , f1 ), (a2 , f2 ) ∈ V × `∞ (V × [0, ∞)) we have that |a1 − a2 | +

2 d  X 1 1 −√ , √ az bz z=1

=

p

Since N/n − → q, from Theorem 4.4 of [3] ( d )! ( d )! X X N d , ζˆz,t (a) : t ≥ 0 − → q, ζz,t (a) : t ≥ 0 . n z=1 z=1



≤C

z=1

z=1

almost surely, uniformly for a ∈ V and t ≥ 0. The statement of the theorem then follows from theorem 5.1 of [3].

N n

!2 ζz,t (b)

for some arbitrary constant C. For each z ∈ Nd , since az and bz are sufficiently close to qz ∈ (0, 1), it follows that there exists some δ > 0 such that az ∧ bz > δ. Therefore, 

ζˆt (a) → ζt (a)



d X

ζz,t (a) −

2 ζz,t (a).

z=1

d X

9

d X

ζz,t (a) −

!2 ≤ C |a − b|2 ,

ζz,t (b)

z=1

z=1

and so, by Kolmogorov’s continuity criterion (cf. [9]), the desired result follows.

|f1 (a, t) − f2 (a, t)| < δ

sup a∈V,t≥0

The above lemma, along with the argument immediately preceding, gives the following.

for some δ > 0, then sup |g(a1 , f1 )(t) − g(a2 , f2 )(t)| t≥0

P Pd n Theorem 9 Let d z=1 ζz,t (·) and z=1 ζz,t (·) be defined as in Corollary 7, then

= sup |f1 (a1 , t) − f2 (a2 , t)| t≥0

≤ sup |f1 (a1 , t) − f1 (a2 , t)| + sup |f1 (a2 , t) − f2 (a2 , t)| . t≥0

t≥0

Therefore, g is a continuous at any (a, f ) such that f is continuous at a, uniformly in t. It thus follows P from the continuous mapping theorem (cf. [17]) that if a 7→ d z=1 ζz,t (a) is continuous, uniformly in t, then ! ! d d X N Xˆ d , ζz − → g q, ζz . (B.11) g n z=1 z=1 Lemma 8 If {ζt (a) : t ≥ 0} is defined as in Corollary 7, then the map a 7→

d X

d X

ζˆz,t



z=1

N n



d X

d

− →

(B.12)

P Corollary 10 If ζˆ = d z=1 ζz,τz (q), then ζˆ ∼ N (0, σ 2 ), where

ζz,t (a)

z=1

is continuous for a ∈ V , uniformly in t ≥ 0.

σt2

=

d X

qz φ2z,0

−

z=1

=

d X z=1

d X

−

d X

Z

τz

qz 0

z=1

ζz,t (b)

qz φz,0

z=1

Proof For any a, b ∈ V , it follows that ζz,t (a) −

!2

d X

z=1

d X

in DR [0, ∞), as n → ∞.

ζz,t (q),

z=1

dSz,t Cz,t−



φz,t Sz,t

z=1

 qz

1 1 −√ √ az bz



Proof Note that when t = τz , we have

 Z ˜z,t∧τ − × φz,t∧τz B z

t∧τz

0

 ˜z,s φz,s . dB

1 ζz,τz (q) = ζz,τ + z

√

Z

τz

qz 0

˜z,t φz,t , dB

2

10

Aaron Heuser et al.

which are independent and normally distributed, implying that ζˆ is also normally distributed. Furthermore Eζˆ2 =

Z d  X √ 1 ζz,τ + qz z z=1

τz

˜z,t φz,t dB

2

0

 Z d X √ 1 + ζz,τ + q z z

τz

˜z,t φz,t dB



0

z,w=1 z6=w

 Z √ 1 + × ζw,τ q z w

τw

˜w,t φw,t dB



0

Z d  X
2 1 + Eq = E ζz,τ z z

τz

˜z,t φz,t dB

2 

0

z=1

−

d X

1 Eζz,τ ζ1 z w,τw

z,w=1 z6=w

=

d X

qz (1 −

qz )φ2z,0

Z

τz

dΛz,t

+ qz 0

z=1

−

d X

φ2z,t

!

Sz,t Cz,t−

qz qw φz,0 φw,0 ,

z,w=1 z6=w

which after recombining the final terms, gives the desired result. The above result completes the proof of Theorem 1, and thus completes the body of this paper.