Why Allow Effects to Differ by Recurrence Order?
A covariate can have a different association with the first recurrence than with later recurrences. For example, a preventive treatment might be associated with a lower first-event hazard but a weaker association among people who have already had an event. A model that forces one coefficient to apply at every order cannot represent that pattern. An event-specific covariate effect assigns a distinct coefficient to the covariate at each recurrence order, while retaining the risk-set and estimand structure of the chosen recurrent-event model.
The recurrence-order distinction concerns the covariate effect, not necessarily the baseline hazard. A model can already allow each order its own baseline hazard and still constrain the covariate coefficient to be shared. Conversely, separate coefficients let the covariate association vary by order. These are different modeling choices.
Modeling the Order-Specific Effects
Write \(\lambda_{ik}(t)\) for the event hazard for subject \(i\) at recurrence order \(k\), under the risk-set definition of the selected model. For a PWP model, the risk set for order \(k\) includes subjects who have experienced the required earlier recurrences; its hazard is conditional on reaching that order. A WLW formulation instead gives event-order-specific marginal hazards. The coefficients must be interpreted according to that underlying model, not as if these two risk sets were interchangeable.
A general event-specific proportional-hazards specification is
Here \(\mathcal{H}_i(t)\) represents the history needed to define the model and risk set, and \(\lambda_{0k}(t)\) is the baseline hazard for order \(k\). Additional covariates can be included in the linear predictor. The event-specific effect is \(\beta_k\); the order-specific baseline hazard \(\lambda_{0k}(t)\) is a separate feature. This formulation assumes that, within each order, the covariate’s hazard ratio is constant over the model’s time scale. It does not assume that the hazard ratio is the same across orders.
For \(K\) orders, one convenient parameterization uses order 1 as the reference:
The interaction \(x_i I(k=j)\) allows the covariate slope at order \(j\) to differ from the order-1 slope. Thus \(\beta_1\) is the effect at order 1, and \(\delta_j\) is the change in the log hazard ratio from order 1 to order \(j\). Alternatively, estimate a separate \(\beta_k\) for each order. These parameterizations describe the same type of model when all required parameters are estimable.
Worked Examples
Worked Example: Interpreting an Event-Order Interaction
Hypothetical data and fitted model. Suppose 240 participants, 120 assigned to each of two treatment groups, are followed for recurrent events. The analysis records 84 first events, 57 second events, and 31 third events. A model with separate baseline hazards by order and robust standard errors estimates the treatment log hazard ratio at order 1 as \(\hat\beta_1=0.20\), with order-2 interaction \(\hat\delta_2=0.45\) and order-3 interaction \(\hat\delta_3=0.10\). Treatment is coded \(x_i=1\) and control \(x_i=0\).
The estimated treatment coefficient at order 1 is \(0.20\), so the estimated hazard ratio is
At order 2, the coefficient is \(0.20+0.45=0.65\), giving \(\exp(0.65)=1.916\). At order 3, it is \(0.20+0.10=0.30\), giving \(\exp(0.30)=1.350\). Under the fitted model, the estimated treatment-to-control hazard ratios are therefore 1.221, 1.916, and 1.350 at orders 1, 2, and 3, respectively. These estimates describe conditional hazards under this model; they do not establish that treatment caused the differences.
The order-2 interaction is itself a log hazard-ratio contrast. Exponentiating it gives \(\exp(0.45)=1.568\): the order-2 treatment hazard ratio is estimated to be 1.568 times the order-1 treatment hazard ratio. This comparison is distinct from either order-specific treatment-versus-control comparison.
Worked Example: Testing Whether Two Order-Specific Effects Differ
Hypothetical fitted results. For a binary exposure in a two-order model, suppose \(\hat\beta_1=0.30\) and \(\hat\beta_2=-0.10\), with robust standard errors 0.18 and 0.20. Suppose the estimated correlation between the two coefficient estimates is 0.25. Consider the null hypothesis \(H_0:\beta_1-\beta_2=0\) against a two-sided alternative.
The estimated contrast is \(0.30-(-0.10)=0.40\). Its variance must account for the covariance between the estimates. The covariance is \(0.25(0.18)(0.20)=0.009\), so
Thus \(SE(\hat\beta_1-\hat\beta_2)=\sqrt{0.0544}=0.2332\), and the Wald statistic is \(z=0.40/0.2332=1.715\). Its two-sided large-sample p-value is approximately 0.086. At the 0.05 level, this is not sufficient evidence to reject equal effects; it is not proof that the effects are equal.
A 95% Wald confidence interval for the log hazard-ratio contrast is \(0.40\pm1.96(0.2332)=(-0.0571,0.8571)\). Exponentiating the endpoints gives an interval of approximately \((0.945,2.36)\) for the ratio of the order-1 to order-2 hazard ratios. The interval includes 1, consistent with the two-sided test. The covariance term matters: treating the estimates as independent would give a different standard error and test.
Worked Example: A Precise First-Order Estimate and an Uncertain Later-Order Estimate
Hypothetical fitted results. In a study with 180 participants, a binary exposure is present in 90 participants. A model estimates a log hazard ratio of \(0.25\) (standard error \(0.15\)) for the first event and \(0.80\) (standard error \(0.50\)) for the third event. The third-order risk set contains far fewer observations than the first-order risk set. Assume the model converged and the third-order coefficient is estimable.
For the first event, the estimated hazard ratio is \(\exp(0.25)=1.284\). The approximate 95% log-scale interval is \(0.25\pm1.96(0.15)=(-0.044,0.544)\), which becomes \((0.957,1.723)\) on the hazard-ratio scale.
For the third event, the hazard ratio is \(\exp(0.80)=2.226\). The log-scale interval is \(0.80\pm1.96(0.50)=(-0.180,1.780)\); exponentiating gives approximately \((0.835,5.93)\). The point estimate is larger at order 3, but its interval is wide and includes 1. This is evidence of substantial uncertainty, not a precise finding that the exposure has a strong third-order effect. Fewer later-order events commonly reduce precision.
Estimation, Assumptions, and Comparisons
Fit the event-specific model using the estimation procedure appropriate to the selected recurrent-event framework. In a Cox-type model, the event-order indicators and covariate-by-order interactions enter the linear predictor, while order-specific baseline hazards can be handled by stratification. For a marginal model, retain its appropriate estimating equations and variance method; when robust inference is used, calculate contrasts from the robust covariance matrix. The interaction changes the covariate effect, not the definition of the risk set.
Interpretation relies on several conditions. The event order and risk sets must be coded according to the target model; the covariate and outcome records must be aligned with the subject and event order; and censoring must be adequately handled under the model’s assumptions. The proportional-effects assumption applies within each order for the specification above. The covariate is treated as time-fixed here. If later values of a covariate change during follow-up, that requires a different formulation. In addition, the model must have enough information at each order to estimate its coefficients; a nominally fitted coefficient can be unstable when events or exposed subjects are sparse.
A shared-effects model constrains \(\beta_1=\cdots=\beta_K\), which can be more precise when a common effect is plausible, but cannot describe order-specific associations. An event-specific model relaxes that constraint and answers a more detailed question, usually at the cost of more parameters and wider uncertainty. Fitting completely separate models by order can also produce order-specific estimates, but a pooled interaction model makes the contrasts explicit and can conveniently retain a common set of adjustment covariates. The appropriate choice depends on the estimand and scientific question, not only on which model gives narrower intervals.
To compare effects, test contrasts such as \(H_0:\beta_k-\beta_l=0\), or jointly test equality across all orders. For a Wald test, the variance of each contrast comes from the fitted covariance matrix; a multi-contrast test uses the full covariance matrix of the contrast vector. The test concerns equality of modeled log hazard ratios. It does not test whether baseline hazards are equal, nor does it establish a difference in cumulative event burden.
Specify whether the effects are conditional or marginal, what counts as each event order, and which risk set applies.
Use separate order-specific slopes or a reference-order slope plus covariate-by-order interactions.
Confirm that event-order coding, risk sets, convergence, and event counts support the requested coefficients.
Give order-specific hazard ratios with uncertainty, and calculate comparisons using the covariance between estimates.
Common Pitfalls and Practical Limits
- Assuming different baseline hazards imply different covariate effects. Stratifying on event order allows baseline hazards to vary; it does not, by itself, allow slopes to vary. Add an interaction or estimate separate slopes to do that.
- Comparing point estimates by eye. Two estimates that look different need not differ statistically. Test their contrast using its standard error, including covariance.
- Reading an interaction as an order-specific hazard ratio. In a reference-order parameterization, \(\delta_k\) is a difference in log effects. The order-\(k\) coefficient is \(\beta_1+\delta_k\).
- Ignoring the selected model’s estimand. A conditional PWP hazard and a marginal WLW hazard have different risk-set interpretations. Event-specific coefficients do not make those interpretations identical.
- Over-interpreting unstable later-order estimates. Sparse events can produce large standard errors, wide intervals, and sensitivity to a small number of observations. More coefficients do not automatically mean more reliable understanding.
Event-specific effects matter when the scientific question concerns how an association changes across recurrence order. They also make the precision-versus-flexibility trade-off visible: a shared coefficient may be efficient if its constraint is reasonable, while order-specific coefficients reveal variation only to the extent the data support it. Report the model, risk-set interpretation, order-specific estimates, and uncertainty together.
Check Your Understanding
Answer each question, then compare with the solutions.
- In a reference-order model, \(\hat\beta_1=-0.20\) and \(\hat\delta_2=0.30\). What are the estimated hazard ratios at orders 1 and 2, and what does \(\hat\delta_2\) represent?
- Why is the standard error of \(\hat\beta_1-\hat\beta_2\) generally not obtained by simply adding the two standard errors?
- A model has a separate baseline hazard at each event order but one shared covariate coefficient. Does it estimate event-specific covariate effects? Explain.
- What is one likely consequence of having very few observed third events when estimating a third-order coefficient?
Solutions. (1) The order-1 hazard ratio is \(\exp(-0.20)=0.819\). The order-2 coefficient is \(-0.20+0.30=0.10\), so its hazard ratio is \(\exp(0.10)=1.105\). The interaction is the order-2 minus order-1 log hazard-ratio contrast; its exponent, \(\exp(0.30)=1.350\), is the ratio of those hazard ratios. (2) The estimates may be correlated, so the contrast variance is \(\mathrm{Var}(\hat\beta_1)+\mathrm{Var}(\hat\beta_2)-2\mathrm{Cov}(\hat\beta_1,\hat\beta_2)\). (3) No. The baseline hazards vary, but the shared coefficient still constrains the covariate effect to be equal across orders. (4) The coefficient may be imprecise, with a large standard error and a wide confidence interval; its point estimate should not be treated as conclusive evidence of a large effect.