Tutorials › Joint Modeling of Longitudinal and Time-to-Event Data › Worked Linear-Marker Weibull Prediction

Dynamic Prediction · Tutorial 543 of 697

Worked Linear-Marker Weibull Prediction

Learn to calculate a patient’s horizon-specific conditional event risk by integrating Weibull survival over the random-effects distribution informed by the patient’s marker history.

Advanced 11 min read

What You'll Learn

  • Specify a linear mixed marker model linked to a Weibull event hazard.
  • Distinguish risk conditional on survival to a landmark from marginal survival.
  • Integrate conditional survival over a patient-specific random-effects distribution.
  • Use Simpson’s rule to approximate the cumulative hazard integrals in a worked prediction.
  • Update random-effects weights for survival to the landmark.
  • Diagnose plug-in, conditioning, and numerical integration errors.

From a Marker History to a Conditional Risk

A dynamic prediction at landmark time \(s\) asks about a patient who is known to be event-free at \(s\), and whose longitudinal marker history has been observed up to that time. For a future horizon \(w>0\), the target is the probability of an event during \((s,s+w]\), conditional on that information. This is not generally the same as evaluating a hazard at the most recent observed marker value: the prediction averages future survival over uncertainty about the patient’s latent marker trajectory.

This tutorial works through that calculation for a linear mixed model with a random intercept and a Weibull baseline hazard. We treat the fitted model parameters as known to focus on integration over the patient’s random effect. In an applied analysis, uncertainty in estimated parameters may also be propagated, for example by the posterior-draw approach covered in the earlier material on Bayesian predictions.

Model and assumptions: For patient \(i\), let \(Y_{ij}\) be the marker observed at time \(t_{ij}\). The longitudinal model is \(Y_{ij}=\beta_0+\beta_1t_{ij}+b_i+\epsilon_{ij}\), where \(b_i\sim N(0,\sigma_b^2)\) and the measurement errors \(\epsilon_{ij}\sim N(0,\sigma^2)\) are independent of one another and of \(b_i\). Conditional on \(b_i\), the event hazard is \(h_i(t\mid b_i)=\lambda\nu t^{\nu-1}\exp\{\gamma^\mathsf{T}x_i+\alpha m_i(t;b_i)\}\), where \(m_i(t;b_i)=\beta_0+\beta_1t+b_i\). Here \(\lambda>0\), \(\nu>0\), and \(m_i(t;b_i)\) is the latent marker mean. We assume the specified marker and hazard models are adequate, that event risk depends on the marker through this shared random effect, and that censoring is conditionally non-informative for the prediction calculation.

The baseline hazard is Weibull: \(h_0(t)=\lambda\nu t^{\nu-1}\), so the baseline cumulative hazard is \(H_0(t)=\lambda t^\nu\). Conditional on \(b_i\), the hazard is known as a function of time and the event-free probability from time 0 to \(u\) is \(\exp\{-H_i(u\mid b_i)\}\), where \(H_i(u\mid b_i)=\int_0^u h_i(t\mid b_i)\,dt\). The random effect makes this cumulative hazard patient-specific.

Prediction target: Let \(\mathcal{Y}_i(s)\) denote patient \(i\)’s marker observations through \(s\). The conditional event risk through horizon \(s+w\) is \(P(s<T_i\le s+w\mid T_i>s,\mathcal{Y}_i(s),x_i)\). If \(p_i(b\mid T_i>s,\mathcal{Y}_i(s),x_i)\) is the random-effects density conditional on the patient’s marker history and event-free status at \(s\), then this risk equals $$ 1-\int \exp\{-[H_i(s+w\mid b)-H_i(s\mid b)]\}\,p_i(b\mid T_i>s,\mathcal{Y}_i(s),x_i)\,db. $$

The density in this expression must reflect both the observed marker history and survival to the landmark. Up to a normalizing constant, it is proportional to the marker likelihood given \(b\), multiplied by the random-effects prior density and the event-free probability through \(s\). Thus, survival to \(s\) is part of the information used to form the prediction, not an optional adjustment made after calculating risk.

Why the Formula Averages Survival

For a fixed value \(b\), conditional survival from \(s\) to \(s+w\) is \(\exp\{-[H_i(s+w\mid b)-H_i(s\mid b)]\}\). The patient’s random effect is not observed, so conditional survival given the patient’s information is the average of this quantity over the conditional distribution of \(b\). The event risk is one minus that average.

1
Condition on a random effect.
For each possible \(b\), the hazard model gives the cumulative hazard increment from \(s\) to \(s+w\), and hence conditional survival over that interval.
2
Average over patient-specific uncertainty.
Weight the conditional survival probabilities by \(p_i(b\mid T_i>s,\mathcal{Y}_i(s),x_i)\) and integrate over \(b\).
3
Convert survival to risk.
Subtract the integrated conditional survival from 1. The result is the event probability in the stated future interval, conditional on the information at the landmark.

