Tesis Maestria 2026
Gabriela Carro
## New names:
## • `` -> `...44`table1=table(sob1$cicat)
table2=table(sob1$AMPU)
table3=table(sob1$fallecio)
table1##
## 0 1
## 189 216table2##
## 0 1
## 348 57table3##
## 0 1
## 391 14prop.table(table1)##
## 0 1
## 0.467 0.533prop.table(table2)##
## 0 1
## 0.859 0.141prop.table(table3)##
## 0 1
## 0.9654 0.0346#pasar variables categóricas a factor y chr a factor
names(sob1)<- tolower(names (sob1))
sob1 <- sob1 %>%
mutate(across(19:35, as.factor))
#crear cicat_factor
sob1$cicat_factor <- factor(sob1$cicat, levels = c(0, 1), labels = c( "No Cicatrizo", "Cicatrizo"))
table(sob1$cicat_factor)##
## No Cicatrizo Cicatrizo
## 189 216#crear ampu_factor
sob1$ampu_factor <- factor(sob1$ampu, levels = c(0, 1), labels = c( "No Amputado", "Amputado"))
table(sob1$ampu_factor)##
## No Amputado Amputado
## 348 57#crear fallecio_factor
sob1$fallecio_factor <- factor(sob1$fallecio, levels = c(0, 1), labels = c( "No Fallecido", "Fallecido"))
table(sob1$fallecio_factor)##
## No Fallecido Fallecido
## 391 14#paso variables chr a factor
sob1 <- sob1 %>% mutate_if(is.character, as.factor)#recategorizacion de variables
#ISQuEMIA
#isq_factor
sob1 <- sob1 %>%
mutate(isq_factor = factor(
if_else(isq %in% c(2, 3), "SI", "NO/leve")
))
#ISQuEMIA2
#isq_factor2
library(dplyr)
sob1 <- sob1 %>%
mutate(isq_factor2 = factor(
case_when(
isq == 0 ~ "NO",
isq == 1 ~ "LEVE",
isq == 2 ~ "MODERADA",
isq == 3 ~ "GRAVE",
TRUE ~ NA_character_
),
levels = c("NO", "LEVE", "MODERADA", "GRAVE")
))
table(sob1$isq)##
## 0 1 2 3
## 205 31 47 122table(sob1$isq_factor)##
## NO/leve SI
## 236 169table(sob1$isq_factor2)##
## NO LEVE MODERADA GRAVE
## 205 31 47 122#localizacion
sob1 <- sob1 %>%
mutate(loca_factor = factor(if_else(local %in% c(2, 3), "Met/Tar", "Fal")),
loca_factor = factor(loca_factor, levels = rev(levels(loca_factor))))
table(sob1$local)##
## 1 2 3
## 195 134 76table(sob1$loca_factor)##
## Met/Tar Fal
## 210 195#topog
sob1<- sob1 %>%
mutate(topo_factor = factor(if_else(topo %in% c(1,2), "dor/plan/lat", "Mas de dos")))
sob1 <- sob1 %>%
mutate(topo_factor = factor(if_else(topo %in% c(1, 2), "dor/plan/lat", "Mas de dos")),
topo_factor = factor(topo_factor, levels = rev(levels(topo_factor))))
table(sob1$topo)##
## 1 2 3
## 152 103 150table(sob1$topo_factor)##
## Mas de dos dor/plan/lat
## 150 255# zon
sob1 <- sob1 %>%
mutate(zon_factor = factor(if_else(zon %in% c(2, 3), "Dos o mas", "Una")), zon_factor = factor(zon_factor, levels = rev(levels(zon_factor))))
table(sob1$zon)##
## 1 2 3
## 242 89 74table(sob1$zon_factor)##
## Una Dos o mas
## 242 163#infe
sob1 <- sob1 %>%
mutate(infe_factor = case_when(
inf == 0 ~ "No",
inf %in% c(1, 2) ~ "Leve/Mod",
inf == 3 ~ "Grave",
TRUE ~ NA_character_ # conserva los NA originales
)) %>%
mutate(infe_factor = factor(infe_factor, levels = c("No", "Leve/Mod", "Grave")))
sob1$infe_factor <- factor (sob1$infe_factor, levels= c ("Grave","Leve/Mod","No"))
table(sob1$infe)## Warning: Unknown or uninitialised column: `infe`.## < table of extent 0 >table(sob1$infe_factor)##
## Grave Leve/Mod No
## 64 255 85#edema
sob1 <- sob1 %>%
mutate(ede_factor = case_when(
ede == 2 ~ "Unilateral",
ede %in% c(0, 1) ~ "Sin/local",
ede == 3 ~ "Bilateral",
TRUE ~ NA_character_ # conserva los NA originales
)) %>%
mutate(ede_factor = factor(ede_factor, levels = c("Sin/local", "Unilateral", "Bilateral")))
sob1$ede_factor <- factor (sob1$ede_factor, levels= c ("Bilateral","Unilateral","Sin/local"))
table(sob1$ede)##
## 0 1 2 3
## 142 98 124 41table(sob1$ede_factor)##
## Bilateral Unilateral Sin/local
## 41 124 240#neur
sob1 <- sob1 %>%
mutate(neur_factor = factor(
if_else(neu %in% c(1, 2, 3), "Leve/mod/charcot", "No"),
levels = c("No", "Leve/mod/charcot")
))
table(sob1$neur)## Warning: Unknown or uninitialised column: `neur`.## < table of extent 0 >table(sob1$neur_factor)##
## No Leve/mod/charcot
## 4 401#tisu (prof)
sob1 <- sob1 %>%
mutate(prof_factor = factor(if_else(prof %in% c(2, 3), "parcial/total", "Piel")))
table(sob1$prof)##
## 1 2 3
## 31 133 241table(sob1$prof_factor)##
## parcial/total Piel
## 374 31#area
sob1$area_factor <- with(sob1, ifelse(area == 1, "Menor a 10",
ifelse(area %in% c(1, 2), "10 a 40",
ifelse(area == 3, "Mas de 40", NA))))
sob1$area_factor <- factor(sob1$area_factor, levels = c("Menor a 10", "10 a 40", "Mas de 40"))
table(sob1$area)##
## 1 2 3
## 187 157 61table(sob1$area_factor)##
## Menor a 10 10 a 40 Mas de 40
## 187 157 61#cica(fase)
sob1 <- sob1 %>%
mutate(fase_factor = factor(if_else(fase %in% c(2, 3), "Granu/infla", "Epit")))
sob1$fase_factor <- factor (sob1$fase_factor, levels= c ("Granu/infla","Epit"))
table(sob1$fase)##
## 1 2 3
## 18 51 336table(sob1$fase_factor)##
## Granu/infla Epit
## 387 18hist(sob1$dias_out)
vis_miss(sob1)
#Tabla 1
library(tableone)
#Las variables van directamente sin el nombre del dataframe
catvars = c( "sexo","tbq", "iam", "hta", "dial", "local","topo", "inf", "ede", "neu", "prof", "area", "zon", "fase")
sob1 <- sob1 %>% mutate_if(is.character, as.factor)
#edad y enojo son continuas
vars = c("dias_ev","san_elian", "sexo","tbq", "iam", "hta", "dial", "local","topo", "inf", "ede", "neu", "prof", "area", "zon", "fase" )
tabla_1 <- CreateTableOne(vars = vars, strata = "isq_factor", factorVars = catvars, data = sob1)
print(tabla_1)## Stratified by isq_factor
## NO/leve SI p test
## n 236 169
## dias_ev (mean (SD)) 47.59 (118.96) 54.08 (104.37) 0.570
## san_elian (mean (SD)) 16.92 (3.97) 20.19 (3.41) <0.001
## sexo = M (%) 190 (80.5) 133 (78.7) 0.748
## tbq (%) 0.185
## 0 106 (45.1) 62 (36.9)
## 1 41 (17.4) 39 (23.2)
## 2 88 (37.4) 67 (39.9)
## iam = 1 (%) 32 (13.6) 51 (30.4) <0.001
## hta = 1 (%) 143 (60.6) 129 (76.8) 0.001
## dial = 1 (%) 10 ( 4.2) 9 ( 5.3) 0.785
## local (%) <0.001
## 1 101 (42.8) 94 (55.6)
## 2 99 (41.9) 35 (20.7)
## 3 36 (15.3) 40 (23.7)
## topo (%) 0.002
## 1 104 (44.1) 48 (28.4)
## 2 60 (25.4) 43 (25.4)
## 3 72 (30.5) 78 (46.2)
## inf (%) 0.700
## 0 54 (22.9) 31 (18.3)
## 1 19 ( 8.1) 13 ( 7.7)
## 2 127 (53.8) 96 (56.8)
## 3 35 (14.8) 29 (17.2)
## 4 1 ( 0.4) 0 ( 0.0)
## ede (%) 0.050
## 0 70 (29.7) 72 (42.6)
## 1 64 (27.1) 34 (20.1)
## 2 78 (33.1) 46 (27.2)
## 3 24 (10.2) 17 (10.1)
## neu (%) <0.001
## 0 1 ( 0.4) 3 ( 1.8)
## 1 6 ( 2.5) 19 (11.2)
## 2 200 (84.7) 140 (82.8)
## 3 29 (12.3) 7 ( 4.1)
## prof (%) 0.022
## 1 25 (10.6) 6 ( 3.6)
## 2 79 (33.5) 54 (32.0)
## 3 132 (55.9) 109 (64.5)
## area (%) 0.063
## 1 120 (50.8) 67 (39.6)
## 2 86 (36.4) 71 (42.0)
## 3 30 (12.7) 31 (18.3)
## zon (%) 0.026
## 1 154 (65.3) 88 (52.1)
## 2 46 (19.5) 43 (25.4)
## 3 36 (15.3) 38 (22.5)
## fase (%) 0.010
## 1 15 ( 6.4) 3 ( 1.8)
## 2 36 (15.3) 15 ( 8.9)
## 3 185 (78.4) 151 (89.3)kableone(tabla_1)| NO/leve | SI | p | test | |
|---|---|---|---|---|
| n | 236 | 169 | ||
| dias_ev (mean (SD)) | 47.59 (118.96) | 54.08 (104.37) | 0.570 | |
| san_elian (mean (SD)) | 16.92 (3.97) | 20.19 (3.41) | <0.001 | |
| sexo = M (%) | 190 (80.5) | 133 (78.7) | 0.748 | |
| tbq (%) | 0.185 | |||
| 0 | 106 (45.1) | 62 (36.9) | ||
| 1 | 41 (17.4) | 39 (23.2) | ||
| 2 | 88 (37.4) | 67 (39.9) | ||
| iam = 1 (%) | 32 (13.6) | 51 (30.4) | <0.001 | |
| hta = 1 (%) | 143 (60.6) | 129 (76.8) | 0.001 | |
| dial = 1 (%) | 10 ( 4.2) | 9 ( 5.3) | 0.785 | |
| local (%) | <0.001 | |||
| 1 | 101 (42.8) | 94 (55.6) | ||
| 2 | 99 (41.9) | 35 (20.7) | ||
| 3 | 36 (15.3) | 40 (23.7) | ||
| topo (%) | 0.002 | |||
| 1 | 104 (44.1) | 48 (28.4) | ||
| 2 | 60 (25.4) | 43 (25.4) | ||
| 3 | 72 (30.5) | 78 (46.2) | ||
| inf (%) | 0.700 | |||
| 0 | 54 (22.9) | 31 (18.3) | ||
| 1 | 19 ( 8.1) | 13 ( 7.7) | ||
| 2 | 127 (53.8) | 96 (56.8) | ||
| 3 | 35 (14.8) | 29 (17.2) | ||
| 4 | 1 ( 0.4) | 0 ( 0.0) | ||
| ede (%) | 0.050 | |||
| 0 | 70 (29.7) | 72 (42.6) | ||
| 1 | 64 (27.1) | 34 (20.1) | ||
| 2 | 78 (33.1) | 46 (27.2) | ||
| 3 | 24 (10.2) | 17 (10.1) | ||
| neu (%) | <0.001 | |||
| 0 | 1 ( 0.4) | 3 ( 1.8) | ||
| 1 | 6 ( 2.5) | 19 (11.2) | ||
| 2 | 200 (84.7) | 140 (82.8) | ||
| 3 | 29 (12.3) | 7 ( 4.1) | ||
| prof (%) | 0.022 | |||
| 1 | 25 (10.6) | 6 ( 3.6) | ||
| 2 | 79 (33.5) | 54 (32.0) | ||
| 3 | 132 (55.9) | 109 (64.5) | ||
| area (%) | 0.063 | |||
| 1 | 120 (50.8) | 67 (39.6) | ||
| 2 | 86 (36.4) | 71 (42.0) | ||
| 3 | 30 (12.7) | 31 (18.3) | ||
| zon (%) | 0.026 | |||
| 1 | 154 (65.3) | 88 (52.1) | ||
| 2 | 46 (19.5) | 43 (25.4) | ||
| 3 | 36 (15.3) | 38 (22.5) | ||
| fase (%) | 0.010 | |||
| 1 | 15 ( 6.4) | 3 ( 1.8) | ||
| 2 | 36 (15.3) | 15 ( 8.9) | ||
| 3 | 185 (78.4) | 151 (89.3) |
#variables categorizadas
#Las variables van directamente sin el nombre del dataframe
catvars2 = c( "sexo","tbq", "iam", "hta", "dial", "loca_factor","topo_factor", "infe_factor", "ede_factor", "neur_factor", "prof_factor", "area_factor", "zon_factor", "fase_factor")
sob1 <- sob1 %>% mutate_if(is.character, as.factor)
#edad y enojo son continuas
#vars2 = c("dias_ev","san_elian", "sexo","tbq", "iam", "hta", "dial", "loca_factor","topo_factor", "infe_factor", "ede_factor", "neur_factor", "prof_factor", "area_factor", "zon_factor", "fase_factor" )
#tabla_2 <- CreateTableOne(vars = vars2, strata = "isq_factor", factorVars = catvars2, data = sob1)
#print(tabla_2)
#kableone(tabla_2)#Elemento surv tiempo al evento cicatrización
tiempo_evento <- Surv(sob1$dias_out,sob1$cicat)
#evaluacion de la cantidad de eventos por tiempo de seguimiento a intervalos de a 30 dias
sob1 %>%
filter(cicat == 1) %>%
mutate(intervalo = cut(dias_out,
breaks = seq(0, max(dias_out), by = 30))) %>%
count(intervalo) %>%
ggplot(aes(x = intervalo, y = n)) +
geom_col(fill = "darkorange") +
labs(x = "Intervalos de tiempo (30 dias)",
y = "Eventos",
title = "Eventos por intervalo de seguimiento") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
#Modelo de cox
#Efecto de la isquemia sobre el tiempo a la cicatrización sin eventos competitivos ajustado por edad y edema
cox_model <- coxph(tiempo_evento ~ isq_factor+ edad + ede_factor , data = sob1)
summary(cox_model)## Call:
## coxph(formula = tiempo_evento ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 405, number of events= 216
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.74425 0.47509 0.17021 -4.37 0.000012 ***
## edad -0.01267 0.98741 0.00781 -1.62 0.10
## ede_factorUnilateral -0.06738 0.93484 0.25814 -0.26 0.79
## ede_factorSin/local 0.13589 1.14555 0.24062 0.56 0.57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.475 2.105 0.340 0.663
## edad 0.987 1.013 0.972 1.003
## ede_factorUnilateral 0.935 1.070 0.564 1.550
## ede_factorSin/local 1.146 0.873 0.715 1.836
##
## Concordance= 0.628 (se = 0.022 )
## Likelihood ratio test= 36.8 on 4 df, p=0.0000002
## Wald test = 33.6 on 4 df, p=0.0000009
## Score (logrank) test = 35.3 on 4 df, p=0.0000004tab_model(cox_model,
show.ci = TRUE, # Mostrar intervalo de confianza
show.se = TRUE, # Mostrar error estándar
transform = "exp", # Para mostrar HR en lugar de log(HR)
dv.labels = "Modelo de Cox - Datos drugs")| Modelo de Cox - Datos drugs | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| isq factorSI | 0.48 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.99 | 0.01 | 0.00 – Inf | 0.105 |
| ede factor [Unilateral] | 0.93 | 0.24 | 0.00 – Inf | 0.794 |
| ede factorSin/local | 1.15 | 0.28 | 0.00 – Inf | 0.572 |
| Observations | 405 | |||
| R2 Nagelkerke | 0.087 | |||
confint(cox_model)## 2.5 % 97.5 %
## isq_factorSI -1.078 -0.41065
## edad -0.028 0.00263
## ede_factorUnilateral -0.573 0.43857
## ede_factorSin/local -0.336 0.60749#Supuesto de proporcionalidad
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph <- cox.zph(cox_model)
test_ph## chisq df p
## isq_factor 2.038 1 0.15
## edad 0.397 1 0.53
## ede_factor 1.610 2 0.45
## GLOBAL 3.921 4 0.42#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
#gráficos loglog
##loglog isq_factor
ggsurvplot(
survfit(tiempo_evento ~ isq_factor, data = sob1),
fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons"
)## 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 every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.## Ignoring unknown labels:
## • fill : "Strata"
## • linetype : "1"
## Ignoring unknown labels:
## • fill : "Strata"
## • linetype : "1"
#ggsurvplot(
# survfit(tiempo_evento ~ infe_factor, data = sob1),
# fun = "cloglog",
# xlab = "Tiempo (dias)",
# ylab = "Log(-Log(S(t)))",
# palette = "simpsons"
#)
ggsurvplot(
survfit(tiempo_evento ~ ede_factor, data = sob1),
fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons"
)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.## Ignoring unknown labels:
## • fill : "Strata"
## • linetype : "1"
## Ignoring unknown labels:
## • fill : "Strata"
## • linetype : "1"
#Gráfico de Kaplan Meier
# Graficar las curvas de Kaplan-Meier estratificadas por isq
km_global <- survfit(tiempo_evento ~ 1, data = sob1)#ponemos 1 para marcar que es crudo, sin estratificar por ninguna variable
#6b. Ajustar el modelo Kaplan-Meier estratificado por "beck_cat"
km_fit <- survfit(tiempo_evento ~ isq_factor,data = sob1)
ggsurvplot(
km_fit,
data = sob1,
fun = "event", # 👈 ESTO hace la curva ascendente
risk.table = TRUE,
pval = TRUE,
conf.int = TRUE,
cumcensor = TRUE,
ggtheme = theme_minimal(),
palette = "simpsons",
xlab = "Tiempo de seguimiento (dias)",
ylab = "Probabilidad acumulada del evento"
)## Ignoring unknown labels:
## • linetype : "1"
## Ignoring unknown labels:
## • linetype : "1"
## Ignoring unknown labels:
## • linetype : "1"
## Ignoring unknown labels:
## • linetype : "1"
## Ignoring unknown labels:
## • colour : "Strata"
## Ignoring unknown labels:
## • colour : "Strata"
#logrank
logrank <- survdiff(tiempo_evento ~ isq_factor, data = sob1)
logrank## Call:
## survdiff(formula = tiempo_evento ~ isq_factor, data = sob1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## isq_factor=NO/leve 236 157 117 13.7 30.6
## isq_factor=SI 169 59 99 16.2 30.6
##
## Chisq= 30.6 on 1 degrees of freedom, p= 0.00000003#Residuos de Martingale y Deviance: Adecuacion del modelo
# Martingale y Deviance
ggcoxdiagnostics(cox_model, type = "martingale", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggcoxdiagnostics(cox_model, type = "deviance", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
library(psych)## Warning: package 'psych' was built under R version 4.4.3##
## Adjuntando el paquete: 'psych'## The following objects are masked from 'package:scales':
##
## alpha, rescale## The following object is masked from 'package:ggsurvfit':
##
## %+%## The following object is masked from 'package:Hmisc':
##
## describe## The following objects are masked from 'package:ggplot2':
##
## %+%, alphatable(sob1$estado)##
## ampu AMPU cerro CERRO fallecio FALLECIO persiste persisTE
## 1 56 6 210 8 6 10 1
## perSISTE PERsiste PERSISTE
## 15 7 85sob1 <- sob1 %>%
mutate(
estado = tolower(trimws(estado))
)
sob1 <- sob1 %>%
mutate(
estado_unif = case_when(
estado %in% c("ampu") ~ "amputacion",
estado %in% c("cerro") ~ "cerro_herida",
estado %in% c("murio", "fallecio") ~ "fallecio",
estado %in% c("persiste") ~ "persiste",
estado %in% c("perdido") ~ "perdido_seguimiento",
TRUE ~ NA_character_
)
)
table(sob1$estado_unif)##
## amputacion cerro_herida fallecio persiste
## 57 216 14 118describeBy(
x = sob1$edad,
group = sob1$estado_unif,
mat = TRUE
)## item group1 vars n mean sd median trimmed mad min max range
## X11 1 amputacion 1 57 62.9 8.38 63.0 62.7 7.41 42 88 46
## X12 2 cerro_herida 1 216 55.1 10.90 55.5 55.1 11.12 22 89 67
## X13 3 fallecio 1 14 62.5 19.74 64.0 63.1 19.27 30 88 58
## X14 4 persiste 1 118 58.5 10.01 57.5 58.5 9.64 32 89 57
## skew kurtosis se
## X11 0.3347 0.632 1.110
## X12 -0.0242 0.378 0.741
## X13 -0.2246 -1.416 5.277
## X14 0.0962 -0.043 0.922#Discriminación y calibración
#C de harrell
summary(cox_model)$concordance## C se(C)
## 0.6283 0.0219#bootstrap
#Indice D de Somer (Dxy), Metricas Q D U, slope y R2 con paquete `rms`
dd <- datadist(sob1); options(datadist = "dd")## Warning in datadist(sob1): ...44 is constantcox_metrics <- cph(tiempo_evento ~ edad + isq_factor + ede_factor,
data = sob1, x = TRUE, y = TRUE, surv = TRUE)
# Validamos métricas con bootstrap
invisible(capture.output({ #esta linea es para evitar que se imprima el output del bootstrap)
cox_rms <- validate(cox_metrics, method = "boot", B = 200)
}))
# Revisamos las metricas validadas (en la columna index.corrected)
cox_rms## index.orig training test optimism index.corrected n
## Dxy 0.2566 0.2655 0.2440 0.0215 0.2351 200
## R2 0.0872 0.0968 0.0808 0.0160 0.0712 200
## Slope 1.0000 1.0000 0.9200 0.0800 0.9200 200
## D 0.0167 0.0188 0.0153 0.0035 0.0132 200
## U -0.0009 -0.0009 0.0007 -0.0016 0.0007 200
## Q 0.0176 0.0197 0.0147 0.0051 0.0125 200
## g 0.5073 0.5380 0.4817 0.0562 0.4510 200round(cox_rms[,5],2)## Dxy R2 Slope D U Q g
## 0.24 0.07 0.92 0.01 0.00 0.01 0.45# Calibración. Evaluación Gráfica
invisible(capture.output({cal <- calibrate(cox_metrics, method = "boot", B = 200, u = 120)
}))
plot(cal,
subtitles = FALSE,
lwd = 2,
lty = 2,
main = "Grafico de Calibracion",
cex.main = 1.2,
cex.axis = 0.7,
cex.lab = 0.9,
xlab = "Probabilidad Predicha de Cicatrizacion",
ylab = "Cicatrizacion Observada")
La curva azul está sistemáticamente por encima de la diagonal.
Eso sugiere que el modelo está sobreestimando la probabilidad de cicatrización.
Porque ignora que:
muchos pacientes mueren
muchos se amputan
Y por lo tanto, inflan la probabilidad teórica de cicatrizar.
#evento cica tomando evento competitivo am y muerte
sob1 <- sob1 %>%
mutate(
evento_cica = case_when(
estado_unif == "cerro_herida" ~ 1,
estado_unif %in% c( "fallecio", "amputacion") ~ 2,
estado_unif %in% c("perdido_seguimiento", "persiste") ~ 0,
TRUE ~ NA_real_
)
)
table(sob1$evento_cica)##
## 0 1 2
## 118 216 71#Opcion con objeto creado
#si evento es 1, es porque se amputo o persiste
#si evento es 2, es porque se murió o se perdió
tiempo_cicat <- Surv(sob1$dias_out, sob1$evento_cica==1)
tiempo_compet <- Surv(sob1$dias_out, sob1$evento_cica==2)
cox_cicat <- coxph(tiempo_cicat ~
isq_factor + edad+
ede_factor , data = sob1)
summary(cox_cicat)## Call:
## coxph(formula = tiempo_cicat ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 405, number of events= 216
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.74425 0.47509 0.17021 -4.37 0.000012 ***
## edad -0.01267 0.98741 0.00781 -1.62 0.10
## ede_factorUnilateral -0.06738 0.93484 0.25814 -0.26 0.79
## ede_factorSin/local 0.13589 1.14555 0.24062 0.56 0.57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.475 2.105 0.340 0.663
## edad 0.987 1.013 0.972 1.003
## ede_factorUnilateral 0.935 1.070 0.564 1.550
## ede_factorSin/local 1.146 0.873 0.715 1.836
##
## Concordance= 0.628 (se = 0.022 )
## Likelihood ratio test= 36.8 on 4 df, p=0.0000002
## Wald test = 33.6 on 4 df, p=0.0000009
## Score (logrank) test = 35.3 on 4 df, p=0.0000004cox_compet <- coxph(tiempo_compet ~
isq_factor + edad+
ede_factor , data = sob1)
summary(cox_compet)## Call:
## coxph(formula = tiempo_compet ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 405, number of events= 71
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.2783 3.5907 0.3007 4.25 0.000021 ***
## edad 0.0270 1.0274 0.0124 2.18 0.029 *
## ede_factorUnilateral 0.3001 1.3500 0.4248 0.71 0.480
## ede_factorSin/local -0.2945 0.7449 0.4145 -0.71 0.477
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 3.591 0.279 1.992 6.47
## edad 1.027 0.973 1.003 1.05
## ede_factorUnilateral 1.350 0.741 0.587 3.10
## ede_factorSin/local 0.745 1.342 0.331 1.68
##
## Concordance= 0.709 (se = 0.031 )
## Likelihood ratio test= 42.8 on 4 df, p=0.00000001
## Wald test = 36.8 on 4 df, p=0.0000002
## Score (logrank) test = 42.2 on 4 df, p=0.00000002#CIF
cifivhx <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor, data=sob1)
cifivhx## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 82 0.501 0.035 0.432, 0.567
## NO/leve 200 22 0.742 0.034 0.669, 0.802
## NO/leve 300 3 0.893 0.029 0.819, 0.938
## SI 100 68 0.186 0.032 0.128, 0.252
## SI 200 23 0.398 0.044 0.311, 0.483
## SI 300 5 0.495 0.053 0.386, 0.594## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 82 0.076 0.018 0.046, 0.116
## NO/leve 200 22 0.076 0.018 0.046, 0.116
## NO/leve 300 3 0.076 0.018 0.046, 0.116
## SI 100 68 0.283 0.036 0.216, 0.355
## SI 200 23 0.341 0.039 0.266, 0.417
## SI 300 5 0.355 0.041 0.276, 0.434## • Tests## outcome statistic df p.value
## 1 55.0 1.00 <0.001
## 2 40.0 1.00 <0.001#Gráfico CIF
library(tidycmprsk)
library(ggsurvfit)
ggcuminc(
cifivhx,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO/leve" = "#56B4E9", "SI" = "#E69F00")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c(
"Evento de interes (cicatrizacion)",
"Evento competitivo (muerte / amputacion mayor)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_step()`).
#Modelo de causa específica Cuando hacés predicción con este modelo, lo que obtenés es:
“Probabilidad de cicatrizar si los eventos competitivos fueran censura”
Es decir:
Trata muerte y amputación como si fueran pérdidas no informativas.
La probabilidad predicha no es la CIF real.
Entonces tu gráfico 1 está calibrando:
Riesgo derivado del modelo causa-específica (no la incidencia acumulada real).
#modelo de causa específica
modelo_CSC <- CSC(
formula = Hist(dias_out, evento_cica) ~ isq_factor+edad + ede_factor, data = sob1,
cause = 1 # Aclara que CHD es el evento de interés
)
# Ver resultados. Ojo! No funciona si hacés summary(modelo)
# Vamos a tener el modelo de interés y debajo el modelo competitivo
modelo_CSC## CSC(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor, data = sob1, cause = 1)
##
## Right-censored response of a competing.risks model
##
## No.Observations: 405
##
## Pattern:
##
## Cause event right.censored
## 1 216 0
## 2 71 0
## unknown 0 118
##
##
## ----------> Cause: 1
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 405, number of events= 216
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.74425 0.47509 0.17021 -4.37 0.000012 ***
## edad -0.01267 0.98741 0.00781 -1.62 0.10
## ede_factorUnilateral -0.06738 0.93484 0.25814 -0.26 0.79
## ede_factorSin/local 0.13589 1.14555 0.24062 0.56 0.57
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.475 2.105 0.340 0.663
## edad 0.987 1.013 0.972 1.003
## ede_factorUnilateral 0.935 1.070 0.564 1.550
## ede_factorSin/local 1.146 0.873 0.715 1.836
##
## Concordance= 0.628 (se = 0.022 )
## Likelihood ratio test= 36.8 on 4 df, p=0.0000002
## Wald test = 33.6 on 4 df, p=0.0000009
## Score (logrank) test = 35.3 on 4 df, p=0.0000004
##
##
##
## ----------> Cause: 2
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 405, number of events= 71
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.2783 3.5907 0.3007 4.25 0.000021 ***
## edad 0.0270 1.0274 0.0124 2.18 0.029 *
## ede_factorUnilateral 0.3001 1.3500 0.4248 0.71 0.480
## ede_factorSin/local -0.2945 0.7449 0.4145 -0.71 0.477
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 3.591 0.279 1.992 6.47
## edad 1.027 0.973 1.003 1.05
## ede_factorUnilateral 1.350 0.741 0.587 3.10
## ede_factorSin/local 0.745 1.342 0.331 1.68
##
## Concordance= 0.709 (se = 0.031 )
## Likelihood ratio test= 42.8 on 4 df, p=0.00000001
## Wald test = 36.8 on 4 df, p=0.0000002
## Score (logrank) test = 42.2 on 4 df, p=0.00000002#Modelo de Fine and Gray Probabilidad real acumulada de cicatrización considerando muerte y amputación como eventos competitivos.
Tu segundo gráfico está calibrando:
Probabilidad real acumulada observada vs predicha (CIF).
#modelo de Fine and gray
modelo_fg <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg##
## Right-censored response of a competing.risks model
##
## No.Observations: 405
##
## Pattern:
##
## Cause event right.censored
## 1 216 0
## 2 71 0
## unknown 0 118
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.995 0.370 0.1699 -5.858 0.0000000047
## edad -0.015 0.985 0.0077 -1.952 0.0510000000
## ede_factorUnilateral -0.167 0.846 0.2521 -0.663 0.5100000000
## ede_factorSin/local 0.176 1.193 0.2349 0.751 0.4500000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.370 2.706 0.265 0.516
## edad 0.985 1.015 0.970 1.000
## ede_factorUnilateral 0.846 1.182 0.516 1.387
## ede_factorSin/local 1.193 0.838 0.753 1.890
##
## Num. cases = 405
## Pseudo Log-likelihood = -1111
## Pseudo likelihood ratio test = 67.1 on 4 df,
##
## Convergence: TRUE#Comparación Causa específica y FyG
#score comparado con ipa
score_comp <-Score(
list("FG model" = modelo_fg, "CSC model" = modelo_CSC),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1,
times = c(17:200),
cause = 1,
metrics ="auc",
summary = "IPA"
)IPA = Index of Prediction Accuracy (Índice de Precisión Predictiva)
Es una medida que evalúa qué tan bien un modelo predice el riesgo en el tiempo, y se basa en el Brier score. Brier score: Mide el error cuadrático medio entre:lo que el modelo predice (probabilidad de evento a tiempo t) y lo que realmente ocurrió (valor esperado 0 a 0.12) IPA = 0 → el modelo no mejora al modelo nulo
IPA > 0 → mejora predictiva
IPA = 1 → predicción perfecta
IPA < 0 → peor que no usar modelo 😅
En la práctica:
0.05–0.10 ya puede ser clínicamente relevante
0.15 suele ser buen desempeño
#Extraigo los AUC por tiempo y modelo
score <- as.data.frame(score_comp$AUC$score)
# Gráfico
ggplot(score, aes(x = times, y = AUC, color = model)) +
geom_line(linewidth = 1.2) +
geom_point(size = 3) +
geom_hline(yintercept = 0.5, linetype = "dashed", alpha = 0.2) +
geom_hline(yintercept = 0.7, linetype = "dotted", alpha = 0.2) +
scale_color_manual(values = c("FG model" = "#E74C3C", "CSC model" = "#3498DB")) +
scale_y_continuous(limits = c(0.45, 0.85), breaks = seq(0.5, 0.8, 0.1)) +
labs(
title = "Discriminacion de Modelos en el Tiempo",
subtitle = "Linea punteada = Azar (0.5), Linea punteada = Adecuada (0.7)",
x = "Dias de seguimiento",
y = "AUC (Area bajo la curva ROC)",
color = "Modelo"
) +
theme_minimal() +
theme(legend.position = "bottom")
#Interpretamos que año a año la capacidad discriminativa de los modelos se mantiene alta, tanto en CSC como en FGR# Calibracion global
s <- Score(list("FGR" = modelo_fg),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(s, models = "FGR" )
La curva negra está mucho más cerca de la diagonal.
Eso es esperable porque:
Fine–Gray modela directamente la CIF
Está alineado con lo que realmente ocurre clínicamente
Además te muestra:
AUC = 70
Brier ≈ 0.20
Eso es desempeño predictivo razonable.