Numerical Estimation Of Joint Endpoint Models

Why Numerical Integration Is Needed

In a joint model, shared random effects can represent subject-level features that influence more than one endpoint. Conditional on those effects, the endpoint likelihoods may be straightforward to calculate. But the likelihood for the observed data must account for the fact that the effects are unobserved. For nonlinear endpoint models, integrating them out is usually not available in closed form. Numerical integration is therefore part of fitting the model, not merely a final computational detail.

This tutorial derives the marginal likelihood for two endpoints and compares three ways to evaluate its integral: adaptive Gaussian quadrature, the Laplace approximation, and Monte Carlo integration. The discussion assumes the model structure and cross-endpoint association have already been specified. It focuses on evaluating the resulting likelihood accurately enough for parameter estimation.

Conditional and Marginal Likelihoods

Consider subjects $i=1,\ldots,n$, with two endpoints observed at one or more occasions $j$. Let $C_{ij}$ denote the concentration used by the PK/PD model at occasion $j$, treated here as known. For illustration, let endpoint 1 be binary and endpoint 2 continuous. A scalar shared random effect $b_i$ can enter both endpoint models:

\begin{aligned} \operatorname{logit}\{P(Y_{i1j}=1\mid b_i)\} &=\alpha_1+ \frac{E_{\max,1}C_{ij}}{EC_{50,1}+C_{ij}}+\lambda_1 b_i,\\ Y_{i2j}\mid b_i &\sim N\left( \alpha_2+\frac{E_{\max,2}C_{ij}}{EC_{50,2}+C_{ij}}+\lambda_2 b_i,\, \sigma_2^2 \right),\\ b_i&\sim N(0,\omega^2). \end{aligned}

The two endpoints have different observation models and different exposure-response parameters, but the same subject-level $b_i$ contributes to each. The example assumes that, conditional on $b_i$ and the model parameters, all endpoint observations are independent across endpoints and occasions, and that the continuous residuals are independent with common variance $\sigma_2^2$. It also assumes subjects are independent and random effects are drawn from the stated normal distribution. These assumptions define the likelihood being evaluated; they should be revisited if the data-generating structure requires additional dependence or residual variation.

Definition: Let $\mathbf y_i$ denote all observed endpoint data for subject $i$, $\boldsymbol\theta$ the fixed-effect and residual parameters, and $\mathbf b_i$ the subject's shared random-effect vector. The conditional likelihood is $L_i(\boldsymbol\theta,\mathbf b_i)=p(\mathbf y_i\mid\mathbf b_i,\boldsymbol\theta)$. The marginal likelihood contribution is the conditional likelihood averaged over the random-effect distribution: $$ L_i^{\mathrm{marg}}(\boldsymbol\theta) =\int L_i(\boldsymbol\theta,\mathbf b_i) \,p(\mathbf b_i\mid\boldsymbol\theta)\,d\mathbf b_i. $$

Under the conditional independence assumptions above, the subject-level conditional likelihood factors into contributions from the two endpoints:

L_i(\boldsymbol\theta,b_i) = \left\{\prod_j p(Y_{i1j}\mid b_i,\boldsymbol\theta)\right\} \left\{\prod_j p(Y_{i2j}\mid b_i,\boldsymbol\theta)\right\}.

The marginal likelihood contribution is the integral of this product over $b_i$, not the product of two separately integrated endpoint likelihoods. In general, the shared effect makes those two operations different. For independent subjects, the full log likelihood is $\ell(\boldsymbol\theta)=\sum_i\log L_i^{\mathrm{marg}}(\boldsymbol\theta)$. Estimation maximizes this function over the unknown parameters.

Deriving the Integral and Its Quadrature Form

The integral follows directly from the law of total probability. For a small region around an observed outcome, its probability is obtained by averaging the conditional probability over possible random effects. In density notation, the same argument gives:

\begin{aligned} p(\mathbf y_i\mid\boldsymbol\theta) &=\int p(\mathbf y_i,b_i\mid\boldsymbol\theta)\,db_i\\ &=\int p(\mathbf y_i\mid b_i,\boldsymbol\theta) \,p(b_i\mid\boldsymbol\theta)\,db_i. \end{aligned}

The second line uses the factorization of a joint density into a conditional density and a marginal density. If $b_i\sim N(0,\omega^2)$, substitute $b_i=\omega z$, where $z\sim N(0,1)$. Then the contribution can be written as an expectation:

L_i^{\mathrm{marg}}(\boldsymbol\theta) =E_{Z}\left[L_i(\boldsymbol\theta,\omega Z)\right], \qquad Z\sim N(0,1).

Gaussian quadrature approximates this expectation by evaluating the conditional likelihood at selected nodes and combining those values with weights. If $z_q$ and $w_q$ are standard-normal quadrature nodes and weights, normalized so that the weights sum to one, then:

L_i^{\mathrm{marg}}(\boldsymbol\theta) \approx \sum_{q=1}^{Q}w_q L_i(\boldsymbol\theta,\omega z_q).

