R-tutorials Datasets Courses 

R-tutorial: survival analysis

Table of Contents

Open R-markdown version of this file

This tutorial shows some basic tools for survival analysis using R. In particular, how to obtain the Kaplan-Meier graph and how to fit a univariate and a multiple Cox regression model. It is also shown how to export the results in a publishable table format.

For the purpose of illustration we use the German breast cancer data which are for example shipped with the R-package pec:

Loading data for tutorial

We load the data from the 'pec' package:

library(data.table)
library(pec)
data(GBSG2,package="pec")
setDT(GBSG2)
GBSG2
     horTh age menostat tsize tgrade pnodes progrec estrec time cens
  1:    no  70     Post    21     II      3      48     66 1814    1
  2:   yes  56     Post    12     II      7      61     77 2018    1
  3:   yes  58     Post    35     II      9      52    271  712    1
  4:   yes  59     Post    17     II      4      60     29 1807    1
  5:    no  73     Post    35     II      1      26     65  772    1
 ---                                                                
682:    no  49      Pre    30    III      3       1     84  721    0
683:   yes  53     Post    25    III     17       0      0  186    0
684:    no  51      Pre    25    III      5      43      0  769    1
685:    no  52     Post    23     II      3      15     34  727    1
686:    no  55     Post    23     II      9     116     15 1701    1

Then, we make sure that tumor grade is a factor:

GBSG2$tgrade <- as.factor(as.character(GBSG2$tgrade))
GBSG2$Age <- cut(GBSG2$age,c(0,40,60,81),labels=c("21-40","41-60","61-80"))

and we categorize age into 4 categories:

GBSG2$Age <- cut(GBSG2$age,c(-Inf,40,50,65,Inf),labels=c("<40","40-50","51-65",">65"))

We also load the R-packages survival and prodlim and Publish

library(prodlim)
library(survival)
library(Publish)

Estimation of the median follow-up time

When reporting the length of follow-up usually the aim is to describe how long the study was able to observe the enrolled subjects. Naive estimates using e.g., the observed times of all subjects tend to underestimate follow-up because of early events. Schemper and Smith (Control Clin Trials, 1996,17:343–346) suggested to estimate potential follow-up based on the Kaplan-Meier for the censored times. An approximation of this so-called reverse Kaplan-Meier can be obtained by reversing the event indicator so that the outcome of interest becomes being censored. However, when there are ties in the data in the sense that some subjects have exactly the same event times then this would not be exactly correct. Here is how to calculate the reverse Kaplan-Meier estimate which is always exactly correct and how to obtain median follow-up with 95% confidence limits and inter quartile range (IQR):

quantile(prodlim(Hist(time,cens)~1,data=GBSG2,reverse=TRUE))
Quantiles of the potential follow up time distribution based on the Kaplan-Meier method 
applied to the censored times reversing the roles of event status and censored.

Table of quantiles and corresponding confidence limits:
     q quantile lower upper
1 0.00       NA    NA    NA
2 0.25     1976  1869  2027
3 0.50     1645  1570  1714
4 0.75     1100   974  1195
5 1.00        8     8    17
Median time (IQR):1645.00 (1100.00;1976.00)

The median potential follow-up time of the GBSG2 study was 1645 days (IQR: [1100 days;1714 days]). This means that 50% of the patients would have been observed for at least 1645 days had there been no events.

Kaplan-Meier

Kaplan-Meier graph

km0 <- prodlim(Hist(time,cens)~1,data=GBSG2)
plot(km0)

km-gbsg2.png

after some fine tuning:

km0 <- prodlim(Hist(time,cens)~1,data=GBSG2)
par(mar=c(7,7,5,5),                          # margin of figure
    mgp=c(4,1,0))                            # move axis label away from figure
plot(km0,
     xlab="Years",                           # label for x-axis
     ylab="Absolute risk of death",          # label for y-axis
     type="cuminc",                          # increasing risks = 1-survival instead of decreasing survival
     axis1.at=seq(0,2900,365.25),            # time grid for x-axis
     axis1.labels=0:7,                       # time labels for x-axis
     axis2.las=2,                            # rotate labels of y-axis
     atrisk.dist=1,                          # adjust numbers below the figure
     atrisk.labels="Number of \npatients: ") # labels for numbers below figure

