Non-Parametric Estimation

Non-parametric survival analysis is the attempt to capture the distribution of survival data without making any assumptions about the shape of the distribution. That is, non-parametric analysis, unlike parametric analysis, does not assume that the survival data was Weibull distributed or that it was Normally distributed etc. Concretely, non-parametric estimation does not attempt to estimate the parameters of a distribution, therefore “non-parametric.” Parametric analysis is covered in more detail in the Parametric Estimation section but it is important to contrast non-parametric estimation against what it is not. So what exactly is non-parametric analysis?

Survival analysis is using statistics to answer the question ‘what is the probability that the thing survived to a particular time?’ Non-parametric analysis answers this by estimating the probability from the proportion failed up to a given time. This can be done by either estimating the probability of surviving a particular segment (the Kaplan-Meier approach) or by estimating the hazard rate and accumulating it (the Nelson-Aalen and Fleming-Harrington approach). When the data are too incomplete for either — left or interval censoring, or right truncation — the Turnbull estimator finds the most likely curve directly.

The price of making no assumption about shape is that a non-parametric estimate can only say something about times at which data were actually seen: it is a step function that changes only at observed values, it cannot extrapolate beyond the largest observation, and it cannot say anything below the earliest time at which items were being watched. Parametric models trade an assumption for the ability to interpolate and extrapolate; the non-parametric estimate is the yardstick against which that assumption is judged.

This page is the theory. For runnable examples of everything described here see Non-Parametric SurPyval Modelling.

The xrd format. Non-parametric estimation is well understood by appreciating the data format used to estimate the CDF. Specifically, the ‘xrd’ format and particularly understanding the r and d sets of that format.

The number of components at risk, \(r\), at a given time, \(x\), is the number of things at risk just prior to time \(x\). The number of deaths, \(d\), is the number of the at risk items that died (or failed) at time \(x\). So for completely observed data the number at risk counts down for every death. So r would count down, e.g. 6, 5, 4, 3… for each death, 1, 1, 1, 1, … So in this example there were 6 items at risk and one death at the first time. Then, because there was 1 death at the first time the number of items at risk has decreased to 5, therefore for the next death there are only 5 at risk. This continues further until there are no more items at risk because they have all died, i.e. there is 1 at risk and 1 death.

This can be extended to more than one death. For example, the risk set could be 8, 6, 5, 3, 2, 1. with an accompanying death set of 2, 1, 2, 1, 1, 1. In this example there were times where there were 2 deaths and therefore the number at risk decreased by 2 after that number of deaths.

So a complete example of this format is:

x = [1, 2, 3, 4, 5, 6]
r = [7, 5, 4, 3, 2, 1]
d = [2, 1, 1, 1, 1, 1]

This format for data is not how survival data is usually provided in text books or papers. Survival data is usually displayed with the simple list of failure times such as “1, 3, 6, 7, 10, 16”. The first step surpyval does for non-parametric analysis is to transform data into the xrd format. All the fit() methods for surpyval take as input the xcnt format, see more at the Types of Data docs. So if you provide surpyval with the data “1, 2, 3, 4, 5, 6” it will assume that each of them are one death, and then create the risk set from the death counts resulting in the xrd format from above. (If you already have data in the xrd format you can skip the conversion: every estimator except Turnbull has a from_xrd(x, r, d) method.)

Right censoring enters through the risk set. A right censored item is one we stopped watching while it was still working. It contributes no death, but it was at risk up to the time it was censored, so it is counted in \(r\) at every time up to and including its censoring time, and removed afterwards. In particular, an item censored at exactly the same time as a death is counted as at risk for that death: the convention is that failures at a given time happen just before censorings at that time. For example the data 1, 2, 2+, 3 (where + marks a censored value) becomes \(x = (1, 2, 3)\), \(r = (4, 3, 1)\) and \(d = (1, 1, 1)\). The censored item is in the risk set at time 2, and gone by time 3. This is the only thing censoring changes, which is why the Kaplan-Meier, Nelson-Aalen and Fleming-Harrington estimators handle right censoring with no modification at all.

Left truncation (delayed entry) also enters through the risk set. A left truncated item is one that only came under observation at some entry time \(t_l\); had it failed before then we would never have known it existed. Such an item cannot be counted as at risk before it entered. SurPyval uses the standard \((t_l, x]\) convention (the same one as R’s survival package and lifelines): an item is at risk at time \(t\) if \(t_l < t \leq x\). An item entering at exactly the time of an event is therefore not at risk for that event, and an observed value equal to its own entry time is rejected as invalid (it would have a zero-length observation window). For example, with values 2, 3, 3, 4, 5, 6 and entry times 0, 0, 1, 1, 2, 2 the risk set is \(r = (4, 5, 3, 2, 1)\) at \(x = (2, 3, 4, 5, 6)\): the two items entering at time 2 are not at risk for the failure at 2, but they are at risk by time 3, so the risk set grows from 4 to 5. With delayed entry the risk set need not decrease, which is the whole mechanism by which truncation is handled.

Two consequences are worth stating. First, the estimate is really of survival conditional on surviving to the earliest entry time — nothing can be learned about what happened before anyone was being watched. Second, if only a handful of items have entered at early times the risk set there is small, and a single early failure can move the estimate a long way (with one item at risk and one failure, the Kaplan-Meier estimate drops to zero and stays there). Always look at \(r\) when working with delayed entry.

What does not fit into the xrd format. Left censoring (we only know the failure was before some time), interval censoring (it was between two times) and right truncation (items that fail after some time are never seen) cannot be expressed as a count of deaths at a time with a known risk set: the death times are not known, or the number at risk is not known. SurPyval raises an error if you pass such data to the Kaplan-Meier, Nelson-Aalen or Fleming-Harrington estimators. The Turnbull estimator, described below, handles them by estimating the \(r\) and \(d\) sets rather than counting them.

Given we now understand the format of the data we can estimate the probability of survival to some time with non-parametric methods. The first method we will visit is the Kaplan-Meier.

Kaplan-Meier Estimation

Kaplan-Meier [KM] is a very popular method for estimating survival curves for populations. The insight for this method is that for each time there is a death, we can estimate the probability of having survived since the previous deaths. Using the data from above as an example, at time 1, there are 7 items at risk and there are 2 deaths. We can therefore say that the probability of surviving this period was (7 - 2)/7, i.e. 5/7. Then the next time there is a death, the probability of having survived that extra time is (5 - 1)/5, i.e. 4/5.

