📊 statistics

A fast and stable algorithm for non-parametric maximum likelihood estimation of survival functions for left-truncated and interval-censored data

This paper introduces a fast and stable product-limit style EM algorithm combined with a modified iterative convex minorant step to efficiently compute the non-parametric maximum likelihood estimator for survival functions using left-truncated and interval-censored data, demonstrating superior convergence and scalability over existing methods.

Original authors: Zachary Waller, Adele H. Marshall, Frank Kee, Felicity Lamrock

Published 2026-08-09

📖 4 min read☕ Coffee break read

CC BY 4.0

Original authors: Zachary Waller, Adele H. Marshall, Frank Kee, Felicity Lamrock

Original paper licensed under CC BY 4.0 (https://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are a detective trying to figure out exactly when a specific event happens in a group of people, like when a secret club member finally decides to quit. But there's a catch: you don't get to see the moment they quit. You only get to peek at them at random times. Sometimes you look, and they are still there; the next time you look, they are gone. You know they quit sometime between those two peeks, but you don't know the exact second. This is called "interval censoring."

Now, add a second twist. Imagine you only start watching these people after they have already been in the club for a while. If someone quit before you started your watch, you never even knew they existed. This is called "left truncation." It's like trying to guess the lifespan of a tree, but you only start measuring it once it's already ten feet tall, and you only check on it every few years to see if it's still standing.

Scientists who study survival—like how long patients stay healthy or how long machines keep working—face this exact puzzle. They need a mathematical way to draw a map of time that shows the probability of an event happening, even when their data is full of these "I don't know exactly when" gaps and "I wasn't watching from the start" holes. The problem is that the old maps they used were incredibly slow to draw and sometimes got stuck in a loop, unable to find the best answer. If you want to know how confident you can be in these maps, you have to redraw them thousands of times, which makes the old, slow methods impossible to use for big, complicated real-world problems.

This paper introduces a new, super-fast detective tool called the "Product-Limit" (PL) algorithm, which works like a clever shortcut to solve this puzzle. The authors, researchers from Queen's University Belfast, realized that instead of treating the "missing time" as a messy mystery, they could reorganize the math to look more like a famous, simple method used for easier data. They call this a "reparameterization," which is just a fancy way of saying they changed the way they asked the question to make it easier to answer.

Think of the old way of solving this as trying to fill a giant, leaky bucket by pouring water in one drop at a time, hoping you eventually get it full. It works, but it takes forever, and if the bucket has a big hole (heavy truncation), the water might never stay in. The new PL algorithm is like realizing you can just plug the hole first and then pour the water in a steady stream. By treating the "start watching" times and the "stop watching" times as exact moments (which they are), and only using the complex math for the "I don't know exactly when" gaps, the new method skips the slow, repetitive steps.

The researchers tested this new tool against nine other existing methods using computer simulations. They created thousands of fake scenarios with different levels of missing data and "late starts." The results were clear: the new PL algorithm, especially when combined with a second step called "ICM," was dramatically faster and more stable than the others. In some tests, it was hundreds of times faster. While the old methods sometimes gave up or got stuck in a loop, the new method kept marching forward, finding the best possible map every time.

To prove it worked in the real world, the team applied their new algorithm to a famous dataset about older adults losing their ability to do daily tasks (like bathing or dressing). This data is tricky because the study only started watching people once they were already 65, and they only checked in every few years. The old methods took over 20 seconds to draw a map and sometimes got stuck after a million tries without finding the best answer. The new PL-ICM algorithm did the same job in a fraction of a second (0.003 seconds for women and 0.002 seconds for men) and found a more accurate map.

The paper suggests that this new approach is a game-changer for handling messy survival data. It doesn't just speed things up; it solves problems that other methods simply couldn't handle, allowing scientists to draw clearer, more reliable pictures of how time affects events, even when their data is full of gaps and late starts. The authors are confident that this method is ready to be used for complex studies, potentially helping researchers understand everything from disease progression to machine failure much more quickly and accurately.

Technical Summary: A Fast and Stable Algorithm for Non-Parametric Maximum Likelihood Estimation of Survival Functions for Left-Truncated and Interval-Censored Data

Problem Statement

The paper addresses the computational challenges associated with calculating the Non-Parametric Maximum Likelihood Estimator (NPMLE) for survival functions when data is both left-truncated (LT) and interval-censored (IC).

Interval Censoring: Events are known only to have occurred within a specific time interval $(L, R]$ rather than at an exact time.
Left Truncation: Observation of an event is conditional on the event occurring after a specific entry time $T$ (e.g., a participant enters a study only after reaching a certain age).

While the Kaplan-Meier estimator provides an analytical solution for right-censored and left-truncated data, the NPMLE for LTIC data generally lacks a closed-form solution. Existing algorithms, such as Turnbull's self-consistency algorithm, often suffer from extremely slow convergence or numerical instability, particularly when heavy truncation is present. Furthermore, standard methods for calculating confidence intervals (e.g., bootstrapping or likelihood-ratio tests) require repeated NPMLE calculations, making computational efficiency critical.

Methodology

1. Reparameterization as a Product-Limit (PL)

The authors propose a novel Expectation-Maximization (EM) algorithm, termed the Product-Limit (PL) algorithm, which reparameterizes the likelihood function.

Inner-Intervals: The algorithm operates on "inner-intervals" (or Turnbull-intervals), which are the disjoint intervals where the survival function is allowed to change.
Interval Hazard ( $h_j$ ): Instead of estimating the probability mass directly, the likelihood is reparameterized in terms of an "interval hazard" $h_j = P(X \in Q_j | X > l_j)$ , where $Q_j$ is the $j$ -th inner-interval.
Analytical Handling of Truncation: This reparameterization allows left-truncation and right-censoring to be handled analytically, similar to the Kaplan-Meier estimator. The survival function is expressed as a product: $S(r_j) = \prod (1 - h_k)$ .
Iterative Updates: The algorithm iteratively updates $h_j$ using a ratio of expected events to a risk-set estimate ( $d_j / R_j$ ). The authors implement an $O(n)$ per-iteration complexity using recursive sums, similar to efficient implementations of Turnbull's algorithm for interval-censored data.

2. Hybrid PL-ICM Algorithm

To further improve performance, the PL algorithm is combined with a modified Iterative Convex Minorant (ICM) step.

Mechanism: The hybrid approach alternates between steps of the PL-EM algorithm and steps of the ICM algorithm (adapted for LTIC data by Pan, 1999).
Stability: The ICM step utilizes a Newton-like update with a Levenberg-Marquardt adjustment to ensure the Hessian is negative-definite, followed by a Weighted Pool-Adjacent-Violators Algorithm (WPAVA) to enforce monotonicity of the cumulative hazard.
Strategy: The authors employ 10 PL steps for every 1 ICM step, leveraging the low computational cost of the PL updates to accelerate convergence.

3. Handling Non-Existence of Maximizers

The paper addresses scenarios where the NPMLE does not exist under traditional parameterizations (e.g., when the survival function drops to zero prematurely due to heavy truncation).

Buffering: A small buffer $\eta$ is introduced to prevent $h_j$ from reaching exactly 1, ensuring numerical stability.
Conditional Estimation: If $\hat{h}_j = 1$ for some $j < m$ , the algorithm identifies that the maximizer exists only in the limit. The authors propose conditioning the survival estimate on survival past time $r_j$ to obtain a valid estimate, a process that is straightforward under the PL parameterization but difficult under Turnbull's.

Key Results

Simulation Study

The authors compared the proposed PL and PL-ICM algorithms against nine other methods (including Turnbull, Shen, Yu, Breslow, ICM, and Quasi-Newton) across simulations with 0%, 25%, and 50% truncation rates.

Speed: The PL-ICM algorithm was consistently the fastest, often converging in fewer than 100 iterations and milliseconds of CPU time. The PL algorithm was the fastest pure EM method.
Stability: Under heavy truncation (50%), other EM algorithms (Turnbull, Shen, Yu) showed increased iteration counts or failure to converge. The Yu algorithm frequently left the parameter space, resulting in non-monotonic survival functions.
Accuracy: The PL and PL-ICM methods achieved the highest log-likelihoods, matching or exceeding the performance of other stable methods while requiring significantly fewer iterations.

Real-World Application (MHCPS)

The algorithms were applied to the Massachusetts Health Care Panel Study (MHCPS) data, a dataset known for causing convergence issues in previous studies.

Convergence: The proposed PL-ICM algorithm converged rapidly (77–99 iterations) to a log-likelihood of -435.67 (males) and -597.41 (females).
Comparison:
- Turnbull/Shen: Failed to converge within 1 million iterations for females and stopped prematurely for males, yielding lower log-likelihoods.
- ICM: Converged for males but failed for females within the iteration limit.
- Quasi-Newton: Converged faster than ICM but failed to reach the maximum likelihood achieved by PL-ICM.
Insight: The PL-ICM algorithm identified a specific time point (age 65.3) where the estimated hazard was 1, indicating a precipitous drop in survival. This allowed for a principled conditional survival estimate, whereas the Turnbull algorithm obscured this feature due to its slow convergence and "ghost" observation mechanism.

Significance and Claims

The paper claims that the proposed PL-ICM algorithm represents a significant advancement in the analysis of LTIC data for the following reasons:

Computational Efficiency: It is the first algorithm to combine a product-limit reparameterization with an ICM step specifically for LTIC data, achieving linear time complexity per iteration ( $O(n)$ ) and drastically reducing the number of iterations required for convergence compared to existing methods.
Stability under Heavy Truncation: Unlike Turnbull's algorithm, which adds "ghost" observations that can degrade convergence, the PL approach accounts for truncation analytically. This makes it robust in scenarios with heavy truncation where other methods fail or become unstable.
Handling Non-Existence: The reparameterization allows for the detection of cases where the NPMLE does not exist in the standard parameterization (i.e., when the survival function drops to zero). The method facilitates the calculation of conditional survival estimates in these cases, providing a solution where traditional methods yield non-identifiable results.
Practical Utility: By significantly reducing the computational cost of calculating the NPMLE, the algorithm makes resampling methods (like bootstrapping for confidence intervals) and likelihood-ratio tests feasible for complex LTIC datasets, a task previously considered impractical.

The authors conclude that their method offers a fast, stable, and conceptually appealing alternative to existing algorithms, enabling the analysis of larger and more complex LTIC datasets than previously practical.