Asymptotic convergence in distribution of the area bounded by prevalence-weighted Kaplan-Meier curves using empirical process modeling

TL;DR

This paper analyzes the asymptotic distribution of the area between prevalence-weighted Kaplan-Meier survival curves, providing theoretical insights into their large sample behavior for population comparisons.

Contribution

It extends empirical process modeling to the area between weighted Kaplan-Meier curves, offering new theoretical results for population-level survival analysis.

Findings

Derived the asymptotic distribution of the area statistic

Applied empirical process theory to weighted survival curves

Supported by NIH-SSA medical condition prioritization case

Abstract

The Kaplan-Meier product-limit estimator is a simple and powerful tool in time to event analysis. An extension exists for populations stratified into cohorts where a population survival curve is generated by weighted averaging of cohort-level survival curves. For making population-level comparisons using this statistic, we analyze the statistics of the area between two such weighted survival curves. We derive the large sample behavior of this statistic based on an empirical process of product-limit estimators. This estimator was used by an interdisciplinary NIH-SSA team in the identification of medical conditions to prioritize for adjudication in disability benefits processing.

Equations

$${\theta }_t^{(i)}=\sum _{z=1}^d\underset{S_{z,t}^{(i)}}{\underbrace{\mathbb{P}\left\{T^{\gamma }>t\mid z^{\gamma }=z,\gamma \in {\Gamma }^{(i)}\right\}}}\underset{q_z^{(i)}}{\underbrace{\mathbb{P}\left\{z^{\gamma }=z\mid \gamma \in {\Gamma }^{(i)}\right\}}},$$

$${\hat{\theta }}_t^{(i)}=\sum _{z=1}^d{\hat{q}}_z^{(i)}{\hat{S}}_{z,t}^{(i)},$$

$${\hat{q}}_z^{(i)}=\frac{n_z^{(i)}}{n^{(i)}}.$$

$$\hat{\Theta }=\sqrt{\frac{n^{(1)}{n}^{(2)}}{n^{(1)}+n^{(2)}}}{\int }_0^{\tau }\mathrm{d}t\left({\hat{\theta }}_t^{(1)}-{\hat{\theta }}_t^{(2)}\right),$$

$${\sigma }^2$$

$$-\sum _{z=1}^d{\int }_0^{{\tau }_z}\mathrm{d}{S}_{z,t}{W}_{z,t}\times {\left(\frac{{\varphi }_{z,t}}{S_{z,t}}\right)}^2,$$

$$W_{z,t}=\left(\frac{p^{(1)}{C}_{z,t-}^{(1)}{q}_z^{(2)}+p^{(2)}{C}_{z,t-}^{(2)}{q}_z^{(1)}}{C_{z,t-}^{(1)}{C}_{z,t-}^{(2)}}\right).$$

$${\hat{\sigma }}^2$$

$$-\sum _{z=1}^d{\int }_0^{{\tau }_z}\mathrm{d}{\hat{S}}_{z,t}{\hat{W}}_{z,t}\times {\left(\frac{{\hat{\varphi }}_{z,t}}{{\hat{S}}_{z,t}}\right)}^2,$$

$${\hat{\varphi }}_{z,t}={\int }_t^{{\tau }_z}\mathrm{d}s{\hat{S}}_{z,s},$$

$${\hat{W}}_{z,t}=\left(\frac{p^{(1)}{\hat{C}}_{z,t-}^{(1)}{\hat{q}}_z^{(2)}+p^{(2)}{\hat{C}}_{z,t-}^{(2)}{\hat{q}}_z^{(1)}}{{\hat{C}}_{z,t-}^{(1)}{\hat{C}}_{z,t-}^{(2)}}\right).$$

$${\sigma }^2$$

$$-\sum _{z=1}^d{\int }_0^{{\tau }_z}\mathrm{d}{S}_{z,t}{W}_{z,t}\times {\left(\frac{{\varphi }_{z,t}}{S_{z,t}}\right)}^2,$$

$$W_{z,t}=\left(\frac{p^{(1)}{C}_{z,t-}^{(1)}{q}_z^{(2)}+p^{(2)}{C}_{z,t-}^{(2)}{q}_z^{(1)}}{C_{z,t-}^{(1)}{C}_{z,t-}^{(2)}}\right).$$

$${\hat{\sigma }}^2$$

$$-\sum _{z=1}^d{\int }_0^{{\tau }_z}\mathrm{d}{\hat{S}}_{z,t}{\hat{W}}_{z,t}\times {\left(\frac{{\hat{\varphi }}_{z,t}}{{\hat{S}}_{z,t}}\right)}^2,$$

$${\hat{\varphi }}_{z,t}={\int }_t^{{\tau }_z}\mathrm{d}s{\hat{S}}_{z,s},$$

$${\hat{W}}_{z,t}=\left(\frac{p^{(1)}{\hat{C}}_{z,t-}^{(1)}{\hat{q}}_z^{(2)}+p^{(2)}{\hat{C}}_{z,t-}^{(2)}{\hat{q}}_z^{(1)}}{{\hat{C}}_{z,t-}^{(1)}{\hat{C}}_{z,t-}^{(2)}}\right).$$

$${\mathbb{P}}_n$$

$${\delta }_x(A)$$

$${\mathcal{G}}_{n}h$$

$${\Vert x\Vert }_{\mathcal{H}}$$

$${\ell }^{\infty }(\mathcal{H})$$

$$\begin{array}{cc}0& \text{dormancy}\\ 1& \text{activity}\\ 2& \text{inactivity}\\ 3& \text{censored}.\end{array}$$

$$n=\sum _{z=1}^d{n}_z$$

$$\sum _{j=0}^3{\mu }_{j,z,t}^{n_z}=n_z$$

$$n_z={\mu }_{2,z,t}^{n_z}+{\mu }_{3,z,t}^{n_z}\forall t>{\tau }_z$$

