🧠 Introduction

In survival analysis, a competing risk is an event that either precludes the occurrence of the primary event of interest or modifies the probability of its occurrence. For example, when studying cardiovascular death, non-cardiovascular death is a competing risk — because once it occurs, the individual can no longer experience the event of interest.

In this simulation, we model: - The event of interest: cardiovascular death (event code = 1) - The competing event: non-cardiovascular death (event code = 2) - The covariates include: - age (continuous) - sex (binary: 0 = female, 1 = male) - treatment (binary: A or B), where treatment B reduces the risk of cardiovascular death

📦 Load required packages

# install.packages(c("riskRegression", "survival", "cmprsk", "data.table", "prodlim"))

library(survival)
library(cmprsk)
library(riskRegression)
## riskRegression version 2023.12.21
library(data.table)
library(prodlim)  # Needed for Hist()

🧪 Simulate data with competing risks

set.seed(123)

n <- 1000
age <- rnorm(n, mean = 65, sd = 10)
sex <- rbinom(n, 1, 0.5)  # 0 = female, 1 = male
treatment <- rbinom(n, 1, 0.5)  # 0 = A, 1 = B

lambda1 <- 0.03 * exp(0.03 * (age - 65) + 0.5 * sex - 1.0 * treatment)  # CV death
lambda2 <- 0.015 * exp(0.02 * (age - 65) - 0.3 * sex)  # Non-CV death

T1 <- rexp(n, rate = lambda1)
T2 <- rexp(n, rate = lambda2)
C <- runif(n, 0, 10)

time <- pmin(T1, T2, C)
status <- ifelse(time == T1, 1, ifelse(time == T2, 2, 0))

data <- data.frame(
  time = time,
  status = status,
  age = age,
  sex = sex,
  treatment = factor(treatment, levels = c(0, 1), labels = c("A", "B"))
)

📊 Cumulative Incidence Function (CIF) curves

cif_simple <- cuminc(ftime = data$time, fstatus = data$status)

plot(cif_simple,
     xlab = "Time (years)",
     ylab = "Cumulative incidence",
     lwd = 2,
     col = c("blue", "red"))

legend("topleft",
       legend = c("Cardiovascular death", "Non-cardiovascular death"),
       col = c("blue", "red"),
       lwd = 2)

Explanation:
The cumulative incidence function (CIF) plots show the probability over time of experiencing either cardiovascular death (blue) or non-cardiovascular death (red), while properly accounting for competing events. This provides a more accurate picture than Kaplan-Meier when more than one type of event is possible.

⚙️ Cause-Specific Cox Models

# Cox model for cardiovascular death
cox_cause1 <- coxph(Surv(time, status == 1) ~ age + sex + treatment, data = data)
summary(cox_cause1)
## Call:
## coxph(formula = Surv(time, status == 1) ~ age + sex + treatment, 
##     data = data)
## 
##   n= 1000, number of events= 114 
## 
##                 coef exp(coef)  se(coef)      z Pr(>|z|)    
## age         0.041952  1.042844  0.009356  4.484 7.33e-06 ***
## sex         0.320136  1.377316  0.189747  1.687 0.091570 .  
## treatmentB -0.683679  0.504757  0.198840 -3.438 0.000585 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## age           1.0428     0.9589    1.0239    1.0621
## sex           1.3773     0.7260    0.9496    1.9978
## treatmentB    0.5048     1.9812    0.3418    0.7453
## 
## Concordance= 0.635  (se = 0.028 )
## Likelihood ratio test= 35.96  on 3 df,   p=8e-08
## Wald test            = 36.12  on 3 df,   p=7e-08
## Score (logrank) test = 36.69  on 3 df,   p=5e-08
# Cox model for non-cardiovascular death
cox_cause2 <- coxph(Surv(time, status == 2) ~ age + sex + treatment, data = data)
summary(cox_cause2)
## Call:
## coxph(formula = Surv(time, status == 2) ~ age + sex + treatment, 
##     data = data)
## 
##   n= 1000, number of events= 58 
## 
##                coef exp(coef) se(coef)      z Pr(>|z|)  
## age         0.03412   1.03471  0.01329  2.567   0.0103 *
## sex        -0.32392   0.72331  0.26688 -1.214   0.2249  
## treatmentB -0.18669   0.82970  0.26417 -0.707   0.4798  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##            exp(coef) exp(-coef) lower .95 upper .95
## age           1.0347     0.9665    1.0081     1.062
## sex           0.7233     1.3825    0.4287     1.220
## treatmentB    0.8297     1.2053    0.4944     1.392
## 
## Concordance= 0.607  (se = 0.036 )
## Likelihood ratio test= 8.44  on 3 df,   p=0.04
## Wald test            = 8.47  on 3 df,   p=0.04
## Score (logrank) test = 8.49  on 3 df,   p=0.04

