Worked Likelihood for a Current-Value Random-Slope Model

Why the Random Slope Changes the Likelihood

In a shared-random-effects joint model, longitudinal measurements and event times are conditionally connected through subject-specific latent effects. With a random intercept alone, the latent biomarker trajectory shifts up or down but does not change its rate of change. A random slope allows both the level and the time trend to vary by subject. When the event hazard depends on the subject’s current fitted biomarker value, that changing trajectory enters the hazard and its cumulative integral. The likelihood must therefore use the same random-effects value consistently in the longitudinal and survival components.

This tutorial derives that likelihood for a Gaussian longitudinal model and a Weibull baseline hazard. The earlier material on the random-intercept Weibull model establishes the basic survival contribution; here the key additional step is evaluating the current-value association along a subject-specific linear trajectory.

Model and Assumptions

Definition: For subject $i$, let $y_{ij}$ be the longitudinal measurement at time $t_{ij}$, for $j=1,\ldots,n_i$. Let $T_i$ be the observed event or censoring time, and let $\delta_i=1$ indicate an observed event and $\delta_i=0$ right censoring. The longitudinal model is $$ y_{ij}=x_{ij}^{T}\beta+z_{ij}^{T}b_i+\epsilon_{ij}, \qquad b_i=\begin{pmatrix}b_{0i}\\b_{1i}\end{pmatrix}\sim N_2(0,D), \qquad \epsilon_{ij}\sim N(0,\sigma^2). $$ For a linear random-intercept-and-slope trajectory, $z_{ij}=(1,t_{ij})^T$, so the fitted current value is $m_i(t;b_i)=\beta_0+\beta_1t+b_{0i}+b_{1i}t$. Conditional on $b_i$, the longitudinal errors are independent. The survival hazard is $h_i(t\mid b_i)=h_0(t)\exp\{\gamma^Tw_i+\alpha m_i(t;b_i)\}$, where $w_i$ is a vector of survival covariates and $\alpha$ is the current-value association parameter. For a Weibull baseline, $h_0(t)=\lambda\rho t^{\rho-1}$, with cumulative baseline hazard $H_0(t)=\lambda t^\rho$, for $\lambda>0$, $\rho>0$, and $t>0$.

The random-effects covariance matrix $D$ must be positive definite, and $\sigma^2>0$. The event process is assumed conditionally independent of the longitudinal measurement errors given $b_i$ and the included covariates. Censoring is assumed non-informative conditional on the modeled information, so the censoring mechanism contributes no additional parameter-dependent likelihood term. The visit times are treated as fixed or as an ignorable process. Measurements are modeled as Gaussian with constant residual variance; dependence not captured by the random effects would require a different longitudinal model.

The current value $m_i(t;b_i)$ is the latent mean trajectory, not the noisy observed measurement. Since $\alpha m_i(t;b_i)$ appears in an exponential, its units must be dimensionless: if the biomarker is measured in units $U$, $\alpha$ has units $U^{-1}$. The baseline hazard and cumulative hazard must use the same time unit as $t$.

Deriving the Conditional and Marginal Likelihood

Given $b_i$, the longitudinal density is a product of normal densities. The event-time contribution follows from the hazard and cumulative hazard. Specifically, the conditional survival function is $S_i(t\mid b_i)=\exp\{-H_i(t\mid b_i)\}$, where $H_i(t\mid b_i)=\int_0^t h_i(u\mid b_i)\,du$. An observed event contributes its conditional hazard at $T_i$ multiplied by survival to $T_i$; a right-censored subject contributes survival to $T_i$ only. These cases can be written together using $\delta_i$.

L_i(b_i)= \left\{\prod_{j=1}^{n_i} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\!\left[-\frac{(y_{ij}-m_i(t_{ij};b_i))^2}{2\sigma^2}\right]\right\} \{h_i(T_i\mid b_i)\}^{\delta_i} \exp\{-H_i(T_i\mid b_i)\}.

The subject’s marginal likelihood contribution is the conditional likelihood averaged over the random-effects distribution. This integration is essential: treating an estimated or most likely $b_i$ as if it were known does not give the marginal likelihood.

L_i=\int_{\mathbb R^2} L_i(b_i)\, \phi_2(b_i;0,D)\,db_i, \qquad \ell_i=\log L_i, \qquad \ell=\sum_i\ell_i,

where $\phi_2(\,\cdot\,;0,D)$ is the bivariate normal density. For the Weibull baseline, the cumulative hazard is

H_i(T_i\mid b_i)= \int_0^{T_i}\lambda\rho u^{\rho-1} \exp\{\gamma^Tw_i+\alpha m_i(u;b_i)\}\,du.

