The Multimarker E-Step In EM Estimation

Why the E-Step Is Subject-Specific

In a shared-random-effects multimarker joint model, a subject’s markers and event time provide joint information about that subject’s latent trajectory. The E-step must therefore average over the random effects conditional on both the observed marker data and the observed event or censoring time. Using marker data alone would omit information carried by the survival outcome; using survival data alone would discard the observed longitudinal measurements.

This tutorial focuses on deriving and calculating those subject-level conditional expectations. The integration is over the full vector of random effects shared across markers, not over one independently fitted random effect per marker. The earlier material on quadrature for multimarker random-effects integrals supplies numerical integration methods; here the emphasis is on the distribution being integrated and on the quantities the EM algorithm needs.

Definition: For subject \(i\), let \(Y_i\) denote all observed longitudinal measurements across markers, \(T_i\) the observed event or censoring time, and \(\delta_i\) the event indicator, where \(\delta_i=1\) denotes an event and \(\delta_i=0\) denotes censoring. Let \(b_i\) be the stacked vector of that subject’s shared random effects. At parameter value \(\theta^{(r)}\), the E-step distribution is the conditional distribution of \(b_i\) given \((Y_i,T_i,\delta_i)\), calculated under \(\theta^{(r)}\). Its conditional mean is \(m_i^{(r)}=E[b_i\mid Y_i,T_i,\delta_i;\theta^{(r)}]\), and its conditional raw second moment is \(S_i^{(r)}=E[b_i b_i^\mathsf{T}\mid Y_i,T_i,\delta_i;\theta^{(r)}]\). The conditional covariance is \(S_i^{(r)}-m_i^{(r)}m_i^{(r)\mathsf{T}}\), not \(S_i^{(r)}\).

Model Contributions to the Conditional Distribution

For clarity, write a general marker model as

Y_i \mid b_i \sim f_Y(Y_i\mid b_i;\theta), \qquad b_i\sim N(0,D).

Here \(Y_i\) stacks the measurements from all markers, and the design matrices and fixed effects may differ by marker. The marker-specific residual model is part of \(f_Y\). A common specification assumes conditionally independent measurement errors across markers given \(b_i\); other residual covariance structures can be used, but must be included in the specified \(f_Y\). The joint model’s hazard, conditional on the random effects and covariates, is written \(h_i(t\mid b_i;\theta)\). The conditional survival contribution through \(T_i\) is

f_T(T_i,\delta_i\mid b_i;\theta) = \{h_i(T_i\mid b_i;\theta)\}^{\delta_i} \exp\left\{-\int_0^{T_i}h_i(u\mid b_i;\theta)\,du\right\}.

This expression assumes the event-time model is correctly specified and that censoring is noninformative conditional on the modeled quantities. It accommodates an event and right censoring: when \(\delta_i=0\), there is no hazard factor at \(T_i\), but the survival probability through \(T_i\) remains. The longitudinal and event-time observations are conditionally independent given \(b_i\) and the modeled covariates, as encoded by their product in the joint likelihood.

Let \(\phi(b_i;0,D)\) denote the multivariate normal random-effects density. Bayes’ rule gives the E-step density, up to its normalizing constant:

p_i^{(r)}(b) = \frac{ f_Y(Y_i\mid b;\theta^{(r)}) f_T(T_i,\delta_i\mid b;\theta^{(r)}) \phi(b;0,D^{(r)}) }{ \int f_Y(Y_i\mid u;\theta^{(r)}) f_T(T_i,\delta_i\mid u;\theta^{(r)}) \phi(u;0,D^{(r)})\,du }.

The denominator is the subject-specific observed-data likelihood contribution. It normalizes the density over the entire stacked random-effects vector. This formula also makes clear why the posterior distribution generally is not Gaussian: even if the prior is Gaussian and the marker model is Gaussian, the survival likelihood can be a nonlinear function of \(b\).

What the E-Step Must Average

The EM algorithm forms the expected complete-data log likelihood under the current conditional distributions. For independent subjects, its E-step objective is

