Chapter 3 Wednesday: Competing risk

3.1 Key (operative) concepts

  1. Patient are exposed simultaneously to \(k (\geq2)\) causes

  2. Effect Free Survival (EFS) is univariate, i.e. only the First Observed Event (FOE) is considered and of interest

  3. The interest is not in the survival model

“At \(\infty\) all individuals will not die in the ICU”

  1. Type of observed time
  • Cansored (conventionally coded w/ \(0\))
  • Failure w/ a FOE different from the last absorbing one (coded w/ \(1 -- k-1\))
  • Failure w/ the FOE as the last absorbing event (coded w/ \(k\))
  • \(T_k = min\Bigl(\tilde{T_k}^{D_k} | D_k \in \{\textrm{causes of failure for }k\} \Bigr)\)
  1. Cumulative Incidence Function (CIF) do not require independence between causes

  2. In competing risk, K-M is biased, i.e. overestimates the CIF (because it the independence assumption is violated)

  3. Tests
  • w/out competing risk: log-rank
  • w/ competing risk: modified \(\chi^2\) (Gray 1988)
  1. Regression strategies for competing risk
  • Case Specific Hazard Ratio (CS-HR) — Cox, useful for clinical interests (present it for each competing risk taken singularly)
  • Subdistribution Hazard Ratio (SHR) — Fine-Gray, useful for administrative] interests (present it for the global risk considering the contribution of each competing one)

Test the proportional hazard assumption for SHR

There are formulas for the sample size calculation when considering competing risk

3.2 Data manipulation

set.seed(171004)
data(mgus, package = 'survival')
# ?mgus

mgus_df <- as_tibble(mgus)
dd <- datadist(mgus_df)

mgus_df
## # A tibble: 241 x 12
##       id   age    sex  dxyr   pcdx pctime futime death   alb creat   hgb
##  * <dbl> <dbl> <fctr> <dbl> <fctr>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
##  1     1    78 female    68   <NA>     NA    748     1   2.8   1.2  11.5
##  2     2    73 female    66     LP   1310   6751     1    NA    NA    NA
##  3     3    87   male    68   <NA>     NA    277     1   2.2   1.1  11.2
##  4     4    86   male    69   <NA>     NA   1815     1   2.8   1.3  15.3
##  5     5    74 female    68   <NA>     NA   2587     1   3.0   0.8   9.8
##  6     6    81   male    68   <NA>     NA    563     1   2.9   0.9  11.5
##  7     7    72 female    68   <NA>     NA   1135     1   3.0   0.8  13.5
##  8     8    79 female    69   <NA>     NA   2016     1   3.1   0.8  15.5
##  9     9    85 female    70   <NA>     NA   2422     1   3.2   1.0  12.4
## 10    10    58   male    65   <NA>     NA   6155     1   3.5   1.0  14.8
## # ... with 231 more rows, and 1 more variables: mspike <dbl>
  1. Find number of patient w/ malignancy (AKA transition), death (w/out malignancy) and Free of Events.
mgus_df <- mgus_df %>% 
  mutate(
    malignancy = !is.na(pcdx)
  )

mgus_df %>% 
  group_by(malignancy, death) %>% 
  summarise(n = n())
## # A tibble: 4 x 3
## # Groups:   malignancy [?]
##   malignancy death     n
##        <lgl> <dbl> <int>
## 1      FALSE     0    14
## 2      FALSE     1   163
## 3       TRUE     0     2
## 4       TRUE     1    62

Patients w/ malignancy as a FOE are 64; patients which experiment death as FOE are 163, while the ones FoE are 14. 163.

  1. Find the indicator for censored, malignancy and death (indicator)
  2. Find the time-to-event to use in the models (time_t)
mgus_df <- mgus_df %>% 
  mutate(
    indicator = if_else(malignancy, 1, 2 * death),
    time_t    = pmin(futime, pctime, na.rm = TRUE)
)

mgus_df
## # A tibble: 241 x 15
##       id   age    sex  dxyr   pcdx pctime futime death   alb creat   hgb
##    <dbl> <dbl> <fctr> <dbl> <fctr>  <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
##  1     1    78 female    68   <NA>     NA    748     1   2.8   1.2  11.5
##  2     2    73 female    66     LP   1310   6751     1    NA    NA    NA
##  3     3    87   male    68   <NA>     NA    277     1   2.2   1.1  11.2
##  4     4    86   male    69   <NA>     NA   1815     1   2.8   1.3  15.3
##  5     5    74 female    68   <NA>     NA   2587     1   3.0   0.8   9.8
##  6     6    81   male    68   <NA>     NA    563     1   2.9   0.9  11.5
##  7     7    72 female    68   <NA>     NA   1135     1   3.0   0.8  13.5
##  8     8    79 female    69   <NA>     NA   2016     1   3.1   0.8  15.5
##  9     9    85 female    70   <NA>     NA   2422     1   3.2   1.0  12.4
## 10    10    58   male    65   <NA>     NA   6155     1   3.5   1.0  14.8
## # ... with 231 more rows, and 4 more variables: mspike <dbl>,
## #   malignancy <lgl>, indicator <dbl>, time_t <dbl>
  1. Estimate the naive K-M and the cumulative incidence functions
