The M-Step for Population Parameters

What the M-Step Optimizes

In an NLME model, the latent individual effects make direct likelihood calculations difficult because they must be integrated out. The EM approach instead treats those effects as missing data: the E-step computes conditional expectations under the current parameter values, and the M-step uses those expectations to update population parameters. This tutorial focuses on that maximization: which updates can be derived explicitly, which usually require numerical optimization, and how to verify that an update has actually maximized the intended objective.

Definition: Let $y_i$ denote the observations for subject $i$, and let $\eta_i$ be that subject’s latent random-effects vector. Write $\phi$ for all population parameters, including fixed effects, random-effect covariance parameters, and residual parameters. Given the current value $\phi^{(r)}$, the M-step maximizes $$ Q(\phi\mid\phi^{(r)})= \sum_{i=1}^{N} E_{\eta_i\mid y_i,\phi^{(r)}}\left[ \log p(y_i,\eta_i\mid\phi) \right]. $$ The conditional distribution inside the expectation is computed at $\phi^{(r)}$ and held fixed while maximizing over $\phi$.

A common NLME specification has $\eta_i\sim N(0,\Omega)$, with observations conditionally distributed around a structural model $f(t_{ij},\eta_i;\beta)$. Here $\beta$ denotes fixed effects, $\Omega$ is the random-effect covariance matrix, and residual parameters determine the conditional observation distribution. The usual likelihood factorization assumes independent subjects and, conditional on $\eta_i$ and the population parameters, the specified within-subject observation model. The M-step formulas below also require the relevant conditional expectations to exist and the proposed parameter values to satisfy the model’s constraints.

The word “expected” is essential. For example, the random-effect contribution to $Q$ involves $E(\eta_i\eta_i^\mathsf{T}\mid y_i,\phi^{(r)})$, not merely the outer product of the conditional mean $E(\eta_i\mid y_i,\phi^{(r)})$. Replacing a latent effect by its conditional mean generally discards conditional uncertainty and changes the M-step.

Gaussian Random-Effect Covariance Updates

Suppose the random effects are independent across subjects and $\eta_i\sim N(0,\Omega)$, where $\Omega$ is an unrestricted positive-definite $q\times q$ covariance matrix. The random-effect part of the complete-data log-likelihood, omitting constants, is $$ -\frac{N}{2}\log|\Omega|-\frac{1}{2}\sum_{i=1}^{N}\eta_i^\mathsf{T}\Omega^{-1}\eta_i. $$ Taking the E-step expectation replaces each quadratic form by its conditional expectation. The identity $\eta_i^\mathsf{T}\Omega^{-1}\eta_i=\operatorname{tr}(\Omega^{-1}\eta_i\eta_i^\mathsf{T})$ then gives a matrix expression that can be maximized directly.

Q_{\Omega}(\Omega)= -\frac{N}{2}\log|\Omega| -\frac{1}{2}\operatorname{tr}\left(\Omega^{-1}S\right)+C, \qquad S=\sum_{i=1}^{N}E(\eta_i\eta_i^\mathsf{T}\mid y_i,\phi^{(r)}).

For positive-definite $S$, differentiating with respect to $\Omega$ and setting the derivative to zero yields $\Omega^{-1}S\Omega^{-1}=N\Omega^{-1}$, so the unrestricted maximum is $\Omega^{(r+1)}=S/N$. Equivalently, the update averages the conditional second-moment matrices. If $\Omega$ is constrained, for example to be diagonal or to have a particular covariance structure, that unrestricted update may not satisfy the constraint; the constrained M-step must then maximize $Q_\Omega$ over the allowed matrices.

Worked Example: Updating a Random-Intercept Variance

Consider a hypothetical model with one random intercept $b_i\sim N(0,\tau^2)$ for each of $N=3$ subjects. Suppose the E-step at the current parameter values gives $E(b_i^2\mid y_i,\phi^{(r)})$ equal to $0.30,\ 0.25,$ and $0.50$. These are conditional second moments, not just squared conditional means. The variance update is $$ \tau_{\text{new}}^2= \frac{1}{3}\sum_{i=1}^{3}E(b_i^2\mid y_i,\phi^{(r)}) =\frac{0.30+0.25+0.50}{3} =\frac{1.05}{3}=0.350. $$

Thus the M-step estimate of the random-intercept variance is $0.350$ in squared random-intercept units. Its square root, if needed to report the random-intercept standard deviation, is $\sqrt{0.350}\approx0.592$ in random-intercept units. The calculation is valid for the stated zero-mean Gaussian random-intercept model with an unconstrained positive variance; it would not by itself enforce an additional variance floor or other structural constraint.

Fixed Effects: Explicit and Numerical Updates

