An EM Algorithm for Latent Class Joint Models

Why Use EM for a Latent Class Joint Model?

In a latent class joint model, each subject belongs to one of a fixed number of unobserved classes. The class affects both the longitudinal measurements and the event-time distribution. Because class membership is not observed, the likelihood sums over possible classes for each subject. That sum can make direct maximization inconvenient. The expectation-maximization (EM) algorithm handles it by alternating between estimating each subject’s class probabilities and fitting the model components using those probabilities as weights.

The weights are not hard assignments: a subject may contribute partly to several classes. This preserves uncertainty about class membership during estimation. The algorithm is an optimization method for the specified likelihood; by itself, it does not establish that the chosen number of classes or model form is appropriate.

Model and Assumptions

Suppose there are \(n\) independent subjects and \(G\) classes. Subject \(i\) has longitudinal data \(y_i=(y_{i1},\ldots,y_{im_i})\), event or censoring time \(t_i\), event indicator \(\delta_i\), and possibly baseline covariates \(x_i\). Here \(\delta_i=1\) denotes an event and \(\delta_i=0\) denotes right censoring. Let \(Z_{ig}=1\) if subject \(i\) belongs to class \(g\), and \(Z_{ig}=0\) otherwise. Exactly one \(Z_{ig}\) is 1 for each subject.

Let \(\pi_{ig}(\alpha)=P(Z_{ig}=1\mid x_i)\) be the class probability, with \(\sum_{g=1}^G\pi_{ig}(\alpha)=1\). Write \(f_g(y_i\mid x_i;\theta_g)\) for the longitudinal-data contribution and \(h_g(t_i,\delta_i\mid x_i;\psi_g)\) for the event-data contribution in class \(g\). Their product is the class-specific contribution \(L_{ig}\). The model assumes that, conditional on class and included covariates, the joint contribution factors into these longitudinal and event components. Subjects are independent, and censoring is noninformative for the event model conditional on the modeled information.

For a right-censored event time with class-specific hazard \(\lambda_g(t\mid x_i;\psi_g)\) and cumulative hazard \(\Lambda_g(t\mid x_i;\psi_g)\), the event contribution is

The observed-data likelihood and log likelihood are

where \(\Theta\) collects all model parameters. The EM algorithm treats the unobserved indicators \(Z_{ig}\) as missing data. The model assumptions must be considered before applying the algorithm: EM maximizes the likelihood that was specified, not a more general model that might better describe dependence or censoring.

The E-Step: Class Responsibilities

At iteration \(r\), evaluate the class probabilities and class-specific contributions at the current estimates, denoted \(\pi_{ig}^{(r)}\) and \(L_{ig}^{(r)}\). Bayes’ rule gives the conditional probability that subject \(i\) is in class \(g\), given the observed data and current estimates. This probability, also called the responsibility, is

For each subject, \(0\leq \tau_{ig}^{(r)}\leq 1\) and \(\sum_g\tau_{ig}^{(r)}=1\). The E-step does not change the model parameters. It computes the expected latent class indicators under the current estimates.

The M-Step: Maximize the Expected Complete-Data Log Likelihood

If class membership were observed, the complete-data log likelihood would be

The EM algorithm replaces each unobserved \(Z_{ig}\) by its conditional expectation from the E-step. The resulting function to maximize is

When the parameter blocks are separate as written, the M-step maximizes \(Q\) with respect to class-membership, longitudinal, and survival parameters. Each block uses the same responsibilities, but a different weighted log likelihood.

Updates for Class Proportions and Longitudinal Parameters

If class probabilities are constant across subjects, write them as \(\pi_g\), where \(\sum_g\pi_g=1\). The class-proportion part of \(Q\) is \(\sum_g(\sum_i\tau_{ig})\log\pi_g\). Maximizing subject to the sum-to-one constraint gives

Thus the updated proportion is the average responsibility for that class, not the fraction of subjects assigned to it by a hard rule. If class membership is modeled using baseline covariates, for example a multinomial logistic model, the class-membership M-step is a responsibility-weighted multinomial regression. Its coefficients generally require numerical maximization rather than this closed-form update.

