Piecewise Exponential Hazards for MCMC

Why Divide Follow-Up Time into Intervals?

In a Bayesian joint model, the survival likelihood depends on the hazard over each subject's observed follow-up. A piecewise exponential baseline represents that hazard with a separate constant rate on each of several time intervals. This gives a flexible stepwise baseline and, under a common prior specification, makes the baseline-hazard update especially simple within MCMC.

The simplification does not remove the joint-model structure. Conditional on a subject's random effects and the other model parameters, the subject-specific relative hazard can still vary with time through the longitudinal trajectory. The survival likelihood must then account for that variation when calculating exposure. The key is to separate the interval-specific baseline rate from the time-varying relative hazard.

Definition: Fix interval boundaries \(0=a_0<a_1<\cdots<a_K\), and define interval \(k\) as \([a_{k-1},a_k)\). A piecewise-constant baseline hazard has \(h_0(t)=\lambda_k\) for \(t\) in interval \(k\), where \(\lambda_k>0\). Each \(\lambda_k\) is a hazard rate, with units of inverse time.

For subject \(i\), let \(T_i\) be the event time, \(C_i\) the censoring time, \(\tau_i=\min(T_i,C_i)\) the observed follow-up time, and \(\delta_i=I(T_i\le C_i)\) the event indicator. Conditional on the longitudinal random effects \(b_i\) and the other parameters, write the hazard as \(h_i(t\mid b_i)=h_0(t)q_i(t)\), where \(q_i(t)=\exp\{\eta_i(t)\}\) is the positive relative hazard. For example, \(\eta_i(t)\) may include fixed covariates and an association term involving the subject's fitted longitudinal trajectory at time \(t\). The derivation below does not require \(q_i(t)\) to be constant.

Event and Exposure Contributions

For right-censored survival data, an observed event contributes the hazard at the observed event time; every subject contributes survival through their observed follow-up. The survival contribution is the exponential of minus the cumulative hazard. Thus the likelihood contribution for subject \(i\), conditional on the model parameters and \(b_i\), is

L_i=\{h_i(\tau_i\mid b_i)\}^{\delta_i} \exp\left\{-\int_0^{\tau_i}h_i(u\mid b_i)\,du\right\}.

This expression treats a censored observation as survival through \(\tau_i=C_i\), without an event factor. Define the event count in interval \(k\) and the subject's relative-hazard-weighted exposure in that interval by

D_k=\sum_i \delta_i I\{a_{k-1}\le T_i<a_k\}, \qquad A_{ik}=\int_{a_{k-1}}^{\min(\tau_i,a_k)} I\{\tau_i\ge a_{k-1}\}q_i(u)\,du,

with \(A_{ik}=0\) if \(\tau_i<a_{k-1}\). For an event at a boundary, use the stated half-open convention: an event at \(a_k\) is assigned to the interval beginning at \(a_k\), provided that interval is in the modeled follow-up. A single boundary point has zero duration and therefore does not change the exposure integral.

The total weighted exposure is \(E_k=\sum_i A_{ik}\). Substituting \(h_i(t\mid b_i)=\lambda_kq_i(t)\) within interval \(k\) gives the survival likelihood, as a function of the baseline rates,

L(\boldsymbol{\lambda}\mid\text{other parameters, data}) \ \propto\ \prod_{k=1}^K \lambda_k^{D_k}\exp(-\lambda_kE_k).

The event factors \(q_i(T_i)\) are omitted from this proportionality because they do not depend on \(\boldsymbol{\lambda}\). The exposure \(E_k\) is generally not just total person-time in interval \(k\): it weights time at risk by the relative hazard \(q_i(t)\). When \(q_i(t)=1\), it reduces to ordinary person-time.

Worked Example: Building Event Counts and Weighted Exposure

Suppose follow-up is divided into \([0,2)\) and \([2,5)\). Consider three hypothetical subjects, conditional on their other parameters. Subject 1 has an event at time 1.5 and \(q_1(t)=2\). Subject 2 is censored at time 3 and has \(q_2(t)=0.5\). Subject 3 has an event at time 4 and has \(q_3(t)=1.5\). These relative hazards are constant in time for this example.

Subject 1 contributes \(2(1.5)=3\) units of weighted exposure to the first interval, and one event there. Subject 2 contributes \(0.5(2)=1\) to the first interval and \(0.5(3-2)=0.5\) to the second; censoring at time 3 creates no event. Subject 3 contributes \(1.5(2)=3\) to the first interval and \(1.5(4-2)=3\) to the second, with one event in the second interval.

Therefore \(D_1=1\), \(D_2=1\), \(E_1=3+1+3=7\), and \(E_2=0.5+3=3.5\). The baseline-rate part of the likelihood is

