Maximum Likelihood and EM Estimation for Cure Models

From Latent Susceptibility to an Estimable Likelihood

A mixture cure model represents the population as a mixture of two latent groups: people who are susceptible to the event and people who will never experience it. For a person whose follow-up ends without an event, group membership is unknown. The likelihood must therefore account for both explanations: the person may be cured, or may be susceptible but have an event time later than their censoring time. Maximum likelihood estimation handles this uncertainty directly; the expectation-maximization algorithm (EM) handles it by repeatedly using the posterior probability of susceptibility.

Let $Z_i=1$ indicate that individual $i$ is susceptible and $Z_i=0$ indicate cure. Write $\pi_i=P(Z_i=1\mid x_i)$ for the susceptibility probability given incidence covariates $x_i$; the cure probability is $1-\pi_i$. Conditional on $Z_i=1$, let the event time have density $f(t\mid w_i;\theta)$ and survival function $S(t\mid w_i;\theta)$, where $w_i$ are latency covariates and $\theta$ denotes latency parameters. A cured individual has no event at any finite follow-up time.

Definition: For right-censored data, $Y_i=\min(T_i,C_i)$ is the observed time and $\delta_i=I(T_i\leq C_i)$ indicates an observed event. Conditional on covariates, the observed-data likelihood contribution is $$ L_i(\gamma,\theta)= \left\{\pi_i f(Y_i\mid w_i;\theta)\right\}^{\delta_i} \left\{(1-\pi_i)+\pi_i S(Y_i\mid w_i;\theta)\right\}^{1-\delta_i}. $$ Here the incidence model specifies $\pi_i$, possibly through parameters $\gamma$. The expression assumes independent individuals and censoring independent of latent status and, among susceptible individuals, event time, conditional on modeled covariates. Thus the censoring-process likelihood factors out and is omitted.

An observed event requires susceptibility, so its contribution is the susceptibility probability times the conditional event density. For a censored observation, a cured individual contributes $1-\pi_i$, while a susceptible individual contributes $\pi_i S(Y_i)$. Adding these mutually exclusive possibilities gives the censored contribution. The full likelihood is $\prod_i L_i$, or, more conveniently for computation, the log likelihood $\ell(\gamma,\theta)=\sum_i\log L_i(\gamma,\theta)$.

Worked Likelihood Contributions

Worked Example: One Right-Censored Observation

Suppose a hypothetical individual is censored at 3 years. The fitted susceptibility probability is $\pi=0.60$, and susceptible event times follow an exponential distribution with rate $\lambda=0.20$ per year. The exponential survival function is $S(t)=e^{-\lambda t}$. The likelihood contribution is $$ L=(1-\pi)+\pi S(3)=0.40+0.60e^{-0.20(3)} =0.40+0.60(0.5488)=0.7293. $$

Among individuals with these model values who remain event-free to 3 years, the posterior probability of susceptibility is the susceptible component divided by the total: $$ q=P(Z=1\mid Y=3,\delta=0)= \frac{\pi S(3)}{(1-\pi)+\pi S(3)} =\frac{0.3293}{0.7293}=0.4515. $$

The posterior probability of cure is the other component, $1-q=0.5485$, or 54.85%. This is a conditional model probability for an individual with the specified follow-up and covariates; censoring does not establish that the individual is cured.

For observed events, the posterior susceptibility probability is 1, provided the modeled event density is positive. For a censored individual $i$, the general E-step probability is $$ q_i=P(Z_i=1\mid Y_i,\delta_i=0,x_i,w_i)= \frac{\pi_iS(Y_i\mid w_i;\theta)} {(1-\pi_i)+\pi_iS(Y_i\mid w_i;\theta)}. $$ Longer event-free follow-up generally lowers this probability when susceptible survival decreases with time. This posterior probability is the key quantity passed from the E step to the M step.

Deriving the EM Updates

If susceptibility status were observed, the complete-data likelihood would be simpler. For each individual it includes a Bernoulli contribution for $Z_i$, and, for susceptible individuals, the appropriate event or survival contribution. Its log likelihood is $$ \ell_c(\gamma,\theta)= \sum_i\left[ Z_i\log(\pi_i)+(1-Z_i)\log(1-\pi_i) +Z_i\{\delta_i\log f(Y_i\mid w_i;\theta) +(1-\delta_i)\log S(Y_i\mid w_i;\theta)\} \right]. $$

