Computing the Marginal Likelihood Score

Why the Marginal Score Is a Posterior Average

In a joint model, each subject’s longitudinal measurements and event-time outcome are connected through latent random effects. The likelihood contribution observed for that subject integrates those effects out. Differentiating the resulting integral gives the marginal likelihood score: the rate of change of the observed-data log likelihood with respect to a model parameter. It is central to direct likelihood optimization, and its calculation must account for uncertainty about each subject’s random effects rather than treating them as known.

The key result is that the marginal score is an average of the complete-data score, with the average taken under the random-effects distribution conditional on that subject’s observed data. Numerical integration can approximate this posterior average using the quadrature methods introduced in the earlier material on numerical integration.

Definition: Let $D_i$ denote all observed data for subject $i$, including longitudinal measurements and the event-time or censoring information. Let $b_i$ be that subject’s random-effects vector, and let $\theta$ collect model parameters. Write $f_i(D_i,b_i;\theta)$ for the complete-data likelihood contribution, including the random-effects density. The marginal contribution is $$ L_i(\theta)=\int f_i(D_i,b_i;\theta)\,db_i. $$ The complete-data score is $s_i^c(b_i;\theta)=\partial\log f_i(D_i,b_i;\theta)/\partial\theta$. The posterior density of the random effects is $p(b_i\mid D_i;\theta)=f_i(D_i,b_i;\theta)/L_i(\theta)$.

For example, under conditional independence of longitudinal and event-time outcomes given $b_i$, a complete-data contribution can be written as the product of the longitudinal density, the event-time likelihood, and the random-effects density. With right censoring, the event-time factor has the usual form $h_i(T_i\mid b_i;\theta)^{\delta_i}S_i(T_i\mid b_i;\theta)$, where $T_i$ is follow-up time, $\delta_i$ indicates an observed event, $h_i$ is the conditional hazard, and $S_i$ is the conditional survival function. The product form assumes the relevant conditional independence and censoring assumptions used to specify the likelihood. Any parameter dependence in any factor—including the random-effects density—contributes to the complete-data score.

Deriving the Score Identity

Assume that $f_i$ is differentiable in $\theta$, its derivative can be moved inside the integral, and $0<L_i(\theta)<\infty$. These regularity conditions justify differentiating the marginal likelihood integral. The differentiation rule and the identity $\partial f_i/\partial\theta=f_i\,\partial\log f_i/\partial\theta$ yield the result directly.

\begin{aligned} \frac{\partial \log L_i(\theta)}{\partial\theta} &=\frac{1}{L_i(\theta)} \frac{\partial}{\partial\theta}\int f_i(D_i,b_i;\theta)\,db_i\\ &=\frac{1}{L_i(\theta)} \int f_i(D_i,b_i;\theta) \frac{\partial\log f_i(D_i,b_i;\theta)}{\partial\theta}\,db_i\\ &=\int p(b_i\mid D_i;\theta)\,s_i^c(b_i;\theta)\,db_i\\ &=E\!\left[s_i^c(b_i;\theta)\mid D_i;\theta\right]. \end{aligned}

For a vector parameter, the equality holds component by component. For independent subjects, the sample log likelihood is $\ell(\theta)=\sum_i\log L_i(\theta)$, so its score is the sum of these subject-specific posterior expectations. This is an identity at the same parameter value $\theta$: the complete-data score and the posterior distribution must both be evaluated there.

Theorem: Under the differentiation and integrability conditions above, the observed-data score for subject $i$ is the posterior expectation of the complete-data score: $$ s_i(\theta)=\frac{\partial\log L_i(\theta)}{\partial\theta} =E\!\left[s_i^c(b_i;\theta)\mid D_i;\theta\right]. $$

Evaluating the Posterior Expectation by Quadrature

A quadrature rule replaces an integral by a weighted sum. Suppose a rule provides nodes $b_{ik}$ and nonnegative integration contributions $A_{ik}$ that include the likelihood and the appropriate random-effects integration weights. Then the marginal likelihood and its derivative numerator are approximated by $\sum_k A_{ik}$ and $\sum_k A_{ik}s_i^c(b_{ik};\theta)$, respectively. Dividing the two gives the numerical posterior average.