$${\hat{\Lambda }}_{z,t}^{n_z}$$

$${\hat{S}}_{z,t}^{n_z}$$

$$B_{z,t}^{n_z}=\sqrt{n_z}\frac{{\hat{S}}_{z,t}^{n_z}-S_{z,t}}{S_{z,t}}$$

$${⟨B_{z,t}^{n_z},B_{w,t}^{n_z}⟩}_t={\delta }_{zw}{n}_z{\int }_0^{t\wedge {\tau }_z}\mathrm{d}{\Lambda }_{z,s}{\left(\frac{{\hat{S}}_{z,s-}^{n_z}}{S_{z,s}}\right)}^2\frac{1_{\left\{{\mu }_{1,z,s-}^{n_z}>0\right\}}}{{\mu }_{1,z,s-}^{n_z}},$$

$$\underset{z\in {\mathbb{N}}_d}{\inf }{\mu }_{1,z,{\tau }_z-}^{n_z}$$

$$c<S_{z,{\tau }_z}\le S_{z,t}\le 1,t\ge 0,$$

Full text

{fmtext}

1 Introduction

Survival analysis addresses the classical statistical problem of determining characteristics of the waiting time until 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 Kaplan and Meier [1958]. Their product-limit estimator Gill [1980], Gill et al. [1983], Van Der Vaart et al. [1996], Shorack and Wellner [2009] 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 are only partially observed because they fall outside the observation window. It was not, however, designed to account for varying subpopulations that may yield non-homogeneity 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. Hougaard [1984]), and extended by Aalen (cf. Aalen [1994]), 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. Lin and Ying [1993]), 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.

These complexities arose in the identification of new disorders to incorporate into the United States Social Security Administration (SSA)’s Compassionate Allowances (CAL) initiative. The CAL initiative seeks to identify candidate medical conditions for fast-tracking in the processing of disability applications. The intent of this initiative is to prioritize applicants who are most likely to die in the time-course of usual case processing so that they may receive benefits while still living.

At its inception, the CAL initiative identified conditions based on the counsel of expert opinion Rasch et al. [2014]. The SSA in collaboration with the National Institutes of Health (NIH) sought to expand the list of CAL conditions systematically, using a data-based approach. Using in-part the survival estimator described in this manuscript, the NIH identified 24 conditions for inclusion into the list of conditions Rasch et al. [2014].

The methodology used in CAL is related to that of the work of Pepe and Fleming (cf. Pepe and Fleming [1989, 1991]), where 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 Pepe and Fleming [1991] 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.

2 Statistical method

Suppose $\Gamma^{(1)}$ and $\Gamma^{(2)}$ are disjoint populations of individuals where each individual belongs to exactly one of $d$ distinct cohorts labeled $z\in\mathbb{Z}_{d}$. For randomly selected individuals $\gamma\in\Gamma^{(i)}$ within population $i$, we desire to understand the statistics of the event time $T^{\gamma}$ under the assumption that survival is conditional on cohort $z^{\gamma}$ and population.

One representation of the marginal survival probability for members of population $i$, $\theta^{(i)}_{t}=\mathbb{P}{\left\{T^{\gamma}>t\mid\gamma\in\Gamma^{(i)}\right\}},$ is found by conditioning on cohort,

[TABLE]

where $S^{(i)}_{z,t}$ represents the survival function for individuals of cohort $z$ in population $i$, where each individual’s cohort membership is known.

We use this representation of the survival probability as motivation to formulate an estimator for the population-average survival functions

[TABLE]

where $\hat{q}^{(i)}_{z}$ and $\hat{S}^{(i)}_{z,t}$ are estimators of the cohort prevalence and cohort-wise survival respectively. This weighted Kaplan Meier method has appeared previously in the literature Murray [2001], and has been empirically validated against the pure Kaplan Meier method Zare et al. [2014], where the weighting procedure was found to reduce the bias in the construction of survival curves. The asymptotic convergence of the product-limit estimator and weighted variants is well established Cai [1998], Pepe and Fleming [1991]. We use this survival curve reconstruction method as a base in constructing a new statistic for comparing populations. The focus of this manuscript is not the properties of this survival estimator but rather the asymptotic convergence of its bounding area and the use of such a quantity for evaluating a null hypothesis.

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 individuals within each cohort, $n_{z}^{(i)},$ is counted, from which an estimator of the cohort distribution is obtained,

[TABLE]

In turn, we assume that the cohort-level survival functions $\hat{S}^{(i)}_{z,t}$ are estimated independently using the product-limit estimator. Note that since the product limit estimator is not a linear functional of sampled lifetimes, $\hat{\theta}_{t}^{(i)}$ 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 $\hat{S}$ and all instances of weighted sums of product limit estimators using the Greek letter $\hat{\theta}.$

With these elements in place, we define our test statistic

[TABLE]

where $\tau=\inf\left\{\tau_{z}:z\in\mathbb{Z}_{d}\right\}$, and $\tau_{z}$ denotes the time at which cohort $z$ is censored in observations. Note that in the absence of censoring this statistic is equivalent to comparison of mean lifetimes between the two populations Pepe and Fleming [1989]. We state here the main result of the paper – the large sample behavior of this statistic within a null-hypothesis statistical testing framework.

Theorem 1**.**

Let $C^{(i)}_{z,t}$ denote the probability that a $z$-type individual has not yet been censored at time $t\geq 0$ (the survival probability relative to the occurrence of censoring), and $q^{(i)}_{z}$ denote the probability that an individual in population $i$ is of cohort $z$, and let $p^{(i)}=n^{(i)}/(n^{(1)}+n^{(2)}).$ Suppose that $\theta^{(1)}_{t}=\theta^{(2)}_{t}.$ Then $\hat{\Theta}\xrightarrow{d}N(0,\sigma^{2})$, as $n^{(i)}\to\infty$, with

[TABLE]

