Germán Rodríguez
Survival Analysis Princeton University

Competing Risks: Cox Models

We continue our competing risk analysis of supreme court tenure by introducing two covariates: age, how old the justice was when appointed, and year, the calendar year of appointment.

We obviously expect survival to improve over time and decline with the age of the justice at appointment, but it is not clear whether the trend over time is to retire earlier or later.

Cause-Specific Hazards

We start with Cox regressions looking at overall tenure in the court

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

. stset tenure, fail(event) // event > 0

     failure event:  event != 0 & event < .
obs. time interval:  (0, tenure]
 exit on or before:  failure

------------------------------------------------------------------------------
        112  total observations
          0  exclusions
------------------------------------------------------------------------------
        112  observations remaining, representing
        103  failures in single-record/single-failure 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

. stcox age year, efron

         failure _d:  event
   analysis time _t:  tenure

Iteration 0:   log likelihood = -385.80479
Iteration 1:   log likelihood = -371.44771
Iteration 2:   log likelihood = -371.38331
Iteration 3:   log likelihood =  -371.3833
Refining estimates:
Iteration 0:   log likelihood =  -371.3833

Cox regression -- Efron method for ties

No. of subjects =          112                  Number of obs    =         112
No. of failures =          103
Time at risk    =         1858
                                                LR chi2(2)       =       28.84
Log likelihood  =    -371.3833                  Prob > chi2      =      0.0000

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |   1.085238   .0188957     4.70   0.000     1.048828    1.122912
        year |   .9939418   .0017327    -3.49   0.000     .9905515    .9973436
------------------------------------------------------------------------------
> library(foreign)

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

> library(survival)

> coxb <- coxph(Surv(tenure, event > 0) ~ age + year, data = sc)

> coef(summary(coxb))
             coef exp(coef)    se(coef)         z     Pr(>|z|)
age   0.081799727 1.0852384 0.017411527  4.698022 2.626938e-06
year -0.006076662 0.9939418 0.001743268 -3.485788 4.906890e-04

We see that the risk of leaving the court at any length of service is 8% higher for every year of age at appointment, and about half a percent lower per calendar year. This, of course, confounds trends in mortality with trends in retirement.

To look at the risk of death we use death as the event indicator, treating everything else as censoring

. streset , fail(event == 1) // just death
-> stset tenure, failure(event==1)

     failure event:  event == 1
obs. time interval:  (0, tenure]
 exit on or before:  failure

------------------------------------------------------------------------------
        112  total observations
          0  exclusions
------------------------------------------------------------------------------
        112  observations remaining, representing
         49  failures in single-record/single-failure 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

. stcox age year, efron

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

Iteration 0:   log likelihood = -185.71011
Iteration 1:   log likelihood = -173.76436
Iteration 2:   log likelihood = -173.72112
Iteration 3:   log likelihood = -173.72111
Refining estimates:
Iteration 0:   log likelihood = -173.72111

Cox regression -- Efron method for ties

No. of subjects =          112                  Number of obs    =         112
No. of failures =           49
Time at risk    =         1858
                                                LR chi2(2)       =       23.98
Log likelihood  =   -173.72111                  Prob > chi2      =      0.0000

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |   1.071606   .0255427     2.90   0.004     1.022694    1.122856
        year |   .9881385   .0026446    -4.46   0.000     .9829687    .9933355
------------------------------------------------------------------------------
> coxd <- coxph(Surv(tenure, event==1) ~ age + year, data = sc)

> coef(summary(coxd))
            coef exp(coef)    se(coef)         z     Pr(>|z|)
age   0.06915817 1.0716057 0.023835916  2.901427 3.714671e-03
year -0.01193245 0.9881385 0.002676374 -4.458439 8.255857e-06

We see that the risk of death while in the court is about 7% higher per year of age at appointment, and declines just over one percent per calendar year of appointment, reflecting general improvements in mortality over time.

To look at the risk of retirement we change the event of interest

. streset , fail(event == 2) // just retirement
-> stset tenure, failure(event==2)

     failure event:  event == 2
obs. time interval:  (0, tenure]
 exit on or before:  failure

------------------------------------------------------------------------------
        112  total observations
          0  exclusions
------------------------------------------------------------------------------
        112  observations remaining, representing
         54  failures in single-record/single-failure 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

. stcox age year, efron

         failure _d:  event == 2
   analysis time _t:  tenure

Iteration 0:   log likelihood = -202.65368
Iteration 1:   log likelihood = -193.95438
Iteration 2:   log likelihood = -193.69369
Iteration 3:   log likelihood = -193.69336
Refining estimates:
Iteration 0:   log likelihood = -193.69336

Cox regression -- Efron method for ties

No. of subjects =          112                  Number of obs    =         112
No. of failures =           54
Time at risk    =         1858
                                                LR chi2(2)       =       17.92
Log likelihood  =   -193.69336                  Prob > chi2      =      0.0001

------------------------------------------------------------------------------
          _t | Haz. Ratio   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
         age |    1.10455   .0282053     3.89   0.000     1.050629    1.161238
        year |   .9996366   .0024224    -0.15   0.881        .9949    1.004396
------------------------------------------------------------------------------
> coxr <- coxph(Surv(tenure, event == 2) ~ age + year, data = sc)

> coef(summary(coxr))
              coef exp(coef)    se(coef)          z     Pr(>|z|)
