## 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)

> 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)

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.