Checking Fine-Gray Fits with Residuals and Time Interactions

Why Check Time Constancy?

A Fine-Gray model with a constant covariate coefficient assumes that the covariate’s association with the subdistribution hazard is the same throughout follow-up. That assumption can be useful, but it is not guaranteed by fitting the model. If the association changes over time, a single coefficient compresses that change into one average effect and can give a misleading summary.

This tutorial uses the residuals and time-interaction features of cmprsk::crr to investigate that issue. Residual patterns are diagnostic evidence, not a definitive test. A time-by-covariate term gives a direct model-based check, but its conclusion depends on the chosen time function and on the rest of the model being reasonably specified.

Definition: Let $T$ be event or censoring time and $J$ the event type. For cause 1, the cumulative incidence function is $F_1(t)=P(T\leq t,J=1)$, and the cause-1 subdistribution hazard is $\lambda_1^*(t)=f_1(t)/\{1-F_1(t-)\}$, where $f_1(t)$ is the cause-1 event density. A proportional subdistribution-hazards model specifies $\lambda_1^*(t\mid x)=\lambda_{10}^*(t)\exp(\beta x)$. Its proportionality assumption is that $\beta$ does not vary with time.

The Fine-Gray risk-set construction accounts for competing events in estimating this model. The coefficient describes a multiplicative association with the subdistribution hazard; it is not, in general, a risk ratio or a direct ratio of cumulative incidence probabilities. The model also requires appropriately coded event types, independent subjects, and censoring that is independent of the event process under the analysis assumptions (or otherwise adequately handled). There must be enough observed cause-1 events across follow-up to estimate the effects. A residual plot cannot repair incorrect event coding, dependent censoring, or inadequate data support.

Residuals as a Diagnostic

A fitted crr object contains a residual matrix, available as fit$res, and the corresponding unique cause-1 event times, available as fit$uftime. The residual rows correspond to those event times; the columns correspond to model covariates. These are estimating-equation residual contributions, not ordinary outcome-minus-prediction residuals. Consequently, a point’s sign should not be read as a simple statement that the subject’s observed risk was above or below a predicted probability.

For a constant-coefficient model, inspect each covariate’s residuals against event time. A roughly patternless cloud is compatible with proportionality, although it does not prove it. A persistent trend or other systematic temporal structure can suggest that a coefficient changes over time. The plot can be noisy, especially with few cause-1 events, ties, or sparse late follow-up. A smooth line is a visual aid; it is not a calibrated hypothesis test.

Worked Example: Plotting a Treatment Residual

Suppose a hypothetical study follows 300 people for relapse (cause 1), death before relapse (cause 2), or censoring (code 0). There are 90 cause-1 events, 75 cause-2 events, and 135 censored observations. Treatment is coded 1 for the new treatment and 0 for control; age is represented in decades from its mean. Assume independent participants, valid event coding, and censoring plausibly independent of the event process after the covariates are considered. These assumptions make the model and diagnostic interpretable; they do not establish that the proportionality assumption holds.

Fit the constant-effect model and plot the treatment residuals against their event times:

\texttt{X <- cbind(treatment = treatment, age10 = age10)}\\ \texttt{fit0 <- crr(ftime, fstatus, cov1 = X, failcode = 1, cencode = 0)}\\ \texttt{plot(fit0\$uftime, fit0\$res[, "treatment"],}\\ \texttt{\ \ xlab = "Time (days)", ylab = "Treatment residual")}\\ \texttt{abline(h = 0)}\\ \texttt{lines(lowess(fit0\$uftime, fit0\$res[, "treatment"]), col = 2)}

Here each plotted point uses one unique cause-1 event time and its treatment-column residual. Suppose the hypothetical plot has residuals fluctuating around zero early in follow-up, followed by a smooth curve that tends to rise late in follow-up. The trend is a reason to investigate time variation, not proof of it: random fluctuation, an omitted covariate, sparse late events, or other model problems could also produce structure. Check the event-time distribution and the number at risk over time before interpreting a late pattern.

Adding a Time-by-Covariate Term

In crr, cov1 holds covariates with constant coefficients. cov2 holds covariates whose coefficients vary according to functions returned by tf. Each column of cov2 is paired with a column of the time-function matrix. Thus, to retain a treatment main effect and add a treatment-by-time term, include treatment in both matrices: once in cov1 for its main effect, and once in cov2 for its interaction.

