R-tutorial: survival analysis
Table of Contents
- Loading data for tutorial
- Estimation of the median follow-up time
- Kaplan-Meier
- Cox regression
- Competing risks
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)
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
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
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")
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)
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")
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"))
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)
Stratified Aalen-Johansen
library(riskRegression) data(Melanoma) aj2 <- prodlim(Hist(time,status)~sex,data=Melanoma) plot(aj2,cause=1)
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