The approximation is deterministic for fixed nodes and weights. Increasing $Q$ generally improves accuracy, but the rate of improvement depends on how well the chosen nodes cover the region where the integrand contributes most. Adaptive quadrature shifts and rescales the nodes to concentrate them near the important region for each subject. In multiple dimensions, it is common to use a multivariate Gaussian proposal centered near the integrand's mode; quadrature evaluations are then adjusted for the ratio of the random-effect density to that proposal density.

Worked Example: Three-Node Gaussian Quadrature

Suppose one hypothetical subject has two binary endpoint observations, $Y_{i1}=1$ and $Y_{i2}=0$, with conditional probabilities $P(Y_{i1}=1\mid b)=\operatorname{logit}^{-1}(-0.2+b)$ and $P(Y_{i2}=1\mid b)=\operatorname{logit}^{-1}(0.5+1.5b)$. Let $b\sim N(0,0.5^2)$. The conditional likelihood for these observations is $L_i(b)=P(Y_{i1}=1\mid b)\{1-P(Y_{i2}=1\mid b)\}$. Use the three standard-normal nodes $z=(-1.2247,0,1.2247)$, with weights $(1/6,2/3,1/6)$.

Transform each node to the random-effect scale using $b=0.5z$. The resulting $b$ values are approximately $-0.6124,0,0.6124$. Substitution into the two logistic probabilities gives conditional likelihoods of approximately $0.1854,0.1700,0.1173$, respectively. For example, at $b=0$, the likelihood is $\operatorname{logit}^{-1}(-0.2)\{1-\operatorname{logit}^{-1}(0.5)\} =(0.4502)(0.3775)=0.1700$.

The three-node marginal likelihood approximation is therefore \[ \frac{1}{6}(0.1854)+\frac{2}{3}(0.1700)+\frac{1}{6}(0.1173) =0.1638 \] to four significant figures. This is an approximation to the probability of the observed pair after averaging over the shared random effect. Using two separately integrated endpoint probabilities instead would not compute this shared-effect likelihood.

Laplace Approximation

The Laplace approximation replaces the log integrand by a quadratic expansion near its mode. For a $d$-dimensional random effect, define $g_i(\mathbf b)=\log L_i(\boldsymbol\theta,\mathbf b)+\log p(\mathbf b\mid\boldsymbol\theta)$. If $\widehat{\mathbf b}_i$ is a local maximum of $g_i$, and $H_i=-\nabla^2 g_i(\widehat{\mathbf b}_i)$ is positive definite, then the approximation is:

L_i^{\mathrm{marg}}(\boldsymbol\theta) \approx \exp\{g_i(\widehat{\mathbf b}_i)\} (2\pi)^{d/2}|H_i|^{-1/2}.

To justify the formula, use a second-order Taylor expansion around the mode: $g_i(\mathbf b)\approx g_i(\widehat{\mathbf b}_i) -\tfrac12(\mathbf b-\widehat{\mathbf b}_i)^\mathsf{T} H_i(\mathbf b-\widehat{\mathbf b}_i)$. Exponentiating gives a constant times a multivariate Gaussian kernel. Integrating that kernel over $\mathbf b$ contributes the factor $(2\pi)^{d/2}|H_i|^{-1/2}$. The approximation is most credible when the log integrand is smooth and nearly quadratic over the region that contributes materially to the integral.

Worked Example: Calculating a Laplace Approximation

For a hypothetical subject in a nonlinear two-endpoint fit, suppose numerical optimization of the log integrand, including the random-effect density, gives a mode $\widehat b_i=0.20$, a log-integrand value $g_i(\widehat b_i)=-2.10$, and curvature $H_i=-g_i''(\widehat b_i)=4.00$. These are illustrative optimizer outputs, not observations from a real study. For one random effect, $d=1$, so:

\[ L_i^{\mathrm{marg}}\approx \exp(-2.10)\sqrt{\frac{2\pi}{4.00}} =(0.1225)(1.2533)=0.1534. \]

The curvature is positive as required for $H_i=-g_i''(\widehat b_i)>0$ at a local maximum of $g_i$. The approximation is the area under the local Gaussian approximation to the full integrand; it is not simply the integrand's value at the mode. A flat or asymmetric integrand may have a materially different area from this local approximation.

Monte Carlo Integration

Monte Carlo integration estimates the expectation by drawing random effects from their distribution and averaging the corresponding conditional likelihoods. With independent draws $\mathbf b_i^{(m)}$, $m=1,\ldots,M$, from $p(\mathbf b_i\mid\boldsymbol\theta)$:

\widehat L_i^{\mathrm{marg}} =\frac{1}{M}\sum_{m=1}^{M} L_i(\boldsymbol\theta,\mathbf b_i^{(m)}).

This estimator is unbiased for the marginal likelihood when the draws come from the stated random-effect distribution. Its accuracy is assessed by its Monte Carlo standard error, which decreases at an approximate rate proportional to $1/\sqrt{M}$. If draws instead come from a proposal distribution, each likelihood value must be multiplied by the appropriate density ratio; otherwise the estimator targets the wrong integral.

