Worked Joint Model for Recurrent Hospitalizations and Death

Why Model the Three Processes Together?

A disease marker measured repeatedly may describe changing health, while hospitalizations recur and death prevents any later measurement or hospitalization. An analysis that treats these processes as unrelated can miss dependence: for example, an unobserved tendency toward a higher marker may also be associated with more hospitalizations and greater mortality. A joint model represents that dependence explicitly, while preserving the distinction between recurrent hospitalizations and the mutually exclusive causes of terminal death.

The model below is a worked specification, not a claim that one specification suits every application. Its effects are conditional on a person’s random effects and measured covariates. This conditioning matters: conditional hazard ratios are not generally equal to population-average hazard ratios after the random effects have been integrated out.

A Joint Model for the Marker and Event Processes

Definition: For person \(i\), let \(Y_{ij}\) be the disease-marker measurement at time \(t_{ij}\), \(m_i(t)\) the latent marker trajectory, \(N_i(t)\) the cumulative number of hospitalizations by time \(t\), and \(T_i\) the time to death or censoring. If death occurs, \(J_i\) denotes its cause, one of \(1,\ldots,K\). Let \(b_i=(b_{0i},b_{1i})^\mathsf{T}\) be person-specific random effects. A joint model specifies the marker distribution and the recurrent-event and cause-specific death processes together, with shared random effects and/or the latent marker linking them.

One simple longitudinal submodel is a linear mixed model:

Y_{ij}=m_i(t_{ij})+\epsilon_{ij}, \qquad m_i(t)=\beta_0+\beta_1t+b_{0i}+b_{1i}t.

Here \(\epsilon_{ij}\) are measurement errors, often assumed independent with \(\epsilon_{ij}\sim N(0,\sigma^2)\), and \(b_i\sim N(0,\Sigma)\). The event processes use the latent trajectory \(m_i(t)\), not the noisy observed value \(Y_{ij}\). This avoids treating measurement error as if it were true short-term biological change.

For a recurrent hospitalization process, a possible intensity is

\lambda_{Ri}(t\mid b_i)=Y_{R0}(t) \exp\{\alpha_R m_i(t)+\kappa_R b_{0i}\}, \qquad t<T_i.

For death from cause \(k\), use a cause-specific hazard

\lambda_{Dki}(t\mid b_i)=\lambda_{Dk0}(t) \exp\{\gamma_k^\mathsf{T}X_i+\alpha_{Dk}m_i(t)+\kappa_{Dk}b_{0i}\}, \qquad k=1,\ldots,K.

The \(Y_{R0}(t)\) and \(\lambda_{Dk0}(t)\) terms are baseline rates; \(X_i\) is a vector of measured covariates. The coefficients \(\alpha_R\) and \(\alpha_{Dk}\) link the latent marker to hospitalization and cause-\(k\) death, while the \(\kappa\) coefficients link the shared random intercept to those processes. Other specifications can share both random effects or use a different association structure. An intensity is a conditional instantaneous event rate; a cause-specific death hazard is the instantaneous rate of death from cause \(k\) among those still alive just before time \(t\).

Assumptions for this specification: The marker model and event intensities are correctly specified; the random-effects distribution and measurement-error assumptions are adequate; and, conditional on covariates and \(b_i\), the marker measurements, recurrent-event process, and terminal-event process have the stated dependence structure. Death ends the at-risk period for hospitalization and for every death cause. Censoring is assumed independent of the event processes conditional on the modeled information. If censoring depends on unmodeled health or event history, this likelihood does not account for that dependence.

Conditional Likelihood and Its Interpretation

Conditional on \(b_i\), the recurrent-event contribution up to the terminal time \(T_i\) has a survival factor for having no hospitalization between events and one intensity factor for each observed hospitalization. If \(t_{i1},\ldots,t_{in_i}\) are the hospitalization times, this contribution is

\left\{\prod_{\ell=1}^{n_i}\lambda_{Ri}(t_{i\ell}\mid b_i)\right\} \exp\left\{-\int_0^{T_i}\lambda_{Ri}(u\mid b_i)\,du\right\}.

If death occurs from cause \(J_i\), the terminal-process contribution is the cause-specific hazard at death multiplied by survival from all death causes up to that time:

\lambda_{D J_i i}(T_i\mid b_i) \exp\left\{-\int_0^{T_i}\sum_{k=1}^{K}\lambda_{Dki}(u\mid b_i)\,du\right\}.

For censoring rather than observed death, the terminal contribution is the exponential survival term alone. The complete conditional likelihood also includes the longitudinal measurement density. The subject-specific likelihood is obtained by integrating the product of these conditional contributions over the distribution of \(b_i\). Thus, shared random effects induce dependence even when the processes are conditionally independent given those effects.