where for $0\leq t\wedge\tau_{z}$, where $\tau_{z}$ is the time at which samples of cohort $z$ are censored, $\phi_{z,t}=\int_{t}^{\tau_{z}}\!\mathrm{d}s\,S_{z,s}$, $\phi_{z}\equiv\phi_{z,0}$, $S_{z,t}$ is the survival function for the pooled data of cohort $z$, and

[TABLE]

Note that this quantity is well-defined since by definition of $\tau_{z}$, ${C}^{(z)}_{z,t}>0$ for all $t\leq\tau_{z}$. The variance $\sigma^{2}$ may be consistently estimated by

[TABLE]

where for $0\leq t\wedge\tau_{z}$, $\hat{S}_{z,t}$ is the product-limit estimate of the pooled data for cohort $z$,

[TABLE]

$\hat{C}^{(i)}_{z,t}$* is the product-limit estimate associated to the event of censoring for cohort $z$ within population $i$, $\hat{\phi}_{z}\equiv\hat{\phi}_{z,0},$ and*

[TABLE]

Note that this quantity is also well-defined since $\hat{C}^{(z)}_{z,t}>0$ for all $t\leq\tau_{z}$. In Appendix A, we provide a proof of Theorem 1 in an empirical process framework. Note that since survival estimates $\hat{\theta}$ and $\hat{S}$ are step functions, all integrals are exactly computable.

3 Numerical investigation

A computational implementation of the test statistic $\hat{\Theta}$ and weighted survival estimators is available in the form of a package for R. This package also contains a class to handle arithmetic involving right-continuous piecewise linear functions. In the appendices we have provided source code that may be used for installing and invoking this package.

Here, we present a computational investigation of the weighted survival curve estimator and the corresponding test statistic. Using simulations, we investigated the statistical power of $\hat{\Theta}$, contrasted with that of existing non-parametric methods. Using a real dataset, we demonstrate the computation of $\hat{\Theta}$, $\hat{\theta}_{t}$, and evaluate Type-I error.

3.1 Evaluating statistical power through simulations