\widehat{s}_i(\theta) = \frac{\sum_{k=1}^{K} A_{ik}\,s_i^c(b_{ik};\theta)} {\sum_{k=1}^{K} A_{ik}}.

The normalized contributions $A_{ik}/\sum_j A_{ij}$ act as discrete posterior weights. If the random effects are Gaussian, quadrature may use a standard-normal transformation, a Cholesky transformation for correlated effects, or an adaptive rule centered near the posterior mass. In every case, the integration weights, transformation Jacobian, and random-effects density must be handled consistently. The numerator and denominator must approximate integrals over the same posterior target. A common multiplicative constant shared by all $A_{ik}$ cancels in the ratio.

Worked Example: A Scalar Survival-Model Score

Consider one hypothetical subject with scalar random effect $b\sim N(0,1)$. At the parameter value $\alpha=0$, suppose the part of the complete-data likelihood that depends on $b$ and $\alpha$ is $$ q_\alpha(b)=\exp\{(0.5+\alpha)b\} \exp\{-0.8\exp((0.5+\alpha)b)\}. $$ Other factors not depending on $b$ or $\alpha$ cancel from the posterior weights. Its complete-data score at $\alpha=0$ is $$ s_\alpha^c(b)= \left.\frac{\partial\log q_\alpha(b)}{\partial\alpha}\right|_{\alpha=0} =b\{1-0.8\exp(0.5b)\}. $$ Use the three-point standard-normal Gauss-Hermite rule, whose nodes are $-\sqrt{3},0,\sqrt{3}$ and whose normalized weights are $1/6,2/3,1/6$.

At the three nodes, direct substitution gives $q_0(-\sqrt{3})\approx0.300436$, $q_0(0)\approx0.449329$, and $q_0(\sqrt{3})\approx0.354896$. Multiplying by the quadrature weights gives the integration contributions $A_k\approx0.0500726,0.299553,0.0591493$, with sum $0.408775$. Thus the approximate posterior weights are $0.12249,0.73281,0.14470$.

The complete-data scores at the nodes are approximately $-1.1492,0,-1.5622$. Therefore the three-node marginal score is $$ 0.12249(-1.1492)+0.73281(0)+0.14470(-1.5622) \approx -0.3668. $$ The result is a score per subject with respect to $\alpha$, not an estimated value of $\alpha$. Its negative sign means the subject’s marginal log likelihood is locally decreasing as $\alpha$ increases at $\alpha=0$, under this hypothetical contribution.

Worked Example: A Gaussian Posterior Average

Suppose a hypothetical observation $y=1$ follows $Y\mid b\sim N(\beta+b,1)$, with $b\sim N(0,1)$, and evaluate the score for $\beta$ at $\beta=0$. The complete-data score is $s_\beta^c(b)=y-\beta-b=1-b$. Completing the square in the product of the prior and conditional likelihood gives $b\mid y,\beta=0\sim N(0.5,0.5)$.

Use a three-point standard-normal rule transformed to this posterior: the nodes are $0.5+\sqrt{0.5}(-\sqrt{3},0,\sqrt{3})$, or approximately $-0.724745,0.500000,1.724745$. The corresponding scores $1-b$ are $1.724745,0.500000,-0.724745$. Using weights $1/6,2/3,1/6$, $$ \frac{1}{6}(1.724745)+\frac{2}{3}(0.500000) +\frac{1}{6}(-0.724745)=0.500000. $$ The posterior average score is $0.5$ per unit of $\beta$. Here the transformed quadrature integrates the linear score exactly; this calculation illustrates the identity, not a claim that adaptive quadrature is exact in general.

Vector Scores and a Quadrature Check

A joint model usually has several parameters, including longitudinal fixed effects, survival parameters, random-effects covariance parameters, and association parameters. At each quadrature node, calculate the full complete-data score vector, then average each component with the same normalized posterior weights. The entries can differ in units; for example, a score for a parameter measured in years has inverse-year units, while a score for a dimensionless association parameter is dimensionless.

Worked Example: Averaging Two Score Components

