Assembling a Joint Likelihood From a Study Dataset

From a Study Record to a Likelihood Contribution

A joint model combines repeated measurements of a marker with a time-to-event outcome. To assemble its likelihood from a dataset, each subject’s record must contribute both the observed marker information and the event or censoring information. In a shared random-effects model, these two parts are linked through the same latent subject-specific random effects. Consequently, the observed-data contribution is not generally the product of a separately marginalized marker likelihood and a separately marginalized survival likelihood.

For subject \(i\), let \(Y_{ij}\) denote the marker value recorded at time \(t_{ij}\), for measurement occasions \(j\) in the observed set \(\mathcal O_i\). Let \(x_i\) and \(z_i\) denote the subject’s recorded covariates for the longitudinal and event models, respectively. Let \(T_i\) be the event time and \(C_i\) the censoring time. The observed follow-up time is \(U_i=\min(T_i,C_i)\), and the event indicator is \(\delta_i=I(T_i\le C_i)\): \(\delta_i=1\) means an event was observed at \(U_i\), and \(\delta_i=0\) means the observation was right-censored there. Write \(b_i\) for the subject’s random-effects vector and \(\theta\) for all model parameters.

Definition: A subject’s conditional likelihood contribution, \(L_i(b_i;\theta)\), is the joint density or probability of that subject’s observed marker values and observed event-time information, conditional on the random effects and covariates. The observed-data contribution is obtained by integrating \(L_i(b_i;\theta)\) against the random-effects density \(f_b(b_i;\theta)\).

For a standard basic joint model, denote the longitudinal conditional density at occasion \(j\) by \(f_Y(Y_{ij}\mid b_i,x_i,t_{ij};\theta)\), the event hazard by \(h_i(t\mid b_i,z_i;\theta)\), and its cumulative hazard through time \(u\) by \(H_i(u\mid b_i,z_i;\theta)=\int_0^u h_i(s\mid b_i,z_i;\theta)\,ds\). The hazard has units of inverse time; the cumulative hazard is dimensionless. Covariates may enter either component according to the specified model.

L_i(b_i;\theta) = \left\{\prod_{j\in\mathcal O_i} f_Y(Y_{ij}\mid b_i,x_i,t_{ij};\theta)\right\} h_i(U_i\mid b_i,z_i;\theta)^{\delta_i} \exp\{-H_i(U_i\mid b_i,z_i;\theta)\}.

The event part follows from the event-time density \(h_i(u\mid b_i,z_i;\theta)\exp\{-H_i(u\mid b_i,z_i;\theta)\}\) when an event occurs, and the survival probability \(\exp\{-H_i(u\mid b_i,z_i;\theta)\}\) when follow-up ends in right censoring. The exponent \(\delta_i\) selects the hazard factor only for an observed event. This expression omits the censoring-time density, which is appropriate when censoring is noninformative for the event process under the analysis and its distribution does not involve the parameters being estimated.

Assumptions Behind the Assembly

The factorization above is a model statement, not an automatic property of any dataset. It uses conditional independence: given \(b_i\) and the relevant covariates, the marker measurements have the specified longitudinal density, and the event process has the specified hazard; the marker and event components are conditionally independent apart from their shared dependence on \(b_i\). Within the longitudinal component, the product over \(j\) also presumes conditional independence of the measurement errors given the modeled trajectory and random effects. A model with residual serial correlation would replace that product by an appropriate joint longitudinal density.

The random-effects density \(f_b(b_i;\theta)\) must be specified, commonly as a multivariate normal density. Once the conditional likelihood is assembled, the subject’s observed-data contribution is

L_i(\theta) = \int L_i(b_i;\theta) f_b(b_i;\theta)\,db_i, \qquad L(\theta)=\prod_{i=1}^{n}L_i(\theta),

where the second equality assumes independent subjects conditional on the shared model parameters. The integral sums or averages the conditional contributions over the random-effects distribution, rather than choosing a single estimated value of \(b_i\). If marker values are missing, \(\mathcal O_i\) contains only the observed occasions. This observed-value product is appropriate when the missingness mechanism can be ignored for the likelihood analysis, for example under an appropriate missing-at-random assumption conditional on included information. A nonignorable missingness mechanism requires modeling the missingness process as well.

Why the Integral Keeps the Components Joint

The longitudinal and event terms both depend on \(b_i\). That dependence makes them informative about the same latent subject characteristics. Conditional on a particular \(b_i\), their contributions multiply; marginally, they must be integrated together. In general, an integral of a product is not the product of the integrals:

\int L_{Yi}(b_i)L_{Ti}(b_i)f_b(b_i)\,db_i \ne \left\{\int L_{Yi}(b_i)f_b(b_i)\,db_i\right\} \left\{\int L_{Ti}(b_i)f_b(b_i)\,db_i\right\}.

