Curva de supervivencia a Roya naranja
Julián Estivent Serna Triana
15/10/2022
Carga de datos
datos_xp <- read_excel("C:/Users/admin/Downloads/Dataset Sancarlos (1).xlsx")
#View(datos_xp)
data<-data.frame(datos_xp$`Tiempo de evaluacion`, datos_xp$Variedad, datos_xp$`Edad (meses)`, datos_xp$`Severidad RN`, datos_xp$`Reaccion RN`, datos_xp$`ID RN`)
datos1 <- data%>%
group_by(datos_xp..Tiempo.de.evaluacion., datos_xp.Variedad) %>%
summarise_at(vars(datos_xp..ID.RN., datos_xp..Reaccion.RN., datos_xp..Edad..meses..),
list(name = mean))
datos<- na.omit(datos1)Transformando datos
t.ev<- datos$datos_xp..Tiempo.de.evaluacion.
edad<- datos$datos_xp..Edad..meses.._name
Variedad <- datos$datos_xp.Variedad
ID<- datos$datos_xp..ID.RN._name
R <- datos$datos_xp..Reaccion.RN._name
EstadoID <- as.integer(with(datos, ifelse(ID < 500.1, "0",
ifelse(ID >= 500.1 & ID < 512.01, "1", "2"))))
Estado <- as.integer(if_else(datos$datos_xp..Reaccion.RN._name> 5, "2", "1"))
sup.caf1<-data.frame(t.ev, edad, Variedad, ID, R, Estado, EstadoID)
sup.caf<- na.omit(sup.caf1)
#View(sup.caf)Creando surv frame para caña
caf.sur<-Surv(sup.caf$t.ev, sup.caf$Estado)Kaplan - Meier
cafkm_fit <- survfit(caf.sur ~ 1)
summary(cafkm_fit, times = c(1, 2, 3, 4, 5 , 6 ,7 , 8, 9, 10, 11, 12, 13, 14, 15, 16))## Call: survfit(formula = caf.sur ~ 1)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 180 3 0.983 0.00954 0.965 1.000
## 2 159 2 0.971 0.01282 0.946 0.996
## 3 144 2 0.957 0.01579 0.927 0.989
## 4 132 2 0.943 0.01859 0.907 0.980
## 5 119 2 0.927 0.02139 0.886 0.970
## 6 105 2 0.909 0.02436 0.863 0.958
## 7 90 5 0.859 0.03180 0.799 0.924
## 8 72 5 0.799 0.03922 0.726 0.880
## 9 57 5 0.729 0.04666 0.643 0.827
## 10 43 3 0.678 0.05183 0.584 0.788
## 11 27 1 0.653 0.05567 0.553 0.772
## 12 19 0 0.653 0.05567 0.553 0.772
## 13 19 0 0.653 0.05567 0.553 0.772
## 14 13 1 0.603 0.07050 0.479 0.758
## 15 7 0 0.603 0.07050 0.479 0.758
## 16 2 0 0.603 0.07050 0.479 0.758plot(cafkm_fit,xlab="Tiempo de supervivencia (Tiempo de evaluación)",ylab="Proporción de plantas vivas")
title("Curva de supervivenvia de proporción de individuos vs tiempo")
abline(h = 0.5, col='red')
abline(h = c(0.25, 0.75), col='blue')
abline(v = c(9), col='blue')
\[- El\ tiempo\ requerido\ para\ que\ el\
25\%\ de\ las\ plantas\ pasen\ a\ ser\ susceptibles\\ es\ de\ 9\
tiempos\ de\ muestreo.\] ## Intervalo de confianza
plot(cafkm_fit,xlab="Tiempo de evaluación",ylab="Proporción de plantas resistentes")
title("Curva de supervivenvia de proporción de individuos vs tiempo (zona de confianza)", cex.main = 1 )
abline(h = 0.75, col='red')
abline(v = 9, col='red')
points(c(9, 9), c(0.88, 0.65), pch =16, col='blue')
points(c(8, 16), c(0.75, 0.75), pch =16, col='red')
\[La\ zona\ de\ confianza\ indica\ una\
proporción\ de\ susceptibilidad\ entre\ 12\%\ y\ 35\%\ a\ los\ 8\
tiempos\ de\ evaluación,\\ desde\ los\ 8\ tiempos\ en\ adelante\ la\
proporción\ de\ plantas\ susceptibles\ puede\ ser\ del\ 75%. \]
## Modelado no parámetrico
ifit <- survfit(Surv(sup.caf$t.ev, sup.caf$Estado)~EstadoID, data=sup.caf)
ifit## Call: survfit(formula = Surv(sup.caf$t.ev, sup.caf$Estado) ~ EstadoID,
## data = sup.caf)
##
## n events median 0.95LCL 0.95UCL
## EstadoID=0 135 0 NA NA NA
## EstadoID=1 13 1 NA NA NA
## EstadoID=2 32 32 7 6 9summary(ifit)## Call: survfit(formula = Surv(sup.caf$t.ev, sup.caf$Estado) ~ EstadoID,
## data = sup.caf)
##
## EstadoID=0
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
##
## EstadoID=1
## time n.risk n.event survival std.err lower 95% CI
## 6.000 9.000 1.000 0.889 0.105 0.706
## upper 95% CI
## 1.000
##
## EstadoID=2
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 1 32 3 0.9062 0.0515 0.81068 1.000
## 2 29 2 0.8438 0.0642 0.72688 0.979
## 3 27 2 0.7812 0.0731 0.65038 0.938
## 4 25 2 0.7188 0.0795 0.57870 0.893
## 5 23 2 0.6562 0.0840 0.51070 0.843
## 6 21 1 0.6250 0.0856 0.47789 0.817
## 7 20 5 0.4688 0.0882 0.32415 0.678
## 8 15 5 0.3125 0.0819 0.18692 0.522
## 9 10 5 0.1562 0.0642 0.06985 0.350
## 10 5 3 0.0625 0.0428 0.01633 0.239
## 11 2 1 0.0312 0.0308 0.00454 0.215
## 14 1 1 0.0000 NaN NA NAsurvp<- ggsurvplot(ifit, conf.int=TRUE, pval=TRUE, risk.table=TRUE,
legend.labs=c("Resistente", "Intermedio", "Susceptible"), legend.title="Reacción",
palette=c("dodgerblue2", "orchid2", "brown2"),
title="Supervivencia de roya naranja en caña de acuerdo al ID",
risk.table.height=.35
)
survp
survp2<- survp$plot + geom_hline(yintercept = 0.5)+
geom_hline(yintercept = 0.82)+
geom_hline(yintercept = 0.48)+
geom_vline(xintercept = 7)+
geom_vline(xintercept = 9)+
geom_vline(xintercept = 6)
survp2
\[\small-Las\ variedades\ susceptibles\
alcanzan\ el\ 50\%\ de\ supervivencia\\ \small a\ los\ 7\ tiempos\ de\
evaluación,\ según\ el\ intervalo\ de\ confianza\ ese\ porcentaje\ de\
supervivencia\ se\ puede\ alcanzar\ entre\ los\ 6\ y\ 9\
tiempos.\] ## Intervalos de confianza
library(km.ci)
a<-km.ci(cafkm_fit, conf.level=0.95, tl=NA, tu=NA, method="loghall")
plot(a, lty=2, lwd=2,xlab = 'Tiempo de evaluación', col = 'red')
lines(cafkm_fit, lwd=2, lty=1, col = 'green')
lines(cafkm_fit, lwd=1, lty=4, conf.int=T, col = 'black')
linetype<-c(1, 2, 4)
title("Intervalos de confianza para distintos estimadores ")
legend(x = "bottomright", .7, c("Kaplan-Meier", "Hall-Wellner", "Pointwise"),
lty = linetype,
col = c('red', 'green', 'black'))
abline(h = 0.7, col='maroon3', lwd=2)
abline(v = 10, col='maroon3', lwd=2)
## Regresión de Cox - Estimador de riesgo acumulado
aalen.fit<- survfit(coxph(caf.sur~1), type="aalen")
plot(aalen.fit,xlab="Tiempo de evaluación",ylab="Proporción de plantas resistentes",col="red",lwd=1,lty=1)
title("Comparación Kaplan-Meier / Nelson-Aalen")
lines(cafkm_fit, lwd=1, lty=1)
legend(10, .2, c("Nelson-Aalen","Kaplan-Meier"), lty=c(1,1),
col=c("red","black"))
sum_aalen.fit = summary(aalen.fit)
png(file="C:/Users/admin/Desktop/Diseno Experimentos/Rn personal/g5.png",
width=600, height=350)
barplot(sum_aalen.fit$time, cumsum(sum_aalen.fit$n.event))
title("Gráfico de barras de muertes acumuladas en el tiempo ")
dev.off()## png
## 2mod_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
## -5.141 -1.530 -0.164 1.707 3.981
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.7875 1.6671 -0.472 0.647
## sum_aalen.fit$time 2.7806 0.2180 12.757 1.64e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.831 on 10 degrees of freedom
## Multiple R-squared: 0.9421, Adjusted R-squared: 0.9363
## F-statistic: 162.7 on 1 and 10 DF, p-value: 1.64e-07plot(sum_aalen.fit$time, cumsum(sum_aalen.fit$n.event),xlab="Tiempo de evaluación",ylab="Individuos muertos", pch = 16)
abline(mod_suv, col = "red", lwd = 1.8)
title("Función de supervivencia")
\[\small De\ acuerdo\ al\ R^2(0.93),\ la\
ocurrencia\ de\ plantas\ susceptibles\ se\ puede\ dar\ como\\ \small
producto\ de\ una\ regresión\ lineal.\]
survdiff(caf.sur~EstadoID,sup.caf)## Call:
## survdiff(formula = caf.sur ~ EstadoID, data = sup.caf)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## EstadoID=0 135 0 24.81 24.815 105.365
## EstadoID=1 13 1 2.22 0.668 0.747
## EstadoID=2 32 32 5.97 113.543 145.053
##
## Chisq= 146 on 2 degrees of freedom, p= <2e-16# Prueba de log-rank or Mantel-Haenszel
survdiff(caf.sur~EstadoID,sup.caf, rho = 0)## Call:
## survdiff(formula = caf.sur ~ EstadoID, data = sup.caf, rho = 0)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## EstadoID=0 135 0 24.81 24.815 105.365
## EstadoID=1 13 1 2.22 0.668 0.747
## EstadoID=2 32 32 5.97 113.543 145.053
##
## Chisq= 146 on 2 degrees of freedom, p= <2e-16# Prueba de Peto & Peto modification of the Gehan-Wilcoxon test
survdiff(caf.sur~EstadoID,sup.caf, rho = 1)## Call:
## survdiff(formula = caf.sur ~ EstadoID, data = sup.caf, rho = 1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## EstadoID=0 135 0.000 21.61 21.606 99.490
## EstadoID=1 13 0.927 2.00 0.575 0.706
## EstadoID=2 32 27.994 5.32 96.761 137.668
##
## Chisq= 138 on 2 degrees of freedom, p= <2e-16Modelo parámetrico
par.wei<-survreg(caf.sur~1,dist="w")
par.wei## Call:
## survreg(formula = caf.sur ~ 1, dist = "w")
##
## Coefficients:
## (Intercept)
## 2.953463
##
## Scale= 0.5523486
##
## Loglik(model)= -144.3 Loglik(intercept only)= -144.3
## n= 180kappa<-par.wei$scale
lambda<-exp(-par.wei$coeff[1])
zeit<-seq(from=0,to=20,length.out=16)
s<-exp(-(lambda*zeit)^kappa)
h<-lambda^kappa *kappa*zeit^(kappa-1)
par(mfrow=c(2,1))
plot(zeit,h,xlab="Tiempo de evaluación",ylab="h(t)", pch = 16, cex = 0.6, las = 1)
plot(zeit,s,xlab="Tiempo de evaluación",ylab="s(t)", pch = 16, cex = 0.6, las = 1)
ufit <- survfit(Surv(sup.caf$t.ev, sup.caf$Estado)~Variedad, data=sup.caf)
ufit## Call: survfit(formula = Surv(sup.caf$t.ev, sup.caf$Estado) ~ Variedad,
## data = sup.caf)
##
## n events median 0.95LCL 0.95UCL
## Variedad=CC 00-3257 5 0 NA NA NA
## Variedad=CC 01-1228 7 0 NA NA NA
## Variedad=CC 01-1940 1 1 14.0 NA NA
## Variedad=CC 01-678 6 1 10.0 NA NA
## Variedad=CC 01-746 6 0 NA NA NA
## Variedad=CC 03-154 8 0 NA NA NA
## Variedad=CC 03-469 5 3 8.5 7 NA
## Variedad=CC 05-430 2 0 NA NA NA
## Variedad=CC 06-489 5 0 NA NA NA
## Variedad=CC 09-066 6 1 NA NA NA
## Variedad=CC 09-535 5 0 NA NA NA
## Variedad=CC 09-702 9 9 6.0 4 NA
## Variedad=CC 10-450 9 0 NA NA NA
## Variedad=CC 10-476 6 0 NA NA NA
## Variedad=CC 11-595 9 0 NA NA NA
## Variedad=CC 11-600 12 0 NA NA NA
## Variedad=CC 11-606 7 0 NA NA NA
## Variedad=CC 16-3205 9 0 NA NA NA
## Variedad=CC 84-75 6 4 8.5 7 NA
## Variedad=CC 93-4181 9 0 NA NA NA
## Variedad=CC 93-4418 4 0 NA NA NA
## Variedad=CC 97-7170 9 7 7.0 4 NA
## Variedad=CC 97-7565 6 4 9.0 8 NA
## Variedad=CC00-3257 1 0 NA NA NA
## Variedad=CC01-1228 1 0 NA NA NA
## Variedad=CC01-1940 1 1 1.0 NA NA
## Variedad=CC01-678 1 0 NA NA NA
## Variedad=CC01-746 1 0 NA NA NA
## Variedad=CC03-154 1 0 NA NA NA
## Variedad=CC03-469 1 0 NA NA NA
## Variedad=CC05-430 1 0 NA NA NA
## Variedad=CC06-489 1 0 NA NA NA
## Variedad=CC09-535 1 0 NA NA NA
## Variedad=CC09-702 1 1 1.0 NA NA
## Variedad=CC10-450 1 0 NA NA NA
## Variedad=CC11-595 1 0 NA NA NA
## Variedad=CC11-600 1 0 NA NA NA
## Variedad=CC11-606 1 0 NA NA NA
## Variedad=CC84-75 1 0 NA NA NA
## Variedad=CC85-92 1 0 NA NA NA
## Variedad=CC93-4181 1 0 NA NA NA
## Variedad=CC93-4418 1 0 NA NA NA
## Variedad=CC97-7170 1 1 1.0 NA NA
## Variedad=V7151 9 0 NA NA NAu<- ggsurvplot(ufit, legend="none")
u1<-ggsurvplot(ufit)
g9 <- u$plot +
geom_text_repel(
aes(label = Variedad), data = ufit,
fontface ="plain", color = "black", size = 3, max.overlaps = 3, min.segment.length = 0
)
u1
g9