\lambda_1^1\exp(-7\lambda_1)\lambda_2^1\exp(-3.5\lambda_2).

The event relative-hazard factors are \(q_1(1.5)=2\) and \(q_3(4)=1.5\); they multiply this expression but do not change its dependence on \(\lambda_1,\lambda_2\). The censored subject contributes exposure despite having no event factor.

Gamma Full Conditional Updates

Suppose the interval-specific baseline hazards have independent Gamma priors parameterized by shape and rate:

\lambda_k\sim\operatorname{Gamma}(r_k,s_k), \qquad p(\lambda_k)\propto \lambda_k^{r_k-1}\exp(-s_k\lambda_k), \quad \lambda_k>0.

Here \(r_k\) is the prior shape and \(s_k\) is the prior rate, not the scale. Conditional on all other parameters—including the random effects and the quantities determining \(q_i(t)\)—the likelihood kernel for \(\lambda_k\) is \(\lambda_k^{D_k}\exp(-E_k\lambda_k)\). Multiplying prior and likelihood kernels gives the full conditional:

p(\lambda_k\mid\text{rest, data}) \propto \lambda_k^{r_k+D_k-1}\exp\{-(s_k+E_k)\lambda_k\}, \qquad \lambda_k\mid\text{rest, data} \sim \operatorname{Gamma}(r_k+D_k,\ s_k+E_k).

This is a Gibbs update: draw directly from the displayed conditional distribution rather than using a Metropolis proposal for \(\lambda_k\). The update is conditional on the current values of all parameters that determine \(q_i(t)\); it is not a marginal posterior distribution for \(\lambda_k\). The joint-model longitudinal likelihood does not add terms to this conditional unless it depends directly on \(\lambda_k\).

Worked Example: A Conditional Gamma Update

Use the hypothetical data from the previous example. Suppose the prior for each interval hazard is \(\operatorname{Gamma}(2,1)\), with shape 2 and rate 1. For interval 1, \(D_1=1\) and \(E_1=7\). Therefore

\lambda_1\mid\text{rest, data}\sim \operatorname{Gamma}(2+1,1+7)=\operatorname{Gamma}(3,8).

For interval 2, \(D_2=1\) and \(E_2=3.5\), so

\lambda_2\mid\text{rest, data}\sim \operatorname{Gamma}(2+1,1+3.5)=\operatorname{Gamma}(3,4.5).

For a Gamma distribution parameterized by shape and rate, the mean is shape divided by rate. The conditional means are \(3/8=0.375\) and \(3/4.5=0.667\) time\(^{-1}\), respectively. These are means of the conditional distributions at the current values of the other parameters, not summaries that account for their posterior uncertainty.

When the Relative Hazard Changes Within an Interval

The piecewise exponential assumption applies to the baseline hazard, not necessarily to each subject's full hazard. If \(q_i(t)\) varies over time, the exposure integral must retain that variation. For instance, if \(q_i(t)=\exp\{\alpha m_i(t)\}\), where \(m_i(t)\) is the longitudinal trajectory, then \(A_{ik}\) is the integral of that exponential over the subject's at-risk portion of interval \(k\). Replacing it by interval duration times a single arbitrary value generally changes the likelihood.

When the integral has no convenient closed form, it can be evaluated by numerical quadrature conditional on the current trajectory and parameters. The Gamma full conditional remains valid if the resulting \(E_k\) is calculated accurately and does not depend on \(\lambda_k\). If the relative-hazard predictor itself depends on a baseline hazard parameter, this conditional derivation no longer applies as written.

Worked Example: Integrating a Time-Varying Relative Hazard

Suppose a hypothetical subject is at risk from time 0 to 2, with one interval covering that span, and has \(q_i(t)=\exp(0.2t)\). There is no event before the end of follow-up. The weighted exposure is not simply \(2q_i(0)\); it is

A_{i1}=\int_0^2\exp(0.2t)\,dt =\frac{\exp(0.4)-1}{0.2} \approx\frac{1.4918-1}{0.2} =2.459.

The contribution to the hazard-rate likelihood is therefore \(\exp(-2.459\lambda_1)\). Treating \(q_i(t)\) as its value at time 0 would instead use exposure \(2\), while using its value at time 2 would give \(2\exp(0.4)\approx2.984\). Neither endpoint substitution equals the required integral. This check matters because an exposure error directly changes the posterior rate parameter \(s_1+E_1\).

Interval Boundaries, Diagnostics, and Limitations

Boundary choices affect the model, not merely the bookkeeping. They determine which event count informs each \(\lambda_k\), how weighted exposure is partitioned, and where the baseline hazard is allowed to jump. The interval structure should be specified and recorded clearly; changing cutpoints during analysis changes the likelihood and prior-to-parameter correspondence.

