Tesis Maestría en efectividad clínica
Carro Gabriela
2026-02-18
#CARGA DE LA BASE DE DATOS
#carga de datos
library(readxl)
sob1 <- read_excel("C:/Users/Administrador/Downloads/TESIS_may.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
## 170 236prop.table(table1)##
## 0 1
## 0.419 0.581#PROPORCION AMPUTADOS
table2=table(sob1$ampu)
table2##
## 0 1
## 348 58prop.table(table2)##
## 0 1
## 0.857 0.143#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# Número de hombres
n_hombres <- sum(sob1$sexo == "M", na.rm = TRUE)
# Total de pacientes
n_total <- sum(!is.na(sob1$sexo))
# Proporción e IC95%
prop.test(n_hombres, n_total)##
## 1-sample proportions test with continuity correction
##
## data: n_hombres out of n_total, null probability 0.5
## X-squared = 143, df = 1, p-value <0.0000000000000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.755 0.835
## sample estimates:
## p
## 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 100 80.6 78.5 89.3 66 1 382 381 1.15 0.79 4#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
## 170 236#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
## 348 58#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
## 58 236 14 98sob1 <- 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
## 58 236 14 98table(sob1$estado)##
## ampu cerro fallecio persiste
## 58 236 14 98table(sob1$cicat_factor)##
## No Cicatrizo Cicatrizo
## 170 236##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 58 406 0.143 0.0174 0.109 0.177 14.3 10.9-17.7%
## 2 cerro_herida 236 406 0.581 0.0245 0.533 0.629 58.1 53.3-62.9%
## 3 fallecio 14 406 0.0345 0.00906 0.0167 0.0522 3.4 1.7-5.2%
## 4 persiste 98 406 0.241 0.0212 0.200 0.283 24.1 20-28.3%tabla_estado## # A tibble: 4 × 7
## estado_unif n N_total prop se li ls
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 amputacion 58 406 0.143 0.0174 0.109 0.177
## 2 cerro_herida 236 406 0.581 0.0245 0.533 0.629
## 3 fallecio 14 406 0.0345 0.00906 0.0167 0.0522
## 4 persiste 98 406 0.241 0.0212 0.200 0.283#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 100 80.6 78.5 89.3 66 1 382 381 1.15 0.79 4describeBy(
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 58 41.0 41.0 22.5 34.4 17.8 1 211 210
## X12 2 cerro_herida 1 236 105.1 74.0 86.5 94.8 62.3 9 330 321
## X13 3 fallecio 1 14 43.1 40.4 24.0 39.9 25.2 5 120 115
## X14 4 persiste 1 98 132.1 95.1 114.5 123.3 103.0 6 382 376
## skew kurtosis se
## X11 1.834 3.677 5.38
## X12 1.128 0.624 4.82
## X13 0.661 -1.253 10.80
## X14 0.693 -0.286 9.60tapply(sob1$dias_out,
sob1$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 14.0 56.2
##
## $cerro_herida
## 25% 75%
## 47 138
##
## $fallecio
## 25% 75%
## 11.2 69.8
##
## $persiste
## 25% 75%
## 60 199#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 94.2 72.7 77 84.8 62.3 5 349 344 1.15 0.91 4.72describeBy(
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 12 34.4 36.9 20.5 28.5 13.3 7 121 114
## X12 2 cerro_herida 1 168 93.4 65.2 77.0 84.1 53.4 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 51 119.3 91.3 115.0 110.1 108.2 6 349 343
## skew kurtosis se
## X11 1.379 0.327 10.66
## X12 1.225 0.994 5.03
## X13 0.708 -1.298 9.44
## X14 0.630 -0.343 12.79tapply(sob1_per_noisq$dias_out,
sob1_per_noisq$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 13.2 33.2
##
## $cerro_herida
## 25% 75%
## 42 121
##
## $fallecio
## 25% 75%
## 10.8 34.8
##
## $persiste
## 25% 75%
## 45.5 185.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 109 90.1 86 97 86 1 382 381 1.02 0.27 6.93describeBy(
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 68 134.1 85.9 113.5 127.5 76.4 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 47 146.0 98.1 114.0 137.7 84.5 7 382 375
## skew kurtosis se
## X11 1.847 3.806 6.21
## X12 0.722 -0.468 10.42
## X13 0.134 -1.992 16.58
## X14 0.710 -0.496 14.31tapply(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%
## 67.2 186.2
##
## $fallecio
## 25% 75%
## 11.8 99.0
##
## $persiste
## 25% 75%
## 70 216#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 = 1007, p-value = 0.2
## alternative hypothesis: true location shift is not equal to 0#Tabla incidencia + intervalos de confianza
#ISQ
tabla_estado_i <- sob1_per_isq %>%
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_i## # A tibble: 4 × 7
## estado_unif n N_total prop se li ls
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 amputacion 46 169 0.272 0.0342 0.205 0.339
## 2 cerro_herida 68 169 0.402 0.0377 0.328 0.476
## 3 fallecio 8 169 0.0473 0.0163 0.0153 0.0794
## 4 persiste 47 169 0.278 0.0345 0.211 0.346#Tabla incidencia + intervalos de confianza
#NO ISQ
tabla_estado_ni <- sob1_per_noisq %>%
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_ni## # A tibble: 4 × 7
## estado_unif n N_total prop se li ls
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 amputacion 12 237 0.0506 0.0142 0.0227 0.0785
## 2 cerro_herida 168 237 0.709 0.0295 0.651 0.767
## 3 fallecio 6 237 0.0253 0.0102 0.00532 0.0453
## 4 persiste 51 237 0.215 0.0267 0.163 0.268#diferencia entre cicatrizados, amputados y fallecios en isq y no isq
library(DescTools)##
## Adjuntando el paquete: 'DescTools'## The following objects are masked from 'package:psych':
##
## AUC, ICC, SD## The following objects are masked from 'package:caret':
##
## MAE, RMSE## The following objects are masked from 'package:Hmisc':
##
## %nin%, Label, Mean, Quantile## The following object is masked from 'package:data.table':
##
## %like%library(emmeans)## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'##
## Adjuntando el paquete: 'emmeans'## The following object is masked from 'package:rms':
##
## contrastlibrary(expss)## Cargando paquete requerido: maditr##
## To modify variables or add new variables:
## let(mtcars, new_var = 42, new_var2 = new_var*hp) %>% head()##
## Adjuntando el paquete: 'maditr'## The following object is masked from 'package:DescTools':
##
## %like%## The following objects are masked from 'package:dplyr':
##
## between, coalesce, first, last## The following object is masked from 'package:purrr':
##
## transpose## The following objects are masked from 'package:data.table':
##
## copy, dcast, let, melt## The following object is masked from 'package:readr':
##
## cols## Registered S3 methods overwritten by 'expss':
## method from
## [.labelled Hmisc
## as.data.frame.labelled base
## print.labelled Hmisc##
## Use 'expss_output_rnotebook()' to display tables inside R Notebooks.
## To return to the console output, use 'expss_output_default()'.##
## Adjuntando el paquete: 'expss'## The following object is masked from 'package:naniar':
##
## is_na## The following objects are masked from 'package:broom.helpers':
##
## contains, vars, where## The following object is masked from 'package:ggpubr':
##
## compare_means## The following objects are masked from 'package:gtsummary':
##
## contains, vars, where## The following objects are masked from 'package:stringr':
##
## fixed, regex## The following objects are masked from 'package:dplyr':
##
## compute, contains, na_if, recode, vars, where## The following objects are masked from 'package:purrr':
##
## keep, modify, modify_if, when## The following objects are masked from 'package:tidyr':
##
## contains, nest## The following object is masked from 'package:ggplot2':
##
## vars## The following objects are masked from 'package:data.table':
##
## copy, likelibrary(descr)
chisq.test(sob1$ampu_factor,sob1$isq_factor)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: sob1$ampu_factor and sob1$isq_factor
## X-squared = 38, df = 1, p-value = 0.0000000008chisq.test(sob1$cicat_factor,sob1$isq_factor)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: sob1$cicat_factor and sob1$isq_factor
## X-squared = 37, df = 1, p-value = 0.000000001chisq.test(sob1$fallecio_factor,sob1$isq_factor)##
## Pearson's Chi-squared test with Yates' continuity correction
##
## data: sob1$fallecio_factor and sob1$isq_factor
## X-squared = 0.9, df = 1, p-value = 0.4#Densidad de incidencia y razón de densidades
#cicatrizacion totaldensidad de incidencia evento/1000 personas/dia
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.0058tiempo1_meses <- tiempo1_total / 30.44
tasa1_mensual <- eventos / tiempo1_meses
tasa1_mensual## [1] 0.176tasa1_1000 <- tasa1_dia * 1000
tasa1_1000## [1] 5.8#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.0037tiempo2_meses <- tiempo2_total / 30.44
tasa2_mensual <- eventos / tiempo2_meses
tasa2_mensual## [1] 0.113tasa2_1000 <- tasa2_dia * 1000
tasa2_1000## [1] 3.7#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.00752tiempo3_meses <- tiempo3_total / 30.44
tasa3_mensual <- eventos / tiempo3_meses
tasa3_mensual## [1] 0.229tasa3_1000 <- tasa3_dia * 1000
tasa3_1000## [1] 7.52# 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.00142tiempo1_meses <- tiempo1_total / 30.44
tasa1_mensual <- eventos / tiempo1_meses
tasa1_mensual## [1] 0.0434tasa1_1000 <- tasa1_dia * 1000
tasa1_1000## [1] 1.42#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.0025tiempo2_meses <- tiempo2_total / 30.44
tasa2_mensual <- eventos / tiempo2_meses
tasa2_mensual## [1] 0.0761tasa2_1000 <- tasa2_dia * 1000
tasa2_1000## [1] 2.5#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.000537tiempo3_meses <- tiempo3_total / 30.44
tasa3_mensual <- eventos / tiempo3_meses
tasa3_mensual## [1] 0.0164tasa3_1000 <- tasa3_dia * 1000
tasa3_1000## [1] 0.537#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.000344tiempo1_meses <- tiempo1_total / 30.44
tasa1_mensual <- eventos / tiempo1_meses
tasa1_mensual## [1] 0.0105tasa1_1000 <- tasa1_dia * 1000
tasa1_1000## [1] 0.344#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.000435tiempo2_meses <- tiempo2_total / 30.44
tasa2_mensual <- eventos / tiempo2_meses
tasa2_mensual## [1] 0.0132tasa2_1000 <- tasa2_dia * 1000
tasa2_1000## [1] 0.435#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.000269tiempo3_meses <- tiempo3_total / 30.44
tasa3_mensual <- eventos / tiempo3_meses
tasa3_mensual## [1] 0.00818tasa3_1000 <- tasa3_dia * 1000
tasa3_1000## [1] 0.269#COMPORTAMIENTO DE LAS VARIABLES
#comportamiento de las variables
#edad % e IC
#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
## 69 168xci <- 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 = 41, df = 1, p-value = 0.0000000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.646 0.765
## sample estimates:
## p
## 0.709#porcentaje de amputacion en no isq/leve
table(sob1_per_noisq$ampu_factor)##
## No Amputado Amputado
## 225 12xa <- 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 = 190, df = 1, p-value <0.0000000000000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.0276 0.0890
## sample estimates:
## p
## 0.0506#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
## 101 68xcii <- 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 = 6, df = 1, p-value = 0.01
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.329 0.481
## sample estimates:
## p
## 0.402#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 168 126 13.9 30.7
## isq_factor=SI 169 68 110 16.0 30.7
##
## Chisq= 30.7 on 1 degrees of freedom, p= 0.00000003summary(km_fit)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## isq_factor=NO/leve 237 237 237 168 122 6.9 91 84
## isq_factor=SI 169 169 169 68 196 11.8 170 153
## 0.95UCL
## isq_factor=NO/leve 109
## isq_factor=SI 239#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 239#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= 236
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.6761 0.5086 0.1607 -4.21 0.000026 ***
## edad -0.0156 0.9845 0.0075 -2.08 0.038 *
## ede_factorUnilateral -0.2211 0.8016 0.1482 -1.49 0.136
## ede_factorBilateral -0.2344 0.7911 0.2338 -1.00 0.316
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.509 1.97 0.371 0.697
## edad 0.985 1.02 0.970 0.999
## ede_factorUnilateral 0.802 1.25 0.600 1.072
## ede_factorBilateral 0.791 1.26 0.500 1.251
##
## Concordance= 0.627 (se = 0.021 )
## Likelihood ratio test= 39.3 on 4 df, p=0.00000006
## Wald test = 36.5 on 4 df, p=0.0000002
## Score (logrank) test = 38 on 4 df, p=0.0000001tab_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.51 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.038 |
| ede factor [Unilateral] | 0.80 | 0.12 | 0.00 – Inf | 0.136 |
| ede factorBilateral | 0.79 | 0.18 | 0.00 – Inf | 0.316 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.093 | |||
confint(cox_model)## 2.5 % 97.5 %
## isq_factorSI -0.9910 -0.361092
## edad -0.0303 -0.000899
## ede_factorUnilateral -0.5115 0.069301
## ede_factorBilateral -0.6926 0.223788##Supuesto de proporcionalidad
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph <- cox.zph(cox_model)
test_ph## chisq df p
## isq_factor 3.263 1 0.071
## edad 0.149 1 0.700
## ede_factor 1.823 2 0.402
## GLOBAL 5.381 4 0.250#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.6269 0.0209#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.2538 0.2616 0.2417 0.0199 0.2339 200
## R2 0.0926 0.1032 0.0858 0.0173 0.0752 200
## Slope 1.0000 1.0000 0.9258 0.0742 0.9258 200
## D 0.0163 0.0185 0.0150 0.0035 0.0128 200
## U -0.0008 -0.0009 0.0006 -0.0014 0.0006 200
## Q 0.0171 0.0193 0.0144 0.0049 0.0123 200
## g 0.5060 0.5324 0.4803 0.0521 0.4539 200round(cox_rms[,5],2)## Dxy R2 Slope D U Q g
## 0.23 0.08 0.93 0.01 0.00 0.01 0.45##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
## 98 236 72#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 = 51, df = 2, p-value = 0.00000000001chisq.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 = 38, df = 1, p-value = 0.0000000008chisq.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= 236
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.6761 0.5086 0.1607 -4.21 0.000026 ***
## edad -0.0156 0.9845 0.0075 -2.08 0.038 *
## ede_factorUnilateral -0.2211 0.8016 0.1482 -1.49 0.136
## ede_factorBilateral -0.2344 0.7911 0.2338 -1.00 0.316
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.509 1.97 0.371 0.697
## edad 0.985 1.02 0.970 0.999
## ede_factorUnilateral 0.802 1.25 0.600 1.072
## ede_factorBilateral 0.791 1.26 0.500 1.251
##
## Concordance= 0.627 (se = 0.021 )
## Likelihood ratio test= 39.3 on 4 df, p=0.00000006
## Wald test = 36.5 on 4 df, p=0.0000002
## Score (logrank) test = 38 on 4 df, p=0.0000001cox_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= 72
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.2297 3.4202 0.2950 4.17 0.000031 ***
## edad 0.0264 1.0267 0.0124 2.13 0.033 *
## ede_factorUnilateral 0.6125 1.8450 0.2531 2.42 0.016 *
## ede_factorBilateral 0.2554 1.2909 0.4148 0.62 0.538
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 3.42 0.292 1.918 6.10
## edad 1.03 0.974 1.002 1.05
## ede_factorUnilateral 1.85 0.542 1.124 3.03
## ede_factorBilateral 1.29 0.775 0.573 2.91
##
## Concordance= 0.709 (se = 0.03 )
## Likelihood ratio test= 41.4 on 4 df, p=0.00000002
## Wald test = 36.2 on 4 df, p=0.0000003
## Score (logrank) test = 41 on 4 df, p=0.00000003##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 4 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 236 0
## 2 72 0
## unknown 0 98
##
##
## ----------> Cause: 1
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 406, number of events= 236
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.6761 0.5086 0.1607 -4.21 0.000026 ***
## edad -0.0156 0.9845 0.0075 -2.08 0.038 *
## ede_factorUnilateral -0.2211 0.8016 0.1482 -1.49 0.136
## ede_factorBilateral -0.2344 0.7911 0.2338 -1.00 0.316
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.509 1.97 0.371 0.697
## edad 0.985 1.02 0.970 0.999
## ede_factorUnilateral 0.802 1.25 0.600 1.072
## ede_factorBilateral 0.791 1.26 0.500 1.251
##
## Concordance= 0.627 (se = 0.021 )
## Likelihood ratio test= 39.3 on 4 df, p=0.00000006
## Wald test = 36.5 on 4 df, p=0.0000002
## Score (logrank) test = 38 on 4 df, p=0.0000001
##
##
##
## ----------> Cause: 2
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 406, number of events= 72
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.2297 3.4202 0.2950 4.17 0.000031 ***
## edad 0.0264 1.0267 0.0124 2.13 0.033 *
## ede_factorUnilateral 0.6125 1.8450 0.2531 2.42 0.016 *
## ede_factorBilateral 0.2554 1.2909 0.4148 0.62 0.538
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 3.42 0.292 1.918 6.10
## edad 1.03 0.974 1.002 1.05
## ede_factorUnilateral 1.85 0.542 1.124 3.03
## ede_factorBilateral 1.29 0.775 0.573 2.91
##
## Concordance= 0.709 (se = 0.03 )
## Likelihood ratio test= 41.4 on 4 df, p=0.00000002
## Wald test = 36.2 on 4 df, p=0.0000003
## Score (logrank) test = 41 on 4 df, p=0.00000003#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 236 0
## 2 72 0
## unknown 0 98
##
##
## 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.9254 0.396 0.15978 -5.79 0.000000007
## edad -0.0169 0.983 0.00754 -2.24 0.025000000
## ede_factorUnilateral -0.3528 0.703 0.14683 -2.40 0.016000000
## ede_factorBilateral -0.2780 0.757 0.23411 -1.19 0.240000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.396 2.52 0.290 0.542
## edad 0.983 1.02 0.969 0.998
## ede_factorUnilateral 0.703 1.42 0.527 0.937
## ede_factorBilateral 0.757 1.32 0.479 1.198
##
## Num. cases = 406
## Pseudo Log-likelihood = -1223
## Pseudo likelihood ratio test = 70.1 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 55.3 0.000000000000106 1
## 2 38.6 0.000000000509064 1
## Estimates and Variances:
## $est
## 100 200 300
## NO/leve 1 0.501 0.7403 0.8732
## SI 1 0.180 0.3981 0.5220
## NO/leve 2 0.075 0.0804 0.0804
## SI 2 0.281 0.3336 0.3432
##
## $var
## 100 200 300
## NO/leve 1 0.001158 0.001028 0.000839
## SI 1 0.000964 0.001788 0.002283
## NO/leve 2 0.000308 0.000334 0.000334
## SI 2 0.001252 0.001449 0.001501##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.496 0.1210 0.259 0.734
## 2: FG model 18 0.518 0.1078 0.306 0.729
## 3: FG model 19 0.518 0.1078 0.306 0.729
## 4: FG model 20 0.518 0.1078 0.306 0.729
## 5: FG model 21 0.554 0.0896 0.378 0.729
## ---
## 180: FG model 196 0.725 0.0295 0.668 0.783
## 181: FG model 197 0.723 0.0298 0.664 0.781
## 182: FG model 198 0.723 0.0298 0.664 0.781
## 183: FG model 199 0.727 0.0298 0.669 0.785
## 184: FG model 200 0.726 0.0300 0.667 0.785##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= 236
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.6755 0.5089 0.1607 -4.20 0.000026 ***
## edad -0.0156 0.9845 0.0075 -2.08 0.037 *
## ede_factorUnilateral -0.2212 0.8016 0.1482 -1.49 0.136
## ede_factorBilateral -0.2359 0.7899 0.2338 -1.01 0.313
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.509 1.97 0.371 0.697
## edad 0.984 1.02 0.970 0.999
## ede_factorUnilateral 0.802 1.25 0.600 1.072
## ede_factorBilateral 0.790 1.27 0.500 1.249
##
## Concordance= 0.627 (se = 0.021 )
## Likelihood ratio test= 39.3 on 4 df, p=0.00000006
## Wald test = 36.5 on 4 df, p=0.0000002
## Score (logrank) test = 38 on 4 df, p=0.0000001tab_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.51 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.037 |
| ede factor [Unilateral] | 0.80 | 0.12 | 0.00 – Inf | 0.136 |
| ede factorBilateral | 0.79 | 0.18 | 0.00 – Inf | 0.313 |
| Observations | 390 | |||
| R2 Nagelkerke | 0.096 | |||
confint(cox_model_s)## 2.5 % 97.5 %
## isq_factorSI -0.9904 -0.360555
## edad -0.0303 -0.000936
## ede_factorUnilateral -0.5116 0.069236
## ede_factorBilateral -0.6940 0.222312# 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 3.25 1 0.071
## edad 0.15 1 0.698
## ede_factor 1.80 2 0.407
## GLOBAL 5.34 4 0.254#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 236 0
## 2 72 0
## unknown 0 82
##
##
## 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.9268 0.396 0.15999 -5.79 0.0000000069
## edad -0.0169 0.983 0.00755 -2.24 0.0250000000
## ede_factorUnilateral -0.3568 0.700 0.14719 -2.42 0.0150000000
## ede_factorBilateral -0.2792 0.756 0.23397 -1.19 0.2300000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.396 2.53 0.289 0.542
## edad 0.983 1.02 0.969 0.998
## ede_factorUnilateral 0.700 1.43 0.525 0.934
## ede_factorBilateral 0.756 1.32 0.478 1.196
##
## Num. cases = 390
## Pseudo Log-likelihood = -1224
## Pseudo likelihood ratio test = 70.4 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" )
Análisis de sensibilidad: todos los pacientes seguidos menos de 30 días se consideran cicatrizados
#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
## 58 252 14 82tiempo_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= 236
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.6761 0.5086 0.1607 -4.21 0.000026 ***
## edad -0.0156 0.9845 0.0075 -2.08 0.038 *
## ede_factorUnilateral -0.2211 0.8016 0.1482 -1.49 0.136
## ede_factorBilateral -0.2344 0.7911 0.2338 -1.00 0.316
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.509 1.97 0.371 0.697
## edad 0.985 1.02 0.970 0.999
## ede_factorUnilateral 0.802 1.25 0.600 1.072
## ede_factorBilateral 0.791 1.26 0.500 1.251
##
## Concordance= 0.627 (se = 0.021 )
## Likelihood ratio test= 39.3 on 4 df, p=0.00000006
## Wald test = 36.5 on 4 df, p=0.0000002
## Score (logrank) test = 38 on 4 df, p=0.0000001tab_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.51 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.038 |
| ede factor [Unilateral] | 0.80 | 0.12 | 0.00 – Inf | 0.136 |
| ede factorBilateral | 0.79 | 0.18 | 0.00 – Inf | 0.316 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.093 | |||
confint(cox_model_s1)## 2.5 % 97.5 %
## isq_factorSI -0.9910 -0.361092
## edad -0.0303 -0.000899
## ede_factorUnilateral -0.5115 0.069301
## ede_factorBilateral -0.6926 0.223788# 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 3.263 1 0.071
## edad 0.149 1 0.700
## ede_factor 1.823 2 0.402
## GLOBAL 5.381 4 0.250#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 252 0
## 2 72 0
## unknown 0 82
##
##
## 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.9101 0.402 0.15361 -5.92 0.0000000031
## edad -0.0165 0.984 0.00706 -2.34 0.0190000000
## ede_factorUnilateral -0.3171 0.728 0.14216 -2.23 0.0260000000
## ede_factorBilateral -0.2842 0.753 0.22585 -1.26 0.2100000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.402 2.48 0.298 0.544
## edad 0.984 1.02 0.970 0.997
## ede_factorUnilateral 0.728 1.37 0.551 0.962
## ede_factorBilateral 0.753 1.33 0.483 1.172
##
## Num. cases = 406
## Pseudo Log-likelihood = -1319
## Pseudo likelihood ratio test = 71.6 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
## 98 58 250tiempo_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= 58
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.5163 4.5554 0.3475 4.36 0.000013 ***
## edad 0.0225 1.0227 0.0138 1.63 0.102
## ede_factorUnilateral 0.6438 1.9038 0.2799 2.30 0.021 *
## ede_factorBilateral 0.1270 1.1354 0.4858 0.26 0.794
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 4.56 0.220 2.306 9.00
## edad 1.02 0.978 0.996 1.05
## ede_factorUnilateral 1.90 0.525 1.100 3.30
## ede_factorBilateral 1.14 0.881 0.438 2.94
##
## Concordance= 0.735 (se = 0.029 )
## Likelihood ratio test= 40.8 on 4 df, p=0.00000003
## Wald test = 33.4 on 4 df, p=0.000001
## Score (logrank) test = 39.9 on 4 df, p=0.00000004tab_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.56 | 1.58 | 0.00 – Inf | <0.001 |
| edad | 1.02 | 0.01 | 0.00 – Inf | 0.102 |
| ede factor [Unilateral] | 1.90 | 0.53 | 0.00 – Inf | 0.021 |
| ede factorBilateral | 1.14 | 0.55 | 0.00 – Inf | 0.794 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.119 | |||
confint(cox_model_a)## 2.5 % 97.5 %
## isq_factorSI 0.83532 2.1973
## edad -0.00448 0.0494
## ede_factorUnilateral 0.09521 1.1925
## ede_factorBilateral -0.82524 1.0792# 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 0.435 1 0.51
## edad 0.104 1 0.75
## ede_factor 2.601 2 0.27
## GLOBAL 3.433 4 0.49#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 58 0
## 2 250 0
## unknown 0 98
##
##
## 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.6325 5.12 0.3518 4.640 0.0000035
## edad 0.0204 1.02 0.0114 1.795 0.0730000
## ede_factorUnilateral 0.6436 1.90 0.2824 2.279 0.0230000
## ede_factorBilateral 0.1417 1.15 0.4696 0.302 0.7600000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 5.12 0.195 2.568 10.20
## edad 1.02 0.980 0.998 1.04
## ede_factorUnilateral 1.90 0.525 1.094 3.31
## ede_factorBilateral 1.15 0.868 0.459 2.89
##
## Num. cases = 406
## Pseudo Log-likelihood = -317
## Pseudo likelihood ratio test = 45.7 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
## 98 58 250#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 88 0.049 0.014 0.026, 0.082
## NO/leve 200 25 0.054 0.015 0.029, 0.089
## NO/leve 300 4 0.054 0.015 0.029, 0.089
## SI 100 75 0.243 0.034 0.180, 0.311
## SI 200 29 0.281 0.036 0.213, 0.353
## SI 300 9 0.290 0.037 0.221, 0.364## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 88 0.528 0.034 0.459, 0.592
## NO/leve 200 25 0.767 0.031 0.699, 0.821
## NO/leve 300 4 0.900 0.027 0.830, 0.942
## SI 100 75 0.218 0.033 0.156, 0.286
## SI 200 29 0.451 0.043 0.365, 0.532
## SI 300 9 0.575 0.047 0.477, 0.661## • Tests## outcome statistic df p.value
## 1 38.4 1.00 <0.001
## 2 49.7 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 4 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 38.4 0.00000000058891 1
## 2 49.7 0.00000000000181 1
## Estimates and Variances:
## $est
## 100 200 300
## NO/leve 1 0.0486 0.0539 0.0539
## SI 1 0.2431 0.2809 0.2905
## NO/leve 2 0.5278 0.7668 0.8996
## SI 2 0.2176 0.4508 0.5747
##
## $var
## 100 200 300
## NO/leve 1 0.000205 0.000232 0.000232
## SI 1 0.001132 0.001297 0.001355
## NO/leve 2 0.001155 0.000969 0.000748
## SI 2 0.001110 0.001839 0.002240infecció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
## 98 236 72#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 50 0
## 2 11 0
## unknown 0 25
##
##
## 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.0884 0.337 0.4521 -2.408 0.016
## edad -0.0191 0.981 0.0151 -1.266 0.210
## ede_factorUnilateral 0.1720 1.188 0.3806 0.452 0.650
## ede_factorBilateral -0.3670 0.693 0.3766 -0.974 0.330
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.337 2.970 0.139 0.817
## edad 0.981 1.019 0.953 1.011
## ede_factorUnilateral 1.188 0.842 0.563 2.504
## ede_factorBilateral 0.693 1.443 0.331 1.449
##
## Num. cases = 86
## Pseudo Log-likelihood = -179
## Pseudo likelihood ratio test = 18.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 156 0
## 2 36 0
## unknown 0 63
##
##
## 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.9354 0.392 0.17639 -5.303 0.00000011
## edad -0.0146 0.986 0.00916 -1.593 0.11000000
## ede_factorUnilateral 0.0432 1.044 0.17459 0.247 0.80000000
## ede_factorBilateral -0.0817 0.922 0.28361 -0.288 0.77000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.392 2.548 0.278 0.555
## edad 0.986 1.015 0.968 1.003
## ede_factorUnilateral 1.044 0.958 0.742 1.470
## ede_factorBilateral 0.922 1.085 0.529 1.607
##
## Num. cases = 255
## Pseudo Log-likelihood = -729
## Pseudo likelihood ratio test = 42.4 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 30 0
## 2 25 0
## unknown 0 10
##
##
## 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.5095 0.221 0.6402 -2.358 0.01800
## edad -0.0131 0.987 0.0267 -0.492 0.62000
## ede_factorUnilateral -1.5021 0.223 0.4144 -3.625 0.00029
## ede_factorBilateral -0.4494 0.638 1.2472 -0.360 0.72000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.221 4.52 0.0630 0.775
## edad 0.987 1.01 0.9367 1.040
## ede_factorUnilateral 0.223 4.49 0.0988 0.502
## ede_factorBilateral 0.638 1.57 0.0554 7.353
##
## Num. cases = 65
## Pseudo Log-likelihood = -101
## Pseudo likelihood ratio test = 23.9 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 13 0.643 0.071 0.487, 0.763
## NO/leve 200 7 0.701 0.071 0.538, 0.816
## NO/leve 250 4 0.789 0.075 0.594, 0.897
## NO/leve 300 1 0.963 0.091 0.006, 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 7 0.305 0.095 0.137, 0.493
## SI 250 3 0.305 0.095 0.137, 0.493
## SI 300 3 0.305 0.095 0.137, 0.493## • 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 13 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 7 0.312 0.090 0.151, 0.488
## SI 250 3 0.312 0.090 0.151, 0.488
## SI 300 3 0.312 0.090 0.151, 0.488## • Tests## outcome statistic df p.value
## 1 14.6 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" = "#56a", "SI" = "#E69F05")
) +
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 6 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 45 0.579 0.044 0.487, 0.659
## NO/leve 200 10 0.828 0.038 0.739, 0.889
## NO/leve 300 1 0.919 0.040 0.793, 0.970
## SI 100 51 0.209 0.041 0.135, 0.294
## SI 200 17 0.470 0.054 0.361, 0.571
## SI 300 5 0.606 0.058 0.483, 0.708## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 45 0.051 0.019 0.022, 0.096
## NO/leve 200 10 0.051 0.019 0.022, 0.096
## NO/leve 300 1 0.051 0.019 0.022, 0.096
## SI 100 51 0.230 0.042 0.154, 0.315
## SI 200 17 0.286 0.046 0.200, 0.377
## SI 300 5 0.286 0.046 0.200, 0.377## • Tests## outcome statistic df p.value
## 1 40.9 1.00 <0.001
## 2 22.3 1.00 <0.001ggcuminc(
cifivhxasleve,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO/leve" = "#56a", "SI" = "#E69F05")
) +
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 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 18 0.278 0.076 0.142, 0.431
## NO/leve 150 12 0.417 0.084 0.253, 0.573
## NO/leve 200 8 0.503 0.086 0.326, 0.656
## NO/leve 250 3 0.657 0.084 0.466, 0.794
## NO/leve 300 2 0.657 0.084 0.466, 0.794
## SI 50.0 16 0.000 0.000 NA, NA
## SI 100 11 0.040 0.040 0.003, 0.175
## SI 150 5 0.223 0.091 0.078, 0.415
## SI 200 5 0.223 0.091 0.078, 0.415
## SI 250 3 0.223 0.091 0.078, 0.415
## SI 300 1 0.353 0.115 0.146, 0.569## • 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 18 0.222 0.070 0.103, 0.370
## NO/leve 150 12 0.250 0.074 0.122, 0.401
## NO/leve 200 8 0.250 0.074 0.122, 0.401
## NO/leve 250 3 0.250 0.074 0.122, 0.401
## NO/leve 300 2 0.250 0.074 0.122, 0.401
## SI 50.0 16 0.448 0.095 0.261, 0.619
## SI 100 11 0.485 0.096 0.291, 0.654
## SI 150 5 0.534 0.100 0.324, 0.705
## SI 200 5 0.534 0.100 0.324, 0.705
## SI 250 3 0.582 0.103 0.358, 0.752
## SI 300 1 0.582 0.103 0.358, 0.752## • Tests## outcome statistic df p.value
## 1 6.02 1.00 0.014
## 2 6.94 1.00 0.008ggcuminc(
cifivhxasgrave,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO/leve" = "#56a", "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= 236
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factor2LEVE 0.02260 1.02286 0.22067 0.10 0.91843
## isq_factor2MODERADA -0.61232 0.54209 0.21908 -2.79 0.00519 **
## isq_factor2GRAVE -0.71444 0.48947 0.19619 -3.64 0.00027 ***
## edad -0.01525 0.98487 0.00757 -2.01 0.04410 *
## ede_factorUnilateral -0.21845 0.80377 0.15034 -1.45 0.14623
## ede_factorBilateral -0.23840 0.78789 0.23405 -1.02 0.30839
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factor2LEVE 1.023 0.978 0.664 1.576
## isq_factor2MODERADA 0.542 1.845 0.353 0.833
## isq_factor2GRAVE 0.489 2.043 0.333 0.719
## edad 0.985 1.015 0.970 1.000
## ede_factorUnilateral 0.804 1.244 0.599 1.079
## ede_factorBilateral 0.788 1.269 0.498 1.246
##
## Concordance= 0.627 (se = 0.021 )
## Likelihood ratio test= 39.5 on 6 df, p=0.0000006
## Wald test = 36.5 on 6 df, p=0.000002
## Score (logrank) test = 38.1 on 6 df, p=0.000001tab_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.02 | 0.23 | 0.00 – Inf | 0.918 |
| isq factor2MODERADA | 0.54 | 0.12 | 0.00 – Inf | 0.005 |
| isq factor2GRAVE | 0.49 | 0.10 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.044 |
| ede factor [Unilateral] | 0.80 | 0.12 | 0.00 – Inf | 0.146 |
| ede factorBilateral | 0.79 | 0.18 | 0.00 – Inf | 0.308 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.093 | |||
confint(cox_modelisq2)## 2.5 % 97.5 %
## isq_factor2LEVE -0.4099 0.455110
## isq_factor2MODERADA -1.0417 -0.182930
## isq_factor2GRAVE -1.0990 -0.329908
## edad -0.0301 -0.000403
## ede_factorUnilateral -0.5131 0.076220
## ede_factorBilateral -0.6971 0.220323#
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 74 0.506 0.037 0.433, 0.576
## NO 200 23 0.725 0.035 0.650, 0.787
## NO 300 4 0.856 0.033 0.776, 0.909
## LEVE 100 14 0.472 0.095 0.281, 0.641
## LEVE 200 2 0.830 0.082 0.592, 0.936
## LEVE 300 0 0.968 0.076 0.026, 1.00
## MODERADA 100 25 0.231 0.065 0.118, 0.366
## MODERADA 200 13 0.484 0.080 0.322, 0.628
## MODERADA 300 3 0.720 0.084 0.518, 0.848
## GRAVE 100 50 0.161 0.035 0.099, 0.235
## GRAVE 200 16 0.369 0.051 0.271, 0.467
## GRAVE 300 6 0.429 0.056 0.318, 0.536## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO 100 74 0.082 0.020 0.049, 0.126
## NO 200 23 0.088 0.021 0.053, 0.133
## NO 300 4 0.088 0.021 0.053, 0.133
## LEVE 100 14 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 25 0.157 0.055 0.068, 0.280
## MODERADA 200 13 0.157 0.055 0.068, 0.280
## MODERADA 300 3 0.157 0.055 0.068, 0.280
## GRAVE 100 50 0.328 0.044 0.245, 0.414
## GRAVE 200 16 0.404 0.047 0.311, 0.495
## GRAVE 300 6 0.418 0.048 0.323, 0.511## • Tests## outcome statistic df p.value
## 1 57.1 3.00 <0.001
## 2 52.2 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 42 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 4.556 3 0.21
## edad 0.162 1 0.69
## ede_factor 1.871 2 0.39
## GLOBAL 6.758 6 0.34#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 236 0
## 2 72 0
## unknown 0 98
##
##
## 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.0946 1.099 0.19578 0.483 0.630000000
## isq_factor2MODERADA -0.5544 0.574 0.19569 -2.833 0.004600000
## isq_factor2GRAVE -1.1171 0.327 0.20060 -5.569 0.000000026
## edad -0.0152 0.985 0.00767 -1.988 0.047000000
## ede_factorUnilateral -0.3524 0.703 0.14942 -2.358 0.018000000
## ede_factorBilateral -0.2993 0.741 0.23589 -1.269 0.200000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factor2LEVE 1.099 0.91 0.749 1.613
## isq_factor2MODERADA 0.574 1.74 0.391 0.843
## isq_factor2GRAVE 0.327 3.06 0.221 0.485
## edad 0.985 1.02 0.970 1.000
## ede_factorUnilateral 0.703 1.42 0.525 0.942
## ede_factorBilateral 0.741 1.35 0.467 1.177
##
## Num. cases = 406
## Pseudo Log-likelihood = -1221
## Pseudo likelihood ratio test = 75.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 236 0
## 2 72 0
## unknown 0 98
##
##
## 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.0157 0.362 0.16348 -6.21 0.00000000052
## edad -0.0156 0.985 0.00767 -2.03 0.04200000000
## ede_factorUnilateral -0.3588 0.698 0.14745 -2.43 0.01500000000
## ede_factorBilateral -0.2834 0.753 0.23679 -1.20 0.23000000000
## loca_factorMet/Tar -0.3748 0.687 0.13493 -2.78 0.00550000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.362 2.76 0.263 0.499
## edad 0.985 1.02 0.970 0.999
## ede_factorUnilateral 0.698 1.43 0.523 0.933
## ede_factorBilateral 0.753 1.33 0.474 1.198
## loca_factorMet/Tar 0.687 1.45 0.528 0.896
##
## Num. cases = 406
## Pseudo Log-likelihood = -1219
## Pseudo likelihood ratio test = 77.9 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 236 0
## 2 72 0
## unknown 0 98
##
##
## 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.8390 0.432 0.15696 -5.35 0.00000009
## edad -0.0165 0.984 0.00721 -2.29 0.02200000
## ede_factorUnilateral -0.2289 0.795 0.14773 -1.55 0.12000000
## ede_factorBilateral -0.2940 0.745 0.24033 -1.22 0.22000000
## topo_factorlat/Mas de dos -0.5055 0.603 0.13627 -3.71 0.00021000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.432 2.31 0.318 0.588
## edad 0.984 1.02 0.970 0.998
## ede_factorUnilateral 0.795 1.26 0.595 1.062
## ede_factorBilateral 0.745 1.34 0.465 1.194
## topo_factorlat/Mas de dos 0.603 1.66 0.462 0.788
##
## Num. cases = 406
## Pseudo Log-likelihood = -1217
## Pseudo likelihood ratio test = 83.5 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 236 0
## 2 72 0
## unknown 0 98
##
##
## 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.8479 0.428 0.15714 -5.396 0.000000068
## edad -0.0168 0.983 0.00755 -2.227 0.026000000
## ede_factorUnilateral -0.1148 0.892 0.15764 -0.728 0.470000000
## ede_factorBilateral -0.2972 0.743 0.24798 -1.198 0.230000000
## zon_factorDos o mas -0.7700 0.463 0.15194 -5.068 0.000000400
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.428 2.33 0.315 0.583
## edad 0.983 1.02 0.969 0.998
## ede_factorUnilateral 0.892 1.12 0.655 1.214
## ede_factorBilateral 0.743 1.35 0.457 1.208
## zon_factorDos o mas 0.463 2.16 0.344 0.624
##
## Num. cases = 406
## Pseudo Log-likelihood = -1209
## Pseudo likelihood ratio test = 97.7 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 236 0
## 2 72 0
## unknown 0 98
##
##
## 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.8645 0.421 0.1650 -5.24 0.00000016
## edad -0.0179 0.982 0.0075 -2.39 0.01700000
## ede_factorUnilateral -0.3773 0.686 0.1466 -2.57 0.01000000
## ede_factorBilateral -0.3017 0.740 0.2358 -1.28 0.20000000
## neur_factorMod/charcot 0.4788 1.614 0.3624 1.32 0.19000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.421 2.37 0.305 0.582
## edad 0.982 1.02 0.968 0.997
## ede_factorUnilateral 0.686 1.46 0.514 0.914
## ede_factorBilateral 0.740 1.35 0.466 1.174
## neur_factorMod/charcot 1.614 0.62 0.793 3.284
##
## Num. cases = 406
## Pseudo Log-likelihood = -1222
## Pseudo likelihood ratio test = 72.4 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 236 0
## 2 72 0
## unknown 0 98
##
##
## 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.907 0.404 0.15942 -5.69 0.000000013
## edad -0.016 0.984 0.00763 -2.10 0.035000000
## ede_factorUnilateral -0.327 0.721 0.14749 -2.21 0.027000000
## ede_factorBilateral -0.265 0.767 0.22811 -1.16 0.240000000
## prof_factorparcial/total -0.409 0.665 0.30045 -1.36 0.170000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.404 2.48 0.295 0.552
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.721 1.39 0.540 0.963
## ede_factorBilateral 0.767 1.30 0.491 1.199
## prof_factorparcial/total 0.665 1.50 0.369 1.198
##
## Num. cases = 406
## Pseudo Log-likelihood = -1222
## Pseudo likelihood ratio test = 73 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 236 0
## 2 72 0
## unknown 0 98
##
##
## 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.8960 0.408 0.15958 -5.615 0.00000002000
## edad -0.0171 0.983 0.00773 -2.206 0.02700000000
## ede_factorUnilateral -0.0358 0.965 0.15221 -0.235 0.81000000000
## ede_factorBilateral -0.3379 0.713 0.24922 -1.356 0.18000000000
## area_factor10 a 40 -0.3376 0.713 0.13998 -2.412 0.01600000000
## area_factorMas de 40 -1.5685 0.208 0.25161 -6.234 0.00000000046
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.408 2.45 0.299 0.558
## edad 0.983 1.02 0.968 0.998
## ede_factorUnilateral 0.965 1.04 0.716 1.300
## ede_factorBilateral 0.713 1.40 0.438 1.163
## area_factor10 a 40 0.713 1.40 0.542 0.939
## area_factorMas de 40 0.208 4.80 0.127 0.341
##
## Num. cases = 406
## Pseudo Log-likelihood = -1201
## Pseudo likelihood ratio test = 115 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 236 0
## 2 72 0
## unknown 0 98
##
##
## 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.8795 0.415 0.16113 -5.46 0.000000048
## edad -0.0183 0.982 0.00749 -2.44 0.014000000
## ede_factorUnilateral -0.3372 0.714 0.14910 -2.26 0.024000000
## ede_factorBilateral -0.2542 0.776 0.23516 -1.08 0.280000000
## fase_factorGranu/infla -0.7832 0.457 0.37916 -2.07 0.039000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.415 2.41 0.303 0.569
## edad 0.982 1.02 0.968 0.996
## ede_factorUnilateral 0.714 1.40 0.533 0.956
## ede_factorBilateral 0.776 1.29 0.489 1.230
## fase_factorGranu/infla 0.457 2.19 0.217 0.961
##
## Num. cases = 406
## Pseudo Log-likelihood = -1220
## Pseudo likelihood ratio test = 76.4 on 5 df,
##
## Convergence: TRUE