To be clear, this is the chance of survival between each death. Therefore the chance of surviving up to a given time is the chance of surviving each segment. Therefore the probability of surviving up to any given time is the probability of surviving through all the previous segments. The probability of surviving multiple outcomes is the multiplication of each of the survival probabilities. Surviving through three sections is equal to the probability that I survive the first, then multiply this by the probability of surviving the second, then multiplying this result with the probability of surviving the third. So continuing our example from above, the probability of surviving the first two segments is (5/7) x (4/5) = 4/7.

Therefore using the at risk count, r, and the death count, d, can be used to estimate the segment survival probabilities and the survival probability to any point can be found by multiplying these probabilities. Formally, this has the following formula:

\[R(x) = \prod_{i:x_{i} \leq x}^{} \left ( 1 - \frac{d_{i} }{r_{i}} \right )\]

where \(x_i\) are the distinct observed times, \(r_i\) the number at risk just before \(x_i\) and \(d_i\) the number of failures at \(x_i\). This is why it is also called the product-limit estimator. Because each factor is a ratio of counts the estimate is a step function that drops only at failure times; censored values change later factors (through \(r\)) but never cause a drop themselves. For the example above the estimate is 5/7, 4/7, 3/7, 2/7, 1/7 and finally 0: with no censoring the Kaplan-Meier is exactly one minus the empirical CDF.

Why this is the “right” answer. The Kaplan-Meier is the non-parametric maximum likelihood estimator (NPMLE) for right censored and left truncated data: of all distributions, the one putting probability mass \(d_i / r_i \times R(x_{i-1})\) at each failure time is the one that makes the observed data most likely. An equivalent and very intuitive construction is Efron’s redistribute-to-the-right algorithm [Efron1967np]: start with mass \(1/N\) on every observation, then, working from left to right, take the mass of each censored observation and share it equally among all observations to its right. A censored item has, after all, failed at some later time, and with no other information it is equally likely to be any of the later ones. The mass left on the failures is exactly the Kaplan-Meier. Keep this picture in mind: it is the same idea the Turnbull estimator generalises.

Greenwood’s variance

The Kaplan-Meier is an estimate, so it has uncertainty. Each factor \(1 - d_i/r_i\) is an estimated binomial proportion, and the uncertainty of the product is most easily found on the log scale, where the product becomes a sum. Working with the cumulative hazard \(H(x) = -\ln R(x)\), Greenwood’s formula [Greenwood1926np] is

\[\widehat{Var}\left(\hat{H}(x)\right) = \sum_{i:x_{i} \leq x} \frac{d_{i}}{r_{i} \left ( r_{i} - d_{i} \right )}\]

and, by the delta method, \(\widehat{Var}(\hat{R}(x)) \approx \hat{R}(x)^2 \, \widehat{Var}(\hat{H}(x))\). Two things are visible in the formula. The terms grow as the risk set shrinks, so the estimate is least certain in the right-hand tail where few items remain. And when \(d_i = r_i\) (everyone still at risk fails, which happens at the last value of a data set with no right censoring) the term is undefined: the estimate has reached 0 and Greenwood’s formula has nothing to say there. SurPyval stores the cumulative variance of \(\hat{H}\) for every fitted model in the attribute greenwood (the name is historical; for the Nelson-Aalen and Fleming-Harrington estimators it holds their own variance, given below).

From a variance to confidence bounds

All the estimators on this page produce confidence bounds the same way, from \(\hat{\sigma}^2(x) = \widehat{Var}(\hat{H}(x))\). Write \(z\) for the standard normal quantile: \(z = \Phi^{-1}(1 - \alpha/2)\) for a two-sided interval and \(z = \Phi^{-1}(1 - \alpha)\) for a one-sided bound, where \(\alpha\) is alpha_ci (0.05 by default).

The naive (“normal”, or plain Greenwood) interval is \(\hat{R} \pm z \hat{R} \hat{\sigma}\). It is symmetric, but a survival probability is not: near 0 or 1 this interval spills outside \([0, 1]\) and it tends to undercover in small samples.

The default (“exp”) interval instead applies the normal approximation to \(\ln \hat{H}(x) = \ln(-\ln \hat{R}(x))\), which is unbounded in both directions and much closer to normally distributed. Its standard error is \(\hat{\sigma}/\hat{H}\) (delta method again), and transforming back gives

\[\left[ \hat{R}(x)^{\exp\left(z\hat{\sigma}(x)/\hat{H}(x)\right)},\; \hat{R}(x)^{\exp\left(-z\hat{\sigma}(x)/\hat{H}(x)\right)} \right]\]

which always lies within \([0, 1]\). This is the “log(-log)” or “exponential Greenwood” interval, and it is what cb() and plot() use unless told otherwise (bound_type='normal' selects the plain interval). Bounds on \(F = 1 - R\) and on \(H = -\ln R\) are the corresponding transformations of the bounds on \(R\). A few practical details of the implementation:

  • The bounds are asymptotic, so only the normal (dist='z') statistic is offered; for small samples, or for the Turnbull estimator, use the bootstrap (below).

  • Where the variance is undefined because the estimate has reached zero (Kaplan-Meier with no right censoring at the last value) the lower bound is set to 0 and the upper bound to the last finite upper bound, so bounds can still be drawn to the last observation.

  • Where the variance is still zero (before the first failure, e.g. when the first values are right censored) the 'exp' bounds are \([1, 1]\): the estimate is 1 and the formula has no uncertainty to report.

  • Bounds are only reported within the range of the data; below the first or above the last observed value cb() returns nan.

Pointwise bounds, bands and the bootstrap