Worked Example: Moving a Boundary Changes the Likelihood

A hypothetical subject has an event at time 2.8 and a constant relative hazard \(q_i(t)=1.5\). Compare boundaries \([0,2,4]\) with \([0,3,4]\). With the first choice, the event lies in interval 2. The subject contributes \(1.5(2)=3\) units of weighted exposure to interval 1 and \(1.5(2.8-2)=1.2\) to interval 2. The baseline-rate likelihood contribution is

1.5\lambda_2\exp(-3\lambda_1-1.2\lambda_2).

With the boundary moved to time 3, the event lies in interval 1, and the subject contributes \(1.5(2.8)=4.2\) units of exposure there and none to interval 2. The contribution becomes

1.5\lambda_1\exp(-4.2\lambda_1).

The data and event time are unchanged, but the model assigns the event and risk exposure to different baseline parameters. This illustrates why boundary placement can affect posterior estimates, particularly with few events.

Useful checks include verifying that each event is assigned to exactly one interval, that censored subjects contribute exposure through their censoring times, and that exposure is zero after follow-up ends. For each interval, check that \(D_k\) is a nonnegative integer and \(E_k\ge0\). Recompute selected exposure integrals independently, especially when \(q_i(t)\) varies. A negative exposure, an event assigned to multiple intervals, or exposure continuing after an event or censoring time indicates an implementation error.

A structurally correct calculation can still fit poorly. Very wide intervals impose a nearly constant baseline over long spans; very narrow intervals can leave intervals with no events and make estimates sensitive to their priors. Compare posterior baseline-hazard summaries across reasonable, prespecified cutpoint choices, and examine whether important conclusions change. This is a sensitivity analysis, not a way to choose boundaries solely because they give a preferred result.

The piecewise exponential model is often compared with a Weibull baseline, which imposes a smooth parametric shape, and spline baselines, which allow smooth variation subject to modeling and prior choices. The piecewise model is transparent and has the convenient conditional update derived here, but its hazard can jump at cutpoints. The alternatives avoid fixed step changes in different ways, but do not automatically provide the same intervalwise Gamma update. None removes the need to check whether the baseline representation is adequate for the data.

Common pitfall: Treating every subject's interval exposure as raw time at risk, even when \(q_i(t)\) varies. The correct exposure is the integral of the relative hazard over time at risk. Raw person-time is valid only when \(q_i(t)=1\) throughout that exposure, or when the model's relative hazard is otherwise accounted for in the calculation.

Practical MCMC Workflow

A reliable implementation separates event allocation, exposure calculation, and the hazard update. At each iteration, use the current longitudinal and association parameters to evaluate \(q_i(t)\), compute each \(A_{ik}\), aggregate \(E_k\) and \(D_k\), and draw each \(\lambda_k\) from its full conditional when the prior assumptions permit. Then use the updated hazards in the other parameter updates required by the joint model.

Fix the time partition.
Record ordered cutpoints and the boundary convention; confirm that all modeled event times fall in the covered intervals.

Allocate events and calculate exposure.
Compute \(D_k\) from observed events and \(A_{ik}\) by integrating the current \(q_i(t)\) only over time at risk.

Update the baseline rates.
With independent Gamma(shape, rate) priors, draw \(\lambda_k\) from \(\operatorname{Gamma}(r_k+D_k,s_k+E_k)\).

Validate the implementation.
Check interval membership, exposure totals, units, and selected likelihood calculations against hand-worked cases.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a right-censored subject, which survival-likelihood contributions remain, and which event contribution is absent?
An interval has \(D_k=2\), \(E_k=5\), and prior \(\operatorname{Gamma}(3,2)\) in shape-rate form. What is the full conditional distribution of \(\lambda_k\)?
If \(q_i(t)=1\) throughout an interval, what does its weighted exposure reduce to?
A subject's event occurs exactly at an interior cutpoint \(a_k\). Under the half-open interval convention used here, which interval receives the event?
Why can moving a cutpoint change the posterior even when the observed follow-up times stay fixed?

Solutions. (1) The subject contributes survival through the censoring time, represented by the cumulative-hazard exponential, but no event hazard factor. (2) Multiplication of the likelihood kernel by the prior gives \(\operatorname{Gamma}(3+2,2+5)=\operatorname{Gamma}(5,7)\). (3) The integral of \(q_i(t)=1\) over time at risk is ordinary time at risk in that interval. (4) The event is assigned to the interval beginning at \(a_k\), namely \([a_k,a_{k+1})\), if that interval is modeled. (5) The cutpoint changes which baseline rate receives the event and how exposure is divided across rates; consequently it changes the likelihood as a function of the interval-specific hazards.