tuned-km-gbsg2.png

Stratified Kaplan-Meier

km1 <- prodlim(Hist(time,cens)~tgrade,data=GBSG2)
par(mar=c(7,7,5,5), mgp=c(3,1,0))
plot(km1,
     atrisk.labels=paste("Tumor grade: ",c("I","II","III"),": "),
     atrisk.title="",
     xlab="Years",  # label for x-axis
     axis1.at=seq(0,2900,365.25), # time grid for x-axis
     axis1.labels=0:7, # time labels for x-axis
     legend.x="bottomleft", # positition of legend
     legend.cex=0.8, # font size of legend
     legend.title="Tumor Grade\n", # 
     logrank=TRUE) # show log-rank p-value

km-strata-gbsg2.png

Log-rank test

survdiff(Surv(time,cens)~horTh,data=GBSG2)
Call:
survdiff(formula = Surv(time, cens) ~ horTh, data = GBSG2)

            N Observed Expected (O-E)^2/E (O-E)^2/V
horTh=no  440      205      180      3.37      8.56
horTh=yes 246       94      119      5.12      8.56

 Chisq= 8.6  on 1 degrees of freedom, p= 0.003

Median survival

overall Kaplan-Meier

km0 <- prodlim(Hist(time,cens)~1,data=GBSG2)
quantile(km0)
Quantiles of the event time distribution based on the   method.

Table of quantiles and corresponding confidence limits:
     q quantile lower upper
1 0.00       NA    NA    NA
2 0.25       NA  2456    NA
3 0.50     1807  1528  2018
4 0.75      727   612   801
5 1.00       72    72   120
Median time (IQR):1807.00 (727.00;--)

stratified Kaplan-Meier

km1 <- prodlim(Hist(time,cens)~tgrade,data=GBSG2)
quantile(km1)

Quantiles of the event time distribution based on the method.

Table of quantiles and corresponding confidence limits:

tgrade=I

q quantile lower upper 1 0.00 NA NA NA 2 0.25 NA NA NA 3 0.50 NA 1990 NA 4 0.75 1459 991 NA 5 1.00 476 476 662 Median time (IQR):– (1459.00;–)

Table of quantiles and corresponding confidence limits:

tgrade=II

q quantile lower upper 1 0.00 NA NA NA 2 0.25 NA 2456 NA 3 0.50 1730 1481 2018 4 0.75 745 646 842 5 1.00 72 72 171 Median time (IQR):1730.00 (745.00;–)

Table of quantiles and corresponding confidence limits:

tgrade=III

q quantile lower upper 1 0.00 NA NA NA 2 0.25 NA 2034 NA 3 0.50 1337 956 2034 4 0.75 476 403 571 5 1.00 98 98 177 Median time (IQR):1337.00 (476.00;–)

Quantiles of event time distribution

tgrade=I

Median time (lower.95; upper.95): – (1990;–) IQR (time): (1459;–)

tgrade=II

Median time (lower.95; upper.95): 1730 (1481;2018) IQR (time): (745;–)

tgrade=III

Median time (lower.95; upper.95): 1337 (956;2034) IQR (time): (476;–)

Tabulate results

Overall Kaplan-Meier

km0 <- prodlim(Hist(time,cens)~1,data=GBSG2)
publish(km0,times=seq(0,2900,365.25),org=TRUE)
Interval No. at risk No. of events No. lost to follow-up Survival probability CI.95
0.0– 0.0 686 0 0 100.0 0.0–100.0
0.0– 365.2 686 56 28 91.6 89.4– 93.7
365.2– 730.5 602 109 35 74.6 71.3– 78.0
730.5–1095.8 458 59 68 64.3 60.5– 68.1
1095.8–1461.0 331 39 64 55.9 51.8– 60.0
1461.0–1826.2 228 22 85 49.2 44.7– 53.7
1826.2–2191.5 121 11 74 41.9 36.3– 47.4
2191.5–2556.8 36 3 30 34.3 24.7– 43.8

Stratified Kaplan-Meier

