Descripción estadística de datos de supervivencia
Julián Estivent Serna Triana
19/6/2022
Transformando datos
t.exp<- (lung$time)/30
t.plant<- lung$age
hibrido<- lung$sex
severidad <- lung$ph.ecog
estado<- lung$status
N17 <- lung$meal.cal
P17<- lung$wt.loss
sup.caf<-data.frame(hibrido,t.exp,t.plant, estado, severidad, N17, P17 )
#View(sup.caf)Creando data frame como si fuera para cultivo de café
caf.sur<-Surv(sup.caf$t.exp, sup.caf$estado)modelado no- parámetrico
caf.fit<-survfit(caf.sur~1)
summary(caf.fit)## Call: survfit(formula = caf.sur ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.167 228 1 0.9956 0.00438 0.9871 1.000
## 0.367 227 3 0.9825 0.00869 0.9656 1.000
## 0.400 224 1 0.9781 0.00970 0.9592 0.997
## 0.433 223 2 0.9693 0.01142 0.9472 0.992
## 0.500 221 1 0.9649 0.01219 0.9413 0.989
## 0.867 220 1 0.9605 0.01290 0.9356 0.986
## 1.000 219 1 0.9561 0.01356 0.9299 0.983
## 1.033 218 1 0.9518 0.01419 0.9243 0.980
## 1.767 217 2 0.9430 0.01536 0.9134 0.974
## 1.800 215 1 0.9386 0.01590 0.9079 0.970
## 1.967 214 1 0.9342 0.01642 0.9026 0.967
## 2.000 213 2 0.9254 0.01740 0.8920 0.960
## 2.033 211 1 0.9211 0.01786 0.8867 0.957
## 2.067 210 1 0.9167 0.01830 0.8815 0.953
## 2.167 209 2 0.9079 0.01915 0.8711 0.946
## 2.367 207 1 0.9035 0.01955 0.8660 0.943
## 2.633 206 1 0.8991 0.01995 0.8609 0.939
## 2.700 205 2 0.8904 0.02069 0.8507 0.932
## 2.933 203 2 0.8816 0.02140 0.8406 0.925
## 3.067 201 1 0.8772 0.02174 0.8356 0.921
## 3.100 199 1 0.8728 0.02207 0.8306 0.917
## 3.167 198 2 0.8640 0.02271 0.8206 0.910
## 3.500 196 1 0.8596 0.02302 0.8156 0.906
## 3.567 194 2 0.8507 0.02362 0.8056 0.898
## 3.667 192 1 0.8463 0.02391 0.8007 0.894
## 3.867 191 1 0.8418 0.02419 0.7957 0.891
## 3.933 190 1 0.8374 0.02446 0.7908 0.887
## 4.067 189 1 0.8330 0.02473 0.7859 0.883
## 4.367 188 1 0.8285 0.02500 0.7810 0.879
## 4.400 187 2 0.8197 0.02550 0.7712 0.871
## 4.500 185 1 0.8153 0.02575 0.7663 0.867
## 4.733 184 1 0.8108 0.02598 0.7615 0.863
## 4.800 183 1 0.8064 0.02622 0.7566 0.859
## 4.833 182 2 0.7975 0.02667 0.7469 0.852
## 4.900 180 1 0.7931 0.02688 0.7421 0.848
## 5.100 179 1 0.7887 0.02710 0.7373 0.844
## 5.200 178 2 0.7798 0.02751 0.7277 0.836
## 5.433 176 3 0.7665 0.02809 0.7134 0.824
## 5.533 173 2 0.7577 0.02845 0.7039 0.816
## 5.567 171 1 0.7532 0.02863 0.6991 0.811
## 5.667 170 1 0.7488 0.02880 0.6944 0.807
## 5.833 167 1 0.7443 0.02898 0.6896 0.803
## 5.867 165 1 0.7398 0.02915 0.6848 0.799
## 5.900 164 1 0.7353 0.02932 0.6800 0.795
## 5.967 162 2 0.7262 0.02965 0.6704 0.787
## 6.000 160 1 0.7217 0.02981 0.6655 0.783
## 6.033 159 2 0.7126 0.03012 0.6559 0.774
## 6.067 157 1 0.7081 0.03027 0.6511 0.770
## 6.100 156 1 0.7035 0.03041 0.6464 0.766
## 6.200 154 1 0.6989 0.03056 0.6416 0.761
## 6.300 152 1 0.6943 0.03070 0.6367 0.757
## 6.467 149 1 0.6897 0.03085 0.6318 0.753
## 6.567 147 1 0.6850 0.03099 0.6269 0.749
## 6.633 145 1 0.6803 0.03113 0.6219 0.744
## 6.700 144 2 0.6708 0.03141 0.6120 0.735
## 6.733 142 1 0.6661 0.03154 0.6071 0.731
## 6.900 139 1 0.6613 0.03168 0.6020 0.726
## 6.933 138 1 0.6565 0.03181 0.5970 0.722
## 7.000 137 1 0.6517 0.03194 0.5920 0.717
## 7.067 135 1 0.6469 0.03206 0.5870 0.713
## 7.267 134 1 0.6421 0.03218 0.5820 0.708
## 7.400 132 1 0.6372 0.03231 0.5769 0.704
## 7.433 130 1 0.6323 0.03243 0.5718 0.699
## 7.533 126 1 0.6273 0.03256 0.5666 0.694
## 7.633 125 1 0.6223 0.03268 0.5614 0.690
## 7.667 124 1 0.6172 0.03280 0.5562 0.685
## 7.967 121 2 0.6070 0.03304 0.5456 0.675
## 8.167 117 1 0.6019 0.03316 0.5402 0.670
## 8.200 116 1 0.5967 0.03328 0.5349 0.666
## 8.900 112 1 0.5913 0.03341 0.5294 0.661
## 8.933 111 1 0.5860 0.03353 0.5239 0.656
## 8.967 110 1 0.5807 0.03364 0.5184 0.651
## 9.000 108 1 0.5753 0.03376 0.5128 0.645
## 9.433 104 1 0.5698 0.03388 0.5071 0.640
## 9.467 103 1 0.5642 0.03400 0.5014 0.635
## 9.500 101 2 0.5531 0.03424 0.4899 0.624
## 9.533 99 1 0.5475 0.03434 0.4841 0.619
## 9.600 98 1 0.5419 0.03444 0.4784 0.614
## 9.700 97 1 0.5363 0.03454 0.4727 0.608
## 9.767 94 1 0.5306 0.03464 0.4669 0.603
## 10.033 91 1 0.5248 0.03475 0.4609 0.597
## 10.100 89 1 0.5189 0.03485 0.4549 0.592
## 10.167 87 1 0.5129 0.03496 0.4488 0.586
## 10.200 86 1 0.5070 0.03506 0.4427 0.581
## 10.333 85 2 0.4950 0.03523 0.4306 0.569
## 10.667 82 1 0.4890 0.03532 0.4244 0.563
## 10.967 81 1 0.4830 0.03539 0.4183 0.558
## 11.233 79 1 0.4768 0.03547 0.4121 0.552
## 11.333 78 1 0.4707 0.03554 0.4060 0.546
## 11.500 77 1 0.4646 0.03560 0.3998 0.540
## 11.600 76 1 0.4585 0.03565 0.3937 0.534
## 11.667 75 1 0.4524 0.03569 0.3876 0.528
## 11.700 74 1 0.4463 0.03573 0.3815 0.522
## 11.767 73 2 0.4340 0.03578 0.3693 0.510
## 12.033 70 1 0.4278 0.03581 0.3631 0.504
## 12.100 69 2 0.4154 0.03583 0.3508 0.492
## 12.133 67 1 0.4092 0.03582 0.3447 0.486
## 12.367 65 2 0.3966 0.03581 0.3323 0.473
## 12.900 60 1 0.3900 0.03582 0.3258 0.467
## 13.000 59 1 0.3834 0.03582 0.3193 0.460
## 13.133 58 1 0.3768 0.03580 0.3128 0.454
## 14.200 55 1 0.3700 0.03580 0.3060 0.447
## 14.267 54 1 0.3631 0.03579 0.2993 0.440
## 14.300 53 1 0.3563 0.03576 0.2926 0.434
## 14.433 52 1 0.3494 0.03573 0.2860 0.427
## 14.733 51 1 0.3426 0.03568 0.2793 0.420
## 14.800 50 1 0.3357 0.03561 0.2727 0.413
## 15.000 48 1 0.3287 0.03555 0.2659 0.406
## 15.167 47 1 0.3217 0.03548 0.2592 0.399
## 15.233 46 1 0.3147 0.03539 0.2525 0.392
## 15.333 44 1 0.3076 0.03530 0.2456 0.385
## 15.767 43 1 0.3004 0.03520 0.2388 0.378
## 15.900 42 1 0.2933 0.03508 0.2320 0.371
## 17.300 39 1 0.2857 0.03498 0.2248 0.363
## 17.333 38 1 0.2782 0.03485 0.2177 0.356
## 17.467 37 2 0.2632 0.03455 0.2035 0.340
## 17.767 34 1 0.2554 0.03439 0.1962 0.333
## 18.333 32 1 0.2475 0.03423 0.1887 0.325
## 18.600 30 1 0.2392 0.03407 0.1810 0.316
## 18.900 28 1 0.2307 0.03391 0.1729 0.308
## 19.133 27 1 0.2221 0.03371 0.1650 0.299
## 19.433 26 1 0.2136 0.03348 0.1571 0.290
## 20.433 24 1 0.2047 0.03325 0.1489 0.281
## 20.800 23 1 0.1958 0.03297 0.1407 0.272
## 21.367 22 1 0.1869 0.03265 0.1327 0.263
## 21.433 21 1 0.1780 0.03229 0.1247 0.254
## 21.800 20 1 0.1691 0.03188 0.1169 0.245
## 21.833 19 1 0.1602 0.03142 0.1091 0.235
## 22.900 18 1 0.1513 0.03090 0.1014 0.226
## 22.967 17 1 0.1424 0.03034 0.0938 0.216
## 23.500 16 1 0.1335 0.02972 0.0863 0.207
## 23.567 15 1 0.1246 0.02904 0.0789 0.197
## 24.267 14 1 0.1157 0.02830 0.0716 0.187
## 24.367 13 1 0.1068 0.02749 0.0645 0.177
## 24.500 12 1 0.0979 0.02660 0.0575 0.167
## 25.500 10 1 0.0881 0.02568 0.0498 0.156
## 26.367 9 1 0.0783 0.02462 0.0423 0.145
## 27.133 7 1 0.0671 0.02351 0.0338 0.133
## 29.433 4 1 0.0503 0.02285 0.0207 0.123Valores únicos
plot(caf.fit,xlab="Tiempo de supervivencia (meses)",ylab="Proporción de plantas vivas")
title("Curva de supervivenvia de proporción de individuos vs tiempo")
abline(h = 0.5, col='red')
abline(v = 10.32, col='red')
abline(h = c(0.25, 0.75), col='blue')
abline(v = c(5.65, 18.5), col='blue')
\[El\ tiempo\ requerido\ para\ una\
supervivencia\ del\ 75\%\ es\ de\ 5.65\ meses,\\ para\ que\ mueran\ el\
50\%\ de\ las\ plantas\ 10.3\ meses,\\ y\ la\ mortalidad\ del\ 75\%\ se\
alcanza\ a\ los\ 18.5\ meses\ después\ de\ siembra.\] ###
Intervalo
plot(caf.fit,xlab="Tiempo de supervivencia (meses)",ylab="Proporción de plantas vivas")
title("Curva de supervivenvia de proporción de individuos vs tiempo (zona de confianza)", cex.main = 1 )
abline(h = 0.5, col='red')
abline(v = 10.32, col='red')
points(c(10.32, 10.32), c(0.44, 0.58), pch =16, col='blue')
points(c(9.5, 12.2), c(0.5, 0.5), pch =16, col='red')
\[La\ zona\
de\ confianza\ indica\ una\ probabilidad\ de\ supervivencia\ entre\ el\ 44\%\ y\ 58\%\ a\ los\ 10.3\ meses,\\ mientras\ que\ entre\ los\ 9.5\ y\ 12.2\ meses\ la\ probabilidad\ de\ supervivencia\ es\ del\ 50\%\]
## Modelado no parámetrico
caf.fit.strata<-survfit(caf.sur~severidad,sup.caf)
plot(caf.fit.strata, lty = 1:4,col=1:4,xlab="Tiempo de supervivencia (meses)",ylab="Proporción de plantas vivas",
lwd=3)
title("Curvas de supervivencia estratificadas por severidad ")
legend(20, .99, c("Severidad=0", "Severidad=1","Severidad=2","Severidad=3"), lty = 1:4,col=1:4, lwd=3)
abline(h = 0.5)
points(c(6.5, 10.32, 13.2), c(0.5, 0.5, 0.5 ), pch =16, col='red')
\[\small El\ grado\ de\ severidad\ 3\ lleva\
la\ proporción\ de\ plantas\ vivas\ a\ 0\%\ antes\ de\ los\ 5\ meses,\\
\small para\ la\ severidad\ de\ grado\ 2\ el\ 50\%\ de\ plantas\
muertas\ se\ alcanza\ a los\ 6.5\ meses,\\ \small para\ el\ grado\ 1\ a\
los\ 10.3\ meses\ y\ para\ el\ grado\ 0\ a\ los\ 13.2 meses.\\ \small
Exceptuando\ el\ mayor\ y\ el\ menor\ grado\ de\ severidad,\ los\
grados\ 1\ y\ 2\ caen\ dentro\ de\ la\ zona\ de\ confianza\]
caf.fit.strata<-survfit(caf.sur~hibrido,sup.caf)
plot(caf.fit.strata, col=1:2, xlab = 'meses', lwd=3)
title("Curvas de supervivencia estratificadas por Híbrido ")
legend(20, .9, c("Híbrido 1", "Híbrido 2"), col=1:2, lwd=3)
abline(h = 0.5)
abline(h = 0.25)
abline(h = 0.08)
points(c(9, 14.2, 15.3, 23, 25.5 ), c(0.5, 0.5, 0.25, 0.25, 0.08), pch =16, col='red')
\[\small El\ híbrido\ 2\ tiene\ mejor\
proporción\ de\ supervivencia\ a\ lo\ largo\ del\ tiempo,\ el\ 50\%\ se\
alcanza\ a\ los\ 9\ meses\ en\ el\ híbrido\ 1\ y\ a\ los\ 14.2\ en\ el\
híbrido\ 2.\\ \small El\ 75\%\ de\ plantas\ muertas\ para\ el\ híbrido\
1\ se\ presenta\ a\ los\ 15.3\ meses\ y\ para\ el\ hibrido\ 2\ a\ los\
23\ meses.\\ \small Finalmente\ ambos\ convergen\ a\ los\ 25.5\ meses,\
dónde\ alcanzan\ una\ proporción\ del\ 92\%\ de\ plantas\ muertas\\
\small en\ el\ que\ el\ híbrido\ 2\ se\ mantiene\ y\ el\ híbrido\ 1\
sigue\ aumentando\ hasta\ el\ fin\ del\ experimento \] ##
Intervalos de Confianza para el estimador Kaplan-Meier
library(km.ci)## Warning: package 'km.ci' was built under R version 4.1.3a<-km.ci(caf.fit, conf.level=0.95, tl=NA, tu=NA, method="loghall")
plot(a, lty=2, lwd=2,xlab = 'meses', col = 'red')
lines(caf.fit, lwd=2, lty=1, col = 'green')
lines(caf.fit, lwd=1, lty=4, conf.int=T, col = 'black')
linetype<-c(1, 2, 4)
title("Intervalos de confianza para distintos estimadores ")
legend(x = "topright", .9, c("Kaplan-Meier", "Hall-Wellner", "Pointwise"),
lty = linetype,
col = c('red', 'green', 'black'))
abline(h = 0.5, col='maroon3', lwd=2)
abline(v = 10.32, col='maroon3', lwd=2)
\[El~ estimador~ Hall-Wellner~ tiene~
intervalos~ más~ angostos,\\ por~ lo~ que~ podría~ dar~ resultados~ más~
precisos~ para~ describir~ la~ curva~ de~ supervivencia\] ##
Estimador de riesgo acumulado
aalen.fit<- survfit(coxph(caf.sur~1), type="aalen")
sum_aalen.fit = summary(aalen.fit)
plot(aalen.fit,xlab = 'meses', col="red",lwd=1,lty=1)
title("Estimadores de riesgo acumulado ")
lines(caf.fit, lwd=1, lty=1)
legend("topright", .9,
c("Nelson-Aalen", "Kaplan-Meier"),
lty=c(1,1),
col=c("red", "black"))
\[Ambos~ estimadores~ son~ muy~
similares\]
barplot(sum_aalen.fit$time, cumsum(sum_aalen.fit$n.event))
title("Gráfico de barras de las muertes acumuladas en el tiempo ")
mod_suv = lm(cumsum(sum_aalen.fit$n.event) ~ sum_aalen.fit$time)
summary(mod_suv)##
## Call:
## lm(formula = cumsum(sum_aalen.fit$n.event) ~ sum_aalen.fit$time)
##
## Residuals:
## Min 1Q Median 3Q Max
## -49.044 -11.535 4.049 12.868 20.208
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.1780 2.1715 10.21 <2e-16 ***
## sum_aalen.fit$time 6.5187 0.1773 36.76 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14.43 on 137 degrees of freedom
## Multiple R-squared: 0.908, Adjusted R-squared: 0.9073
## F-statistic: 1352 on 1 and 137 DF, p-value: < 2.2e-16plot(sum_aalen.fit$time, cumsum(sum_aalen.fit$n.event),xlab="Tiempo de supervivencia (meses)",ylab="Plantas muertas", pch = 16)
abline(mod_suv)
\[La~ ocurrencia~ de~ muertes~ se~ puede~ estimar~ como~ producto~ de~ una~ regresión~ lineal,\\
ya~ que~ el~ r^2~ es~ de~ 0.90\]
survdiff(caf.sur~severidad,sup.caf)## Call:
## survdiff(formula = caf.sur ~ severidad, data = sup.caf)
##
## n=227, 1 observation deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## severidad=0 63 37 54.153 5.4331 8.2119
## severidad=1 113 82 83.528 0.0279 0.0573
## severidad=2 50 44 26.147 12.1893 14.6491
## severidad=3 1 1 0.172 3.9733 4.0040
##
## Chisq= 22 on 3 degrees of freedom, p= 7e-05# Prueba de log-rank or Mantel-Haenszel
survdiff(caf.sur~hibrido,sup.caf, rho = 0)## Call:
## survdiff(formula = caf.sur ~ hibrido, data = sup.caf, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## hibrido=1 138 112 91.6 4.55 10.3
## hibrido=2 90 53 73.4 5.68 10.3
##
## Chisq= 10.3 on 1 degrees of freedom, p= 0.001# Preuba de Peto & Peto modification of the Gehan-Wilcoxon test
survdiff(caf.sur~hibrido,sup.caf, rho = 1)## Call:
## survdiff(formula = caf.sur ~ hibrido, data = sup.caf, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## hibrido=1 138 70.4 55.6 3.95 12.7
## hibrido=2 90 28.7 43.5 5.04 12.7
##
## Chisq= 12.7 on 1 degrees of freedom, p= 4e-04survdiff(caf.sur~hibrido + severidad,sup.caf)## Call:
## survdiff(formula = caf.sur ~ hibrido + severidad, data = sup.caf)
##
## n=227, 1 observation deleted due to missingness.
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## hibrido=1, severidad=0 36 28 33.051 0.772 0.986
## hibrido=1, severidad=1 71 54 43.318 2.634 3.636
## hibrido=1, severidad=2 29 28 14.416 12.799 14.128
## hibrido=1, severidad=3 1 1 0.172 3.973 4.004
## hibrido=2, severidad=0 27 9 21.101 6.940 8.020
## hibrido=2, severidad=1 42 28 40.210 3.707 4.999
## hibrido=2, severidad=2 21 16 11.731 1.553 1.693
##
## Chisq= 32.9 on 6 degrees of freedom, p= 1e-05\[-Se~ rechaza~ la~ H_0~ de~ igualdad~ de~ supervivencia~ en~ los~ 4~ grupos~ funcionales\\ -Se~ rechaza~ la~ H_0~ de~ igualdad~ de~ supervivencia~ entre~ los~ tipos~ de~ híbridos~\\ - El~ híbrido~ 2~ tiene~ mejor~ supervivencia~ de ~acuerdo~ a ~las~ pruebas~ y~ los~ estimadores,\\ por~ lo~ que~ sería~ el~ mejor~ a~ utilizar~ para~ llevar~ a~ cabo~ el~ cultivo\]
caf.fit.strata<-survfit(caf.sur~hibrido,sup.caf)
plot(caf.fit.strata, conf.int = 0.95,
col=1:2, xlab = 'meses', lwd=1)
title("Curvas de supervivencia estratificadas por híbrido y estimadores ")
legend("topright", .9, c("Hibrido 1", "Hibrido 2"), col=1:2, lwd=3)
# abline(v = 400)
abline(h = 0.25, col = 'maroon2', lwd=2)
abline(h = 0.5, col = 'maroon2')
abline(v = c(7, 10.66), col = 1)
abline(v = c(11.66, 18.3), col = 2)
\[-Para~ el~ hibrido~ el~ 50\%~ de~
supervivencia~ se~ da~ entre~ los~ 7~ y~ 10.6~ meses\\ - Para ~el~
hibrido~ 2~ el~ 50\%~ de~ supervivencia~ se~ da~ entre~ 11.6~ y~ 18.3~
meses\]
Modelo Parametrico
par.wei<-survreg(caf.sur~1,dist="w")
par.wei## Call:
## survreg(formula = caf.sur ~ 1, dist = "w")
##
## Coefficients:
## (Intercept)
## 2.633707
##
## Scale= 0.7593936
##
## Loglik(model)= -592.7 Loglik(intercept only)= -592.7
## n= 228kappa<-par.wei$scale
lambda<-exp(-par.wei$coeff[1])
zeit<-seq(from=0,to=36.6,length.out=33.3)
s<-exp(-(lambda*zeit)^kappa)
h<-lambda^kappa *kappa*zeit^(kappa-1)
par(mfrow=c(2,1))
plot(zeit,h,xlab="Meses",ylab="h(t)", pch = 16, cex = 0.6, las = 1)
plot(zeit,s,xlab="Meses",ylab="s(t)", pch = 16, cex = 0.6, las = 1)