Q(\theta\mid\theta^{(r)}) = \sum_i E_{p_i^{(r)}}\left[ \log f_Y(Y_i\mid b_i;\theta) +\log f_T(T_i,\delta_i\mid b_i;\theta) +\log\phi(b_i;0,D) \right].

The expectation is taken using the conditional distribution at \(\theta^{(r)}\); the parameter inside the complete-data log likelihood is the candidate \(\theta\) being maximized. For Gaussian random-effects terms, the update for \(D\) depends on \(S_i^{(r)}\), the conditional raw second moment. The longitudinal part may also require conditional means and second moments. But the event-time part often contains an exponential of a function of \(b_i\), through the hazard and cumulative hazard. In that case, the event-time expectation generally cannot be recovered from \(m_i^{(r)}\) and \(S_i^{(r)}\) alone: it must be integrated over the conditional distribution, or approximated using a suitable numerical method.

Calculating the Moments by Normalized Quadrature

Suppose numerical integration evaluates the subject’s unnormalized posterior contribution at nodes \(b_{i1},\ldots,b_{iM}\). Define

q_{im} = f_Y(Y_i\mid b_{im};\theta^{(r)}) f_T(T_i,\delta_i\mid b_{im};\theta^{(r)}) \phi(b_{im};0,D^{(r)})\,a_{im}, \qquad w_{im}=\frac{q_{im}}{\sum_{\ell=1}^M q_{i\ell}},

where \(a_{im}\) is the quadrature weight, including any required transformation factor. The normalized weights are nonnegative and sum to one. They approximate the posterior probabilities associated with the quadrature nodes. The corresponding moment approximations are

m_i^{(r)}\approx\sum_{m=1}^M w_{im}b_{im}, \qquad S_i^{(r)}\approx\sum_{m=1}^M w_{im}b_{im}b_{im}^{\mathsf T}, \qquad \operatorname{Cov}(b_i\mid Y_i,T_i,\delta_i)\approx S_i^{(r)}-m_i^{(r)}m_i^{(r)\mathsf T}.

These are posterior-weighted averages, not averages under the original random-effects distribution. In particular, the survival contribution can shift the conditional mean and change the conditional covariance. Integrating in the full random-effects dimension can be computationally demanding, so node accuracy and stability matter.

Worked Example: Conditional Mean and Covariance From Three Nodes

Suppose a two-dimensional random-effects vector has three quadrature nodes with normalized posterior weights \(0.2\), \(0.5\), and \(0.3\), at \(b_1=(-1,0)^\mathsf T\), \(b_2=(0,1)^\mathsf T\), and \(b_3=(1,0)^\mathsf T\). These values are hypothetical and illustrate the moment calculation.

The conditional mean is the weighted node average:

m = 0.2\begin{pmatrix}-1\\0\end{pmatrix} +0.5\begin{pmatrix}0\\1\end{pmatrix} +0.3\begin{pmatrix}1\\0\end{pmatrix} = \begin{pmatrix}0.1\\0.5\end{pmatrix}.

For the raw second moment, first take each node’s outer product and then weight it:

S = 0.2\begin{pmatrix}1&0\\0&0\end{pmatrix} +0.5\begin{pmatrix}0&0\\0&1\end{pmatrix} +0.3\begin{pmatrix}1&0\\0&0\end{pmatrix} = \begin{pmatrix}0.5&0\\0&0.5\end{pmatrix}.

Subtracting the mean’s outer product gives the conditional covariance:

S-mm^\mathsf T = \begin{pmatrix}0.5&0\\0&0.5\end{pmatrix} - \begin{pmatrix}0.01&0.05\\0.05&0.25\end{pmatrix} = \begin{pmatrix}0.49&-0.05\\-0.05&0.25\end{pmatrix}.

The off-diagonal covariance is \(-0.05\), even though the raw second-moment off-diagonal entry is zero. The distinction matters: the raw second moment is the quantity used in the unrestricted covariance update below, while the covariance describes uncertainty around the conditional mean.

