Retinal Thickness and Diabetic Retinopathy Progression

Why Link Retinal Thickness to Progression?

Suppose a hypothetical study follows people with diabetic retinopathy, recording optical imaging measurements of retinal thickness at several visits and the time to first sight-threatening progression. The scans describe a changing biological process, while progression is an event that may occur between visits. A participant who progresses may stop contributing routine scans, so the available measurements can become fewer precisely for people whose disease is advancing. An analysis that treats the observed scans and event times as unrelated may therefore miss important dependence.

A joint model addresses this by representing each participant’s underlying thickness trajectory and relating that latent trajectory to the event hazard. The model does not assert that retinal thickness causes progression. It estimates an association conditional on specified covariates and assumptions. Nor does it presume that thicker retina always means more severe disease: thickness can reflect different pathological processes, so the clinical direction must be justified for the particular measurement and endpoint.

Model and Assumptions

Let \(Y_{ij}\) be participant \(i\)'s measured retinal thickness, in micrometers, at time \(t_{ij}\), measured in years from baseline. Let \(T_i\) be time to first sight-threatening diabetic retinopathy progression, and let \(C_i\) be the censoring time. The observed follow-up time is \(\tilde T_i=\min(T_i,C_i)\), and \(\delta_i=1\) if progression is observed and \(0\) if follow-up is censored. Define \(m_i(t)\) as the participant's latent, or underlying, mean thickness at time \(t\).

One simple longitudinal submodel is a linear trajectory with participant-specific intercept and slope:

Here \(\beta_0\) is the population intercept, \(\beta_1\) is the population change in thickness per year, and \(b_{0i}\) and \(b_{1i}\) are participant-specific deviations. A standard specification assumes \(\boldsymbol b_i=(b_{0i},b_{1i})^\mathsf{T}\sim N(\boldsymbol 0,D)\) and independent measurement errors \(\epsilon_{ij}\sim N(0,\sigma^2)\), conditional on the random effects and covariates. Covariates can be added to either submodel as appropriate.

For illustration, use a proportional hazards event submodel whose association is with the current latent thickness, centered at 250 micrometers and scaled per 10 micrometers:

The baseline hazard is \(h_0(t)\), \(\boldsymbol Z_i\) contains event-model covariates, and \(\alpha\) is the log hazard ratio for a 10-micrometer higher current latent thickness, holding the included covariates fixed. Centering changes the reference level, not the association estimate. The form assumes proportional hazards conditional on the covariates and random effects. It also assumes the chosen current-value association is an adequate description; a slope association or a different trajectory may be needed if scientific reasoning supports it.

Key assumptions include a correctly specified longitudinal mean and random-effects distribution; measurement errors that are conditionally independent with the specified variance structure; and censoring that is noninformative conditional on the modeled information. The visit schedule must also be ignorable for the likelihood used: for example, visits must not occur more often because of unmodeled worsening thickness. If visits or censoring remain informative after conditioning, they may need their own model. The event definition and time origin must be consistent across participants.

Likelihood: How the Two Processes Connect

Conditional on \(\boldsymbol b_i\), the longitudinal observations contribute their measurement densities. The event contribution is the hazard at the observed time if progression occurred, multiplied by the probability of remaining event-free up to that time. Write the conditional survival function as \(S_i(t\mid\boldsymbol b_i)=\exp\{-\int_0^t h_i(u\mid\boldsymbol b_i)\,du\}\). The participant's likelihood contribution is then integrated over the random-effects distribution:

This expression makes the connection explicit: conditional on a participant's random effects, the two data components are modeled together; the same random effects enter both components; and the likelihood averages over their unknown values. For a censored participant, \(\delta_i=0\), so the event factor is survival through the censoring time, not a density for an event at that time.

The factorization follows by conditional independence: given \(\boldsymbol b_i\) and covariates, the joint density of the observed scans and event record is the scan-density product times the event-record contribution. Multiplying by the density of \(\boldsymbol b_i\) and integrating gives the marginal likelihood because the random effects are unobserved. If the random effects affect both thickness and the hazard, integrating them jointly preserves the dependence between the two observed processes. Fitting the two components separately generally does not.

Worked Likelihood Contributions

Diagnostics, Alternatives, and Limitations

A joint model depends on both submodels being credible. Inspect longitudinal residuals against fitted thickness and time for curvature, changing spread, or systematic patterns. Compare observed and model-predicted trajectories, and examine whether the assumed random-effects distribution and measurement-error structure are plausible. For the event component, assess whether the proportional-hazards assumption is reasonable and compare predicted and observed event patterns over follow-up. Check sensitivity to influential participants, sparse scans, and alternative trajectory forms. A model that fits event times but poorly describes scan trajectories, or vice versa, is not adequate simply because it converges.

Also ask whether the visit schedule could be informative. If participants are scanned more frequently when clinicians suspect worsening disease, a model treating visit times as fixed and ignorable may not capture all relevant dependence. Similarly, a single current-value association can be inadequate if the rate of thickness change, a delayed effect, or another aspect of the trajectory is more relevant. These alternatives require scientific justification and careful model checking rather than selection solely by statistical fit.

Joint modeling is not a cure for poor measurement or an ill-defined endpoint. Retinal thickness may vary by retinal region, instrument, and image quality; those sources of variation require appropriate measurement protocols and model terms. Sight-threatening progression should be defined in advance, with a clear time origin and consistent event ascertainment. If death or another event prevents progression from being observed, treating it as ordinary censoring requires additional assumptions; a competing-risks analysis may be more appropriate. Finally, model-based associations and predictions should not be described as causal effects without a causal design and assumptions that support that interpretation.