Direct Maximum Likelihood Optimization

Why Optimize the Integrated Likelihood Directly?

In a shared-random-effects joint model, the longitudinal and event-time outcomes for a subject are conditionally modeled given that subject’s latent random effects. Those effects are not observed, so maximum likelihood for the observed data requires integrating them out. Direct maximum likelihood optimization treats this integrated likelihood as the objective and searches over the model parameters themselves, rather than iterating through an EM expectation and maximization step.

This approach is attractive when the marginal likelihood can be evaluated accurately at proposed parameter values and a suitable numerical optimizer is available. It also has practical challenges: covariance parameters must remain valid throughout the search, and an optimizer can converge to a maximum of an inaccurately evaluated likelihood. The integration method and its numerical tolerances therefore matter as much as the choice of optimizer.

Definition: Let \(Y_i\) denote subject \(i\)’s longitudinal data, \(T_i\) the observed event or censoring time, and \(\delta_i\) the event indicator. Let \(b_i\) be the subject’s random-effects vector and \(\theta\) the model parameter vector, including fixed effects, random-effects covariance parameters, survival parameters, and any longitudinal-survival association parameter. The subject’s marginal likelihood contribution is \(L_i(\theta)=\int f_i(Y_i,T_i,\delta_i\mid b_i;\theta)\,p(b_i;\theta)\,db_i\). If subjects are independent, the observed-data likelihood is \(L(\theta)=\prod_{i=1}^n L_i(\theta)\), and the marginal log likelihood is \(\ell(\theta)=\sum_{i=1}^n \log L_i(\theta)\).

For example, under right censoring, \(\delta_i=1\) indicates an observed event and \(\delta_i=0\) indicates right censoring. The conditional event-time contribution has the form \(h_i(T_i\mid b_i;\theta)^{\delta_i}S_i(T_i\mid b_i;\theta)\), where \(h_i\) and \(S_i\) are the conditional hazard and survival functions. It is combined with the conditional longitudinal density before integration. This likelihood formulation assumes the specified conditional models and random-effects distribution are appropriate, subjects are independent, and censoring is noninformative under the model (or is otherwise modeled). The random-effects integral must also be finite.

Enforcing Valid Covariance Parameters

An unconstrained numerical optimizer may propose negative variances or an invalid correlation. A covariance matrix for random effects must be positive definite (or, in some models, positive semidefinite). Rather than repeatedly rejecting invalid proposals, reparameterize so that ordinary unconstrained optimizer coordinates map to valid parameters.

For a random-effect standard deviation \(\sigma\), use \(\sigma=\exp(\lambda)\), so every finite \(\lambda\) gives \(\sigma>0\). Equivalently, a variance can be written as \(\exp(\gamma)\). For a correlation \(\rho\) in a two-dimensional random-effects covariance matrix, set \(\rho=\tanh(\zeta)\); any finite \(\zeta\) gives \(-1<\rho<1\). If the model permits a zero variance or a correlation exactly at a boundary, these values are limiting cases, not attained by finite transformed coordinates.

D= \begin{pmatrix} \sigma_1^2 & \rho\sigma_1\sigma_2\\ \rho\sigma_1\sigma_2 & \sigma_2^2 \end{pmatrix}, \qquad \sigma_j=\exp(\lambda_j),\quad \rho=\tanh(\zeta).

For a higher-dimensional random-effects vector, a standard alternative is the Cholesky parameterization \(D=CC^\mathsf{T}\), with \(C\) lower triangular and its diagonal entries constrained to be positive by exponentiating unconstrained coordinates. This produces a positive-definite covariance matrix for finite transformed coordinates. Fixed effects and many association parameters can remain on their original, unconstrained scales; positive hazard or scale parameters can likewise be represented with logarithms.

Worked Example: Transforming a Random Intercept-Slope Covariance

Suppose a hypothetical model uses a random intercept and random slope with standard deviations \(0.8\) outcome units and \(1.2\) outcome units per year, and correlation \(-0.35\). Calculate unconstrained coordinates and verify the resulting covariance matrix is valid.

The transformed coordinates are \(\lambda_1=\log(0.8)\approx-0.2231\), \(\lambda_2=\log(1.2)\approx0.1823\), and \(\zeta=\operatorname{atanh}(-0.35)\approx-0.3654\). Transforming back gives \(\exp(\lambda_1)=0.8\), \(\exp(\lambda_2)=1.2\), and \(\tanh(\zeta)=-0.35\).

The covariance is \(D=\begin{pmatrix}0.8^2 & (-0.35)(0.8)(1.2)\\(-0.35)(0.8)(1.2)&1.2^2\end{pmatrix} =\begin{pmatrix}0.64&-0.336\\-0.336&1.44\end{pmatrix}\). Its determinant is \(0.64(1.44)-(-0.336)^2=0.9216-0.112896=0.808704>0\). Its first diagonal entry is positive, so this symmetric two-by-two matrix is positive definite. The off-diagonal entry has units of outcome units squared per year.