The interval above is pointwise: at any single time \(x\) it covers the true \(R(x)\) with probability \(1-\alpha\). It does not mean that the whole true curve lies between the two bound curves with that probability; a curve has many opportunities to escape somewhere. If the question is “is this whole curve (for example a fitted Weibull) consistent with the data?”, you need a simultaneous confidence band, which is wider. SurPyval provides the two classical bands via band(): the Hall-Wellner band [HallWellner1980np], whose width is proportional to \((1 + N\hat{\sigma}^2(x))/\sqrt{N}\) with \(N\) the total number of items, and Nair’s equal-precision band [Nair1984np], which is the pointwise interval with a larger critical value. The critical values come from the supremum of a (standardised) Brownian bridge over the range of the data, computed numerically (to about 1e-8, rather than read from a table or simulated) so that they are exact for whatever range the data give; the bands are defined only between the first and last points with a positive, finite variance. Near either end of the data the estimate rests on a handful of failures, or of items at risk, and its error is far from the normal the bands assume. Two things then matter. The first is the scale the band is applied on: SurPyval applies it, by default, to \(\arcsin\sqrt{\hat{R}(x)}\), whose half width is \(\tfrac12 c\,\hat{\sigma}(x)\sqrt{\hat{R}(x)/(1 - \hat{R}(x))}\) for a band of half width \(c\,\hat{\sigma}(x)\) on \(\hat{R}\) itself (the delta method), cut at 0 and \(\pi/2\) [BorganLiestol1990np]. The second is the range: the equal-precision band’s boundary, \(c\sqrt{a(1-a)}\) for the Brownian bridge in \(a = N\hat{\sigma}^2/(1 + N\hat{\sigma}^2)\), makes its critical value grow without limit as the range reaches \(a = 0\) or 1, and the theory is for a range \([t_L, t_U]\) inside the data. By default it covers the times with \(0.1 \le a \le 0.9\) (Klein and Moeschberger tabulate it for \(a_L\) from 0.02 and \(a_U\) to 0.98), and x_range=(t_L, t_U) sets any range. In simulation (n = 40 to 400) it then covers 94% to 96% for a nominal 95%; over the first to the last event it covered about 93% on the arcsine scale, 87% to 89% on the log(-log) scale (bound_type='exp', the pointwise default) and 83% untransformed ('normal'), most misses at the first few events. The Hall-Wellner band, whose boundary in \(a\) is flat, covers about 95% over the first to the last event. See [KleinMoeschberger2003np] (section 4.4) for the derivations.

The asymptotic formulas all rest on large-sample theory for counts that were actually observed. When that is doubtful (small samples, or Turnbull estimates built from expected counts) the non-parametric bootstrap is the robust alternative: bootstrap_cb() resamples the original observations with replacement, refits the same estimator to each resample and reports the percentile interval of the refitted curves at each requested time.

Quantiles, the median and the mean

Because the estimate is a step function, its quantile is defined as the smallest observed value at which the estimated CDF reaches the requested probability, \(\hat{q}(p) = \min\{x_i : \hat{F}(x_i) \geq p\}\), and the median is \(\hat{q}(0.5)\). SurPyval allows \(10^{-9}\) of round-off in that comparison, so that a curve which reaches \(p\) exactly on paper, such as the Kaplan-Meier of 1 to 30 at 15, is not put a step late. With right censoring the estimate may never reach \(p\) (the curve stops above \(1-p\)), in which case the quantile is undefined and SurPyval returns nan rather than guess. A confidence interval for a quantile is found by the Brookmeyer-Crowley method [BrookmeyerCrowley1982np]: invert the pointwise bounds, i.e. take all times at which the confidence interval for \(R(x)\) contains \(1 - p\).

In SurPyval quantile_cb(p) returns, for each \(p\), the first observed time at which the lower bound of \(R\) has fallen to \(1 - p\) (the lower end of the interval) and the first observed time at which the upper bound of \(R\) has fallen below \(1 - p\) (the upper end). With heavy censoring the upper bound may never get that low, and the upper end is then nan: the interval is open to the right, the data being consistent with the quantile lying beyond the last observation.

The mean is the area under the survival curve. With right censoring the curve does not reach zero, so the area to infinity is unknown; what can be estimated is the area up to a horizon \(\tau\), the restricted mean (see Restricted mean survival time below). SurPyval’s mean() integrates the step function from 0 to \(\tau\), defaulting \(\tau\) to the largest observed value; if the estimate reaches zero this is the ordinary mean of the estimated distribution. A \(\tau\) beyond the last observation holds the curve at its final value out to \(\tau\), which is an extrapolation, not an estimate.

No failures at all: success-run testing

A common reliability test is to run \(n\) items for a fixed duration (or demand) and hope that none fail. If none do, the Kaplan-Meier estimate is 1 with zero variance, which is true of the sample but useless as a statement about the population. The question is instead: what is the lowest reliability consistent with seeing \(n\) successes in a row? If the true probability of success is \(R\) and trials are independent, the chance of \(n\) successes is \(R^n\). The smallest \(R\) for which that chance is still at least \(\alpha\) is the lower \(1-\alpha\) confidence bound

\[R_L = \alpha^{1/n}.\]

For example, 59 successes demonstrate \(R \geq 0.95\) at 95% confidence (\(0.05^{1/59} \approx 0.9505\)), and 22 successes demonstrate \(R \geq 0.90\) at 90% confidence. This is surpyval.success_run(n, confidence=...) (or alpha=...).

Nelson-Aalen Estimation

The Nelson-Aalen estimator [NA], instead of finding the probability, estimates the cumulative hazard function, and given that we know the relationship between the cumulative hazard function and the reliability function, the Nelson-Aalen cumulative hazard estimate can be converted to a survival curve.

A word on names, because the literature is not consistent. The survival curve \(e^{-\hat{H}}\) built from the Nelson-Aalen cumulative hazard is called the Breslow estimator by some authors and, confusingly, the Fleming-Harrington estimator by others (R’s survfit has used type="fleming-harrington" for it, and type="fh2" for the tie-corrected version). In SurPyval, NelsonAalen is \(e^{-\hat{H}}\) with the plain Nelson-Aalen \(\hat{H}\), and FlemingHarrington is \(e^{-\hat{H}}\) with the tie-corrected \(\hat{H}\) described in Fleming-Harrington Estimation below. The two are identical when there are no tied failures.

The first step in computing the NA estimate is to convert your data to the x, r, d format. Once in this format the instantaneous hazard rate is found by:

\[h(x) = \frac{d_{x} }{r_{x}}\]

This estimate of the instantaneous hazard rate is the proportion of deaths/failures at a value, x, among those at risk there. Strictly, \(d_x/r_x\) is the increment of the cumulative hazard at \(x\) (a probability of failing at \(x\) given survival to just before it), not a rate per unit time. Then to find the cumulative hazard rate for any x we simply take the sum of the instantaneous hazard rates for all the values below x. Mathematically:

\[H(x) = \sum_{i:x_{i} \leq x}^{} \frac{d_{i} }{r_{i}}\]

Then, since we know that the reliability, or survival function, is related to the cumulative hazard function, we can easily compute it.

\[R(x) = e^{-H(x)}\]

So we now have the survival/reliability function. One benefit of the Nelson-Aalen estimator is that it does not estimate a probability of 0 for the highest value (in a completely observed data set). This means that for a completely observed data set the whole estimation can be plotted on a transformed y-axis. For this reason SurPyval uses the Nelson-Aalen as the default plotting position.