The right-hand side effectively averages the two components separately and loses their shared mixing variable. Likewise, replacing \(b_i\) by an estimated value \(\hat b_i\) gives a conditional plug-in contribution \(L_i(\hat b_i;\theta)\), not the integrated observed-data likelihood. These are different operations with different statistical meanings.

A short derivation makes the event factor’s role explicit. Conditional on \(b_i\), survival through time \(u\) is \(S_i(u\mid b_i)=\exp\{-H_i(u\mid b_i)\}\). Differentiating gives the event-time density:

-\frac{d}{du}S_i(u\mid b_i) = h_i(u\mid b_i)\exp\{-H_i(u\mid b_i)\}.

Thus an event observed at \(u\) contributes this density, while censoring at \(u\) contributes \(S_i(u\mid b_i)\). Multiplying by the observed longitudinal densities gives the conditional subject likelihood; averaging that entire product over \(f_b\) gives the observed-data contribution.

Worked Example: Translating a Mixed Dataset Record

Suppose a hypothetical subject has covariates \(x_i\) and \(z_i\), marker values \(Y_{i1}=4.8\) at \(t_{i1}=0\) years and \(Y_{i2}=5.1\) at \(t_{i2}=1\) year, and no marker recorded at a scheduled 2-year visit. The subject has an observed event at \(U_i=2.4\) years, so \(\delta_i=1\). Suppose the model specifies conditionally independent normal marker errors with standard deviation \(\sigma=0.45\) around mean \(m_i(t;b_i,x_i)\).

The observed set is \(\mathcal O_i=\{1,2\}\); the missing 2-year marker does not contribute a normal density under an ignorable missingness assumption. For each observed value, the normal density is

f_Y(Y_{ij}\mid b_i,x_i,t_{ij}) = \frac{1}{0.45\sqrt{2\pi}} \exp\left\{-\frac{[Y_{ij}-m_i(t_{ij};b_i,x_i)]^2}{2(0.45)^2}\right\}.

The event indicator is one, so this record contributes a hazard at 2.4 years as well as survival through that time. Its full conditional contribution is

L_i(b_i;\theta) = f_Y(4.8\mid b_i,x_i,0) f_Y(5.1\mid b_i,x_i,1) h_i(2.4\mid b_i,z_i) \exp\{-H_i(2.4\mid b_i,z_i)\}.

The event factor has units of inverse years, while the two marker densities each have units inverse to the marker’s measurement unit. The resulting likelihood contribution is a density in the observed continuous outcomes. To obtain the observed-data contribution, integrate this entire expression over \(f_b(b_i;\theta)\). The absent 2-year value is not imputed as if it had been observed.

Numerical Contributions at Fixed Random Effects

Worked Example: An Observed Event and Two Marker Densities

Consider a hypothetical subject with two marker residuals \(0.0526\) and \(0.2572\) marker units, conditional residual standard deviation \(0.45\) marker units, and an event at \(U_i=0.0481\) years. At the random-intercept value \(b_i=-0.2\), suppose the hazard is constant through the event time and equals \(0.15\exp(b_i)\) per year. Assume the two marker observations are conditionally independent normal measurements.

For a normal residual \(r\), the density is \(\exp\{-r^2/[2(0.45)^2]\}/[0.45\sqrt{2\pi}]\). Substitution gives densities approximately \(0.8805\) and \(0.7529\) per marker unit, so their product is approximately \(0.6630\) per squared marker unit.

The event hazard at this random-effects value is

h_i=0.15\exp(-0.2)=0.122810\ \text{per year}.

Because the hazard is constant over this interval, the cumulative hazard through the event time is \(H_i=0.1228096(0.0481)\approx0.005907\). The event-time factor is therefore

h_i\exp(-H_i) \approx 0.1228096\exp(-0.005907) \approx 0.12209\ \text{per year}.

Thus the conditional subject contribution is approximately \(0.6630(0.12209)=0.08094\), in the corresponding joint density units. If \(b_i\sim N(0,1)\), its density at \(-0.2\) is \(\phi(-0.2)\approx0.3910\). The integrand \(L_i(b_i;\theta)f_b(b_i;\theta)\) at this point is approximately \(0.08094(0.3910)=0.03165\). This is one value of the integrand, not the subject’s integrated likelihood.

Worked Example: A Right-Censored Subject

For another hypothetical subject, suppose the same two residuals and residual standard deviation apply, but follow-up ends by right censoring at \(U_i=1.315\) years. Let \(\delta_i=0\). At \(b_i=-0.3\), assume a constant hazard \(0.2\exp(b_i)\) per year. The marker-density product remains approximately \(0.6630\) per squared marker unit.