Whether fixed effects have a closed-form update depends on how they enter the conditional observation model. A linear Gaussian mean with fixed residual variance produces a quadratic function of the fixed effects, so its maximizer is a weighted least-squares solution. In contrast, a nonlinear structural model typically makes $Q$ nonlinear in its fixed effects. The E-step expectations can still be used to define the objective, but finding its maximum then requires an appropriate numerical optimizer.

Worked Example: A Closed-Form Fixed-Effect Update

Suppose the hypothetical observation model is $y_{ij}=\beta x_{ij}+b_i+\epsilon_{ij}$, with $\epsilon_{ij}\sim N(0,\sigma^2)$. Treat $\sigma^2$ as fixed during this part of the M-step. The $Q$-terms involving $\beta$ are proportional to $$ -\frac{1}{2\sigma^2}\sum_{i,j} E\left[(y_{ij}-\beta x_{ij}-b_i)^2\mid y_i,\phi^{(r)}\right]. $$ Differentiating with respect to $\beta$, setting the derivative to zero, and solving gives $$ \beta_{\text{new}}= \frac{\sum_{i,j}x_{ij}\{y_{ij}-E(b_i\mid y_i,\phi^{(r)})\}} {\sum_{i,j}x_{ij}^2}, $$ provided the denominator is positive. The conditional variance of $b_i$ contributes to the expected squared residual, but does not change this derivative when it is independent of $\beta$.

For a small hypothetical data set, let the three pairs $(x_{ij},y_{ij})$ be $(1,2.0),(2,4.1),(1,3.8)$. Suppose the corresponding conditional means of the subject effects for these observations are $0.2,-0.1,0.4$, respectively. Then the numerator is $1(2.0-0.2)+2(4.1-(-0.1))+1(3.8-0.4)=1.8+8.4+3.4=13.6$. The denominator is $1^2+2^2+1^2=6$, giving $\beta_{\text{new}}=13.6/6\approx2.267$. This update assumes the stated linear mean, common fixed $\sigma^2$, and nonzero design information; it is not a general formula for nonlinear fixed effects.

Worked Example: Numerical Maximization for a Nonlinear Parameter

Consider a hypothetical nonlinear M-step in which a one-dimensional profile of the expected complete-data log-likelihood, after omitting constants, is $Q(x)=6.1x-1.3x^2-5x^4$, with $x=e^{-k}$ and $k>0$. Thus the admissible range is $0<x<1$. Differentiation gives $Q'(x)=6.1-2.6x-20x^3$, so an interior stationary point satisfies $$ 20x^3+2.6x-6.1=0. $$

The left side is strictly increasing for $x>0$, since its derivative is $60x^2+2.6>0$. Its value at zero is $-6.1$, and at one it is $16.5$, so there is exactly one root in the admissible interval. Numerical root finding gives $x\approx0.6090$. Since $Q'(x)$ changes from positive to negative at this root, it is the maximum over $0<x<1$. Transforming back, $k=-\log(0.6090)\approx0.4959$, or $k\approx0.496$ in inverse units of whatever time scale defines $k$. In an actual fit, this root would be found by a numerical optimizer or root solver, and its solution and constraints would need to be checked rather than assumed.

Residual-Parameter Updates

For independent Gaussian residuals with common variance $\sigma^2$, write $r_{ij}(\eta_i)=y_{ij}-f(t_{ij},\eta_i;\beta)$, and let $M$ be the total number of observations. The residual contribution to $Q$, omitting constants, is $$ -\frac{M}{2}\log(\sigma^2) -\frac{1}{2\sigma^2}\sum_{i,j}E\{r_{ij}(\eta_i)^2\mid y_i,\phi^{(r)}\}. $$ Differentiating with respect to $\sigma^2$ yields the closed-form update: the average expected squared residual. The residuals must be squared before taking their conditional expectation; squaring a conditional mean residual alone is generally too small because it omits conditional variance.

Worked Example: Updating a Common Residual Variance

Suppose a hypothetical data set has four observations and the E-step provides the expected squared residuals $0.64,\ 1.44,\ 2.25,$ and $0.49$, each in squared response units. Under independent Gaussian residuals with a common variance, the update is $$ \sigma_{\text{new}}^2= \frac{0.64+1.44+2.25+0.49}{4} =\frac{4.82}{4}=1.205. $$

The updated residual variance is $1.205$ squared response units, and the corresponding residual standard deviation is $\sqrt{1.205}\approx1.098$ response units. This is the unconstrained update under the common-variance Gaussian assumption. If the residual model instead includes proportional error, variance covariates, correlated errors, or another distribution, its contribution to $Q$ changes and this simple average need not maximize it.