How it relates to the Kaplan-Meier. Compare the factor each estimator applies at a failure time: Kaplan-Meier multiplies by \(1 - d_i/r_i\), Nelson-Aalen by \(e^{-d_i/r_i}\). Since \(e^{-u} \geq 1 - u\) for all \(u\), the Nelson-Aalen survival estimate is always at least as large as the Kaplan-Meier. When \(d_i/r_i\) is small (large risk set, few ties) the two factors are almost equal and so are the estimates; they separate when a large fraction of the risk set fails at once — with many ties, or in the tail where the risk set is small. In the worked example above, the survival at time 1 is 0.714 by Kaplan-Meier and 0.751 by Nelson-Aalen, and at the last value 0 versus 0.077.

Variance. The variance of \(\hat{H}\) used by SurPyval for the Nelson-Aalen estimator is Aalen’s (Poisson-type) estimator [Aalen1978np],

\[\widehat{Var}\left(\hat{H}(x)\right) = \sum_{i:x_{i} \leq x} \frac{d_{i}}{r_{i}^{2}},\]

which Klein [Klein1991np] recommends for its small-sample performance. Unlike Greenwood’s formula it stays finite when \(d_i = r_i\), so Nelson-Aalen bounds are defined right up to the last observation. The bounds themselves are built exactly as described in From a variance to confidence bounds.

Hazard rates and densities from a step function

A non-parametric estimate of \(H\) is a step function, and the derivative of a step function is zero almost everywhere and infinite at the steps, so the hazard rate \(h(x) = dH/dx\) cannot be read off it directly. SurPyval offers two approximations. hf() differences the cumulative hazard between the points you ask for, so the answer is an increment over your grid, not a rate per unit time, and it depends on the grid you choose. (Two conventions of hf() are worth knowing: the first point has nothing before it to difference from and repeats the second point’s value, and a zero increment, i.e. a gap in the grid with no failure, is replaced by the previous non-zero one. A single point returns the jump of the step it falls in.) df() is the drop in the survival over the same step, the probability the estimate puts there (\(S(1 - e^{-\Delta H})\) for a step of \(\Delta H\) starting at \(S\), close to \(h e^{-H}\) for small steps), so it inherits the same grid dependence and conventions; it stays finite where a Kaplan-Meier estimate falls to zero, whose last hf() is infinite. The better estimate is smoothed_hf(), which spreads each jump \(\Delta \hat{H}(x_i)\) over a neighbourhood with an Epanechnikov kernel of bandwidth \(b\),

\[\hat{h}(t) = \frac{1}{b} \sum_{i} K\left( \frac{t - x_i}{b} \right) \Delta \hat{H}(x_i), \qquad K(u) = \tfrac{3}{4}(1 - u^2) \text{ for } |u| \leq 1,\]

with a correction near the ends of the data where part of the kernel falls outside the observed range (the sum is divided by the share of the kernel’s mass that lies inside the range). The bandwidth trades bias (large \(b\) flattens real features) against variance (small \(b\) is noisy); the default of one eighth of the observed range is only a starting point. The estimate is nan outside the observed range, and an infinite jump (a Kaplan-Meier falling to zero at the last failure) is left out of the sum. See [KleinMoeschberger2003np] (section 6.2).

Fleming-Harrington Estimation

The Fleming-Harrington estimator [FH], uses the same principle as the Nelson-Aalen estimator. That is, it finds the cumulative hazard function and then converts that to the reliability/survival estimate. However, the NA estimate assumes, for any given step, that the number of items at risk is the same for each of the tied deaths; the FH estimate changes this. If \(d\) items fail “at the same time” they did not really fail at the same instant — the ties are an artefact of recording to finite precision — so the FH estimator imagines them failing one after another, with the risk set shrinking by one after each. Mathematically, the hazard rate is calculated with:

\[h(x) = \frac{1}{r_{x}} + \frac{1}{r_{x} - 1} + \frac{1}{r_{x} - 2} + ... + \frac{1}{r_{x} - (d_{x} - 1)}\]

Which can be summarised as:

\[h(x) = \sum_{j=0}^{d_x-1} \frac{1}{r_x - j}\]

The cumulative hazard rate therefore becomes:

\[H(x) = \sum_{i:x_{i} \leq x} \sum_{j=0}^{d_i-1} \frac{1}{r_i - j}\]

and, as for the Nelson-Aalen, \(R(x) = e^{-H(x)}\). You can see that the cumulative hazard rate will be slightly higher than the NA estimate since:

\[\underbrace{\frac{1}{r_{x}} + ... + \frac{1}{r_{x}}}_{d_x \text{ terms}} \leq \frac{1}{r_{x}} + \frac{1}{r_{x} - 1} + ... + \frac{1}{r_{x} - (d_{x} - 1)}\]

with equality exactly when there is a single death/failure at that time (\(d_x = 1\)). So the Fleming-Harrington estimate is identical to the Nelson-Aalen when there are no tied failures, and differs from it only at tied failure times. There is also an ordering between all three estimators: since each term \(e^{-1/(r-j)} \geq 1 - 1/(r-j)\) and \(\prod_{j=0}^{d-1}(1 - \frac{1}{r-j}) = 1 - \frac{d}{r}\), at every step

\[R_{KM}(x) \leq R_{FH}(x) \leq R_{NA}(x).\]

In the worked example the survival at time 1 (two tied failures among seven) is 0.714, 0.734 and 0.751 respectively. The Fleming-Harrington and Nelson-Aalen estimates are particularly useful for small samples, see [FH].

The variance of \(\hat{H}\) uses the same tie-splitting as the estimator,

\[\widehat{Var}\left(\hat{H}(x)\right) = \sum_{i:x_{i} \leq x} \sum_{j=0}^{d_i-1} \frac{1}{(r_i - j)^2},\]

which is the variance used by R’s survfit with ctype=2, and reduces to the Nelson-Aalen variance \(\sum d_i/r_i^2\) when there are no ties.

The Turnbull estimator (below) produces expected, and therefore fractional, death counts. SurPyval extends the tie-splitting ladder to a fractional \(d\) by taking the whole terms \(1/r, 1/(r-1), \ldots\) for the first \(\lceil d \rceil - 1\) failures and a pro-rata share of the next term for the remainder, so the Fleming-Harrington estimator can be applied to a Turnbull ladder. For example \(d = 2.5\) gives \(1/r + 1/(r-1) + 0.5/(r-2)\), and any \(d \leq 1\) gives \(d/r\), the Nelson-Aalen term. This extension is SurPyval’s own choice, not something from [FH]. Counts that differ from a whole number only by floating-point round-off (the EM can produce \(1 + 2\times 10^{-16}\)) are treated as that whole number, so they do not add a spurious extra term with a near-zero risk set.

