Fitting Fine-Gray Models with cmprsk

What a Fine–Gray Fit Estimates

A Fine–Gray model relates covariates to the cumulative incidence function (CIF) for one event type in the presence of competing events. It is useful when the scientific question concerns the probability of experiencing a particular event by a given time, while recognizing that other event types can prevent that event from occurring. This tutorial focuses on fitting the model with R’s cmprsk::crr function and obtaining coefficient estimates and standard errors.

Let \(T\) be the time to the first event, \(J\) the event type, and \(k\) the event type of interest. The CIF for cause \(k\), conditional on covariates \(z\), is \(F_k(t\mid z)=P(T\leq t,J=k\mid z)\). Unlike a cause-specific Cox model, the Fine–Gray model is built directly around this CIF through the cause-\(k\) subdistribution hazard.

The conceptual risk set for this hazard differs from an ordinary Cox risk set: people who experience a competing event are retained in the subdistribution risk set, with weighting to account for censoring. Consequently, crr is not simply a Cox regression with competing events coded as censored. Its estimation procedure uses censoring adjustments, and the resulting coefficient is interpreted on the subdistribution hazard scale.

Data and Model Assumptions

For each person, prepare a follow-up time, a status code, and a row of covariates. The follow-up time is the time to the first event or censoring. The status code must distinguish the event of interest, each competing event, and censoring. In the call to crr, the arguments failcode and cencode identify the target event and the censoring code; every other non-censoring status represents a competing event.

These assumptions are substantive, not guaranteed by successful code execution. In particular, an estimate from a model with a time-varying covariate effect is a potentially misleading summary if interpreted as a constant SHR. The next tutorial addresses residuals and time interactions; here, the priority is to specify and fit the intended model correctly.

The crr function takes a covariate matrix rather than an R formula. Construct the matrix explicitly, check factor reference categories, and remove the intercept. Rows in the time, status, and covariate inputs must refer to the same people in the same order. Handle missingness consistently before fitting; do not allow the three inputs to silently represent different subsets.

Fitting and Reading the Output

A minimal call has the structure cmprsk::crr(ftime, fstatus, cov1, failcode, cencode). The returned object includes coefficient estimates in fit$coef and an estimated coefficient covariance matrix in fit$var. The estimated standard error for coefficient \(j\) is the square root of the corresponding diagonal element of that matrix:

A large-sample Wald statistic for testing \(H_0:\beta_j=0\) is \(z=\hat\beta_j/\widehat{SE}(\hat\beta_j)\). A two-sided p-value is calculated from the standard normal distribution as \(2P(Z\geq |z|)\), where \(Z\sim N(0,1)\). A corresponding approximate 95% confidence interval on the coefficient scale is \(\hat\beta_j\pm1.96\widehat{SE}(\hat\beta_j)\); exponentiating its endpoints gives a confidence interval for the SHR. These are large-sample Wald calculations, so they can be unreliable with sparse events or imprecise estimates.

The crr output also contains other components, including residual-related quantities. This fitting tutorial uses the coefficient and covariance components to obtain estimates and standard errors; it does not treat successful convergence as evidence that model assumptions hold.

Coefficient Scale, Units, and Coding Choices

A coefficient is tied to the units and coding of its covariate. A one-year age effect may be inconvenient to communicate, while a ten-year contrast may be more useful. If age has coefficient \(\hat\beta_{\text{age}}\), the model-implied SHR for a ten-year increase is \(\exp(10\hat\beta_{\text{age}})\), with other covariates held fixed. Its confidence interval is obtained by exponentiating \(10\) times the coefficient-scale interval endpoints. State the contrast and units rather than reporting a bare exponentiated coefficient.

For a factor, the interpretation depends on its reference category and contrast coding. Check the levels and columns in the design matrix before fitting. A mislabeled reference group does not change the underlying fit, but it changes how each coefficient must be described. Similarly, switching the target cause by changing failcode changes the model’s estimand; it is not a cosmetic change to the same analysis.

Why an SHR Is Not a Risk Ratio

The distinction follows directly from the model. Define the cumulative subdistribution hazard \(H_k^*(t\mid z)=\int_0^t\lambda_k^*(u\mid z)\,du\). Since \(\lambda_k^*=-d\log(1-F_k)/dt\), integration gives \(1-F_k(t\mid z)=\exp\{-H_k^*(t\mid z)\}\). Under the proportional model, \(H_k^*(t\mid z)=H_{k0}^*(t)\exp(\beta^\mathsf{T}z)\), so:

The coefficient multiplies the cumulative subdistribution hazard, but the CIF is a nonlinear transformation that also depends on time and the baseline hazard. Therefore, even if the SHR is constant over time, the ratio of CIFs generally is not. This derivation also explains why a positive coefficient corresponds to a larger modeled CIF at a fixed time, all else equal, while the magnitude of the risk difference depends on time and baseline risk.

A cause-specific Cox model instead models the instantaneous event rate among people still free of every event type. It can be appropriate for questions about that cause-specific rate, but it is not a direct substitute for a Fine–Gray fit when the target is the subdistribution hazard linked to the CIF. Neither coefficient should be described as a causal effect without additional design and causal assumptions.

Diagnostics and Common Failure Modes

Before interpreting output, check that the observed status values are exactly the intended codes, that failcode identifies the desired event, and that cencode identifies censoring. A wrong code can produce a numerically valid fit to the wrong event. Inspect the design matrix for an intercept, duplicated columns, unexpected missing values, and factor levels with no observations. Confirm that every column has a clear interpretation and that the reference groups are as intended.

Warnings, non-finite estimates, or very large standard errors can indicate sparse event information, collinearity, or a coding problem. A fit may also be numerically successful despite poor proportional subdistribution hazards or an unsuitable covariate form. Assess the proportionality assumption and residual behavior rather than assuming that software output validates the model. When censoring varies systematically across groups, consider whether censoring strata are needed and whether the specified cengroup structure is scientifically defensible.

• Wrong: “I coded competing events as censored, so the target event is handled correctly.” Correction: Supply their distinct event codes; with the selected failcode, they are treated as competing events by crr.
• Wrong: “The SHR is a risk ratio.” Correction: It is a ratio of subdistribution hazards; calculate or predict CIFs separately to discuss absolute risk.
• Wrong: “The function returned coefficients, so proportionality holds.” Correction: Estimation output does not establish model adequacy; use appropriate diagnostics.
• Wrong: “A coefficient for treatment means the treatment caused the difference.” Correction: It describes an adjusted model association unless a causal design and assumptions justify causal interpretation.