Análise de sobrevivência
Juan Pablo Edwards
25/10/2014
Example 1: lung
data(melanom); attach(melanom)
melanom=melanom[ order(melanom[,1],melanom[,4]), ]
msurv <- with(melanom, Surv(days, status==1))
mean(msurv)## [1] 1076.539mean(msurv[,1])## [1] 2152.8summary(msurv)## time status
## Min. : 10 Min. :0.000
## 1st Qu.:1525 1st Qu.:0.000
## Median :2005 Median :0.000
## Mean :2153 Mean :0.278
## 3rd Qu.:3042 3rd Qu.:1.000
## Max. :5565 Max. :1.000(fit.bysex <- survfit(Surv(days, status == 1) ~ sex, data=melanom))## Call: survfit(formula = Surv(days, status == 1) ~ sex, data = melanom)
##
## records n.max n.start events median 0.95LCL 0.95UCL
## sex=1 126 126 126 28 NA NA NA
## sex=2 79 79 79 29 NA 2388 NAplot(fit.bysex,conf.int=TRUE, col=c("black","grey"), lty=1:2)
pl2= ggsurv(fit.bysex)My example for phytopathology
status: 0 sadio ; 1 sintoma ; 2 esporul ; 3 morto por MA ; 4 morto por outro
df <- expand.grid(isol=1:20, placas=1:5, folhas=1:3)
df <- df[order(df$isol,df$placas),]; head(df)## isol placas folhas
## 1 1 1 1
## 101 1 1 2
## 201 1 1 3
## 21 1 2 1
## 121 1 2 2
## 221 1 2 3df$dias <- sample(c(1:14), replace=TRUE, size=nrow(df))
df$status <-sample(c(0,1), replace=TRUE, size=nrow(df))
attach(df); df[1:15,]## The following object is masked from melanom:
##
## status## isol placas folhas dias status
## 1 1 1 1 13 0
## 101 1 1 2 5 1
## 201 1 1 3 4 1
## 21 1 2 1 4 1
## 121 1 2 2 5 0
## 221 1 2 3 5 1
## 41 1 3 1 2 0
## 141 1 3 2 3 1
## 241 1 3 3 5 1
## 61 1 4 1 7 1
## 161 1 4 2 5 0
## 261 1 4 3 7 1
## 81 1 5 1 13 1
## 181 1 5 2 5 1
## 281 1 5 3 3 1Apenas o isolado 1 Jeito tradicional
isol1=df[1:15,] #isol1=edit(df[1:15,])
boxplot(isol1$dias); points(mean(isol1$dias), col=2);abline(h=mean(isol1$dias), lty=3)
stripchart(dias ~ isol, data = isol1,
vertical = TRUE, method = "jitter",
pch = 21, col = "Grey70", bg = "Grey90",
add = TRUE) 
Com analise de sobreviencia
(msurv <- with(isol1, Surv(dias, status==1)))## [1] 13+ 5 4 4 5+ 5 2+ 3 5 7 5+ 7 13 5 3isol1$surv=msurv
(fit.isol1 <- survfit(Surv(dias, status == 1) ~ isol, data=isol1))## Call: survfit(formula = Surv(dias, status == 1) ~ isol, data = isol1)
##
## records n.max n.start events median 0.95LCL 0.95UCL
## 15 15 15 11 5 5 NAplot(fit.isol1); abline(v=12, h=0.5, col=2, lty=3)
como seria do jeito tradicional
trad = round((isol_mean <- aggregate(dias ~ isol*placas, data = df, mean)),1)
(trad <- trad[order(trad$isol),])## isol placas dias
## 1 1 1 7.3
## 21 1 2 4.7
## 41 1 3 3.3
## 61 1 4 6.3
## 81 1 5 7.0
## 2 2 1 6.0
## 22 2 2 9.0
## 42 2 3 9.7
## 62 2 4 10.3
## 82 2 5 5.3
## 3 3 1 6.3
## 23 3 2 9.0
## 43 3 3 5.7
## 63 3 4 8.7
## 83 3 5 3.0
## 4 4 1 6.7
## 24 4 2 9.0
## 44 4 3 11.0
## 64 4 4 6.0
## 84 4 5 6.7
## 5 5 1 4.7
## 25 5 2 6.0
## 45 5 3 6.3
## 65 5 4 9.3
## 85 5 5 5.0
## 6 6 1 7.0
## 26 6 2 9.3
## 46 6 3 7.7
## 66 6 4 6.7
## 86 6 5 4.7
## 7 7 1 8.3
## 27 7 2 10.0
## 47 7 3 8.0
## 67 7 4 9.0
## 87 7 5 4.3
## 8 8 1 5.3
## 28 8 2 12.0
## 48 8 3 11.7
## 68 8 4 8.0
## 88 8 5 9.3
## 9 9 1 8.0
## 29 9 2 9.3
## 49 9 3 9.7
## 69 9 4 3.7
## 89 9 5 7.0
## 10 10 1 6.7
## 30 10 2 8.7
## 50 10 3 6.0
## 70 10 4 7.3
## 90 10 5 10.0
## 11 11 1 7.0
## 31 11 2 6.0
## 51 11 3 3.3
## 71 11 4 8.0
## 91 11 5 7.7
## 12 12 1 9.0
## 32 12 2 4.7
## 52 12 3 7.0
## 72 12 4 8.0
## 92 12 5 3.7
## 13 13 1 8.7
## 33 13 2 7.3
## 53 13 3 8.3
## 73 13 4 9.7
## 93 13 5 11.7
## 14 14 1 7.3
## 34 14 2 7.3
## 54 14 3 12.7
## 74 14 4 5.7
## 94 14 5 7.7
## 15 15 1 4.7
## 35 15 2 13.3
## 55 15 3 11.0
## 75 15 4 12.7
## 95 15 5 6.0
## 16 16 1 8.0
## 36 16 2 10.0
## 56 16 3 4.3
## 76 16 4 9.3
## 96 16 5 6.7
## 17 17 1 5.3
## 37 17 2 10.3
## 57 17 3 6.3
## 77 17 4 9.0
## 97 17 5 9.7
## 18 18 1 11.0
## 38 18 2 10.7
## 58 18 3 3.3
## 78 18 4 10.3
## 98 18 5 7.7
## 19 19 1 6.7
## 39 19 2 7.7
## 59 19 3 6.7
## 79 19 4 7.3
## 99 19 5 6.0
## 20 20 1 6.3
## 40 20 2 9.7
## 60 20 3 6.3
## 80 20 4 9.3
## 100 20 5 6.7round((isol_mean <- aggregate(dias ~ isol, data = df, mean)),2)## isol dias
## 1 1 5.73
## 2 2 8.07
## 3 3 6.53
## 4 4 7.87
## 5 5 6.27
## 6 6 7.07
## 7 7 7.93
## 8 8 9.27
## 9 9 7.53
## 10 10 7.73
## 11 11 6.40
## 12 12 6.47
## 13 13 9.13
## 14 14 8.13
## 15 15 9.53
## 16 16 7.67
## 17 17 8.13
## 18 18 8.60
## 19 19 6.87
## 20 20 7.67E com analise de sobrevivencia?
(fit.isol <- survfit(Surv(dias, status == 1) ~ isol, data=df))## Call: survfit(formula = Surv(dias, status == 1) ~ isol, data = df)
##
## records n.max n.start events median 0.95LCL 0.95UCL
## isol=1 15 15 15 11 5 5 NA
## isol=2 15 15 15 6 11 7 NA
## isol=3 15 15 15 4 NA 12 NA
## isol=4 15 15 15 5 14 NA NA
## isol=5 15 15 15 6 8 7 NA
## isol=6 15 15 15 8 7 5 NA
## isol=7 15 15 15 9 11 9 NA
## isol=8 15 15 15 4 13 12 NA
## isol=9 15 15 15 11 9 6 NA
## isol=10 15 15 15 10 11 10 NA
## isol=11 15 15 15 10 7 5 NA
## isol=12 15 15 15 7 10 8 NA
## isol=13 15 15 15 6 14 10 NA
## isol=14 15 15 15 8 11 9 NA
## isol=15 15 15 15 7 14 11 NA
## isol=16 15 15 15 8 12 6 NA
## isol=17 15 15 15 7 12 9 NA
## isol=18 15 15 15 5 13 13 NA
## isol=19 15 15 15 8 9 6 NA
## isol=20 15 15 15 9 11 7 NA(test.isol <- survdiff(Surv(dias, status == 1) ~ isol, data=df))## Call:
## survdiff(formula = Surv(dias, status == 1) ~ isol, data = df)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## isol=1 15 11 4.61 8.869071 9.99306
## isol=2 15 6 6.85 0.106125 0.12070
## isol=3 15 4 5.63 0.473824 0.53592
## isol=4 15 5 8.46 1.415030 1.83070
## isol=5 15 6 6.29 0.013605 0.01797
## isol=6 15 8 7.25 0.077095 0.10015
## isol=7 15 9 8.26 0.066109 0.08534
## isol=8 15 4 10.20 3.768955 4.86690
## isol=9 15 11 7.34 1.821085 2.22442
## isol=10 15 10 7.95 0.526163 0.64792
## isol=11 15 10 4.56 6.478453 7.15514
## isol=12 15 7 5.10 0.704414 0.78687
## isol=13 15 6 8.86 0.922084 1.12478
## isol=14 15 8 8.31 0.011218 0.01446
## isol=15 15 7 10.93 1.414592 1.84568
## isol=16 15 8 6.63 0.283727 0.32292
## isol=17 15 7 7.59 0.045697 0.05533
## isol=18 15 5 9.58 2.188442 2.83530
## isol=19 15 8 5.65 0.980516 1.10381
## isol=20 15 9 8.94 0.000446 0.00061
##
## Chisq= 34.5 on 19 degrees of freedom, p= 0.0159Fit a stratified model
cox.model= coxph(Surv(dias, status)~ isol, df)
summary(cox.model)## Call:
## coxph(formula = Surv(dias, status) ~ isol, data = df)
##
## n= 300, number of events= 149
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isol -0.009131 0.990910 0.014535 -0.628 0.53
##
## exp(coef) exp(-coef) lower .95 upper .95
## isol 0.9909 1.009 0.9631 1.02
##
## Concordance= 0.539 (se = 0.029 )
## Rsquare= 0.001 (max possible= 0.991 )
## Likelihood ratio test= 0.39 on 1 df, p=0.5298
## Wald test = 0.39 on 1 df, p=0.5299
## Score (logrank) test = 0.39 on 1 df, p=0.5298plot(fit.isol,conf.int=TRUE, col=c(1:20))
pl2= ggsurv(fit.isol)+
scale_x_continuous( breaks=seq(0,15,by=1)) +
scale_y_continuous(breaks = seq(0, 1, by = .25))
pl2
med.surv <- data.frame(time = c(4,4,5,5,9,9), quant = c(.5,0,.5,0,.5,0),
iso = c('1', '1', '2', '2', '3', '3')
)
pl2 + geom_line(data = med.surv, aes(time, quant, group = iso),
col = 'darkblue', linetype = 3) 
msurv <- with(df, Surv(dias, status==1))
msurv[1:15]## [1] 13+ 5 4 4 5+ 5 2+ 3 5 7 5+ 7 13 5 3mean(msurv)## [1] 4.063333mean(msurv[,2])## [1] 0.4966667summary(msurv)## time status
## Min. : 1.00 Min. :0.0000
## 1st Qu.: 4.00 1st Qu.:0.0000
## Median : 7.00 Median :0.0000
## Mean : 7.63 Mean :0.4967
## 3rd Qu.:11.25 3rd Qu.:1.0000
## Max. :14.00 Max. :1.0000