M-Step Updates for Longitudinal Fixed Effects

Why Posterior Random-Effect Means Enter the Update

In a longitudinal mixed model, each subject’s repeated measurements share latent random effects. The EM algorithm treats those effects as missing data: the E-step averages the complete-data log likelihood over their posterior distribution, and the M-step maximizes that average. For the Gaussian longitudinal submodel, the fixed-effects update is a least-squares calculation with a specific adjustment: the posterior mean of each subject’s random effects is subtracted from that subject’s response vector.

This is not the same as fitting a regression after substituting estimated random effects as if they were known. In the M-step, the posterior distribution is fixed at the current EM iterate, and the update maximizes the expected complete-data criterion. The distinction matters because posterior uncertainty contributes to other parameter updates, even though it does not change the fixed-effects formula under the assumptions below.

Definition: For subject $i=1,\ldots,n$, let $y_i$ be the vector of $n_i$ longitudinal measurements, $X_i$ the fixed-effects design matrix, and $Z_i$ the random-effects design matrix. The Gaussian mixed submodel is $$ y_i\mid b_i \sim N(X_i\beta+Z_i b_i,\ R_i), $$ where $\beta$ is the fixed-effects vector, $b_i$ is the latent random-effects vector, and $R_i$ is the residual covariance matrix. At EM iteration $k$, write $m_i^{(k)}=E(b_i\mid \text{observed data};\theta^{(k)})$ for the E-step posterior mean. The conditioning includes the observed event-time information when it is part of the joint model.

For the update derived here, assume subjects are conditionally independent given their random effects, $R_i$ is positive definite and does not depend on $\beta$, and the posterior distribution in the E-step is evaluated at the current parameter vector $\theta^{(k)}$. The design matrices and the posterior moments are therefore fixed while maximizing with respect to $\beta$. If $R_i=\sigma^2 I$ for every subject, the common residual variance factor does not affect the maximizer.

Deriving the Fixed-Effects Update

Ignoring terms that do not depend on $\beta$, the Gaussian longitudinal log likelihood conditional on the random effects is $$ \ell_L(\beta;b)= -\frac{1}{2}\sum_{i=1}^n (y_i-X_i\beta-Z_i b_i)^\mathsf{T}R_i^{-1}(y_i-X_i\beta-Z_i b_i). $$ The longitudinal contribution to the EM auxiliary function is the conditional expectation of this expression under the E-step posterior. Define $r_i=y_i-Z_i b_i$. Expanding the quadratic shows that its $\beta$-dependent terms are $2\beta^\mathsf{T}X_i^\mathsf{T}R_i^{-1}r_i-\beta^\mathsf{T}X_i^\mathsf{T}R_i^{-1}X_i\beta$. Taking the posterior expectation replaces $b_i$ in the linear term by $m_i^{(k)}$; the quadratic term in $\beta$ is unchanged.

Q_L(\beta\mid\theta^{(k)}) = C + \sum_{i=1}^n \left[ \beta^\mathsf{T}X_i^\mathsf{T}R_i^{-1}(y_i-Z_i m_i^{(k)}) -\frac{1}{2}\beta^\mathsf{T}X_i^\mathsf{T}R_i^{-1}X_i\beta \right],

Here $C$ does not depend on $\beta$. Differentiating with respect to $\beta$ and setting the gradient to zero gives the normal equations. If the matrix on the left is nonsingular, the unique maximizer is:

\left(\sum_{i=1}^n X_i^\mathsf{T}R_i^{-1}X_i\right)\beta^{(k+1)} = \sum_{i=1}^n X_i^\mathsf{T}R_i^{-1}(y_i-Z_i m_i^{(k)}), $$ $$ \beta^{(k+1)} = \left(\sum_{i=1}^n X_i^\mathsf{T}R_i^{-1}X_i\right)^{-1} \left(\sum_{i=1}^n X_i^\mathsf{T}R_i^{-1}(y_i-Z_i m_i^{(k)})\right).