Using simulations, we explored the statistical power of the test statistic $\hat{\Theta}$ in a case where populations are difficult to distinguish based purely on mean survival time.

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 $\lambda^{-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 $\lambda=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 $q_{2}$. 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 $q_{2}\approx 0.245$. For this reason, we chose values near $0.25$ for our simulations.

To compare the reweighted Kaplan-Meier estimator (Eq. 2.2) to the standard Kaplan-Meier estimator, we estimated survival for the admixed population for $q_{2}=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 the 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 $q_{2}.$ For samples of size $n^{(i)}\in\left\{30,50,100,200,1000\right\}$ taken from each population, we performed $1000$ null hypothesis statistical tests using our method, the log-rank method Berty et al. [2010], and the standard Kaplan-Meier Wilcoxon signed-rank difference-of-mean methods Wilcoxon [1945], Schoenfeld [1981]. The power of the test, or the proportion of times that the null hypothesis was correctly rejected, is shown in Fig. 4.

3.2 Evaluating Type-I error in a real world example

We applied the survival estimator and statistic to NCCTG Lung Cancer data Loprinzi et al. [1994] available within the survival package for R. We compared the survival between male ($n^{(1)}=136$) and female ($n^{(2)}=90$) cancer patients, organized by ECOG performance score ($z\in\{0,1,2\})$ as cohort. Using males as population 1 and females as population 2, we arrived at the test-statistic estimate: $\hat{\Theta}=-961$, with 95% asymptotic confidence interval: $(-1527,-396)$, which would support rejection ($P\approx 0.0009$) of the null hypothesis ($\hat{\theta}_{t}^{(1)}=\hat{\theta}_{t}^{(2)}$) at $\alpha=0.05$. For reference, both the Wilcoxon ($P\approx 0.0012$) and log-rank ($P\approx 0.0015$) tests referenced in Fig. 4 also rejected the null hypothesis.

In theory, the Type-I error is set by the significance level at study design. Whether a statistic controls Type-I error correctly depends on accurate evaluation of its sampling distribution. In the case of $\hat{\Theta}$, our main result is that the sampling distribution for this estimator converges asymptotically in distribution to a Gaussian with a definite variance. However, small-sample behavior is not guaranteed. To evaluate Type-I error, we used the same dataset, restricted to male patients. For each of $n\in\{40,80,136\}$, we sampled without replacement the $n$ male patients split into two groups so that $n^{(1)}=n^{(2)}=n/2$, and compared survival between the two random groups. Repeating this procedure $10,000$ times, we generated the observed distribution of $P$-values, presented in Fig. 6 in log-scale. The distributions computed using the three methods are similar. The three methods all rejected $H_{0}$ approximately $5\%$ of the time except for the case of $\hat{\Theta}$ at $n=40$, which rejected $H_{0}$ approximately $6\%$ of the time. Essentially, asymptotic convergence as defined by the accurate evaluation of $\alpha=0.05$ Type-I error occurs somewhere in between $40$ and $80$ samples for this particular dataset.

Probing deeper, we examined the sampling distributions of $\hat{\Theta}$ for each of $n\in\{50,60,70\}$, in each instance compared to the Gaussian distribution stated in Thm. 1, where the approximation is computed using only the first sample of size $n$. The results for these simulations are shown in Fig. 7, where it is seen that the sampling distribution of $\hat{\Theta}$ is approximately the same as the computed asymptotic Gaussian distribution, which is traced out in red.

The R code used to compute these examples is available in Appendix BB.3.

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 statistical 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 the survival curves generated from the two methods are distinct yet similar. In the second and third rows of Fig. 3, one sees that this reweighted estimator has comparable performance to the pure Kaplan-Meier estimator at large sample sizes. Asymptotically, both estimators converge to the true survival function, with variance converging to zero. At small sample sizes, there are differences. The reweighted estimator has reduced variance at the cost of larger bias, in a time-dependent manner. It also appears to have smaller variance at the cost of larger error at earlier times. This error at earlier times is mitigated by decreased error at later times (better reconstruction of tails), however, the estimator variance is lower at all times. Hence, dependent on costs, for small samples, this reweighted estimator may be preferable to the pure Kaplan-Meier estimator.

In simulations of the test statistic derived from the reweighted survival estimator, we saw superior performance compared to existing methods. In Fig. 4, it is seen that in all cases, the test statistic $\hat{\Theta}$ was better at distinguishing between the two populations than either the Wilcoxon signed-rank test or the log-rank test. The relatively-high statistical power of this statistic is due to tighter variation in the test-statistic. In nearly all cases $(>99.5\%)$, the estimator variance for the tested method was less than that of the other two tests (not shown).

This manuscript derives the asymptotic convergence in distribution of the $\hat{\Theta}$ statistic. Numerically, we demonstrated convergence of the statistic in Figs. 6 and 7, where we verified that the asymptotic approximation respects Type-I error at $\alpha=0.05$ and where we observe good match between the sampling distribution of $\hat{\Theta}$ and the asymptotic Gaussian distribution provided by Thm. 1.

A variant of this method was used in Rasch et al. [2014] 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 entry times and a heuristic “survival surface” was generated in order to get a single overall picture of the survivability of a given disorder. The censoring parameters $\tau_{z}$ varied due to the finite sampling window and the fact that more-recent cohorts are not observed for as long a time period as older cohorts, as depicted in Fig 1b. It was also expected that survival by cohort would vary due to differences in health care administration and treatment between entry cohorts. The use of the empirical prevalences ($\hat{q}_{z}$) allowed the accounting for variability in disability application volume by sufferers of given disorders, conditional on entry date.

We note that a strong limitation of the presented method lies in its framing in terms of null hypothesis statistical testing. The $\hat{\Theta}$ statistic only provides a $P$-value, as opposed to other tests such as the log-rank test which provide hazard ratios as well. As a trade-off for statistical power, one is sacrificing interpretability in the form of effect sizes.

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. [2014]. 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 construct 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 Chang et al. [2014] in order to infer true continuously-indexed surfaces.

\ethics

This section does not apply.

\dataccess

All data in this manuscript is simulated, with R source code provided in Appendix B.2.

\competing

The authors declare no competing interests.

\disclaimer

This section does not apply.

\aucontribute

A.H. and M.H. developed the statistical method. A.H., M.H, and J.C.C. wrote the proof. A.H., M.H, and J.C.C. performed the simulations. J.C.C generated the figures. A.H., M.H, and J.C.C. wrote the manuscript. All authors gave final approval for publication.

\funding

This work is supported by the Intramural Research Program of the National Institutes of Health Clinical Center and the US Social Security Administration. \ack The authors would like to thank Dr. Leighton Chan and Dr. Elizabeth Rasch for insightful discussions, guidance and support, Dr. Pei-Shu Ho for help obtaining data.

Appendix A Proof of the main theorem

To prove the main theorem, we use an empirical process modeling framework to develop the asymptotic properties of first deterministically proportionally-weighted Kaplan-Meier estimators. We then replace the deterministic proportions with estimates given by the sample prevalences of the cohorts. Here, we restate the main theorem and prove it through a series of lemmata.

Theorem 1**.**

Let $C^{(i)}_{z,t}$ denote the probability that a $z$-type individual has not yet been censored at time $t\geq 0$ (the survival probability relative to the occurrence of censoring), and $q^{(i)}_{z}$ denote the probability that an individual in population $i$ is of cohort $z$, and let $p^{(i)}=n^{(i)}/(n^{(1)}+n^{(2)}).$ Suppose that $\theta^{(1)}_{t}=\theta^{(2)}_{t}.$ Then $\hat{\Theta}\xrightarrow{d}N(0,\sigma^{2})$, as $n^{(i)}\to\infty$, with

[TABLE]

where for $0\leq t\wedge\tau_{z}$, where $\tau_{z}$ is the time at which samples of cohort $z$ are censored, $\phi_{z,t}=\int_{t}^{\tau_{z}}\!\mathrm{d}s\,S_{z,s}$, $\phi_{z}\equiv\phi_{z,0}$, $S_{z,t}$ is the survival function for the pooled data of cohort $z$, and

[TABLE]

The variance $\sigma^{2}$ may be consistently estimated by

[TABLE]

where for $0\leq t\wedge\tau_{z}$, $\hat{S}_{z,t}$ is the product-limit estimate of the pooled data for cohort $z$,

[TABLE]

$\hat{C}^{(i)}_{z,t}$* is the product-limit estimate associated to the event of censoring for cohort $z$ within population $i$, $\hat{\phi}_{z}\equiv\hat{\phi}_{z,0},$ and*

[TABLE]

Overview of Proof of Theorem 1.

To prove the main theorem, we turn to the modeling framework that we present in  A.2. In general, we proceed by first assuming fixed sample proportions and then extending results to random proportions as given by empirical prevalence (Eq 2.3). The convergence of $\hat{\Theta}$ follows directly from corollary 10 and Eq. A.20. The consistency of $\hat{\sigma}^{2}$ follows from theorem 4.2.2 of Gill [1980], which provides for weak convergence of the product limit estimator to a Gaussian process, and the Glivenko-Cantelli theorem.

∎

A.1 Preliminaries and Notation

Given any pair of random elements $X,Y$, we denote equality in a distributional sense by $X\approx Y$. Let $\mathbb{P}$ be a probability measure on the measurable space $\left(X,\mathcal{A}\right)$. The empirical measure generated by the sample of random elements $x_{1},\dots,x_{n}$, $n\in\mathbb{N}$ is given by

[TABLE]

where for any $x\in X$, and any $A\in\mathcal{A}$,

[TABLE]

Note that alternatively, when needed, one may write $\delta_{x}(A)$ as the indicator function $1_{A}(x)$ on the set $A$. Furthermore, in the case that $A=\left\{k\right\}$, $k\in\mathbb{Z}$, and $x\in\mathbb{Z}$, we write $\delta_{x}(A)\equiv\delta_{x,k}$.

Given $\mathcal{H}$, a class of measurable functions $h:X\to\mathbb{R}$, the empirical measure generates the map $\mathcal{H}\to\mathbb{R}$ given by $h\mapsto\mathbb{P}_{n}h$, where for any signed measure $Q$ and measurable function $h$, we use the notation $Qh=\int\!\mathrm{d}Q\,h$. Furthermore, define the $\mathcal{H}$-indexed empirical process $\mathcal{G}_{n}$ by

[TABLE]

and with the empirical process, identify the signed measure $\mathcal{G}_{n}=n^{-1/2}\sum_{i=1}^{n}\left(\delta_{x_{i}}-\mathbb{P}\right)$.

Note that for a measurable function $h$, from the law of large numbers and the central limit theorem, it follows that $\mathbb{P}_{n}h\displaystyle\xrightarrow{a.s.}\mathbb{P}h$, and $\mathcal{G}_{n}h\xrightarrow{d}N{\left(0,\mathbb{P}\left(h-\mathbb{P}h\right)^{2}\right)}$, provided $\mathbb{P}h$ exists and $\mathbb{P}h^{2}<\infty$, and where “$\xrightarrow{d}$” denotes convergence in distribution. In addition to the preceding notation, given the elements $f$, and $f_{n}$, $n\in\mathbb{N}$, we also denote respectively, convergence in probability and in distribution, of $f_{n}$ to $f$, by $f_{n}\xrightarrow{p}f$.

For any map $x:\mathcal{H}\to\mathbb{R}^{k}$, $k\in\mathbb{N}$, define the uniform norm $\left\|x\right\|_{\mathcal{H}}$ by

[TABLE]

and in the case that $\mathcal{H}\subset\mathbb{R}$, write $\left\|\cdot\right\|_{\mathcal{H}}\equiv\left\|\cdot\right\|_{\infty}$. A class $\mathcal{H}$ for which $\left\|\mathbb{P}_{n}-\mathbb{P}\right\|_{\mathcal{H}}\to 0$ is called a $\mathbb{P}$-Glivenko-Cantelli class. Denote by $\ell^{\infty}(\mathcal{H})$ the class of uniformly bounded functions on $\mathcal{H}$. That is, for a general $k\in\mathbb{N}$,

[TABLE]

If for some tight Borel measurable element $\mathcal{G}\in\ell^{\infty}(\mathcal{H})$, $\mathcal{G}_{n}\xrightarrow{d}\mathcal{G}$, in $\ell^{\infty}(\mathcal{H})$, we say that $\mathcal{H}$ is a $\mathbb{P}$-Donsker class.

A.2 Empirical process framework

To prove Theorem 1, we turn to an empirical modeling framework that will provide us the asymptotic statistics of the weighted product limit estimator. Consider a closed particle-system, such that according to a predefined set of characteristics, the system can be subdivided into mutually exclusive subsystems.

Each particle corresponds to the observed state of a particular individual in a fixed population cohort. Note that we will restrict this discussion to only a single population of particles. These arguments will extend to multiple populations as mentioned in this manuscript by treating separate populations as independent.

At any given time $t\geq 0$, each particle will have exactly one associated state $x$ in the set $\mathbb{Z}_{4}$, referring respectively to states of

[TABLE]

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=1),$ and that particles are considered dormant for all $t<0$.