Worked Example: Monte Carlo Estimate and Simulation Error

Suppose four independent random-effect draws from the correct distribution give conditional likelihood values $0.12, 0.20, 0.16,$ and $0.24$ for a hypothetical subject. The Monte Carlo estimate is \[ \widehat L_i^{\mathrm{marg}}=\frac{0.12+0.20+0.16+0.24}{4} =\frac{0.72}{4}=0.1800. \]

The sample standard deviation of these four likelihood values is \[ s=\sqrt{\frac{(-0.06)^2+(0.02)^2+(-0.02)^2+(0.06)^2}{4-1}} =\sqrt{\frac{0.0080}{3}}=0.05164. \] Thus the estimated Monte Carlo standard error of the mean is $SE=s/\sqrt{4}=0.02582$. This is large relative to the estimate of $0.1800$, so four draws provide a rough calculation, not a precise integral. Increasing the number of draws can reduce simulation error, but the adequacy should be assessed from the resulting error and stability of the fit.

Choosing and Checking a Numerical Method

The methods differ in how they spend computation and what can make them inaccurate. Adaptive quadrature is often effective for low-dimensional random effects when the integrand is concentrated in a region that can be located and scaled. The Laplace approximation can be computationally efficient, especially when many subject-level integrals must be evaluated during optimization, but it relies on a good local quadratic approximation. Monte Carlo is flexible and extends naturally to higher dimensions, but introduces simulation variability and may require many draws for stable likelihood comparisons.

Method	Strength	Main numerical concern	Useful check
Adaptive quadrature	Accurate deterministic integration for low-dimensional effects	Insufficient nodes or poorly located/scaled quadrature region	Increase quadrature order and compare estimates
Laplace approximation	Fast approximation from the mode and curvature	Skewed, multimodal, or non-quadratic log integrand	Compare with higher-order quadrature or Monte Carlo where feasible
Monte Carlo integration	Flexible in higher dimensions and not tied to a local quadratic form	High simulation error or poor coverage of influential regions	Increase draws; assess standard errors and run-to-run stability

Several checks apply regardless of method. The log likelihood and parameter estimates should be stable when numerical settings are tightened. For quadrature, compare successive node counts rather than assuming a customary order is sufficient. For Laplace, verify that the mode was found and that the negative Hessian is positive definite; a non-positive curvature matrix means the stated approximation is not justified at that point. For Monte Carlo, track simulation error relative to the changes in log likelihood that influence optimization, and check whether repeated runs give materially different estimates.

In multiple dimensions, tensor-product quadrature can require a number of evaluations that grows rapidly with both dimension and node count. Monte Carlo avoids that particular growth but can be inefficient if most prior draws fall in regions contributing little to the integral. Adaptive proposals can help, but then density-ratio weighting must be implemented correctly. All three methods can struggle when a likelihood is weakly informative about random effects or has several important modes; numerical precision does not resolve model identifiability.

Common Pitfall and Practical Implications

A common mistake is to evaluate each endpoint's marginal likelihood separately and multiply the results, reasoning that this is equivalent to integrating the joint model. With shared random effects, the correct expression is $\int L_{i1}(b)L_{i2}(b)p(b)\,db$, whereas separate marginalization gives $\{\int L_{i1}(b)p(b)\,db\}\{\int L_{i2}(b)p(b)\,db\}$. These are generally unequal: the first uses the same latent subject effect in both endpoint contributions, while the second effectively averages each endpoint independently. The correction is to form the full conditional joint likelihood first and integrate it over the shared effects.

A second practical pitfall is treating a converged optimizer as proof that the integral is accurate. An optimizer can converge precisely to the maximum of a numerically inaccurate likelihood. Numerical integration settings should therefore be checked at and around the fitted solution. The method's limitations also matter: quadrature accuracy depends on node placement and dimension, Laplace accuracy depends on local shape, and Monte Carlo accuracy depends on adequate sampling. When feasible, agreement between methods with different approximation errors is reassuring, although it does not by itself establish that the endpoint model is appropriate.

Check Your Understanding

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

Why is the marginal likelihood for shared effects an integral of the product of the endpoint likelihoods?
For standard-normal Monte Carlo draws, what quantity is averaged to estimate a subject's marginal likelihood?
What does a non-positive-definite $H_i$ imply for the Laplace approximation as stated?
Why should quadrature accuracy be checked by changing the node count?

Solutions. (1) Conditional on the same subject effect, the joint likelihood is the product of the endpoint likelihoods under conditional independence; averaging that joint likelihood over the shared effect gives the marginal likelihood. (2) Average the conditional joint likelihood $L_i(\boldsymbol\theta,\omega Z^{(m)})$ across draws $Z^{(m)}\sim N(0,1)$. (3) The required local Gaussian integral is not justified by the stated formula at that point, because the negative Hessian is not positive definite; the mode or curvature calculation needs investigation. (4) A finite quadrature rule is an approximation, and stability across increasing node counts provides evidence that the numerical integral is adequately resolved.