A useful time function for positive event times is $\log(t/180)$, when time is measured in days and 180 days is a clinically interpretable reference. The model for binary treatment $x$ is then

\lambda_1^*(t\mid x,z) = \lambda_{10}^*(t) \exp\left\{\beta x+\gamma x\log(t/180)+\theta z\right\},

where $z$ is another covariate. At $t=180$ days, $\log(t/180)=0$, so the treatment subdistribution-hazard ratio is $\exp(\beta)$. At another time $t$, it is $\exp\{\beta+\gamma\log(t/180)\}$. A nonzero $\gamma$ means that this ratio changes with time under the selected log-time form.

Worked Example: Fitting and Interpreting a Log-Time Interaction

Continue with the hypothetical data above. Assume that the residual plot motivates checking whether the treatment effect changes over time. The treatment interaction is coded as follows; time is in days, so the function is defined for positive event times.

\texttt{X1 <- cbind(treatment = treatment, age10 = age10)}\\ \texttt{X2 <- cbind(treatment = treatment)}\\ \texttt{fit1 <- crr(ftime, fstatus, cov1 = X1, cov2 = X2,}\\ \texttt{\ \ tf = function(uft) cbind(log(uft / 180)),}\\ \texttt{\ \ failcode = 1, cencode = 0)}

Suppose, for illustration, the fitted treatment main-effect estimate is $\hat{\beta}=-0.18$, the time-interaction estimate is $\hat{\gamma}=0.32$, and the estimated standard error of $\hat{\gamma}$ is $0.14$. These are stipulated hypothetical fit results, not results from a real study. For one interaction coefficient, the Wald statistic is

z=\frac{\hat{\gamma}}{SE(\hat{\gamma})} =\frac{0.32}{0.14}=2.286.

A two-sided normal-reference p-value is approximately $0.0223$. At the 0.05 significance level, this is evidence against $H_0:\gamma=0$, conditional on the model and the specified time function. It is not the probability that $H_0$ is true. The approximate 95% Wald interval for $\gamma$ is $0.32\pm1.96(0.14)=[0.0456,0.5944]$, in log subdistribution-hazard-ratio units per unit increase in $\log(t/180)$.

At 180 days, the estimated treatment subdistribution-hazard ratio is $\exp(-0.18)=0.835$. At 540 days, $\log(540/180)=\log(3)=1.0986$, so the estimated log ratio is $-0.18+0.32(1.0986)=0.1716$, and the ratio is $\exp(0.1716)=1.187$. The fitted association thus changes from below 1 at 180 days to above 1 at 540 days. These are model-based subdistribution-hazard ratios, not estimates that treatment causes a change in cumulative incidence. For inference about the ratio at either time, its standard error must account for the covariance between $\hat{\beta}$ and $\hat{\gamma}$; the interaction’s standard error alone is not enough.

In an actual analysis, inspect the coefficient names and variance matrix in the fitted object to identify the interaction coefficient and its standard error. The test above is appropriate for one prespecified time interaction. With several interaction terms, test the full block using a Wald statistic $\hat{\boldsymbol{\gamma}}^\mathsf{T}\widehat{\operatorname{Var}}(\hat{\boldsymbol{\gamma}})^{-1}\hat{\boldsymbol{\gamma}}$, compared approximately with a chi-squared distribution with degrees of freedom equal to the number of interaction coefficients.

Choosing and Checking the Time Scale

Time functions should be chosen before examining many alternatives, and their units must be consistent. With $\log(t/180)$, both $t$ and 180 must be expressed in the same units. The reference time makes the main effect interpretable at that time; it does not change the fitted model when the same interaction is parameterized consistently.

Worked Example: Converting the Reference Time Correctly

If time is recorded in days, 540 days gives $\log(540/180)=\log(3)=1.0986$. Converting to years using 365.25 days per year gives $540/365.25=1.47844\ldots$, or about 1.478 years. The reference time must also be converted: $180/365.25=0.4928\ldots$ years. The same time multiplier is therefore

\log\left(\frac{1.47844\ldots}{180/365.25}\right) = \log(3) \approx 1.0986.

Using years in the numerator but retaining 180 as a days-based reference would instead calculate $\log(1.478/180)\approx-4.802$. That arithmetic is correct for the mismatched expression, but the expression is not the intended time function. It can produce a badly mis-scaled interaction. Convert the reference and all event times together, or retain the original time unit.