km1 <- prodlim(Hist(time,cens)~tgrade,data=GBSG2)
publish(km1,times=seq(0,2900,365.25),org=TRUE)

tgrade=I

Interval No. at risk No. of events No. lost to follow-up Survival probability CI.95
0.0– 0.0 81 0 0 100.0 0.0–100.0
0.0– 365.2 81 0 4 100.0 0.0–100.0
365.2– 730.5 77 5 7 93.1 87.2– 98.9
730.5–1095.8 65 6 8 83.8 74.9– 92.6
1095.8–1461.0 51 5 10 74.4 63.4– 85.5
1461.0–1826.2 36 0 15 74.4 63.4– 85.5
1826.2–2191.5 21 2 14 62.0 43.8– 80.2
2191.5–2556.8 5 0 5 NA NA– NA

tgrade=II

Interval No. at risk No. of events No. lost to follow-up Survival probability CI.95
0.0– 0.0 444 0 0 100.0 0.0–100.0
0.0– 365.2 444 33 19 92.3 89.8– 94.8
365.2– 730.5 392 69 19 75.8 71.7– 79.9
730.5–1095.8 304 44 42 64.1 59.4– 68.8
1095.8–1461.0 218 26 37 55.7 50.6– 60.8
1461.0–1826.2 155 20 54 47.0 41.4– 52.6
1826.2–2191.5 81 7 46 40.5 33.9– 47.1
2191.5–2556.8 28 3 22 30.7 18.9– 42.5

tgrade=III

Interval No. at risk No. of events No. lost to follow-up Survival probability CI.95
0.0– 0.0 161 0 0 100.0 0.0–100.0
0.0– 365.2 161 23 5 85.3 79.7– 90.8
365.2– 730.5 133 35 9 62.5 54.9– 70.2
730.5–1095.8 89 9 18 55.3 47.2– 63.4
1095.8–1461.0 62 8 17 47.2 38.6– 55.9
1461.0–1826.2 37 2 16 43.5 34.0– 53.0
1826.2–2191.5 19 2 14 35.4 22.6– 48.2
2191.5–2556.8 3 0 3 NA NA– NA

Cox regression

Univariate Cox regression

cox1 <- coxph(Surv(time,cens)~tgrade,data=GBSG2)
summary(cox1)
Call:
coxph(formula = Surv(time, cens) ~ tgrade, data = GBSG2)

  n= 686, number of events= 299 

            coef exp(coef) se(coef)     z  Pr(>|z|)    
tgradeII  0.8718    2.3912   0.2460 3.543  0.000395 ***
tgradeIII 1.1537    3.1701   0.2615 4.411 0.0000103 ***
---
codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

          exp(coef) exp(-coef) lower .95 upper .95
tgradeII      2.391     0.4182     1.476     3.873
tgradeIII     3.170     0.3155     1.899     5.293

Concordance= 0.577  (se = 0.015 )
Rsquare= 0.035   (max possible= 0.995 )
Likelihood ratio test= 24.23  on 2 df,   p=0.000005
Wald test            = 19.75  on 2 df,   p=0.00005
Score (logrank) test = 21.1  on 2 df,   p=0.00003

Comparison of Kaplan-Meier and Cox regression

library(pec)
km.grade <- prodlim(Surv(time,cens)~tgrade,data=GBSG2)
cox.grade <- coxph(Surv(time,cens)~tgrade,data=GBSG2,x=TRUE,y=TRUE)
newdata <- data.frame(tgrade=c("I","II","III"))
## first show Kaplan-Meier without confidence limits
plot(km.grade, lty=1, lwd=3,
      col=c("darkgreen","darkorange","red"), confint=FALSE)
## now add survival estimates based on Cox regression 
plotPredictSurvProb(cox.grade, lty=2,
                    col=c("darkgreen","darkorange","red"),
                    add=TRUE, sort(unique(GBSG2$time)),
                    newdata=newdata)
mtext("Comparison of univariate  Cox regression and stratified Kaplan-Meier")

kaplan-meier-vs-cox-survival-curves.png

Multiple Cox regression

