Why Integration Is Challenging for Discrete Markers
A non-Gaussian longitudinal model specifies the probability of each observed marker value conditional on a subject’s random effects. The joint model’s likelihood for that subject is not just this conditional probability: the random effects are unobserved, so their values must be integrated out. With discrete markers, the conditional likelihood can change rapidly across the random-effects distribution. Repeated binary observations, for example, can give much greater weight to a narrow range of random-intercept values than the normal random-effects distribution does on its own.
This tutorial focuses on whether numerical integration is adequate for that marginal likelihood. A fit can appear to converge while still depending materially on the number or placement of quadrature points. In particular, a stable optimizer is not evidence by itself that the integral has been approximated accurately.
This expression relies on a correctly specified conditional marker model, a correctly specified random-effects distribution, and the stated conditional-independence structure. The event-time contribution also relies on the model’s censoring assumptions; in the usual setup, censoring is noninformative conditional on the modeled covariates and random effects. If the random-effects distribution or dependence assumptions are wrong, more accurate integration will not repair the model.
From a Normal Expectation to Gauss–Hermite Quadrature
For numerical work it is useful to express the integral as an expectation. If \(b_i\sim N_q(0,D)\), then $$ L_i(\theta)=E_{b_i}\left[f(y_i\mid b_i;\theta)f_T(T_i,\delta_i\mid b_i;\theta)\right]. $$ Let \(D=LL^\mathsf{T}\) be a Cholesky factorization and let \(z\sim N_q(0,I)\), so that \(b_i=Lz\). Gauss–Hermite quadrature approximates an integral against the weight \(\exp(-x^2)\). Changing variables from a standard normal coordinate \(z_k\) to \(x_k=z_k/\sqrt{2}\) gives the approximation $$ L_i(\theta)\approx \pi^{-q/2}\sum_{r_1=1}^{m}\cdots\sum_{r_q=1}^{m} \left(\prod_{k=1}^{q} w_{r_k}\right) f(y_i\mid b_{r};\theta)f_T(T_i,\delta_i\mid b_{r};\theta), \qquad b_r=\sqrt{2}Lx_r. $$ Here \(x_r\) and \(w_r\) are the one-dimensional Gauss–Hermite nodes and weights, and \(m\) is the number of nodes per dimension. This is a weighted sum of conditional likelihoods at selected random-effects values. Increasing \(m\) usually improves the approximation when the integrand is smooth and well covered by the nodes, but it does not guarantee adequacy in every model.
Worked Example: A Binary Marker Creates a Concentrated Integrand
Worked Example: A Binary Marker Creates a Concentrated Integrand
Suppose, purely hypothetically, that one subject has 12 binary marker observations, all equal to 1. Use a logistic random-intercept model with \(p(b)=\operatorname{expit}(b)=1/(1+\exp(-b))\), and assume \(b\sim N(0,1)\). For these marker observations, the integrand in the marginal probability is $$ g(b)=p(b)^{12}\phi(b;0,1). $$ The assumptions for this calculation are a conditional Bernoulli model with independent observations given \(b\), the stated logistic link, and a standard normal random intercept.
At \(b=0\), \(p(0)=0.5\) and \(\phi(0;0,1)=0.39894\). Therefore $$ g(0)=(0.5)^{12}(0.39894)=0.0000974. $$ At \(b=1.76\), \(p(1.76)\approx0.85321\), \(p(1.76)^{12}\approx0.1488\), and \(\phi(1.76;0,1)\approx0.08478\). Thus $$ g(1.76)\approx(0.85321)^{12}(0.08478)\approx(0.1488)(0.08478)\approx0.01262. $$ The ratio is \(g(1.76)/g(0)\approx0.01262/0.0000974\approx129.6\). Although \(b=1.76\) is away from the random-effects mean, the marker likelihood makes the integrand there about 129.6 times its value at zero.
The log integrand, up to a constant, is \(12\log p(b)-b^2/2\). Its first derivative is \(12\{1-p(b)\}-b\), and its second derivative is \(-12p(b)\{1-p(b)\}-1\). At \(b=1.76\), the first derivative is approximately \(12(1-0.85321)-1.76\approx0.0015\), close to zero; the mode is about 1.76. The second derivative there is approximately \(-12(0.85321)(1-0.85321)-1\approx-2.50\), indicating local concentration around that mode. This calculation explains why nodes concentrated around zero alone can miss important integrand mass. It does not, by itself, tell us how many nodes suffice for a particular fitted model.
Worked Example: Checking Likelihood and Estimate Stability
Worked Example: Checking Likelihood and Estimate Stability
Consider a hypothetical joint-model fit to repeated binary markers and event times. The model has one normal random intercept and an association parameter \(\alpha\) linking that random effect to the event submodel. The following illustrative software output compares the same model and data fitted with increasing one-dimensional Gauss–Hermite node counts. The reported log likelihood is the maximized log likelihood under that quadrature rule.
| Nodes per dimension | Maximized log likelihood | Estimate of \(\alpha\) | Estimate of random-intercept SD |
|---|---|---|---|
| 15 | -412.638 | 0.714 | 0.936 |
| 30 | -412.521 | 0.728 | 0.961 |
| 60 | -412.520 | 0.728 | 0.962 |
From 15 to 30 nodes, the maximized log likelihood changes by \(-412.521-(-412.638)=0.117\), and \(\alpha\) changes by \(0.728-0.714=0.014\). Those shifts suggest that 15 nodes may not be adequate for stable estimation. From 30 to 60 nodes, the log likelihood changes by \(0.001\), while the displayed estimates of \(\alpha\) and the random-intercept standard deviation change by 0.000 and 0.001, respectively. On these summaries, 30 nodes appear much more adequate than 15, and the comparison provides little evidence of further change at 60 nodes.
The rounded table is not enough to establish that every inferential quantity is stable. A careful check also compares unrounded parameter estimates, standard errors, and relevant derived quantities, and verifies that the optimizer converged at each node count. A small change in log likelihood alone cannot rule out a meaningful shift in a weakly identified parameter. These values are an illustration of the diagnostic calculation, not results from an actual study.
Worked Example: The Cost of More Dimensions
Worked Example: The Cost of More Dimensions
Suppose a hypothetical model has three random effects and uses a tensor-product Gauss–Hermite rule with \(m=15\) nodes per dimension. The number of node combinations evaluated for one subject is \(15^3=3{,}375\). Increasing to 30 nodes per dimension requires \(30^3=27{,}000\) combinations, eight times as many because \((30/15)^3=8\). Increasing to 60 requires \(60^3=216{,}000\) combinations, or 64 times the original number because \((60/15)^3=64\).
This calculation checks the computational cost under the tensor-product assumption: every node in each dimension is combined with every node in the other dimensions. It shows why a high node count that is feasible for one random effect may be impractical for several. In a joint model, each combination may require evaluating multiple marker and event contributions, and this work is repeated over subjects and optimizer iterations. A researcher should therefore report the rule used and the adequacy checks, not just the final parameter estimates.
Diagnostics, Alternatives, and Limitations
A useful diagnostic is a sequence of fits with increasing node counts, using the same data, model, parameterization, and optimization tolerances. Compare maximized log likelihoods and estimates, including their standard errors. Examine the parameters that matter for the scientific interpretation, not just the fixed effects. Also check whether estimates are on boundaries, whether the optimizer reports convergence, and whether starting values or optimization settings materially change the result. When feasible, repeat the comparison with a different integration strategy to identify sensitivity to the node placement.
Fixed Gauss–Hermite quadrature places nodes according to the normal random-effects distribution. If the subject-specific likelihood shifts the important mass away from that distribution’s center, or is sharply concentrated, a modest fixed rule can approximate the integral poorly. Adaptive Gauss–Hermite quadrature recenters and rescales nodes using the local shape of the integrand, which can be more efficient for a concentrated, approximately unimodal integrand. Adaptation is not a universal solution: an inaccurate mode or curvature estimate, a multimodal integrand, or difficult tails can still cause problems.
Monte Carlo integration instead approximates the expectation by averaging conditional likelihood values at random draws from the random-effects distribution. It avoids a deterministic tensor grid and can be useful as dimension increases, but introduces simulation error. Repeated runs or increased simulation size are needed to assess whether that error affects the likelihood and estimates. These methods are alternatives for approximating the same marginal likelihood; none removes the need to check numerical stability.
A common mistake is to conclude that integration is adequate because the optimizer converged or because the estimates look plausible. That reasoning confuses optimization of an approximate likelihood with accuracy of the approximation itself. The correction is to increase the quadrature resolution, compare likelihood and parameter results, and examine relevant standard errors and diagnostics. Another mistake is to keep increasing the node count without considering the model dimension: a tensor grid’s cost grows as \(m^q\), so a more suitable integration strategy may be needed.
Finally, quadrature stability is a numerical check, not evidence that the model assumptions are true. It cannot validate the chosen marker distribution, random-effects distribution, event model, or conditional-independence assumptions. Nor does a stable rounded table prove that all quantities are stable. Report enough detail for readers to understand the integration rule and the checks used.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- Why can 12 repeated binary observations all equal to 1 shift the important integrand mass away from \(b=0\), even when the random-effects distribution is centered at zero?
- A one-dimensional fit has maximized log likelihoods of \(-250.10\), \(-250.02\), and \(-250.02\) at 10, 20, and 40 nodes. The estimate of a key association parameter changes from 0.40 to 0.46 to 0.46. What does this suggest, and what else should be checked?
- For a tensor-product rule with four random effects and 12 nodes per dimension, how many node combinations are evaluated for one subject?
- How does adaptive Gauss–Hermite quadrature differ from fixed Gauss–Hermite quadrature, and what problem can remain even with adaptation?
Solutions. (1) The conditional likelihood \(p(b)^{12}\) is much larger at positive \(b\) than at zero; its product with the normal density can therefore peak away from the normal mean, as the worked calculation shows. (2) The change from 10 to 20 nodes is \(0.08\) log-likelihood units and 0.06 for the association estimate, while the 20- and 40-node results agree at the displayed precision. This suggests 10 nodes may be inadequate and 20 may be sufficient for these reported quantities. Check unrounded estimates, standard errors, optimizer convergence, and other parameters or derived quantities before deciding. (3) There are \(12^4=20{,}736\) combinations. (4) Adaptive quadrature recenters and rescales nodes using the integrand’s local mode and curvature; it can still perform poorly if those local features are estimated badly or the integrand is multimodal or has difficult tails.