Because event observations force $Z_i=1$, their expected susceptibility indicator is 1. For censored observations, its conditional expectation is $q_i$ above. In the E step, calculate these expectations using the current parameter values. In the M step, maximize the expected complete-data log likelihood, replacing each $Z_i$ by its current expectation $q_i$. The resulting objective separates into incidence and latency terms: $$ Q(\gamma,\theta)= \sum_i\{q_i\log(\pi_i)+(1-q_i)\log(1-\pi_i)\} +\sum_iq_i\{\delta_i\log f(Y_i\mid w_i;\theta) +(1-\delta_i)\log S(Y_i\mid w_i;\theta)\}. $$ For an observed event set $q_i=1$; for a censored observation use the posterior probability. Maximize the first sum for the incidence model and the second for the latency model. With logistic incidence, the first is a weighted Bernoulli log likelihood with fractional outcomes $q_i$; standard logistic-regression optimization can maximize it. The latency update is a weighted survival log likelihood, with the same $q_i$ weights. In general these maximizations may require numerical optimization.

Justification for likelihood ascent: Let $Z$ denote all missing susceptibility indicators, and let the current parameters be $\psi^{(r)}=(\gamma^{(r)},\theta^{(r)})$. The observed log-likelihood difference between a proposed $\psi$ and the current value can be written $$ \ell(\psi)-\ell(\psi^{(r)} )=Q(\psi\mid\psi^{(r)})-Q(\psi^{(r)}\mid\psi^{(r)}) +\operatorname{KL}\!\left\{ P_{\psi^{(r)}}(Z\mid Y,\delta)\,\middle\|\, P_{\psi}(Z\mid Y,\delta) \right\}. $$ The Kullback–Leibler divergence $\operatorname{KL}$ is nonnegative. Thus, if the M step increases or maximizes $Q$, the observed log likelihood cannot decrease. This guarantees ascent, not that EM reaches the global maximum.

A Complete Exponential-Latency EM Iteration

Worked Example: Updating Susceptibility and the Event Rate

Consider four hypothetical independent individuals with no covariates. Two experience events at 2 and 5 years; two are right-censored at 3 and 8 years. Assume a common susceptibility probability $\pi$, exponential susceptible event times with rate $\lambda$ per year, and noninformative censoring. Start the EM algorithm at $\pi^{(0)}=0.60$ and $\lambda^{(0)}=0.20$ per year. These are valid interior starting values, and all observed times have positive event densities or survival probabilities under the model.

E step. The two observed events have $q=1$. For the censored observation at 3 years, $$ q_3=\frac{0.60e^{-0.20(3)}}{0.40+0.60e^{-0.20(3)}}=0.4515. $$ For the one censored at 8 years, $$ q_8=\frac{0.60e^{-0.20(8)}}{0.40+0.60e^{-0.20(8)}} =\frac{0.1211}{0.5211}=0.2324. $$ The posterior cure probabilities for these two individuals are respectively $1-0.4515=0.5485$ and $1-0.2324=0.7676$.

M step for incidence. With a common susceptibility probability, maximize the Bernoulli part of $Q$. Its maximizer is the average expected susceptibility: $$ \pi^{(1)}=\frac{1+1+q_3+q_8}{4} =\frac{2+0.4515+0.2324}{4}=0.6710. $$

M step for latency. For an exponential distribution, $\log f(y)=\log\lambda-\lambda y$ and $\log S(y)=-\lambda y$. Thus the latency part of $Q$ is $$ 2\log\lambda-\lambda\{2+5+3q_3+8q_8\}. $$ Differentiating with respect to $\lambda$ gives $2/\lambda-(7+3q_3+8q_8)$. Setting this equal to zero yields $$ \lambda^{(1)}=\frac{2}{7+3(0.4515)+8(0.2324)} =\frac{2}{10.214}\approx0.1958\text{ per year}. $$ The two observed events contribute their full event-time exposure; censored observations contribute exposure weighted by their posterior susceptibility.

Checking the Update and Convergence

Worked Example: Verifying the Observed Likelihood

For the same four observations, the observed log likelihood at any $(\pi,\lambda)$ is $$ \ell(\pi,\lambda)=2\log(\pi\lambda)-7\lambda +\log\{1-\pi+\pi e^{-3\lambda}\} +\log\{1-\pi+\pi e^{-8\lambda}\}. $$ The first term accounts for the two event densities and their combined event times, 2 and 5 years; the last two terms are the censored mixture contributions.