Interpretation: Holding \(b_i\), covariates, and time fixed, increasing the latent marker by one unit multiplies the hospitalization intensity by \(\exp(\alpha_R)\) and the cause-\(k\) death hazard by \(\exp(\alpha_{Dk})\). These are conditional associations, not necessarily causal effects. A positive coefficient says that a higher latent marker is associated with a higher conditional rate or hazard.

Worked Example: Evaluating Marker-Linked Rates

Suppose, hypothetically, that the fitted marker trajectory for one person is \(m_i(t)=1.2+0.15t+0.2-0.05t\), with time in years and marker units on a specified clinical scale. At year 2, \(m_i(2)=1.2+0.30+0.20-0.10=1.60\). The person’s random intercept is \(b_{0i}=0.20\). For illustration, take constant baseline rates \(Y_{R0}=0.008\) hospitalizations per year, \(\lambda_{D10}=0.004\) per year, and \(\lambda_{D20}=0.002\) per year. Let \((\alpha_R,\kappa_R)=(0.40,0.50)\), \((\alpha_{D1},\kappa_{D1})=(0.50,0.25)\), and \((\alpha_{D2},\kappa_{D2})=(0.20,-0.25)\), with no additional covariate terms.

The conditional hospitalization intensity at year 2 is

0.008\exp\{0.40(1.60)+0.50(0.20)\} =0.008\exp(0.74) \approx 0.01677\ \text{hospitalizations per year}.

For cause 1, the conditional death hazard is \(0.004\exp\{0.50(1.60)+0.25(0.20)\}=0.004\exp(0.85)\approx0.009359\) per year. For cause 2, it is \(0.002\exp\{0.20(1.60)-0.25(0.20)\}=0.002\exp(0.27)\approx0.002620\) per year. These rates apply at that time conditional on the specified trajectory and random effect. The assumptions used in this calculation are the stated model forms, constant baseline rates, and the supplied fitted values; it does not establish that these assumptions are adequate for actual data.

Worked Conditional Event Likelihood

The likelihood calculation makes explicit how recurrent events and competing death causes both contribute. Under constant rates, the hospitalization and terminal survival factors are exponentials of integrated rates. Each observed event contributes its corresponding instantaneous rate at its event time.

Worked Example: One Hospitalization Before Cause-Specific Death

Consider hypothetical conditional rates held constant for two years: hospitalization intensity \(0.02195\) per year, cause-1 death hazard \(0.06376\) per year, and cause-2 death hazard \(0.01000\) per year. Suppose one hospitalization occurs at year \(0.8\), followed by observed cause-1 death at year 2. Assume the processes are conditionally independent given the random effects and that these rates are constant over the interval.

The event-rate factors are the hospitalization intensity at the hospitalization time and the cause-1 hazard at death. Their product is

(0.02195)(0.06376)=0.001399532\approx0.001400\ \text{year}^{-2}.

The integrated hospitalization intensity is \(2(0.02195)=0.04390\). The integrated total death hazard, summing over both competing causes, is \(2(0.06376+0.01000)=0.14752\). The conditional event-time likelihood contribution is therefore

(0.02195)(0.06376) \exp\{-0.04390-0.14752\} =0.001399532\exp(-0.19142) \approx0.001156\ \text{year}^{-2}.

The exponent includes both the absence of additional hospitalizations and survival from all causes of death until year 2. The competing cause-2 hazard contributes to survival even though cause 1 is observed. The units are year\(^{-2}\) because the contribution contains two event-time density factors, each with units year\(^{-1}\). This is a conditional likelihood contribution, not a probability by itself; event times are continuous.

Shared Effects and Conditional Hazard Ratios

The marker and shared random effect can both contribute to event risk. A coefficient’s interpretation must therefore hold the other terms in the linear predictor fixed. For example, changing \(b_{0i}\) affects the event intensities through the \(\kappa\) coefficients, while a change in \(m_i(t)\) acts through the \(\alpha\) coefficients.

Worked Example: Comparing Two Conditional Effects

Suppose the hospitalization marker coefficient is \(\alpha_R=0.40\), and the cause-1 death marker coefficient is \(\alpha_{D1}=0.50\), as in the preceding example. For two otherwise identical people whose latent marker values differ by \(0.5\) units at a given time, with the same random effects, the hospitalization intensity ratio is

\exp\{0.40(0.5)\}=\exp(0.20)\approx1.221.