Key result: If \(\eta\) is an unconstrained coordinate and \(\theta=g(\eta)\) is a one-to-one transformation onto the interior of the parameter space, maximizing \(\ell(g(\eta))\) over \(\eta\) gives the same maximum-likelihood estimate on the original scale as maximizing \(\ell(\theta)\) over that interior. No Jacobian factor is added to a frequentist likelihood. A Jacobian is used when transforming a probability density, as in a change of variables for a Bayesian prior or posterior density, not simply to re-express the parameter over which a likelihood is maximized.

The reason is direct: the transformation changes the coordinates used to search, not the likelihood assigned to a particular model parameter. Every admissible interior value of \(\theta\) corresponds to one \(\eta\), and the objective at that point is the same \(\ell(\theta)\). Thus the ordering of candidate parameter values by likelihood is unchanged. The transformed objective may have different numerical curvature, however, so the optimizer’s path and convergence behavior can differ.

Numerical Integration Inside the Objective

In general, the integral defining \(L_i(\theta)\) has no closed form. A numerical method approximates it at each proposed \(\theta\). For instance, quadrature replaces an integral by a weighted sum of integrand evaluations. Importance sampling instead estimates the integral with weighted evaluations drawn from a proposal distribution. The earlier material on importance sampling for random-effects integration covers that integration strategy; for direct optimization, the essential point is that the same marginal likelihood objective must be evaluated consistently at every parameter proposal.

Write \(\widetilde L_i(\theta)\) for a numerically computed approximation to \(L_i(\theta)\). The optimizer then works with \(\widetilde\ell(\theta)=\sum_i\log\widetilde L_i(\theta)\). A usable approximation should keep subject contributions positive and finite and be accurate enough that numerical error does not change which candidate appears better. In particular, a small apparent improvement in \(\widetilde\ell\) is not persuasive if it is of the same size as changes caused by tightening integration tolerances or increasing integration resolution.

Worked Example: Comparing Two Candidate Fits

Suppose a hypothetical three-subject joint model is evaluated by numerical integration at two candidate parameter vectors, A and B. The computed subject likelihood contributions are A: \(0.018,\ 0.011,\ 0.024\), and B: \(0.020,\ 0.010,\ 0.025\). Assume both candidates use the same likelihood definition and integration settings.

The corresponding log likelihoods, calculated from the unrounded contributions, are \(\ell_A=\log(0.018)+\log(0.011)+\log(0.024)\approx-12.2569\) and \(\ell_B=\log(0.020)+\log(0.010)+\log(0.025)\approx-12.2061\). The difference should be computed from the unrounded contributions: \(\ell_B-\ell_A=\log(0.020/0.018)+\log(0.010/0.011)+\log(0.025/0.024)\approx0.0508723\), or about \(0.0509\).

Thus B has the larger approximate log likelihood, by about \(0.0509\), and is preferred between these two candidates by this criterion. This comparison does not establish that B is the global maximum or that the difference is meaningful relative to integration error. Subtracting the displayed rounded log likelihoods gives \(0.0508\); the slight discrepancy is rounding, which is why comparisons should use unrounded values.

Running and Checking the Optimization

A direct fit typically supplies the optimizer with transformed coordinates, maps them to the original model parameters, evaluates the numerically integrated log likelihood, and returns that objective. The optimizer searches for a maximum, often by minimizing the negative log likelihood. It is important to check the final result on both scales: transformed coordinates reveal whether the search has run far toward a boundary, while original-scale estimates are the quantities needed to interpret the fitted model.

Choose valid coordinates.
Use logarithms for positive parameters and a Cholesky or correlation transformation for covariance parameters. Specify which parameters are unrestricted and which have boundaries.

Evaluate the marginal objective.
At each proposal, compute all subject integrals using the stated numerical method, then sum their log contributions. Check that the contributions are positive and finite.

Optimize and restart.
Run the optimizer from multiple plausible starting points. A small change in its objective or transformed coordinates is useful evidence of numerical convergence, but it is not proof of a global maximum.

Check numerical stability.
Re-evaluate the fitted parameters with tighter integration tolerances or a more accurate integration setting. Confirm that the objective and important estimates remain stable at the precision needed for the analysis.

Worked Example: Distinguishing Optimization Convergence From Stability

Consider a hypothetical model fit in transformed coordinates. The negative log likelihood at the starting point is \(48.31\). The optimizer reports successive values \(45.82,\ 44.76,\ 44.73\), and terminates with a small change in its final iterations. A second starting point reaches a negative log likelihood of \(44.729\).

Since maximization of the log likelihood is equivalent to minimization of its negative, these values improve from \(48.31\) to \(44.73\). The two starts differ at the reported precision by \(44.730-44.729=0.001\), which is encouraging evidence that both reached a similar solution. It is not enough on its own to show that numerical integration is accurate.