Unlike the random-intercept case, the exponent includes a time-varying linear trajectory when $b_{1i}$ or $\beta_1$ is nonzero. Except in special cases, this integral has no elementary closed form and must be evaluated numerically for each random-effects value considered during likelihood computation.

Worked Example: Evaluating the Survival Term Along a Trajectory

Consider a hypothetical subject with $\beta_0=1.8$, $\beta_1=0.12$ biomarker units per year, $b_{0i}=0.2$ biomarker units, and $b_{1i}=-0.05$ biomarker units per year. There are no survival covariates. Let $\alpha=0.6$ per biomarker unit, $\lambda=0.04$ per year$^{1.5}$, $\rho=1.5$, and suppose an event occurs at $T_i=4$ years.

The trajectory at any time $t$ is $m_i(t;b_i)=1.8+0.12t+0.2-0.05t=2.0+0.07t$. Thus its value at the event time is $m_i(4;b_i)=2.28$ biomarker units. The conditional event hazard is

h_i(4\mid b_i)=0.04(1.5)4^{0.5}\exp(0.6\times2.28) =0.12\exp(1.368)\approx0.4712\ \text{per year}.

For this trajectory, $H_i(4\mid b_i)=\int_0^4 0.06u^{0.5}\exp\{0.6(2+0.07u)\}\,du$. This integral can be evaluated by numerical quadrature; it is approximately $1.1762$. The conditional event-time density contribution is therefore about $0.4712\exp(-1.1762)=0.1454$ per year. The event density is not just the hazard: it includes the probability of surviving event-free up to year 4.

Combining Longitudinal and Survival Contributions

A crucial bookkeeping rule is that every term for a subject at a fixed $b_i$ must be evaluated using that same trajectory. The longitudinal means at the visit times, the hazard at the event or censoring time, and the cumulative hazard along the interval all depend on the same $b_{0i}$ and $b_{1i}$. The following example makes this consistency visible.

Worked Example: One Conditional Subject Likelihood

Suppose a hypothetical subject has observations $y_{i1}=2.3$ at $t_{i1}=0$ years and $y_{i2}=2.1$ at $t_{i2}=2$ years. Use $\beta_0=2.0$, $\beta_1=-0.2$ biomarker units per year, $b_{0i}=0.3$, and $b_{1i}=0.1$ biomarker units per year. Let $\sigma=0.4$ biomarker units. The subject has an event at year 3. For the survival model, take $\lambda=0.05$ per year$^{1.4}$, $\rho=1.4$, $\alpha=0.5$ per biomarker unit, and no survival covariates.

The single trajectory is $m_i(t;b_i)=2.0+0.3+(-0.2+0.1)t=2.3-0.1t$. It gives $m_i(0;b_i)=2.3$, $m_i(2;b_i)=2.1$, and $m_i(3;b_i)=2.0$. Therefore both observed longitudinal residuals are zero, and their joint conditional density is

\left(\frac{1}{\sqrt{2\pi}(0.4)}\right)^2 \exp\!\left[-\frac{0^2+0^2}{2(0.4)^2}\right] \approx(0.997356)^2\approx0.9947.

At year 3, the hazard is $0.05(1.4)3^{0.4}\exp(0.5\times2.0)\approx0.2953$ per year. The cumulative hazard must use $m_i(u;b_i)=2.3-0.1u$ throughout $0\leq u\leq3$:

H_i(3\mid b_i)=\int_0^3 0.07u^{0.4}\exp\{0.5(2.3-0.1u)\}\,du \approx0.6741.

Thus the conditional event-time density is approximately $0.2953\exp(-0.6741)=0.1505$ per year. Multiplying by the longitudinal density gives $L_i(b_i)\approx0.9947(0.1505)=0.1497$, with units determined by the two longitudinal densities and the event-time density. This is a conditional contribution at the specified $b_i$; it is not yet the marginal contribution for the subject.

Integrating Over the Random Effects

The marginal likelihood averages contributions such as the one just calculated over plausible random intercepts and slopes under $N_2(0,D)$. Numerical quadrature or Monte Carlo methods can approximate this integral. As a compact illustration, a two-node product quadrature will be used below. It is deliberately coarse: the point is to show exactly how the conditional contributions are combined, not to recommend two nodes for a final fit.

Worked Example: A Product-Quadrature Approximation

Use the same fixed effects, residual standard deviation, Weibull parameters, and association parameter as in the preceding example, but suppose the observed values are $y_{i1}=2.1$ at year 0 and $y_{i2}=1.5$ at year 2. Let $b_{0i}\sim N(0,0.3^2)$ and $b_{1i}\sim N(0,0.1^2)$, independently. Approximate the expectation over the random effects with nodes $b_{0i}=\pm0.3$, $b_{1i}=\pm0.1$, each receiving product weight $1/4$. These are the two-point standard-normal nodes $\pm1$, rescaled by the respective standard deviations.