Why the Event Observation Changes the Weights

A survival outcome changes the E-step distribution through both the hazard at an observed event time and the cumulative hazard accumulated up to the observed time. The following discrete example shows the calculation explicitly.

Worked Example: An Event Updates a Two-Point Random-Effects Distribution

Let the scalar random effect take values \(b=-1\) and \(b=1\), each with prior probability \(0.5\). Suppose the hypothetical longitudinal likelihoods at those values are \(f_Y(Y\mid -1)=0.5\) and \(f_Y(Y\mid 1)=1\). For a subject with an event at \(T=2\), let the hazard be \(h(t\mid b)=0.2\exp(0.5b)\), constant over time. The event-time contribution is \(h(2\mid b)\exp\{-2h(2\mid b)\}\).

At \(b=-1\), the hazard is \(0.2e^{-0.5}\approx0.121306\), so the event-time contribution is

0.121306\exp\{-2(0.121306)\}\approx0.0952.

At \(b=1\), the hazard is \(0.2e^{0.5}\approx0.329744\), and the event-time contribution is

0.2e^{0.5}\exp\{-2(0.2e^{0.5})\} \approx 0.329744\exp(-0.659489) \approx 0.1705.

Multiplying the event contributions by the longitudinal likelihoods gives approximately \(0.0476\) for \(b=-1\) and \(0.1705\) for \(b=1\). The equal prior probabilities cancel when normalizing. Thus

w_-=\frac{0.0476}{0.0476+0.1705}\approx0.2182, \qquad w_+=\frac{0.1705}{0.0476+0.1705}\approx0.7818.

The conditional mean is \(E[b\mid Y,T,\delta=1]=(-1)w_-+(1)w_+\approx0.5636\). Since \(b^2=1\) at either node, the conditional raw second moment is \(1\), and the conditional variance is \(1-(0.5636)^2\approx0.6824\). This result is specific to the stated longitudinal likelihoods and hazard. It demonstrates how the joint evidence reweights the two possible random effects; it is not the mean under the prior distribution.

A Covariance Update From the Raw Second Moments

If the random effects are independent across subjects with common distribution \(N(0,D)\), and \(D\) is unrestricted, the part of \(Q\) involving \(D\) is

Q_D = -\frac{n}{2}\log|D| -\frac{1}{2}\sum_i \operatorname{tr}\left(D^{-1}S_i^{(r)}\right) +\text{constant}.

To see why the update uses raw second moments, the Gaussian log density contains \(b_i^\mathsf T D^{-1}b_i\). Its conditional expectation is

E[b_i^\mathsf T D^{-1}b_i\mid Y_i,T_i,\delta_i] = \operatorname{tr}\left(D^{-1}E[b_i b_i^\mathsf T\mid Y_i,T_i,\delta_i]\right) = \operatorname{tr}(D^{-1}S_i^{(r)}).

Differentiating \(Q_D\) with respect to \(D\) gives \(-\tfrac{n}{2}D^{-1}+\tfrac{1}{2}D^{-1}(\sum_i S_i^{(r)})D^{-1}\). Setting this equal to zero and multiplying on both sides by \(D\) yields the maximizer

D^{(r+1)}=\frac{1}{n}\sum_i S_i^{(r)}.

This update presumes a zero-mean Gaussian random-effects distribution and no additional constraints on \(D\). If the covariance structure is restricted, such as by tying entries or imposing a structured form, the M-step must maximize \(Q_D\) subject to that structure instead.

Worked Example: Updating a Two-Marker Random-Effects Covariance

For two hypothetical subjects, suppose the E-step returns

S_1= \begin{pmatrix}0.5&0.6\\0.6&1.7\end{pmatrix}, \qquad S_2= \begin{pmatrix}1.4&-0.2\\-0.2&0.8\end{pmatrix}.

With an unrestricted covariance and \(n=2\), average the two raw second-moment matrices entry by entry:

D^{(r+1)} = \frac{1}{2}(S_1+S_2) = \begin{pmatrix}0.95&0.20\\0.20&1.25\end{pmatrix}.