Turnbull Estimation

The Turnbull estimator is a remarkable non-parametric estimation method for data that can handle arbitrary censoring and truncation [TB]. The Turnbull estimator can be found with a procedure of finding the most likely survival curve from the data, for that reason it is also known as the Non-Parametric Maximum Likelihood Estimator (NPMLE). The Kaplan-Meier is also a non-parametric maximum likelihood estimator, so is there a contradiction? No: for data that the Kaplan-Meier can handle (exact, right censored and left truncated observations), the Turnbull NPMLE is the Kaplan-Meier. Turnbull’s contribution is the generalisation to data that the Kaplan-Meier cannot handle. (Whether SurPyval’s Turnbull output equals the Kaplan-Meier depends on an option; see Why a Turnbull fit need not equal the Kaplan-Meier.)

Every observation as an interval

The Turnbull estimate is really an estimate of the observed failures given censoring, and then the ‘ghost’ failures (as Turnbull describes it) due to truncation. Turnbull’s estimate converts all failures to interval failures regardless of the censoring. This is because a left censored point is equivalent to an intervally censored observation in the interval -Inf to x, and a right censored point is equivalent to an intervally censored observation in the interval x to Inf. An exactly observed failure at \(x\) is the degenerate interval \([x, x]\).

SurPyval uses the standard \((l, r]\) convention for censoring intervals: an observation censored in \((l, r]\) failed strictly after \(l\) and at or before \(r\). So an interval whose right end coincides with an exactly observed failure time can have failed at that time, and a left censored observation at \(x\), i.e. \((-\infty, x]\), can have failed at \(x\) itself. This matches Turnbull’s paper, and the \((t_l, x]\) convention used for truncation windows.

Next the time axis is cut into pieces at every distinct endpoint: every observed value, every interval end and every truncation time, with an extra zero-width piece \([x, x]\) at each exactly observed time to hold the mass of that failure. Let \(p_j\) be the (unknown) probability that a failure falls in piece \(j\); the survival curve is determined by these masses, so estimating the \(p_j\) is estimating the curve. Then for all the pieces between negative infinity and infinity we find how many failures happened in each. This value need not be a whole number since a single observation could have failed across several pieces.

The self-consistency (EM) algorithm

If we knew the curve we could say where each observation most likely failed; if we knew where each failed we could estimate the curve by counting. Turnbull’s algorithm alternates between these two, which makes it an instance of the Expectation-Maximisation (EM) algorithm. Suppose observation \(i\) (which may represent \(n_i\) identical items) has censoring interval \(A_i\) and truncation window \(B_i\), and write \(\alpha_{ij} = 1\) if piece \(j\) lies inside \(A_i\) (the observation could have failed there) and 0 otherwise, and \(\beta_{ij} = 1\) if piece \(j\) lies inside \(B_i\) (a failure there would have been observable) and 0 otherwise.

E-step, observed failures. Given the current masses \(p\), the expected number of observation \(i\)’s failures that fell in piece \(j\) is its count shared out over the pieces it could have failed in, in proportion to how likely each piece is:

\[\mu_{ij}(p) = n_i \frac{\alpha_{ij}p_j}{\sum_{k=1}^{m}\alpha_{ik}p_k}\]

where \(m\) is the number of pieces. For a right censored observation this is exactly the redistribute-to-the-right idea from the Kaplan-Meier section; for an exact observation all the count lands on its own zero-width piece.

E-step, ghosts. If an observation is truncated, it was only a possible observation among others that would have been seen had the observation not been limited. If the probability of falling inside the window is \(P(B_i) = \sum_k \beta_{ik}p_k\), then for every item seen there were, on average, \((1 - P(B_i))/P(B_i)\) unseen ‘ghost’ items whose failures fell outside it. Spreading them over the pieces outside the window gives:

\[\nu_{ij}(p) = n_i \frac{(1 - \beta_{ij})p_j}{\sum_{k=1}^{m}\beta_{ik}p_k}\]

M-step. We can then estimate the probability of failure in each piece as the total expected failures in that piece divided by the total expected failures:

\[p_j = \frac{\sum_{i=1}^{N} \left( \mu_{ij} + \nu_{ij} \right)}{M(p)} \qquad\text{where}\qquad M(p) = \sum_{i=1}^{N} \sum_{j=1}^{m} \left( \mu_{ij} + \nu_{ij} \right)\]

Using this estimate of the masses, it can be input to the start of this procedure and it done again. This can then be repeated over and over until the values do not change; a solution that reproduces itself this way is called self-consistent. At this point we have reached the NPMLE estimate of the survival function!

It helps to see the M-step in xrd terms. Call \(d_j = \sum_i (\mu_{ij} + \nu_{ij})\) the expected deaths in piece \(j\) and \(r_j = \sum_{k \geq j} d_k\) the expected number at risk just before it. The Kaplan-Meier on this expected ladder, \(\prod (1 - d_j/r_j)\), telescopes to a curve whose mass in piece \(j\) is exactly \(d_j / M\). The M-step is therefore just a Kaplan-Meier applied to expected counts, which is why a fitted Turnbull model has r and d attributes like every other estimator, with fractional values.

In SurPyval the iteration stops when the largest change in any \(p_j\) falls below tol (default \(10^{-10}\)), or after max_iter iterations (default 1000), with a warning if the tolerance was not reached. The iteration starts from equal masses on every piece (under truncation, on every piece some observation could have failed in). With the Kaplan-Meier option and no truncation it starts instead on Turnbull’s innermost intervals, the runs of pieces that begin where some observation’s possible failure pieces begin and end where some observation’s end, with no other beginning or end in between. Every other piece is dominated: any observation that could have failed there could also have failed in an innermost interval, so the NPMLE gives it no mass [TB]. Starting there spares the EM from draining that mass away, which it does only slowly, leaving residues of around \(10^{-9}\) that the confidence bounds are sensitive to (see below). Under left truncation the dominance argument fails in general (moving mass changes the windows’ probabilities), but with only exact and right censored observations the delayed-entry Kaplan-Meier is an NPMLE with its mass on the innermost intervals, so the iteration starts there too; that also makes it return the Kaplan-Meier where the data leave the maximum not unique (see below). EM is reliable but can be slow, especially when many observations are censored far to the right; raising max_iter is the first thing to try. Under truncation SurPyval also intersects each observation’s support with its own truncation window (an observed failure cannot have happened where it would not have been observed) and confines the mass to pieces that at least one observation could have failed in, which keeps the iteration away from meaningless solutions.

