Worked Likelihood for a Random-Intercept Weibull Model

Why the Random Intercept Must Be Integrated Out

In a shared-random-effects joint model, a subject’s biomarker measurements and event time are linked through a latent subject-specific effect. The observed-data likelihood must account for uncertainty about that effect. It is therefore not obtained by estimating a separate intercept for each subject and then treating that estimate as known: it averages the conditional likelihood over the random-intercept distribution.

This tutorial works through that marginal likelihood for Gaussian biomarker measurements and a Weibull event hazard. The focus is the likelihood contribution and its numerical evaluation at specified parameter values, not the full optimization of all model parameters.

Definition: For subject \(i\), let \(y_{ij}\) be the biomarker 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, while \(\delta_i=0\) indicates right censoring. The shared random intercept is \(b_i\sim N(0,\tau^2)\). Conditional on \(b_i\), the biomarker model is \(y_{ij}=\beta_0+\beta_1t_{ij}+b_i+\epsilon_{ij}\), where the errors are independent \(N(0,\sigma^2)\). The latent mean trajectory is \(m_i(t;b_i)=\beta_0+\beta_1t+b_i\).

The event hazard conditional on \(b_i\) is specified as

h_i(t\mid b_i)=\kappa\rho t^{\rho-1} \exp\{\gamma z_i+\alpha m_i(t;b_i)\},\qquad \kappa>0,\quad \rho>0.

Here \(\kappa\rho t^{\rho-1}\) is the Weibull baseline hazard, \(z_i\) is a fixed covariate, \(\gamma\) is its log-hazard coefficient, and \(\alpha\) is the association parameter for the latent biomarker level. The cumulative hazard is

H_i(T_i\mid b_i)=\int_0^{T_i} h_i(u\mid b_i)\,du.

When \(\beta_1\ne 0\), the association term changes over time, and this integral generally needs numerical evaluation. If \(\beta_1=0\), then \(m_i(t;b_i)=\beta_0+b_i\) is constant and the cumulative hazard has the simpler form \(\kappa T_i^\rho\exp\{\gamma z_i+\alpha(\beta_0+b_i)\}\).

Constructing the Marginal Likelihood

The likelihood construction uses three assumptions. First, conditional on \(b_i\), the longitudinal errors are independent with the stated Gaussian distribution. Second, conditional on \(b_i\) and the covariates, the longitudinal measurements and event process are independent. Third, subjects are independent, with their random intercepts independently drawn from the same \(N(0,\tau^2)\) distribution. For right censoring, censoring must be non-informative for the event process conditional on the variables in the model; its likelihood contribution can then be omitted.

For a fixed value \(b\), the longitudinal density is

f_i(\mathbf y_i\mid b)= \prod_{j=1}^{n_i} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{[y_{ij}-(\beta_0+\beta_1t_{ij}+b)]^2}{2\sigma^2}\right\}.

For an observed event at \(T_i\), the conditional time-to-event density is the hazard at that time multiplied by the conditional survival probability. For a censored observation, only the survival probability enters. Both cases are captured by \(h_i(T_i\mid b)^{\delta_i}\exp\{-H_i(T_i\mid b)\}\). Multiplying by the random-intercept density gives the conditional-data integrand.

L_i(\theta)= \int_{-\infty}^{\infty} f_i(\mathbf y_i\mid b)\, h_i(T_i\mid b)^{\delta_i} \exp\{-H_i(T_i\mid b)\} \frac{1}{\sqrt{2\pi\tau^2}} \exp\left(-\frac{b^2}{2\tau^2}\right)\,db.

The observed-data likelihood is \(L(\theta)=\prod_i L_i(\theta)\), and its log is \(\ell(\theta)=\sum_i\log L_i(\theta)\). The integral is essential: it averages the full conditional contribution, including both data sources, over plausible values of \(b\). In general, integrating the longitudinal and event components separately and multiplying their marginal likelihoods would discard their dependence through the shared effect.

Justification: The density of the observed data for subject \(i\) is obtained by conditioning on the unobserved \(b_i\), multiplying the conditional longitudinal and event contributions under conditional independence, and then integrating with respect to the density of \(b_i\). Thus \(L_i(\theta)=\int f(\mathbf y_i,T_i,\delta_i\mid b)\,f(b)\,db\). This is the law of total probability applied to a continuous latent variable.

Worked Example: Writing a Censored Subject’s Contribution

Suppose a hypothetical subject has a single measurement \(y=1.0\) at time \(t=0\), is censored at \(T=1\), and has \(z=0\). Consider the parameter values \(\beta_0=1\), \(\beta_1=0\), \(\sigma^2=0.25\), \(\tau=0.5\), \(\kappa=0.1\), \(\rho=1.5\), and \(\alpha=0.5\). Since \(z=0\), the term \(\gamma z\) is zero. The observed event indicator is \(\delta=0\).