At the starting values, substituting $\pi=0.60$ and $\lambda=0.20$ gives $$ \ell(0.60,0.20)=2\log(0.12)-1.40+\log(0.7293)+\log(0.5211) \approx-6.608. $$ At the updated values, using $\pi=0.6710$ and $\lambda=0.1958$ gives $$ \begin{aligned} \ell(0.6710,0.1958) &=2\log(0.6710\times0.1958)-7(0.1958)\\ &\quad+\log(0.3290+0.6710e^{-3(0.1958)}) +\log(0.3290+0.6710e^{-8(0.1958)})\\ &\approx-6.541. \end{aligned} $$ The observed log likelihood has increased by about $0.0672$, as expected for an EM step. The updated parameters have not necessarily maximized the observed likelihood; further E and M steps are needed until a stated convergence criterion is met.

A practical implementation should track the observed log likelihood, parameter changes, and the size of the EM update. A substantial decrease in observed likelihood signals an implementation error, an incorrectly maximized M step, or a numerical problem. Tiny changes in likelihood alone can be misleading if parameters are still moving appreciably, so convergence checks commonly consider both likelihood and parameter changes. Multiple starting values are useful because different starts can lead to different stationary points.

Direct Maximization, Limitations, and Common Mistakes

Direct maximum likelihood maximizes the observed-data log likelihood, including the mixture terms for censored observations. It can converge quickly, but its objective may be difficult to optimize and its gradients require careful handling. EM replaces that problem with simpler weighted maximizations and naturally incorporates uncertainty about susceptibility. Its main costs are that it may converge slowly near a maximum and that monotone likelihood ascent does not establish a global maximum. Neither approach removes the need to check model assumptions, parameter constraints, and sensitivity to initialization.

Do not label every censored individual cured. A censoring time supplies incomplete information. Use the posterior probability $q_i$, not a deterministic cure classification.
Do not use susceptible survival alone for a censored observation. Its likelihood contribution is $(1-\pi_i)+\pi_iS(Y_i)$, because cure is also possible.
Do not set a censored subject’s latency weight to 1. The expected complete-data contribution is weighted by $q_i$, which may be far below 1.
Do not treat convergence as proof of identification or a global optimum. Limited follow-up can leave the cure fraction and long-term latency poorly distinguished; a stable algorithm can still estimate them imprecisely or settle at a local solution.

The parametric latency model must be appropriate for susceptible event times, and incidence and latency components must be estimable from the available covariates and follow-up. The introductory material on identifiability and follow-up explains why observing enough of the event-time tail matters: EM is an estimation algorithm, not a remedy for insufficient information. A fitted cure fraction should therefore be interpreted with the model and follow-up limits in view.

Check Your Understanding

Use the definitions and update equations above. Solutions follow the questions.

Write the likelihood contribution for a censored individual with observed time $y$, susceptibility probability $\pi$, and susceptible survival $S(y)$.
For $\pi=0.50$, exponential rate $\lambda=0.25$ per year, and censoring at 4 years, calculate the posterior susceptibility and cure probabilities.
An EM E step gives $q=1,1,0.40,0.20$ for four subjects. Two events occur at 1 and 3 years; the other two are censored at 2 and 4 years. Calculate the no-covariate incidence and exponential-rate updates.
Why does an increase in observed log likelihood at every EM iteration not prove that the algorithm found the global maximum?

Solutions. (1) The contribution is $(1-\pi)+\pi S(y)$, the sum of the cured and susceptible-survivor possibilities. (2) $S(4)=e^{-0.25(4)}=e^{-1}=0.3679$. Thus $q=0.50(0.3679)/\{0.50+0.50(0.3679)\}=0.2690$, and the posterior cure probability is $1-q=0.7310$. (3) The incidence update is $(1+1+0.40+0.20)/4=0.6500$. The latency exposure is $1+3+2(0.40)+4(0.20)=5.60$ years, so $\lambda_{\text{new}}=2/5.60=0.3571$ per year. (4) EM guarantees nondecreasing likelihood under valid E and M steps, but it may converge to a local maximum or a stationary point; likelihood ascent alone does not compare all possible solutions.