Diagnosing Event-Order Heterogeneity

Why Diagnose Event-Order Heterogeneity?

A model that gives every recurrence the same baseline pattern and the same covariate effect can be concise and precise. But it may hide a meaningful feature of the data: the first event may occur under a different pattern of risk from later events, or a covariate may be more strongly associated with later recurrences. Diagnosing these differences helps determine whether a shared-effect specification is a reasonable summary.

“Event-order heterogeneity” means that a recurrence-related quantity varies with the event number $k$. It is not the same as a covariate effect changing over follow-up time. The earlier material on proportional effects addresses whether a covariate effect changes with time; here the question is whether the pattern differs across event orders. If the apparent difference also changes over time, an event-order diagnostic alone cannot describe that pattern.

Definition: For subject $i$, let $k$ index event order, such as the first or second recurrence. An event-order-specific hazard model can be written $$ \lambda_{ik}(t\mid X_i)=\lambda_{0k}(t)\exp(X_i^\mathsf{T}\beta_k) $$ while the subject is in the risk set appropriate to the chosen model. Here $\lambda_{0k}(t)$ is the baseline hazard for event order $k$, and $\beta_k$ is its covariate-effect vector. Baseline heterogeneity means the $\lambda_{0k}(t)$ patterns differ; covariate-effect heterogeneity means the $\beta_k$ differ.

This notation describes a family of models, not a claim that every recurrent-event analysis has the same risk set or estimand. For example, ordered-event conditional models and marginal event-order models define who is at risk differently. Use the risk-set definition and interpretation belonging to the selected model; comparisons across event orders are meaningful only when the time scale and event-order estimand are clearly stated.

What to Compare and Under Which Assumptions

A useful diagnostic separates two questions. First, do the baseline recurrence patterns differ by order? Second, do the covariate effects differ by order? A model may have different baselines but a common covariate effect, or similar baselines but different covariate effects. Allowing event-order-specific baselines does not by itself require event-order-specific covariate coefficients.

For covariate effects, the shared-effect model imposes $\beta_1=\cdots=\beta_K$. The event-specific model relaxes this restriction. For a single covariate, the null hypothesis of a common effect is $H_0:\beta_1=\cdots=\beta_K$. For several covariates, the null states equality across orders for every coefficient being tested. These hypotheses concern log hazard ratios; exponentiating a coefficient gives a hazard ratio for a one-unit increase in that covariate, under the model’s other conditions.

Before using these comparisons, check that the analysis has a coherent event-order definition and time scale, that the risk sets match the selected estimand, and that censoring and entry are handled appropriately. The model must also represent the covariate form adequately. A proportional-effects assumption is relevant when interpreting a constant $\beta_k$ over time; evidence against it can make an event-order comparison misleading if the time pattern differs across orders. Use subject-level robust covariance estimates or a subject-level resampling method when repeated events create within-subject dependence. Robust variance addresses uncertainty from clustering; it does not correct a wrong risk set, time scale, or mean-model form.

Baseline curves are often compared visually using event-order-specific estimated cumulative baseline hazards, $\widehat{\Lambda}_{0k}(t)$. A cumulative hazard summarizes hazard over time; its slope reflects the local hazard pattern. Curves that separate, converge, or cross suggest different patterns, but do not alone establish statistical evidence. Estimates at later event orders can be unstable when few subjects reach those orders. Show numbers at risk or otherwise disclose the support behind each curve.

There is not a universal likelihood-ratio test for equality of arbitrary, unspecified baseline hazard functions. A practical inferential comparison can be a prespecified contrast of estimated cumulative hazards at a fixed time, using its estimated covariance, or a joint resampling procedure for the curves. A single-time contrast does not test equality over the whole follow-up period. If many times or orders are inspected, pointwise intervals or p-values do not automatically control the overall error rate.

Testing Covariate-Effect Differences

For fitted event-order-specific coefficients, a Wald test compares the estimated contrasts with their estimated covariance. For two event orders and one covariate, define $\widehat{\Delta}=\hat{\beta}_1-\hat{\beta}_2$. Its standard error uses both variances and their covariance:

SE(\widehat{\Delta})= \sqrt{\widehat{\operatorname{Var}}(\hat{\beta}_1) +\widehat{\operatorname{Var}}(\hat{\beta}_2) -2\widehat{\operatorname{Cov}}(\hat{\beta}_1,\hat{\beta}_2)}.