The matrix is nonsingular when the fixed-effects design has full column rank after accounting for the positive-definite weights. The Hessian of $Q_L$ with respect to $\beta$ is $-\sum_i X_i^\mathsf{T}R_i^{-1}X_i$, which is negative definite under this condition. Thus the solution is a maximum, not merely a stationary point.

Key result: In the Gaussian longitudinal submodel, the fixed-effects M-step is generalized least squares of $y_i-Z_i m_i^{(k)}$ on $X_i$, using residual precision matrices $R_i^{-1}$. For $R_i=\sigma^2 I$, this reduces to ordinary least squares on the responses adjusted by the E-step posterior random-effect means.

The derivation also explains why the posterior covariance of $b_i$ does not appear in this particular update. The expected quadratic includes a term involving $E(b_i b_i^\mathsf{T}\mid\text{data})$, but that term is independent of $\beta$, so it disappears when differentiating with respect to $\beta$. This does not mean that posterior uncertainty is irrelevant to EM; it generally matters for updating covariance parameters.

Worked Calculations

Worked Example: A Random-Intercept Model With Equal Residual Variance

Suppose the hypothetical longitudinal model is $y_{ij}=\beta+b_i+\varepsilon_{ij}$, with a common residual variance and a random intercept $b_i$. At the current E-step, suppose the posterior means are $m_1^{(k)}=0.5$ and $m_2^{(k)}=-1$. The observed measurements are $4,6$ for subject 1 and $2,5$ for subject 2. The design for $\beta$ is an intercept, so the update is the mean of the four adjusted measurements.

For subject 1, subtract the posterior mean $0.5$: the adjusted values are $4-0.5=3.5$ and $6-0.5=5.5$. For subject 2, subtract $-1$: the adjusted values are $2-(-1)=3$ and $5-(-1)=6$. Therefore, $$ \beta^{(k+1)}=\frac{3.5+5.5+3+6}{4}=\frac{18}{4}=4.5. $$

The update is $4.5$ response units. The calculation uses the posterior means supplied by the E-step, not the unobserved true random intercepts.

Worked Example: A Weighted Intercept Update

Consider three hypothetical subjects with one longitudinal measurement each and a random intercept. Let the observations be $y=(5,3,8)$, the E-step posterior means be $m^{(k)}=(1,-0.5,0.5)$, and the residual variances be $1,4,1$, respectively. Thus $R_i$ is the scalar residual variance for subject $i$, and the fixed-effects design is an intercept.

The adjusted responses are $5-1=4$, $3-(-0.5)=3.5$, and $8-0.5=7.5$. Generalized least squares weights these values by the inverse residual variances, $1$, $1/4$, and $1$. Hence, $$ \beta^{(k+1)} =\frac{1(4)+(1/4)(3.5)+1(7.5)}{1+1/4+1} =\frac{4+0.875+7.5}{2.25} =\frac{12.375}{2.25} =5.5. $$

The updated intercept is $5.5$ response units. The middle subject receives less weight because its residual variance is four times as large as the other two residual variances.

Worked Example: Updating an Intercept and Time Slope

Suppose a hypothetical model has a fixed intercept and a fixed linear time slope, as well as a subject-specific random intercept. The four measurements are: subject A has $y=2,4$ at times $0,1$, and subject B has $y=1,4$ at times $0,1$. The E-step posterior means of the random intercepts are $0.5$ for A and $-0.5$ for B. Assume a common residual variance.

Subtracting each posterior mean gives adjusted responses $1.5,3.5$ for A and $1.5,4.5$ for B. With columns $1$ and time, the design matrix has rows $(1,0),(1,1),(1,0),(1,1)$. The normal-equation components are $$ X^\mathsf{T}X= \begin{pmatrix}4&2\\2&2\end{pmatrix}, \qquad X^\mathsf{T}(y-Zm^{(k)})= \begin{pmatrix}11\\8\end{pmatrix}. $$

