Análise de Sobrevivência Lista 1
Mauro Citro Lalucci
2023-05-31
Pacotes utilizados nesta lista
library(survminer)## Warning: package 'survminer' was built under R version 4.2.3## Carregando pacotes exigidos: ggplot2## Warning: package 'ggplot2' was built under R version 4.2.3## Carregando pacotes exigidos: ggpubr## Warning: package 'ggpubr' was built under R version 4.2.3library(survival)##
## Attaching package: 'survival'## The following object is masked from 'package:survminer':
##
## myeloma1 - Exercício 1
Os dados da tabela abaixo referem-se aos tempos de sobrevivência (em dias) de pacientes com câncer submetidos à radioterapia (o símbolo + indica censura).
Tabela: Tempos de pacientes submetidos à radioterapia. \[ \begin{matrix} 7 & 34 & 42 & 63 & 64 & 74+ & 83 & 84\\ 91 & 108 & 112 & 129 & 133 & 133 & 139 & 140\\ 140 & 146 & 149 & 154 & 160 & 165 & 173 & 176\\ 185+ & 218 & 225 & 241 & 248 & 273 & 277 & 279+\\ 297 & 319+ & 405 & 417 & 420 & 440 & 523+ & 583\\ 594 & 1101 & 1116+ & 1146 & 1226+ & 1349+ & 1412+ & 1417 \end{matrix} \]
cens <- as.numeric(c(1,1,1,1,1,0,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,0,1,1,1,1,1,1,0,1,0,1,1,1,1,0,1,1,1,0,1,0,0,0,1))
# 1 = normal; 0 = censura
tempo <- as.numeric(c(7, 34, 42, 63, 64, 74, 83, 84, 91, 108, 112, 129, 133, 133, 139,
140, 140, 146, 149, 154, 160, 160, 165, 173, 176, 185, 218,
225, 241, 248, 273, 277, 279, 297, 319, 405, 417, 420, 440,
523, 583, 594, 1101, 1116, 1146, 1226, 1349, 1412, 1417))Obtenha: o estimador de Kaplan-Meier, de Nelson-Aalen, o tempo mediano, o tempo médio e construa a curva TTT.
1.1 - Estimador de Kaplan-Meier
surv <- Surv(tempo, cens)
ekm <- survfit(surv ~ 1)
tabela_ekm <- summary(ekm)
plot(ekm,conf.int=F, xlab="Tempo (dias)", ylab="S(t) estimada")
title('Kaplan-Meier')
O gráfico acima nos mostra a função de sobrevivência para o grupo
analisado.
1.2 - Estimador de Nelson-Aalen
ena <- survfit(surv ~ 1, type = "fleming-harrington")
summary(ena)## Call: survfit(formula = surv ~ 1, type = "fleming-harrington")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 7 49 1 0.9798 0.0200 0.94138 1.000
## 34 48 1 0.9596 0.0280 0.90629 1.000
## 42 47 1 0.9394 0.0339 0.87523 1.000
## 63 46 1 0.9192 0.0387 0.84633 0.998
## 64 45 1 0.8990 0.0428 0.81885 0.987
## 83 43 1 0.8783 0.0466 0.79165 0.975
## 84 42 1 0.8577 0.0498 0.76533 0.961
## 91 41 1 0.8370 0.0528 0.73974 0.947
## 108 40 1 0.8163 0.0553 0.71475 0.932
## 112 39 1 0.7957 0.0577 0.69028 0.917
## 129 38 1 0.7750 0.0598 0.66628 0.901
## 133 37 2 0.7342 0.0632 0.62024 0.869
## 139 35 1 0.7135 0.0647 0.59734 0.852
## 140 34 2 0.6728 0.0671 0.55328 0.818
## 146 32 1 0.6521 0.0682 0.53126 0.800
## 149 31 1 0.6314 0.0691 0.50951 0.782
## 154 30 1 0.6107 0.0698 0.48804 0.764
## 160 29 2 0.5700 0.0709 0.44671 0.727
## 165 27 1 0.5493 0.0713 0.42593 0.708
## 173 26 1 0.5285 0.0715 0.40540 0.689
## 176 25 1 0.5078 0.0717 0.38511 0.670
## 218 23 1 0.4862 0.0718 0.36403 0.649
## 225 22 1 0.4646 0.0718 0.34322 0.629
## 241 21 1 0.4430 0.0716 0.32269 0.608
## 248 20 1 0.4214 0.0713 0.30244 0.587
## 273 19 1 0.3998 0.0708 0.28247 0.566
## 277 18 1 0.3782 0.0702 0.26279 0.544
## 297 16 1 0.3553 0.0696 0.24197 0.522
## 405 14 1 0.3308 0.0690 0.21979 0.498
## 417 13 1 0.3063 0.0681 0.19811 0.474
## 420 12 1 0.2818 0.0669 0.17695 0.449
## 440 11 1 0.2573 0.0654 0.15634 0.423
## 583 9 1 0.2302 0.0639 0.13367 0.397
## 594 8 1 0.2032 0.0618 0.11192 0.369
## 1101 7 1 0.1761 0.0592 0.09114 0.340
## 1146 5 1 0.1442 0.0564 0.06700 0.310
## 1417 1 1 0.0531 0.0570 0.00647 0.435plot(ena,conf.int=F, xlab="Tempo (dias)", ylab="S(t) estimada")
title('Nelson-Aalen')
Temos a função de sobrevivência dos dados observados estimada utilizando o estimador de Nelson-Aalen
1.3 - Tempo mediano
Kaplan-Meier \[ \frac{218-176}{0,479-0,501}=\frac{Md-176}{0,5-0,501}\rightarrow Md=174,0909\mbox{ dias} \]
Nelson-Aalen \[ \frac{218-176}{0,4862-0,5078}=\frac{Md-176}{0,5-0,5078}\rightarrow Md=191,1667\mbox{ dias} \]
Função para Obter o valor mediano da sobrevivência com mais precisão
calcular_Md <- function(survfit_obj) {
# Obtendo o resumo do objeto survfit
summary_obj <- summary(survfit_obj)
# Extraindo as estimativas de sobrevivencia
survival <- summary_obj$surv
time <- summary_obj$time
tabela_sobrevivencia <- data.frame(Tempo = time,
Sobrevida = survival)
prob_acima <- tail(subset(tabela_sobrevivencia, Sobrevida > 0.5)$Sobrevida, 1)
prob_abaixo <- head(subset(tabela_sobrevivencia, Sobrevida < 0.5)$Sobrevida, 1)
tempo_acima <- tail(subset(tabela_sobrevivencia, Sobrevida > 0.5)$Tempo, 1)
tempo_abaixo <- head(subset(tabela_sobrevivencia, Sobrevida < 0.5)$Tempo, 1)
# Calculando Md
resultado <- (tempo_abaixo - tempo_acima) / (prob_abaixo - prob_acima)
Md <- resultado * (0.5 - prob_acima) + tempo_acima
return(Md)
}paste(calcular_Md(ekm),"Mediana de Kaplan-Meier em dias")## [1] "178.286931818181 Mediana de Kaplan-Meier em dias"paste(calcular_Md(ena),"Mediana de Nelson-Aalen em dias")## [1] "191.183699986637 Mediana de Nelson-Aalen em dias"1.4 - Tempo médio
t<- tempo[cens==1]
tj<-c(0,as.numeric(levels(as.factor(t))))
surv<-c(1,as.numeric(levels(as.factor(ekm$surv))))
surv<-sort(surv, decreasing=T)
k<-length(tj)-1
prod<-matrix(0,k,1)
for(j in 1:k){
prod[j]<-(tj[j+1]-tj[j])*surv[j]
}
TMekm<-sum(prod)
paste(TMekm,"Tempo médio de Kaplan-Meier")## [1] "427.11553620437 Tempo médio de Kaplan-Meier" t<- tempo[cens==1]
tj<-c(0,as.numeric(levels(as.factor(t))))
surv<-c(1,as.numeric(levels(as.factor(ena$surv))))
surv<-sort(surv, decreasing=T)
k<-length(tj)-1
prod<-matrix(0,k,1)
for(j in 1:k){
prod[j]<-(tj[j+1]-tj[j])*surv[j]
}
TMena<-sum(prod)
paste(TMena,"Tempo médio de Nelson-Aalen")## [1] "439.8365558102 Tempo médio de Nelson-Aalen"1.5 - Construir a curva TTT
TTTplot <- function(time,censur=NULL)
{
n <- length(time)
oo <- order(time)
if(is.null(censur)) censur = rep(1,n)
sorttime <- time[oo]
sortcens <- cens[oo]
ttt <- numeric(n)
ttt[1] <- n*sorttime[1]
for(i in 2:n)
{
ttt[i] <- ttt[i-1]+(n-i+1)*(sorttime[i]-sorttime[i-1])
}
ttt <- data.frame(time=sorttime,censur = sortcens,
rank=1:n,"cum total time"=ttt)
colnames(ttt) <- c("time","censur","rank","Cum.Total.Time")#,"i/n","TTT")
TTT <- ttt[ttt$censur==1,]
nfail <- dim(TTT)[1]
TTT <- cbind(TTT[,-3],TTT$Cum.Total.Time/TTT$Cum.Total.Time[nfail],
(1:nfail)/nfail)
colnames(TTT) <- c("time","censur","Cum.Total.Time","TTT","i/n")
W <- sum(TTT$"Cum.Total.Time"[1:(nfail-1)])/TTT$"Cum.Total.Time"[nfail]
Z <- (W-(nfail-1)/2)/sqrt((nfail-1)/12)
print(TTT)
cat("Barlow-Proschan's test\n")
cat(paste("W=",round(W,4)," k=",nfail,"\n",sep=""))
cat(paste("Z=",round(Z,4)," p.value=",round(pnorm(-abs(Z)),4),"\n",sep=""))
plot(TTT$"i/n",TTT$TTT,xlab="",ylab="TTT-plot",xlim=c(0,1),
ylim=c(0,1),pch=20, main="",sub=paste("nº. falhas = ",nfail," nº. censuras = ", n-nfail,sep=""),type="l",lwd=1)
lines(c(0,1),c(0,1),lty=2,lwd=1)
# return(list(TTT=TTT,W=W,Z=Z))
}
TTTplot(tempo,cens)## time censur Cum.Total.Time TTT i/n
## 1 7 1 343 0.01952191 0.025
## 2 34 1 1639 0.09328401 0.050
## 3 42 1 2015 0.11468412 0.075
## 4 63 1 2981 0.16966420 0.100
## 5 64 1 3026 0.17222538 0.125
## 7 83 1 3853 0.21929425 0.150
## 8 84 1 3895 0.22168469 0.175
## 9 91 1 4182 0.23801935 0.200
## 10 108 1 4862 0.27672168 0.225
## 11 112 1 5018 0.28560046 0.250
## 12 129 1 5664 0.32236767 0.275
## 13 133 1 5812 0.33079112 0.300
## 14 133 1 5812 0.33079112 0.325
## 15 139 1 6022 0.34274331 0.350
## 16 140 1 6056 0.34467843 0.375
## 17 140 1 6056 0.34467843 0.400
## 18 146 1 6248 0.35560615 0.425
## 19 149 1 6341 0.36089926 0.450
## 20 154 1 6491 0.36943654 0.475
## 21 160 1 6665 0.37933978 0.500
## 22 160 1 6665 0.37933978 0.525
## 23 165 1 6800 0.38702334 0.550
## 24 173 1 7008 0.39886170 0.575
## 25 176 1 7083 0.40313034 0.600
## 27 218 1 8058 0.45862265 0.625
## 28 225 1 8212 0.46738759 0.650
## 29 241 1 8548 0.48651110 0.675
## 30 248 1 8688 0.49447923 0.700
## 31 273 1 9163 0.52151394 0.725
## 32 277 1 9235 0.52561184 0.750
## 34 297 1 9557 0.54393853 0.775
## 36 405 1 11091 0.63124644 0.800
## 37 417 1 11247 0.64012521 0.825
## 38 420 1 11283 0.64217416 0.850
## 39 440 1 11503 0.65469550 0.875
## 41 583 1 12873 0.73266932 0.900
## 42 594 1 12961 0.73767786 0.925
## 43 1101 1 16510 0.93966989 0.950
## 45 1146 1 16750 0.95332954 0.975
## 49 1417 1 17570 1.00000000 1.000
## Barlow-Proschan's test
## W=16.29 k=40
## Z=-1.7806 p.value=0.0375
2 - Exercício 2
Os dados representados na tabela abaixo representam o tempo (em dias) até a morte de pacientes com câncer de ovário tratados na Mayo Clinic. O símbolo + indica censura.
\[ 1-\mbox{Tumor Grande}\\ \begin{matrix} 28 & 89 & 175 & 195 & 309\\ 377+ & 393+ & 421+ & 447+ & 462\\ 709+ & 744+ & 770+ & 1106+ & 1206+ \end{matrix}\\ 2-\mbox{Tumor Pequeno}\\ \begin{matrix} 34 & 88 & 137 & 199 & 280\\ 291 & 299+ & 300+ & 309 & 351\\ 358 & 369 & 369 & 370 & 375\\ 382 & 392 & 429+ & 451 & 1119+ \end{matrix} \]
tempos <- c(28,89,175,195,309,377,393,421,447,462,709,744,770,1106,1206,34,88,137,199,280,291,299,300,309,351,358,369,369,370,375,382,392,429,451,1119)
# 1 = normal; 0 = censura
cens <- c(1,1,1,1,1,0,0,0,0,1,0,0,0,0,0,1,1,1,1,1,1,0,0,1,1,1,1,1,1,1,1,1,0,1,0)2.1 - Estimador de Kaplan-Meier
# Armazenando os grupos, sendo grupo 1 = 15 obs; grupo 2 = 20 obs
grupos=c(rep(1,15),rep(2,20))
# Montando o grafico
ekm=survfit(Surv(tempos,cens)~grupos)
summary(ekm)## Call: survfit(formula = Surv(tempos, cens) ~ grupos)
##
## grupos=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 28 15 1 0.933 0.0644 0.815 1.000
## 89 14 1 0.867 0.0878 0.711 1.000
## 175 13 1 0.800 0.1033 0.621 1.000
## 195 12 1 0.733 0.1142 0.540 0.995
## 309 11 1 0.667 0.1217 0.466 0.953
## 462 6 1 0.556 0.1434 0.335 0.922
##
## grupos=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 34 20 1 0.9500 0.0487 0.8591 1.000
## 88 19 1 0.9000 0.0671 0.7777 1.000
## 137 18 1 0.8500 0.0798 0.7071 1.000
## 199 17 1 0.8000 0.0894 0.6426 0.996
## 280 16 1 0.7500 0.0968 0.5823 0.966
## 291 15 1 0.7000 0.1025 0.5254 0.933
## 309 12 1 0.6417 0.1093 0.4596 0.896
## 351 11 1 0.5833 0.1139 0.3979 0.855
## 358 10 1 0.5250 0.1165 0.3399 0.811
## 369 9 2 0.4083 0.1162 0.2338 0.713
## 370 7 1 0.3500 0.1133 0.1856 0.660
## 375 6 1 0.2917 0.1084 0.1408 0.604
## 382 5 1 0.2333 0.1012 0.0997 0.546
## 392 4 1 0.1750 0.0912 0.0630 0.486
## 451 2 1 0.0875 0.0769 0.0156 0.489plot(ekm,xlab="Tempo (dias)", ylab="S(t) estimada", lty=c(2,1))
legend(1,0.3,c("Tumor Grande","Tumor Pequeno"),lty=c(2,1),lwd=1,bty="o")
2.2 - Estimador de Nelson-Aalen
ena=survfit(coxph(Surv(tempos,cens)~1,method="breslow"))
summary(ena)## Call: survfit(formula = coxph(Surv(tempos, cens) ~ 1, method = "breslow"))
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 28 35 1 0.972 0.0278 0.919 1.000
## 34 34 1 0.944 0.0387 0.871 1.000
## 88 33 1 0.915 0.0467 0.828 1.000
## 89 32 1 0.887 0.0531 0.789 0.998
## 137 31 1 0.859 0.0584 0.752 0.982
## 175 30 1 0.831 0.0629 0.716 0.964
## 195 29 1 0.803 0.0668 0.682 0.945
## 199 28 1 0.775 0.0701 0.649 0.925
## 280 27 1 0.746 0.0730 0.616 0.904
## 291 26 1 0.718 0.0755 0.585 0.883
## 309 23 2 0.659 0.0802 0.519 0.836
## 351 21 1 0.628 0.0821 0.486 0.811
## 358 20 1 0.597 0.0836 0.454 0.786
## 369 19 2 0.538 0.0852 0.394 0.733
## 370 17 1 0.507 0.0857 0.364 0.706
## 375 16 1 0.476 0.0858 0.334 0.678
## 382 14 1 0.443 0.0860 0.303 0.648
## 392 13 1 0.411 0.0856 0.273 0.618
## 451 8 1 0.362 0.0881 0.225 0.584
## 462 7 1 0.314 0.0886 0.181 0.546plot(ena,conf.int=F, xlab="t", ylab="S(t)", lty=1, lwd=1)
2.3 - Tempo mediano
Função para Obter o valor mediano da sobrevivência com mais precisão
calcular_Md(ekm)## [1] 360.3571calcular_Md(ena)## [1] 371.11982.4 O Tempo Médio
t<- tempos[cens==1]
tj<-c(0,as.numeric(levels(as.factor(t))))
surv<-c(1,as.numeric(levels(as.factor(ekm$surv))))
surv<-sort(surv, decreasing=T)
k<-length(tj)-1
prod<-matrix(0,k,1)
for(j in 1:k){
prod[j]<-(tj[j+1]-tj[j])*surv[j]
}
TMekm<-sum(prod)
paste(TMekm,"Tempo médio de Kaplan-Meier")## [1] "317.769444444444 Tempo médio de Kaplan-Meier" t<- tempos[cens==1]
tj<-c(0,as.numeric(levels(as.factor(t))))
surv<-c(1,as.numeric(levels(as.factor(ena$surv))))
surv<-sort(surv, decreasing=T)
k<-length(tj)-1
prod<-matrix(0,k,1)
for(j in 1:k){
prod[j]<-(tj[j+1]-tj[j])*surv[j]
}
TMena<-sum(prod)
paste(TMena,"Tempo médio de Nelson-Aalen")## [1] "342.331818683236 Tempo médio de Nelson-Aalen"2.5 Comparando as curvas por Log-rank e Peto
Vamos comparar as curvas pelo teste de Log-rank
survdiff(formula=Surv(tempos[1:35],cens[1:35])~grupos[1:35],rho=0)## Call:
## survdiff(formula = Surv(tempos[1:35], cens[1:35]) ~ grupos[1:35],
## rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## grupos[1:35]=1 15 6 11.3 2.51 5.57
## grupos[1:35]=2 20 16 10.7 2.67 5.57
##
## Chisq= 5.6 on 1 degrees of freedom, p= 0.02Vamos comparar as curvas pelo teste de Peto
survdiff(formula = Surv(tempos, cens) ~ grupos, rho = 1)## Call:
## survdiff(formula = Surv(tempos, cens) ~ grupos, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## grupos=1 15 4.66 7.40 1.013 2.74
## grupos=2 20 10.60 7.87 0.953 2.74
##
## Chisq= 2.7 on 1 degrees of freedom, p= 0.1Exercício 3
Dados de tempo até a ocorrência de cio em ovelhas submetidas a três tratamentos. Tabela: Tempo até a ocorrência de cio em ovelhas de três amostras, os tempos de sobrevivência estão em dias.
\[ 1-\mbox{Tratamento A}\\ \begin{matrix} 57 & 60+ & 60+ & 60+ & 48\\ 40 & 49 & 45 & 25 & 60+ & 51\\ \end{matrix}\\ 2-\mbox{Tratamento A}\\ \begin{matrix} 32 & 60+ & 60+ & 22 & 60+\\ 54 & 60+ & 60+ & 60+ & 51\\ \end{matrix}\\ 3-\mbox{Tratamento C}\\ \begin{matrix} 28 & 32 & 60+ & 41 & 60+\\ 54 & 59 & 40 & 32 & 60+ & 60+\\ \end{matrix} \] # 3.1 Kaplan-Meier (estratificado)
tempos <- c(57,60,60,60,48,40,49,45,25,60,51,32,60,60,22,60,54,60,60,60,51,28,32,60,41,60,54,59,40,32,60,60)
# 1 = normal; 0 = censura
cens <- c(1,0,0,0,1,1,1,1,1,0,1,1,0,0,1,0,1,0,0,0,1,1,1,0,1,0,1,1,1,1,0,0)Agora vamos plotar o gráfico com os dados coletados para Kaplan-Meier com estratificação.
# Armazenando os grupos, sendo grupo 1 = 15 obs; grupo 2 = 20 obs
grupos=c(rep(1,11),rep(2,10),rep(3,11))
# Montando o grafico
ekm=survfit(Surv(tempos,cens)~grupos)
summary(ekm)## Call: survfit(formula = Surv(tempos, cens) ~ grupos)
##
## grupos=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 25 11 1 0.909 0.0867 0.754 1.000
## 40 10 1 0.818 0.1163 0.619 1.000
## 45 9 1 0.727 0.1343 0.506 1.000
## 48 8 1 0.636 0.1450 0.407 0.995
## 49 7 1 0.545 0.1501 0.318 0.936
## 51 6 1 0.455 0.1501 0.238 0.868
## 57 5 1 0.364 0.1450 0.166 0.795
##
## grupos=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 22 10 1 0.9 0.0949 0.732 1.000
## 32 9 1 0.8 0.1265 0.587 1.000
## 51 8 1 0.7 0.1449 0.467 1.000
## 54 7 1 0.6 0.1549 0.362 0.995
##
## grupos=3
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 28 11 1 0.909 0.0867 0.754 1.000
## 32 10 2 0.727 0.1343 0.506 1.000
## 40 8 1 0.636 0.1450 0.407 0.995
## 41 7 1 0.545 0.1501 0.318 0.936
## 54 6 1 0.455 0.1501 0.238 0.868
## 59 5 1 0.364 0.1450 0.166 0.795plot(ekm,xlab="Tempo (semanas)", ylab="S(t) estimada", lty=c(2,1))
legend(1,0.3,c("Tratamento A","Tratamento B","Tratamento C"),lty=c(2,1),lwd=1,bty="o")
Nelson-Aalen
Agora vamos fazer com os mesmos dados, o estimador de Nelson-Aalen no R.
ena=survfit(coxph(Surv(tempos,cens)~1,method="breslow"))
summary(ena)## Call: survfit(formula = coxph(Surv(tempos, cens) ~ 1, method = "breslow"))
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 22 32 1 0.969 0.0303 0.912 1.000
## 25 31 1 0.938 0.0421 0.859 1.000
## 28 30 1 0.908 0.0508 0.813 1.000
## 32 29 3 0.818 0.0670 0.697 0.961
## 40 26 2 0.758 0.0745 0.625 0.919
## 41 24 1 0.727 0.0776 0.590 0.896
## 45 23 1 0.696 0.0802 0.555 0.872
## 48 22 1 0.665 0.0824 0.522 0.848
## 49 21 1 0.634 0.0842 0.489 0.823
## 51 20 2 0.574 0.0863 0.427 0.771
## 54 18 2 0.513 0.0871 0.368 0.716
## 57 16 1 0.482 0.0872 0.338 0.688
## 59 15 1 0.451 0.0870 0.309 0.658plot(ena,conf.int=F, xlab="t", ylab="S(t)", lty=1, lwd=1)
Tempo Mediano
Vamos verificar a situação da descritiva dos nossos dados para Kaplan-Meier.
summary(ekm)## Call: survfit(formula = Surv(tempos, cens) ~ grupos)
##
## grupos=1
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 25 11 1 0.909 0.0867 0.754 1.000
## 40 10 1 0.818 0.1163 0.619 1.000
## 45 9 1 0.727 0.1343 0.506 1.000
## 48 8 1 0.636 0.1450 0.407 0.995
## 49 7 1 0.545 0.1501 0.318 0.936
## 51 6 1 0.455 0.1501 0.238 0.868
## 57 5 1 0.364 0.1450 0.166 0.795
##
## grupos=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 22 10 1 0.9 0.0949 0.732 1.000
## 32 9 1 0.8 0.1265 0.587 1.000
## 51 8 1 0.7 0.1449 0.467 1.000
## 54 7 1 0.6 0.1549 0.362 0.995
##
## grupos=3
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 28 11 1 0.909 0.0867 0.754 1.000
## 32 10 2 0.727 0.1343 0.506 1.000
## 40 8 1 0.636 0.1450 0.407 0.995
## 41 7 1 0.545 0.1501 0.318 0.936
## 54 6 1 0.455 0.1501 0.238 0.868
## 59 5 1 0.364 0.1450 0.166 0.795Para calcular o Tempo Mediano para o grupo 1, vamos utilizar a fórmula:
\[ \frac{51-49}{0,455-0,545}=\frac{Md-49}{0,5-0,545}\rightarrow Md=48\mbox{ dias} \]
Para calcular o Tempo Mediano para o grupo 3, vamos utilizar a fórmula:
\[ \frac{51-49}{0,455-0,545}=\frac{Md-49}{0,5-0,545}\rightarrow Md=48\mbox{ dias} \] Função para Obter o valor mediano da sobrevivência com mais precisão
calcular_Md(ekm)## [1] 46calcular_Md(ena)## [1] 55.29998O Tempo Médio
Vamos cálcular o tempo médio pelo indicador Kaplan-Meier.
t<- tempos[cens==1]
tj<-c(0,as.numeric(levels(as.factor(t))))
surv<-c(1,as.numeric(levels(as.factor(ekm$surv))))
surv<-sort(surv, decreasing=T)
k<-length(tj)-1
prod<-matrix(0,k,1)
for(j in 1:k){
prod[j]<-(tj[j+1]-tj[j])*surv[j]
}
tm<-sum(prod)
tm## [1] 47.68182Vamos cálcular o tempo médio pelo indicador Nelson-Aalen.
t<- tempos[cens==1]
tj<-c(0,as.numeric(levels(as.factor(t))))
surv<-c(1,as.numeric(levels(as.factor(ena$surv))))
surv<-sort(surv, decreasing=T)
k<-length(tj)-1
prod<-matrix(0,k,1)
for(j in 1:k){
prod[j]<-(tj[j+1]-tj[j])*surv[j]
}
tm<-sum(prod)
tm## [1] 49.81581Comparando as curvas por Log-rank e Peto
Vamos comparar as curvas pelo teste de Log-rank
Grupo 1 vs grupo 2
survdiff(formula=Surv(tempos[1:21],cens[1:21])~grupos[1:21],rho=0)## Call:
## survdiff(formula = Surv(tempos[1:21], cens[1:21]) ~ grupos[1:21],
## rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## grupos[1:21]=1 11 7 5.38 0.488 0.978
## grupos[1:21]=2 10 4 5.62 0.467 0.978
##
## Chisq= 1 on 1 degrees of freedom, p= 0.3Grupo 1 vs grupo 2
survdiff(formula=Surv(tempos[21:32],cens[21:32])~grupos[21:32],rho=0)## Call:
## survdiff(formula = Surv(tempos[21:32], cens[21:32]) ~ grupos[21:32],
## rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## grupos[21:32]=2 1 1 0.644 0.1966 0.228
## grupos[21:32]=3 11 7 7.356 0.0172 0.228
##
## Chisq= 0.2 on 1 degrees of freedom, p= 0.6Grupo 1 vs grupo 3
survdiff(formula = Surv(c(tempos[1:11], tempos[22:32]), c(cens[1:11],
cens[22:32])) ~ c(grupos[1:11], grupos[22:32]), rho = 0)## Call:
## survdiff(formula = Surv(c(tempos[1:11], tempos[22:32]), c(cens[1:11],
## cens[22:32])) ~ c(grupos[1:11], grupos[22:32]), rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## c(grupos[1:11], grupos[22:32])=1 11 7 7.26 0.00922 0.0196
## c(grupos[1:11], grupos[22:32])=3 11 7 6.74 0.00993 0.0196
##
## Chisq= 0 on 1 degrees of freedom, p= 0.9Vamos comparar as curvas pelo teste de Peto
survdiff(formula = Surv(tempos, cens) ~ grupos, rho = 1)## Call:
## survdiff(formula = Surv(tempos, cens) ~ grupos, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## grupos=1 11 4.97 4.56 0.0362 0.0723
## grupos=2 10 3.09 4.53 0.4560 0.9241
## grupos=3 11 5.34 4.31 0.2466 0.4792
##
## Chisq= 1 on 2 degrees of freedom, p= 0.6