# Using survival
cuminc(mgus_df$time_t, mgus_df$indicator) %>%                  # ?cmprsk::cuminc
  plot(                                                  # ?cmprsk:::plot.cuminc
    main    = 'Cumulative Incidence Estimates curves',
    col     = c('blue', 'red'),
    xlab    = 'Days',
    curvlab = c('Transition', 'Death'),
    wh      = c(1, 1)                                          # legend position
)

survfit(Surv(time_t, malignancy) ~ 1,        # using `rms::npsurv()` is the same
  data = mgus_df
  ) %>% 
  lines(                          # Use `lines()` to draw over the previous plot
    fun       = 'event',                            # plot the cumulative events
    conf.int   = FALSE,
    col       = 'black',
    lty       = 3
)
legend(x = 1, y = 0.86,
  legend = 'Naive K-M (Overestimation of Transitions)',
  col    = 'black',
  lty    = 3,
  bty    = 'n' # remove box arround the legend (because we have to add an entry)
)

3.3 Simulation of Competing risk

  1. Specify two cause-specific exponential hazard \(\lambda_1(t)\) and \(\lambda_2(t)\) of means \(0.8\) and \(1.2\). (Set sample size as you like.)
n        <- 1e4
lambda_1 <- 0.8
lambda_2 <- 1.2
  1. Simulate survival times \(T\) based on the all causes hazard \(\lambda_.(t) = \lambda_1(t) + \lambda_2(t)\).
lambda    <- lambda_1 + lambda_2 
surv_time <- rexp(n,
  rate = 1 / lambda
)
  1. Generate Bernoulli \(B(p)\) random variables, w/ \(p = \lambda_1(t) / \lambda_.(t)\), i.e. is the probability of occurrence of the event of type 1.
p_cens     <- lambda_1 / lambda
transition <- rbinom(n,
  size = 1,
  prob = p_cens
) %>%
  as.logical     # Set as logical to use the variable for conditional statements
  1. Simulate uniform censoring times over \([0, 1]\).
censor_time <- runif(n,
  min = 0,
  max = 1
)
  1. Estimate the Cumulative Incidence of each competing event, w/ and w/out censoring; discuss the results.
# create the dataset
sim_data <- data_frame(
  id         = seq_len(n),
  transition = transition,
  surv_t     = surv_time,
  cens_t     = censor_time,
  time_t     = pmin(surv_t, cens_t),
  status     = case_when(
    time_t == cens_t ~ 0L,              # All the censored patients has status 0
    transition       ~ 1L,    # Among the other, the ones which has a transition
                              # have state 1
    TRUE             ~ 2L     # All the other were dead (before the end of f-up)
  )
)

# Explore a (random) sample of three cases for each staus
sim_data %>% 
  group_by(status) %>% 
  sample_n(3)
## # A tibble: 9 x 6
## # Groups:   status [3]
##      id transition      surv_t    cens_t      time_t status
##   <int>      <lgl>       <dbl>     <dbl>       <dbl>  <int>
## 1  6759      FALSE 6.844765752 0.5175084 0.517508388      0
## 2   906       TRUE 1.378642827 0.5063402 0.506340163      0
## 3  3568       TRUE 0.623047318 0.5752797 0.575279657      0
## 4  2196       TRUE 0.506380841 0.9344626 0.506380841      1
## 5  8119       TRUE 0.157231928 0.9007847 0.157231928      1
## 6  5854       TRUE 0.742379822 0.8771130 0.742379822      1
## 7  1685      FALSE 0.539729742 0.9635529 0.539729742      2
## 8  4550      FALSE 0.001042022 0.3692372 0.001042022      2
## 9  6893      FALSE 0.410360153 0.9322769 0.410360153      2
# Using survival
cuminc(sim_data$time_t, sim_data$status) %>%                  # ?cmprsk::cuminc
  plot(                                                  # ?cmprsk:::plot.cuminc
    main    = 'Cumulative Incidence Estimates curves',
    col     = c('blue', 'red'),
    xlab    = 'Time (normalized [0, 1])',
    curvlab = c('Transition', 'Event'),
    wh      = c(0.01, 1)                                          # legend position
)