Writing the fixed-effects update as $(\beta_0,\beta_1)^\mathsf{T}$, the equations are $4\beta_0+2\beta_1=11$ and $2\beta_0+2\beta_1=8$. Subtracting the second equation from the first gives $2\beta_0=3$, so $\beta_0=1.5$. Substitution gives $2(1.5)+2\beta_1=8$, hence $\beta_1=2.5$. The update is an intercept of $1.5$ response units and a slope of $2.5$ response units per time unit.

Scope, Diagnostics, and Alternatives

These equations update the fixed effects in the Gaussian longitudinal contribution. They are the full fixed-effects M-step when the other complete-data likelihood terms do not depend directly on $\beta$, as in a shared-random-effects specification whose survival component depends on $b_i$ but not directly on $\beta$. In a joint model where the event hazard depends on the longitudinal trajectory $X_i(t)\beta+Z_i(t)b_i$, the survival contribution can also depend on $\beta$. In that case, the displayed generalized least-squares solution maximizes the longitudinal contribution only; it is not, by itself, the maximizer of the full auxiliary function. The full update must account for the additional survival term.

A practical implementation should check the dimensions of $X_i$, $Z_i$, and $R_i$, verify that each $R_i$ is positive definite, and confirm that the fixed-effects information matrix $\sum_i X_i^\mathsf{T}R_i^{-1}X_i$ has full rank. A singular or nearly singular matrix may indicate redundant fixed-effect columns, insufficient variation in time or covariates, or severe numerical scaling problems. Large changes in the update can also signal unstable posterior means or poorly conditioned design matrices; inspect the E-step calculations and the conditioning of the normal equations rather than silently accepting an implausible estimate.

The main alternative within this Gaussian calculation is ordinary least squares versus generalized least squares. Ordinary least squares is appropriate when residual variances are common and the within-subject residual covariance is proportional to the identity; generalized least squares uses $R_i^{-1}$ to account for specified heteroscedasticity or residual correlation. Neither alternative justifies treating the latent $b_i$ as observed. The E-step posterior mean is the correct quantity for this fixed-effects update under the stated model, while an individual conditional mode or a plug-in estimate is not generally the same quantity.

The update also relies on a linear mean $X_i\beta+Z_i b_i$ and residual covariance that does not depend on $\beta$. Nonlinear mean structures, parameter-dependent residual covariance, or survival terms involving $\beta$ can prevent a closed-form update. In such cases, the expected complete-data criterion still defines the M-step target, but numerical maximization may be needed.

Common pitfall: Subtracting the posterior random-effect mean is not an optional correction, and one must not subtract the posterior variance as though it were another response adjustment. The mean enters the linear term in $\beta$; the posterior second moment contributes only terms that do not involve $\beta$ in this update. Also, the inverse residual variance is the weight: a variance of $4$ gives weight $1/4$, not $4$.

Check Your Understanding

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

In an equal-variance random-intercept model, one subject has measurements $6$ and $8$, and the current posterior mean of the intercept is $2$. What adjusted measurements enter the fixed-intercept update?
For one observation with residual variance $9$, what generalized least-squares weight does it receive?
Why does the posterior covariance of the random effects not appear in the Gaussian fixed-effects update derived here?
If the survival hazard depends directly on $X_i(t)\beta+Z_i(t)b_i$, is the longitudinal generalized least-squares update necessarily the full M-step? Explain.

Solutions. (1) The adjusted values are $6-2=4$ and $8-2=6$. (2) The weight is the inverse variance, $1/9$. (3) Its contribution to the expected complete-data quadratic is independent of $\beta$; differentiating $Q_L$ with respect to $\beta$ removes that term, leaving the posterior mean in the normal equations. (4) No. The survival contribution then also depends on $\beta$, so maximizing the longitudinal term alone need not maximize the full auxiliary function.

M-Step Updates for Longitudinal Fixed Effects

What You'll Learn