Competing Risks:
Singer and Willet (2003) Applied Longitudinal Data Analysis, Oxford, analyze data on the length of service of Supreme Court Justices. The dataset is available in Stata format from UCLA at http://www.ats.ucla.edu/stat/examples/alda. You have to download a zip file and then extract judges.dta. I have saved a copy in the course website, so you can
. use http://data.princeton.edu/pop509/judges.dta
. Desc
Contains data from http://data.princeton.edu/pop509/judges.dta
obs: 109
vars: 7
size: 3,488 (99.9% of memory free)
-------------------------------------------------------------------------------
storage display value
variable name type format label variable label
-------------------------------------------------------------------------------
id float %9.0g
tenure float %9.0g
dead float %9.0g
retire float %9.0g
leave float %9.0g
age float %9.0g
year float %9.0g
-------------------------------------------------------------------------------
Sorted by:
. Tab dead retire
| retire
dead | 0 1 | Total
-----------+----------------------+----------
0 | 9 53 | 62
1 | 47 0 | 47
-----------+----------------------+----------
Total | 56 53 | 109
As you can see we have information on 109 judges, of whom 53 have retired, 47 have died, and 9 were serving in 1999 when data collection ended (and indeed 9 would usually be censored).
There are two predictors of interest: age of the judge when nominated, and year when nominated. Following Singer and Willet we will treat both variables linearly. You may want to explore if this assumption is consistent with the data.
We first look at overall risk of leaving the court:
. stset tenure, fail(leave)
failure event: leave != 0 & leave < .
Obs. Time interval: (0, tenure]
exit on or before: failure
------------------------------------------------------------------------------
109 total obs.
0 exclusions
------------------------------------------------------------------------------
109 obs. Remaining, representing
100 failures in single record/single failure data
1740 total analysis time at risk, at risk from t = 0
earliest observed entry t = 0
last observed exit t = 36
. sts graph
failure _d: leave
analysis time _t: tenure
. graph export jsurv.png, replace width(400)
(note: file jsurv.png not found)
(file jsurv.png written in PNG format)
The median stay in the court is about 15 years (you can verify this using
sts list to obtain more details). The survival curve looks
pretty linear. What does that say about the overall risk of leaving?
Next we fit a proportional hazards model:
. stcox age year, efron
failure _d: leave
analysis time _t: tenure
Iteration 0: log likelihood = -369.81014
Iteration 1: log likelihood = -356.01736
Iteration 2: log likelihood = -355.88236
Iteration 3: log likelihood = -355.88231
Refining estimates:
Iteration 0: log likelihood = -355.88231
Cox regression -- Efron method for ties
No. of subjects = 109 Number of obs = 109
No. of failures = 100
Time at risk = 1740
LR chi2(2) = 27.86
Log likelihood = -355.88231 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. Z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
age | 1.089874 .0192433 4.87 0.000 1.052803 1.128251
year | .9947939 .0018698 -2.78 0.005 .9911359 .9984655
------------------------------------------------------------------------------
. estimates store leave
We see that the risk of leaving the court at any duration of service increases about 9.0% per year of age when the judge is nominated. The risk also declines about half a percent per year, indicating a trend for judges to stay longer over time.
See the Cox log for various diagnostics you could apply to these data.
The Risk of Retiring
We now look at the risk of retiring by resetting the data so only
retirement counts as a "failure".
This can be done using streset, so we only specify
what we are changing:
. streset, fail(retire)
-> stset tenure, failure(retire)
failure event: retire != 0 & retire < .
Obs. Time interval: (0, tenure]
exit on or before: failure
------------------------------------------------------------------------------
109 total obs.
0 exclusions
------------------------------------------------------------------------------
109 obs. Remaining, representing
53 failures in single record/single failure data
1740 total analysis time at risk, at risk from t = 0
earliest observed entry t = 0
last observed exit t = 36
. stcox age year, efron
failure _d: retire
analysis time _t: tenure
Iteration 0: log likelihood = -197.70313
Iteration 1: log likelihood = -187.16616
Iteration 2: log likelihood = -186.74899
Iteration 3: log likelihood = -186.74827
Refining estimates:
Iteration 0: log likelihood = -186.74827
Cox regression -- Efron method for ties
No. of subjects = 109 Number of obs = 109
No. of failures = 53
Time at risk = 1740
LR chi2(2) = 21.91
Log likelihood = -186.74827 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. Z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
age | 1.111944 .0286727 4.12 0.000 1.057143 1.169586
year | 1.000714 .0026451 0.27 0.787 .9955435 1.005912
------------------------------------------------------------------------------
The risk of retirement increases 11.2% per year of age at nomination; it does not appear to have changed over time; if anything it has increased slightly, but the difference is far from significant.
I would like to save the baseline cumulative hazard, which by default is evaluated when all predictors are zero. In this case that is far outside the range of the data and doesn't make sense. I will center age on 50 years and calendar year on 1900, so the baseline will refer to a 50-year old judge nominated in 1900. This doesn't change the coefficients at all, as you will see when I refit the model below.
Confusingly, the stcurve command can plot the baseline
survival, cumulative hazard and (smoothed) hazard, but it re-scales
it to the mean of all predictors. If you are going to use this
command you don't need to scale things, but I will do my own plot.
. gen agec = age - 50
. gen yearc = year - 1900
. stcox agec yearc, efron basech(Hret)
failure _d: retire
analysis time _t: tenure
Iteration 0: log likelihood = -197.70313
Iteration 1: log likelihood = -187.16616
Iteration 2: log likelihood = -186.74899
Iteration 3: log likelihood = -186.74827
Refining estimates:
Iteration 0: log likelihood = -186.74827
Cox regression -- Efron method for ties
No. of subjects = 109 Number of obs = 109
No. of failures = 53
Time at risk = 1740
LR chi2(2) = 21.91
Log likelihood = -186.74827 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. Z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
agec | 1.111944 .0286727 4.12 0.000 1.057143 1.169586
yearc | 1.000714 .0026451 0.27 0.787 .9955435 1.005912
------------------------------------------------------------------------------
. label var Hret "Cum hazard of retirement"
. estimates store retire
Risk of Death
We now reset things so death is the only cause of failure and fit a proportional hazards model.
. streset, fail(dead)
-> stset tenure, failure(dead)
failure event: dead != 0 & dead < .
Obs. Time interval: (0, tenure]
exit on or before: failure
------------------------------------------------------------------------------
109 total obs.
0 exclusions
------------------------------------------------------------------------------
109 obs. Remaining, representing
47 failures in single record/single failure data
1740 total analysis time at risk, at risk from t = 0
earliest observed entry t = 0
last observed exit t = 36
. stcox agec yearc, basech(Hdead)
failure _d: dead
analysis time _t: tenure
Iteration 0: log likelihood = -176.71579
Iteration 1: log likelihood = -167.55139
Iteration 2: log likelihood = -167.5303
Iteration 3: log likelihood = -167.5303
Refining estimates:
Iteration 0: log likelihood = -167.5303
Cox regression -- Breslow method for ties
No. of subjects = 109 Number of obs = 109
No. of failures = 47
Time at risk = 1740
LR chi2(2) = 18.37
Log likelihood = -167.5303 Prob > chi2 = 0.0001
------------------------------------------------------------------------------
_t | Haz. Ratio Std. Err. Z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
agec | 1.067512 .0254077 2.74 0.006 1.018857 1.11849
yearc | .9889157 .002853 -3.86 0.000 .9833398 .9945232
------------------------------------------------------------------------------
. label var Hdead "Cum hazard of death"
. estimates store dead
The risk of death increases 6.9% per year of age at nomination, and has declined about one percent per calendar year, a notable 69% in the course of a century, reflecting general increases in survival.
We see that the reason why judges are staying longer in the court is simply because they are living longer.
Death vs Retirement
Here's a summary of the estimates:
. estimates table retire dead, eform t
----------------------------------------
Variable | retire dead
-------------+--------------------------
agec | 1.1119443 1.0675115
| 4.12 2.74
yearc | 1.0007144 .98891572
| 0.27 -3.86
----------------------------------------
legend: b/t
Here's a plot of the two cumulative hazards.
. twoway (scatter Hret _t if _t<36, c(J) sort ) /// > (scatter Hdead _t if _t<36, c(J) sort ) . graph export compHaz.png, replace width(400) (file compHaz.png written in PNG format)
For the first 15 years or so the two hazards are very similar, so leaving the bench is just as likely to be by death or by retirement. In later years the risk of retirement grows much higher than the risk of death, so judges are more likely to leave by retirement than die on the job. (I exclude the last time point, 36 years, because both cumulative hazards are very large and it obscures the comparison at lower durations.)
Exercise: consider fitting piece-wise exponential models using suitable duration categories, so you can estimate the two hazards as well as the probability that if a judges leaves it is by dead or retirement given age at nomination, year of nomination, and length of service.