cox2 <- coxph(Surv(time,cens)~tgrade+age+tsize+pnodes+progrec+estrec,data=GBSG2)
summary(cox2)
Call:
coxph(formula = Surv(time, cens) ~ tgrade + age + tsize + pnodes + 
    progrec + estrec, data = GBSG2)

  n= 686, number of events= 299 

                coef  exp(coef)   se(coef)      z        Pr(>|z|)    
tgradeII   0.6459107  1.9077236  0.2490306  2.594         0.00949 ** 
tgradeIII  0.8231461  2.2776542  0.2676730  3.075         0.00210 ** 
age       -0.0028977  0.9971065  0.0061744 -0.469         0.63885    
tsize      0.0071504  1.0071760  0.0038631  1.851         0.06417 .  
pnodes     0.0496877  1.0509429  0.0075636  6.569 0.0000000000505 ***
progrec   -0.0022732  0.9977294  0.0005736 -3.963 0.0000739478425 ***
estrec     0.0001392  1.0001392  0.0004428  0.314         0.75322    
---
codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

          exp(coef) exp(-coef) lower .95 upper .95
tgradeII     1.9077     0.5242    1.1710    3.1081
tgradeIII    2.2777     0.4390    1.3479    3.8489
age          0.9971     1.0029    0.9851    1.0092
tsize        1.0072     0.9929    0.9996    1.0148
pnodes       1.0509     0.9515    1.0355    1.0666
progrec      0.9977     1.0023    0.9966    0.9989
estrec       1.0001     0.9999    0.9993    1.0010

Concordance= 0.684  (se = 0.018 )
Rsquare= 0.131   (max possible= 0.995 )
Likelihood ratio test= 95.93  on 7 df,   p=<2e-16
Wald test            = 104.5  on 7 df,   p=<2e-16
Score (logrank) test = 111.2  on 7 df,   p=<2e-16

Plotting individualized survival probabilities

The following code extracts the survival probabilities for specific combinations of the risk factors: for all three categories of tumor grade and fixed age, tumor size and number of positive lymph nodes.

library(pec)
data(GBSG2)
cox2 <- coxph(Surv(time,cens)~tgrade+age+tsize+pnodes,data=GBSG2,x=TRUE,y=TRUE)
newdata <- data.frame(tgrade=c("I","II","III"),age=50,tsize=30,pnodes=8)
plotPredictSurvProb(cox2,
                    sort(unique(GBSG2$time)),
                    newdata=newdata,
                    col=c("darkgreen","darkorange","red"),
                    legend.title="Tumor grade",
                    legend.legend=c("I","II","III"))