Explanation:
In cause-specific models, each event type is modeled separately using a Cox proportional hazards model. Other events are treated as censored. These models estimate the instantaneous risk of each event and are useful for understanding etiological effects of covariates.

Interpretation: We compare the effect of Treatment B on the cause-specific hazards for:

Cardiovascular Death (status == 1)

  • Hazard Ratio: exp(coef) = 0.505
  • 95% CI: [0.342, 0.745]
  • p-value: 0.0006 (statistically significant)

Interpretation:
Treatment B is associated with a 49.5% reduction in the cause-specific hazard of cardiovascular death compared to Treatment A. This effect is statistically significant, suggesting a protective benefit of Treatment B for this specific cause.

Non-Cardiovascular Death (status == 2)

  • Hazard Ratio: exp(coef) = 0.830
  • 95% CI: [0.494, 1.392]
  • p-value: 0.48 (not statistically significant)

Interpretation:
For non-cardiovascular death, Treatment B is associated with a 17% reduction in the hazard, but this effect is not statistically significant. We cannot conclude that Treatment B has any meaningful impact on this cause of death.

🔢 Fine-Gray Subdistribution Hazard Model

dataDT <- as.data.table(data)

fg_model <- FGR(Hist(time, status) ~ age + sex + treatment, data = dataDT, cause = 1)
summary(fg_model)
## Competing Risks Regression
## 
## Call:
## FGR(formula = Hist(time, status) ~ age + sex + treatment, data = dataDT, 
##     cause = 1)
## 
##               coef exp(coef) se(coef)     z p-value
## age         0.0401     1.041   0.0085  4.72 2.4e-06
## sex         0.3212     1.379   0.1863  1.72 8.5e-02
## treatmentB -0.6682     0.513   0.1986 -3.36 7.7e-04
## 
##            exp(coef) exp(-coef)  2.5% 97.5%
## age            1.041      0.961 1.024 1.058
## sex            1.379      0.725 0.957 1.986
## treatmentB     0.513      1.951 0.347 0.757
## 
## Num. cases = 1000
## Pseudo Log-likelihood = -712 
## Pseudo likelihood ratio test = 34.3  on 3 df,

Explanation:
The Fine-Gray model estimates the subdistribution hazard for cardiovascular death (event 1), directly modeling the cumulative incidence rather than the instantaneous hazard. Unlike cause-specific Cox models, this model keeps individuals with competing events in the risk set, adjusting their weights, and is useful when estimating absolute risk.

Interpretation:

The summary(fg_model) output includes estimated coefficients (coef) and their exponentiated values (exp(coef)), which are interpreted as subdistribution hazard ratios (SHR).

For example:

📝 Conclusion

This document demonstrates how to simulate competing risks data, visualize cumulative incidence, and fit both cause-specific and Fine-Gray models. Each method offers a different perspective: Cox models are ideal for causal inference, while Fine-Gray models are better suited for risk prediction.