At the four nodes, the trajectory values at years $0$ and $2$, respectively, are $(2.3,2.1)$, $(2.3,1.7)$, $(1.7,1.5)$, and $(1.7,1.1)$. The corresponding values at year 3 are $2.0,\ 1.4,\ 1.4,$ and $0.8$. Thus each survival contribution uses the same trajectory that produced its two longitudinal means.

The longitudinal density at each node is $0.9974^2$ times the exponential of minus half the sum of squared standardized residuals. The four values are approximately $0.2850,\ 0.7747,\ 0.6030,$ and $0.3659$. Using the Weibull hazard and cumulative-hazard formulas along each node’s own trajectory, the conditional event-time density values at year 3 are approximately $0.1505,\ 0.1238,\ 0.1328,$ and $0.1063$ per year. Their products with the longitudinal densities are approximately $0.04289,\ 0.09591,\ 0.08008,$ and $0.03890$. The quadrature approximation to the marginal subject likelihood is therefore

L_i\approx\frac{0.04289+0.09591+0.08008+0.03890}{4} \approx0.0644.

This calculation illustrates why substituting a single estimated random-effects value is not equivalent to integration: the likelihood is an average of the joint contributions, not the joint contribution at the average or best-fitting node. In an actual analysis, quadrature resolution should be sufficient for the integral and the resulting estimates to be stable.

Censoring, Diagnostics, and Alternatives

If the subject in the second example were right-censored at year 3 rather than having an observed event, then $\delta_i=0$. The survival contribution would be $\exp\{-H_i(3\mid b_i)\}$, with no hazard multiplier. Censoring does not mean the subject had no survival information: survival through the censoring time still contributes. Conversely, treating a censored subject as if an event occurred inserts a hazard factor that is not supported by the observed data.

Several checks help detect a faulty implementation. For selected random-effects values, calculate the trajectory at every visit and at the event time independently, then confirm that those values match the inputs used in the longitudinal and survival terms. Check that $H_i(t\mid b_i)\geq0$, is nondecreasing in $t$, and equals zero at time zero. For $\alpha=0$, the survival component should no longer depend on the biomarker trajectory. For $\delta_i=0$, the hazard-at-$T_i$ factor must disappear. Compare numerical integration at increasing resolutions; a substantial shift in the likelihood or estimates signals inadequate integration or numerical instability. A very small marginal likelihood can be valid, but overflow in $\exp\{\alpha m_i(t;b_i)\}$ or underflow in $\exp\{-H_i(t\mid b_i)\}$ calls for stable log-scale computation.

A random-intercept model is a simpler alternative when subject-specific departures in trajectory slope are not needed; it imposes $b_{1i}=0$. That restriction can be inadequate if individual trends vary, but an unnecessary random slope can be difficult to estimate with sparse or short follow-up. A current-value association also differs from linking the hazard to the current slope or to a cumulative history of the biomarker: those are different scientific models, not interchangeable computational shortcuts. The current-value model’s limitation is that it summarizes the trajectory’s relationship with event risk through its value at time $t$, and it relies on the chosen longitudinal mean structure and baseline hazard being adequate.

Common Mistake

Common pitfall: Using one set of random effects to compute the longitudinal means and a different, inconsistent trajectory to compute the survival contribution. For example, if a model gives $m_i(t)=1.5$ throughout follow-up, its fitted means at every visit must also be $1.5$; it cannot simultaneously produce visit means $2.30$ and $2.25$. The correction is to define $m_i(t;b_i)$ once and evaluate it at all visit and survival times using the same $b_i$. Also do not replace the required integral over $b_i$ by the likelihood at a single random-effects estimate.

Check Your Understanding

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

If $\delta_i=0$, which part of the conditional survival contribution remains?
For $m_i(t)=1.2+0.3t$, what is the fitted current value at $t=2$, and what value must be used in the longitudinal mean at a visit at $t=2$?
When $\alpha=0$, which part of the conditional likelihood becomes independent of $b_i$ through the biomarker association?
Why does evaluating $L_i(b_i)$ at one estimated random-effects value not generally equal the marginal likelihood?

Solutions. (1) The survival contribution is $\exp\{-H_i(T_i\mid b_i)\}$; there is no hazard multiplier. (2) The current value is $1.2+0.3(2)=1.8$, and the fitted longitudinal mean uses $1.8$ at that visit, before comparing it with the observed measurement. (3) The survival hazard and cumulative hazard no longer depend on the current biomarker value through $\alpha m_i(t;b_i)$; the longitudinal likelihood can still depend on $b_i$. (4) The marginal likelihood averages the conditional likelihood over the full Gaussian distribution of the random effects, and in general an average of a function is not its value at one estimated point.