The diagonal entries are the updated random-effects variances, and \(0.20\) is the updated covariance in the units of the two random effects. The matrix is positive definite because its first diagonal entry is positive and its determinant is \(0.95(1.25)-(0.20)^2=1.1475>0\). Averaging conditional covariances alone would generally give the wrong update: it would omit variation in the subject-specific conditional means.

Diagnostics, Alternatives, and Limitations

Numerical accuracy in the E-step affects every subsequent M-step. Useful checks include verifying that normalized posterior weights sum to one, that computed covariances are symmetric and positive semidefinite up to numerical tolerance, and that the observed-data log likelihood behaves consistently with the EM iterations. A covariance with a materially negative eigenvalue, or a large change in moments when integration settings are modestly refined, signals a numerical or implementation problem.

Posterior weights concentrated on a few nodes can indicate that the current quadrature rule does not adequately resolve the conditional distribution. Increasing the number or adapting the placement of nodes can help, but high random-effects dimension makes tensor-product quadrature expensive. Computation on the log scale is often important: when likelihood contributions are very small, subtracting the largest log weight before exponentiating helps prevent underflow without changing normalized weights.

Quadrature approximates the defining integrals directly and can be effective when the dimension is manageable and the nodes capture the posterior mass. A Laplace approximation instead uses a local Gaussian approximation around a posterior mode; it can be less costly in higher dimensions but may be inaccurate for skewed or multimodal conditional distributions. Monte Carlo EM approximates the conditional expectations with posterior draws and is useful in more difficult dimensions, but introduces simulation error that must be controlled. None of these methods removes the need to use the joint longitudinal and event-time contribution in the conditional distribution.

A common mistake is to call \(E[b_i b_i^\mathsf T\mid\text{data}]\) the conditional covariance, then use it as though it measured variability around the posterior mean. It is the raw second moment; the covariance is obtained by subtracting the mean’s outer product. A related mistake is to compute \(E[b_i\mid Y_i]\) and treat it as the joint-model E-step mean. The correct distribution conditions on \(Y_i,T_i,\delta_i\), so the event and censoring information must enter through the event-time likelihood contribution.

The E-step does not make random effects observed or eliminate modeling assumptions. Its moments are conditional on the fitted model, including its marker residual structure, hazard specification, and censoring assumptions. Poorly identified marker associations or a weakly supported covariance structure can still produce unstable estimates, even when the conditional integration is implemented correctly.

Check Your Understanding

Answer each question, then compare your reasoning with the solutions.

In a two-dimensional model, \(m=(0.2,-0.1)^\mathsf T\) and \(S=\begin{pmatrix}0.50&0.08\\0.08&0.30\end{pmatrix}\). What is the conditional covariance?
For a censored subject with \(\delta_i=0\), which parts of the event-time likelihood remain in the E-step contribution?
Why does the unrestricted update for \(D\) average \(S_i\) rather than the conditional covariance matrices?
Suppose normalized quadrature weights do not sum to one because the unnormalized values are \(0.2,0.3,0.5\). What should be done before calculating a conditional mean?

Solutions. (1) Subtract \(mm^\mathsf T=\begin{pmatrix}0.04&-0.02\\-0.02&0.01\end{pmatrix}\) from \(S\). The conditional covariance is \(\begin{pmatrix}0.46&0.10\\0.10&0.29\end{pmatrix}\). (2) The hazard-at-\(T_i\) factor is absent, but the survival factor \(\exp\{-\int_0^{T_i}h_i(u\mid b_i)\,du\}\) remains and weights the random effects. (3) The Gaussian quadratic term involves \(E[b_i b_i^\mathsf T\mid\text{data}]=S_i\), which includes both conditional covariance and the outer product of the conditional mean. (4) Normalize the values by their sum, \(1.0\), giving weights \(0.2,0.3,0.5\); then calculate the weighted mean. In general, divide every nonnegative unnormalized weight by their total before computing posterior expectations.