Let $d\in\mathbb{N}$ and $\tau\in(0,\infty)$ be fixed. We will assume the existence of a collection of individuals $\Gamma$, assumed to be infinite in size, where each individual $\gamma\in\Gamma$ exhibits a càdlàg path-valued state $x_{t}^{\gamma}$, for $t\geq 0$. For each $\gamma\in\Gamma$, $x^{\gamma}_{t}$ is determined by the individuals particle type $z^{\gamma}$ and a random jump time $\xi^{\gamma}$. The particle type $z^{\gamma}$ is distributed in the population through the probability mass $\mathbb{P}(z^{\gamma}=z)=q_{z}$, where $\mathbf{q}=(q_{1},\dots,q_{d})\in(0,1)^{d}$ satisfies $\sum_{z=1}^{d}q_{z}=1$. Let $S_{t}=(S_{1,t},\dots,S_{d,t})$ be the survival vector $S_{z,t}=\mathbb{P}\left\{T_{z}>t\right\}$, which is assumed continuous for $t\geq 0$. Suppose that it is desired to understand the event probabilities for randomly selected $\gamma\in\Gamma$, unconditional on subgroup membership. We assume that members of each cohort are in the inactive (0) state at times $t<0$.

Given a random sample $\gamma_{1},\dots,\gamma_{n}$, $n\in\mathbb{N}$ of individuals, let $\mathbf{n}=(n_{1},\dots,n_{d})$ and

[TABLE]

where $n_{z}$ is the random number of drawn individuals of cohort $z$. In considering the event time probabilities of each subgroup, the random number of 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 $n_{z}=a_{z}n$, $a_{z}\in(0,1)$, of individuals of cohort $z$, and let $\mu_{j,z,t}^{n_{z}}\geq 0$ be the number of the cohort $z$ individuals who are in state $j\in\mathbb{Z}_{4}$ at time $t$, so that

[TABLE]

is conserved. Also, we assume that there exists $\tau_{z}<\infty$ when all particles either become inactive or censored so that $\tau_{z}$ is the infimum time where the condition

[TABLE]

holds.

For the sample of size $n_{z}$, we denote the $z$-type cumulative hazard by $\Lambda_{z,t}$ and respectively define the $z$-type cumulative hazard and survival estimates by

[TABLE]

Define further

[TABLE]