A logarithmic interaction is undefined at time zero, but crr evaluates the time function at cause-1 event times. Those times must be positive. If zero-time cause-1 events are possible, investigate the time recording and choose a defensible, prespecified transformation that is defined at every event time; do not silently add an arbitrary constant. Also inspect the fitted range: a model may estimate an interaction even where late follow-up is supported by very few events.

How to Use the Evidence

Fit and inspect the constant-effect model.
Check status coding, covariates, convergence, event-time support, and plots of the relevant columns of fit$res against fit$uftime.

Specify a focused time interaction.
Choose a time function that addresses the observed pattern and makes sense for the application. Include the covariate in cov1 for its main effect and in cov2 for its varying effect.

Test, interpret, and diagnose again.
Assess the interaction coefficient or block, report effects at meaningful times, and re-examine residual patterns and late-time support. Explain the chosen function and avoid treating a non-significant test as proof of proportionality.

Residual plots and time interactions answer related but distinct questions. The residual plot can reveal unanticipated patterns without first selecting a functional form. An interaction test quantifies evidence for one specified form, such as a linear trend in log time. A significant interaction does not establish that this is the correct temporal shape; a non-significant interaction does not rule out other forms. If a log-linear change is implausible, a prespecified alternative such as a piecewise effect or a flexible function may be considered, provided the data support the added complexity.

These diagnostics do not replace other model checks. A temporal pattern can arise from an omitted or misspecified covariate, influential observations, or coding errors, as well as from a genuinely time-varying association. Examine whether the pattern is driven by a small number of late events. The cause-specific Cox model addresses a different hazard estimand; switching models is not simply a correction for nonproportionality in a Fine-Gray fit. The present checks concern the Fine-Gray subdistribution-hazard model.

Common Mistakes and Limitations

Calling residuals direct prediction errors. The crr residuals are estimating-equation contributions at event times. Their plot is diagnostic; a point is not a subject-level difference between observed and predicted cumulative incidence.
Using mixed time units. The time argument and its reference must use the same units. A numerically computable interaction can still be scientifically and mathematically mis-specified if they do not.
Reading a smooth curve as a formal test. A visually rising or falling trend is suggestive, not a p-value. A formal interaction test depends on the chosen function and model assumptions.
Overinterpreting a non-significant interaction. Sparse events or low power can obscure time variation. “Did not detect evidence” is more accurate than “proved proportional hazards.”
Confusing a subdistribution-hazard ratio with a risk ratio. The ratio is not generally the ratio of cumulative incidence at a given time. Report the estimand precisely and give time-specific effects only with appropriate uncertainty calculations.
Searching many functions without accounting for selection. Trying numerous transformations and reporting only the most favorable result inflates false-positive risk. Prefer a prespecified scientific rationale and describe exploratory searches as such.

Key takeaway: Use the time-indexed residuals from crr to look for temporal structure, then use a scientifically motivated cov2-by-tf term to quantify a specified departure from proportional subdistribution hazards. Check units, event support, model assumptions, and the distinction between a subdistribution hazard and cumulative incidence.

Check Your Understanding

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

What does a systematic trend in a Fine-Gray residual plot suggest, and why is it not conclusive by itself?
For a model with $\beta=-0.10$ and $\gamma=0.20$ multiplying $x\log(t/100)$, calculate the fitted subdistribution-hazard ratio for a one-unit increase in $x$ at 100 and 200 time units.
If an interaction estimate is $0.24$ with standard error $0.12$, calculate its Wald statistic and approximate two-sided p-value. State the conclusion at $\alpha=0.05$.
Why is $\log(1.478/180)$ not the appropriate time multiplier for 540 days when time has been converted to years?

Solutions. (1) A trend suggests that a constant coefficient may not describe the time pattern, but noise, sparse events, omitted covariates, or other misspecification can also produce structure; the plot is not a formal test. (2) At 100 units the multiplier is zero, so the ratio is $\exp(-0.10)=0.905$. At 200 units, $\log(2)=0.6931$, giving $\exp\{-0.10+0.20(0.6931)\}=\exp(0.0386)=1.039$. (3) $z=0.24/0.12=2.00$, with two-sided $p\approx0.0455$. This is evidence against a zero interaction at the 0.05 level under the model and test assumptions. (4) It mixes years in the numerator with a reference of 180 days. The reference must be converted to $180/365.25$ years; then the multiplier is $\log\{1.47844\ldots/(180/365.25)\}=\log(3)\approx1.0986$.