Why a Turnbull fit need not equal the Kaplan-Meier

Once the EM has produced the expected ladder \((r_j, d_j)\), SurPyval lets you choose how to turn it into a survival curve with turnbull_estimator: 'Kaplan-Meier', 'Nelson-Aalen' or 'Fleming-Harrington' (the default). Under truncation the EM always iterates with the Kaplan-Meier (self-consistency) update above and the chosen estimator is applied to the converged ladder; without truncation the chosen estimator’s curve is also used within the iteration, i.e. the masses \(p_j\) of the next E-step are the drops of the chosen estimator’s curve. The Nelson-Aalen and Fleming-Harrington curves never reach zero, so in that case the drops add up to less than one and a small share of the right censored items’ expected failures is placed beyond the last finite value (it shows up in \(r\) but not in the reported \(d\)). Either way:

  • Only turnbull_estimator='Kaplan-Meier' gives the NPMLE. With it, on exact, right censored and left truncated data, Turnbull agrees with KaplanMeier for both the curve and the confidence bounds, to within the EM’s tolerance (around \(10^{-9}\)): with censoring or truncation the Turnbull fit is iterated to convergence rather than computed in closed form, so it is not equal to the last digit.

  • The Nelson-Aalen and Fleming-Harrington options are \(e^{-H}\) constructions on the expected ladder. They are not trying to maximise the likelihood, and they give a (usually slightly) higher survival curve, by the ordering shown in the Fleming-Harrington section. They are offered for the same reasons those estimators are preferred elsewhere: they do not drop to zero at the last failure and they behave better in the far tail.

  • Because the ladder is an expected one — right censored items appear as fractional failures spread over later pieces, and truncation adds ghost failures — a Turnbull fit with the Fleming-Harrington option is not in general equal to FlemingHarrington.fit on the same data either, even where the latter can be used.

For example, with values 2, 3, 3, 4, 5, 6 and entry times 0, 0, 1, 1, 2, 2 the Turnbull estimate of survival at 2 is 0.750 with the Kaplan-Meier option (identical to KaplanMeier), 0.765 with Fleming-Harrington and 0.779 with Nelson-Aalen. None of these is a mistake; they are different estimators. When comparing a Turnbull fit with a Kaplan-Meier fit, compare like with like.

Uncertainty in a Turnbull estimate

The confidence bounds from cb() for a Turnbull model use the variance formula of the chosen estimator, but the choice of ladder matters:

  • With no truncation and only exact and right censored observations (where Turnbull reduces to the Kaplan-Meier) the variance is computed from the observed counts, i.e. the ordinary Greenwood ladder. The EM’s expected ladder would treat redistributed censored mass as extra observed failures and give intervals that are too narrow.

  • With truncation, the variance ladder also uses observed counts: exact failures count once at their time, censored items leave the risk set at their censoring time, and each item is only at risk inside its own truncation window. Ghost failures, which are needed to get the estimate right, are not data and are excluded. Only genuinely interval or left censored items, whose failure time is unknown, are spread over their possible pieces.

  • With interval or left censoring and no truncation, the expected ladder is used as though it were observed. This ignores the uncertainty in where the censored failures were allocated, so the bounds are an approximation.

In the last case, and whenever you want calibrated intervals for a Turnbull estimate, use bootstrap_cb(). The simultaneous bands of band() rest on theory for right censored data and should not be used with interval censored Turnbull estimates.

The variance at each value uses the same pieces as the estimate there: the expected failures in a piece \((x_k, x_{k+1}]\) enter both at \(x_{k+1}\), where the estimate drops, so where the estimate is still 1 the bounds are \([1, 1]\). Expected counts that are equal, or zero, up to rounding error are treated as exactly so: at a last value where every item still at risk fails the variance is undefined, just as for exact data. One sensitivity remains. The width of the log(-log) interval on its own scale is roughly \(1/\sqrt{D}\) for \(D\) expected failures so far, so where the ladder holds only a tiny fraction of a failure the bounds are \([0, 1]\) even though the estimate is 1 to eight decimal places. That happens where the EM leaves a small residue of mass on a piece the NPMLE gives none, which it can with the Nelson-Aalen and Fleming-Harrington options or under truncation (the Kaplan-Meier option without truncation avoids it, as described above). bound_type='normal' or bootstrap_cb() avoid it too.

What the data cannot tell you: non-identifiability

The NPMLE is a delicate object, and it is worth being clear about what it does not determine.

Mass inside an interval. The likelihood only depends on how much mass lies in each piece, not on where within the piece it lies. Within a piece where the estimate drops, every curve between the value at the left end and the value at the right end fits the data equally well. SurPyval’s step function places each piece’s drop at the piece’s right end, which is a convention, not an estimate. The fitted model keeps both ends of that range: R_upper[k] is the estimate at x[k] (before the drop in the piece \((x_k, x_{k+1}]\), and the same as R[k]) and R_lower[k] the estimate at x[k+1] (after it), so every curve lying between them fits equally well. Evaluate Turnbull curves and bounds at piece boundaries (model.bounds, which also includes \(\pm\infty\) and any truncation times) when precision matters.

Mass outside the observation windows. With truncation the estimate is conditional: with left truncation it describes survival given survival to the earliest entry time; with right truncation it describes the distribution given failure before the latest truncation time. There is no information below the smallest left truncation time or above the largest right truncation time, and in effect the estimate assumes that the first value, if left truncated, had a 100% chance of observation, and likewise for the last value if right truncated. Only a parametric model, by assuming a shape, can extrapolate into those regions.

A likelihood with no maximum. Under truncation each observation’s likelihood is a ratio, the mass of the pieces it could have failed in (its support) over the mass of its truncation window. The supremum of such a product can lie where some windows carry no mass at all, a limit the parameter space does not contain: the NPMLE then does not exist, the EM climbs towards the boundary without settling, and raising max_iter does not help. Or the likelihood can be flat in some direction, so the maximum exists but is not unique. Which applies is a property of the data, not of the fit, and SurPyval decides it before the EM runs; the answer is model.npmle, one of "exists", "not unique", "does not exist" or "undetermined".