At a given \(b\), the measurement density is \(f(y=1\mid b)=(1/\sqrt{2\pi(0.25)})\exp\{-(1-(1+b))^2/[2(0.25)]\} =0.7979\exp(-2b^2)\). Because \(\beta_1=0\), the cumulative hazard at time 1 is \(H(1\mid b)=0.1(1^{1.5})\exp\{0.5(1+b)\}=0.1\exp\{0.5(1+b)\}\).

For a censored observation, the conditional event contribution is the survival probability, \(\exp\{-H(1\mid b)\}\), with no hazard multiplier. At \(b=0\), \(H=0.1\exp(0.5)=0.1649\), so the integrand excluding the random-effect density is \(0.7979\exp(-0.1649)=0.6766\). This is not yet the marginal likelihood contribution: it is the conditional-data contribution at one particular random-intercept value. It must still be averaged over \(b\sim N(0,0.25)\).

Numerical Evaluation by Gauss-Hermite Quadrature

The integrand is usually not available in a form that can be integrated analytically. A direct numerical approach is to transform the integral into an expectation over the normal random effect and approximate that expectation by weighted evaluations of the integrand. Three-point Gauss-Hermite quadrature for \(b\sim N(0,\tau^2)\) uses nodes \(0,\pm\sqrt{3}\tau\) and weights \(2/3,1/6,1/6\). Thus, for a function \(q\),

E\{q(b)\}\ \approx\ \frac{2}{3}q(0)+\frac{1}{6}q(\sqrt{3}\tau)+\frac{1}{6}q(-\sqrt{3}\tau).

The approximation is exact when \(q\) is a polynomial of degree at most five, but the likelihood integrand here is not generally such a polynomial. The quadrature is an approximation whose adequacy should be checked by increasing the number of nodes when computing final estimates.

Worked Example: Marginal Likelihood for the Censored Subject

Continue with the hypothetical subject and parameters above. Here \(\tau=0.5\), so the three quadrature nodes are \(b=0\) and \(b=\pm\sqrt{3}(0.5)=\pm0.8660\). Define \(q(b)=f(y=1\mid b)\exp\{-H(1\mid b)\}\). Each \(q(b)\) is the conditional-data density; its weighted average approximates the integral with respect to the random-intercept distribution.

At \(b=0\), the calculations from the previous example give \(q(0)=0.6766\). At \(b=0.8660\), the measurement density is \(0.7979\exp[-2(0.8660)^2]=0.1780\). The cumulative hazard is \(0.1\exp[0.5(1+0.8660)]=0.2542\), so the survival probability is \(\exp(-0.2542)=0.7755\), giving \(q(0.8660)=0.1780(0.7755)=0.1381\).

At \(b=-0.8660\), the measurement density is again \(0.1780\). The cumulative hazard is \(0.1\exp[0.5(1-0.8660)]=0.1069\), and the survival probability is \(\exp(-0.1069)=0.8986\). Thus \(q(-0.8660)=0.1780(0.8986)=0.1600\).

The quadrature estimate of the marginal contribution is therefore \((2/3)(0.6766)+(1/6)(0.1381)+(1/6)(0.1600)=0.5007\), to four significant figures. It is a density contribution, not a probability for the subject’s exact measurement value. Its numerical value depends on the measurement units and on the model parameters.

Worked Example: Changing Censoring to an Observed Event

Keep the same hypothetical measurement and parameter values, but now suppose the subject experiences the event at \(T=1\), so \(\delta=1\). The longitudinal density and cumulative hazard at each node are unchanged. The event contribution now includes \(h(1\mid b)\), where \(h(1\mid b)=0.1(1.5)(1^{0.5})\exp\{0.5(1+b)\}=0.15\exp\{0.5(1+b)\}\).

At \(b=0\), the hazard is \(0.15\exp(0.5)=0.2473\). Multiplying the measurement density, hazard, and survival probability gives \(q(0)=0.7979(0.2473)\exp(-0.1649)=0.1673\). At \(b=0.8660\), the hazard is \(0.15\exp(0.9330)=0.3813\), so \(q(0.8660)=0.1780(0.3813)\exp(-0.2542)=0.05265\). At \(b=-0.8660\), the hazard is \(0.15\exp(0.0670)=0.1604\), giving \(q(-0.8660)=0.1780(0.1604)\exp(-0.1069)=0.02566\).

The marginal event contribution is approximately \((2/3)(0.1673)+(1/6)(0.05265)+(1/6)(0.02566)=0.1246\). Compared with the censored case, the integrand includes the hazard at the observed event time. Omitting that factor would incorrectly use a survival contribution in place of an event-time density.

