EM Estimation for Shared-Random-Effects Models

Why use EM for shared-random-effects models?

In a shared-parameter model, a latent subject-specific effect helps describe both the outcome history and the dropout process. The observed-data likelihood averages over that unobserved effect, which can make direct maximization difficult. The expectation-maximization (EM) algorithm instead alternates between calculating the posterior distribution of the latent effects under the current parameters and maximizing an expected complete-data log-likelihood.

The algorithm does not make dropout ignorable. The dropout contribution remains in the likelihood and therefore in the posterior distribution calculated at each E-step. Nor does EM create information about parameters that the observed data cannot identify: the identifiability conditions discussed in the earlier material on shared effects still apply.

Model and complete-data likelihood

Consider \(n\) independent subjects. Subject \(i\) has a \(q\)-dimensional random effect \(b_i\), with \(b_i\sim N_q(0,D)\). At observed outcome times \(j\in\mathcal O_i\), suppose \(Y_{ij}=x_{ij}^{T}\beta+z_{ij}^{T}b_i+\varepsilon_{ij}\), where the errors are independent \(N(0,\sigma^2)\) conditional on \(b_i\). For a discrete-time dropout process, let \(d_{it}=1\) indicate dropout in interval \(t\), and let \(r_{it}=1\) indicate the subject is at risk at the start of that interval. Use \(\operatorname{logit}(h_{it})=w_{it}^{T}\gamma+a_{it}^{T}b_i\), where \(h_{it}=P(d_{it}=1\mid r_{it}=1,b_i)\). The dropout contribution is included only for intervals in which \(r_{it}=1\).

For example, subject \(i\)'s complete-data log-likelihood, omitting constants that do not depend on \(\theta\), is

Here \(m_i\) is the number of observed outcomes for subject \(i\). The term \(-q\log(2\pi)/2\), as well as the normal-outcome constant, is omitted. A different outcome or dropout model changes the corresponding likelihood terms, but not the EM principle.

The E-step and M-step

At iteration \(k\), the E-step uses the current estimate \(\theta^{(k)}\) to form each subject's posterior distribution of \(b_i\), conditional on that subject's observed outcomes and dropout history. In this model its density is proportional to the product of the observed-outcome density, the dropout-history likelihood, and the \(N_q(0,D^{(k)})\) density. The E-step computes expectations under that posterior; common required quantities include \(E[b_i\mid\mathrm{obs}_i]\), \(E[b_ib_i^{T}\mid\mathrm{obs}_i]\), and expectations of the dropout log-likelihood terms.

In the M-step, maximize \(Q(\theta\mid\theta^{(k)})=\sum_i E_{\theta^{(k)}}[\ell_{c,i}(\theta)\mid\mathrm{obs}_i]\). For the Gaussian outcome component, the regression update has a closed form when the design matrices are fixed. The residual variance update uses the posterior expectation of the squared residual, not just the squared residual evaluated at the posterior mean. The covariance update uses posterior second moments. For a logistic dropout model, the expected log-likelihood generally has no closed-form maximizing coefficient, so numerical maximization is used.

Why the EM step increases likelihood

Let \(Z\) denote all the latent random effects, and let \(p_\theta(y_{\mathrm{obs}},Z)\) be the joint density of the observed data and latent effects. Set \(q(Z)=p_{\theta^{(k)}}(Z\mid y_{\mathrm{obs}})\). For any candidate \(\theta\), writing the observed likelihood as an integral against \(q\) and applying Jensen's inequality gives

The second expectation on the right does not depend on the candidate \(\theta\). The first is \(Q(\theta\mid\theta^{(k)})\). At \(\theta=\theta^{(k)}\), the inequality is an equality because \(q\) is that parameter value's posterior. Consequently, an exact M-step that does not decrease \(Q\) cannot decrease the observed-data log-likelihood. With approximate integration or an incomplete numerical maximization, this guarantee may not hold exactly, so the observed likelihood should be monitored.

Worked calculations

General M-step formulas and practical checks

For the Gaussian outcome model, stack the observed outcomes for subject \(i\) into \(Y_i\), and let \(X_i\) and \(Z_i\) be the corresponding fixed- and random-effect design matrices. If \(\sum_i X_i^{T}X_i\) is invertible, the update is

The denominator in the variance update is the number of observed outcome measurements under this likelihood, not the number of subjects. If \(D\) is unrestricted and the random-effects mean is fixed at zero, its update is

This update uses second moments, including uncertainty about the latent effects. Constraints such as a diagonal \(D\) require the corresponding constrained maximization rather than blindly retaining all estimated off-diagonal entries.

For the dropout coefficients, the M-step maximizes \(\sum_i E[\sum_{t:r_{it}=1}\{d_{it}\eta_{it}-\log(1+\exp(\eta_{it}))\}\mid\mathrm{obs}_i]\). The expectation is over each subject's posterior distribution of \(b_i\). It can be evaluated by quadrature or another integration method, then maximized numerically. The dropout update is not an ordinary logistic regression on estimated random effects: treating posterior means as known covariates discards posterior uncertainty and does not, in general, maximize this expected likelihood.

A practical E-step often needs numerical integration because the posterior combines the random-effects density with outcome and dropout terms. The earlier material on maximum likelihood estimation with adaptive quadrature develops a direct-likelihood integration approach; in EM, integration instead supplies posterior expectations for \(Q\). Check that increasing the integration accuracy does not materially change those expectations or the resulting updates. Also verify that probabilities are valid, covariance estimates are positive semidefinite, and the observed-data log-likelihood does not fall materially after an iteration. A substantial decrease can signal inaccurate integration, a failed M-step optimization, or an implementation error.

Comparison, limitations, and common mistakes

EM is useful when complete-data terms have convenient updates and posterior moments can be calculated reliably. Its limitations include the need for repeated E-steps, potential sensitivity to starting values in a non-concave problem, and slow progress in some regions of the parameter space. An increase in likelihood does not establish that the model fits well, that parameters are identifiable, or that the shared-effects assumptions are scientifically appropriate. Those questions require separate assessment.

Common mistake: “I can update the variance using \((Y_i-X_i^T\beta-Z_i^TE[b_i\mid\mathrm{obs}_i])^2\), because the posterior mean is the best estimate of \(b_i\).” This replaces a latent random effect by a point estimate. The correct EM update averages the squared residual over its full posterior distribution: \(E[(Y_i-X_i^T\beta-Z_i^Tb_i)^2\mid\mathrm{obs}_i]\). The difference includes posterior variance and can materially affect \(\sigma^2\) and \(D\).