#SETWD
setwd("~/MAESTRIA EPIDEMIOLOGIA ICESI/ensayos_clinicos/analisis_sobrevida")
library(readxl)
## Warning: package 'readxl' was built under R version 4.5.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.5.3
##
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.5.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.3
library(randomizr)
## Warning: package 'randomizr' was built under R version 4.5.3
library(survival)
## Warning: package 'survival' was built under R version 4.5.3
library(survminer)
## Warning: package 'survminer' was built under R version 4.5.3
## Cargando paquete requerido: ggpubr
## Warning: package 'ggpubr' was built under R version 4.5.2
##
## Adjuntando el paquete: 'survminer'
## The following object is masked from 'package:survival':
##
## myeloma
library(haven)
## Warning: package 'haven' was built under R version 4.5.2
library(stats)
library(epitools)
##
## Adjuntando el paquete: 'epitools'
## The following object is masked from 'package:survival':
##
## ratetable
library(table1)
## Warning: package 'table1' was built under R version 4.5.3
##
## Adjuntando el paquete: 'table1'
## The following objects are masked from 'package:base':
##
## units, units<-
library(ggfortify)
## Warning: package 'ggfortify' was built under R version 4.5.3
library(haven)
EMR <- read_dta("EMR.dta")
View(EMR)
##Explorar base de datos
str(EMR)
## tibble [126 × 11] (S3: tbl_df/tbl/data.frame)
## $ ID : chr [1:126] "1" "2" "3" "4" ...
## ..- attr(*, "label")= chr "ID"
## ..- attr(*, "format.stata")= chr "%9s"
## $ sexo : dbl+lbl [1:126] 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 1, 0, 1, 1, ...
## ..@ label : chr "Genero"
## ..@ format.stata: chr "%10.0g"
## ..@ labels : Named num [1:2] 0 1
## .. ..- attr(*, "names")= chr [1:2] "Masculino" "Femenino"
## $ talla : num [1:126] 167 168 169 155 173 174 162 165 155 NA ...
## ..- attr(*, "label")= chr "Talla"
## ..- attr(*, "format.stata")= chr "%8.0g"
## $ peso : num [1:126] 63 79 65 48 68 69 54 62 48 NA ...
## ..- attr(*, "label")= chr "Peso"
## ..- attr(*, "format.stata")= chr "%8.0g"
## $ edad : num [1:126] 75 51 58 76 63 62 73 77 61 19 ...
## ..- attr(*, "label")= chr "Edad"
## ..- attr(*, "format.stata")= chr "%8.0g"
## $ diaagnosticodeingreso : chr [1:126] "SHOCK SEPTICO DE ORIGEN ABDOMINAL" "FALLA RESPIRATORIA POR OBSTRUCCION DE VIA AEREA SECUNDARIA A LESIO TUMORAL MAS FALLA RESPIRATORIA POR NEUMONIA ASPIRATIVA" "PCR (AESP)" "MASA TIROIDEA EN ESTUDIO CA ANAPLASICO DETIROIDES OBSTRUCCION PARCIAL VAS DISFAGIA SECUNDARIA" ...
## ..- attr(*, "label")= chr "Dx Ingreso"
## ..- attr(*, "format.stata")= chr "%155s"
## $ Clasificaciondiagnostica: dbl+lbl [1:126] 1, 1, 1, 1, 1, 1, 2, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, ...
## ..@ label : chr "ClasiDx"
## ..@ format.stata: chr "%9.0g"
## ..@ labels : Named num [1:2] 1 2
## .. ..- attr(*, "names")= chr [1:2] "Medico" "Quiurgico"
## $ apacheii : num [1:126] 29 35 22 31 24 19 5 28 19 17 ...
## ..- attr(*, "label")= chr "ApacheII"
## ..- attr(*, "format.stata")= chr "%8.0g"
## $ Horasdedestete : num [1:126] 4 30 6 6 1 6 5 6 5 0.2 ...
## ..- attr(*, "label")= chr "Horasdedestete"
## ..- attr(*, "format.stata")= chr "%10.0g"
## $ Tratamiento : dbl+lbl [1:126] 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, ...
## ..@ label : chr "Tratamiento"
## ..@ format.stata: chr "%12.0g"
## ..@ labels : Named num [1:2] 0 1
## .. ..- attr(*, "names")= chr [1:2] "Convencional" "EMR"
## $ extubacion : num [1:126] 1 1 0 1 1 0 0 1 1 1 ...
## ..- attr(*, "format.stata")= chr "%9.0g"
## - attr(*, "label")= chr "An\xe1lsis Ecc"
summary(EMR)
## ID sexo talla peso
## Length:126 Min. :0.0000 Min. :145.0 Min. :45.00
## Class :character 1st Qu.:0.0000 1st Qu.:156.0 1st Qu.:49.00
## Mode :character Median :0.0000 Median :163.0 Median :59.00
## Mean :0.4365 Mean :162.9 Mean :58.14
## 3rd Qu.:1.0000 3rd Qu.:169.0 3rd Qu.:65.00
## Max. :1.0000 Max. :185.0 Max. :80.00
## NA's :1 NA's :1
## edad diaagnosticodeingreso Clasificaciondiagnostica apacheii
## Min. :18.00 Length:126 Min. :1.000 Min. : 5.00
## 1st Qu.:44.50 Class :character 1st Qu.:1.000 1st Qu.:18.00
## Median :62.00 Mode :character Median :1.000 Median :24.00
## Mean :57.49 Mean :1.262 Mean :25.65
## 3rd Qu.:71.00 3rd Qu.:2.000 3rd Qu.:34.00
## Max. :92.00 Max. :2.000 Max. :66.00
##
## Horasdedestete Tratamiento extubacion
## Min. : 0.000 Min. :0.0000 Min. :0.0000
## 1st Qu.: 2.000 1st Qu.:0.0000 1st Qu.:1.0000
## Median : 5.000 Median :0.0000 Median :1.0000
## Mean : 9.069 Mean :0.4921 Mean :0.8175
## 3rd Qu.:11.000 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :72.000 Max. :1.0000 Max. :1.0000
##
head(EMR)
## # A tibble: 6 × 11
## ID sexo talla peso edad diaagnosticodeingreso Clasificaciondiagnos…¹
## <chr> <dbl+lbl> <dbl> <dbl> <dbl> <chr> <dbl+lbl>
## 1 1 0 [Mascu… 167 63 75 SHOCK SEPTICO DE ORI… 1 [Medico]
## 2 2 0 [Mascu… 168 79 51 FALLA RESPIRATORIA P… 1 [Medico]
## 3 3 0 [Mascu… 169 65 58 PCR (AESP) 1 [Medico]
## 4 4 1 [Femen… 155 48 76 MASA TIROIDEA EN EST… 1 [Medico]
## 5 5 0 [Mascu… 173 68 63 FIEBRE DE ORIGEN DES… 1 [Medico]
## 6 6 0 [Mascu… 174 69 62 INSUFICIENCIARESPIRA… 1 [Medico]
## # ℹ abbreviated name: ¹Clasificaciondiagnostica
## # ℹ 4 more variables: apacheii <dbl>, Horasdedestete <dbl>,
## # Tratamiento <dbl+lbl>, extubacion <dbl>
prop.table(table(EMR$extubacion))*100
##
## 0 1
## 18.25397 81.74603
#Crear objetos tratamiento y extubación
EMR$Tratamiento_f <- factor(EMR$Tratamiento,
levels = c(0, 1),
labels = c("Convencional", "EMR"))
EMR$extubacion_f <- factor(EMR$extubacion,
levels = c(0, 1),
labels = c("No", "Si"))
#revisar NA
table(EMR$Tratamiento_f, useNA = "ifany")
##
## Convencional EMR
## 64 62
table(EMR$extubacion_f, useNA = "ifany")
##
## No Si
## 23 103
#Tabla 2x2
tab <- table(EMR$Tratamiento_f, EMR$extubacion_f)
tab
##
## No Si
## Convencional 12 52
## EMR 11 51
prop.table(tab, margin = 1)
##
## No Si
## Convencional 0.1875000 0.8125000
## EMR 0.1774194 0.8225806
#calcular asociación
epitools::riskratio(tab)
## $data
##
## No Si Total
## Convencional 12 52 64
## EMR 11 51 62
## Total 23 103 126
##
## $measure
## risk ratio with 95% C.I.
## estimate lower upper
## Convencional 1.000000 NA NA
## EMR 1.012407 0.8584392 1.19399
##
## $p.value
## two-sided
## midp.exact fisher.exact chi.square
## Convencional NA NA NA
## EMR 0.8874054 1 0.8835695
##
## $correction
## [1] FALSE
##
## attr(,"method")
## [1] "Unconditional MLE & normal approximation (Wald) CI"
promedio <- tapply(EMR$Horasdedestete, EMR$Tratamiento_f, mean, na.rm = TRUE)
mediana <- tapply(EMR$Horasdedestete, EMR$Tratamiento_f, median, na.rm = TRUE)
data.frame(
Tratamiento = names(promedio),
Promedio = as.numeric(promedio),
Mediana = as.numeric(mediana)
)
## Tratamiento Promedio Mediana
## 1 Convencional 8.784375 5
## 2 EMR 9.362903 5
#Según el tipo de tratamiento, el promedio de horas de destete fue de 8.78 horas en el grupo con tratamiento convencional y de 9.36 horas en el grupo EMR. La mediana fue de 5 horas en ambos grupos, es decir que el 50% de los pacientes de ambos grupos se destetaron alrededor de este tiempo. Esto sugiere que el tiempo típico de destete fue similar entre tratamientos, aunque el grupo EMR presentó un promedio ligeramente mayor.
#Discriminado por tipo de tratamiento y por extubacion exitosa o no exitosa
promedio <- aggregate(Horasdedestete ~ Tratamiento_f + extubacion_f,
data = EMR,
FUN = mean,
na.rm = TRUE)
mediana <- aggregate(Horasdedestete ~ Tratamiento_f + extubacion_f,
data = EMR,
FUN = median,
na.rm = TRUE)
names(promedio)[3] <- "Promedio"
names(mediana)[3] <- "Mediana"
resultado <- merge(promedio, mediana,
by = c("Tratamiento_f", "extubacion_f"))
resultado
## Tratamiento_f extubacion_f Promedio Mediana
## 1 Convencional No 12.666667 7
## 2 Convencional Si 7.888462 4
## 3 EMR No 9.500000 5
## 4 EMR Si 9.333333 5
#Discriminando por tipo de tratamiento y extubación exitosa, el grupo Convencional sin extubación exitosa presentó el mayor ztiempo de destete. En el grupo EMR, los tiempos de destete fueron similares entre pacientes con y sin extubación exitosa, con promedios cercanos a 9.5 y 9.33 horas, y medianas de 5 horas.
#KMgeneral
km_general <- survfit(Surv(Horasdedestete, extubacion) ~ 1, data = EMR)
km_general
## Call: survfit(formula = Surv(Horasdedestete, extubacion) ~ 1, data = EMR)
##
## n events median 0.95LCL 0.95UCL
## [1,] 126 103 5 4 7
summary(km_general)
## Call: survfit(formula = Surv(Horasdedestete, extubacion) ~ 1, data = EMR)
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.0 126 3 0.9762 0.0136 0.9499 1.000
## 0.2 122 1 0.9682 0.0157 0.9380 0.999
## 1.0 120 16 0.8391 0.0330 0.7769 0.906
## 2.0 103 16 0.7088 0.0409 0.6330 0.794
## 3.0 86 8 0.6428 0.0432 0.5634 0.733
## 4.0 75 9 0.5657 0.0450 0.4839 0.661
## 5.0 66 8 0.4971 0.0456 0.4152 0.595
## 6.0 56 7 0.4350 0.0456 0.3542 0.534
## 7.0 46 4 0.3972 0.0454 0.3175 0.497
## 8.0 42 3 0.3688 0.0450 0.2904 0.468
## 9.0 38 2 0.3494 0.0447 0.2719 0.449
## 10.0 35 2 0.3294 0.0443 0.2531 0.429
## 11.0 33 4 0.2895 0.0432 0.2161 0.388
## 12.0 29 2 0.2695 0.0425 0.1979 0.367
## 13.0 26 1 0.2592 0.0421 0.1885 0.356
## 15.0 24 2 0.2376 0.0412 0.1690 0.334
## 18.0 21 1 0.2262 0.0408 0.1589 0.322
## 19.0 20 2 0.2036 0.0397 0.1389 0.298
## 22.0 16 2 0.1782 0.0386 0.1165 0.272
## 25.0 12 1 0.1633 0.0382 0.1033 0.258
## 29.0 10 1 0.1470 0.0377 0.0889 0.243
## 30.0 9 3 0.0980 0.0341 0.0495 0.194
## 31.0 6 1 0.0817 0.0321 0.0378 0.176
## 40.0 5 1 0.0653 0.0295 0.0269 0.159
## 49.0 3 1 0.0436 0.0265 0.0132 0.144
## 59.0 2 1 0.0218 0.0203 0.0035 0.136
## 72.0 1 1 0.0000 NaN NA NA
#Probabilidad Extubacion no exitosa
plot(km_general,
xlab = "Horas de destete",
ylab = "Probabilidad de no extubación exitosa",
main = "Kaplan-Meier: muestra total")
#Probabilidad Extubacion exitosa
plot(km_general,
fun = "event",
xlab = "Horas de destete",
ylab = "Probabilidad acumulada de extubación exitosa",
main = "Extubación exitosa acumulada: muestra total")
library(survival)
#ajustar Kaplan-Meier por tratamiento
km_trat <- survfit(Surv(Horasdedestete, extubacion) ~ Tratamiento_f, data = EMR)
#ver resumen
km_trat
## Call: survfit(formula = Surv(Horasdedestete, extubacion) ~ Tratamiento_f,
## data = EMR)
##
## n events median 0.95LCL 0.95UCL
## Tratamiento_f=Convencional 64 52 5 4 8
## Tratamiento_f=EMR 62 51 6 4 11
summary(km_trat)
## Call: survfit(formula = Surv(Horasdedestete, extubacion) ~ Tratamiento_f,
## data = EMR)
##
## Tratamiento_f=Convencional
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0.0 64 1 0.9844 0.0155 0.9545 1.000
## 0.2 63 1 0.9688 0.0217 0.9270 1.000
## 1.0 62 11 0.7969 0.0503 0.7042 0.902
## 2.0 51 6 0.7031 0.0571 0.5996 0.824
## 3.0 45 5 0.6250 0.0605 0.5170 0.756
## 4.0 38 5 0.5428 0.0627 0.4327 0.681
## 5.0 33 3 0.4934 0.0632 0.3839 0.634
## 6.0 29 5 0.4083 0.0627 0.3022 0.552
## 7.0 21 1 0.3889 0.0627 0.2836 0.533
## 8.0 20 2 0.3500 0.0621 0.2472 0.496
## 9.0 17 1 0.3294 0.0618 0.2281 0.476
## 10.0 16 1 0.3088 0.0613 0.2093 0.456
## 11.0 15 2 0.2677 0.0596 0.1730 0.414
## 15.0 12 1 0.2454 0.0587 0.1535 0.392
## 22.0 8 2 0.1840 0.0579 0.0994 0.341
## 25.0 6 1 0.1533 0.0558 0.0752 0.313
## 30.0 5 2 0.0920 0.0474 0.0335 0.253
## 40.0 3 1 0.0613 0.0403 0.0169 0.222
## 59.0 1 1 0.0000 NaN NA NA
##
## Tratamiento_f=EMR
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 62 2 0.9677 0.0224 0.92475 1.000
## 1 58 5 0.8843 0.0411 0.80725 0.969
## 2 52 10 0.7143 0.0587 0.60807 0.839
## 3 41 3 0.6620 0.0616 0.55157 0.795
## 4 37 4 0.5904 0.0645 0.47658 0.731
## 5 33 5 0.5010 0.0660 0.38696 0.649
## 6 27 2 0.4639 0.0661 0.35079 0.613
## 7 25 3 0.4082 0.0655 0.29800 0.559
## 8 22 1 0.3896 0.0651 0.28080 0.541
## 9 21 1 0.3711 0.0646 0.26379 0.522
## 10 19 1 0.3516 0.0641 0.24592 0.503
## 11 18 2 0.3125 0.0626 0.21096 0.463
## 12 16 2 0.2734 0.0606 0.17709 0.422
## 13 13 1 0.2524 0.0595 0.15904 0.401
## 15 12 1 0.2314 0.0581 0.14141 0.379
## 18 11 1 0.2103 0.0565 0.12422 0.356
## 19 10 2 0.1683 0.0525 0.09133 0.310
## 29 5 1 0.1346 0.0516 0.06346 0.286
## 30 4 1 0.1010 0.0485 0.03940 0.259
## 31 3 1 0.0673 0.0424 0.01957 0.231
## 49 2 1 0.0337 0.0319 0.00526 0.215
## 72 1 1 0.0000 NaN NA NA
library(survival)
library(survminer)
library(ggplot2)
km_trat <- survfit(Surv(Horasdedestete, extubacion) ~ Tratamiento_f, data = EMR)
g <- ggsurvplot(
km_trat,
data = EMR,
fun = "event",
conf.int = TRUE,
conf.int.style = "ribbon",
conf.int.alpha = 0.10,
censor = FALSE,
palette = c("black", "forestgreen"),
size = 1.1,
break.time.by = 10,
xlab = "Horas de destete",
ylab = "Probabilidad acumulada de extubación exitosa",
title = "Probabilidad acumulada de extubación exitosa según tratamiento",
legend.title = "Tratamiento",
legend.labs = c("Convencional", "EMR"),
risk.table = TRUE,
risk.table.title = "Número a riesgo",
risk.table.height = 0.18,
risk.table.y.text = FALSE,
tables.theme = theme_cleantable(),
ggtheme = theme_classic()
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
## Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
g$plot <- g$plot +
theme(
plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
axis.title = element_text(size = 12),
axis.text = element_text(size = 10),
legend.title = element_text(size = 10),
legend.text = element_text(size = 10)
)
g$table <- g$table +
theme(
plot.title = element_text(face = "bold", size = 10),
axis.title.x = element_blank(),
axis.text.x = element_text(size = 9),
panel.grid = element_blank()
)
g
## Ignoring unknown labels:
## • colour : "Tratamiento"
#Log rank test
logrank_trat <- survdiff(Surv(Horasdedestete, extubacion) ~ Tratamiento_f, data = EMR)
logrank_trat
## Call:
## survdiff(formula = Surv(Horasdedestete, extubacion) ~ Tratamiento_f,
## data = EMR)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## Tratamiento_f=Convencional 64 52 51 0.0214 0.0475
## Tratamiento_f=EMR 62 51 52 0.0210 0.0475
##
## Chisq= 0 on 1 degrees of freedom, p= 0.8
#No es estadisticamente significativo
cox_trat <- coxph(Surv(Horasdedestete, extubacion) ~ Tratamiento_f, data = EMR)
summary(cox_trat)
## Call:
## coxph(formula = Surv(Horasdedestete, extubacion) ~ Tratamiento_f,
## data = EMR)
##
## n= 126, number of events= 103
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Tratamiento_fEMR -0.0468 0.9543 0.1984 -0.236 0.814
##
## exp(coef) exp(-coef) lower .95 upper .95
## Tratamiento_fEMR 0.9543 1.048 0.6469 1.408
##
## Concordance= 0.515 (se = 0.03 )
## Likelihood ratio test= 0.06 on 1 df, p=0.8
## Wald test = 0.06 on 1 df, p=0.8
## Score (logrank) test = 0.06 on 1 df, p=0.8
cox_edad <- coxph(Surv(Horasdedestete, extubacion) ~ edad, data = EMR)
summary(cox_edad)
## Call:
## coxph(formula = Surv(Horasdedestete, extubacion) ~ edad, data = EMR)
##
## n= 126, number of events= 103
##
## coef exp(coef) se(coef) z Pr(>|z|)
## edad 0.002444 1.002447 0.004875 0.501 0.616
##
## exp(coef) exp(-coef) lower .95 upper .95
## edad 1.002 0.9976 0.9929 1.012
##
## Concordance= 0.496 (se = 0.037 )
## Likelihood ratio test= 0.25 on 1 df, p=0.6
## Wald test = 0.25 on 1 df, p=0.6
## Score (logrank) test = 0.25 on 1 df, p=0.6
cox_apache <- coxph(Surv(Horasdedestete, extubacion) ~ apacheii, data = EMR)
summary(cox_apache)
## Call:
## coxph(formula = Surv(Horasdedestete, extubacion) ~ apacheii,
## data = EMR)
##
## n= 126, number of events= 103
##
## coef exp(coef) se(coef) z Pr(>|z|)
## apacheii 0.011689 1.011758 0.008771 1.333 0.183
##
## exp(coef) exp(-coef) lower .95 upper .95
## apacheii 1.012 0.9884 0.9945 1.029
##
## Concordance= 0.537 (se = 0.037 )
## Likelihood ratio test= 1.74 on 1 df, p=0.2
## Wald test = 1.78 on 1 df, p=0.2
## Score (logrank) test = 1.78 on 1 df, p=0.2
EMR$Clasificaciondiagnostica_f <- factor(
EMR$Clasificaciondiagnostica,
levels = c(1, 2),
labels = c("Medico", "Quirurgico")
)
cox_clasificacion <- coxph(
Surv(Horasdedestete, extubacion) ~ Clasificaciondiagnostica_f,
data = EMR
)
summary(cox_clasificacion)
## Call:
## coxph(formula = Surv(Horasdedestete, extubacion) ~ Clasificaciondiagnostica_f,
## data = EMR)
##
## n= 126, number of events= 103
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Clasificaciondiagnostica_fQuirurgico 0.06501 1.06717 0.22522 0.289 0.773
##
## exp(coef) exp(-coef) lower .95 upper .95
## Clasificaciondiagnostica_fQuirurgico 1.067 0.9371 0.6863 1.659
##
## Concordance= 0.505 (se = 0.027 )
## Likelihood ratio test= 0.08 on 1 df, p=0.8
## Wald test = 0.08 on 1 df, p=0.8
## Score (logrank) test = 0.08 on 1 df, p=0.8
#Proporcionalidad del modelo por tratamiento
proporcionalidad_trat <- cox.zph(cox_trat)
proporcionalidad_trat
## chisq df p
## Tratamiento_f 0.356 1 0.55
## GLOBAL 0.356 1 0.55
plot(proporcionalidad_trat)
# Proporcionalidad del modelo por edad
proporcionalidad_edad <- cox.zph(cox_edad)
proporcionalidad_edad
## chisq df p
## edad 0.601 1 0.44
## GLOBAL 0.601 1 0.44
plot(proporcionalidad_edad)
#Proporcionalidad del modelo por Apache II
proporcionalidad_apache <- cox.zph(cox_apache)
proporcionalidad_apache
## chisq df p
## apacheii 0.000639 1 0.98
## GLOBAL 0.000639 1 0.98
plot(proporcionalidad_apache)