Thus, the person with the higher latent marker has a hospitalization intensity about 22.1% higher at that time, conditional on the model and other terms. For cause-1 death, the corresponding ratio is \(\exp\{0.50(0.5)\}=\exp(0.25)\approx1.284\), or about 28.4% higher conditional hazard. These comparisons are not statements about lifetime risks or causal effects.

Alternatively, compare two people with the same marker and covariates but random intercepts differing by \(0.30\). The hospitalization intensity ratio from the random-effect term is \(\exp\{0.50(0.30)\}=\exp(0.15)\approx1.162\); the cause-1 hazard ratio is \(\exp\{0.25(0.30)\}=\exp(0.075)\approx1.078\). Shared effects thus permit a latent propensity to be associated with multiple processes, with the strength and direction specified separately for each process.

Estimation, Diagnostics, and Alternatives

In practice, the random effects are unobserved. Estimation therefore combines the marker density and event-process likelihood and integrates over the random-effects distribution, commonly using numerical quadrature or simulation-based methods. The examples above condition on specified random effects to make the event contributions transparent; they do not replace this integration in population-level estimation.

Several checks are important before interpreting fitted parameters:

Assess the marker trajectory. Residual patterns, nonlinear trends, and changing variability can indicate that a linear trajectory or constant error variance is inadequate. Sparse measurements can make an individual trajectory uncertain, which in turn affects event-risk estimates.
Check event-process fit. Compare observed and model-predicted recurrent-event counts and timing, and assess whether the chosen baseline intensity captures changes over time or by event history. Check cause-specific event patterns separately; a good fit for total mortality can conceal poor fit for one cause.
Examine dependence assumptions. The chosen shared random effects and marker links are substantive assumptions. Sensitivity to alternative link terms or random-effects distributions can reveal that conclusions depend strongly on the assumed dependence structure.
Check censoring and event recording. Informative loss to follow-up, missed hospitalizations, or uncertain cause of death can invalidate the assumed likelihood unless represented in the model.

A pair of separate models may be simpler, but it does not generally account for dependence induced by shared latent health or informative termination of recurrence at death. A model that uses the marker as an observed time-varying covariate is another alternative, but it does not automatically correct for measurement error or account for the marker’s latent trajectory. The joint model is useful when this shared structure is scientifically plausible and supported by the data. It is not automatically superior: a complex model may be weakly identified when measurements or events are sparse, and its conclusions can be sensitive to baseline-rate and random-effects assumptions.

Common Pitfall

Common mistake: Interpreting \(\exp(\alpha_{Dk})\) as the population-average death-risk ratio, or as a causal effect of changing the marker. Correction: It is a conditional cause-specific hazard ratio for a one-unit difference in the latent marker, holding the modeled random effects and covariates fixed. Marginal comparisons require accounting for the distribution of random effects and the fact that people with different histories remain at risk for different lengths of time. The model alone does not establish causation.

A second frequent error is to use only the hazard for the observed death cause in the survival exponent. Survival up to death requires survival from every competing cause, so the exponent contains the integral of \(\sum_k\lambda_{Dki}(t\mid b_i)\). The observed cause enters through the event-rate factor at \(T_i\), not by removing the other causes from the preceding survival contribution.

Key takeaway: A joint model links a marker trajectory, recurrent hospitalization intensity, and cause-specific death hazards through explicitly stated association terms and shared random effects. Correct interpretation depends on distinguishing conditional effects from marginal ones and including all relevant event-free survival contributions in the likelihood.

Check Your Understanding

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

If \(\alpha_R=0.30\), what is the conditional hospitalization intensity ratio for a one-unit higher latent marker, holding other terms fixed?
Why does the terminal survival contribution include every cause-specific death hazard when cause 1 is observed?
With a constant hospitalization intensity of \(0.02\) per year over 3 years and one hospitalization, what is the hospitalization-process contribution if the intensity at the event time is also \(0.02\) per year?
What does a positive shared-random-intercept coefficient in the cause-2 hazard imply, conditional on the model?

Solutions. (1) The ratio is \(\exp(0.30)\approx1.350\): about a 35.0% higher conditional intensity. (2) The person must remain alive from all causes up to the observed death time; the cause observed at death determines the event factor, while all causes contribute to preceding survival. (3) The contribution is \(0.02\exp\{-3(0.02)\}=0.02\exp(-0.06)\approx0.01884\) per year, a one-event time-density contribution. (4) Among otherwise comparable individuals at the same modeled marker and covariates, a higher random intercept multiplies the conditional cause-2 hazard according to the exponential of its coefficient times the intercept difference; it does not by itself establish a causal effect or a population-average ratio.