mtext("Individualized survival curves from multiple Cox regression
age=50, tumor size= 30, no. positive lymph nodes=8",line=1.5,cex=1.3)

multiple-cox-survival-curves.jpg

Stratified

Allowing different baseline hazards between subgroups defined by menopause status:

cox3 <- coxph(Surv(time,cens)~tgrade+age+strata(menostat)+tsize+pnodes+progrec+estrec,data=GBSG2,x=TRUE)
summary(cox3)
Call:
coxph(formula = Surv(time, cens) ~ tgrade + age + strata(menostat) + 
    tsize + pnodes + progrec + estrec, data = GBSG2, x = TRUE)

  n= 686, number of events= 299 

                coef  exp(coef)   se(coef)      z        Pr(>|z|)    
tgradeII   0.6206848  1.8602015  0.2492349  2.490        0.012761 *  
tgradeIII  0.8050319  2.2367678  0.2676438  3.008        0.002631 ** 
age       -0.0102152  0.9898368  0.0092931 -1.099        0.271670    
tsize      0.0072755  1.0073020  0.0038726  1.879        0.060287 .  
pnodes     0.0502917  1.0515778  0.0076171  6.602 0.0000000000404 ***
progrec   -0.0022259  0.9977766  0.0005729 -3.885        0.000102 ***
estrec     0.0001704  1.0001704  0.0004505  0.378        0.705218    
---
codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

          exp(coef) exp(-coef) lower .95 upper .95
tgradeII     1.8602     0.5376    1.1413    3.0319
tgradeIII    2.2368     0.4471    1.3237    3.7795
age          0.9898     1.0103    0.9720    1.0080
tsize        1.0073     0.9928    0.9997    1.0150
pnodes       1.0516     0.9510    1.0360    1.0674
progrec      0.9978     1.0022    0.9967    0.9989
estrec       1.0002     0.9998    0.9993    1.0011

Concordance= 0.687  (se = 0.025 )
Rsquare= 0.131   (max possible= 0.99 )
Likelihood ratio test= 96.29  on 7 df,   p=<2e-16
Wald test            = 105.3  on 7 df,   p=<2e-16
Score (logrank) test = 111.7  on 7 df,   p=<2e-16

Interaction (effect modification).

Testing if tumor grade effect depends on age:

cox4 <- coxph(Surv(time,cens)~tgrade+Age+strata(menostat)+tsize+pnodes+progrec+estrec,data=GBSG2)
cox5 <- coxph(Surv(time,cens)~tgrade*Age+strata(menostat)+tsize+pnodes+progrec+estrec,data=GBSG2)
anova(cox4,cox5,test="chisq")
Analysis of Deviance Table
 Cox model: response is  Surv(time, cens)
 Model 1: ~ tgrade + Age + strata(menostat) + tsize + pnodes + progrec + estrec
 Model 2: ~ tgrade * Age + strata(menostat) + tsize + pnodes + progrec + estrec
   loglik  Chisq Df P(>|Chi|)
1 -1535.1                    
2 -1531.0 8.2495  6    0.2204

Possible ways to conclude:

  • The tumor grade effect on the hazard rate is not significantly different between age groups (p=0.47).
  • Within the limitations of this study we were not able to identify a significant difference of the tumor grade effect on the the hazard rate between age groups (p=0.47).

Tabulate results

Multiple Cox regression

Apply publish function to the Cox fit:

cox3 <- coxph(Surv(time,cens)~tgrade+age+strata(menostat)+tsize+pnodes+progrec+estrec,data=GBSG2)
publish(cox3,org=TRUE)
Variable Units HazardRatio CI.95 p-value
tgrade I Ref    
  II 1.86 [1.14;3.03] 0.01276
  III 2.24 [1.32;3.78] 0.00263
age   0.99 [0.97;1.01] 0.27167
tsize   1.01 [1.00;1.01] 0.06029
pnodes   1.05 [1.04;1.07] < 0.001
progrec   1.00 [1.00;1.00] < 0.001
estrec   1.00 [1.00;1.00] 0.70522

With interaction term

cox5 <- coxph(Surv(time,cens)~tgrade * menostat+tsize+age+pnodes+progrec+estrec,data=GBSG2)
publish(cox5,org=TRUE)
Variable Units HazardRatio CI.95 p-value
tsize   1.01 [1.00;1.01] 0.0691
age   0.99 [0.97;1.01] 0.2426
pnodes   1.05 [1.04;1.07] <0.001
progrec   1.00 [1.00;1.00] <0.001
estrec   1.00 [1.00;1.00] 0.7903
tgrade(I): menostat(Pre vs Post)   0.53 [0.19;1.45] 0.2125
tgrade(II): menostat(Pre vs Post)   0.84 [0.57;1.24] 0.3771
tgrade(III): menostat(Pre vs Post)   0.81 [0.47;1.40] 0.4435
menostat(Post): tgrade(II vs I)   1.58 [0.87;2.88] 0.1330
menostat(Post): tgrade(III vs I)   1.92 [1.00;3.67] 0.0483
menostat(Pre): tgrade(II vs I)   2.53 [1.10;5.83] 0.0297
menostat(Pre): tgrade(III vs I)   2.95 [1.23;7.10] 0.0156

Note: the main effects do usually not have an interpretation in the presence of interaction (effect modification).

Forest plots: showing hazard ratios with confidence intervals

par(cex=4)
cox2 <- coxph(Surv(time,cens)~tgrade+age+tsize+pnodes,data=GBSG2)
rt2 <- do.call("rbind",regressionTable(cox2))
plotConfidence(x=rt2$HazardRatio,
               lower=rt2$Lower,
               upper=rt2$Upper,
               xlim=c(0,5),
               factor.reference.pos=1,
               labels=c("Tumor grade I",
                   "                     II",
                   "                     III",
                   "Age (years)",
                   "Tumor size (mm)",
                   "No. positive\nlymph nodes"),
               title.labels="Factor",
               cex=2,
               xlab.cex=1.3,
               xlab="Hazard ratio")

forest-cox2.jpg

Competing risks

Visualization of the model

library(riskRegression)
data(Melanoma)
levels(Melanoma$event) <- gsub("\\.","\n\n",levels(Melanoma$event))
with(Melanoma,plot(Hist(time,event,cens.code="censored"),arrowLabelStyle="count"))

cr-model-melanoma.png

Aalen-Johansen graph

The graph shows the cumulative incidence of cancer related death since surgery.

library(riskRegression)
data(Melanoma)
aj <- prodlim(Hist(time,status)~1,data=Melanoma)
plot(aj,cause=1)

aalen-johansen-graph.png

Stratified Aalen-Johansen

library(riskRegression)
data(Melanoma)
aj2 <- prodlim(Hist(time,status)~sex,data=Melanoma)
plot(aj2,cause=1)

aalen-johansen-graph-sex.png

Gray's test

library(cmprsk)
with(Melanoma,cuminc(ftime=time,fstatus=status,group=sex))$Tests
       stat        pv df
1 5.8140209 0.0158989  1
2 0.8543656 0.3553203  1

  • The Gray test shows a significantly higher risk of cancer related death for males compared to females (p-value= 0.0159).
  • The Gray test shows no non-significant difference in risks of other cause mortality between males and females (p-value= 0.3553).

Cox regression

The function fits Cox regression models, one for each competing risk.

Same formula for both competing risks

library(riskRegression)
CSC(Hist(time,status)~sex+age+invasion+logthick+ulcer,data=Melanoma)
CSC(formula = Hist(time, status) ~ sex + age + invasion + logthick + 
    ulcer, data = Melanoma)

Right-censored response of a competing.risks model

No.Observations: 205 

Pattern:
         
Cause     event right.censored
  1          57              0
  2          14              0
  unknown     0            134


----------> Cause:  1 

Call:
survival::coxph(formula = survival::Surv(time, status) ~ sex + 
    age + invasion + logthick + ulcer, x = TRUE, y = TRUE)

  n= 205, number of events= 57 

                    coef exp(coef) se(coef)     z Pr(>|z|)   
sexMale         0.416256  1.516273 0.277562 1.500  0.13370   
age             0.010690  1.010747 0.008465 1.263  0.20667   
invasionlevel.1 0.325882  1.385252 0.386854 0.842  0.39957   
invasionlevel.2 0.399668  1.491329 0.536187 0.745  0.45604   
logthick        0.422402  1.525621 0.250260 1.688  0.09144 . 
ulcerpresent    0.951752  2.590245 0.329443 2.889  0.00386 **
---
codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                exp(coef) exp(-coef) lower .95 upper .95
sexMale             1.516     0.6595    0.8801     2.612
age                 1.011     0.9894    0.9941     1.028
invasionlevel.1     1.385     0.7219    0.6490     2.957
invasionlevel.2     1.491     0.6705    0.5214     4.266
logthick            1.526     0.6555    0.9342     2.492
ulcerpresent        2.590     0.3861    1.3581     4.940

Concordance= 0.767  (se = 0.04 )
Rsquare= 0.199   (max possible= 0.937 )
Likelihood ratio test= 45.44  on 6 df,   p=0.00000004
Wald test            = 38.67  on 6 df,   p=0.0000008
Score (logrank) test = 45.79  on 6 df,   p=0.00000003



----------> Cause:  2 

Call:
survival::coxph(formula = survival::Surv(time, status) ~ sex + 
    age + invasion + logthick + ulcer, x = TRUE, y = TRUE)

  n= 205, number of events= 14 

                    coef exp(coef) se(coef)      z Pr(>|z|)    
sexMale          0.19618   1.21675  0.55962  0.351 0.725915    
age              0.09048   1.09469  0.02554  3.543 0.000396 ***
invasionlevel.1 -1.40328   0.24579  0.73843 -1.900 0.057388 .  
invasionlevel.2 -2.08963   0.12373  1.23631 -1.690 0.090987 .  
logthick         0.44436   1.55949  0.39440  1.127 0.259884    
ulcerpresent     0.12655   1.13491  0.65348  0.194 0.846444    
---
codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                exp(coef) exp(-coef) lower .95 upper .95
sexMale            1.2167     0.8219   0.40630     3.644
age                1.0947     0.9135   1.04125     1.151
invasionlevel.1    0.2458     4.0685   0.05781     1.045
invasionlevel.2    0.1237     8.0819   0.01097     1.396
logthick           1.5595     0.6412   0.71989     3.378
ulcerpresent       1.1349     0.8811   0.31529     4.085

Concordance= 0.835  (se = 0.085 )
Rsquare= 0.097   (max possible= 0.481 )
Likelihood ratio test= 20.88  on 6 df,   p=0.002
Wald test            = 16.94  on 6 df,   p=0.01
Score (logrank) test = 18.36  on 6 df,   p=0.005

Different formula for the two competing risks

library(riskRegression)
CSC(list(Hist(time,status)~sex+age+invasion+logthick+ulcer,Hist(time,status)~sex+age),data=Melanoma)
CSC(formula = list(Hist(time, status) ~ sex + age + invasion + 
    logthick + ulcer, Hist(time, status) ~ sex + age), data = Melanoma)

Right-censored response of a competing.risks model

No.Observations: 205 

Pattern:
         
Cause     event right.censored
  1          57              0
  2          14              0
  unknown     0            134


----------> Cause:  1 

Call:
survival::coxph(formula = survival::Surv(time, status) ~ sex + 
    age + invasion + logthick + ulcer, x = TRUE, y = TRUE)

  n= 205, number of events= 57 

                    coef exp(coef) se(coef)     z Pr(>|z|)   
sexMale         0.416256  1.516273 0.277562 1.500  0.13370   
age             0.010690  1.010747 0.008465 1.263  0.20667   
invasionlevel.1 0.325882  1.385252 0.386854 0.842  0.39957   
invasionlevel.2 0.399668  1.491329 0.536187 0.745  0.45604   
logthick        0.422402  1.525621 0.250260 1.688  0.09144 . 
ulcerpresent    0.951752  2.590245 0.329443 2.889  0.00386 **
---
codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

                exp(coef) exp(-coef) lower .95 upper .95
sexMale             1.516     0.6595    0.8801     2.612
age                 1.011     0.9894    0.9941     1.028
invasionlevel.1     1.385     0.7219    0.6490     2.957
invasionlevel.2     1.491     0.6705    0.5214     4.266
logthick            1.526     0.6555    0.9342     2.492
ulcerpresent        2.590     0.3861    1.3581     4.940

Concordance= 0.767  (se = 0.04 )
Rsquare= 0.199   (max possible= 0.937 )
Likelihood ratio test= 45.44  on 6 df,   p=0.00000004
Wald test            = 38.67  on 6 df,   p=0.0000008
Score (logrank) test = 45.79  on 6 df,   p=0.00000003



----------> Cause:  2 

Call:
survival::coxph(formula = survival::Surv(time, status) ~ sex + 
    age, x = TRUE, y = TRUE)

  n= 205, number of events= 14 

           coef exp(coef) se(coef)     z Pr(>|z|)    
sexMale 0.36910   1.44643  0.54712 0.675 0.499914    
age     0.07581   1.07876  0.02155 3.518 0.000434 ***
---
codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

        exp(coef) exp(-coef) lower .95 upper .95
sexMale     1.446     0.6914     0.495     4.227
age         1.079     0.9270     1.034     1.125

Concordance= 0.821  (se = 0.085 )
Rsquare= 0.076   (max possible= 0.481 )
Likelihood ratio test= 16.22  on 2 df,   p=0.0003
Wald test            = 13.96  on 2 df,   p=0.0009
Score (logrank) test = 14.47  on 2 df,   p=0.0007
Last update: 23 Nov 2018 by Thomas Alexander Gerds.