survfit(Surv(time_t, transition) ~ 1,        # using `rms::npsurv()` is the same
  data = sim_data
  ) %>% 
  lines(                          # Use `lines()` to draw over the previous plot
    fun       = 'event',                            # plot the cumulative events
    conf.int   = FALSE,
    col       = 'black',
    lty       = 3
)
legend(x = 0.01, y = 0.86,
  legend = 'Naive K-M (overestimation of transitions)',
  col    = 'black',
  lty    = 3,
  bty    = 'n' # remove box arround the legend (because we have to add an entry)
)

3.4 Estimation of the effect of sex on MGUS incidence

  1. Compare the results of Cox cause specific hazard model…

For clinical questions, i.e. cause specific risk to experiment the event w/out taking into account the other couse(s)

dd <- datadist(mgus_df)

cox_sex <- cph(Surv(time_t, malignancy) ~ sex,
  data = mgus_df,
  x    = TRUE,
  y    = TRUE
)

summary(cox_sex)                  # this is good for a clean view of the effects
##              Effects              Response : Surv(time_t, malignancy) 
## 
##  Factor            Low High Diff. Effect   S.E.    Lower 0.95 Upper 0.95
##  sex - female:male 2   1    NA    0.049342 0.25103 -0.44268   0.54136   
##   Hazard Ratio     2   1    NA    1.050600      NA  0.64232   1.71830
cox_sex                    # Here there are more informations (and the p-values)
## Cox Proportional Hazards Model
##  
##  cph(formula = Surv(time_t, malignancy) ~ sex, data = mgus_df, 
##      x = TRUE, y = TRUE)
##  
##                      Model Tests       Discrimination    
##                                           Indexes        
##  Obs       241    LR chi2      0.04    R2       0.000    
##  Events     64    d.f.            1    Dxy     -0.039    
##  Center -0.028    Pr(> chi2) 0.8441    g        0.024    
##                   Score chi2   0.04    gr       1.025    
##                   Pr(> chi2) 0.8441                      
##  
##           Coef    S.E.   Wald Z Pr(>|Z|)
##  sex=male -0.0493 0.2510 -0.20  0.8442  
## 
Predict(cox_sex) %>%      # It is necessary to have the predictions for the plot
  plot

cox_sex_death <- cph(Surv(time_t, indicator == 2) ~ sex,
  data = mgus_df,
  x    = TRUE,
  y    = TRUE
)

summary(cox_sex_death)
##              Effects              Response : Surv(time_t, indicator == 2) 
## 
##  Factor            Low High Diff. Effect   S.E.    Lower 0.95 Upper 0.95
##  sex - female:male 2   1    NA    -0.44221 0.16183 -0.75939   -0.12502  
##   Hazard Ratio     2   1    NA     0.64262      NA  0.46795    0.88248
cox_sex_death
## Cox Proportional Hazards Model
##  
##  cph(formula = Surv(time_t, indicator == 2) ~ sex, data = mgus_df, 
##      x = TRUE, y = TRUE)
##  
##                      Model Tests       Discrimination    
##                                           Indexes        
##  Obs       241    LR chi2      7.65    R2       0.031    
##  Events    163    d.f.            1    Dxy      0.124    
##  Center 0.2514    Pr(> chi2) 0.0057    g        0.218    
##                   Score chi2   7.59    gr       1.243    
##                   Pr(> chi2) 0.0059                      
##  
##           Coef   S.E.   Wald Z Pr(>|Z|)
##  sex=male 0.4422 0.1618 2.73   0.0063  
## 
Predict(cox_sex_death) %>% 
  plot

  1. …to those of the Fine and Gray model

For administrative questions, i.e. overall risk of experiment each event taking into account the competing risk

mgus_num <- mgus_df %>% 
  mutate(sex = as.numeric(sex))

crr(                                      # We do not have a plot method for crr
  ftime   = mgus_num$time_t,
  fstatus = mgus_num$indicator,
  cov1    = mgus_num$sex
) %>%
  summary
## Competing Risks Regression
## 
## Call:
## crr(ftime = mgus_num$time_t, fstatus = mgus_num$indicator, cov1 = mgus_num$sex)
## 
##                 coef exp(coef) se(coef)     z p-value
## mgus_num$sex1 -0.339     0.713    0.249 -1.36    0.17
## 
##               exp(coef) exp(-coef)  2.5% 97.5%
## mgus_num$sex1     0.713        1.4 0.437  1.16
## 
## Num. cases = 241
## Pseudo Log-likelihood = -341 
## Pseudo likelihood ratio test = 1.83  on 1 df,

References

Gray, Robert J. 1988. “A Class of K-Sample Tests for Comparing the Cumulative Incidence of a Competing Risk.” The Annals of Statistics. JSTOR, 1141–54.