Encoding Joint Likelihoods in INLA

Why joint-likelihood encoding requires careful indexing

A joint model can be statistically well specified and still be encoded incorrectly. The danger is especially clear when a subject contributes several longitudinal measurements and several survival intervals: the data rows do not correspond one-to-one, but both components must refer to the same subject-level latent effects. A shifted index can silently attach one subject’s survival information to another subject’s longitudinal trajectory.

INLA represents models through likelihood contributions and a latent Gaussian field. For this tutorial, the essential task is to define the two likelihood components, show how their linear predictors use common subject effects, and ensure that each observation is mapped to the correct subject. The earlier material on latent Gaussian representations and piecewise survival likelihoods supplies the model ingredients; here the emphasis is on combining and indexing them.

Specify the likelihood components before arranging the data

For a Gaussian longitudinal model, a useful example is a random intercept and random slope:

Here \(t_{ij}\) is measurement time, and \(b_{0i}\) and \(b_{1i}\) are subject-specific departures in intercept and slope. Their joint distribution is Gaussian, with a specified covariance structure. For a piecewise-constant hazard, write the interval-specific hazard as

The \(w_i\) are survival covariates, \(\lambda_{0k}\) is the baseline hazard in interval \(k\), and \(\alpha_0,\alpha_1\) specify how the shared effects enter the survival predictor. This is a shared-random-effects association specification; it is not automatically identical to associating the hazard with the entire longitudinal mean trajectory. The survival contribution over interval \(k\), as a likelihood in the model parameters, is proportional to

A Poisson representation uses \(d_{ik}\) with mean \(\mu_{ik}\), so that \(\log(\mu_{ik})=\log(E_{ik})+\log(\lambda_{ik})\). The Poisson probability mass is proportional to the piecewise-exponential survival likelihood, not generally equal to it: their ratio includes the data-only factor \(E_{ik}^{d_{ik}}/d_{ik}!\). Since this factor does not depend on model parameters, it does not change likelihood-based parameter inference. It must not, however, be described as equality of the two probability contributions.

An INLA encoding therefore needs the right response for each likelihood, its own predictor structure, and the right subject index wherever a shared effect is used. A missing longitudinal response on a survival row is structural, not a value to be imputed; similarly, a survival response is absent on a longitudinal row.

How the common latent field enters an INLA data layout

A practical layout keeps track of the likelihood and subject for every observation. One common strategy is a combined response with a separate column for each likelihood: the longitudinal column contains \(y_{ij}\) on longitudinal rows and is missing on survival rows, while the survival column contains \(d_{ik}\) on survival rows and is missing on longitudinal rows. The likelihood family assigned to each column must match its response. In a stacked construction, the same principle applies even if separate stacks are assembled first: each row’s response, covariates, exposure, and indices must remain together.

For the survival rows, the linear predictor for the Poisson mean includes \(\log(E_{ik})\) as an offset. Equivalently, the hazard predictor is \(\log(\lambda_{0k})+w_i^\mathsf{T}\gamma+\alpha_0b_{0i}+\alpha_1b_{1i}\), with the exposure added to form the log Poisson mean. Rows with zero exposure should not be included as event-time information; an event row must have positive exposure.

The model formula or equivalent design matrices must make the subject mapping explicit. Separate likelihood components may require separate predictor-specific covariates or design matrices. A reliable implementation should not assume that a variable’s position in one data frame automatically indexes the correct latent effect in another. If a shared effect is expressed through a copied or otherwise scaled latent component in an INLA construction, verify that the copy refers to the intended underlying field and that its scale enters the intended predictor. The exact syntax depends on the chosen INLA stack and model specification; the statistical requirement does not.

Justification, diagnostics, and alternatives

The conditional joint likelihood follows from the stated conditional independence assumption. Given the latent effects and parameters, each longitudinal measurement contributes its Gaussian density, and each survival interval contributes its piecewise-exponential likelihood factor. Therefore, apart from factors depending only on observed survival exposures and event indicators, the joint likelihood is the product

where \(f_{\mathrm{G}}\) is the Gaussian density and \(\theta\) denotes model parameters. This product is the reason separate likelihood components can represent a joint model: the shared latent effects connect their predictors, while each observed datum contributes through its appropriate sampling model. It is a likelihood conditional on the latent effects, not a claim that the two data sources are marginally independent after integrating over those effects.

Several checks catch common encoding failures before interpreting an approximation:

• Subject-index audit. For each row, compare the subject identifier with the index supplied to every shared latent effect. Check both the first and last subject and at least one subject with multiple visits and intervals.
• Row and response audit. Count observed longitudinal responses and event-interval rows separately. Confirm that missing entries occur only in the response column for the other likelihood, rather than being accidental missing data.
• Exposure audit. Check that interval exposures are nonnegative, have the time units used by the hazard, and sum to each subject’s observed follow-up. Confirm that \(\log(E_{ik})\) is used as the offset where required.
• Predictor audit. Evaluate a few predictors by hand at fixed latent-effect values. A shared effect should change both intended components, and should not change predictors for unrelated subjects.
• Likelihood audit. For a small data set, calculate selected Gaussian and survival contributions directly. Compare the Poisson and piecewise-exponential contributions allowing for the exposure- and event-dependent proportionality factor.

A common mistake is to reuse the longitudinal row number as the survival subject index. That is wrong because visit rows and interval rows have different lengths and ordering. The correction is to build each index from the subject identifier using the same subject-to-index map, then check the map against the data rows. Another mistake is to call the Poisson mass exactly equal to the event-time probability contribution. It is instead proportional to the piecewise-exponential likelihood as a function of model parameters, with a data-only factor.

A combined multi-likelihood response can be concise, but separate stacks or explicitly constructed design matrices may make predictor-specific structure easier to inspect. Either can encode the same model if the response-to-family assignments, offsets, and shared-effect indices agree. Fitting separate longitudinal and survival models is not an equivalent alternative: it does not generally retain the dependence induced by shared latent effects. A two-stage analysis can be useful for different purposes, but it does not encode this same joint likelihood. INLA also has limits: the model must fit the supported latent-Gaussian and likelihood structure, and incorrect data assembly is not repaired by numerical approximation. More complex association structures may require additional modeling machinery rather than a simple shared index.