EM Computation of Interval-Censored Distributions

Why Use EM to Compute the Turnbull Estimate?

With interval-censored data, the event time is not observed exactly. Each subject instead contributes an interval known to contain the event time, or a one-sided interval for censoring. The Turnbull estimator represents the event-time distribution by probability masses on support cells compatible with the observed intervals. Its likelihood is straightforward to write down, but maximizing it directly can involve many overlapping compatibility patterns.

The expectation-maximization (EM) algorithm provides a systematic way to compute a nonparametric maximum likelihood estimate (NPMLE). Its key idea is to assign each observation a fractional amount of probability across the support cells it could have come from. These fractions depend on the current mass estimates; summing the allocations and renormalizing gives the next estimates.

Definition: Let \(J_1,\ldots,J_m\) denote the support cells used in the Turnbull representation, with probability masses \(p_1,\ldots,p_m\), where \(p_j\geq 0\) and \(\sum_{j=1}^m p_j=1\). For observation \(i\), let \(C_i\) be its event-time set under the study’s interval and censoring conventions. Its compatibility set is \(A_i=\{j:J_j\text{ is compatible with }C_i\}\). The probability assigned to observation \(i\) is \(\sum_{j\in A_i}p_j\).

The support cells are those produced by the Turnbull support construction, such as the maximal-intersection representation. Compatibility must be determined using the same endpoint convention as the observation likelihood: an event interval, a left-censored observation, and a right-censored observation can imply different inclusion rules at their boundaries. Once the cells and compatibility sets have been specified, the computation uses only those sets and the masses.

For \(n\) independent observations, the observed-data likelihood and log-likelihood, apart from any factors not involving the masses, are

L(\mathbf p)=\prod_{i=1}^n\left(\sum_{j\in A_i}p_j\right), \qquad \ell(\mathbf p)=\sum_{i=1}^n\log\left(\sum_{j\in A_i}p_j\right).

This likelihood is maximized over the probability simplex. The EM computation described here assumes independent observations, a valid support representation, and a nonempty compatibility set for every observation. Initialize with masses that sum to one and are strictly positive on every support cell being considered. This ensures each observation has a positive fitted probability at the start and allows the algorithm to allocate mass to any compatible cell.

Deriving the Fractional-Allocation Update

Imagine, temporarily, that each observation also carried an unobserved label indicating which compatible support cell generated its event time. If the labels were observed, maximizing the likelihood would amount to assigning probability to each cell in proportion to the number of labels assigned there. Since the labels are missing, EM replaces each label count by its conditional expected value under the current masses.

For current masses \(\mathbf p^{(r)}\), the conditional probability that observation \(i\) belongs to cell \(j\) is zero if \(j\notin A_i\). If \(j\in A_i\), it is proportional to \(p_j^{(r)}\), normalized over the compatible cells:

w_{ij}^{(r)} = \begin{cases} \dfrac{p_j^{(r)}}{\sum_{k\in A_i}p_k^{(r)}}, & j\in A_i,\\[6pt] 0, & j\notin A_i. \end{cases}

These are fractional allocations: for each observation, they are nonnegative and sum to one. The expected number assigned to cell \(j\) is \(\sum_i w_{ij}^{(r)}\). Dividing by \(n\) converts those expected counts into updated probability masses:

EM update: For \(j=1,\ldots,m\),

p_j^{(r+1)} = \frac{1}{n}\sum_{i=1}^n w_{ij}^{(r)} = \frac{1}{n}\sum_{i:j\in A_i} \frac{p_j^{(r)}}{\sum_{k\in A_i}p_k^{(r)}}.

The update automatically satisfies the probability constraint. In fact,

\sum_{j=1}^m p_j^{(r+1)} = \frac{1}{n}\sum_{i=1}^n\sum_{j=1}^m w_{ij}^{(r)} = \frac{1}{n}\sum_{i=1}^n 1 =1.

This derivation is also the computational interpretation of the Turnbull self-consistency equations: each observation contributes its current conditional allocation to the cells it can occupy. If an iteration assigns zero mass to a cell, later iterations keep that mass at zero. Thus a strictly positive initialization avoids prematurely excluding a cell that might receive positive mass at a maximum.

Worked Example: One EM Update From Three Compatibility Sets

Consider three support cells and three hypothetical independent observations with compatibility sets \(A_1=\{1,2\}\), \(A_2=\{2,3\}\), and \(A_3=\{1,2,3\}\). Start at \(p_1^{(0)}=p_2^{(0)}=p_3^{(0)}=1/3\). The broad third observation is compatible with every cell.

For observation 1, the two compatible masses are equal, so its allocations are \(w_{11}^{(0)}=1/2\) and \(w_{12}^{(0)}=1/2\). Observation 2 allocates \(1/2\) to cell 2 and \(1/2\) to cell 3. Observation 3 allocates \(1/3\) to each cell. The expected counts are therefore