Suppose a more accurate integration setting gives a final negative log likelihood of \(44.728\), while the random-effects standard deviations change from \(0.81\) and \(1.17\) to \(0.81\) and \(1.18\), respectively. The objective changes by \(0.002\), and the displayed standard deviations are stable to about \(0.01\). That supports stability at these reporting precisions; it does not prove exact integration or global optimality. The analyst should also examine whether any transformed covariance coordinate is extreme, since a very large absolute correlation coordinate can indicate a fitted correlation close to its boundary.

How Direct Optimization Compares With EM

EM and direct optimization target the same observed-data likelihood when they are correctly implemented for the same model. EM uses latent random effects to construct an auxiliary function and alternates an E-step with an M-step. Direct optimization instead evaluates the marginal likelihood and searches over its parameters. EM can be convenient when conditional expectations and M-step updates are manageable, while direct optimization can be attractive when a reliable integrated objective is available and the full parameter vector is practical to optimize.

EM’s likelihood improvement per iteration can be modest near a maximum, and the algorithm may require many iterations. Direct optimization may converge in fewer objective evaluations but can be sensitive to starting values, parameter scaling, and numerical integration noise. Monte Carlo EM adds simulation error to the E-step, whereas direct optimization with deterministic quadrature has quadrature error instead. These approaches are not interchangeable guarantees: all require careful convergence and likelihood checks.

Approach	What is optimized or updated?	Important numerical concern
Direct marginal-likelihood optimization	The numerically integrated observed-data log likelihood	Integration error can distort the objective and its apparent maximum
EM	An auxiliary function based on the latent random effects	Slow progress near convergence; E- and M-steps must be implemented correctly
Monte Carlo EM	A simulation-based approximation to the E-step followed by an M-step	Monte Carlo variation can obscure small changes in the objective

Failure Modes and Limitations

A converged optimizer has found a point satisfying its stopping criteria for the objective it was given. That statement alone does not show that the point is a global maximum, that the integrated likelihood was evaluated accurately, or that the statistical model is adequate. Common warning signs include materially different solutions from different starts, unstable estimates as integration accuracy changes, invalid or nonfinite subject contributions, and covariance parameters driven close to a boundary.

Boundary estimates. A transformed standard deviation can approach zero only as its log coordinate tends to negative infinity. A correlation can approach \(-1\) or \(1\) only as its transformed coordinate diverges. Extreme coordinates can signal a genuine boundary estimate, weak identification, or a numerical problem; investigate rather than interpreting them automatically as ordinary interior estimates.
Integration noise mistaken for improvement. If changing integration accuracy changes the objective as much as the optimizer’s reported improvement, the maximum is not yet numerically established.
Poor parameter scaling. Parameters on very different scales can make an optimizer’s steps inefficient. Transformation and sensible scaling can help, but do not fix a misspecified likelihood.
Over-reliance on one start. A local search can settle at a local maximum or a flat region. Multiple starts and comparison of final likelihoods help reveal this risk.

Direct optimization also has a computational cost: many likelihood evaluations may each require a multidimensional integral. Its limitations include dependence on the accuracy and stability of that integration, possible difficulties near parameter boundaries, and the inability of numerical convergence checks to verify the scientific assumptions of the joint model. The optimizer finds a maximum for the specified likelihood; it does not establish that the likelihood’s assumptions are true.

Common pitfall: “The optimizer converged, so the reported estimate is the maximum-likelihood estimate.” This confuses an algorithm’s stopping rule with statistical and numerical validation. Check multiple starting values, inspect boundary behavior, and verify that the integrated objective and estimates are stable to increased integration accuracy before treating the solution as a reliable maximum.

Check Your Understanding

Answer each question, then compare with the brief solutions.

A random-effect standard deviation is optimized using \(\sigma=\exp(\lambda)\). What values can \(\sigma\) take for finite \(\lambda\), and how is zero approached?
For a two-dimensional covariance matrix, \(\sigma_1=0.5\), \(\sigma_2=2\), and \(\rho=0.4\). What is the off-diagonal covariance?
Why is a Jacobian not added when maximizing the likelihood in transformed coordinates?
Two starting values give nearly identical final objective values, but tighter integration settings change the objective substantially. What conclusion is justified?
Give one advantage and one numerical concern for direct optimization compared with EM.

Solutions. (1) Every finite \(\lambda\) gives \(\sigma>0\); \(\sigma\) approaches zero as \(\lambda\) tends to negative infinity. (2) The covariance is \(\rho\sigma_1\sigma_2=(0.4)(0.5)(2)=0.4\), in the product units of the two random effects. (3) The transformation changes the coordinates of the parameter search, not the likelihood value assigned to each parameter; a Jacobian is for transforming a probability density. (4) The optimizer appears reproducible for the approximate objective, but integration stability is not established; improve integration accuracy and reassess the fit. (5) Direct optimization works on the integrated observed-data objective and may avoid many EM iterations; its maximum can be distorted by integration error or sensitivity to starting values.