For a hypothetical subject, suppose numerical integration gives unnormalized posterior contributions $A_1=0.12$, $A_2=0.30$, and $A_3=0.08$. Their sum is $0.50$, so the normalized weights are $0.24,0.60,0.16$. Suppose the complete-data score components at the same nodes are $s_\beta^c=(1,0,-2)$ and $s_\alpha^c=(-0.5,1,2)$. These are hypothetical quadrature outputs from a specified joint model.

For the first component, $$ \widehat{s}_\beta=0.24(1)+0.60(0)+0.16(-2)=-0.08. $$ For the second, $$ \widehat{s}_\alpha=0.24(-0.5)+0.60(1)+0.16(2)=0.80. $$ The approximate subject score is therefore $(-0.08,0.80)$. Each component uses the same posterior weights, but its own complete-data score values.

Diagnostics, Alternatives, and Limitations

The posterior-average form is often more convenient than differentiating a numerically computed marginal likelihood: it requires evaluating complete-data scores at integration nodes and averaging them. But it does not remove integration error. Increase the number of quadrature points, compare with an alternative suitable rule, or vary the adaptive integration settings. The marginal likelihood and each score component should stabilize to a tolerance appropriate for the optimization task. A stable likelihood alone is not enough; a score can remain sensitive to integration error.

A useful independent check is a central finite difference of the marginal log likelihood: $$ \frac{\ell(\theta+h e_j)-\ell(\theta-h e_j)}{2h}, $$ where $e_j$ selects parameter component $j$. Compare this with the integrated score component for several reasonable step sizes $h$. Large discrepancies can indicate an incorrect complete-data derivative, omitted parameter dependence in the random-effects density, inconsistent likelihood and score quadrature, or a finite-difference step that is too large or too small. Finite differences are a diagnostic, not automatically a better answer: they can be noisy when each likelihood evaluation has integration error.

If quadrature is applied after a parameter-dependent transformation, the exact integral is unchanged only when the transformation and its Jacobian are handled correctly. Differentiating an approximate sum whose nodes or weights change with $\theta$ can also introduce numerical effects not represented by simply averaging node scores. A robust implementation defines the target integral clearly, includes all density and Jacobian terms, and checks the resulting score against independent perturbations of the marginal likelihood.

Omitting the prior score: if a parameter changes the random-effects distribution, its density is part of $f_i$ and must be differentiated.
Using prior weights as posterior weights: the data likelihood changes the random-effects weights; multiply by the subject likelihood before normalizing.
Forgetting normalization: the score is a weighted average, not the unnormalized sum in its numerator.
Mixing parameter values: evaluate the complete-data score and posterior weights at the same $\theta$.

Direct numerical differentiation of the marginal likelihood is a practical alternative and an important validation check, but repeated likelihood evaluations may be costly and sensitive to integration error. The posterior-average calculation uses analytic complete-data derivatives when available. Its limitations are the usual ones for numerical integration: accuracy may deteriorate with higher-dimensional or sharply concentrated posterior distributions, and a finite quadrature rule only approximates the integral. A score computed this way is also not itself a standard error or an uncertainty interval.

Key takeaway: Integrate the complete-data score against the subject-specific random-effects posterior, using the same integration target and parameter value in the score and posterior weights. The normalized weighted sum is the marginal likelihood score.

Check Your Understanding

Answer each question, then compare your work with the solutions.

Why must the random-effects density be included in the complete-data score when its parameters are part of $\theta$?
Suppose two quadrature contributions are $0.2$ and $0.8$, with complete-data scores $3$ and $-1$. What is the approximate marginal score?
In the Gaussian example, what is the posterior mean of $b$, and how does it determine the score for $\beta$ at $\beta=0$?
Name two checks that can reveal a numerical or implementation problem in a computed marginal score.

Solutions. (1) The complete-data likelihood includes the random-effects density, so differentiating its log contributes to the score whenever that density depends on the parameter. (2) The contributions sum to $1$, so the score is $0.2(3)+0.8(-1)= -0.2$. (3) The posterior mean is $0.5$; since the complete-data score is $1-b$, its posterior expectation is $1-0.5=0.5$. (4) Increase or vary quadrature settings and check score stability; also compare the score with a central finite difference of the marginal log likelihood.