Germán Rodríguez
Survival Analysis Princeton University

Competing Risks: The Fine and Gray Model

We conclude our competing risk analysis of supreme court tenure by fitting the Fine and Gray model, which focuses on the sub-hazard of each risk.

This model can be fit in Stata using the stcrreg command with death as the event of interest as retirement as the competing event.

This model can be fit in R using the competing risk regression function `crr()`. The function requires a model matrix instead of a model formula; in our example the matrix is simple, for more general helps type ?model.matrix. I specify the event of interest and the censoring code, although they happen to coincide with the defaults.

. use http://data.princeton.edu/pop509/justices
(Supreme court justices 1-Jan-2016)

. streset , fail(event==1)
-> stset tenure, id(id) failure(event==1)

                id:  id
     failure event:  event == 1
obs. time interval:  (tenure[_n-1], tenure]
 exit on or before:  failure

------------------------------------------------------------------------------
        112  total observations
          0  exclusions
------------------------------------------------------------------------------
        112  observations remaining, representing
        112  subjects
         49  failures in single-failure-per-subject data
       1858  total analysis time at risk and under observation
                                                at risk from t =         0
                                     earliest observed entry t =         0
                                          last observed exit t =        36

. stcrreg age year, compete(event==2)

         failure _d:  event == 1
   analysis time _t:  tenure
                 id:  id

Iteration 0:   log pseudolikelihood = -210.85503  
Iteration 1:   log pseudolikelihood = -207.50164  
Iteration 2:   log pseudolikelihood = -207.43579  
Iteration 3:   log pseudolikelihood = -207.43578  

Competing-risks regression                       No. of obs       =        112
                                                 No. of subjects  =        112
Failure event  : event == 1                      No. failed       =         49
Competing event: event == 2                      No. competing    =         54
                                                 No. censored     =          9

                                                 Wald chi2(2)     =      16.11
Log pseudolikelihood = -207.43578                Prob > chi2      =     0.0003

                                   (Std. Err. adjusted for 112 clusters in id)
------------------------------------------------------------------------------
             |               Robust
          _t |        SHR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |   1.012261   .0182725     0.68   0.500     .9770735    1.048715
        year |   .9909425   .0023241    -3.88   0.000     .9863978    .9955081
------------------------------------------------------------------------------
> library(foreign)

> sc <- read.dta("http://data.princeton.edu/pop509/justices.dta")

> library(cmprsk)

> fg <- crr(sc$tenure, sc$event, sc[,c("age", "year")], 
+   failcode = 1, cencode = 0)

> summary(fg)
Competing Risks Regression

Call:
crr(ftime = sc$tenure, fstatus = sc$event, cov1 = sc[, c("age", 
    "year")], failcode = 1, cencode = 0)

        coef exp(coef) se(coef)      z p-value
age   0.0122     1.012  0.01797  0.678 5.0e-01
year -0.0091     0.991  0.00233 -3.897 9.7e-05

     exp(coef) exp(-coef)  2.5% 97.5%
age      1.012      0.988 0.977 1.049
year     0.991      1.009 0.986 0.995

Num. cases = 112
Pseudo Log-likelihood = -207 
Pseudo likelihood ratio test = 14.3  on 2 df,

We see that the probability that a judge will die while in the court doesn't appear to vary with age at appointment but decreases with calendar year of appointment, so judges appointed more recently are less likely to die on the job.

Predicted Cumulative Incidence

We predict the CIF of death for a justice appointed at age 55 in 1950 and in 2000

. stcurve, cif at1(age = 55 year = 1950) at2(age = 55 year = 2000) ///
>         title("Cumulative incidence of death") xtitle(Tenure) ///
>         subtitle("Justice appointed at age 55 in 1950 and 2000") ///
>         text(.20 38 "2000") text(.30 38 "1950") legend(off)

. graph export justices-c55y.png, width(500) replace
(file justices-c55y.png written in PNG format)

CIF Death

> nd <- matrix(c(55, 55, 1950, 2000), 2, 2)

> cifs <- predict(fg, nd)

> cifd <- data.frame( 
+   year = factor(rep(c("1950","2000"),rep(nrow(cifs),2))),
+   time = c(cifs[,1], cifs[,1]), 
+   cif = c(cifs[,2], cifs[,3]))

> ggplot(cifd, aes(time, cif, color=year)) + geom_step()

> ggsave("justices-c55yr.png", width = 500/72, height = 400/72, dpi = 72)

CIF Death

Aside: It is instructive to compute these curves 'by hand'. Because predict sets the covariates to zero, it will be advantageous to center age on 55 and year on 1950, as we did for the Cox model, so the predicted survival is our first scenario.

. gen agec = age - 55

. gen yearc = year - 1950

. stcrreg agec yearc, compete(event==2)

         failure _d:  event == 1
   analysis time _t:  tenure
                 id:  id

Iteration 0:   log pseudolikelihood = -210.85503  
Iteration 1:   log pseudolikelihood = -207.50164  
Iteration 2:   log pseudolikelihood = -207.43579  
Iteration 3:   log pseudolikelihood = -207.43578  

Competing-risks regression                       No. of obs       =        112
                                                 No. of subjects  =        112
Failure event  : event == 1                      No. failed       =         49
Competing event: event == 2                      No. competing    =         54
                                                 No. censored     =          9

                                                 Wald chi2(2)     =      16.11
Log pseudolikelihood = -207.43578                Prob > chi2      =     0.0003

                                   (Std. Err. adjusted for 112 clusters in id)
------------------------------------------------------------------------------
             |               Robust
          _t |        SHR   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        agec |   1.012261   .0182725     0.68   0.500     .9770735    1.048715
       yearc |   .9909425   .0023241    -3.88   0.000     .9863978    .9955081
------------------------------------------------------------------------------

. predict cif1950, basecif

We know that this model is linear in the c-log-log scale, so we transform the CIF for 1950, add 50 times the year coefficient, and then transform back to probabilities.

. gen cif2000 = invcloglog( cloglog(cif1950) + 50 * _b[yearc] )
(2 missing values generated)

. sort _t

. line cif1950 _t, c(J) || line cif2000 _t, c(J) legend(off)

. list _t cif1950 cif2000 in -1

     +-------------------------+
     | _t    cif1950   cif2000 |
     |-------------------------|
112. | 36   .2954049   .199208 |
     +-------------------------+

The last values of the CIF for death by year of appointment are as follows

> group_by(cifd, year) %>% summarize(p = max(cif))
Source: local data frame [2 x 2]

    year         p
  (fctr)     (dbl)
1   1950 0.2954049
2   2000 0.1992080

We see that the probability of dying while serving is declining over time. For a justice appointed at age 55, the probability would be 29.5% if appointed in 1950 and 19.9 if appointed in 2000, assuming of course that current trends continue, so the model applies.

Incidence functions can also be computed from Cox models, which in my view have a clearer interpretation, as shown here. The Cox approach leads to estimated probabilities of death of 27.5% if appointed in 1950 and 17.3% if appointed in 2000, under the proportional hazards assumption. Both methods agree in estimating the overall decline as approximately ten percentage points.