\begin{aligned} N_1^{(0)}&=\frac12+\frac13=\frac56,\\ N_2^{(0)}&=\frac12+\frac12+\frac13=\frac43,\\ N_3^{(0)}&=\frac12+\frac13=\frac56. \end{aligned}

There are \(n=3\) observations, so the updated masses are

(p_1^{(1)},p_2^{(1)},p_3^{(1)}) = \left(\frac{5}{18},\frac49,\frac{5}{18}\right) \approx (0.2778,0.4444,0.2778).

They sum to one. The initial likelihood is \((2/3)(2/3)(1)=4/9\). At the updated masses, the first and second observations each have fitted probability \(13/18\), and the third has probability one. Thus \(L(\mathbf p^{(1)})=(13/18)^2\approx 0.5216\), greater than the initial \(4/9\approx 0.4444\). One iteration has improved the fit, but it has not established that the algorithm has converged.

Iteration, Convergence, and Boundary Estimates

The EM algorithm alternates two calculations: compute the observation-level fractional allocations, then average them to form new masses. It can be implemented directly from a list of compatibility sets. The log-likelihood is often more convenient than the likelihood for monitoring progress, especially when the number of observations is large.

Build and check the support.
Construct the Turnbull support cells and record each observation’s compatible cells, using the intended endpoint rules. Confirm that no compatibility set is empty.

Initialize the masses.
Choose nonnegative masses summing to one, typically strictly positive across the candidate cells. Record the initial log-likelihood.

Allocate and update.
For each observation, divide each compatible current mass by the sum of compatible masses. Average the resulting allocations over observations to get the next masses.

Monitor and stop.
Recompute the log-likelihood and check mass constraints and allocation sums. Stop when the change in log-likelihood and/or the largest change in any mass is below a prespecified numerical tolerance.

Under the stated model, each EM iteration does not decrease the observed-data likelihood. The reason is that the allocation step computes the conditional distribution of the missing cell labels under the current masses, while the update maximizes the expected complete-data log-likelihood over the probability simplex. The general EM likelihood inequality then gives \(L(\mathbf p^{(r+1)})\geq L(\mathbf p^{(r)})\). The log-likelihood has the same monotonicity because the logarithm is increasing.

The log-likelihood is concave in the masses: each term is the logarithm of a nonnegative linear function of \(\mathbf p\). Consequently, any interior stationary point is a global maximum; boundary maxima are also possible. EM is a way to approach a maximum, not a guarantee that all masses will be strictly positive or that convergence will be rapid. The update can become small near a boundary, and the likelihood can be flat in some directions. A small change between iterations is a numerical stopping signal, not proof of a unique estimate.

Worked Example: EM Approaches a Boundary Maximum

Suppose four hypothetical observations have compatibility sets \(A_1=\{1\}\), \(A_2=\{1,2\}\), \(A_3=\{2,3\}\), and \(A_4=\{3\}\). Initialize at \((p_1^{(0)},p_2^{(0)},p_3^{(0)})=(1/4,1/2,1/4)\). The two singleton observations allocate all their mass to cells 1 and 3. Observation 2 allocates \(1/3\) to cell 1 and \(2/3\) to cell 2; observation 3 allocates \(2/3\) to cell 2 and \(1/3\) to cell 3. Expected counts are \(4/3\) for each cell, giving

(p_1^{(1)},p_2^{(1)},p_3^{(1)})=(1/3,1/3,1/3).

At these equal masses, the two non-singleton observations each split their allocation equally. The next expected counts are \(3/2,1,3/2\), so

(p_1^{(2)},p_2^{(2)},p_3^{(2)})=(3/8,1/4,3/8).

The initial likelihood is \((1/4)(3/4)(3/4)(1/4)=9/256\approx0.03516\). After the second update it is \((3/8)(5/8)(5/8)(3/8)=225/4096\approx0.05493\). Both iterations improve the likelihood.

Here the boundary solution can also be checked directly. By symmetry, set \(p_1=p_3=x\) and \(p_2=1-2x\), where \(0\leq x\leq1/2\). The likelihood is \(x^2(1-x)^2\), which increases on this range because its derivative is \(2x(1-x)(1-2x)\geq0\). Its maximum occurs at \(x=1/2\), giving masses \((1/2,0,1/2)\) and likelihood \(1/16=0.0625\). The EM iterates move toward this boundary maximum. The zero middle mass is an estimate, not a computational failure.

Worked Example: Allocations for a Broad and a Right-Censored Observation

Suppose three support cells represent ordered event-time regions, and two hypothetical observations have compatibility sets \(A_1=\{1,2,3\}\) and \(A_2=\{2,3\}\). The second set could arise from a right-censored observation whose event time is known to be later than a threshold. Begin with masses \((0.2,0.3,0.5)\).

