Tesis 17 febrero
Carro Gabriela
2026-02-18
#CARGA DE LA BASE DE DATOS
#carga de datos
library(readxl)
sob1 <- read_excel("C:/Users/Administrador/Downloads/TESI.xlsx")## New names:
## • `` -> `...44`View(sob1)names(sob1)<- tolower(names (sob1)) #PASA A MINUSCULA LOS NOMBRES
#PROPORCION CICATRIZADOS
table1=table(sob1$cicat)
table1##
## 0 1
## 183 223prop.table(table1)##
## 0 1
## 0.451 0.549#PROPORCION AMPUTADOS
table2=table(sob1$ampu)
table2##
## 0 1
## 349 57prop.table(table2)##
## 0 1
## 0.86 0.14#PROPORCION FALLECIDOS
table3=table(sob1$fallecio)
table3##
## 0 1
## 392 14prop.table(table3)##
## 0 1
## 0.9655 0.0345#tabla por sexo
table4=table(sob1$sexo)
table4##
## F M
## 82 324prop.table(table4)##
## F M
## 0.202 0.798#promedio de edad
describe(sob1$edad)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 406 57.4 11.1 57 57.4 10.4 22 89 67 0 0.35 0.55#dias de seguimiento
describe(sob1$dias_out)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 406 96.2 75.5 77.5 86.4 64.5 1 362 361 1.15 0.94 3.75#Descripción de los eventos
#pasar variables categóricas a factor
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
## 183 223#paso variables chr a factor
sob1 <- sob1 %>% mutate_if(is.character, as.factor)
#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
## 349 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
## 392 14#variable unificada
sob1 <- sob1 %>%
mutate(
estado = tolower(trimws(estado))
)
table(sob1$estado)##
## ampu cerro fallecio persiste
## 57 223 14 112sob1 <- 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 223 14 112##Intervalo de confianza de los eventos
#intervalos de confianza
tabla_estado <- sob1 %>%
dplyr::filter(!is.na(estado_unif)) %>%
dplyr::group_by(estado_unif) %>%
dplyr::summarise(
n = n()
) %>%
dplyr::mutate(
N_total = sum(n),
prop = n / N_total,
se = sqrt(prop * (1 - prop) / N_total),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
tabla_estado %>%
dplyr::mutate(
porcentaje = round(prop * 100, 1),
IC95 = paste0(
round(li * 100, 1), "-",
round(ls * 100, 1), "%"
)
)## # A tibble: 4 × 9
## estado_unif n N_total prop se li ls porcentaje IC95
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 amputacion 57 406 0.140 0.0172 0.107 0.174 14 10.7-17.4%
## 2 cerro_herida 223 406 0.549 0.0247 0.501 0.598 54.9 50.1-59.8%
## 3 fallecio 14 406 0.0345 0.00906 0.0167 0.0522 3.4 1.7-5.2%
## 4 persiste 112 406 0.276 0.0222 0.232 0.319 27.6 23.2-31.9%tabla_estado## # A tibble: 4 × 7
## estado_unif n N_total prop se li ls
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 amputacion 57 406 0.140 0.0172 0.107 0.174
## 2 cerro_herida 223 406 0.549 0.0247 0.501 0.598
## 3 fallecio 14 406 0.0345 0.00906 0.0167 0.0522
## 4 persiste 112 406 0.276 0.0222 0.232 0.319#Días de seguimiento
#dias de seguimiento
#en pacientes censurados
#dias de seguimiento
describe(sob1$dias_out)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 406 96.2 75.5 77.5 86.4 64.5 1 362 361 1.15 0.94 3.75describeBy(
x = sob1$dias_out,
group = sob1$estado_unif,
mat = TRUE
)## item group1 vars n mean sd median trimmed mad min max range
## X11 1 amputacion 1 57 39.6 39.9 22 32.8 16.3 1 211 210
## X12 2 cerro_herida 1 223 99.7 70.5 82 89.7 57.8 9 330 321
## X13 3 fallecio 1 14 43.1 40.4 24 39.9 25.2 5 120 115
## X14 4 persiste 1 112 124.7 83.7 114 117.6 83.8 6 362 356
## skew kurtosis se
## X11 1.982 4.5398 5.28
## X12 1.234 1.0772 4.72
## X13 0.661 -1.2525 10.80
## X14 0.718 0.0353 7.91tapply(sob1$dias_out,
sob1$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 14 48
##
## $cerro_herida
## 25% 75%
## 46 128
##
## $fallecio
## 25% 75%
## 11.2 69.8
##
## $persiste
## 25% 75%
## 63.8 175.8#ISQuEMIA
#isq_factor
sob1 <- sob1 %>%
mutate(isq_factor = factor(
if_else(isq %in% c(2, 3), "SI", "NO/leve")
))
#dias seguimiento en subgrupos
sob1_per_isq <- subset(sob1, isq_factor == "SI")
sob1_per_noisq <- subset(sob1, isq_factor == "NO/leve")
#no isquemicos persistente seguimiento
describe(sob1_per_noisq$dias_out)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 237 91.1 68.2 77 82.7 62.3 5 337 332 1.11 0.88 4.43describeBy(
x = sob1_per_noisq$dias_out,
group = sob1_per_noisq$estado_unif,
mat = TRUE
)## item group1 vars n mean sd median trimmed mad min max range
## X11 1 amputacion 1 11 26.5 26.1 19.0 20.8 11.9 7 98 91
## X12 2 cerro_herida 1 162 90.6 63.9 71.5 80.9 52.6 15 301 286
## X13 3 fallecio 1 6 25.8 23.1 17.5 25.8 15.6 5 66 61
## X14 4 persiste 1 58 111.7 76.6 112.5 107.2 86.7 6 337 331
## skew kurtosis se
## X11 1.789 2.182 7.88
## X12 1.313 1.319 5.02
## X13 0.708 -1.298 9.44
## X14 0.532 -0.056 10.06tapply(sob1_per_noisq$dias_out,
sob1_per_noisq$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 12.5 26.0
##
## $cerro_herida
## 25% 75%
## 42 116
##
## $fallecio
## 25% 75%
## 10.8 34.8
##
## $persiste
## 25% 75%
## 50.8 167.0#pacientes isquémicos
# isquemicos persistente seguimiento
describe(sob1_per_isq$dias_out)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 169 103 84.3 82 92.2 80.1 1 362 361 1.07 0.52 6.49describeBy(
x = sob1_per_isq$dias_out,
group = sob1_per_isq$estado_unif,
mat = TRUE
)## item group1 vars n mean sd median trimmed mad min max range
## X11 1 amputacion 1 46 42.7 42.1 24.5 35.9 20.8 1 211 210
## X12 2 cerro_herida 1 61 124.1 81.2 108.0 114.9 72.6 9 330 321
## X13 3 fallecio 1 8 56.1 46.9 50.0 56.1 61.5 6 120 114
## X14 4 persiste 1 54 138.7 89.3 116.5 130.2 79.3 7 362 355
## skew kurtosis se
## X11 1.847 3.8060 6.21
## X12 0.896 0.0851 10.40
## X13 0.134 -1.9918 16.58
## X14 0.743 -0.3552 12.15tapply(sob1_per_isq$dias_out,
sob1_per_isq$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 14.2 59.0
##
## $cerro_herida
## 25% 75%
## 63 160
##
## $fallecio
## 25% 75%
## 11.8 99.0
##
## $persiste
## 25% 75%
## 67.2 199.2#comparar dias de seguimiento en pacientes con y sin isquemia cuyo estado sea "persiste"
# Filtrar solo pacientes con estado "persiste"
sob_persiste <- subset(sob1, estado_unif == "persiste")
# Aplicar Wilcoxon
wilcox.test(dias_out ~ isq_factor, data = sob_persiste)##
## Wilcoxon rank sum test with continuity correction
##
## data: dias_out by isq_factor
## W = 1330, p-value = 0.2
## alternative hypothesis: true location shift is not equal to 0#Densidad de incidencia y razón de densidades
#cicatrizacion totaldensidad de incidencia
sob1$evento_cierre <- ifelse(sob1$estado_unif == "cerro_herida", 1, 0)
eventos <- sum(sob1$evento_cierre, na.rm = TRUE) #nro eventos
tiempo1_total <- sum(sob1$dias_out, na.rm = TRUE) #tiempo persona total
tasa1_dia <- eventos / tiempo1_total
tasa1_dia## [1] 0.00571tiempo1_meses <- tiempo1_total / 30.44
tasa1_mensual <- eventos / tiempo1_meses
tasa1_mensual## [1] 0.174tasa1_1000 <- tasa1_dia * 1000
tasa1_1000## [1] 5.71#cicatrizacion isq tiempo por 1000 personas por dia
sob1_per_isq$evento_cierre <- ifelse(sob1_per_isq$estado_unif == "cerro_herida", 1, 0)
eventos <- sum(sob1_per_isq$evento_cierre, na.rm = TRUE) #nro eventos
tiempo2_total <- sum(sob1_per_isq$dias_out, na.rm = TRUE) #tiempo persona total
tasa2_dia <- eventos / tiempo2_total
tasa2_dia## [1] 0.00349tiempo2_meses <- tiempo2_total / 30.44
tasa2_mensual <- eventos / tiempo2_meses
tasa2_mensual## [1] 0.106tasa2_1000 <- tasa2_dia * 1000
tasa2_1000## [1] 3.49#cicatrizacion no isq
sob1_per_noisq$evento_cierre <- ifelse(sob1_per_noisq$estado_unif == "cerro_herida", 1, 0)
eventos <- sum(sob1_per_noisq$evento_cierre, na.rm = TRUE) #nro eventos
tiempo3_total <- sum(sob1_per_noisq$dias_out, na.rm = TRUE) #tiempo persona total
tasa3_dia <- eventos / tiempo3_total
tasa3_dia## [1] 0.0075tiempo3_meses <- tiempo3_total / 30.44
tasa3_mensual <- eventos / tiempo3_meses
tasa3_mensual## [1] 0.228tasa3_1000 <- tasa3_dia * 1000
tasa3_1000## [1] 7.5# densidad de incidencia amputacion total
sob1$evento_ampu <- ifelse(sob1$estado_unif == "amputacion", 1, 0)
eventos <- sum(sob1$evento_ampu, na.rm = TRUE) #nro eventos
tiempo1_total <- sum(sob1$dias_out, na.rm = TRUE) #tiempo persona total
tasa1_dia <- eventos / tiempo1_total
tasa1_dia## [1] 0.00146tiempo1_meses <- tiempo1_total / 30.44
tasa1_mensual <- eventos / tiempo1_meses
tasa1_mensual## [1] 0.0444tasa1_1000 <- tasa1_dia * 1000
tasa1_1000## [1] 1.46#amputacuin isq
sob1_per_isq$evento_ampu <- ifelse(sob1_per_isq$estado_unif == "amputacion", 1, 0)
eventos <- sum(sob1_per_isq$evento_ampu, na.rm = TRUE) #nro eventos
tiempo2_total <- sum(sob1_per_isq$dias_out, na.rm = TRUE) #tiempo persona total
tasa2_dia <- eventos / tiempo2_total
tasa2_dia## [1] 0.00263tiempo2_meses <- tiempo2_total / 30.44
tasa2_mensual <- eventos / tiempo2_meses
tasa2_mensual## [1] 0.0801tasa2_1000 <- tasa2_dia * 1000
tasa2_1000## [1] 2.63#amputacion no isq
sob1_per_noisq$evento_ampu <- ifelse(sob1_per_noisq$estado_unif == "amputacion", 1, 0)
eventos <- sum(sob1_per_noisq$evento_ampu, na.rm = TRUE) #nro eventos
tiempo3_total <- sum(sob1_per_noisq$dias_out, na.rm = TRUE) #tiempo persona total
tasa3_dia <- eventos / tiempo3_total
tasa3_dia## [1] 0.000509tiempo3_meses <- tiempo3_total / 30.44
tasa3_mensual <- eventos / tiempo3_meses
tasa3_mensual## [1] 0.0155tasa3_1000 <- tasa3_dia * 1000
tasa3_1000## [1] 0.509#evento muerte densidad de incidencia
#cicatrizacion totaldensidad de incidencia
sob1$evento_muerte <- ifelse(sob1$estado_unif == "fallecio", 1, 0)
eventos <- sum(sob1$evento_muerte, na.rm = TRUE) #nro eventos
tiempo1_total <- sum(sob1$dias_out, na.rm = TRUE) #tiempo persona total
tasa1_dia <- eventos / tiempo1_total
tasa1_dia## [1] 0.000358tiempo1_meses <- tiempo1_total / 30.44
tasa1_mensual <- eventos / tiempo1_meses
tasa1_mensual## [1] 0.0109tasa1_1000 <- tasa1_dia * 1000
tasa1_1000## [1] 0.358#muerte isq tiempo por 1000 personas por dia. Tasa de incidencia de muerte por 1000 persons por dia
sob1_per_isq$evento_muerte <- ifelse(sob1_per_isq$estado_unif == "fallecio", 1, 0)
eventos <- sum(sob1_per_isq$evento_muerte, na.rm = TRUE) #nro eventos
tiempo2_total <- sum(sob1_per_isq$dias_out, na.rm = TRUE) #tiempo persona total
tasa2_dia <- eventos / tiempo2_total
tasa2_dia## [1] 0.000458tiempo2_meses <- tiempo2_total / 30.44
tasa2_mensual <- eventos / tiempo2_meses
tasa2_mensual## [1] 0.0139tasa2_1000 <- tasa2_dia * 1000
tasa2_1000## [1] 0.458#muerte no isq
sob1_per_noisq$evento_muerte <- ifelse(sob1_per_noisq$estado_unif == "fallecio", 1, 0)
eventos <- sum(sob1_per_noisq$evento_muerte, na.rm = TRUE) #nro eventos
tiempo3_total <- sum(sob1_per_noisq$dias_out, na.rm = TRUE) #tiempo persona total
tasa3_dia <- eventos / tiempo3_total
tasa3_dia## [1] 0.000278tiempo3_meses <- tiempo3_total / 30.44
tasa3_mensual <- eventos / tiempo3_meses
tasa3_mensual## [1] 0.00846tasa3_1000 <- tasa3_dia * 1000
tasa3_1000## [1] 0.278#COMPORTAMIENTO DE LAS VARIABLES
#comportamiento de las variables
#LOCAl
library(dplyr)
tabla_prop <- sob1 %>%
dplyr::group_by(local) %>%
dplyr::summarise(
n = sum(!is.na(cicat)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
library(ggplot2)
ggplot(tabla_prop, aes(x = local, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Localizacion de la lesion",
y = "Proporcion de cicatrizacion (IC95%)"
) +
theme_minimal(base_size = 14)## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#localizacion
sob1 <- sob1 %>%
mutate(loca_factor = factor(if_else(local %in% c(2, 3), "Met/Tar", "Fal")))
table(sob1$local)##
## 1 2 3
## 196 134 76table(sob1$loca_factor)##
## Fal Met/Tar
## 196 210#TOPOGRAFIA
tabla_prop <- sob1 %>%
dplyr::group_by(topo) %>%
dplyr::summarise(
n = sum(!is.na(topo)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
ggplot(tabla_prop, aes(x = topo, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Topografia de la lesion",
y = "Proporcion de cicatrizacion (IC95%)"
) +
theme_minimal(base_size = 14)
#topog
sob1<- sob1 %>%
mutate(topo_factor = factor(if_else(topo %in% c(1), "dor/plan", "lat/Mas de dos")))
table(sob1$topo)##
## 1 2 3
## 152 104 150table(sob1$topo_factor)##
## dor/plan lat/Mas de dos
## 152 254#Nro de zonas
tabla_prop <- sob1 %>%
dplyr::group_by(zon) %>%
dplyr::summarise(
n = sum(!is.na(zon)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
ggplot(tabla_prop, aes(x = zon, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Numero de Zonas de la lesion",
y = "Proporcion de cicatrizacion (IC95%)"
) +
theme_minimal(base_size = 14)
# 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 90 74table(sob1$zon_factor)##
## Una Dos o mas
## 242 164#INFECCION
tabla_prop <- sob1 %>%
dplyr::group_by(inf) %>%
dplyr::summarise(
n = sum(!is.na(inf)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
ggplot(tabla_prop, aes(x = inf, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Grado de infeccion de la lesion",
y = "Proporcion de cicatrizacion (IC95%)"
) +
theme_minimal(base_size = 14)
#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)##
## No Leve/Mod Grave
## 86 255 65#EDEMA
tabla_prop <- sob1 %>%
dplyr::group_by(ede) %>%
dplyr::summarise(
n = sum(!is.na(ede)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
ggplot(tabla_prop, aes(x = ede, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Edema de la lesion",
y = "Proporcion de cicatrizacion (IC95%)"
) +
theme_minimal(base_size = 14)
#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
## 143 98 124 41table(sob1$ede_factor)##
## Sin/local Unilateral Bilateral
## 241 124 41#NEURO
tabla_prop <- sob1 %>%
dplyr::group_by(neu) %>%
dplyr::summarise(
n = sum(!is.na(neu)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
ggplot(tabla_prop, aes(x = neu, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Neuropatia de la lesion",
y = "Proporcion de cicatrizacion (IC95%)"
) +
theme_minimal(base_size = 14)
#neur
sob1 <- sob1 %>%
mutate(neur_factor = factor(
if_else(neu %in% c(2, 3), "Mod/charcot", "No/leve"),
levels = c("No/leve", "Mod/charcot")
))
table(sob1$neu)##
## 0 1 2 3
## 4 25 341 36table(sob1$neur_factor)##
## No/leve Mod/charcot
## 29 377#PROF
tabla_prop <- sob1 %>%
dplyr::group_by(prof) %>%
dplyr::summarise(
n = sum(!is.na(prof)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
ggplot(tabla_prop, aes(x = prof, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Profundidad de la lesion",
y = "Proporcion de cicatrizacion (IC95%)"
) +
theme_minimal(base_size = 14)
#tisu (prof)
sob1 <- sob1 %>%
mutate(prof_factor = factor(if_else(prof %in% c(1), "Piel", "parcial/total")))
sob1 <- sob1 %>% mutate(prof_factor = factor(prof_factor, levels = c("Piel", "parcial/total")))
table(sob1$prof)##
## 1 2 3
## 32 133 241table(sob1$prof_factor)##
## Piel parcial/total
## 32 374#AREA
tabla_prop <- sob1 %>%
dplyr::group_by(area) %>%
dplyr::summarise(
n = sum(!is.na(area)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
ggplot(tabla_prop, aes(x = area, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Area de la lesion",
y = "Proporcion de cicatrizacion (IC95%)"
) +
theme_minimal(base_size = 14)
#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
## 188 157 61table(sob1$area_factor)##
## Menor a 10 10 a 40 Mas de 40
## 188 157 61#FASE de cicat
tabla_prop <- sob1 %>%
dplyr::group_by(fase) %>%
dplyr::summarise(
n = sum(!is.na(fase)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
ggplot(tabla_prop, aes(x = fase, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Area de la lesion",
y = "Proporcion de cicatrizacion (IC95%)"
) +
theme_minimal(base_size = 14)
#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 ("Epit","Granu/infla"))
table(sob1$fase)##
## 1 2 3
## 19 51 336table(sob1$fase_factor)##
## Epit Granu/infla
## 19 387#isq
#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
## 206 31 47 122table(sob1$isq_factor)##
## NO/leve SI
## 237 169table(sob1$isq_factor2)##
## NO LEVE MODERADA GRAVE
## 206 31 47 122#quiero saber el porcentaje de pacientes isquémicos de mi muestra
x <- sum(sob1$isq_factor2 %in% c("LEVE","MODERADA","GRAVE"), na.rm = TRUE)
n <- sum(!is.na(sob1$isq_factor2))
prop.test(x, n)##
## 1-sample proportions test with continuity correction
##
## data: x out of n, null probability 0.5
## X-squared = 0.06, df = 1, p-value = 0.8
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.443 0.542
## sample estimates:
## p
## 0.493#porcentaje de isquemia moderada grave
xi <- sum(sob1$isq_factor2 %in% c("MODERADA","GRAVE"), na.rm = TRUE)
ni <- sum(!is.na(sob1$isq_factor2))
prop.test(xi, ni)##
## 1-sample proportions test with continuity correction
##
## data: xi out of ni, null probability 0.5
## X-squared = 11, df = 1, p-value = 0.0009
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.368 0.466
## sample estimates:
## p
## 0.416#porcentaje de isquemia no/leve
xin <- sum(sob1$isq_factor2 %in% c("NO","LEVE"), na.rm = TRUE)
nin <- sum(!is.na(sob1$isq_factor2))
prop.test(xin, nin)##
## 1-sample proportions test with continuity correction
##
## data: xin out of nin, null probability 0.5
## X-squared = 11, df = 1, p-value = 0.0009
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.534 0.632
## sample estimates:
## p
## 0.584#porcentaje de cicatrizacion en no isq/leve
table(sob1_per_noisq$cicat_factor)##
## No Cicatrizo Cicatrizo
## 75 162xci <- sum(sob1_per_noisq$cicat_factor %in% c("Cicatrizo"), na.rm = TRUE)
nci <- sum(!is.na(sob1_per_noisq$cicat_factor))
prop.test(xci, nci)##
## 1-sample proportions test with continuity correction
##
## data: xci out of nci, null probability 0.5
## X-squared = 31, df = 1, p-value = 0.00000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.620 0.741
## sample estimates:
## p
## 0.684#porcentaje de amputacion en no isq/leve
table(sob1_per_noisq$ampu_factor)##
## No Amputado Amputado
## 226 11xa <- sum(sob1_per_noisq$ampu_factor %in% c("Amputado"), na.rm = TRUE)
na <- sum(!is.na(sob1_per_noisq$ampu_factor))
prop.test(xa, na)##
## 1-sample proportions test with continuity correction
##
## data: xa out of na, null probability 0.5
## X-squared = 193, df = 1, p-value <0.0000000000000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.0246 0.0838
## sample estimates:
## p
## 0.0464#porcentaje de muerte en no isq/leve
table(sob1_per_noisq$fallecio_factor)##
## No Fallecido Fallecido
## 231 6xf <- sum(sob1_per_noisq$fallecio_factor %in% c("Fallecido"), na.rm = TRUE)
nf <- sum(!is.na(sob1_per_noisq$fallecio_factor))
prop.test(xf, nf)##
## 1-sample proportions test with continuity correction
##
## data: xf out of nf, null probability 0.5
## X-squared = 212, df = 1, p-value <0.0000000000000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.0103 0.0569
## sample estimates:
## p
## 0.0253#porcentaje de cicatrizacion en isquemicos
table(sob1_per_isq$cicat_factor)##
## No Cicatrizo Cicatrizo
## 108 61xcii <- sum(sob1_per_isq$cicat_factor %in% c("Cicatrizo"), na.rm = TRUE)
ncii <- sum(!is.na(sob1_per_isq$cicat_factor))
prop.test(xcii, ncii)##
## 1-sample proportions test with continuity correction
##
## data: xcii out of ncii, null probability 0.5
## X-squared = 13, df = 1, p-value = 0.0004
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.290 0.439
## sample estimates:
## p
## 0.361#porcentaje de amputacion en isquemicos
table(sob1_per_isq$ampu_factor)##
## No Amputado Amputado
## 123 46xai <- sum(sob1_per_isq$ampu_factor %in% c("Amputado"), na.rm = TRUE)
nai <- sum(!is.na(sob1_per_isq$ampu_factor))
prop.test(xai, nai)##
## 1-sample proportions test with continuity correction
##
## data: xai out of nai, null probability 0.5
## X-squared = 34, df = 1, p-value = 0.000000005
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.208 0.347
## sample estimates:
## p
## 0.272#porcentaje de muerte en isquemicos
table(sob1_per_isq$fallecio_factor)##
## No Fallecido Fallecido
## 161 8xfi <- sum(sob1_per_isq$fallecio_factor %in% c("Fallecido"), na.rm = TRUE)
nfi<- sum(!is.na(sob1_per_isq$fallecio_factor))
prop.test(xfi, nfi)##
## 1-sample proportions test with continuity correction
##
## data: xfi out of nfi, null probability 0.5
## X-squared = 137, df = 1, p-value <0.0000000000000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.0222 0.0944
## sample estimates:
## p
## 0.0473#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")
#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)
#kableone(tabla_1)
#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")
vars2 = c("edad","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", "am_cont", "ant_pie" )
tabla_2 <- CreateTableOne(
vars = vars2,
strata = "isq_factor",
factorVars = catvars2,
data = sob1,
addOverall = TRUE
)
print(tabla_2)## Stratified by isq_factor
## Overall NO/leve SI
## n 406 237 169
## edad (mean (SD)) 57.41 (11.09) 53.34 (9.97) 63.12 (10.03)
## dias_ev (mean (SD)) 50.26 (112.85) 47.52 (118.71) 54.08 (104.37)
## san_elian (mean (SD)) 18.26 (4.09) 16.88 (3.99) 20.19 (3.41)
## sexo = M (%) 324 (79.8) 191 (80.6) 133 (78.7)
## tbq (%)
## 0 169 (41.8) 107 (45.3) 62 (36.9)
## 1 80 (19.8) 41 (17.4) 39 (23.2)
## 2 155 (38.4) 88 (37.3) 67 (39.9)
## iam = 1 (%) 83 (20.5) 32 (13.6) 51 (30.4)
## hta = 1 (%) 273 (67.4) 144 (60.8) 129 (76.8)
## dial = 1 (%) 19 ( 4.7) 10 ( 4.2) 9 ( 5.3)
## loca_factor = Met/Tar (%) 210 (51.7) 135 (57.0) 75 (44.4)
## topo_factor = lat/Mas de dos (%) 254 (62.6) 133 (56.1) 121 (71.6)
## infe_factor (%)
## No 86 (21.2) 55 (23.2) 31 (18.3)
## Leve/Mod 255 (62.8) 146 (61.6) 109 (64.5)
## Grave 65 (16.0) 36 (15.2) 29 (17.2)
## ede_factor (%)
## Sin/local 241 (59.4) 135 (57.0) 106 (62.7)
## Unilateral 124 (30.5) 78 (32.9) 46 (27.2)
## Bilateral 41 (10.1) 24 (10.1) 17 (10.1)
## neur_factor = Mod/charcot (%) 377 (92.9) 230 (97.0) 147 (87.0)
## prof_factor = parcial/total (%) 374 (92.1) 211 (89.0) 163 (96.4)
## area_factor (%)
## Menor a 10 188 (46.3) 121 (51.1) 67 (39.6)
## 10 a 40 157 (38.7) 86 (36.3) 71 (42.0)
## Mas de 40 61 (15.0) 30 (12.7) 31 (18.3)
## zon_factor = Dos o mas (%) 164 (40.4) 83 (35.0) 81 (47.9)
## fase_factor = Granu/infla (%) 387 (95.3) 221 (93.2) 166 (98.2)
## am_cont = 1 (%) 23 ( 5.7) 5 ( 2.1) 18 (10.7)
## ant_pie (%)
## 0 142 (35.0) 69 (29.1) 73 (43.2)
## 1 263 (64.8) 167 (70.5) 96 (56.8)
## 2 1 ( 0.2) 1 ( 0.4) 0 ( 0.0)
## Stratified by isq_factor
## p test
## n
## edad (mean (SD)) <0.001
## dias_ev (mean (SD)) 0.565
## san_elian (mean (SD)) <0.001
## sexo = M (%) 0.732
## tbq (%) 0.171
## 0
## 1
## 2
## iam = 1 (%) <0.001
## hta = 1 (%) 0.001
## dial = 1 (%) 0.778
## loca_factor = Met/Tar (%) 0.016
## topo_factor = lat/Mas de dos (%) 0.002
## infe_factor (%) 0.479
## No
## Leve/Mod
## Grave
## ede_factor (%) 0.450
## Sin/local
## Unilateral
## Bilateral
## neur_factor = Mod/charcot (%) <0.001
## prof_factor = parcial/total (%) 0.011
## area_factor (%) 0.057
## Menor a 10
## 10 a 40
## Mas de 40
## zon_factor = Dos o mas (%) 0.012
## fase_factor = Granu/infla (%) 0.036
## am_cont = 1 (%) 0.001
## ant_pie (%) 0.010
## 0
## 1
## 2print(tabla_2,
nonnormal = c("dias_ev", "san_elian"),
showAllLevels = TRUE,
quote = FALSE,
noSpaces = TRUE)## Stratified by isq_factor
## level Overall
## n 406
## edad (mean (SD)) 57.41 (11.09)
## dias_ev (median [IQR]) 20.00 [10.00, 45.00]
## san_elian (median [IQR]) 18.00 [15.00, 21.00]
## sexo (%) F 82 (20.2)
## M 324 (79.8)
## tbq (%) 0 169 (41.8)
## 1 80 (19.8)
## 2 155 (38.4)
## iam (%) 0 321 (79.5)
## 1 83 (20.5)
## hta (%) 0 132 (32.6)
## 1 273 (67.4)
## dial (%) 0 387 (95.3)
## 1 19 (4.7)
## loca_factor (%) Fal 196 (48.3)
## Met/Tar 210 (51.7)
## topo_factor (%) dor/plan 152 (37.4)
## lat/Mas de dos 254 (62.6)
## infe_factor (%) No 86 (21.2)
## Leve/Mod 255 (62.8)
## Grave 65 (16.0)
## ede_factor (%) Sin/local 241 (59.4)
## Unilateral 124 (30.5)
## Bilateral 41 (10.1)
## neur_factor (%) No/leve 29 (7.1)
## Mod/charcot 377 (92.9)
## prof_factor (%) Piel 32 (7.9)
## parcial/total 374 (92.1)
## area_factor (%) Menor a 10 188 (46.3)
## 10 a 40 157 (38.7)
## Mas de 40 61 (15.0)
## zon_factor (%) Una 242 (59.6)
## Dos o mas 164 (40.4)
## fase_factor (%) Epit 19 (4.7)
## Granu/infla 387 (95.3)
## am_cont (%) 0 383 (94.3)
## 1 23 (5.7)
## ant_pie (%) 0 142 (35.0)
## 1 263 (64.8)
## 2 1 (0.2)
## Stratified by isq_factor
## NO/leve SI p
## n 237 169
## edad (mean (SD)) 53.34 (9.97) 63.12 (10.03) <0.001
## dias_ev (median [IQR]) 20.00 [7.50, 30.00] 25.00 [14.00, 60.00] 0.020
## san_elian (median [IQR]) 17.00 [14.00, 20.00] 20.00 [18.00, 23.00] <0.001
## sexo (%) 46 (19.4) 36 (21.3) 0.732
## 191 (80.6) 133 (78.7)
## tbq (%) 107 (45.3) 62 (36.9) 0.171
## 41 (17.4) 39 (23.2)
## 88 (37.3) 67 (39.9)
## iam (%) 204 (86.4) 117 (69.6) <0.001
## 32 (13.6) 51 (30.4)
## hta (%) 93 (39.2) 39 (23.2) 0.001
## 144 (60.8) 129 (76.8)
## dial (%) 227 (95.8) 160 (94.7) 0.778
## 10 (4.2) 9 (5.3)
## loca_factor (%) 102 (43.0) 94 (55.6) 0.016
## 135 (57.0) 75 (44.4)
## topo_factor (%) 104 (43.9) 48 (28.4) 0.002
## 133 (56.1) 121 (71.6)
## infe_factor (%) 55 (23.2) 31 (18.3) 0.479
## 146 (61.6) 109 (64.5)
## 36 (15.2) 29 (17.2)
## ede_factor (%) 135 (57.0) 106 (62.7) 0.450
## 78 (32.9) 46 (27.2)
## 24 (10.1) 17 (10.1)
## neur_factor (%) 7 (3.0) 22 (13.0) <0.001
## 230 (97.0) 147 (87.0)
## prof_factor (%) 26 (11.0) 6 (3.6) 0.011
## 211 (89.0) 163 (96.4)
## area_factor (%) 121 (51.1) 67 (39.6) 0.057
## 86 (36.3) 71 (42.0)
## 30 (12.7) 31 (18.3)
## zon_factor (%) 154 (65.0) 88 (52.1) 0.012
## 83 (35.0) 81 (47.9)
## fase_factor (%) 16 (6.8) 3 (1.8) 0.036
## 221 (93.2) 166 (98.2)
## am_cont (%) 232 (97.9) 151 (89.3) 0.001
## 5 (2.1) 18 (10.7)
## ant_pie (%) 69 (29.1) 73 (43.2) 0.010
## 167 (70.5) 96 (56.8)
## 1 (0.4) 0 (0.0)
## Stratified by isq_factor
## test
## n
## edad (mean (SD))
## dias_ev (median [IQR]) nonnorm
## san_elian (median [IQR]) nonnorm
## sexo (%)
##
## tbq (%)
##
##
## iam (%)
##
## hta (%)
##
## dial (%)
##
## loca_factor (%)
##
## topo_factor (%)
##
## infe_factor (%)
##
##
## ede_factor (%)
##
##
## neur_factor (%)
##
## prof_factor (%)
##
## area_factor (%)
##
##
## zon_factor (%)
##
## fase_factor (%)
##
## am_cont (%)
##
## ant_pie (%)
##
## kableone(tabla_2)| Overall | NO/leve | SI | p | test | |
|---|---|---|---|---|---|
| n | 406 | 237 | 169 | ||
| edad (mean (SD)) | 57.41 (11.09) | 53.34 (9.97) | 63.12 (10.03) | <0.001 | |
| dias_ev (mean (SD)) | 50.26 (112.85) | 47.52 (118.71) | 54.08 (104.37) | 0.565 | |
| san_elian (mean (SD)) | 18.26 (4.09) | 16.88 (3.99) | 20.19 (3.41) | <0.001 | |
| sexo = M (%) | 324 (79.8) | 191 (80.6) | 133 (78.7) | 0.732 | |
| tbq (%) | 0.171 | ||||
| 0 | 169 (41.8) | 107 (45.3) | 62 (36.9) | ||
| 1 | 80 (19.8) | 41 (17.4) | 39 (23.2) | ||
| 2 | 155 (38.4) | 88 (37.3) | 67 (39.9) | ||
| iam = 1 (%) | 83 (20.5) | 32 (13.6) | 51 (30.4) | <0.001 | |
| hta = 1 (%) | 273 (67.4) | 144 (60.8) | 129 (76.8) | 0.001 | |
| dial = 1 (%) | 19 ( 4.7) | 10 ( 4.2) | 9 ( 5.3) | 0.778 | |
| loca_factor = Met/Tar (%) | 210 (51.7) | 135 (57.0) | 75 (44.4) | 0.016 | |
| topo_factor = lat/Mas de dos (%) | 254 (62.6) | 133 (56.1) | 121 (71.6) | 0.002 | |
| infe_factor (%) | 0.479 | ||||
| No | 86 (21.2) | 55 (23.2) | 31 (18.3) | ||
| Leve/Mod | 255 (62.8) | 146 (61.6) | 109 (64.5) | ||
| Grave | 65 (16.0) | 36 (15.2) | 29 (17.2) | ||
| ede_factor (%) | 0.450 | ||||
| Sin/local | 241 (59.4) | 135 (57.0) | 106 (62.7) | ||
| Unilateral | 124 (30.5) | 78 (32.9) | 46 (27.2) | ||
| Bilateral | 41 (10.1) | 24 (10.1) | 17 (10.1) | ||
| neur_factor = Mod/charcot (%) | 377 (92.9) | 230 (97.0) | 147 (87.0) | <0.001 | |
| prof_factor = parcial/total (%) | 374 (92.1) | 211 (89.0) | 163 (96.4) | 0.011 | |
| area_factor (%) | 0.057 | ||||
| Menor a 10 | 188 (46.3) | 121 (51.1) | 67 (39.6) | ||
| 10 a 40 | 157 (38.7) | 86 (36.3) | 71 (42.0) | ||
| Mas de 40 | 61 (15.0) | 30 (12.7) | 31 (18.3) | ||
| zon_factor = Dos o mas (%) | 164 (40.4) | 83 (35.0) | 81 (47.9) | 0.012 | |
| fase_factor = Granu/infla (%) | 387 (95.3) | 221 (93.2) | 166 (98.2) | 0.036 | |
| am_cont = 1 (%) | 23 ( 5.7) | 5 ( 2.1) | 18 (10.7) | 0.001 | |
| ant_pie (%) | 0.010 | ||||
| 0 | 142 (35.0) | 69 (29.1) | 73 (43.2) | ||
| 1 | 263 (64.8) | 167 (70.5) | 96 (56.8) | ||
| 2 | 1 ( 0.2) | 1 ( 0.4) | 0 ( 0.0) |
#Creacion del elemento surv tiempo a la 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))
##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 237 162 120 14.4 32.3
## isq_factor=SI 169 61 103 16.9 32.3
##
## Chisq= 32.3 on 1 degrees of freedom, p= 0.00000001summary(km_fit)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## isq_factor=NO/leve 237 237 237 162 120 6.55 91 84
## isq_factor=SI 169 169 169 61 195 11.94 170 153
## 0.95UCL
## isq_factor=NO/leve 109
## isq_factor=SI 275#mediana
km_fit$table[, c("median","0.95LCL","0.95UCL")]## NULLsurv_median(km_fit)## strata median lower upper
## 1 isq_factor=NO/leve 91 84 109
## 2 isq_factor=SI 170 153 275#Modelo de cox. Se realiza para evaluar los supuestos
#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= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71978 0.48686 0.16765 -4.29 0.000018 ***
## edad -0.01610 0.98402 0.00765 -2.11 0.035 *
## ede_factorUnilateral -0.20647 0.81345 0.15253 -1.35 0.176
## ede_factorBilateral -0.15043 0.86034 0.23500 -0.64 0.522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.676
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.860 1.16 0.543 1.364
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007tab_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.49 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.035 |
| ede factor [Unilateral] | 0.81 | 0.12 | 0.00 – Inf | 0.176 |
| ede factorBilateral | 0.86 | 0.20 | 0.00 – Inf | 0.522 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.095 | |||
confint(cox_model)## 2.5 % 97.5 %
## isq_factorSI -1.0484 -0.39118
## edad -0.0311 -0.00111
## ede_factorUnilateral -0.5054 0.09249
## ede_factorBilateral -0.6110 0.31016##Supuesto de proporcionalidad
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph <- cox.zph(cox_model)
test_ph## chisq df p
## isq_factor 1.8077 1 0.18
## edad 0.0327 1 0.86
## ede_factor 1.4240 2 0.49
## GLOBAL 3.6145 4 0.46#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: 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"
##Residuos de Martingale y deviance para la adecuación 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'
##Discriminación y calibración del COX
#C de harrell
summary(cox_model)$concordance## C se(C)
## 0.6301 0.0213#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.2601 0.2692 0.2470 0.0223 0.2379 200
## R2 0.0954 0.1059 0.0891 0.0168 0.0786 200
## Slope 1.0000 1.0000 0.9191 0.0809 0.9191 200
## D 0.0177 0.0201 0.0165 0.0036 0.0142 200
## U -0.0009 -0.0009 0.0006 -0.0015 0.0006 200
## Q 0.0186 0.0210 0.0159 0.0051 0.0136 200
## g 0.5290 0.5619 0.5058 0.0561 0.4729 200round(cox_rms[,5],2)## Dxy R2 Slope D U Q g
## 0.24 0.08 0.92 0.01 0.00 0.01 0.47##Calibración. Evaluación gráfica
# 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.
#Fine and Gray ##Creación de variable con evento competitivo
#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
## 112 223 71#chi2 para comparación de cicatrizacion, muerte y AM en isq y no isq
chisq.test(sob1$isq_factor,sob1$evento_cica)##
## Pearson's Chi-squared test
##
## data: sob1$isq_factor and sob1$evento_cica
## X-squared = 55, df = 2, p-value = 0.000000000001chisq.test(sob1$isq_factor,sob1$evento_ampu)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: sob1$isq_factor and sob1$evento_ampu
## X-squared = 40, df = 1, p-value = 0.0000000003chisq.test(sob1$isq_factor,sob1$evento_muerte)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: sob1$isq_factor and sob1$evento_muerte
## X-squared = 0.9, df = 1, p-value = 0.4#Creación objeto SURV
#modelo de causa específica
#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= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71978 0.48686 0.16765 -4.29 0.000018 ***
## edad -0.01610 0.98402 0.00765 -2.11 0.035 *
## ede_factorUnilateral -0.20647 0.81345 0.15253 -1.35 0.176
## ede_factorBilateral -0.15043 0.86034 0.23500 -0.64 0.522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.676
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.860 1.16 0.543 1.364
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007cox_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= 406, number of events= 71
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.2845 3.6128 0.3009 4.27 0.00002 ***
## edad 0.0268 1.0272 0.0124 2.16 0.030 *
## ede_factorUnilateral 0.5849 1.7948 0.2557 2.29 0.022 *
## ede_factorBilateral 0.2923 1.3395 0.4145 0.71 0.481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 3.61 0.277 2.003 6.52
## edad 1.03 0.974 1.003 1.05
## ede_factorUnilateral 1.79 0.557 1.087 2.96
## ede_factorBilateral 1.34 0.747 0.594 3.02
##
## Concordance= 0.71 (se = 0.031 )
## Likelihood ratio test= 43 on 4 df, p=0.00000001
## Wald test = 36.9 on 4 df, p=0.0000002
## Score (logrank) test = 42.4 on 4 df, p=0.00000001##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 para comparar los IPA y el AUC
#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: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## ----------> Cause: 1
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71978 0.48686 0.16765 -4.29 0.000018 ***
## edad -0.01610 0.98402 0.00765 -2.11 0.035 *
## ede_factorUnilateral -0.20647 0.81345 0.15253 -1.35 0.176
## ede_factorBilateral -0.15043 0.86034 0.23500 -0.64 0.522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.676
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.860 1.16 0.543 1.364
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007
##
##
##
## ----------> Cause: 2
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 406, number of events= 71
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.2845 3.6128 0.3009 4.27 0.00002 ***
## edad 0.0268 1.0272 0.0124 2.16 0.030 *
## ede_factorUnilateral 0.5849 1.7948 0.2557 2.29 0.022 *
## ede_factorBilateral 0.2923 1.3395 0.4145 0.71 0.481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 3.61 0.277 2.003 6.52
## edad 1.03 0.974 1.003 1.05
## ede_factorUnilateral 1.79 0.557 1.087 2.96
## ede_factorBilateral 1.34 0.747 0.594 3.02
##
## Concordance= 0.71 (se = 0.031 )
## Likelihood ratio test= 43 on 4 df, p=0.00000001
## Wald test = 36.9 on 4 df, p=0.0000002
## Score (logrank) test = 42.4 on 4 df, p=0.00000001#Modelo de Fine and Gray Probabilidad real acumulada de cicatrización considerando muerte y amputación como eventos competitivos. 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: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## 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.9719 0.378 0.16781 -5.791 0.000000007
## edad -0.0172 0.983 0.00771 -2.236 0.025000000
## ede_factorUnilateral -0.3437 0.709 0.14962 -2.297 0.022000000
## ede_factorBilateral -0.1976 0.821 0.23140 -0.854 0.390000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.378 2.64 0.272 0.526
## edad 0.983 1.02 0.968 0.998
## ede_factorUnilateral 0.709 1.41 0.529 0.951
## ede_factorBilateral 0.821 1.22 0.521 1.292
##
## Num. cases = 406
## Pseudo Log-likelihood = -1151
## Pseudo likelihood ratio test = 70.8 on 4 df,
##
## Convergence: TRUE#Gray test
library(cmprsk)
cuminc_obj <- cmprsk::cuminc(
ftime = sob1$dias_out,
fstatus = sob1$evento_cica,
group = sob1$isq_factor
)
cuminc_obj## Tests:
## stat pv df
## 1 56.8 0.0000000000000484 1
## 2 40.4 0.0000000002049831 1
## Estimates and Variances:
## $est
## 100 200 300
## NO/leve 1 0.4989 0.7424 0.8903
## SI 1 0.1830 0.3941 0.4955
## NO/leve 2 0.0751 0.0751 0.0751
## SI 2 0.2822 0.3372 0.3500
##
## $var
## 100 200 300
## NO/leve 1 0.001167 0.001115 0.000864
## SI 1 0.000998 0.001915 0.002577
## NO/leve 2 0.000309 0.000309 0.000309
## SI 2 0.001267 0.001483 0.001592##Performance predictiva con Score
# Evaluar performance predictiva con Score. AUC
score_fg <- Score(
object = list("FG model" = modelo_fg),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1,
times = c(17:200),
metrics = "auc",
summary = "IPA",
cens.model = "km",
cause = 1
)
score_fg$AUC$score## model times AUC se lower upper
## <fctr> <int> <num> <num> <num> <num>
## 1: FG model 17 0.503 0.1179 0.272 0.734
## 2: FG model 18 0.530 0.1064 0.322 0.739
## 3: FG model 19 0.530 0.1064 0.322 0.739
## 4: FG model 20 0.530 0.1064 0.322 0.739
## 5: FG model 21 0.574 0.0900 0.397 0.750
## ---
## 180: FG model 196 0.723 0.0313 0.662 0.785
## 181: FG model 197 0.720 0.0316 0.658 0.782
## 182: FG model 198 0.720 0.0316 0.658 0.782
## 183: FG model 199 0.725 0.0317 0.663 0.787
## 184: FG model 200 0.724 0.0320 0.661 0.787##comparacion causa especifica 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"
)##Índice de Precisión Predictiva
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" )
## Gráfico de ALTMAN. Probabilidad acumulada de incidencia de
cicatrización a 180 días
#altman
edades <- seq(30, 85, by = 1)
new_no <- data.frame(
isq_factor = factor(rep(levels(sob1$isq_factor)[1], length(edades)),
levels = levels(sob1$isq_factor)),
edad = edades,
ede_factor = factor(rep("Sin/local", length(edades)),
levels = levels(sob1$ede_factor))
)
new_si <- new_no
new_si$isq_factor <- factor(rep(levels(sob1$isq_factor)[2], length(edades)),
levels = levels(sob1$isq_factor))
pred_no <- predictRisk(modelo_fg,
newdata = new_no,
times = 180,
cause = 1)
pred_si <- predictRisk(modelo_fg,
newdata = new_si,
times = 180,
cause = 1)
plot(edades, pred_no,
type = "l",
lwd = 3,
xlab = "Edad (años)",
ylab = "Probabilidad de cicatrizacion a 180 dias",
ylim = c(0, max(pred_no, pred_si)),
col = "black")
lines(edades, pred_si,
lwd = 3,
lty = 2,
col = "black")
legend("bottomleft",
legend = c("Sin isquemia/isq leve", "Con isquemia moderada/grave"),
lwd = 3,
lty = c(1,2),
bty = "n")
##ANALISIS DE SENSIBILIDAD
#excluyendo del análisis a los pacientes con herida persistente seguidos 30 días o menos
sob1_menos30 <- sob1 %>%
filter(!(estado_unif == "persiste" & dias_out < 31))
tiempo_evento_s <- Surv(sob1_menos30$dias_out,sob1_menos30$cicat)
#modelo de cox
#Efecto de la isquemia sobre el tiempo a la cicatrización sin eventos competitivos ajustado por edad y edema
cox_model_s <- coxph(tiempo_evento_s ~ isq_factor+ edad + ede_factor , data = sob1_menos30)
summary(cox_model_s)## Call:
## coxph(formula = tiempo_evento_s ~ isq_factor + edad + ede_factor,
## data = sob1_menos30)
##
## n= 390, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71920 0.48714 0.16763 -4.29 0.000018 ***
## edad -0.01613 0.98400 0.00764 -2.11 0.035 *
## ede_factorUnilateral -0.20657 0.81337 0.15254 -1.35 0.176
## ede_factorBilateral -0.15199 0.85899 0.23499 -0.65 0.518
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.677
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.859 1.16 0.542 1.362
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007tab_model(cox_model_s,
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.49 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.035 |
| ede factor [Unilateral] | 0.81 | 0.12 | 0.00 – Inf | 0.176 |
| ede factorBilateral | 0.86 | 0.20 | 0.00 – Inf | 0.518 |
| Observations | 390 | |||
| R2 Nagelkerke | 0.099 | |||
confint(cox_model_s)## 2.5 % 97.5 %
## isq_factorSI -1.0478 -0.39065
## edad -0.0311 -0.00115
## ede_factorUnilateral -0.5055 0.09240
## ede_factorBilateral -0.6126 0.30859# Martingale y Deviance
ggcoxdiagnostics(cox_model_s, type = "martingale", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggcoxdiagnostics(cox_model_s, type = "deviance", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph_s <- cox.zph(cox_model_s)
test_ph_s## chisq df p
## isq_factor 1.7970 1 0.18
## edad 0.0335 1 0.85
## ede_factor 1.4034 2 0.50
## GLOBAL 3.5787 4 0.47#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph_s)## 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.
#modelo de Fine and gray
modelo_fg_s <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor
, data = sob1_menos30,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg_s##
## Right-censored response of a competing.risks model
##
## No.Observations: 390
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 96
##
##
## 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_menos30, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9733 0.378 0.16798 -5.79 0.0000000069
## edad -0.0173 0.983 0.00771 -2.24 0.0250000000
## ede_factorUnilateral -0.3477 0.706 0.14997 -2.32 0.0200000000
## ede_factorBilateral -0.1990 0.820 0.23122 -0.86 0.3900000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.378 2.65 0.272 0.525
## edad 0.983 1.02 0.968 0.998
## ede_factorUnilateral 0.706 1.42 0.526 0.948
## ede_factorBilateral 0.820 1.22 0.521 1.289
##
## Num. cases = 390
## Pseudo Log-likelihood = -1151
## Pseudo likelihood ratio test = 71.2 on 4 df,
##
## Convergence: TRUE# Calibracion global
ss <- Score(list("FGR" = modelo_fg_s),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1_menos30,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(ss, models = "FGR" )
#Considerando a los pacientes con herida persistente seguidos 30 días o menoscomo TODOS CICATRIZADOS sob1_sens
sob1_sens <- sob1
sob1_sens <- sob1 %>%
mutate(
estado_unif = ifelse(
estado_unif == "persiste" & dias_out < 31,
"cerro_herida",
estado_unif
)
)
table(sob1_sens$estado_unif)##
## amputacion cerro_herida fallecio persiste
## 57 239 14 96tiempo_evento_s1 <- Surv(sob1_sens$dias_out,sob1_sens$cicat)
#modelo de cox
#Efecto de la isquemia sobre el tiempo a la cicatrización sin eventos competitivos ajustado por edad y edema
cox_model_s1 <- coxph(tiempo_evento_s1 ~ isq_factor+ edad + ede_factor , data = sob1_sens)
summary(cox_model_s1)## Call:
## coxph(formula = tiempo_evento_s1 ~ isq_factor + edad + ede_factor,
## data = sob1_sens)
##
## n= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71978 0.48686 0.16765 -4.29 0.000018 ***
## edad -0.01610 0.98402 0.00765 -2.11 0.035 *
## ede_factorUnilateral -0.20647 0.81345 0.15253 -1.35 0.176
## ede_factorBilateral -0.15043 0.86034 0.23500 -0.64 0.522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.676
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.860 1.16 0.543 1.364
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007tab_model(cox_model_s1,
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.49 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.035 |
| ede factor [Unilateral] | 0.81 | 0.12 | 0.00 – Inf | 0.176 |
| ede factorBilateral | 0.86 | 0.20 | 0.00 – Inf | 0.522 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.095 | |||
confint(cox_model_s1)## 2.5 % 97.5 %
## isq_factorSI -1.0484 -0.39118
## edad -0.0311 -0.00111
## ede_factorUnilateral -0.5054 0.09249
## ede_factorBilateral -0.6110 0.31016# Martingale y Deviance
ggcoxdiagnostics(cox_model_s1, type = "martingale", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggcoxdiagnostics(cox_model_s1, type = "deviance", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph_s1 <- cox.zph(cox_model_s1)
test_ph_s1## chisq df p
## isq_factor 1.8077 1 0.18
## edad 0.0327 1 0.86
## ede_factor 1.4240 2 0.49
## GLOBAL 3.6145 4 0.46#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph_s1)## 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.
#evento cica: agrega como cicatrizados a los perdidos
sob1_sens <- sob1_sens %>%
mutate(
evento_cica = ifelse(
sob1$estado_unif == "persiste" & dias_out < 31,
1,
evento_cica
)
)
#modelo de Fine and gray
modelo_fg_s1 <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor
, data = sob1_sens,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg_s1##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 239 0
## 2 71 0
## unknown 0 96
##
##
## 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_sens, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9534 0.385 0.16094 -5.924 0.0000000031
## edad -0.0168 0.983 0.00718 -2.339 0.0190000000
## ede_factorUnilateral -0.3064 0.736 0.14465 -2.118 0.0340000000
## ede_factorBilateral -0.2083 0.812 0.22283 -0.935 0.3500000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.385 2.59 0.281 0.528
## edad 0.983 1.02 0.970 0.997
## ede_factorUnilateral 0.736 1.36 0.554 0.977
## ede_factorBilateral 0.812 1.23 0.525 1.257
##
## Num. cases = 406
## Pseudo Log-likelihood = -1247
## Pseudo likelihood ratio test = 72.2 on 4 df,
##
## Convergence: TRUE# Calibracion global
s1 <- Score(list("FGR" = modelo_fg_s1),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1_sens,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(s1, models = "FGR" )
#Análisis exploratorio sobre el efecto de la isquemia sobre la AM tomando como competitivos muerte y cicatrizacion
#evento_ampu
sob1 <- sob1 %>%
mutate(
evento_ampu = case_when(
estado_unif %in% c("cerro_herida","fallecio") ~ 2,
estado_unif %in% c("amputacion") ~ 1,
estado_unif %in% c("perdido_seguimiento", "persiste") ~ 0,
TRUE ~ NA_real_
)
)
table(sob1$evento_ampu)##
## 0 1 2
## 112 57 237tiempo_evento_am <- Surv(sob1$dias_out,sob1$ampu)
#modelo de cox
#Efecto de la isquemia sobre el tiempo a la amputacion
cox_model_a <- coxph(tiempo_evento_am ~ isq_factor+ edad + ede_factor , data = sob1)
summary(cox_model_a)## Call:
## coxph(formula = tiempo_evento_am ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 406, number of events= 57
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.5974 4.9404 0.3586 4.45 0.0000084 ***
## edad 0.0232 1.0235 0.0138 1.68 0.093 .
## ede_factorUnilateral 0.6092 1.8389 0.2835 2.15 0.032 *
## ede_factorBilateral 0.1648 1.1792 0.4856 0.34 0.734
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 4.94 0.202 2.446 9.98
## edad 1.02 0.977 0.996 1.05
## ede_factorUnilateral 1.84 0.544 1.055 3.21
## ede_factorBilateral 1.18 0.848 0.455 3.05
##
## Concordance= 0.737 (se = 0.029 )
## Likelihood ratio test= 42.6 on 4 df, p=0.00000001
## Wald test = 33.9 on 4 df, p=0.0000008
## Score (logrank) test = 41.6 on 4 df, p=0.00000002tab_model(cox_model_a,
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 | 4.94 | 1.77 | 0.00 – Inf | <0.001 |
| edad | 1.02 | 0.01 | 0.00 – Inf | 0.093 |
| ede factor [Unilateral] | 1.84 | 0.52 | 0.00 – Inf | 0.032 |
| ede factorBilateral | 1.18 | 0.57 | 0.00 – Inf | 0.734 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.125 | |||
confint(cox_model_a)## 2.5 % 97.5 %
## isq_factorSI 0.89454 2.3003
## edad -0.00385 0.0503
## ede_factorUnilateral 0.05352 1.1648
## ede_factorBilateral -0.78693 1.1165# Martingale y Deviance
ggcoxdiagnostics(cox_model_a, type = "martingale", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggcoxdiagnostics(cox_model_a, type = "deviance", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph_a <- cox.zph(cox_model_a)
test_ph_a## chisq df p
## isq_factor 1.568891 1 0.21
## edad 0.000702 1 0.98
## ede_factor 3.858376 2 0.15
## GLOBAL 5.486547 4 0.24#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph_a)## 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.
#modelo de Fine and gray
modelo_fg_a <- FGR(
formula = Hist(dias_out, evento_ampu) ~ isq_factor + edad+ede_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg_a##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 57 0
## 2 237 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_ampu) ~ isq_factor + edad +
## ede_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI 1.7075 5.52 0.3684 4.635 0.0000036
## edad 0.0213 1.02 0.0115 1.850 0.0640000
## ede_factorUnilateral 0.6072 1.84 0.2865 2.120 0.0340000
## ede_factorBilateral 0.1664 1.18 0.4689 0.355 0.7200000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 5.52 0.181 2.679 11.35
## edad 1.02 0.979 0.999 1.04
## ede_factorUnilateral 1.84 0.545 1.047 3.22
## ede_factorBilateral 1.18 0.847 0.471 2.96
##
## Num. cases = 406
## Pseudo Log-likelihood = -310
## Pseudo likelihood ratio test = 47.3 on 4 df,
##
## Convergence: TRUE# Calibracion global
sa <- Score(list("FGR" = modelo_fg_a),
formula = Hist(dias_out, evento_ampu) ~ 1,
data = sob1,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(sa, models = "FGR" )
table(sob1$evento_ampu)##
## 0 1 2
## 112 57 237#cif
#CIF
cifivhxa <- cuminc(Surv(dias_out, factor(evento_ampu)) ~ isq_factor, data=sob1)
cifivhxa##
## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 86 0.049 0.014 0.026, 0.082
## NO/leve 200 22 0.049 0.014 0.026, 0.082
## NO/leve 300 3 0.049 0.014 0.026, 0.082
## SI 100 71 0.244 0.034 0.181, 0.312
## SI 200 24 0.283 0.036 0.214, 0.356
## SI 300 8 0.296 0.038 0.224, 0.372## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 86 0.525 0.034 0.456, 0.590
## NO/leve 200 22 0.769 0.032 0.698, 0.825
## NO/leve 300 3 0.917 0.028 0.842, 0.957
## SI 100 71 0.221 0.034 0.159, 0.290
## SI 200 24 0.448 0.044 0.360, 0.532
## SI 300 8 0.549 0.051 0.445, 0.642## • Tests## outcome statistic df p.value
## 1 40.4 1.00 <0.001
## 2 50.8 1.00 <0.001ggcuminc(
cifivhxa,
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 (amputacion)",
"Evento competitivo (muerte / cicatrizacion)"
)
) +
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()`).
#altman
#altman
edades <- seq(30, 85, by = 1)
new_no_a <- data.frame(
isq_factor = factor(rep(levels(sob1$isq_factor)[1], length(edades)),
levels = levels(sob1$isq_factor)),
edad = edades,
ede_factor = factor(rep("Sin/local", length(edades)),
levels = levels(sob1$ede_factor))
)
new_si_a <- new_no_a
new_si_a$isq_factor <- factor(rep(levels(sob1$isq_factor)[2], length(edades)),
levels = levels(sob1$isq_factor))
pred_no_a <- predictRisk(modelo_fg_a,
newdata = new_no_a,
times = 180,
cause = 1)
pred_si_a <- predictRisk(modelo_fg_a,
newdata = new_si_a,
times = 180,
cause = 1)
plot(edades, pred_no_a,
type = "l",
lwd = 3,
xlab = "Edad",
ylab = "Probabilidad de amputacion a 180 dias",
ylim = c(0, max(pred_no_a, pred_si_a)),
col = "black")
lines(edades, pred_si_a,
lwd = 3,
lty = 2,
col = "black")
legend("bottomleft",
legend = c("Sin isquemia/isq leve", "Con isquemia moderada/grave"),
lwd = 3,
lty = c(1,2),
bty = "n")
#Gray test
cuminc_obja <- cmprsk::cuminc(
ftime = sob1$dias_out,
fstatus = sob1$evento_ampu,
group = sob1$isq_factor
)
cuminc_obja## Tests:
## stat pv df
## 1 40.4 0.00000000020845 1
## 2 50.8 0.00000000000105 1
## Estimates and Variances:
## $est
## 100 200 300
## NO/leve 1 0.0487 0.0487 0.0487
## SI 1 0.2440 0.2833 0.2961
## NO/leve 2 0.5253 0.7688 0.9167
## SI 2 0.2212 0.4480 0.5494
##
## $var
## 100 200 300
## NO/leve 1 0.000207 0.000207 0.000207
## SI 1 0.001141 0.001322 0.001441
## NO/leve 2 0.001165 0.001055 0.000769
## SI 2 0.001148 0.001970 0.002550infección como modificador del efecto de la isquemia sobre la cicatrización
#modelo de Fine and gray
table(sob1$infe_factor)##
## No Leve/Mod Grave
## 86 255 65table(sob1$evento_cica)##
## 0 1 2
## 112 223 71#no se puede hacer fine and gray. voy a hacer un modelo esratificado para la variable infección
sob1_grave <- subset(sob1, infe_factor == "Grave")
sob1_levemod <- subset(sob1, infe_factor == "Leve/Mod")
sob1_sin <- subset(sob1, infe_factor == "No")
#sin infeccion
modelo_fg_sin <- FGR( formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor , data = sob1_sin, cause = 1) # Aclara que Cicatrizacion es el evento de interés
modelo_fg_sin##
## Right-censored response of a competing.risks model
##
## No.Observations: 86
##
## Pattern:
##
## Cause event right.censored
## 1 48 0
## 2 11 0
## unknown 0 27
##
##
## 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_sin, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -1.2354 0.291 0.4967 -2.487 0.013
## edad -0.0178 0.982 0.0155 -1.149 0.250
## ede_factorUnilateral 0.0648 1.067 0.4263 0.152 0.880
## ede_factorBilateral -0.3101 0.733 0.3725 -0.832 0.410
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.291 3.440 0.110 0.77
## edad 0.982 1.018 0.953 1.01
## ede_factorUnilateral 1.067 0.937 0.463 2.46
## ede_factorBilateral 0.733 1.364 0.353 1.52
##
## Num. cases = 86
## Pseudo Log-likelihood = -172
## Pseudo likelihood ratio test = 19.4 on 4 df,
##
## Convergence: TRUE#infeccion leve/mod
modelo_fg_levemod <- FGR( formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor , data = sob1_levemod, cause = 1) # Aclara que Cicatrizacion es el evento de interés
modelo_fg_levemod##
## Right-censored response of a competing.risks model
##
## No.Observations: 255
##
## Pattern:
##
## Cause event right.censored
## 1 147 0
## 2 36 0
## unknown 0 72
##
##
## 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_levemod, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9114 0.402 0.18345 -4.968 0.00000068
## edad -0.0145 0.986 0.00925 -1.570 0.12000000
## ede_factorUnilateral 0.0232 1.023 0.18241 0.127 0.90000000
## ede_factorBilateral -0.0382 0.962 0.28301 -0.135 0.89000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.402 2.488 0.281 0.576
## edad 0.986 1.015 0.968 1.004
## ede_factorUnilateral 1.023 0.977 0.716 1.463
## ede_factorBilateral 0.962 1.039 0.553 1.676
##
## Num. cases = 255
## Pseudo Log-likelihood = -689
## Pseudo likelihood ratio test = 38 on 4 df,
##
## Convergence: TRUE#infeccion grave
modelo_fg_grave <- FGR( formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor , data = sob1_grave, cause = 1) # Aclara que Cicatrizacion es el evento de interés
modelo_fg_grave##
## Right-censored response of a competing.risks model
##
## No.Observations: 65
##
## Pattern:
##
## Cause event right.censored
## 1 28 0
## 2 24 0
## unknown 0 13
##
##
## 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_grave, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -1.9935 0.136 0.6891 -2.8929 0.00380
## edad -0.0115 0.989 0.0276 -0.4161 0.68000
## ede_factorUnilateral -1.4872 0.226 0.4127 -3.6038 0.00031
## ede_factorBilateral 0.0882 1.092 1.2836 0.0687 0.95000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.136 7.341 0.0353 0.526
## edad 0.989 1.012 0.9364 1.044
## ede_factorUnilateral 0.226 4.425 0.1007 0.507
## ede_factorBilateral 1.092 0.916 0.0882 13.517
##
## Num. cases = 65
## Pseudo Log-likelihood = -87.9
## Pseudo likelihood ratio test = 27.4 on 4 df,
##
## Convergence: TRUE#gráficos de cif para cada categoría de infección
#SIN INFECCIÓN
#cif
#CIF
cifivhxasin <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor, data=sob1_sin)
cifivhxasin## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 37 0.295 0.063 0.180, 0.420
## NO/leve 100 25 0.464 0.069 0.326, 0.592
## NO/leve 150 12 0.645 0.071 0.488, 0.765
## NO/leve 200 7 0.709 0.072 0.541, 0.824
## NO/leve 250 4 0.793 0.074 0.599, 0.901
## NO/leve 300 1 0.963 0.089 0.011, 1.00
## SI 50.0 19 0.097 0.054 0.024, 0.232
## SI 100 13 0.210 0.079 0.082, 0.377
## SI 150 10 0.210 0.079 0.082, 0.377
## SI 200 6 0.269 0.095 0.109, 0.460
## SI 250 3 0.269 0.095 0.109, 0.460
## SI 300 2 0.269 0.095 0.109, 0.460## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 37 0.019 0.019 0.001, 0.088
## NO/leve 100 25 0.037 0.026 0.007, 0.114
## NO/leve 150 12 0.037 0.026 0.007, 0.114
## NO/leve 200 7 0.037 0.026 0.007, 0.114
## NO/leve 250 4 0.037 0.026 0.007, 0.114
## NO/leve 300 1 0.037 0.026 0.007, 0.114
## SI 50.0 19 0.233 0.079 0.100, 0.398
## SI 100 13 0.268 0.083 0.124, 0.437
## SI 150 10 0.312 0.090 0.151, 0.488
## SI 200 6 0.312 0.090 0.151, 0.488
## SI 250 3 0.312 0.090 0.151, 0.488
## SI 300 2 0.312 0.090 0.151, 0.488## • Tests## outcome statistic df p.value
## 1 15.5 1.00 <0.001
## 2 12.1 1.00 <0.001ggcuminc(
cifivhxasin,
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)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_step()`).
#INFECCIon leve/mod
#cif
#CIF
cifivhxasleve <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor, data=sob1_levemod)
cifivhxasleve##
## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 44 0.574 0.044 0.482, 0.655
## NO/leve 200 8 0.824 0.041 0.727, 0.890
## NO/leve 300 1 0.902 0.054 0.723, 0.968
## SI 100 50 0.210 0.041 0.136, 0.296
## SI 200 13 0.486 0.056 0.373, 0.590
## SI 300 5 0.579 0.060 0.452, 0.687## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 44 0.051 0.019 0.022, 0.096
## NO/leve 200 8 0.051 0.019 0.022, 0.096
## NO/leve 300 1 0.051 0.019 0.022, 0.096
## SI 100 50 0.231 0.042 0.154, 0.316
## SI 200 13 0.288 0.046 0.202, 0.380
## SI 300 5 0.288 0.046 0.202, 0.380## • Tests## outcome statistic df p.value
## 1 37.1 1.00 <0.001
## 2 22.2 1.00 <0.001ggcuminc(
cifivhxasleve,
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 (amputacion)",
"Evento competitivo (muerte / cicatrizacion)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_step()`).
#INFECCIon grave
#cif
#CIF
cifivhxasgrave <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor, data=sob1_grave)
cifivhxasgrave##
## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 25 0.111 0.053 0.034, 0.239
## NO/leve 100 17 0.283 0.078 0.145, 0.438
## NO/leve 150 12 0.428 0.086 0.260, 0.586
## NO/leve 200 7 0.527 0.089 0.341, 0.683
## NO/leve 250 1 0.734 0.085 0.524, 0.862
## NO/leve 300 1 0.734 0.085 0.524, 0.862
## SI 50.0 16 0.000 0.000 NA, NA
## SI 100 8 0.047 0.048 0.003, 0.203
## SI 150 5 0.164 0.090 0.038, 0.369
## SI 200 5 0.164 0.090 0.038, 0.369
## SI 250 3 0.164 0.090 0.038, 0.369
## SI 300 1 0.320 0.126 0.106, 0.560## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 25 0.194 0.067 0.084, 0.338
## NO/leve 100 17 0.223 0.071 0.103, 0.372
## NO/leve 150 12 0.223 0.071 0.103, 0.372
## NO/leve 200 7 0.223 0.071 0.103, 0.372
## NO/leve 250 1 0.223 0.071 0.103, 0.372
## NO/leve 300 1 0.223 0.071 0.103, 0.372
## SI 50.0 16 0.448 0.095 0.261, 0.619
## SI 100 8 0.485 0.096 0.291, 0.654
## SI 150 5 0.544 0.104 0.324, 0.719
## SI 200 5 0.544 0.104 0.324, 0.719
## SI 250 3 0.602 0.109 0.360, 0.777
## SI 300 1 0.602 0.109 0.360, 0.777## • Tests## outcome statistic df p.value
## 1 9.28 1.00 0.002
## 2 8.13 1.00 0.004ggcuminc(
cifivhxasgrave,
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)"
)
) +
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()`).
#Análisis secundario
table(sob1$isq_factor2)##
## NO LEVE MODERADA GRAVE
## 206 31 47 122#Efecto de la isquemia sobre el tiempo a la cicatrización sin eventos competitivos ajustado por edad y edema
cox_modelisq2 <- coxph(tiempo_evento ~ isq_factor2+ edad + ede_factor , data = sob1)
summary(cox_modelisq2)## Call:
## coxph(formula = tiempo_evento ~ isq_factor2 + edad + ede_factor,
## data = sob1)
##
## n= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factor2LEVE 0.0302 1.0306 0.2251 0.13 0.89339
## isq_factor2MODERADA -0.7480 0.4733 0.2394 -3.12 0.00178 **
## isq_factor2GRAVE -0.6937 0.4997 0.2008 -3.45 0.00055 ***
## edad -0.0163 0.9838 0.0077 -2.12 0.03421 *
## ede_factorUnilateral -0.2029 0.8163 0.1547 -1.31 0.18953
## ede_factorBilateral -0.1465 0.8638 0.2356 -0.62 0.53413
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factor2LEVE 1.031 0.97 0.663 1.602
## isq_factor2MODERADA 0.473 2.11 0.296 0.757
## isq_factor2GRAVE 0.500 2.00 0.337 0.741
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.816 1.22 0.603 1.105
## ede_factorBilateral 0.864 1.16 0.544 1.371
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.6 on 6 df, p=0.0000003
## Wald test = 37.2 on 6 df, p=0.000002
## Score (logrank) test = 39 on 6 df, p=0.0000007tab_model(cox_modelisq2,
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 factor2LEVE | 1.03 | 0.23 | 0.00 – Inf | 0.893 |
| isq factor2MODERADA | 0.47 | 0.11 | 0.00 – Inf | 0.002 |
| isq factor2GRAVE | 0.50 | 0.10 | 0.00 – Inf | 0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.034 |
| ede factor [Unilateral] | 0.82 | 0.13 | 0.00 – Inf | 0.190 |
| ede factorBilateral | 0.86 | 0.20 | 0.00 – Inf | 0.534 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.096 | |||
confint(cox_modelisq2)## 2.5 % 97.5 %
## isq_factor2LEVE -0.4110 0.47135
## isq_factor2MODERADA -1.2172 -0.27882
## isq_factor2GRAVE -1.0872 -0.30013
## edad -0.0314 -0.00121
## ede_factorUnilateral -0.5061 0.10023
## ede_factorBilateral -0.6082 0.31525#
cifivhxas22 <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor2, data=sob1)
cifivhxas22## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO 100 73 0.503 0.037 0.429, 0.572
## NO 200 20 0.726 0.036 0.647, 0.790
## NO 300 3 0.875 0.035 0.787, 0.928
## LEVE 100 13 0.478 0.096 0.284, 0.649
## LEVE 200 2 0.858 0.086 0.578, 0.958
## LEVE 300 0 0.968 0.064 0.190, 0.999
## MODERADA 100 22 0.240 0.068 0.122, 0.380
## MODERADA 200 12 0.463 0.083 0.297, 0.613
## MODERADA 300 3 0.636 0.106 0.394, 0.802
## GRAVE 100 49 0.162 0.035 0.100, 0.237
## GRAVE 200 12 0.376 0.053 0.273, 0.479
## GRAVE 300 5 0.442 0.058 0.327, 0.551## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO 100 73 0.082 0.020 0.049, 0.126
## NO 200 20 0.082 0.020 0.049, 0.126
## NO 300 3 0.082 0.020 0.049, 0.126
## LEVE 100 13 0.032 0.032 0.002, 0.144
## LEVE 200 2 0.032 0.032 0.002, 0.144
## LEVE 300 0 0.032 0.032 0.002, 0.144
## MODERADA 100 22 0.161 0.057 0.069, 0.287
## MODERADA 200 12 0.161 0.057 0.069, 0.287
## MODERADA 300 3 0.161 0.057 0.069, 0.287
## GRAVE 100 49 0.329 0.044 0.246, 0.415
## GRAVE 200 12 0.407 0.048 0.313, 0.498
## GRAVE 300 5 0.426 0.050 0.327, 0.522## • Tests## outcome statistic df p.value
## 1 57.7 3.00 <0.001
## 2 53.7 3.00 <0.001ggcuminc(
cifivhxas22,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO" = "#56B4E9", "LEVE" = "#E69F00", "MODERADA"="#B58", "GRAVE"="#009E73")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c(
"Evento de interes (cicatrizacion)",
"Evento competitivo (muerte / amputacion)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 34 rows containing missing values or values outside the scale range
## (`geom_step()`).
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph_isq2 <- cox.zph(cox_modelisq2)
test_ph_isq2## chisq df p
## isq_factor2 3.8524 3 0.28
## edad 0.0327 1 0.86
## ede_factor 1.3943 2 0.50
## GLOBAL 5.6890 6 0.46#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph_isq2)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.## 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.
#modelo de Fine and gray
modelo_fg_isq2 <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor2 + edad+ede_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg_isq2##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor2 + edad +
## ede_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factor2LEVE 0.0950 1.100 0.19799 0.480 0.630000000
## isq_factor2MODERADA -0.6968 0.498 0.22574 -3.087 0.002000000
## isq_factor2GRAVE -1.0946 0.335 0.20344 -5.380 0.000000074
## edad -0.0163 0.984 0.00782 -2.082 0.037000000
## ede_factorUnilateral -0.3405 0.711 0.15203 -2.240 0.025000000
## ede_factorBilateral -0.2122 0.809 0.23300 -0.911 0.360000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factor2LEVE 1.100 0.909 0.746 1.621
## isq_factor2MODERADA 0.498 2.007 0.320 0.775
## isq_factor2GRAVE 0.335 2.988 0.225 0.499
## edad 0.984 1.016 0.969 0.999
## ede_factorUnilateral 0.711 1.406 0.528 0.958
## ede_factorBilateral 0.809 1.236 0.512 1.277
##
## Num. cases = 406
## Pseudo Log-likelihood = -1150
## Pseudo likelihood ratio test = 73.1 on 6 df,
##
## Convergence: TRUE# Calibracion global
s_isq2 <- Score(list("FGR" = modelo_fg_isq2),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(s_isq2, models = "FGR" )
#Análisis exploratorio con variables desbalanceadas
#modelo de Fine and gray con localizacion: modifica el coeficiente en 10.4% a 0.34
modelo_fgL <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor+ loca_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fgL##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor + loca_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -1.074 0.342 0.17099 -6.282 0.00000000034
## edad -0.016 0.984 0.00786 -2.039 0.04100000000
## ede_factorUnilateral -0.348 0.706 0.14946 -2.326 0.02000000000
## ede_factorBilateral -0.211 0.810 0.23538 -0.895 0.37000000000
## loca_factorMet/Tar -0.423 0.655 0.13738 -3.081 0.00210000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.342 2.93 0.244 0.478
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.706 1.42 0.527 0.947
## ede_factorBilateral 0.810 1.23 0.511 1.285
## loca_factorMet/Tar 0.655 1.53 0.500 0.857
##
## Num. cases = 406
## Pseudo Log-likelihood = -1146
## Pseudo likelihood ratio test = 80.3 on 5 df,
##
## Convergence: TRUE#modelo de Fine and gray con topo: modifica el coeficiente en a 0.41 en 7.8%
modelo_fgt <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor+ topo_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fgt##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor + topo_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.8826 0.414 0.16409 -5.379 0.000000075
## edad -0.0171 0.983 0.00736 -2.318 0.020000000
## ede_factorUnilateral -0.2019 0.817 0.15061 -1.340 0.180000000
## ede_factorBilateral -0.2156 0.806 0.24127 -0.894 0.370000000
## topo_factorlat/Mas de dos -0.5182 0.596 0.13943 -3.716 0.000200000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.414 2.42 0.300 0.571
## edad 0.983 1.02 0.969 0.997
## ede_factorUnilateral 0.817 1.22 0.608 1.098
## ede_factorBilateral 0.806 1.24 0.502 1.293
## topo_factorlat/Mas de dos 0.596 1.68 0.453 0.783
##
## Num. cases = 406
## Pseudo Log-likelihood = -1144
## Pseudo likelihood ratio test = 84 on 5 df,
##
## Convergence: TRUE# modelo de fine and gray con número de zonas: modifica el coeficiente en a 0.41 en 7.8%
modelo_fgNum <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor+ zon_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fgNum##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor + zon_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.8937 0.409 0.16416 -5.444 0.000000052
## edad -0.0173 0.983 0.00771 -2.239 0.025000000
## ede_factorUnilateral -0.0793 0.924 0.16078 -0.493 0.620000000
## ede_factorBilateral -0.2231 0.800 0.24875 -0.897 0.370000000
## zon_factorDos o mas -0.7827 0.457 0.15681 -4.992 0.000000600
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.409 2.44 0.297 0.564
## edad 0.983 1.02 0.968 0.998
## ede_factorUnilateral 0.924 1.08 0.674 1.266
## ede_factorBilateral 0.800 1.25 0.491 1.303
## zon_factorDos o mas 0.457 2.19 0.336 0.622
##
## Num. cases = 406
## Pseudo Log-likelihood = -1138
## Pseudo likelihood ratio test = 97 on 5 df,
##
## Convergence: TRUE# modelo de fine and gray con neuropatía: modifica el coeficiente a 0.40a 5.2%
modelo_fgNeu <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor+ neur_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fgNeu##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor + neur_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9107 0.402 0.17366 -5.24 0.00000016
## edad -0.0183 0.982 0.00772 -2.37 0.01800000
## ede_factorUnilateral -0.3669 0.693 0.14949 -2.45 0.01400000
## ede_factorBilateral -0.2190 0.803 0.23313 -0.94 0.35000000
## neur_factorMod/charcot 0.4820 1.619 0.38609 1.25 0.21000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.402 2.486 0.286 0.565
## edad 0.982 1.018 0.967 0.997
## ede_factorUnilateral 0.693 1.443 0.517 0.929
## ede_factorBilateral 0.803 1.245 0.509 1.269
## neur_factorMod/charcot 1.619 0.618 0.760 3.451
##
## Num. cases = 406
## Pseudo Log-likelihood = -1150
## Pseudo likelihood ratio test = 73 on 5 df,
##
## Convergence: TRUE# modelo de fine and gray con profundidad: no modifica el coeficiente
modelo_fgprof <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor+ prof_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fgprof##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor + prof_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9516 0.386 0.16724 -5.690 0.000000013
## edad -0.0164 0.984 0.00778 -2.106 0.035000000
## ede_factorUnilateral -0.3156 0.729 0.15075 -2.094 0.036000000
## ede_factorBilateral -0.1855 0.831 0.22470 -0.826 0.410000000
## prof_factorparcial/total -0.4398 0.644 0.29320 -1.500 0.130000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.386 2.59 0.278 0.536
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.729 1.37 0.543 0.980
## ede_factorBilateral 0.831 1.20 0.535 1.290
## prof_factorparcial/total 0.644 1.55 0.363 1.144
##
## Num. cases = 406
## Pseudo Log-likelihood = -1149
## Pseudo likelihood ratio test = 74.1 on 5 df,
##
## Convergence: TRUE# modelo de fine and gray con area: no modifica el coeficiente
modelo_fgA <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor+ area_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fgA##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor + area_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.94976 0.387 0.16791 -5.6563 0.00000001500
## edad -0.01758 0.983 0.00785 -2.2409 0.02500000000
## ede_factorUnilateral 0.00781 1.008 0.15143 0.0516 0.96000000000
## ede_factorBilateral -0.26196 0.770 0.24759 -1.0581 0.29000000000
## area_factor10 a 40 -0.33746 0.714 0.14170 -2.3816 0.01700000000
## area_factorMas de 40 -1.59324 0.203 0.25923 -6.1460 0.00000000079
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.387 2.585 0.278 0.538
## edad 0.983 1.018 0.968 0.998
## ede_factorUnilateral 1.008 0.992 0.749 1.356
## ede_factorBilateral 0.770 1.299 0.474 1.250
## area_factor10 a 40 0.714 1.401 0.541 0.942
## area_factorMas de 40 0.203 4.920 0.122 0.338
##
## Num. cases = 406
## Pseudo Log-likelihood = -1130
## Pseudo likelihood ratio test = 113 on 6 df,
##
## Convergence: TRUE# modelo de fine and gray con fase: no modifica el coeficiente a 0.40 (5.2%)
modelo_fgfase <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor+ fase_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fgfase##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor + fase_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9238 0.397 0.16916 -5.461 0.000000047
## edad -0.0186 0.982 0.00766 -2.436 0.015000000
## ede_factorUnilateral -0.3276 0.721 0.15227 -2.151 0.031000000
## ede_factorBilateral -0.1708 0.843 0.23225 -0.735 0.460000000
## fase_factorGranu/infla -0.7773 0.460 0.38419 -2.023 0.043000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.397 2.52 0.285 0.553
## edad 0.982 1.02 0.967 0.996
## ede_factorUnilateral 0.721 1.39 0.535 0.971
## ede_factorBilateral 0.843 1.19 0.535 1.329
## fase_factorGranu/infla 0.460 2.18 0.216 0.976
##
## Num. cases = 406
## Pseudo Log-likelihood = -1148
## Pseudo likelihood ratio test = 77.1 on 5 df,
##
## Convergence: TRUE