This order matters: in general, averaging the hazards or random effects first and then exponentiating does not give the same result as averaging conditional survival. Likewise, using only the marker likelihood to form the random-effects distribution omits information from the fact that the patient has remained event-free through \(s\).

Worked Example: Integrating Risk for a Hypothetical Patient

Suppose a hypothetical fitted model has \(\beta_0=0\), \(\beta_1=0.10\) marker units per time unit, \(\alpha=0.50\) per marker unit, \(\gamma^\mathsf{T}x_i=0\), \(\lambda=0.12\), and \(\nu=1.5\). The hazard is therefore \(0.18\sqrt{t}\exp\{0.50(0.10t+b)\}\). Consider a prediction from landmark \(s=1\) to horizon \(s+w=2\). For this example, suppose numerical integration of the patient’s marker and survival information has produced two random-effects quadrature nodes and weights: \(b=-0.20\) with weight \(0.40\), and \(b=0.20\) with weight \(0.60\). These weights sum to 1 and represent the conditional random-effects distribution at the landmark.

First compute the common time integral. The factor involving time is \(f(t)=0.18\sqrt{t}\exp(0.05t)\), and the random effect contributes the separate multiplier \(\exp(0.50b)\). Composite Simpson integration with step \(h=0.25\) gives \(\int_0^1f(t)\,dt\approx0.1218\) and \(\int_0^2f(t)\,dt\approx0.3587\). For example, through time 1 the Simpson calculation is \(0.25/3\) times \([f(0)+f(1)+4\{f(0.25)+f(0.75)\}+2f(0.50)]\), which is approximately \(0.1218\). Applying Simpson’s rule directly over the prediction interval gives an integral of approximately \(0.2369\).

At \(b=-0.20\), the random-effects multiplier is \(\exp(-0.10)=0.9048\). The cumulative hazard increment is \(0.2369(0.9048)=0.2143\), giving conditional survival \(\exp(-0.2143)=0.8071\). At \(b=0.20\), the multiplier is \(\exp(0.10)=1.1052\), the increment is \(0.2369(1.1052)=0.2618\), and conditional survival is \(\exp(-0.2618)=0.7697\).

Average these survival probabilities using the supplied weights: \(0.40(0.8071)+0.60(0.7697)=0.7847\). The predicted risk is \(1-0.7847=0.2153\), or about 21.53%. In context, under this fitted model and conditional random-effects distribution, the hypothetical patient’s probability of an event between times 1 and 2, given survival to time 1 and the marker history, is approximately 21.5%.

Landmark Survival Changes the Random-Effects Weights

A useful way to see why landmark conditioning matters is to start with weights based only on marker data and then update them for survival. Bayes’ rule gives \(p_i(b\mid T_i>s,\mathcal{Y}_i(s),x_i)\propto P(T_i>s\mid b,x_i)\,p_i(b\mid\mathcal{Y}_i(s),x_i)\). The survival factor favors random effects associated with longer survival, relative to the marker-only distribution. The amount of change depends on the model and the landmark.

Worked Example: Updating Two Weights for Survival to the Landmark

Use the same hypothetical hazard parameters and landmark \(s=1\). Suppose the marker history alone gives equal quadrature weights, \(0.50\) at \(b=-0.20\) and \(0.50\) at \(b=0.20\). The cumulative hazard through time 1 is the integral \(0.1218\) multiplied by \(\exp(0.50b)\).

For \(b=-0.20\), this cumulative hazard is \(0.1218(0.9048)=0.1102\), so survival to time 1 is \(\exp(-0.1102)=0.8957\). For \(b=0.20\), it is \(0.1218(1.1052)=0.1346\), and survival is \(\exp(-0.1346)=0.8741\). The unnormalized weights after incorporating survival are \(0.50(0.8957)=0.4479\) and \(0.50(0.8741)=0.4371\). Their sum is \(0.8850\), so the normalized landmark weights are approximately \(0.5061\) and \(0.4939\).

Using the conditional survivals from time 1 to 2 calculated in the preceding example gives integrated survival \(0.5061(0.8071)+0.4939(0.7697)=0.7886\). The conditional event risk is therefore \(1-0.7886=0.2114\), or 21.14%. This differs slightly from the risk formed using the equal marker-only weights, because the patient’s survival to the landmark shifts the random-effects distribution.

Why a Plug-In Random Effect Can Mislead