The hazard and cumulative hazard are

h_i=0.2\exp(-0.3)\approx0.148164\ \text{per year}, \qquad H_i=0.148163644(1.315)\approx0.19484.

Since \(\delta_i=0\), there is no hazard multiplier at the censoring time. The event component is the survival factor:

\exp(-H_i)=\exp(-0.1948352)\approx0.8230.

The conditional contribution is therefore \(0.6629668(0.8229709)\approx0.5456\), with marker-density units. For \(b_i\sim N(0,1)\), \(\phi(-0.3)\approx0.381388\); the integrand at this random-effects value is approximately \(0.5456(0.381388)=0.2081\). As in the event example, this pointwise value must be integrated over all possible random effects to give the observed-data likelihood.

Worked Example: Approximating a Random-Effects Integral

To see what integration does numerically, consider a deliberately simplified hypothetical subject with no marker observations, right-censoring at one year, and conditional hazard \(h(t\mid b)=\exp(b)\) per year. Let \(B\sim N(0,1)\). The observed-data contribution is the expectation of the conditional survival probability:

L_i=\int_{-\infty}^{\infty}\exp\{-\exp(b)\}\phi(b)\,db.

A three-node Gauss–Hermite approximation for a standard normal expectation uses nodes \(0,\sqrt{3},-\sqrt{3}\) and weights \(2/3,1/6,1/6\). Substituting the function \(g(b)=\exp\{-\exp(b)\}\) gives

L_i\approx \frac{2}{3}\exp(-1) +\frac{1}{6}\exp\{-\exp(\sqrt{3})\} +\frac{1}{6}\exp\{-\exp(-\sqrt{3})\} \approx 0.3855.

The terms are approximately \(0.2453\), \(0.000585\), and \(0.1396\), respectively. Their sum is the quadrature approximation, not an exact evaluation of the integral; more nodes can improve accuracy. In this simple example, the survival probability is averaged over the random-effects distribution. Evaluating it only at \(b=0\) would give \(\exp(-1)\approx0.3679\), which is not the same calculation.

Checks, Alternatives, and Limitations

Before evaluating or optimizing a likelihood, check that each record has been translated consistently. In particular, confirm the time units for \(U_i\), \(h_i\), and \(H_i\); verify that \(\delta_i=1\) only for observed events; and ensure that only observed marker values appear in the longitudinal product. A cumulative hazard must be nonnegative, and the event component must be \(h_i(U_i)\exp\{-H_i(U_i)\}\) for an event but \(\exp\{-H_i(U_i)\}\) for right censoring. For numerical integration, inspect convergence as the number of quadrature nodes or integration accuracy increases.

The integral is often not available in closed form. Gaussian quadrature or Monte Carlo integration can approximate it, while a plug-in random effect or a first-order approximation may be computationally simpler but does not equal exact marginalization. A model with no shared random effects has a different factorization and does not represent the same source of dependence. The appropriate choice depends on the scientific model and computational demands; none of these alternatives repairs a misspecified hazard, longitudinal distribution, censoring assumption, or missingness mechanism.

A common mistake is to multiply a marker likelihood averaged over random effects by a survival likelihood averaged over random effects, or to insert a subject’s estimated random effect into both components and call the result marginal. The correction is to form the joint conditional product first and integrate that product against one shared random-effects density. The assembly also relies on the stated conditional-independence and censoring assumptions; if those do not fit the data-generating process, changing the integration method alone will not fix the likelihood.

Check Your Understanding

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

A subject has two observed markers, a third marker is missing, and the subject is right-censored at \(u\). Which marker and event factors appear in the conditional contribution?
Why can the marker and survival likelihoods not generally be integrated separately over the same random effects?
For a subject with an event indicator \(\delta_i=1\), what event factor appears at observed time \(U_i\)? What changes when \(\delta_i=0\)?
What does a single integrand value \(L_i(b_i;\theta)f_b(b_i;\theta)\) represent, and what additional operation gives the subject’s observed-data contribution?

Solutions. (1) Include the two observed marker densities, omit the missing marker density under an ignorable missingness assumption, and multiply by \(\exp\{-H_i(u\mid b_i,z_i)\}\); right censoring contributes no hazard multiplier. (2) Both components depend on the same \(b_i\), so the observed-data contribution averages their product over one shared random-effects value. In general, the integral of a product is not the product of separate integrals. (3) For \(\delta_i=1\), use \(h_i(U_i\mid b_i,z_i)\exp\{-H_i(U_i\mid b_i,z_i)\}\); for \(\delta_i=0\), use only \(\exp\{-H_i(U_i\mid b_i,z_i)\}\). (4) It is the integrand’s value at one particular random-effects vector. Integrate it over the random-effects space to obtain \(L_i(\theta)\).