The covariance matters because estimates from different event orders may be based on some of the same subjects. For more than two orders or covariates, collect the independent contrasts in a vector $\hat{\boldsymbol{\Delta}}$, with covariance matrix $\widehat{V}$. The joint Wald statistic is $\hat{\boldsymbol{\Delta}}^\mathsf{T}\widehat{V}^{-1}\hat{\boldsymbol{\Delta}}$; under the null and standard large-sample conditions, it is compared with a chi-square distribution having degrees of freedom equal to the number of independent contrasts.

Justification: A contrast is a linear combination of estimates. For $D=\hat{\beta}_1-\hat{\beta}_2$, variance algebra gives $\operatorname{Var}(D)=\operatorname{Var}(\hat{\beta}_1)+\operatorname{Var}(\hat{\beta}_2)-2\operatorname{Cov}(\hat{\beta}_1,\hat{\beta}_2)$. If the joint estimates are approximately normal and the null states $D=0$, then $D/SE(D)$ is approximately standard normal. This yields a two-sided Wald test or, equivalently, a confidence interval for the contrast.

Worked Example: Testing a Two-Order Covariate Contrast

Suppose a hypothetical ordered-recurrence analysis estimates the log hazard ratio for a one-unit increase in a covariate as $\hat{\beta}_1=0.20$ for the first event and $\hat{\beta}_2=0.32$ for the second. The robust standard errors are 0.05 and 0.04, and the estimated covariance is 0.001. Test whether the covariate effects differ, using a two-sided normal-reference Wald test.

The estimated contrast is $\widehat{\Delta}=0.20-0.32=-0.12$; its magnitude is 0.12. Its variance is $0.05^2+0.04^2-2(0.001)=0.0025+0.0016-0.002=0.0021$. Thus $SE(\widehat{\Delta})=\sqrt{0.0021}=0.0458$, and the standardized statistic is $-0.12/0.0458=-2.62$.

The two-sided normal-reference p-value is approximately 0.0088. At the 0.05 level, these data provide evidence that the covariate association differs between the two event orders under the fitted model. The estimated hazard ratios are $\exp(0.20)=1.22$ and $\exp(0.32)=1.38$; their ratio is $\exp(-0.12)=0.887$ for first-order versus second-order effects. This is evidence about model-based associations, not evidence that the covariate causes recurrence.

Worked Diagnostics

Worked Example: Comparing Baseline Patterns at a Fixed Time

In a hypothetical analysis on a common follow-up-time scale, estimated baseline cumulative hazards at 12 months are 0.30 for the first event and 0.48 for the second. Their standard errors are 0.05 and 0.06, and their estimated covariance is 0.0012. Consider the prespecified contrast $\widehat{\Lambda}_{02}(12)-\widehat{\Lambda}_{01}(12)$.

The estimated difference is $0.48-0.30=0.18$. Its variance is $0.06^2+0.05^2-2(0.0012)=0.0036+0.0025-0.0024=0.0037$, giving $SE=\sqrt{0.0037}=0.0608$. The standardized contrast is $0.18/0.0608=2.96$, with a two-sided normal-reference p-value of about 0.0031.

At this specified horizon, the estimated baseline cumulative hazard for the second event is higher by 0.18 on the cumulative-hazard scale, with evidence against equality at 12 months. This is not a test that the baseline hazards are unequal at every time, nor does it mean the second event is 0.18 times as likely. The result is interpretable only with the stated model, event-order risk sets, and time scale.

Worked Example: A Joint Test Across Three Event Orders

Suppose a hypothetical model estimates one covariate’s log hazard ratios as $(0.10, 0.25, 0.40)$ for event orders 1, 2, and 3. For the contrasts $\hat{\beta}_2-\hat{\beta}_1$ and $\hat{\beta}_3-\hat{\beta}_1$, suppose the estimated covariance matrix is $\begin{pmatrix}0.03&0.01\\0.01&0.03\end{pmatrix}$. Test the joint null that the three effects are equal.

The contrast vector is $\hat{\boldsymbol{\Delta}}=(0.15,0.30)^\mathsf{T}$. The inverse covariance matrix is $\begin{pmatrix}37.5&-12.5\\-12.5&37.5\end{pmatrix}$. Therefore the Wald statistic is $W=37.5(0.15^2)-25(0.15)(0.30)+37.5(0.30^2)=3.094$. There are two independent contrasts, so the reference distribution is chi-square with 2 degrees of freedom. The p-value is about 0.213.