Observation 1 is compatible everywhere, so its allocation is \((0.2,0.3,0.5)\). For observation 2, the compatible total mass is \(0.3+0.5=0.8\), giving allocations \(0.3/0.8=0.375\) to cell 2 and \(0.5/0.8=0.625\) to cell 3. Summing allocations and dividing by \(n=2\) gives

(p_1^{(1)},p_2^{(1)},p_3^{(1)}) = \left(\frac{0.2}{2},\frac{0.3+0.375}{2},\frac{0.5+0.625}{2}\right) =(0.1,0.3375,0.5625).

The masses still sum to one. The initial likelihood is \(1(0.8)=0.8\); the updated likelihood is \(1(0.3375+0.5625)=0.9\). This illustrates why the algorithm does not assign an entire observation to one arbitrarily chosen compatible cell: it distributes its contribution in proportion to the current masses.

Diagnostics, Alternatives, and Limitations

A reliable implementation checks both the arithmetic and the statistical setup. For every observation at every iteration, the allocations across compatible cells should sum to one. Updated masses should be nonnegative and sum to one, up to numerical rounding. The observed-data log-likelihood should not decrease beyond small numerical error. A decrease can indicate an incorrect compatibility set, an update coded with the wrong denominator, or a normalization error.

Zero fitted probability: if \(\sum_{k\in A_i}p_k^{(r)}=0\), the allocation for observation \(i\) is undefined. Check whether its compatibility set is empty or whether initialization assigned zero to every compatible cell. A positive initialization over all candidate cells prevents this at the first step.
Unintended zeros: a cell initialized at zero remains at zero. If the support is being explored, do not set candidate masses to zero before establishing that those cells can be excluded.
Apparent nonconvergence: slow progress, especially near a boundary, may require many iterations. Compare both log-likelihood changes and mass changes; a stable likelihood alone can coexist with noticeable movement along a nearly flat direction.
Endpoint mistakes: a cell’s compatibility must follow the event-set convention used in the likelihood. Inconsistent treatment of interval endpoints changes the likelihood and therefore changes the estimate.

Direct constrained optimization of \(\ell(\mathbf p)\) over the probability simplex is the main computational alternative. It can provide a useful check on the EM result and may converge faster with a well-chosen optimizer, but it requires explicit constraints and careful treatment of boundary solutions. EM has simple updates, preserves the probability constraint, and exposes the fractional contributions observation by observation. It can, however, converge slowly. The self-consistency equations give the corresponding fixed-point condition; EM supplies a practical iterative route to satisfying it.

The fitted masses describe how the likelihood represents the distribution over the chosen support cells. They do not recover exact event times inside an interval that was never observed, and a cell may represent a region rather than a uniquely located event time. Likewise, a converged iteration and a high likelihood do not by themselves establish that the distribution is uniquely determined by the observed inspection pattern. Treat the resulting masses as the Turnbull estimate under the stated support and observation conventions.

Common mistake: Assigning each interval-censored observation entirely to the compatible cell with the largest current mass, then renormalizing, is not the EM update. That procedure replaces conditional probabilities with hard assignments and can change the likelihood being optimized. The correction is to allocate fractionally: each compatible cell receives its current mass divided by the total current mass over that observation’s compatibility set.

Check Your Understanding

Use the EM update and the compatibility sets stated in each question.

For one observation compatible with cells 1 and 3, the current masses are \(p_1=0.2\), \(p_2=0.5\), and \(p_3=0.3\). What fractional allocations does the observation contribute?
Two observations have compatibility sets \(A_1=\{1,2\}\) and \(A_2=\{2,3\}\), with current masses \((0.2,0.5,0.3)\). Calculate one EM update.
Why is a strictly positive initialization useful when the candidate support cells are still under consideration?
If the log-likelihood has increased and the largest change in mass is below the chosen tolerance, what can be concluded, and what cannot be concluded?

Solutions. (1) The compatible total is \(0.2+0.3=0.5\). The allocations are \(0.2/0.5=0.4\) to cell 1 and \(0.3/0.5=0.6\) to cell 3; cell 2 receives zero. (2) Observation 1 allocates \(0.2/0.7=2/7\) to cell 1 and \(5/7\) to cell 2. Observation 2 allocates \(0.5/0.8=5/8\) to cell 2 and \(0.3/0.8=3/8\) to cell 3. The expected counts are \(2/7\), \(5/7+5/8=75/56\), and \(3/8\). Divide by two: the updated masses are \((1/7,75/112,3/16)\), approximately \((0.1429,0.6696,0.1875)\), which sum to one. (3) A cell initialized at zero stays at zero under the EM update, so positive starting masses avoid excluding a candidate cell solely because of initialization. (4) The numerical sequence has met the selected stopping criterion and is consistent with convergence to a likelihood maximum. This does not prove a unique maximum, exact recovery of event times within support cells, or that the tolerance is adequate for every inferential purpose.

EM Computation of Interval-Censored Distributions

What You'll Learn