A tempting shortcut is to replace the patient’s random effect by its conditional mean and calculate a single risk. That plug-in prediction is \(1-\exp\{-[H_i(s+w\mid \bar b)-H_i(s\mid\bar b)]\}\), where \(\bar b\) is the conditional mean. It is not generally equal to the integral of conditional risks. The exponential transformation is nonlinear, so variation in \(b\) can affect the integrated risk even when the conditional mean is known exactly.

Worked Example: Comparing Integration With a Mean-Effect Plug-In

Return to the two-node distribution in the first worked example: weights \(0.40\) and \(0.60\) at \(b=-0.20\) and \(b=0.20\). Its conditional mean is \(0.40(-0.20)+0.60(0.20)=0.040\). A plug-in calculation at this mean uses multiplier \(\exp(0.50(0.040))=\exp(0.020)=1.0202\), so the interval cumulative hazard is \(0.2369(1.0202)=0.2417\). The resulting plug-in risk is \(1-\exp(-0.2417)=0.2147\).

The integrated risk from the first example was approximately 0.2153, about 0.0007 higher than the plug-in risk of approximately 0.2147, or 0.07 percentage points. The gap is modest here, but need not be modest for a more variable or skewed conditional random-effects distribution, a stronger association parameter, or a longer horizon. Neither method is a universal bound for the other; the risk should be computed by integrating conditional survival when the target averages over random-effect uncertainty.

Numerical Checks, Alternatives, and Limitations

For a one-dimensional random intercept, Gaussian quadrature or numerical integration over the conditional random-effects density is often practical. With multiple random effects, deterministic integration may become expensive; Monte Carlo integration and Laplace approximation are standard alternatives. Their general computational principles are covered in the earlier material on Monte Carlo integration for risk predictions and Laplace approximation in dynamic prediction. Whatever method is used, the target remains the same: average conditional survival under the random-effects distribution given both marker history and survival to the landmark.

Useful diagnostics include checking that quadrature weights are nonnegative and sum to 1, that the integrated survival lies between 0 and 1, and that the resulting risk also lies between 0 and 1. Repeat the calculation with more quadrature nodes or a finer integration grid: material changes indicate inadequate numerical accuracy. Check time units carefully, since the Weibull hazard and the marker slope must use compatible units. A prediction that changes substantially after a small numerical refinement, or produces negative risk, signals an implementation or convergence problem.

The calculation also has substantive limitations. It depends on the specified longitudinal and event models, the shared-random-effect association, and the assumed random-effects distribution. Treating fitted parameters as fixed omits parameter-estimation uncertainty. Conditional risk is not a causal effect, and a predicted risk is not a guarantee about an individual outcome. In addition, this particular worked model uses a random intercept; richer random-effects structures require integration over a higher-dimensional distribution.

  • Do not omit landmark survival. The random-effects density should condition on \(T_i>s\), not only on the observed marker values.
  • Do not average hazards and then exponentiate. Average the conditional survival probabilities, then subtract from 1.
  • Do not report the wrong interval. The prediction uses the cumulative hazard increment \(H_i(s+w\mid b)-H_i(s\mid b)\), not cumulative hazard from time 0 alone.
  • Do not confuse a plug-in estimate with integration. Evaluating risk at the conditional mean random effect generally differs from averaging risk over its conditional distribution.
Key takeaway: For a linear-marker Weibull prediction, compute the conditional cumulative hazard increment at each random-effect value, convert it to conditional survival, average survival under the patient-specific distribution conditioned on marker history and survival to the landmark, and subtract the result from 1.

Check Your Understanding

Answer each question using the definitions and calculations above. Short solutions follow the questions.

  1. Why does the random-effects distribution for a prediction at time \(s\) condition on survival to \(s\), as well as the marker history?
  2. If integrated conditional survival through the horizon is \(0.73\), what is the conditional event risk, and how should it be interpreted?
  3. For a fixed random effect, why is the cumulative hazard increment over \((s,s+w]\) a difference of cumulative hazards rather than the cumulative hazard through \(s+w\)?
  4. What numerical checks would you use to assess a quadrature-based prediction?

Solutions. (1) Being event-free at \(s\) is observed information that can change the relative plausibility of random effects; Bayes’ rule weights the marker-based distribution by \(P(T_i>s\mid b,x_i)\). (2) The risk is \(1-0.73=0.27\), or 27%; under the fitted model and conditioning information, this is the probability of an event during the specified future interval. (3) Conditional survival from \(s\) onward is determined by accumulated hazard after \(s\), which is \(H_i(s+w\mid b)-H_i(s\mid b)\); using \(H_i(s+w\mid b)\) would include hazard accumulated before the landmark. (4) Check that weights are nonnegative and sum to 1, integrated survival and risk are in \([0,1]\), and the estimate is stable when the number of nodes or grid resolution is increased.

Clinical Biostats