Combining Subjects and Checking the Calculation

Under subject independence, the individual marginal contributions multiply. For likelihood computation, it is usually safer to add their logarithms than to multiply many small numbers directly. The next example combines the two hypothetical subjects just calculated: one censored at time 1 and one with an event at time 1, each with the same single measurement and parameter values.

Worked Example: A Two-Subject Likelihood and Log-Likelihood

Using the three-node approximations, the censored subject contributes approximately \(L_1=0.5007\), and the event subject contributes approximately \(L_2=0.1246\). Therefore the two-subject likelihood is \(L\approx(0.5007)(0.1246)=0.06239\). The corresponding log-likelihood is \(\ell=\log(0.06239)\approx-2.774\).

The negative log-likelihood is not a sign of an error: densities can exceed or fall below one, and a likelihood has no requirement to be a probability for the entire parameter space. The value \(-2.774\) is the log of the approximate likelihood for these two records at the stated parameter values. It is not an estimate of a treatment effect or a measure of evidence by itself; comparisons require evaluating likelihoods under specified parameter values or models.

Several checks help reveal implementation errors. Confirm that \(\sigma^2>0\), \(\tau^2>0\), \(\kappa>0\), and \(\rho>0\), and that all times are on the same scale used to specify \(\kappa\). Verify that \(\delta_i=1\) includes the hazard multiplier and \(\delta_i=0\) does not. Check that the longitudinal residual uses the same \(b\) that enters the event hazard. For nonzero \(\beta_1\), check the numerical cumulative-hazard integral separately; a wrong integral affects both event and censoring contributions. When many subjects are used, compute log contributions and use stable numerical methods rather than multiplying extremely small values.

Alternatives, Limitations, and a Common Pitfall

Approach	What it computes	Main caution
Numerical quadrature of the marginal likelihood	A weighted approximation to the integral over the random intercept	Check convergence as the node count increases, especially if the integrand is concentrated or skewed.
Laplace approximation	An approximation based on the integrand near its mode	Accuracy can suffer when the integrand is not well approximated locally by a Gaussian shape.
Plug in a subject-specific intercept estimate	A conditional contribution evaluated at one estimated value of \(b_i\)	It is not the marginal likelihood and generally ignores uncertainty in \(b_i\).

The quadrature example is deliberately small enough to compute by hand. In actual fitting, three nodes may be inadequate; numerical resolution should be increased until likelihoods and resulting estimates are sufficiently stable for the intended precision. If the model has a time-varying latent mean, numerical evaluation of \(H_i(T_i\mid b)\) may be needed at every quadrature node and parameter value. The method also relies on the specified random-effect distribution and the conditional-independence and censoring assumptions; quadrature cannot repair a misspecified model.

Common mistake: “The likelihood for a subject is the longitudinal density times the survival contribution evaluated at the best-fitting random intercept.” This treats a latent quantity as if it were observed and known. The correction is to integrate the conditional-data contribution over the random-intercept density, as in the displayed formula. A subject-specific estimate can be useful in other summaries, but it is not a substitute for marginalization in the observed-data likelihood.

Key takeaway: For this shared-random-intercept model, each subject contributes an integral of the longitudinal density times the conditional event or censoring contribution, weighted by the Gaussian random-effect density. Numerical integration approximates that integral; subject contributions are then multiplied, or their logarithms added.

Check Your Understanding

Answer each question, then compare with the short solutions below.

For an observed event, which additional factor appears in the conditional time-to-event contribution compared with right censoring?
Why can the longitudinal and event likelihood components not generally be marginalized separately and then multiplied?
With \(\tau=0.4\), what are the three random-intercept nodes for the three-point quadrature rule?
For the model with \(\beta_1=0\), \(\gamma z=0\), \(\kappa=0.2\), \(\rho=2\), and \(\alpha(\beta_0+b)=0\), calculate the cumulative hazard at \(T=3\).
What is the main error in evaluating the conditional-data contribution at a single estimated \(b_i\) and calling it the marginal likelihood?

Solutions. (1) The observed-event contribution includes \(h_i(T_i\mid b_i)\); both event and censoring contributions include \(\exp\{-H_i(T_i\mid b_i)\}\). (2) Both data sources depend on the same random intercept, so separate marginalization would generally lose their dependence. (3) The nodes are \(0\) and \(\pm\sqrt{3}(0.4)=\pm0.6928\). (4) \(H(3)=0.2(3^2)\exp(0)=1.8\). (5) It ignores uncertainty in the random intercept; the marginal likelihood integrates the conditional contribution over its Gaussian distribution.