More generally, a Gaussian residual covariance $R_i(\lambda)$, governed by residual parameters $\lambda$, contributes $$ -\frac{1}{2}\sum_i\left[ \log|R_i(\lambda)|+ E\{r_i^\mathsf{T}R_i(\lambda)^{-1}r_i\mid y_i,\phi^{(r)}\} \right] $$ up to constants. If the covariance has a simple scale parameter, an explicit update may be available. With nonlinear variance functions or multiple coupled residual parameters, the M-step is generally a constrained numerical optimization problem.

Reliable M-Step Computation

A correct M-step maximizes the same $Q(\phi\mid\phi^{(r)})$ defined by the E-step. The conditional expectations are not recomputed at trial values of $\phi$ during this maximization. Recomputing them changes the objective and mixes the two steps. When parameters are separable in $Q$, blocks can be updated individually; otherwise, a blockwise or joint optimizer must account for their dependence.

Specify the objective and constraints.
Write the relevant expected complete-data log-likelihood terms, identify which parameters enter each term, and state constraints such as positive variances or positive-definite covariance matrices.

Use explicit updates where justified.
Apply a closed-form formula only when its model assumptions and parameterization match the objective being maximized.

Optimize the remaining parameters.
For nonlinear or coupled blocks, use a constrained numerical method and evaluate the objective using the fixed E-step expectations.

Check the result.
Confirm feasibility, inspect optimizer convergence and gradients, and compare the new and old $Q$-values. A small gradient at a feasible interior solution supports stationarity, but does not alone establish a global maximum.

Useful failure signals include a negative variance, a non-positive-definite covariance matrix, an optimizer stopping at a constraint without a plausible boundary solution, or a reported solution whose $Q$-value is lower than a feasible starting value. Such outcomes can arise from an incorrect derivative, unstable numerical arithmetic, a poorly scaled parameterization, or inaccurate E-step expectations. In nonlinear models, multiple local optima are also possible, so starting values and objective profiles can matter.

Update type	Typical situation	Main check
Closed-form	Unrestricted Gaussian random-effect covariance or common Gaussian residual variance	Verify moments, denominators, and covariance constraints
Closed-form linear solve	Fixed effects in a linear Gaussian mean with fixed covariance parameters	Check design rank and correct weighting
Numerical optimization	Nonlinear fixed effects or structured residual covariance parameters	Check constraints, convergence, and the value of the intended $Q$

The main alternative is to optimize the observed-data likelihood directly, evaluating or approximating the random-effect integrals at each trial parameter value. That approach does not use an EM M-step, and can be attractive when the observed likelihood and its derivatives are reliably available. The M-step instead maximizes an expected complete-data objective and may simplify some parameter blocks, but nonlinear expectations can make its numerical optimization difficult. Neither a closed-form block update nor a numerical optimizer removes the need to verify that the assumed model and the computed expectations are appropriate.

Using only conditional means. Substituting $E(\eta_i\mid y_i)$ for $\eta_i$ in a quadratic term ignores conditional uncertainty. Use the required conditional second moment, such as $E(\eta_i\eta_i^\mathsf{T}\mid y_i)$.
Applying a familiar formula to a different error model. The residual-variance average applies to independent Gaussian errors with a common variance, not automatically to proportional or correlated errors.
Assuming an optimizer’s status message proves maximization. Check feasibility and the objective; a convergence flag alone does not prove a global maximum.

Key takeaway: The M-step maximizes $Q(\phi\mid\phi^{(r)})$ with the E-step distribution fixed. Gaussian second moments often give explicit variance updates; simple linear mean structures can give explicit fixed-effect updates. Nonlinear or constrained parameter blocks generally require careful numerical optimization.

Check Your Understanding

Answer each question, then compare with the brief solutions below.

For zero-mean Gaussian random effects, which conditional moment determines the covariance update: the conditional mean or the conditional second-moment matrix?
A common Gaussian residual variance has expected squared residuals $1.0,\ 2.0,\ 3.0$ for three observations. What is the M-step variance update?
Why is a nonlinear fixed-effect update not generally obtained by substituting conditional random-effect means into the structural model and solving as if the latent effects were observed?
During a numerical M-step, should the E-step conditional distribution be recomputed at each trial value of the population parameters?
Name two checks that are important after a constrained numerical M-step.

Solutions. (1) The update uses the conditional second-moment matrix $E(\eta_i\eta_i^\mathsf{T}\mid y_i,\phi^{(r)})$, which includes conditional uncertainty. (2) The variance is $(1.0+2.0+3.0)/3=2.00$ squared response units. (3) The expected complete-data objective generally depends on conditional expectations of the complete-data log-likelihood, not only on the likelihood evaluated at a conditional mean; nonlinearities make those operations unequal. (4) No. Hold the E-step distribution at $\phi^{(r)}$ fixed while maximizing $Q(\phi\mid\phi^{(r)})$. (5) For example, verify parameter feasibility and positive definiteness, inspect convergence or gradients, and check that the computed objective is not below its feasible starting value.