and note that $\hat{S}^{n_{z}}_{z,t}=\hat{S}^{n_{z}}_{z,\tau_{z}}$ and ${B}^{n_{z}}_{z,t}={B}^{n_{z}}_{z,\tau_{z}}$ for all $t\geq\tau_{z}$.

From Gill [1980], it follows that $\left\{B^{n_{z}}_{z,t}:t\geq 0\right\}$ is a mean-zero square-integrable martingale with Meyer bracket process

[TABLE]

where $t\wedge\tau_{z}=\min\{t,\tau_{z}\},$ and $\delta_{(\cdot,\cdot)}$ is the Kroenicker delta function.

A.3 Convergence Theorems

In order to guarantee convergence of the estimator, we make the following assumptions (based upon an initially known sample size distribution $\mathbf{n}$).

Assumption 2**.**

We assume that the initial sample is chosen large enough to ensure that individuals of cohort $z$, at state ${1}$ (active), exist at all points $t\in[0,\tau_{z}]$, $z\in\left\{1,\dots,d\right\}$. That is,

[TABLE]

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

[TABLE]

for some constant $c>0$.

Assumption 3**.**

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

[TABLE]

Assumption 4**.**

For each $z\in\left\{1,\dots,d\right\}$ there exists a non-increasing continuous function $m_{z}:[0,\infty)\to(0,1]$ such that

[TABLE]

Note that in the case of fixed censoring, that is, in the case that censoring exists only at time $\tau$, the above is satisfied by $m_{z,t}=S_{z,t}$. In the general case, $m_{z,t}$ can be seen as the probability that an individual of cohort $z$ has not yet left state $1$. That is, $m_{z,t}$ is the probability that an individual has not left due to censoring or death by time $t$, and so $m_{z,t}=S_{z,t}C_{z,t-}$, where $C_{z,t}$ is the probability that censoring has not occurred by time $t$.

To prove the main theorem, we now present a series of lemmata.

Lemma 5**.**

If $\hat{q}$ is defined is in Eq. 2.3 and $\hat{S}_{z,s-}^{n_{z}}$ is defined as in Eq. A.15, then

[TABLE]

as $n\to\infty$, uniformly in $t\geq 0$.

Proof.

It is claimed that to prove the statement of the lemma, it suffices to show that

[TABLE]

uniformly in $t\geq 0$, for each $z=1,\dots,d$.

Indeed, for if the above holds, then

[TABLE]

uniformly in $t\geq 0$. Since the central limit theorem implies that $\sqrt{n}(\hat{q}_{z}-q_{z})\xrightarrow{d}N(0,q_{z}(1-q_{z}))$, each term in the sum would converge in probability to [math], uniformly in $t\geq 0$.

And so, if $\mathbb{E}_{N}$ denotes the expectation given $N$, we have that

[TABLE]

for some arbitrary constant $C$. From Lenglart’s inequality (cf. Lenglart [1977]),

[TABLE]

for any arbitrary $\eta,\epsilon>0$. Therefore, from Assumption 3, since $n_{z}\to\infty$ a.s., the desired result follows. ∎

Turning momentarily to the situation where there are two populations denoted by superscripts $(1)$, and $(2)$, for any $t\geq 0$, define

[TABLE]

noting that setting $\theta_{s}^{(1)}=\theta_{s}^{(2)}$ recovers our test statistic of Eq. 2.4. For a general survival function $\theta$, with respective estimate $\hat{\theta}$, define $\hat{Y}_{t}$ by

[TABLE]

If the process $\hat{Y}$ converges in distribution to some $Y\sim N(0,\sigma^{2})$, since $n^{(i)}/(n^{(1)}+n^{(2)})$ converges to $p^{(i)}$, $i=1,2$, it follows that

[TABLE]

Now we turn to analysis under a single population, dropping the superscripts. Note that $\hat{Y}_{t}=\sum_{z=1}^{d}\hat{Z}_{z,t}$, where

[TABLE]

Therefore, if it can be shown that

[TABLE]

uniformly in $t$, then convergence of $(\hat{Y}_{t}:t\geq 0)$ is dependent only upon the convergence of the $d$-dimensional vector-valued process $\hat{\zeta}(\hat{q})$ given by

[TABLE]