At the 0.05 level, this test does not provide sufficient evidence that the covariate effects differ across these three orders. It does not establish that the effects are equal: the estimates remain uncertain, and the test’s power depends on the information available, especially for the third event.

Reading Plots and Choosing a Response

A useful display places the event-order-specific cumulative baseline hazard curves on the same axes and time scale, with uncertainty bands where available. Inspect whether separation is persistent or limited to a region with sparse risk sets, whether curves cross, and whether apparent differences are large relative to their uncertainty. Alongside the plot, report the number of subjects still contributing to each event-order risk set at relevant times. The curves can be difficult to compare if later-order estimates are based on few subjects or if event-order definitions imply different populations.

For covariate effects, plot the event-order-specific coefficient estimates or hazard ratios with confidence intervals, and include a joint test of equality when that is the scientific question. Overlapping separate confidence intervals are not a formal test of equality: the uncertainty of a difference depends on covariance. Conversely, a small p-value does not describe the size or practical importance of the difference. Give the contrast and its interval or scale-specific estimates.

If order-specific effects appear different, a reasonable model-specification response is to report the evidence and fit or present an event-specific-effects model as a sensitivity analysis. If baseline patterns differ but covariate effects appear similar, a shared covariate effect with event-order-specific baselines may still be appropriate. The earlier material on event-specific covariate effects provides the basic interpretation of separate coefficients. Choosing between these specifications also depends on the estimand, scientific purpose, and support in the data—not just a significance threshold.

The main alternatives answer different questions. A pooled shared-effect model is parsimonious but can obscure order differences. Event-order-specific coefficients allow those associations to differ but use more parameters and may be imprecise for late events. Stratifying baselines by order allows different baseline shapes while retaining a common coefficient, but baseline curves should not be treated as comparable without a common time scale and clear risk-set interpretation. A formal comparison of baseline curves at selected times is narrower than a whole-curve comparison. Model-based tests and plots complement one another; neither substitutes for checking the design and assumptions.

Common Mistakes and Limitations

Dropping the covariance term. Treating two order-specific estimates as independent can give the wrong standard error. Use the joint robust covariance estimate or an appropriate subject-level resampling method.
Equating a visual difference with a test result. Curves that look separated may be uncertain, and a pointwise test cannot establish a global difference across all times. State exactly what contrast or global procedure was used.
Calling a nonsignificant result proof of sameness. Failure to reject equality can reflect limited information, particularly at later orders. Report the estimate and uncertainty, not just the p-value.
Confounding event order with follow-up time. Event order and calendar or analysis time are distinct. A difference across event orders is not automatically a time-varying covariate effect; changing time patterns require their own assessment.
Ignoring changing risk sets. Later-event estimates often rely on fewer eligible subjects. A contrast can be unstable or target a different selected population if risk-set definitions are not kept clear.
Interpreting a baseline cumulative-hazard difference as a probability difference. Cumulative hazard is not itself an event probability. Keep the scale and event-order estimand explicit.

These diagnostics are limited by model choice, sparse later events, and multiplicity when many covariates, orders, or time points are examined. A robust variance estimate protects against some within-subject dependence in inference, but it cannot repair bias from informative censoring, an unsuitable time scale, incorrect covariate specification, or an estimand that does not match the question.

Check Your Understanding

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

In an event-order model, what is the difference between baseline heterogeneity and covariate-effect heterogeneity?
Two order-specific coefficient estimates have standard errors 0.03 and 0.04, with covariance 0.0006. What is the standard error of their difference?
A test of equality across three event-order coefficients gives a p-value of 0.18. What is a valid conclusion at significance level 0.05?
Why should a cumulative baseline-hazard plot show information about the risk sets contributing to later event orders?

Solutions. (1) Baseline heterogeneity means $\lambda_{0k}(t)$ differs across orders; covariate-effect heterogeneity means $\beta_k$ differs. (2) The variance is $0.03^2+0.04^2-2(0.0006)=0.0009+0.0016-0.0012=0.0013$, so the standard error is $\sqrt{0.0013}=0.0361$. (3) At 0.05, the test does not provide sufficient evidence against equal coefficients; it does not prove equality. (4) Later orders may have few eligible subjects, making estimates unstable and apparent curve differences less reliable.