R-TUTORIAL ON COMPETING RISKS

This topic is treated in ASAUR Sections 9.2.1 - 9.2.3. Below we have condensed this considerably, concentrating on the estimation of (cumulative) cause-specific hazards and cumulative incidence functions.  

The data are used in ASAUR Section 9.2.3. Below we have taken out only what is needed for the present tutorial.

Here is a summary of what has been done in the book:

We consider prostate cancer patients ages 80 and over diagnosed with stage T2 poorly differentiated prostate cancer. We define indicator variables “status.other” and “status.prost”, and then select the subset “prostateSurvival.highrisk” as follows, using the “prostate survival” data
from Example 1.4 p.7. 

We will be interested in the survival time 'survTime' (months) and 'status' which is 0 for censoring, 1 for death of prostate cancer, 2 for death of other causes

library(asaur)

prostateSurvival <- within(prostateSurvival, {
+ status.prost <- as.numeric({status == 1})
+ status.other <- as.numeric({status == 2})})
attach(prostateSurvival)
prostateSurvival.highrisk <- prostateSurvival[{{grade ==
"poor"} &
+ {stage=="T2"} & {ageGroup == "80+"}},]
head(prostateSurvival.highrisk)

# The last stament gives output:

> head(prostateSurvival.highrisk)
   grade stage ageGroup survTime status status.other status.prost
13  poor    T2      80+       21      0            0            0
38  poor    T2      80+      105      0            0            0
41  poor    T2      80+        2      1            0            1
47  poor    T2      80+       67      2            1            0
78  poor    T2      80+        2      0            0            0
93  poor    T2      80+       60      2            1            0

tt=prostateSurvival.highrisk$survTime
dd=prostateSurvival.highrisk$status


#Estimation of cumulative hazards:

cumh1 = survfit(Surv(tt,dd==1),type="fh2")
cumh2 = survfit(Surv(tt,dd==2),type="fh2")
par(mfrow=c(1,2))
plot(cumh1,fun="cumhaz",mark.time=F,main="Cum.haz. prostate")
plot(cumh2,fun="cumhaz",mark.time=F,main="Cum.haz. other")


#Estimation of cumulative incidence functions:

library(mstate)
library(ggplot2)

cif = Cuminc(time=tt,status=dd)
plot(x = cif,use.ggplot = TRUE,conf.type = "log",lty = 1:2,conf.int = 0.95)