with $a=(a_{1},\dots,a_{d})\in(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 $\hat{\zeta}_{t}(\hat{q})$, we first focus on convergence of $\hat{\zeta}_{t}(a)$.

Let $\phi_{z,t}=\int_{t}^{\tau_{z}}\!\mathrm{d}s\,S_{z,s}$ and write $\hat{\zeta}_{t}(a)=\hat{\zeta}^{1}_{t}+\hat{\zeta}^{2}_{t}(a)$, where

[TABLE]

and

[TABLE]

Lemma 6**.**

Suppose that $\left\{\hat{\zeta}^{1}_{t}(a):t\geq 0\right\}$ and $\left\{\hat{\zeta}^{2}_{t}(a):t\geq 0\right\}$ are the processes respectively defined by equations (A.23) and (A.24), and that $\tilde{B}$ is the $d$-dimensional mean-zero Gaussian process defined by

[TABLE]

Then $\hat{\zeta}^{1}_{t}\xrightarrow{d}\zeta^{1}_{t}\mathrm{and}\quad\hat{\zeta}^{2}_{t}(a)\xrightarrow{d}\zeta^{2}_{t}(a),$ in the space of compactly supported functions $\mathcal{D}_{\mathbb{R}^{d}}[0,\infty)\,\,\mathrm{as}\,\,n\to\infty,$ for each $a\in V$, where $\zeta^{1}_{t}=(\zeta^{1}_{1,t},\dots,\zeta^{1}_{d,t})$ is the mean-zero square-integrable Gaussian process defined by

[TABLE]

and $\zeta^{2}_{t}(a)=(\zeta^{2}_{1,t}(a),\dots,\zeta^{2}_{d,t})$ is given by

[TABLE]

The processes $\hat{\zeta}^{1}$ and $\hat{\zeta}^{2}(a)$ are independent, and there exist a Skorohod representations such that

[TABLE]

and

[TABLE]

almost surely as $n\to\infty$.

Proof.

To begin note that independence follows immediately from the independence of the respective limiting processes. Since $\mathbf{n}$ is a multinomial random variable, (A.26) follows from the central limit theorem. In the case of $\hat{\zeta}^{2}_{t}(a)$, we first consider $B_{z,t}^{na_{z}}$.

An application of Lenglart’s inequality, very similar to that in the proof of Lemma 5, along with Assumption 3, shows that

[TABLE]

Moreover, from Assumption 4,

[TABLE]

It follows that

[TABLE]

uniformly in $t\geq 0$, and since $m_{z,t}=S_{z,t}C_{z,t-}$,

[TABLE]

Therefore, from theorem 4.2.1 of Gill [1980], $B^{na_{z}}_{z,t}\xrightarrow{d}\tilde{B}_{z,t}$, and there exists a Skorohod representation of $B^{na_{z}}_{z,t}$ such that

[TABLE]

almost surely as $n\to\infty$. Since almost sure convergence of $B_{z,t}^{na_{z}}$ implies almost sure convergence of bounded functionals of $B_{z,t}^{na_{z}}$, the desired convergence of $\hat{\zeta}^{2}(a)$ follows from Theorem 2.1 of Gill et al. [1983]. ∎

Corollary 7**.**

If the process $\hat{\zeta}(a)=\left\{\hat{\zeta}_{t}(a)\right\}$ is defined by equation (A.3), then

[TABLE]

Proof.

From the previous theorem we may assume that $\hat{\zeta}^{1}_{z,t}\to\zeta^{1}_{z,t}$ and $\hat{\zeta}^{2}_{z,t}(a)\to\zeta^{2}_{z,t}(a)$ almost surely, uniformly for $a\in V$ and $t\geq 0$. Therefore

[TABLE]

almost surely, uniformly for $a\in V$ and $t\geq 0$. The statement of the theorem then follows from theorem 5.1 of Billingsley [2013]. ∎

Since $\mathbf{n}/n\xrightarrow{p}\mathbf{q}$, from Theorem 4.4 of Billingsley [2013]

[TABLE]

Define the map $g:V\times\ell^{\infty}(V\times[0,\infty))\to\ell^{\infty}([0,\infty))$ by $g(a,f)=f(a,\cdot)$, then

[TABLE]

Furthermore, if for any $(a_{1},f_{1}),(a_{2},f_{2})\in V\times\ell^{\infty}(V\times[0,\infty))$ we have that

[TABLE]

for some $\delta>0$, then

[TABLE]

Therefore, $g$ is a continuous at any $(a,f)$ such that $f$ is continuous at $a$, uniformly in $t$. It thus follows from the continuous mapping theorem (cf. Van Der Vaart et al. [1996]) that if $a\mapsto\sum_{z=1}^{d}\zeta_{z,t}(a)$ is continuous, uniformly in $t$, then

[TABLE]

Lemma 8**.**

If $\left\{\zeta_{t}(a):t\geq 0\right\}$ is defined as in Corollary 7, then the map

[TABLE]

is continuous for $a\in V$, uniformly in $t\geq 0$.

Proof.

For any $a,b\in V$, it follows that

[TABLE]

Since $S_{\tau_{z}}>0$ for all $z$, from Doob’s martingale inequality (cf. Karatzas and Shreve [2012]),

[TABLE]

for some arbitrary constant $C$. For each $z\in\mathbb{N}_{d}$, since $a_{z}$ and $b_{z}$ are sufficiently close to $q_{z}\in(0,1)$, it follows that there exists some $\delta>0$ such that $a_{z}\wedge b_{z}>\delta$. Therefore,

[TABLE]

Combining the above two results gives

[TABLE]

and so, by Kolmogorov’s continuity criterion (cf. Karatzas and Shreve [2012]), the desired result follows. ∎

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

Theorem 9**.**

Let $\sum_{z=1}^{d}\zeta_{z,t}^{n}(\cdot)$ and $\sum_{z=1}^{d}\zeta_{z,t}(\cdot)$ be defined as in Corollary 7, then

[TABLE]

Corollary 10**.**

If $\hat{\zeta}=\sum_{z=1}^{d}\zeta_{z,\tau_{z}}(q)$, then

[TABLE]

where

[TABLE]

Proof.

Note that when $t=\tau_{z}$, we have

[TABLE]

which are independent and normally distributed, implying that $\hat{\zeta}$ is also normally distributed. Furthermore

[TABLE]

which after recombining the final terms, gives the desired result. ∎

Appendix B Computation

B.1 Installation of R package

The following code installs the R package from github sources

⬇

1*# install devtools if not already installed*

2install.packages("devtools")

3library(devtools)

4install_github("joshchang/calonesurv")

5library(calonesurv)

B.2 Simulation of data used in this manuscript

We simulated draws from the populations mentioned in the main text using the following R code:

⬇

1library(survival)

2library(calonesurv)

3

4p = 0.75

5n = 20

6samples = 10000

7tvals = seq(0,1/explambda*4,by=0.05)

8

9explambda = 1/4

10wshape = 1

11wscale = 5

12

13mix_results <- matrix(nrow = length(tvals), ncol = samples)

14pure_results <- matrix(nrow = length(tvals), ncol = samples)

15

16for(i in 1:samples){

17 n1 = sum(rbinom(n,1,p) )

18 n2 = n-n1

19 #while(n2<2){

20 # n1 = sum(rbinom(n,1,p) )

21 # n2 = n-n1

22 #}

23 t1 = sort(rexp(n1,explambda))

24 t2 = sort(rweibull(n2,wshape,wscale))

25

26 km1 = cadlag(survfit(Surv(t1)~1))

27 km3 = cadlag(survfit(Surv(c(t1,t2))~1))

28 pure <- cadlag(km3$**time**,**c**(1,km3$surv))

29

30 if(n2>0){

31 km2 = cadlag(survfit(Surv(t2)~1))

32 mix <- n1/nkm1 + n2/nkm2

33 mix_results[,i] <- mix(tvals)

34 }

35 else{

36 mix_results[,i] <- pure(tvals)

37 }

38

39 pure_results[,i] <- pure(tvals)

40}

B.3 Real-world example

⬇

1require(coin)

2require(survival)

3source("Theta.R")

4library(data.table)

5library(ggplot2)

6

7surv_data = with(subset(survival::lung,ph.ecog %in% 0:2),

8 data.frame(population = sex,

9 censor = as.numeric(status==1),

10 time = time, cohort = ph.ecog ))

11

12out = Theta_hat(surv_data)

13print(out)

14print(confint(out))

15print(pvalue.Theta_hat(out))

B.3.1 Simulations for examining the sampling distribution of $\hat{\Theta}$

⬇

1

2J = 10000

3N = c(40,80,136)

4

5

6typeone = data.table("method" = rep(0,J*length(N)3), "n" = rep(0,Jlength(N)3), "P" = rnorm(Jlength(N)*3))

7Thetavals = data.table("n" = rep(0,Jlength(N)), "Theta" = rep(0,Jlength(N)))

8

9*# Simulations for evaluating Type-I error*

10

11k = 0

12j = 1

13for(n in N){

14 i = 1

15 while(i< J*3){

16 indices = sample(1:136,n)

17 s_data = subset(surv_data,population==1)[indices,]

18 s_data$population[1:(n/2)] = 2

19

20 tryCatch({

21 out = Theta_hat(s_data)

22 p_Theta_hat = pvalue.Theta_hat(out)

23 p_lr = pvalue(logrank_test(Surv(time,1-censor)~as.factor(population), data = s_data))

24 fit = survfit(Surv(time,1-censor)~as.factor(population), data = s_data)

25 t1 = ten(fit)

26 invisible(comp(t1))

27 p_w = attr(t1,’lrt’)$pNorm[2]

28 Thetavals[j] = list( "n" =n , "Theta" = as.numeric(out))

29 typeone[J3k + i,] = list("method"=1,"n"=n,"P" = p_Theta_hat)

30 typeone[J3k + i+1,] = list("method"=2,"n"=n,"P" = p_w)

31 # logrank is 3

32 typeone[J3k + i+2,] = list("method"=3,"n"=n,"P" = p_lr)

33 i = i + 3

34 j = j + 1

35 },error = function(err) print(err))

36

37 }

38 k = k+1

39}

40pdf("fig6.pdf",family="CM Roman", width=8, height=4)

41typeone$method = **as**.**factor**(typeone$method)

42levels(typeone$method) = c("Our Method", "Wilcoxon signed rank","Log-rank")

43pp = ggplot(typeone, aes(factor(n),P))

44print(pp + geom_violin(aes(fill=factor(method)))+geom_hline(aes(yintercept=0.05),linetype=’dashed’)+ ylim(1e-4,1)+

45 scale_y_continuous(trans=’log2’,breaks = c(0.0005,0.005,0.05,0.5,1),limits = c(1e-4,1)) + theme_bw() + theme(legend.position="top") +

46 ylab("P-value") + xlab("Sample size") + guides(fill=guide_legend(title="Method: ")))

47dev.off()

48

49

50typeone$reject = typeone$P < 0.05

51print(ftable(reject ~ method + n, data = typeone))

52

53

54*# Simulations for evaluating asymptotic convergence in distribution*

55

56Thetavals_2 = data.table("n" = rep(0,Jlength(N)), "Theta" = rep(0,Jlength(N)))

57

58j = 1

59for(n in c(50,60,70)){

60 i = 1

61 while(i< J*3){

62 indices = sample(1:136,n)

63 s_data = subset(surv_data,population==1)[indices,]

64 s_data$population[1:(n/2)] = 2

65 tryCatch({

66 out = Theta_hat(s_data)

67 Thetavals_2[j] = list( "n" =n , "Theta" = as.numeric(out))

68 j = j + 1

69 i = i + 3

70 },error = function(err) print(err))

71

72 }

73}

74

75*# compute the theoretical Gaussian densities*

76

77n = 50

78indices = sample(1:136,n)

79s_data = subset(surv_data,population==1)[indices,]

80s_data$population[1:(n/2)] = 2

81out_50 = Theta_hat(s_data)

82

83n = 60

84indices = sample(1:136,n)

85s_data = subset(surv_data,population==1)[indices,]

86s_data$population[1:(n/2)] = 2

87out_60 = Theta_hat(s_data)

88

89n = 70

90indices = sample(1:136,n)

91s_data = subset(surv_data,population==1)[indices,]

92s_data$population[1:(n/2)] = 2

93out_70 = Theta_hat(s_data)

94

95grid <- with(dd, seq(-750,750, length = 200))

96normaldens =

97 data.frame(

98 Theta = rep(grid,3),

99 n = rep(c(50,60,70),each=length(grid)),

100 density = c(dnorm(grid, mean = 0, sd = sqrt(slot(out_50,’variance’))),

101 dnorm(grid, 0, sqrt(slot(out_60,’variance’))),

102 dnorm(grid, 0, sqrt(slot(out_70,’variance’))))

103 )

104

105

106pdf("fig7.pdf",family="CM Roman", width=8, height=3)

107pp2 = ggplot(Thetavals_2,aes(Theta)) + geom_histogram(aes(x = Theta, y = ..density..),bins = 80,fill = "grey") + facet_wrap(~n) +

108 theme_bw() + xlab(expression(hat(Theta)))+

109 geom_line(aes(y = density), data = normaldens, colour = "red", size=1.1)

110print(pp2)

111dev.off()