For the longitudinal component, the M-step maximizes \(\sum_i\tau_{ig}^{(r)}\log f_g(y_i\mid x_i;\theta_g)\) for each class. Consequently, it is a weighted longitudinal fit, with subject \(i\) receiving weight \(\tau_{ig}^{(r)}\) in class \(g\). The precise update depends on the chosen longitudinal model. For example, if each subject contributes one Gaussian observation \(y_i\) with class-specific mean \(\mu_g\) and variance \(\sigma_g^2\), the updates are the weighted mean and weighted variance:

Updating Survival Parameters

The survival M-step maximizes \(\sum_i\tau_{ig}^{(r)}\log h_g(t_i,\delta_i\mid x_i;\psi_g)\) for each class. For a parametric model, this is the usual survival log likelihood with fractional subject weights. One useful illustration is a class-specific exponential model with constant hazard \(\lambda_g\), for which \(h_g(t_i,\delta_i)=\lambda_g^{\delta_i}\exp(-\lambda_g t_i)\). The weighted log likelihood is

Setting its derivative with respect to \(\lambda_g\) to zero gives the weighted event-count divided by weighted exposure-time update:

This requires positive weighted exposure time. More generally, a chosen parametric or semiparametric survival specification determines the form of the weighted survival maximization; an update for an exponential hazard should not be used for a different hazard model without derivation.

Why the Algorithm Increases the Likelihood

The EM likelihood argument follows from a lower bound. For one subject, put \(a_{ig}=\pi_{ig}L_{ig}\), so that the observed log-likelihood contribution is \(\log\sum_g a_{ig}\). For any probabilities \(q_{ig}\) that sum to 1 and any nonnegative candidate contributions \(a_{ig}\), write

The first inequality holds because it retains only terms with \(q_{ig}>0\); Jensen’s inequality for the concave logarithm gives the second inequality. A positive candidate contribution with \(q_{ig}=0\) is simply omitted, so the bound remains valid. Taking \(q_{ig}=\tau_{ig}^{(r)}\) gives equality at the current parameter values: every positive current contribution has positive posterior responsibility, and the Jensen bound is tight. The lower bound, summed across subjects, differs from \(Q(\Theta\mid\Theta^{(r)})\) only by terms that do not depend on \(\Theta\). Maximizing \(Q\) therefore cannot decrease the observed-data log likelihood. An exact M-step increases or preserves it; a generalized EM step has the same guarantee if it increases \(Q\). This is a monotonicity result, not a guarantee of reaching the global maximum.

Diagnostics, Alternatives, and Limitations

Track the observed-data log likelihood after each full iteration. It should not decrease beyond small numerical rounding error when the E- and M-steps are implemented correctly. Also monitor parameter changes and the largest changes in responsibilities. Use stable log-sum-exp calculations for the denominator in the responsibility formula when likelihood contributions are extremely small; direct multiplication can underflow to zero. Check that responsibilities sum to one and that the M-step uses the same weights consistently across all parameter blocks.

Convergence of the log likelihood alone is not enough to establish a trustworthy fit. A nearly empty class can have a very small total responsibility, making its parameter estimates unstable or poorly determined. Different initial values can lead to different stationary solutions. Run multiple sensible initializations, compare attained log likelihoods, and inspect whether estimates and fitted contributions are numerically plausible. If the likelihood decreases, first check coding, weight normalization, and whether each M-step truly maximizes or increases its weighted objective.

EM is a computational strategy, not a remedy for model misspecification. Its weighted updates rely on the stated class-specific longitudinal and survival likelihoods, class-membership model, and censoring assumptions. It may converge slowly, and monotone improvement does not ensure global optimality. The suitability of the number and form of classes, as well as the substantive interpretation of fitted classes, requires separate assessment.

Common Mistakes

• Treating responsibilities as observed class labels. They are conditional probabilities under the current estimates, and generally lie between 0 and 1.
• Updating parameters with hard assignments. This changes the procedure and discards uncertainty; the EM M-step uses \(\tau_{ig}^{(r)}\) as weights.
• Using the new parameters in the same E-step. Compute responsibilities from the current iteration’s parameters, hold them fixed through the M-step, and then recalculate.
• Assuming every M-step has a closed form. Covariate-dependent class probabilities and many longitudinal or survival models require numerical weighted maximization.
• Calling a nondecreasing likelihood proof a global-optimum proof. The guarantee is monotonicity, not global convergence.