The criterion (derived in full in surpyval/univariate/nonparametric/_turnbull_npmle.py) comes from three results:

  • Existence. Call a set of pieces closed if every observation whose window meets it has a support that meets it. If every non-empty closed set meets every window, the maximum is attained (compactness). The condition is checked exactly in \(O((N + M)\log M)\). It is a sample counterpart of Woodroofe’s condition that every time carrying mass be observable with positive probability [Woodroofe1985np].

  • Exact observations. When every support is a single point this condition is the strong connectivity of the graph of Vardi and Wang [Vardi1985np] [Wang1991np], with an edge from \(a\) to \(b\) when an observation observable at \(a\) failed at \(b\). As they showed, if the windows fall into groups that do not overlap, each strongly connected, the NPMLE exists but is not unique (nothing links the groups’ total masses); if a group is connected but not strongly connected, it does not exist.

  • One-sided truncation, any censoring. With left truncation only, the likelihood is a polynomial in the discrete hazards, which always has a maximiser; the NPMLE exists unless every maximiser has a hazard of one before the last entry, where the survival would be zero before someone is observed. That can happen only at a gap: a time \(t\), before the last entry, such that nobody who had entered by \(t\) is known to have survived past it. If some support containing the gap ends before the last support starts, raising the hazard at the gap always raises the likelihood and the NPMLE does not exist; otherwise the likelihood does not depend on the gap’s hazard at all and the NPMLE exists but is not unique. With no gap it exists. Right truncation alone is the mirror image.

The first case of the last result is the collapse seen with left censoring and staggered entry: a left censored observation whose support reaches below the other observations’ entry times, where their conditional likelihoods cannot see it. It also covers the delayed-entry Kaplan-Meier that falls to zero when every item at risk fails before the next item enters. The second needs a censored observation whose support runs across the gap: a left or interval censored one, or an item right censored before the next item enters (one item entering at 0 and censored at 1, another entering at 5 and failing at 6: nothing in the data says how much probability lies in \((1, 5]\)). With exact and right censored observations only, SurPyval returns the delayed-entry Kaplan-Meier, the maximiser that puts none there. Entering every unit at a common time removes both.

With windows truncated on both sides and censored observations none of these arguments applies in full, and existence can depend on the counts, not just on which observation could have failed where. SurPyval reports "exists" when the first condition holds, "does not exist" when it finds a piece whose mass always raises the likelihood, and otherwise "undetermined"; then the EM’s convergence is the guide. For interval censored data with left truncation, Hudgens [Hudgens2005np] gives a necessary and sufficient condition for existence; the one-sided result above was derived and checked numerically on its own, not transcribed from that paper.

SurPyval warns when the verdict is "does not exist" (the estimate is not identifiable) or "not unique". On 240 simulated samples of 30 left truncated observations with every kind of censoring, the verdict agreed with what the EM does every time: all 63 with "does not exist" drifted to the boundary (the smallest window mass shrinking tenfold or more between 1,000 and 10,000 iterations), the one "not unique" sample settled in different places from different starts with the same likelihood, and the 176 with "exists" did not drift. The share of the fitted mass on pieces that some observation gains from and none pays for, still reported as exploitable_mass, ranged from 0.11 to 0.99 on the samples whose NPMLE does not exist, which is why it no longer decides the warning.

“Exists” rules out the free directions above, not every ambiguity: besides the placement of mass inside a piece, interval censored data can have a ridge of maxima with the same likelihood [GentlemanGeyer1994np]. Two of those 176 samples did.

A collapsed estimate. If the iteration produces a non-finite update, or (under truncation) the survival estimate collapses to essentially zero across the region the data can identify, SurPyval warns and sets degenerate = True on the model. The fitted model also records converged and iters. Treat any Turnbull estimate for which a warning was raised with suspicion: it is telling you that the data do not pin down a unique curve.

When to use it

The Turnbull estimator is the only non-parametric method that can handle left censoring, interval censoring, and right truncation (and arbitrary combinations of censoring and truncation). Left truncation / delayed entry on its own is handled by the Kaplan-Meier, Nelson-Aalen and Fleming-Harrington estimators too, via the tl keyword; it is only the left/interval censoring and right truncation that require Turnbull. Turnbull must therefore be used to supply the plotting positions in the parametric package whenever such data is present. For data the simpler estimators can handle, prefer them: they are exact rather than iterative, and their confidence bounds rest on firmer ground.

Comparing two groups: the log-rank test

Having estimated a survival curve for each of several groups, the natural next question is whether the groups genuinely differ, or whether the separation between the curves is just sampling noise. The log-rank test answers this. At every distinct event time it compares the number of failures observed in each group with the number expected if all groups shared one common survival curve, where the expected count in a group is the total number of failures at that time apportioned by the group’s share of the risk set. Summing the observed-minus-expected differences over all event times, and dividing by their variance, gives a statistic that is \(\chi^2\) distributed with \(k - 1\) degrees of freedom for \(k\) groups (counting, as R’s survdiff does, only groups with a positive expected number of failures; a group never at risk at a failure time carries no information). A small \(p\)-value is evidence the groups differ.

Concretely, at event time \(t\) let \(r_{gt}\) and \(d_{gt}\) be the number at risk and the number of failures in group \(g\), and \(r_t\) and \(d_t\) the pooled totals. Under the null hypothesis the failures at \(t\) are shared out like a draw without replacement from the risk set, so the expected count in group \(g\) is \(E_{gt} = d_t r_{gt}/r_t\) with hypergeometric variance \(d_t \frac{r_t - d_t}{r_t - 1} \frac{r_{gt}}{r_t}\left(1 - \frac{r_{gt}}{r_t}\right)\). The test statistic is built from \(\sum_t w_t (d_{gt} - E_{gt})\) and the corresponding (weighted) covariance, with one group dropped because the differences sum to zero. The test uses exact and right censored data.

The plain log-rank weights every event time equally (\(w_t = 1\)), which makes it most sensitive to proportional differences in hazard. Weighted variants change that emphasis: the Gehan-Breslow (\(w_t = r_t\)) and Tarone-Ware (\(w_t = \sqrt{r_t}\)) weights, and the Fleming-Harrington family \(w_t = \hat{S}(t-)^{\rho}(1 - \hat{S}(t-))^{\gamma}\) (with \(\hat{S}\) the pooled Kaplan-Meier just before \(t\)), up-weight early or late times so the test is more sensitive to differences concentrated there. Because the risk set is largest at the start, the Gehan-Breslow and Tarone-Ware weights emphasise early differences (Gehan more strongly), but they also depend on the censoring pattern, which the Fleming-Harrington weights do not. In the Fleming-Harrington family \(\rho > 0, \gamma = 0\) emphasises early differences (\(\rho = 1\) is close to the Peto-Peto test), \(\rho = 0, \gamma > 0\) late ones, and \(\rho = \gamma = 0\) is the plain log-rank. Choose the weighting before looking at the data: trying several and reporting the smallest \(p\)-value inflates the false-positive rate. See [KleinMoeschberger2003np] (chapter 7). (These Fleming-Harrington weights are a different thing from the Fleming-Harrington estimator above; the same two authors are behind both.)