age   0.0994381242 1.1045501 0.025535539  3.8941072 9.856107e-05
year -0.0003634978 0.9996366 0.002423297 -0.1500014 8.807635e-01

We see an increase in the risk of retirement of about 10% per year of age at appointment, and essentially no trend over time, so the risk of retirement depends only on age and not on calendar year of appointment.

Incidence from Cox Models

The next step will be to calculate cumulative incidence functions based on the Cox models for cause-specific hazards. In the next section we will use the Fine and Gray model.

For convenience we center age on 55 and calendar year on 1950. We start by fitting a Cox model to the risk of death

. gen agec = age - 55

. gen yearc = year - 1950

. stset tenure, fail(event == 1)

     failure event:  event == 1
obs. time interval:  (0, tenure]
 exit on or before:  failure

------------------------------------------------------------------------------
        112  total observations
          0  exclusions
------------------------------------------------------------------------------
        112  observations remaining, representing
         49  failures in single-record/single-failure 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

. quietly stcox agec yearc, efron
> sc <- mutate(sc, agec = age - 55, yearc = year - 1955)

> coxdc <-  coxph(Surv(tenure, event == 1) ~ agec + yearc, data = sc)

Next we predict the baseline hazard, representing a justice appointed at age 55 in 1950, and the hazard for a justice appointed at the same age in 2000 Because we want just one hazard contribution per distinct failure time we set the superfluous values to missing To obtain the hazard we first use survfit() to obtain the survival and then difference minus its log

. predict h0_d, basehc
(63 missing values generated)

. gsort _t-_d

. by _t: replace h0_d = . if _n > 1 // set to missing for ties
(24 real changes made, 24 to missing)

. gen h1_d = h0_d * exp(50 * _b[yearc])
(87 missing values generated)
> sfd <- survfit(coxdc, newdata = list(agec = 0, yearc = 0))

> coxf <- data.frame(time = sfd$time, h0d = diff(c(0, -log(sfd$surv)))) %>%
+   mutate(h1d = h0d * exp(50 * coef(coxdc)["yearc"]))  

The next step is to do exactly the same calculations for the risk of retirement

. quietly streset, fail(event == 2)

. quietly stcox agec yearc, efron

. predict h0_r, basehc
(58 missing values generated)

. gsort _t - _d

. by _t: replace h0_r = . if _n > 1 // set to missing for ties
(27 real changes made, 27 to missing)

. gen h1_r = h0_r * exp(50 * _b[yearc])
(85 missing values generated)
> coxrc <-  coxph(Surv(tenure, event == 2) ~ agec + yearc, data = sc)

> sfr <- survfit(coxrc, newdata = list(agec = 0, yearc = 0))

> coxf <- mutate(coxf, h0r = diff(c(0, -log(sfr$surv))),
+         h1r = h0r * exp(50 * coef(coxrc)["yearc"])) 

The overall risk is just the sum of the risks of death and retirement. The survival function can be estimated as the product of the conditional survival probabilities for each failure time To compute a product in Stata we exponentiate the sum of the logs, after setting all missing values to zero. Observations where both hazards are missing don't contribute and can be dropped

. drop if missing(h0_d) & missing(h0_r)
(60 observations deleted)

. foreach h of varlist h0_d h1_d h0_r h1_r {
  2.         quietly replace `h'=0 if missing(`h')
  3. }

. gen S0 = exp(sum(log(1 - h0_d - h0_r)))

. gen S1 = exp(sum(log(1 - h1_d - h1_r)))
> coxf <- mutate(coxf,
+   S0 = cumprod(1 - h0d - h0r), 
+        S1 = cumprod(1 - h1d - h1r)) 

We can then calculate the CIF by computing the product of the cause-specific hazard at each failure time by the survival function just before that time, and summing.

. gen cif0 = sum(S0[_n-1] * h0_d )

. gen cif1 = sum(S1[_n-1] * h1_d )

. sort _t

. list _t cif0 cif1 in -1

     +--------------------------+
     | _t       cif0       cif1 |
     |--------------------------|
 52. | 36   .2748652   .1726615 |
     +--------------------------+
> last <- nrow(coxf)

> coxf <- mutate(coxf,
+   cif0 = cumsum(c(1, S0)[-last] * h0d),
+   cif1 = cumsum(c(1, S1)[-last] * h1d))

Finally we plot the curves

. line cif0 cif1 _t, connect(J J) ///
>    title("Cumulative incidence of death") xtitle(Tenure) ///
>    subtitle("Justice appointed at age 55 in 1950 and 2000") ///
>    note(Based on Cox model of cause-specific hazards) ///
>    text(.17 38 "2000") text(.27 38 "1950") legend(off)

. 
. graph export cif_h.png, width(500) replace
(file cif_h.png written in PNG format)

CIF Death

> g <- ggplot(coxf, aes(time, cif0)) + geom_line() + geom_line(aes(time, cif1)) +
+ ggtitle("Incidence of Death For Justices Appointed at Age 55 in 1950 and 2000") 

> g + annotate("text", x = 38, y=max(coxf$cif0), label="1950") +
+     annotate("text", x = 38, y=max(coxf$cif1), label="2000") 

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

CIF Death In the next section we will estimate these functions using the Fine and Gray model. The general appearance of the curves is very similar. Both approaches estimate a reduction of about 10 percentage points in a fifty year period in the probability that a justice will die while in the court.

Continue with Fine and Gray