Why Average Over Random Effects?
A joint model links a patient’s longitudinal biomarker and event process through shared random effects. Even after observing a patient’s biomarker history, those effects are generally not known exactly. A dynamic survival prediction must therefore account for the remaining uncertainty about them. The key operation is to calculate survival for each possible random-effects value and average those survival probabilities using the patient-specific distribution of the effects given the observed history.
This average is not generally the same as calculating survival at the mean random effect. Survival is a nonlinear function of the effects, so substituting a single “typical” value can give a different prediction. The distinction matters especially when the patient’s random-effects distribution is diffuse or survival changes sharply across its range.
From Conditional Survival to an Individual Prediction
Suppose the event hazard conditional on \(b_i\) and baseline covariates \(x_i\) is \(h_i(u\mid b_i,x_i)\). Its conditional survival function from the start of follow-up through time \(u\) is
For a patient known to be event-free at \(s\), the conditional probability of surviving beyond \(t\), given the random effects, is the survival ratio \(S_i(t\mid b_i,x_i)/S_i(s\mid b_i,x_i)\). Equivalently, it is the exponential of the negative cumulative hazard from \(s\) to \(t\). The individual prediction averages that conditional probability over the distribution of random effects given the patient’s observed history.
This is an application of the law of total probability: conditional on the information available at \(s\), survival past \(t\) can occur under many possible values of the unobserved random effects. The integral adds the conditional survival probabilities, weighting each value by its conditional probability for this patient. The conditioning in \(q_i\) must match the prediction target. In particular, if the prediction is conditional on being event-free at \(s\), the distribution used in the integral must be appropriate to that condition. The formula does not itself prescribe how to update that distribution from the observed history.
Worked Predictions
Worked Example: Two Possible Random-Effects Values
Hypothetical data. At a landmark, suppose a patient’s history-conditioned distribution assigns probability \(0.60\) to \(b=-1\) and \(0.40\) to \(b=1\). Assume the future hazard is constant over the next two years conditional on \(b\): \(h(u\mid b)=0.08\exp(0.7b)\) per year. Calculate the patient’s probability of surviving those two years.
For a constant hazard, conditional survival over two years is \(\exp\{-2h(b)\}\). At \(b=-1\),
Thus the hazard rounds to \(0.03973\) per year and the survival probability to \(0.9236\). At \(b=1\),
Weight the two conditional probabilities by the patient-specific probabilities:
The predicted probability of surviving the next two years, conditional on the stated history and being event-free at the landmark, is about \(0.8440\). For comparison, the mean of the two random-effects values under these weights is \(-0.2\). Substituting that mean into the hazard gives \(\exp\{-2(0.08e^{-0.14})\}=0.8701\), which is not the integrated prediction. It answers a different, plug-in question.
Worked Example: Integrating With Posterior Quadrature Weights
Hypothetical data. Suppose a numerical approximation to the patient’s history-conditioned distribution uses nodes \(b=-1.2,0,0.4\), with weights \(0.25,0.50,0.25\). The weights sum to one. Assume \(h(u\mid b)=0.05e^{0.8b}\) per year over the next four years. Use these nodes to approximate survival through four years.
At each node, calculate the conditional hazard, multiply it by four years to obtain cumulative hazard, and exponentiate its negative. The results are:
| Node \(b\) | Weight | Hazard per year | Four-year survival |
|---|---|---|---|
| \(-1.2\) | 0.25 | \(0.05e^{-0.96}=0.01914\) | \(e^{-4(0.01914)}=0.9263\) |
| \(0\) | 0.50 | \(0.05000\) | \(e^{-4(0.05000)}=0.8187\) |
| \(0.4\) | 0.25 | \(0.05e^{0.32}=0.06886\) | \(e^{-4(0.06886)}=0.7593\) |
The weighted sum is
The quadrature approximation is therefore about \(0.8308\). The same nodes and weights imply a mean random effect of \(0.25(-1.2)+0.50(0)+0.25(0.4)=-0.2\). Plugging this mean into the hazard gives \(\exp\{-4(0.05e^{-0.16})\}\approx0.8433\), again not the weighted average. A quadrature rule approximates the integral by a weighted sum; its nodes and weights are chosen to represent the distribution and should not be mistaken for observed patient measurements.
Worked Example: Integrating Over an Intercept and Slope
Hypothetical data. Let the patient’s random-effects vector be \(b=(b_0,b_1)\). Suppose a quadrature approximation places equal weight \(0.25\) on four possible vectors: \((0,0)\), \((-0.5,0.4)\), \((0.5,-0.4)\), and \((0,0.6)\). Over each of two successive one-year intervals, assume the hazard is \(0.05e^{b_0}\) in the first interval and \(0.05e^{b_0+b_1}\) in the second. Calculate two-year survival by integrating over the vector.
For any node \((b_0,b_1)\), the two-year conditional survival is \(\exp\{-0.05e^{b_0}-0.05e^{b_0+b_1}\}\). For \((0,0)\), this is \(e^{-0.10}=0.9048\). For \((-0.5,0.4)\), it is \(e^{-0.05(e^{-0.5}+e^{-0.1})}=e^{-0.07557}=0.9272\). For \((0.5,-0.4)\), it is \(e^{-0.05(e^{0.5}+e^{0.1})}=e^{-0.1377}=0.8714\). For \((0,0.6)\), it is \(e^{-0.05(1+e^{0.6})}=e^{-0.1411}=0.8684\).
The individual prediction is the equally weighted average:
This example shows that the integration applies to the entire random-effects vector. When the hazard depends on both an intercept and a slope, the joint distribution of those effects matters; separate averages of the components do not generally reproduce the average of the survival function.
Computing and Checking the Integral
For a low-dimensional random-effects distribution, numerical quadrature can approximate the integral with nodes \(b_k\) and normalized weights \(w_k\):
For higher dimensions, simulation-based integration is often convenient. Draw \(b_i^{(m)}\) from the history-conditioned distribution and average the corresponding conditional survival probabilities. If the draws are independent, the Monte Carlo standard error is estimated by the sample standard deviation of those probabilities divided by the square root of the number of draws. With dependent simulation draws, an appropriate effective sample size or simulation-based error estimate is needed.
Several checks help detect a faulty calculation. The weights must be nonnegative and sum to one, and the final survival prediction must lie between zero and one. For a fixed landmark, predicted survival should not increase as the future time \(t\) increases. Increase the number of quadrature nodes or simulation draws and verify that the prediction is stable to a practically meaningful precision. Inspect whether the numerical method adequately represents regions of the random-effects distribution that contribute substantially to survival; a small number of central nodes can miss important tail behavior.
Common Mistakes, Alternatives, and Limits
- Plugging in the mean effect. This computes survival at one value, not the average over the patient’s uncertainty. Use it only when that plug-in quantity is the intended target or as a rough diagnostic comparison.
- Averaging hazards and then calculating survival. In general, \(\exp\{-\int_s^t E[h_i(u\mid b)\mid\mathcal{H}_i(s)]\,du\}\) is not equal to the expected conditional survival. The target here averages the survival probabilities themselves.
- Using the wrong conditioning distribution. The weights must correspond to the prediction information and event-free status specified in the target. Omitting relevant conditioning or applying it twice can distort the prediction.
- Treating a quadrature approximation as exact. Numerical integration error is distinct from model uncertainty. Check sensitivity to the integration scheme and its resolution.
Direct numerical integration, quadrature, and simulation are alternative ways to approximate the same expectation when they use the same conditional distribution and model. A plug-in prediction is computationally simpler but generally targets a different quantity. All these approaches remain limited by the fitted joint model: an incorrect hazard structure, an inadequate longitudinal model, or shared effects that fail to capture relevant dependence can produce misleading predictions even when the integral is evaluated accurately. A prediction based on fixed fitted parameters also omits uncertainty in those estimates.
Check Your Understanding
Answer each question, then compare your reasoning with the solutions.
- In one sentence, explain why conditional survival must be averaged over a history-conditioned random-effects distribution.
- A patient has two possible random-effects values, \(b=0\) and \(b=1\), with respective conditional probabilities \(0.7\) and \(0.3\). The corresponding survival probabilities to a horizon are \(0.90\) and \(0.60\). What is the integrated prediction?
- Why does survival calculated at the posterior mean random effect generally differ from the integrated prediction?
- What two basic checks should be made on the weights in a quadrature approximation?
- A numerical prediction increases when the future horizon is extended, with the landmark and conditioning history unchanged. What should be investigated?
Solutions. (1) The random effects are not known exactly, so the prediction averages conditional survival over the possible effects weighted by their probabilities given the patient’s information. (2) \(0.7(0.90)+0.3(0.60)=0.81\). (3) The survival function is nonlinear in the random effects, so evaluating it at their mean generally does not equal its expectation. (4) Each weight should be nonnegative, and the weights should sum to one. (5) Check the conditional survival calculations, the cumulative-hazard integration, and the numerical integration or simulation; survival for a fixed landmark should be nonincreasing with the horizon.