When a nuisance factor (a site, a batch) also affects survival, comparing groups while ignoring it can be misleading if the groups are unevenly distributed across its levels. The stratified log-rank accumulates the observed-minus-expected numerators and their variances within each stratum before forming the statistic, so groups are only ever compared against others in the same stratum — the same logic as stratification in a Cox model, applied to the two-sample test. The degrees of freedom are still \(k - 1\), and with Fleming-Harrington weights the pooled Kaplan-Meier behind the weights is computed within each stratum.

Restricted mean survival time

A hazard ratio (from a log-rank test or a Cox model) only has a clean interpretation when the proportional-hazards assumption holds. When it does not — curves that cross, an effect that reverses over time — the restricted mean survival time (RMST) is an assumption-light summary. It is the area under the survival curve up to a horizon \(\tau\),

\[\text{RMST}(\tau) = \int_0^{\tau} S(t)\, dt,\]

which is exactly the average event-free time over the first \(\tau\) units and is always well defined. For a non-parametric estimate the integral is a sum of rectangles under the step function. Its variance is

\[\widehat{Var}\left(\widehat{\text{RMST}}(\tau)\right) = \sum_{i: x_i \leq \tau} A_i^2 \, v_i\]

where \(A_i\) is the area under the estimated curve from \(x_i\) to \(\tau\) and \(v_i\) is the increment at \(x_i\) of the estimator’s variance of \(\hat{H}\) (Greenwood’s increment \(d_i / (r_i(r_i - d_i))\) for the Kaplan-Meier, the corresponding increments for the Nelson-Aalen and Fleming-Harrington). For the Kaplan-Meier this is the standard estimator (see [KleinMoeschberger2003np], section 4.5); using the other estimators’ increments is the natural analogue. A Turnbull model uses the increments of its variance ladder, so with left or interval censoring the standard error inherits the approximation described in Uncertainty in a Turnbull estimate. Intuitively, uncertainty in the hazard at \(x_i\) moves the whole curve after \(x_i\), and hence the area \(A_i\). Comparing two groups by the difference in their RMST gives an effect measured in the natural units of time, with the variance of the difference being the sum of the two (independent) variances, and needs no assumption about the shape of, or relationship between, the two hazards. SurPyval’s intervals are the plain normal ones, \(\widehat{\text{RMST}} \pm z\,\widehat{SE}\), and the \(p\)-value for the difference is the two-sided \(z\)-test of a zero difference. The choice of \(\tau\) matters: it should be a time of practical interest at which both curves are still supported by data.

For worked examples — fitting the Kaplan-Meier, Nelson-Aalen, Fleming-Harrington and Turnbull estimators, confidence bounds and bands, Turnbull diagnostics, plotting positions and success-run testing, and comparing groups with the log-rank test and RMST difference — see the Non-Parametric SurPyval Modelling page.

References

[KM]

Kaplan, E. L., & Meier, P. (1958). Nonparametric estimation from incomplete observations. Journal of the American statistical association, 53(282), 457-481.

[NA]

Nelson, Wayne (1969). Hazard plotting for incomplete failure data. Journal of Quality Technology, 1(1), 27-52.

[FH] (1,2,3)

Fleming, Thomas R and Harrington, David P (1984). Nonparametric estimation of the survival distribution in censored data. Communications in Statistics-Theory and Methods, 13(20), 2469-2486.

[TB] (1,2)

Turnbull, Bruce W (1976). The empirical distribution function with arbitrarily grouped, censored and truncated data. Journal of the Royal Statistical Society: Series B (Methodological), 38(3), 290-295.

[Greenwood1926np]

Greenwood, M. (1926). The natural duration of cancer. Reports on Public Health and Medical Subjects, 33, 1-26. London: His Majesty’s Stationery Office.

[Aalen1978np]

Aalen, O. (1978). Nonparametric inference for a family of counting processes. The Annals of Statistics, 6(4), 701-726.

[Klein1991np]

Klein, J. P. (1991). Small sample moments of some estimators of the variance of the Kaplan-Meier and Nelson-Aalen estimators. Scandinavian Journal of Statistics, 18(4), 333-340.

[Efron1967np]

Efron, B. (1967). The two sample problem with censored data. Proceedings of the Fifth Berkeley Symposium on Mathematical Statistics and Probability, 4, 831-853.

[HallWellner1980np]

Hall, W. J., & Wellner, J. A. (1980). Confidence bands for a survival curve from censored data. Biometrika, 67(1), 133-143.

[Nair1984np]

Nair, V. N. (1984). Confidence bands for survival functions with censored data: a comparative study. Technometrics, 26(3), 265-275.

[BorganLiestol1990np]

Borgan, Ø. and Liestøl, K. (1990). A note on confidence intervals and bands for the survival function based on transformations. Scandinavian Journal of Statistics, 17(1), 35-41.

[BrookmeyerCrowley1982np]

Brookmeyer, R., & Crowley, J. (1982). A confidence interval for the median survival time. Biometrics, 38(1), 29-41.

[KleinMoeschberger2003np] (1,2,3,4)

Klein, J. P., & Moeschberger, M. L. (2003). Survival Analysis: Techniques for Censored and Truncated Data (2nd ed.). Springer.

[Filliben1975np]

Filliben, J. J. (1975). The probability plot correlation coefficient test for normality. Technometrics, 17(1), 111-117.

[Vardi1985np]

Vardi, Y. (1985). Empirical distributions in selection bias models. The Annals of Statistics, 13(1), 178-203.

[Wang1991np]

Wang, M.-C. (1991). Nonparametric estimation from cross-sectional survival data. Journal of the American Statistical Association, 86(413), 130-143.

[Woodroofe1985np]

Woodroofe, M. (1985). Estimating a distribution function with truncated data. The Annals of Statistics, 13(1), 163-177.

[Hudgens2005np]

Hudgens, M. G. (2005). On nonparametric maximum likelihood estimation with interval censoring and left truncation. Journal of the Royal Statistical Society: Series B, 67(4), 573-587.

[GentlemanGeyer1994np]

Gentleman, R., & Geyer, C. J. (1994). Maximum likelihood for interval censored data: consistency and computation. Biometrika, 81(3), 618-623.