Tesis 17 febrero
Carro Gabriela
2026-02-18
#CARGA DE LIBRERIAS
knitr::opts_chunk$set(echo = TRUE)
library(readr)
library(data.table)## Warning: package 'data.table' was built under R version 4.4.3library(tidyverse)## Warning: package 'ggplot2' was built under R version 4.4.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ purrr 1.0.2
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 4.0.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::between() masks data.table::between()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::first() masks data.table::first()
## ✖ lubridate::hour() masks data.table::hour()
## ✖ lubridate::isoweek() masks data.table::isoweek()
## ✖ dplyr::lag() masks stats::lag()
## ✖ dplyr::last() masks data.table::last()
## ✖ lubridate::mday() masks data.table::mday()
## ✖ lubridate::minute() masks data.table::minute()
## ✖ lubridate::month() masks data.table::month()
## ✖ lubridate::quarter() masks data.table::quarter()
## ✖ lubridate::second() masks data.table::second()
## ✖ purrr::transpose() masks data.table::transpose()
## ✖ lubridate::wday() masks data.table::wday()
## ✖ lubridate::week() masks data.table::week()
## ✖ lubridate::yday() masks data.table::yday()
## ✖ lubridate::year() masks data.table::year()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(broom)
library(vcd)## Warning: package 'vcd' was built under R version 4.4.3## Cargando paquete requerido: gridlibrary(sjPlot)##
## Adjuntando el paquete: 'sjPlot'
##
## The following object is masked from 'package:ggplot2':
##
## set_themelibrary(epiR)## Cargando paquete requerido: survival
## Package epiR 2.0.77 is loaded
## Type help(epi.about) for summary information
## Type browseVignettes(package = 'epiR') to learn how to use epiR for applied epidemiological analyseslibrary(gtsummary)## Warning: package 'gtsummary' was built under R version 4.4.3library(tableone)
library(ggpubr)
library(vcd)
library(vcdExtra)## Warning: package 'vcdExtra' was built under R version 4.4.3## Cargando paquete requerido: gnm## Warning: package 'gnm' was built under R version 4.4.3##
## Adjuntando el paquete: 'vcdExtra'
##
## The following object is masked from 'package:dplyr':
##
## summariselibrary(CalibrationCurves)## Warning: package 'CalibrationCurves' was built under R version 4.4.3## Cargando paquete requerido: rms## Warning: package 'rms' was built under R version 4.4.3## Cargando paquete requerido: Hmisc## Warning: package 'Hmisc' was built under R version 4.4.3##
## Adjuntando el paquete: 'Hmisc'
##
## The following objects are masked from 'package:dplyr':
##
## src, summarize
##
## The following objects are masked from 'package:base':
##
## format.pval, unitslibrary(pROC)## Warning: package 'pROC' was built under R version 4.4.3## Type 'citation("pROC")' for a citation.
##
## Adjuntando el paquete: 'pROC'
##
## The following objects are masked from 'package:stats':
##
## cov, smooth, varlibrary(glmnet)## Warning: package 'glmnet' was built under R version 4.4.3## Cargando paquete requerido: Matrix
##
## Adjuntando el paquete: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loaded glmnet 4.1-10library(caret)## Warning: package 'caret' was built under R version 4.4.3## Cargando paquete requerido: lattice
##
## Adjuntando el paquete: 'lattice'
##
## The following object is masked from 'package:gnm':
##
## barley
##
##
## Adjuntando el paquete: 'caret'
##
## The following object is masked from 'package:survival':
##
## cluster
##
## The following object is masked from 'package:purrr':
##
## liftlibrary(broom.helpers)## Warning: package 'broom.helpers' was built under R version 4.4.3##
## Adjuntando el paquete: 'broom.helpers'
##
## The following objects are masked from 'package:gtsummary':
##
## all_categorical, all_continuous, all_contrasts, all_dichotomous,
## all_interaction, all_interceptsoptions(scipen = 999, digits = 3, encoding = 'UTF-8')
library(naniar)## Warning: package 'naniar' was built under R version 4.4.3library(dplyr)
library(survival)
library(survminer)##
## Adjuntando el paquete: 'survminer'
##
## The following object is masked from 'package:survival':
##
## myelomalibrary(MASS)##
## Adjuntando el paquete: 'MASS'
##
## The following object is masked from 'package:gtsummary':
##
## select
##
## The following object is masked from 'package:dplyr':
##
## selectlibrary(rms)
library(sjPlot)
pkgbuild::has_build_tools(debug = TRUE)## Trying to compile a simple C file## Running "C:/PROGRA~1/R/R-44~1.2/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 13.3.0'
## gcc -I"C:/PROGRA~1/R/R-44~1.2/include" -DNDEBUG -I"C:/rtools44/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -mfpmath=sse -msse2 -mstackrealign -c foo.c -o foo.o
## gcc -shared -s -static-libgcc -o foo.dll tmp.def foo.o -LC:/rtools44/x86_64-w64-mingw32.static.posix/lib/x64 -LC:/rtools44/x86_64-w64-mingw32.static.posix/lib -LC:/PROGRA~1/R/R-44~1.2/bin/x64 -lR## ## [1] TRUElibrary(ggsurvfit)## Warning: package 'ggsurvfit' was built under R version 4.4.3library(riskRegression)## Warning: package 'riskRegression' was built under R version 4.4.3## Registered S3 method overwritten by 'riskRegression':
## method from
## nobs.multinom broom
## riskRegression version 2025.05.20library(rms)
library(cmprsk)## Warning: package 'cmprsk' was built under R version 4.4.3library(prodlim)## Warning: package 'prodlim' was built under R version 4.4.3library(naniar)
library(ggplot2)
library(scales)## Warning: package 'scales' was built under R version 4.4.3##
## Adjuntando el paquete: 'scales'
##
## The following object is masked from 'package:purrr':
##
## discard
##
## The following object is masked from 'package:readr':
##
## col_factorlibrary(tidycmprsk)## Warning: package 'tidycmprsk' was built under R version 4.4.3##
## Adjuntando el paquete: 'tidycmprsk'
##
## The following objects are masked from 'package:cmprsk':
##
## crr, cuminc
##
## The following object is masked from 'package:gtsummary':
##
## triallibrary(dplyr)
library(psych)## Warning: package 'psych' was built under R version 4.4.3##
## Adjuntando el paquete: 'psych'
##
## The following objects are masked from 'package:scales':
##
## alpha, rescale
##
## The following object is masked from 'package:ggsurvfit':
##
## %+%
##
## The following object is masked from 'package:Hmisc':
##
## describe
##
## The following objects are masked from 'package:ggplot2':
##
## %+%, alphaCARGA DE LA BASE DE DATOS
#carga de datos
library(readxl)
sob1 <- read_excel("C:/Users/Administrador/Downloads/TESI.xlsx")## New names:
## • `` -> `...44`View(sob1)names(sob1)<- tolower(names (sob1)) #PASA A MINUSCULA LOS NOMBRES
#PROPORCION CICATRIZADOS
table1=table(sob1$cicat)
table1##
## 0 1
## 183 223prop.table(table1)##
## 0 1
## 0.451 0.549#PROPORCION AMPUTADOS
table2=table(sob1$ampu)
table2##
## 0 1
## 349 57prop.table(table2)##
## 0 1
## 0.86 0.14#PROPORCION FALLECIDOS
table3=table(sob1$fallecio)
table3##
## 0 1
## 392 14prop.table(table3)##
## 0 1
## 0.9655 0.0345#tabla por sexo
table4=table(sob1$sexo)
table4##
## F M
## 82 324prop.table(table4)##
## F M
## 0.202 0.798#promedio de edad
describe(sob1$edad)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 406 57.4 11.1 57 57.4 10.4 22 89 67 0 0.35 0.55#dias de seguimiento
describe(sob1$dias_out)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 406 96.2 75.5 77.5 86.4 64.5 1 362 361 1.15 0.94 3.75Descripción de los eventos
#pasar variables categóricas a factor
sob1 <- sob1 %>%
mutate(across(19:35, as.factor))
#crear cicat_factor
sob1$cicat_factor <- factor(sob1$cicat, levels = c(0, 1), labels = c( "No Cicatrizo", "Cicatrizo"))
table(sob1$cicat_factor)##
## No Cicatrizo Cicatrizo
## 183 223#paso variables chr a factor
sob1 <- sob1 %>% mutate_if(is.character, as.factor)
#crear ampu_factor
sob1$ampu_factor <- factor(sob1$ampu, levels = c(0, 1), labels = c( "No Amputado", "Amputado"))
table(sob1$ampu_factor)##
## No Amputado Amputado
## 349 57#crear fallecio_factor
sob1$fallecio_factor <- factor(sob1$fallecio, levels = c(0, 1), labels = c( "No Fallecido", "Fallecido"))
table(sob1$fallecio_factor)##
## No Fallecido Fallecido
## 392 14#variable unificada
sob1 <- sob1 %>%
mutate(
estado = tolower(trimws(estado))
)
table(sob1$estado)##
## ampu cerro fallecio persiste
## 57 223 14 112sob1 <- sob1 %>%
mutate(
estado_unif = case_when(
estado %in% c("ampu") ~ "amputacion",
estado %in% c("cerro") ~ "cerro_herida",
estado %in% c("murio", "fallecio") ~ "fallecio",
estado %in% c("persiste") ~ "persiste",
estado %in% c("perdido") ~ "perdido_seguimiento",
TRUE ~ NA_character_
)
)
table(sob1$estado_unif)##
## amputacion cerro_herida fallecio persiste
## 57 223 14 112Intervalo de confianza de los eventos
#intervalos de confianza
tabla_estado <- sob1 %>%
dplyr::filter(!is.na(estado_unif)) %>%
dplyr::group_by(estado_unif) %>%
dplyr::summarise(
n = n()
) %>%
dplyr::mutate(
N_total = sum(n),
prop = n / N_total,
se = sqrt(prop * (1 - prop) / N_total),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
tabla_estado %>%
dplyr::mutate(
porcentaje = round(prop * 100, 1),
IC95 = paste0(
round(li * 100, 1), "–",
round(ls * 100, 1), "%"
)
)## # A tibble: 4 × 9
## estado_unif n N_total prop se li ls porcentaje IC95
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 amputacion 57 406 0.140 0.0172 0.107 0.174 14 10.7–17.4%
## 2 cerro_herida 223 406 0.549 0.0247 0.501 0.598 54.9 50.1–59.8%
## 3 fallecio 14 406 0.0345 0.00906 0.0167 0.0522 3.4 1.7–5.2%
## 4 persiste 112 406 0.276 0.0222 0.232 0.319 27.6 23.2–31.9%tabla_estado## # A tibble: 4 × 7
## estado_unif n N_total prop se li ls
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl>
## 1 amputacion 57 406 0.140 0.0172 0.107 0.174
## 2 cerro_herida 223 406 0.549 0.0247 0.501 0.598
## 3 fallecio 14 406 0.0345 0.00906 0.0167 0.0522
## 4 persiste 112 406 0.276 0.0222 0.232 0.319#dias de seguimiento
#en pacientes censurados
#dias de seguimiento
describe(sob1$dias_out)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 406 96.2 75.5 77.5 86.4 64.5 1 362 361 1.15 0.94 3.75describeBy(
x = sob1$dias_out,
group = sob1$estado_unif,
mat = TRUE
)## item group1 vars n mean sd median trimmed mad min max range
## X11 1 amputacion 1 57 39.6 39.9 22 32.8 16.3 1 211 210
## X12 2 cerro_herida 1 223 99.7 70.5 82 89.7 57.8 9 330 321
## X13 3 fallecio 1 14 43.1 40.4 24 39.9 25.2 5 120 115
## X14 4 persiste 1 112 124.7 83.7 114 117.6 83.8 6 362 356
## skew kurtosis se
## X11 1.982 4.5398 5.28
## X12 1.234 1.0772 4.72
## X13 0.661 -1.2525 10.80
## X14 0.718 0.0353 7.91tapply(sob1$dias_out,
sob1$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 14 48
##
## $cerro_herida
## 25% 75%
## 46 128
##
## $fallecio
## 25% 75%
## 11.2 69.8
##
## $persiste
## 25% 75%
## 63.8 175.8#ISQuEMIA
#isq_factor
sob1 <- sob1 %>%
mutate(isq_factor = factor(
if_else(isq %in% c(2, 3), "SI", "NO/leve")
))
#dias seguimiento en subgrupos
sob1_per_isq <- subset(sob1, isq_factor == "SI")
sob1_per_noisq <- subset(sob1, isq_factor == "NO/leve")
#no isquemicos persistente seguimiento
describe(sob1_per_noisq$dias_out)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 237 91.1 68.2 77 82.7 62.3 5 337 332 1.11 0.88 4.43describeBy(
x = sob1_per_noisq$dias_out,
group = sob1_per_noisq$estado_unif,
mat = TRUE
)## item group1 vars n mean sd median trimmed mad min max range
## X11 1 amputacion 1 11 26.5 26.1 19.0 20.8 11.9 7 98 91
## X12 2 cerro_herida 1 162 90.6 63.9 71.5 80.9 52.6 15 301 286
## X13 3 fallecio 1 6 25.8 23.1 17.5 25.8 15.6 5 66 61
## X14 4 persiste 1 58 111.7 76.6 112.5 107.2 86.7 6 337 331
## skew kurtosis se
## X11 1.789 2.182 7.88
## X12 1.313 1.319 5.02
## X13 0.708 -1.298 9.44
## X14 0.532 -0.056 10.06tapply(sob1_per_noisq$dias_out,
sob1_per_noisq$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 12.5 26.0
##
## $cerro_herida
## 25% 75%
## 42 116
##
## $fallecio
## 25% 75%
## 10.8 34.8
##
## $persiste
## 25% 75%
## 50.8 167.0#pacientes isquémicos
# isquemicos persistente seguimiento
describe(sob1_per_isq$dias_out)## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 169 103 84.3 82 92.2 80.1 1 362 361 1.07 0.52 6.49describeBy(
x = sob1_per_isq$dias_out,
group = sob1_per_isq$estado_unif,
mat = TRUE
)## item group1 vars n mean sd median trimmed mad min max range
## X11 1 amputacion 1 46 42.7 42.1 24.5 35.9 20.8 1 211 210
## X12 2 cerro_herida 1 61 124.1 81.2 108.0 114.9 72.6 9 330 321
## X13 3 fallecio 1 8 56.1 46.9 50.0 56.1 61.5 6 120 114
## X14 4 persiste 1 54 138.7 89.3 116.5 130.2 79.3 7 362 355
## skew kurtosis se
## X11 1.847 3.8060 6.21
## X12 0.896 0.0851 10.40
## X13 0.134 -1.9918 16.58
## X14 0.743 -0.3552 12.15tapply(sob1_per_isq$dias_out,
sob1_per_isq$estado_unif,
quantile,
probs = c(0.25, 0.75),
na.rm = TRUE)## $amputacion
## 25% 75%
## 14.2 59.0
##
## $cerro_herida
## 25% 75%
## 63 160
##
## $fallecio
## 25% 75%
## 11.8 99.0
##
## $persiste
## 25% 75%
## 67.2 199.2#comparar dias de seguimiento en pacientes con y sin isquemia cuyo estado sea "persiste"
# Filtrar solo pacientes con estado "persiste"
sob_persiste <- subset(sob1, estado_unif == "persiste")
# Aplicar Wilcoxon
wilcox.test(dias_out ~ isq_factor, data = sob_persiste)##
## Wilcoxon rank sum test with continuity correction
##
## data: dias_out by isq_factor
## W = 1330, p-value = 0.2
## alternative hypothesis: true location shift is not equal to 0#Densidad de incidencia y razón de densidades
#cicatrizacion totaldensidad de incidencia
sob1$evento_cierre <- ifelse(sob1$estado_unif == "cerro_herida", 1, 0)
eventos <- sum(sob1$evento_cierre, na.rm = TRUE) #nro eventos
tiempo1_total <- sum(sob1$dias_out, na.rm = TRUE) #tiempo persona total
tasa1_dia <- eventos / tiempo1_total
tasa1_dia## [1] 0.00571tiempo1_meses <- tiempo1_total / 30.44
tasa1_mensual <- eventos / tiempo1_meses
tasa1_mensual## [1] 0.174tasa1_1000 <- tasa1_dia * 1000
tasa1_1000## [1] 5.71#cicatrizacion isq tiempo por 1000 personas por dia
sob1_per_isq$evento_cierre <- ifelse(sob1_per_isq$estado_unif == "cerro_herida", 1, 0)
eventos <- sum(sob1_per_isq$evento_cierre, na.rm = TRUE) #nro eventos
tiempo2_total <- sum(sob1_per_isq$dias_out, na.rm = TRUE) #tiempo persona total
tasa2_dia <- eventos / tiempo2_total
tasa2_dia## [1] 0.00349tiempo2_meses <- tiempo2_total / 30.44
tasa2_mensual <- eventos / tiempo2_meses
tasa2_mensual## [1] 0.106tasa2_1000 <- tasa2_dia * 1000
tasa2_1000## [1] 3.49#cicatrizacion no isq
sob1_per_noisq$evento_cierre <- ifelse(sob1_per_noisq$estado_unif == "cerro_herida", 1, 0)
eventos <- sum(sob1_per_noisq$evento_cierre, na.rm = TRUE) #nro eventos
tiempo3_total <- sum(sob1_per_noisq$dias_out, na.rm = TRUE) #tiempo persona total
tasa3_dia <- eventos / tiempo3_total
tasa3_dia## [1] 0.0075tiempo3_meses <- tiempo3_total / 30.44
tasa3_mensual <- eventos / tiempo3_meses
tasa3_mensual## [1] 0.228tasa3_1000 <- tasa3_dia * 1000
tasa3_1000## [1] 7.5# densidad de incidencia amputacion total
sob1$evento_ampu <- ifelse(sob1$estado_unif == "amputacion", 1, 0)
eventos <- sum(sob1$evento_ampu, na.rm = TRUE) #nro eventos
tiempo1_total <- sum(sob1$dias_out, na.rm = TRUE) #tiempo persona total
tasa1_dia <- eventos / tiempo1_total
tasa1_dia## [1] 0.00146tiempo1_meses <- tiempo1_total / 30.44
tasa1_mensual <- eventos / tiempo1_meses
tasa1_mensual## [1] 0.0444tasa1_1000 <- tasa1_dia * 1000
tasa1_1000## [1] 1.46#amputacuin isq
sob1_per_isq$evento_ampu <- ifelse(sob1_per_isq$estado_unif == "amputacion", 1, 0)
eventos <- sum(sob1_per_isq$evento_ampu, na.rm = TRUE) #nro eventos
tiempo2_total <- sum(sob1_per_isq$dias_out, na.rm = TRUE) #tiempo persona total
tasa2_dia <- eventos / tiempo2_total
tasa2_dia## [1] 0.00263tiempo2_meses <- tiempo2_total / 30.44
tasa2_mensual <- eventos / tiempo2_meses
tasa2_mensual## [1] 0.0801tasa2_1000 <- tasa2_dia * 1000
tasa2_1000## [1] 2.63#amputacion no isq
sob1_per_noisq$evento_ampu <- ifelse(sob1_per_noisq$estado_unif == "amputacion", 1, 0)
eventos <- sum(sob1_per_noisq$evento_ampu, na.rm = TRUE) #nro eventos
tiempo3_total <- sum(sob1_per_noisq$dias_out, na.rm = TRUE) #tiempo persona total
tasa3_dia <- eventos / tiempo3_total
tasa3_dia## [1] 0.000509tiempo3_meses <- tiempo3_total / 30.44
tasa3_mensual <- eventos / tiempo3_meses
tasa3_mensual## [1] 0.0155tasa3_1000 <- tasa3_dia * 1000
tasa3_1000## [1] 0.509#evento muerte densidad de incidencia
#cicatrizacion totaldensidad de incidencia
sob1$evento_muerte <- ifelse(sob1$estado_unif == "fallecio", 1, 0)
eventos <- sum(sob1$evento_muerte, na.rm = TRUE) #nro eventos
tiempo1_total <- sum(sob1$dias_out, na.rm = TRUE) #tiempo persona total
tasa1_dia <- eventos / tiempo1_total
tasa1_dia## [1] 0.000358tiempo1_meses <- tiempo1_total / 30.44
tasa1_mensual <- eventos / tiempo1_meses
tasa1_mensual## [1] 0.0109tasa1_1000 <- tasa1_dia * 1000
tasa1_1000## [1] 0.358#muerte isq tiempo por 1000 personas por dia. Tasa de incidencia de muerte por 1000 persons por dia
sob1_per_isq$evento_muerte <- ifelse(sob1_per_isq$estado_unif == "fallecio", 1, 0)
eventos <- sum(sob1_per_isq$evento_muerte, na.rm = TRUE) #nro eventos
tiempo2_total <- sum(sob1_per_isq$dias_out, na.rm = TRUE) #tiempo persona total
tasa2_dia <- eventos / tiempo2_total
tasa2_dia## [1] 0.000458tiempo2_meses <- tiempo2_total / 30.44
tasa2_mensual <- eventos / tiempo2_meses
tasa2_mensual## [1] 0.0139tasa2_1000 <- tasa2_dia * 1000
tasa2_1000## [1] 0.458#muerte no isq
sob1_per_noisq$evento_muerte <- ifelse(sob1_per_noisq$estado_unif == "fallecio", 1, 0)
eventos <- sum(sob1_per_noisq$evento_muerte, na.rm = TRUE) #nro eventos
tiempo3_total <- sum(sob1_per_noisq$dias_out, na.rm = TRUE) #tiempo persona total
tasa3_dia <- eventos / tiempo3_total
tasa3_dia## [1] 0.000278tiempo3_meses <- tiempo3_total / 30.44
tasa3_mensual <- eventos / tiempo3_meses
tasa3_mensual## [1] 0.00846tasa3_1000 <- tasa3_dia * 1000
tasa3_1000## [1] 0.278COMPORTAMIENTO DE LAS VARIABLES
#comportamiento de las variables
#LOCAl
library(dplyr)
tabla_prop <- sob1 %>%
dplyr::group_by(local) %>%
dplyr::summarise(
n = sum(!is.na(cicat)),
eventos = sum(cicat == 1, na.rm = TRUE),
prop = eventos / n,
se = sqrt(prop * (1 - prop) / n),
li = prop - 1.96 * se,
ls = prop + 1.96 * se
)
library(ggplot2)
ggplot(tabla_prop, aes(x = local, y = prop)) +
geom_col(fill = "#4E79A7", width = 0.6) +
geom_errorbar(aes(ymin = li, ymax = ls),
width = 0.2,
size = 0.8) +
geom_text(aes(label = paste0(round(prop*100,1), "%")),
vjust = -0.5,
size = 5) +
scale_y_continuous(limits = c(0,1),
labels = scales::percent_format()) +
labs(
x = "Localización de la lesión",
y = "Proporción de cicatrización (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 lesión",
y = "Proporcion de cicatrización (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 infección 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 cicatrizacion en no isq/leve
table(sob1_per_noisq$cicat_factor)##
## No Cicatrizo Cicatrizo
## 75 162xci <- sum(sob1_per_noisq$cicat_factor %in% c("Cicatrizo"), na.rm = TRUE)
nci <- sum(!is.na(sob1_per_noisq$cicat_factor))
prop.test(xci, nci)##
## 1-sample proportions test with continuity correction
##
## data: xci out of nci, null probability 0.5
## X-squared = 31, df = 1, p-value = 0.00000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.620 0.741
## sample estimates:
## p
## 0.684#porcentaje de amputacion en no isq/leve
table(sob1_per_noisq$ampu_factor)##
## No Amputado Amputado
## 226 11xa <- sum(sob1_per_noisq$ampu_factor %in% c("Amputado"), na.rm = TRUE)
na <- sum(!is.na(sob1_per_noisq$ampu_factor))
prop.test(xa, na)##
## 1-sample proportions test with continuity correction
##
## data: xa out of na, null probability 0.5
## X-squared = 193, df = 1, p-value <0.0000000000000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.0246 0.0838
## sample estimates:
## p
## 0.0464#porcentaje de muerte en no isq/leve
table(sob1_per_noisq$fallecio_factor)##
## No Fallecido Fallecido
## 231 6xf <- sum(sob1_per_noisq$fallecio_factor %in% c("Fallecido"), na.rm = TRUE)
nf <- sum(!is.na(sob1_per_noisq$fallecio_factor))
prop.test(xf, nf)##
## 1-sample proportions test with continuity correction
##
## data: xf out of nf, null probability 0.5
## X-squared = 212, df = 1, p-value <0.0000000000000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.0103 0.0569
## sample estimates:
## p
## 0.0253#porcentaje de cicatrizacion en isquemicos
table(sob1_per_isq$cicat_factor)##
## No Cicatrizo Cicatrizo
## 108 61xcii <- sum(sob1_per_isq$cicat_factor %in% c("Cicatrizo"), na.rm = TRUE)
ncii <- sum(!is.na(sob1_per_isq$cicat_factor))
prop.test(xcii, ncii)##
## 1-sample proportions test with continuity correction
##
## data: xcii out of ncii, null probability 0.5
## X-squared = 13, df = 1, p-value = 0.0004
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.290 0.439
## sample estimates:
## p
## 0.361#porcentaje de amputacion en isquemicos
table(sob1_per_isq$ampu_factor)##
## No Amputado Amputado
## 123 46xai <- sum(sob1_per_isq$ampu_factor %in% c("Amputado"), na.rm = TRUE)
nai <- sum(!is.na(sob1_per_isq$ampu_factor))
prop.test(xai, nai)##
## 1-sample proportions test with continuity correction
##
## data: xai out of nai, null probability 0.5
## X-squared = 34, df = 1, p-value = 0.000000005
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.208 0.347
## sample estimates:
## p
## 0.272#porcentaje de muerte en isquemicos
table(sob1_per_isq$fallecio_factor)##
## No Fallecido Fallecido
## 161 8xfi <- sum(sob1_per_isq$fallecio_factor %in% c("Fallecido"), na.rm = TRUE)
nfi<- sum(!is.na(sob1_per_isq$fallecio_factor))
prop.test(xfi, nfi)##
## 1-sample proportions test with continuity correction
##
## data: xfi out of nfi, null probability 0.5
## X-squared = 137, df = 1, p-value <0.0000000000000002
## alternative hypothesis: true p is not equal to 0.5
## 95 percent confidence interval:
## 0.0222 0.0944
## sample estimates:
## p
## 0.0473#TABLA 1
library(tableone)
#Las variables van directamente sin el nombre del dataframe
#catvars = c( "sexo","tbq", "iam", "hta", "dial", "local","topo", "inf", "ede", "neu", "prof", "area", "zon", "fase")
#continuas
#vars = c("dias_ev","san_elian", "sexo","tbq", "iam", "hta", "dial", "local","topo", "inf", "ede", "neu", "prof", "area", "zon", "fase" )
#tabla_1 <- CreateTableOne(vars = vars, strata = "isq_factor", factorVars = catvars, data = sob1)
#print(tabla_1)
#kableone(tabla_1)
#variables categorizadas
#Las variables van directamente sin el nombre del dataframe
catvars2 = c( "sexo","tbq", "iam", "hta", "dial", "loca_factor","topo_factor", "infe_factor", "ede_factor", "neur_factor", "prof_factor", "area_factor", "zon_factor", "fase_factor")
vars2 = c("edad","dias_ev","san_elian", "sexo","tbq", "iam", "hta", "dial", "loca_factor","topo_factor", "infe_factor", "ede_factor", "neur_factor", "prof_factor", "area_factor", "zon_factor", "fase_factor", "am_cont", "ant_pie" )
tabla_2 <- CreateTableOne(
vars = vars2,
strata = "isq_factor",
factorVars = catvars2,
data = sob1,
addOverall = TRUE
)
print(tabla_2)## Stratified by isq_factor
## Overall NO/leve SI
## n 406 237 169
## edad (mean (SD)) 57.41 (11.09) 53.34 (9.97) 63.12 (10.03)
## dias_ev (mean (SD)) 50.26 (112.85) 47.52 (118.71) 54.08 (104.37)
## san_elian (mean (SD)) 18.26 (4.09) 16.88 (3.99) 20.19 (3.41)
## sexo = M (%) 324 (79.8) 191 (80.6) 133 (78.7)
## tbq (%)
## 0 169 (41.8) 107 (45.3) 62 (36.9)
## 1 80 (19.8) 41 (17.4) 39 (23.2)
## 2 155 (38.4) 88 (37.3) 67 (39.9)
## iam = 1 (%) 83 (20.5) 32 (13.6) 51 (30.4)
## hta = 1 (%) 273 (67.4) 144 (60.8) 129 (76.8)
## dial = 1 (%) 19 ( 4.7) 10 ( 4.2) 9 ( 5.3)
## loca_factor = Met/Tar (%) 210 (51.7) 135 (57.0) 75 (44.4)
## topo_factor = lat/Mas de dos (%) 254 (62.6) 133 (56.1) 121 (71.6)
## infe_factor (%)
## No 86 (21.2) 55 (23.2) 31 (18.3)
## Leve/Mod 255 (62.8) 146 (61.6) 109 (64.5)
## Grave 65 (16.0) 36 (15.2) 29 (17.2)
## ede_factor (%)
## Sin/local 241 (59.4) 135 (57.0) 106 (62.7)
## Unilateral 124 (30.5) 78 (32.9) 46 (27.2)
## Bilateral 41 (10.1) 24 (10.1) 17 (10.1)
## neur_factor = Mod/charcot (%) 377 (92.9) 230 (97.0) 147 (87.0)
## prof_factor = parcial/total (%) 374 (92.1) 211 (89.0) 163 (96.4)
## area_factor (%)
## Menor a 10 188 (46.3) 121 (51.1) 67 (39.6)
## 10 a 40 157 (38.7) 86 (36.3) 71 (42.0)
## Mas de 40 61 (15.0) 30 (12.7) 31 (18.3)
## zon_factor = Dos o mas (%) 164 (40.4) 83 (35.0) 81 (47.9)
## fase_factor = Granu/infla (%) 387 (95.3) 221 (93.2) 166 (98.2)
## am_cont = 1 (%) 23 ( 5.7) 5 ( 2.1) 18 (10.7)
## ant_pie (%)
## 0 142 (35.0) 69 (29.1) 73 (43.2)
## 1 263 (64.8) 167 (70.5) 96 (56.8)
## 2 1 ( 0.2) 1 ( 0.4) 0 ( 0.0)
## Stratified by isq_factor
## p test
## n
## edad (mean (SD)) <0.001
## dias_ev (mean (SD)) 0.565
## san_elian (mean (SD)) <0.001
## sexo = M (%) 0.732
## tbq (%) 0.171
## 0
## 1
## 2
## iam = 1 (%) <0.001
## hta = 1 (%) 0.001
## dial = 1 (%) 0.778
## loca_factor = Met/Tar (%) 0.016
## topo_factor = lat/Mas de dos (%) 0.002
## infe_factor (%) 0.479
## No
## Leve/Mod
## Grave
## ede_factor (%) 0.450
## Sin/local
## Unilateral
## Bilateral
## neur_factor = Mod/charcot (%) <0.001
## prof_factor = parcial/total (%) 0.011
## area_factor (%) 0.057
## Menor a 10
## 10 a 40
## Mas de 40
## zon_factor = Dos o mas (%) 0.012
## fase_factor = Granu/infla (%) 0.036
## am_cont = 1 (%) 0.001
## ant_pie (%) 0.010
## 0
## 1
## 2print(tabla_2,
nonnormal = c("dias_ev", "san_elian"),
showAllLevels = TRUE,
quote = FALSE,
noSpaces = TRUE)## Stratified by isq_factor
## level Overall
## n 406
## edad (mean (SD)) 57.41 (11.09)
## dias_ev (median [IQR]) 20.00 [10.00, 45.00]
## san_elian (median [IQR]) 18.00 [15.00, 21.00]
## sexo (%) F 82 (20.2)
## M 324 (79.8)
## tbq (%) 0 169 (41.8)
## 1 80 (19.8)
## 2 155 (38.4)
## iam (%) 0 321 (79.5)
## 1 83 (20.5)
## hta (%) 0 132 (32.6)
## 1 273 (67.4)
## dial (%) 0 387 (95.3)
## 1 19 (4.7)
## loca_factor (%) Fal 196 (48.3)
## Met/Tar 210 (51.7)
## topo_factor (%) dor/plan 152 (37.4)
## lat/Mas de dos 254 (62.6)
## infe_factor (%) No 86 (21.2)
## Leve/Mod 255 (62.8)
## Grave 65 (16.0)
## ede_factor (%) Sin/local 241 (59.4)
## Unilateral 124 (30.5)
## Bilateral 41 (10.1)
## neur_factor (%) No/leve 29 (7.1)
## Mod/charcot 377 (92.9)
## prof_factor (%) Piel 32 (7.9)
## parcial/total 374 (92.1)
## area_factor (%) Menor a 10 188 (46.3)
## 10 a 40 157 (38.7)
## Mas de 40 61 (15.0)
## zon_factor (%) Una 242 (59.6)
## Dos o mas 164 (40.4)
## fase_factor (%) Epit 19 (4.7)
## Granu/infla 387 (95.3)
## am_cont (%) 0 383 (94.3)
## 1 23 (5.7)
## ant_pie (%) 0 142 (35.0)
## 1 263 (64.8)
## 2 1 (0.2)
## Stratified by isq_factor
## NO/leve SI p
## n 237 169
## edad (mean (SD)) 53.34 (9.97) 63.12 (10.03) <0.001
## dias_ev (median [IQR]) 20.00 [7.50, 30.00] 25.00 [14.00, 60.00] 0.020
## san_elian (median [IQR]) 17.00 [14.00, 20.00] 20.00 [18.00, 23.00] <0.001
## sexo (%) 46 (19.4) 36 (21.3) 0.732
## 191 (80.6) 133 (78.7)
## tbq (%) 107 (45.3) 62 (36.9) 0.171
## 41 (17.4) 39 (23.2)
## 88 (37.3) 67 (39.9)
## iam (%) 204 (86.4) 117 (69.6) <0.001
## 32 (13.6) 51 (30.4)
## hta (%) 93 (39.2) 39 (23.2) 0.001
## 144 (60.8) 129 (76.8)
## dial (%) 227 (95.8) 160 (94.7) 0.778
## 10 (4.2) 9 (5.3)
## loca_factor (%) 102 (43.0) 94 (55.6) 0.016
## 135 (57.0) 75 (44.4)
## topo_factor (%) 104 (43.9) 48 (28.4) 0.002
## 133 (56.1) 121 (71.6)
## infe_factor (%) 55 (23.2) 31 (18.3) 0.479
## 146 (61.6) 109 (64.5)
## 36 (15.2) 29 (17.2)
## ede_factor (%) 135 (57.0) 106 (62.7) 0.450
## 78 (32.9) 46 (27.2)
## 24 (10.1) 17 (10.1)
## neur_factor (%) 7 (3.0) 22 (13.0) <0.001
## 230 (97.0) 147 (87.0)
## prof_factor (%) 26 (11.0) 6 (3.6) 0.011
## 211 (89.0) 163 (96.4)
## area_factor (%) 121 (51.1) 67 (39.6) 0.057
## 86 (36.3) 71 (42.0)
## 30 (12.7) 31 (18.3)
## zon_factor (%) 154 (65.0) 88 (52.1) 0.012
## 83 (35.0) 81 (47.9)
## fase_factor (%) 16 (6.8) 3 (1.8) 0.036
## 221 (93.2) 166 (98.2)
## am_cont (%) 232 (97.9) 151 (89.3) 0.001
## 5 (2.1) 18 (10.7)
## ant_pie (%) 69 (29.1) 73 (43.2) 0.010
## 167 (70.5) 96 (56.8)
## 1 (0.4) 0 (0.0)
## Stratified by isq_factor
## test
## n
## edad (mean (SD))
## dias_ev (median [IQR]) nonnorm
## san_elian (median [IQR]) nonnorm
## sexo (%)
##
## tbq (%)
##
##
## iam (%)
##
## hta (%)
##
## dial (%)
##
## loca_factor (%)
##
## topo_factor (%)
##
## infe_factor (%)
##
##
## ede_factor (%)
##
##
## neur_factor (%)
##
## prof_factor (%)
##
## area_factor (%)
##
##
## zon_factor (%)
##
## fase_factor (%)
##
## am_cont (%)
##
## ant_pie (%)
##
## kableone(tabla_2)| Overall | NO/leve | SI | p | test | |
|---|---|---|---|---|---|
| n | 406 | 237 | 169 | ||
| edad (mean (SD)) | 57.41 (11.09) | 53.34 (9.97) | 63.12 (10.03) | <0.001 | |
| dias_ev (mean (SD)) | 50.26 (112.85) | 47.52 (118.71) | 54.08 (104.37) | 0.565 | |
| san_elian (mean (SD)) | 18.26 (4.09) | 16.88 (3.99) | 20.19 (3.41) | <0.001 | |
| sexo = M (%) | 324 (79.8) | 191 (80.6) | 133 (78.7) | 0.732 | |
| tbq (%) | 0.171 | ||||
| 0 | 169 (41.8) | 107 (45.3) | 62 (36.9) | ||
| 1 | 80 (19.8) | 41 (17.4) | 39 (23.2) | ||
| 2 | 155 (38.4) | 88 (37.3) | 67 (39.9) | ||
| iam = 1 (%) | 83 (20.5) | 32 (13.6) | 51 (30.4) | <0.001 | |
| hta = 1 (%) | 273 (67.4) | 144 (60.8) | 129 (76.8) | 0.001 | |
| dial = 1 (%) | 19 ( 4.7) | 10 ( 4.2) | 9 ( 5.3) | 0.778 | |
| loca_factor = Met/Tar (%) | 210 (51.7) | 135 (57.0) | 75 (44.4) | 0.016 | |
| topo_factor = lat/Mas de dos (%) | 254 (62.6) | 133 (56.1) | 121 (71.6) | 0.002 | |
| infe_factor (%) | 0.479 | ||||
| No | 86 (21.2) | 55 (23.2) | 31 (18.3) | ||
| Leve/Mod | 255 (62.8) | 146 (61.6) | 109 (64.5) | ||
| Grave | 65 (16.0) | 36 (15.2) | 29 (17.2) | ||
| ede_factor (%) | 0.450 | ||||
| Sin/local | 241 (59.4) | 135 (57.0) | 106 (62.7) | ||
| Unilateral | 124 (30.5) | 78 (32.9) | 46 (27.2) | ||
| Bilateral | 41 (10.1) | 24 (10.1) | 17 (10.1) | ||
| neur_factor = Mod/charcot (%) | 377 (92.9) | 230 (97.0) | 147 (87.0) | <0.001 | |
| prof_factor = parcial/total (%) | 374 (92.1) | 211 (89.0) | 163 (96.4) | 0.011 | |
| area_factor (%) | 0.057 | ||||
| Menor a 10 | 188 (46.3) | 121 (51.1) | 67 (39.6) | ||
| 10 a 40 | 157 (38.7) | 86 (36.3) | 71 (42.0) | ||
| Mas de 40 | 61 (15.0) | 30 (12.7) | 31 (18.3) | ||
| zon_factor = Dos o mas (%) | 164 (40.4) | 83 (35.0) | 81 (47.9) | 0.012 | |
| fase_factor = Granu/infla (%) | 387 (95.3) | 221 (93.2) | 166 (98.2) | 0.036 | |
| am_cont = 1 (%) | 23 ( 5.7) | 5 ( 2.1) | 18 (10.7) | 0.001 | |
| ant_pie (%) | 0.010 | ||||
| 0 | 142 (35.0) | 69 (29.1) | 73 (43.2) | ||
| 1 | 263 (64.8) | 167 (70.5) | 96 (56.8) | ||
| 2 | 1 ( 0.2) | 1 ( 0.4) | 0 ( 0.0) |
#elemento surv tiempo a la cicatrización
tiempo_evento <- Surv(sob1$dias_out,sob1$cicat)
#evaluacion de la cantidad de eventos por tiempo de seguimiento a intervalos de a 30 dias
sob1 %>%
filter(cicat == 1) %>%
mutate(intervalo = cut(dias_out,
breaks = seq(0, max(dias_out), by = 30))) %>%
count(intervalo) %>%
ggplot(aes(x = intervalo, y = n)) +
geom_col(fill = "darkorange") +
labs(x = "Intervalos de tiempo (30 dias)",
y = "Eventos",
title = "Eventos por intervalo de seguimiento") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Kaplan meier
# Graficar las curvas de Kaplan-Meier estratificadas por isq
km_global <- survfit(tiempo_evento ~ 1, data = sob1)#ponemos 1 para marcar que es crudo, sin estratificar por ninguna variable
#6b. Ajustar el modelo Kaplan-Meier estratificado por "beck_cat"
km_fit <- survfit(tiempo_evento ~ isq_factor,data = sob1)
ggsurvplot(
km_fit,
data = sob1,
fun = "event", # 👈 ESTO hace la curva ascendente
risk.table = TRUE,
pval = TRUE,
conf.int = TRUE,
cumcensor = TRUE,
ggtheme = theme_minimal(),
palette = "simpsons",
xlab = "Tiempo de seguimiento (dias)",
ylab = "Probabilidad acumulada del evento"
)## Ignoring unknown labels:
## • linetype : "1"
## Ignoring unknown labels:
## • linetype : "1"
## Ignoring unknown labels:
## • linetype : "1"
## Ignoring unknown labels:
## • linetype : "1"
## Ignoring unknown labels:
## • colour : "Strata"
## Ignoring unknown labels:
## • colour : "Strata"
#logrank
logrank <- survdiff(tiempo_evento ~ isq_factor, data = sob1)
logrank## Call:
## survdiff(formula = tiempo_evento ~ isq_factor, data = sob1)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## isq_factor=NO/leve 237 162 120 14.4 32.3
## isq_factor=SI 169 61 103 16.9 32.3
##
## Chisq= 32.3 on 1 degrees of freedom, p= 0.00000001summary(km_fit)$table## records n.max n.start events rmean se(rmean) median 0.95LCL
## isq_factor=NO/leve 237 237 237 162 120 6.55 91 84
## isq_factor=SI 169 169 169 61 195 11.94 170 153
## 0.95UCL
## isq_factor=NO/leve 109
## isq_factor=SI 275#mediana
km_fit$table[, c("median","0.95LCL","0.95UCL")]## NULLsurv_median(km_fit)## strata median lower upper
## 1 isq_factor=NO/leve 91 84 109
## 2 isq_factor=SI 170 153 275#modelo de cox. Se realiza para evaluar los supuestos
#Efecto de la isquemia sobre el tiempo a la cicatrización sin eventos competitivos ajustado por edad y edema
cox_model <- coxph(tiempo_evento ~ isq_factor+ edad + ede_factor , data = sob1)
summary(cox_model)## Call:
## coxph(formula = tiempo_evento ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71978 0.48686 0.16765 -4.29 0.000018 ***
## edad -0.01610 0.98402 0.00765 -2.11 0.035 *
## ede_factorUnilateral -0.20647 0.81345 0.15253 -1.35 0.176
## ede_factorBilateral -0.15043 0.86034 0.23500 -0.64 0.522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.676
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.860 1.16 0.543 1.364
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007tab_model(cox_model,
show.ci = TRUE, # Mostrar intervalo de confianza
show.se = TRUE, # Mostrar error estándar
transform = "exp", # Para mostrar HR en lugar de log(HR)
dv.labels = "Modelo de Cox - Datos drugs")| Modelo de Cox - Datos drugs | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| isq factorSI | 0.49 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.035 |
| ede factor [Unilateral] | 0.81 | 0.12 | 0.00 – Inf | 0.176 |
| ede factorBilateral | 0.86 | 0.20 | 0.00 – Inf | 0.522 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.095 | |||
confint(cox_model)## 2.5 % 97.5 %
## isq_factorSI -1.0484 -0.39118
## edad -0.0311 -0.00111
## ede_factorUnilateral -0.5054 0.09249
## ede_factorBilateral -0.6110 0.31016#supuesto de proporcionalidad
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph <- cox.zph(cox_model)
test_ph## chisq df p
## isq_factor 1.8077 1 0.18
## edad 0.0327 1 0.86
## ede_factor 1.4240 2 0.49
## GLOBAL 3.6145 4 0.46#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
#gráficos loglog
##loglog isq_factor
ggsurvplot(
survfit(tiempo_evento ~ isq_factor, data = sob1),
fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons"
)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.## Ignoring unknown labels:
## • fill : "Strata"
## • linetype : "1"
## Ignoring unknown labels:
## • fill : "Strata"
## • linetype : "1"
#ggsurvplot(
# survfit(tiempo_evento ~ infe_factor, data = sob1),
# fun = "cloglog",
# xlab = "Tiempo (dias)",
# ylab = "Log(-Log(S(t)))",
# palette = "simpsons"
#)
ggsurvplot(
survfit(tiempo_evento ~ ede_factor, data = sob1),
fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons"
)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.## Ignoring unknown labels:
## • fill : "Strata"
## • linetype : "1"
## Ignoring unknown labels:
## • fill : "Strata"
## • linetype : "1"
#Residuos de Martingale y deviance para la adecuación del modelo
# Martingale y Deviance
ggcoxdiagnostics(cox_model, type = "martingale", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggcoxdiagnostics(cox_model, type = "deviance", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
#Discriminación y calibración del COX
#C de harrell
summary(cox_model)$concordance## C se(C)
## 0.6301 0.0213#bootstrap
#Indice D de Somer (Dxy), Metricas Q D U, slope y R2 con paquete `rms`
dd <- datadist(sob1); options(datadist = "dd")## Warning in datadist(sob1): ...44 is constantcox_metrics <- cph(tiempo_evento ~ edad + isq_factor + ede_factor,
data = sob1, x = TRUE, y = TRUE, surv = TRUE)
# Validamos métricas con bootstrap
invisible(capture.output({ #esta linea es para evitar que se imprima el output del bootstrap)
cox_rms <- validate(cox_metrics, method = "boot", B = 200)
}))
# Revisamos las metricas validadas (en la columna index.corrected)
cox_rms## index.orig training test optimism index.corrected n
## Dxy 0.2601 0.2658 0.2460 0.0198 0.2403 200
## R2 0.0954 0.1049 0.0887 0.0162 0.0792 200
## Slope 1.0000 1.0000 0.9291 0.0709 0.9291 200
## D 0.0177 0.0199 0.0164 0.0035 0.0143 200
## U -0.0009 -0.0009 0.0006 -0.0015 0.0006 200
## Q 0.0186 0.0208 0.0158 0.0050 0.0136 200
## g 0.5290 0.5566 0.5027 0.0539 0.4751 200round(cox_rms[,5],2)## Dxy R2 Slope D U Q g
## 0.24 0.08 0.93 0.01 0.00 0.01 0.48# Calibración. Evaluación Gráfica
invisible(capture.output({cal <- calibrate(cox_metrics, method = "boot", B = 200, u = 120)
}))
plot(cal,
subtitles = FALSE,
lwd = 2,
lty = 2,
main = "Grafico de Calibracion",
cex.main = 1.2,
cex.axis = 0.7,
cex.lab = 0.9,
xlab = "Probabilidad Predicha de Cicatrizacion",
ylab = "Cicatrizacion Observada")
La curva azul está sistemáticamente por encima de la diagonal.
Eso sugiere que el modelo está sobreestimando la probabilidad de cicatrización.
Porque ignora que:
muchos pacientes mueren
muchos se amputan
Y por lo tanto, inflan la probabilidad teórica de cicatrizar.
#evento cica tomando evento competitivo am y muerte
sob1 <- sob1 %>%
mutate(
evento_cica = case_when(
estado_unif == "cerro_herida" ~ 1,
estado_unif %in% c( "fallecio", "amputacion") ~ 2,
estado_unif %in% c("perdido_seguimiento", "persiste") ~ 0,
TRUE ~ NA_real_
)
)
table(sob1$evento_cica)##
## 0 1 2
## 112 223 71#modelo de causa específica
#Opcion con objeto creado
#si evento es 1, es porque se amputo o persiste
#si evento es 2, es porque se murió o se perdió
tiempo_cicat <- Surv(sob1$dias_out, sob1$evento_cica==1)
tiempo_compet <- Surv(sob1$dias_out, sob1$evento_cica==2)
cox_cicat <- coxph(tiempo_cicat ~
isq_factor + edad+
ede_factor , data = sob1)
summary(cox_cicat)## Call:
## coxph(formula = tiempo_cicat ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71978 0.48686 0.16765 -4.29 0.000018 ***
## edad -0.01610 0.98402 0.00765 -2.11 0.035 *
## ede_factorUnilateral -0.20647 0.81345 0.15253 -1.35 0.176
## ede_factorBilateral -0.15043 0.86034 0.23500 -0.64 0.522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.676
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.860 1.16 0.543 1.364
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007cox_compet <- coxph(tiempo_compet ~
isq_factor + edad+
ede_factor , data = sob1)
summary(cox_compet)## Call:
## coxph(formula = tiempo_compet ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 406, number of events= 71
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.2845 3.6128 0.3009 4.27 0.00002 ***
## edad 0.0268 1.0272 0.0124 2.16 0.030 *
## ede_factorUnilateral 0.5849 1.7948 0.2557 2.29 0.022 *
## ede_factorBilateral 0.2923 1.3395 0.4145 0.71 0.481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 3.61 0.277 2.003 6.52
## edad 1.03 0.974 1.003 1.05
## ede_factorUnilateral 1.79 0.557 1.087 2.96
## ede_factorBilateral 1.34 0.747 0.594 3.02
##
## Concordance= 0.71 (se = 0.031 )
## Likelihood ratio test= 43 on 4 df, p=0.00000001
## Wald test = 36.9 on 4 df, p=0.0000002
## Score (logrank) test = 42.4 on 4 df, p=0.00000001#CIF
cifivhx <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor, data=sob1)
cifivhx## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 86 0.499 0.034 0.430, 0.564
## NO/leve 200 22 0.742 0.033 0.670, 0.801
## NO/leve 300 3 0.890 0.029 0.816, 0.936
## SI 100 71 0.183 0.032 0.126, 0.249
## SI 200 24 0.394 0.044 0.308, 0.479
## SI 300 8 0.496 0.051 0.393, 0.590## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 86 0.075 0.018 0.045, 0.114
## NO/leve 200 22 0.075 0.018 0.045, 0.114
## NO/leve 300 3 0.075 0.018 0.045, 0.114
## SI 100 71 0.282 0.036 0.215, 0.353
## SI 200 24 0.337 0.039 0.263, 0.413
## SI 300 8 0.350 0.040 0.273, 0.428## • Tests## outcome statistic df p.value
## 1 56.8 1.00 <0.001
## 2 40.4 1.00 <0.001#Gráfico CIF
library(tidycmprsk)
library(ggsurvfit)
ggcuminc(
cifivhx,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO/leve" = "#56B4E9", "SI" = "#E69F00")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c(
"Evento de interes (cicatrizacion)",
"Evento competitivo (muerte / amputacion mayor)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_step()`).
#Modelo de causa específica Cuando hacés predicción con este modelo, lo que obtenés es:
“Probabilidad de cicatrizar si los eventos competitivos fueran censura”
Es decir:
Trata muerte y amputación como si fueran pérdidas no informativas.
La probabilidad predicha no es la CIF real.
Entonces tu gráfico 1 está calibrando:
Riesgo derivado del modelo causa-específica (no la incidencia acumulada real).
#modelo de causa específica
modelo_CSC <- CSC(
formula = Hist(dias_out, evento_cica) ~ isq_factor+edad + ede_factor, data = sob1,
cause = 1 # Aclara que CHD es el evento de interés
)
# Ver resultados. Ojo! No funciona si hacés summary(modelo)
# Vamos a tener el modelo de interés y debajo el modelo competitivo
modelo_CSC## CSC(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor, data = sob1, cause = 1)
##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## ----------> Cause: 1
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71978 0.48686 0.16765 -4.29 0.000018 ***
## edad -0.01610 0.98402 0.00765 -2.11 0.035 *
## ede_factorUnilateral -0.20647 0.81345 0.15253 -1.35 0.176
## ede_factorBilateral -0.15043 0.86034 0.23500 -0.64 0.522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.676
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.860 1.16 0.543 1.364
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007
##
##
##
## ----------> Cause: 2
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 406, number of events= 71
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.2845 3.6128 0.3009 4.27 0.00002 ***
## edad 0.0268 1.0272 0.0124 2.16 0.030 *
## ede_factorUnilateral 0.5849 1.7948 0.2557 2.29 0.022 *
## ede_factorBilateral 0.2923 1.3395 0.4145 0.71 0.481
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 3.61 0.277 2.003 6.52
## edad 1.03 0.974 1.003 1.05
## ede_factorUnilateral 1.79 0.557 1.087 2.96
## ede_factorBilateral 1.34 0.747 0.594 3.02
##
## Concordance= 0.71 (se = 0.031 )
## Likelihood ratio test= 43 on 4 df, p=0.00000001
## Wald test = 36.9 on 4 df, p=0.0000002
## Score (logrank) test = 42.4 on 4 df, p=0.00000001#Modelo de Fine and Gray Probabilidad real acumulada de cicatrización considerando muerte y amputación como eventos competitivos.
Tu segundo gráfico está calibrando:
Probabilidad real acumulada observada vs predicha (CIF).
#modelo de Fine and gray
modelo_fg <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9719 0.378 0.16781 -5.791 0.000000007
## edad -0.0172 0.983 0.00771 -2.236 0.025000000
## ede_factorUnilateral -0.3437 0.709 0.14962 -2.297 0.022000000
## ede_factorBilateral -0.1976 0.821 0.23140 -0.854 0.390000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.378 2.64 0.272 0.526
## edad 0.983 1.02 0.968 0.998
## ede_factorUnilateral 0.709 1.41 0.529 0.951
## ede_factorBilateral 0.821 1.22 0.521 1.292
##
## Num. cases = 406
## Pseudo Log-likelihood = -1151
## Pseudo likelihood ratio test = 70.8 on 4 df,
##
## Convergence: TRUE#Gray test
library(cmprsk)
cuminc_obj <- cmprsk::cuminc(
ftime = sob1$dias_out,
fstatus = sob1$evento_cica,
group = sob1$isq_factor
)
cuminc_obj## Tests:
## stat pv df
## 1 56.8 0.0000000000000484 1
## 2 40.4 0.0000000002049831 1
## Estimates and Variances:
## $est
## 100 200 300
## NO/leve 1 0.4989 0.7424 0.8903
## SI 1 0.1830 0.3941 0.4955
## NO/leve 2 0.0751 0.0751 0.0751
## SI 2 0.2822 0.3372 0.3500
##
## $var
## 100 200 300
## NO/leve 1 0.001167 0.001115 0.000864
## SI 1 0.000998 0.001915 0.002577
## NO/leve 2 0.000309 0.000309 0.000309
## SI 2 0.001267 0.001483 0.001592¿Poner medianas de cicatrización?
En análisis con riesgos competitivos:
👉 La mediana basada en Kaplan–Meier NO corresponde 👉 La mediana basada en CIF tampoco es tan informativa clínicamente
Porque:
La CIF no es una función de supervivencia clásica
La mediana puede no ser alcanzable
Clínicamente es más interpretable decir:
“La probabilidad acumulada de cicatrización a 180 días fue 62% en pacientes sin isquemia vs 38% en pacientes con isquemia.”
Eso es muchísimo más potente.
*CIF (figura principal)
*Probabilidad absoluta a 180 días
*sHR (crudo y ajustado)
*Sensibilidad con cause-specific (opcional, no central)
# Evaluar performance predictiva con Score. AUC
score_fg <- Score(
object = list("FG model" = modelo_fg),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1,
times = c(17:200),
metrics = "auc",
summary = "IPA",
cens.model = "km",
cause = 1
)
score_fg$AUC$score## model times AUC se lower upper
## <fctr> <int> <num> <num> <num> <num>
## 1: FG model 17 0.503 0.1179 0.272 0.734
## 2: FG model 18 0.530 0.1064 0.322 0.739
## 3: FG model 19 0.530 0.1064 0.322 0.739
## 4: FG model 20 0.530 0.1064 0.322 0.739
## 5: FG model 21 0.574 0.0900 0.397 0.750
## ---
## 180: FG model 196 0.723 0.0313 0.662 0.785
## 181: FG model 197 0.720 0.0316 0.658 0.782
## 182: FG model 198 0.720 0.0316 0.658 0.782
## 183: FG model 199 0.725 0.0317 0.663 0.787
## 184: FG model 200 0.724 0.0320 0.661 0.787#comparacion causa especifica y FyG
#score comparado con ipa
score_comp <-Score(
list("FG model" = modelo_fg, "CSC model" = modelo_CSC),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1,
times = c(17:200),
cause = 1,
metrics ="auc",
summary = "IPA"
)#Extraigo las IPAs por tiempo y modelo
ipa <- as.data.frame(score_comp$Brier$score)[, c("model", "times", "IPA")]
ipa## model times IPA
## 1 Null model 17 0.00000
## 2 Null model 18 0.00000
## 3 Null model 19 0.00000
## 4 Null model 20 0.00000
## 5 Null model 21 0.00000
## 6 Null model 22 0.00000
## 7 Null model 23 0.00000
## 8 Null model 24 0.00000
## 9 Null model 25 0.00000
## 10 Null model 26 0.00000
## 11 Null model 27 0.00000
## 12 Null model 28 0.00000
## 13 Null model 29 0.00000
## 14 Null model 30 0.00000
## 15 Null model 31 0.00000
## 16 Null model 32 0.00000
## 17 Null model 33 0.00000
## 18 Null model 34 0.00000
## 19 Null model 35 0.00000
## 20 Null model 36 0.00000
## 21 Null model 37 0.00000
## 22 Null model 38 0.00000
## 23 Null model 39 0.00000
## 24 Null model 40 0.00000
## 25 Null model 41 0.00000
## 26 Null model 42 0.00000
## 27 Null model 43 0.00000
## 28 Null model 44 0.00000
## 29 Null model 45 0.00000
## 30 Null model 46 0.00000
## 31 Null model 47 0.00000
## 32 Null model 48 0.00000
## 33 Null model 49 0.00000
## 34 Null model 50 0.00000
## 35 Null model 51 0.00000
## 36 Null model 52 0.00000
## 37 Null model 53 0.00000
## 38 Null model 54 0.00000
## 39 Null model 55 0.00000
## 40 Null model 56 0.00000
## 41 Null model 57 0.00000
## 42 Null model 58 0.00000
## 43 Null model 59 0.00000
## 44 Null model 60 0.00000
## 45 Null model 61 0.00000
## 46 Null model 62 0.00000
## 47 Null model 63 0.00000
## 48 Null model 64 0.00000
## 49 Null model 65 0.00000
## 50 Null model 66 0.00000
## 51 Null model 67 0.00000
## 52 Null model 68 0.00000
## 53 Null model 69 0.00000
## 54 Null model 70 0.00000
## 55 Null model 71 0.00000
## 56 Null model 72 0.00000
## 57 Null model 73 0.00000
## 58 Null model 74 0.00000
## 59 Null model 75 0.00000
## 60 Null model 76 0.00000
## 61 Null model 77 0.00000
## 62 Null model 78 0.00000
## 63 Null model 79 0.00000
## 64 Null model 80 0.00000
## 65 Null model 81 0.00000
## 66 Null model 82 0.00000
## 67 Null model 83 0.00000
## 68 Null model 84 0.00000
## 69 Null model 85 0.00000
## 70 Null model 86 0.00000
## 71 Null model 87 0.00000
## 72 Null model 88 0.00000
## 73 Null model 89 0.00000
## 74 Null model 90 0.00000
## 75 Null model 91 0.00000
## 76 Null model 92 0.00000
## 77 Null model 93 0.00000
## 78 Null model 94 0.00000
## 79 Null model 95 0.00000
## 80 Null model 96 0.00000
## 81 Null model 97 0.00000
## 82 Null model 98 0.00000
## 83 Null model 99 0.00000
## 84 Null model 100 0.00000
## 85 Null model 101 0.00000
## 86 Null model 102 0.00000
## 87 Null model 103 0.00000
## 88 Null model 104 0.00000
## 89 Null model 105 0.00000
## 90 Null model 106 0.00000
## 91 Null model 107 0.00000
## 92 Null model 108 0.00000
## 93 Null model 109 0.00000
## 94 Null model 110 0.00000
## 95 Null model 111 0.00000
## 96 Null model 112 0.00000
## 97 Null model 113 0.00000
## 98 Null model 114 0.00000
## 99 Null model 115 0.00000
## 100 Null model 116 0.00000
## 101 Null model 117 0.00000
## 102 Null model 118 0.00000
## 103 Null model 119 0.00000
## 104 Null model 120 0.00000
## 105 Null model 121 0.00000
## 106 Null model 122 0.00000
## 107 Null model 123 0.00000
## 108 Null model 124 0.00000
## 109 Null model 125 0.00000
## 110 Null model 126 0.00000
## 111 Null model 127 0.00000
## 112 Null model 128 0.00000
## 113 Null model 129 0.00000
## 114 Null model 130 0.00000
## 115 Null model 131 0.00000
## 116 Null model 132 0.00000
## 117 Null model 133 0.00000
## 118 Null model 134 0.00000
## 119 Null model 135 0.00000
## 120 Null model 136 0.00000
## 121 Null model 137 0.00000
## 122 Null model 138 0.00000
## 123 Null model 139 0.00000
## 124 Null model 140 0.00000
## 125 Null model 141 0.00000
## 126 Null model 142 0.00000
## 127 Null model 143 0.00000
## 128 Null model 144 0.00000
## 129 Null model 145 0.00000
## 130 Null model 146 0.00000
## 131 Null model 147 0.00000
## 132 Null model 148 0.00000
## 133 Null model 149 0.00000
## 134 Null model 150 0.00000
## 135 Null model 151 0.00000
## 136 Null model 152 0.00000
## 137 Null model 153 0.00000
## 138 Null model 154 0.00000
## 139 Null model 155 0.00000
## 140 Null model 156 0.00000
## 141 Null model 157 0.00000
## 142 Null model 158 0.00000
## 143 Null model 159 0.00000
## 144 Null model 160 0.00000
## 145 Null model 161 0.00000
## 146 Null model 162 0.00000
## 147 Null model 163 0.00000
## 148 Null model 164 0.00000
## 149 Null model 165 0.00000
## 150 Null model 166 0.00000
## 151 Null model 167 0.00000
## 152 Null model 168 0.00000
## 153 Null model 169 0.00000
## 154 Null model 170 0.00000
## 155 Null model 171 0.00000
## 156 Null model 172 0.00000
## 157 Null model 173 0.00000
## 158 Null model 174 0.00000
## 159 Null model 175 0.00000
## 160 Null model 176 0.00000
## 161 Null model 177 0.00000
## 162 Null model 178 0.00000
## 163 Null model 179 0.00000
## 164 Null model 180 0.00000
## 165 Null model 181 0.00000
## 166 Null model 182 0.00000
## 167 Null model 183 0.00000
## 168 Null model 184 0.00000
## 169 Null model 185 0.00000
## 170 Null model 186 0.00000
## 171 Null model 187 0.00000
## 172 Null model 188 0.00000
## 173 Null model 189 0.00000
## 174 Null model 190 0.00000
## 175 Null model 191 0.00000
## 176 Null model 192 0.00000
## 177 Null model 193 0.00000
## 178 Null model 194 0.00000
## 179 Null model 195 0.00000
## 180 Null model 196 0.00000
## 181 Null model 197 0.00000
## 182 Null model 198 0.00000
## 183 Null model 199 0.00000
## 184 Null model 200 0.00000
## 185 FG model 17 -0.00473
## 186 FG model 18 -0.00342
## 187 FG model 19 -0.00342
## 188 FG model 20 -0.00342
## 189 FG model 21 -0.00040
## 190 FG model 22 -0.00420
## 191 FG model 23 -0.00420
## 192 FG model 24 0.00449
## 193 FG model 25 0.00474
## 194 FG model 26 0.00474
## 195 FG model 27 0.00724
## 196 FG model 28 0.01362
## 197 FG model 29 0.01362
## 198 FG model 30 0.01362
## 199 FG model 31 0.01362
## 200 FG model 32 0.01656
## 201 FG model 33 0.02070
## 202 FG model 34 0.02073
## 203 FG model 35 0.03658
## 204 FG model 36 0.04297
## 205 FG model 37 0.04297
## 206 FG model 38 0.04297
## 207 FG model 39 0.04425
## 208 FG model 40 0.04251
## 209 FG model 41 0.04935
## 210 FG model 42 0.06500
## 211 FG model 43 0.07229
## 212 FG model 44 0.06823
## 213 FG model 45 0.06960
## 214 FG model 46 0.06494
## 215 FG model 47 0.07900
## 216 FG model 48 0.07900
## 217 FG model 49 0.07780
## 218 FG model 50 0.07328
## 219 FG model 51 0.07501
## 220 FG model 52 0.07510
## 221 FG model 53 0.07834
## 222 FG model 54 0.07681
## 223 FG model 55 0.07681
## 224 FG model 56 0.07665
## 225 FG model 57 0.07409
## 226 FG model 58 0.07123
## 227 FG model 59 0.06740
## 228 FG model 60 0.06638
## 229 FG model 61 0.06757
## 230 FG model 62 0.06757
## 231 FG model 63 0.06691
## 232 FG model 64 0.06952
## 233 FG model 65 0.06929
## 234 FG model 66 0.07067
## 235 FG model 67 0.07094
## 236 FG model 68 0.06580
## 237 FG model 69 0.06944
## 238 FG model 70 0.07471
## 239 FG model 71 0.08167
## 240 FG model 72 0.08915
## 241 FG model 73 0.08380
## 242 FG model 74 0.08380
## 243 FG model 75 0.08629
## 244 FG model 76 0.08104
## 245 FG model 77 0.08309
## 246 FG model 78 0.09034
## 247 FG model 79 0.09096
## 248 FG model 80 0.09409
## 249 FG model 81 0.09409
## 250 FG model 82 0.09531
## 251 FG model 83 0.09133
## 252 FG model 84 0.09397
## 253 FG model 85 0.09929
## 254 FG model 86 0.09544
## 255 FG model 87 0.09458
## 256 FG model 88 0.09843
## 257 FG model 89 0.10471
## 258 FG model 90 0.10348
## 259 FG model 91 0.12012
## 260 FG model 92 0.11877
## 261 FG model 93 0.11877
## 262 FG model 94 0.11877
## 263 FG model 95 0.12059
## 264 FG model 96 0.12606
## 265 FG model 97 0.10903
## 266 FG model 98 0.11382
## 267 FG model 99 0.11454
## 268 FG model 100 0.11454
## 269 FG model 101 0.11321
## 270 FG model 102 0.11680
## 271 FG model 103 0.11717
## 272 FG model 104 0.11175
## 273 FG model 105 0.11905
## 274 FG model 106 0.12098
## 275 FG model 107 0.12411
## 276 FG model 108 0.12061
## 277 FG model 109 0.11851
## 278 FG model 110 0.11949
## 279 FG model 111 0.12408
## 280 FG model 112 0.12127
## 281 FG model 113 0.12346
## 282 FG model 114 0.12500
## 283 FG model 115 0.12166
## 284 FG model 116 0.12166
## 285 FG model 117 0.11914
## 286 FG model 118 0.11914
## 287 FG model 119 0.13104
## 288 FG model 120 0.13055
## 289 FG model 121 0.12967
## 290 FG model 122 0.13407
## 291 FG model 123 0.13407
## 292 FG model 124 0.13407
## 293 FG model 125 0.13407
## 294 FG model 126 0.13912
## 295 FG model 127 0.14527
## 296 FG model 128 0.14474
## 297 FG model 129 0.14336
## 298 FG model 130 0.14519
## 299 FG model 131 0.14519
## 300 FG model 132 0.14519
## 301 FG model 133 0.14474
## 302 FG model 134 0.14474
## 303 FG model 135 0.14665
## 304 FG model 136 0.14665
## 305 FG model 137 0.14665
## 306 FG model 138 0.14249
## 307 FG model 139 0.14818
## 308 FG model 140 0.15151
## 309 FG model 141 0.15553
## 310 FG model 142 0.16455
## 311 FG model 143 0.16591
## 312 FG model 144 0.16354
## 313 FG model 145 0.16120
## 314 FG model 146 0.16120
## 315 FG model 147 0.16107
## 316 FG model 148 0.15475
## 317 FG model 149 0.15475
## 318 FG model 150 0.16012
## 319 FG model 151 0.15789
## 320 FG model 152 0.15789
## 321 FG model 153 0.14825
## 322 FG model 154 0.14076
## 323 FG model 155 0.14076
## 324 FG model 156 0.14076
## 325 FG model 157 0.13842
## 326 FG model 158 0.13527
## 327 FG model 159 0.13527
## 328 FG model 160 0.13022
## 329 FG model 161 0.13022
## 330 FG model 162 0.13815
## 331 FG model 163 0.13131
## 332 FG model 164 0.13131
## 333 FG model 165 0.13131
## 334 FG model 166 0.13178
## 335 FG model 167 0.13677
## 336 FG model 168 0.13380
## 337 FG model 169 0.13380
## 338 FG model 170 0.12691
## 339 FG model 171 0.12691
## 340 FG model 172 0.12361
## 341 FG model 173 0.12361
## 342 FG model 174 0.12375
## 343 FG model 175 0.12375
## 344 FG model 176 0.12375
## 345 FG model 177 0.12938
## 346 FG model 178 0.12938
## 347 FG model 179 0.12938
## 348 FG model 180 0.12938
## 349 FG model 181 0.14238
## 350 FG model 182 0.14925
## 351 FG model 183 0.14925
## 352 FG model 184 0.14925
## 353 FG model 185 0.14925
## 354 FG model 186 0.16729
## 355 FG model 187 0.16729
## 356 FG model 188 0.16729
## 357 FG model 189 0.16729
## 358 FG model 190 0.16767
## 359 FG model 191 0.15728
## 360 FG model 192 0.16272
## 361 FG model 193 0.16272
## 362 FG model 194 0.15336
## 363 FG model 195 0.15336
## 364 FG model 196 0.14951
## 365 FG model 197 0.14516
## 366 FG model 198 0.14516
## 367 FG model 199 0.14996
## 368 FG model 200 0.14799
## 369 CSC model 17 -0.00356
## 370 CSC model 18 -0.00231
## 371 CSC model 19 -0.00231
## 372 CSC model 20 -0.00231
## 373 CSC model 21 0.00051
## 374 CSC model 22 -0.00269
## 375 CSC model 23 -0.00269
## 376 CSC model 24 0.00469
## 377 CSC model 25 0.00494
## 378 CSC model 26 0.00494
## 379 CSC model 27 0.00778
## 380 CSC model 28 0.01304
## 381 CSC model 29 0.01304
## 382 CSC model 30 0.01304
## 383 CSC model 31 0.01304
## 384 CSC model 32 0.01541
## 385 CSC model 33 0.01976
## 386 CSC model 34 0.02071
## 387 CSC model 35 0.03474
## 388 CSC model 36 0.04022
## 389 CSC model 37 0.04022
## 390 CSC model 38 0.04022
## 391 CSC model 39 0.04123
## 392 CSC model 40 0.03973
## 393 CSC model 41 0.04586
## 394 CSC model 42 0.06016
## 395 CSC model 43 0.06691
## 396 CSC model 44 0.06295
## 397 CSC model 45 0.06413
## 398 CSC model 46 0.06018
## 399 CSC model 47 0.07416
## 400 CSC model 48 0.07416
## 401 CSC model 49 0.07289
## 402 CSC model 50 0.06858
## 403 CSC model 51 0.07003
## 404 CSC model 52 0.07055
## 405 CSC model 53 0.07371
## 406 CSC model 54 0.07208
## 407 CSC model 55 0.07208
## 408 CSC model 56 0.07178
## 409 CSC model 57 0.06994
## 410 CSC model 58 0.06756
## 411 CSC model 59 0.06406
## 412 CSC model 60 0.06334
## 413 CSC model 61 0.06513
## 414 CSC model 62 0.06513
## 415 CSC model 63 0.06473
## 416 CSC model 64 0.06684
## 417 CSC model 65 0.06645
## 418 CSC model 66 0.06762
## 419 CSC model 67 0.06841
## 420 CSC model 68 0.06337
## 421 CSC model 69 0.06649
## 422 CSC model 70 0.07172
## 423 CSC model 71 0.07776
## 424 CSC model 72 0.08477
## 425 CSC model 73 0.07980
## 426 CSC model 74 0.07980
## 427 CSC model 75 0.08228
## 428 CSC model 76 0.07705
## 429 CSC model 77 0.08012
## 430 CSC model 78 0.08769
## 431 CSC model 79 0.08879
## 432 CSC model 80 0.09162
## 433 CSC model 81 0.09162
## 434 CSC model 82 0.09239
## 435 CSC model 83 0.08864
## 436 CSC model 84 0.09199
## 437 CSC model 85 0.09700
## 438 CSC model 86 0.09343
## 439 CSC model 87 0.09284
## 440 CSC model 88 0.09625
## 441 CSC model 89 0.10196
## 442 CSC model 90 0.10070
## 443 CSC model 91 0.11648
## 444 CSC model 92 0.11513
## 445 CSC model 93 0.11513
## 446 CSC model 94 0.11513
## 447 CSC model 95 0.11682
## 448 CSC model 96 0.12232
## 449 CSC model 97 0.10598
## 450 CSC model 98 0.11137
## 451 CSC model 99 0.11169
## 452 CSC model 100 0.11169
## 453 CSC model 101 0.11036
## 454 CSC model 102 0.11358
## 455 CSC model 103 0.11360
## 456 CSC model 104 0.10839
## 457 CSC model 105 0.11496
## 458 CSC model 106 0.11675
## 459 CSC model 107 0.11927
## 460 CSC model 108 0.11674
## 461 CSC model 109 0.11415
## 462 CSC model 110 0.11560
## 463 CSC model 111 0.12093
## 464 CSC model 112 0.11754
## 465 CSC model 113 0.11959
## 466 CSC model 114 0.12072
## 467 CSC model 115 0.11744
## 468 CSC model 116 0.11744
## 469 CSC model 117 0.11639
## 470 CSC model 118 0.11639
## 471 CSC model 119 0.12790
## 472 CSC model 120 0.12711
## 473 CSC model 121 0.12674
## 474 CSC model 122 0.13119
## 475 CSC model 123 0.13119
## 476 CSC model 124 0.13119
## 477 CSC model 125 0.13119
## 478 CSC model 126 0.13596
## 479 CSC model 127 0.14210
## 480 CSC model 128 0.14209
## 481 CSC model 129 0.14146
## 482 CSC model 130 0.14330
## 483 CSC model 131 0.14330
## 484 CSC model 132 0.14330
## 485 CSC model 133 0.14305
## 486 CSC model 134 0.14305
## 487 CSC model 135 0.14552
## 488 CSC model 136 0.14552
## 489 CSC model 137 0.14552
## 490 CSC model 138 0.14182
## 491 CSC model 139 0.14727
## 492 CSC model 140 0.15046
## 493 CSC model 141 0.15436
## 494 CSC model 142 0.16376
## 495 CSC model 143 0.16558
## 496 CSC model 144 0.16319
## 497 CSC model 145 0.16072
## 498 CSC model 146 0.16072
## 499 CSC model 147 0.16056
## 500 CSC model 148 0.15449
## 501 CSC model 149 0.15449
## 502 CSC model 150 0.15959
## 503 CSC model 151 0.15741
## 504 CSC model 152 0.15741
## 505 CSC model 153 0.14825
## 506 CSC model 154 0.14090
## 507 CSC model 155 0.14090
## 508 CSC model 156 0.14090
## 509 CSC model 157 0.13858
## 510 CSC model 158 0.13631
## 511 CSC model 159 0.13631
## 512 CSC model 160 0.13040
## 513 CSC model 161 0.13040
## 514 CSC model 162 0.13816
## 515 CSC model 163 0.13086
## 516 CSC model 164 0.13086
## 517 CSC model 165 0.13086
## 518 CSC model 166 0.13164
## 519 CSC model 167 0.13683
## 520 CSC model 168 0.13373
## 521 CSC model 169 0.13373
## 522 CSC model 170 0.12685
## 523 CSC model 171 0.12685
## 524 CSC model 172 0.12338
## 525 CSC model 173 0.12338
## 526 CSC model 174 0.12372
## 527 CSC model 175 0.12372
## 528 CSC model 176 0.12372
## 529 CSC model 177 0.13027
## 530 CSC model 178 0.13027
## 531 CSC model 179 0.13027
## 532 CSC model 180 0.13027
## 533 CSC model 181 0.14351
## 534 CSC model 182 0.15017
## 535 CSC model 183 0.15017
## 536 CSC model 184 0.15017
## 537 CSC model 185 0.15017
## 538 CSC model 186 0.16840
## 539 CSC model 187 0.16840
## 540 CSC model 188 0.16840
## 541 CSC model 189 0.16840
## 542 CSC model 190 0.16954
## 543 CSC model 191 0.15779
## 544 CSC model 192 0.16312
## 545 CSC model 193 0.16312
## 546 CSC model 194 0.15315
## 547 CSC model 195 0.15315
## 548 CSC model 196 0.14899
## 549 CSC model 197 0.14428
## 550 CSC model 198 0.14428
## 551 CSC model 199 0.15025
## 552 CSC model 200 0.14889# Gráfico
library(ggplot2)
ggplot(ipa, aes(x = times, y = IPA, color = model)) +
geom_line(linewidth = 1.2) +
geom_point(size = 3) +
geom_text(aes(label = percent(IPA, accuracy = 0.1)),
vjust = -0.8, size = 3, show.legend = FALSE) +
scale_y_continuous(labels = percent_format(accuracy = 0.1),
limits = c(0, NA)) +
scale_color_manual(
values = c(
"Null model" = "gray40",
"FG model" = "#2ECC71",
"CSC model" = "#3498DB"
),
drop = FALSE # <- muy importante para que no elimine niveles
) +
labs(
title = "(IPA) en el tiempo",
x = "Dias de seguimiento", y = "IPA", color = "Modelo"
) +
theme_minimal() +
theme(legend.position = "bottom", plot.title.position = "plot")+ scale_x_continuous(
limits = c(30, 200),
breaks = seq(30, 200, by = 10)
)## Warning: Removed 39 rows containing missing values or values outside the scale range
## (`geom_line()`).## Warning: Removed 39 rows containing missing values or values outside the scale range
## (`geom_point()`).## Warning: Removed 39 rows containing missing values or values outside the scale range
## (`geom_text()`).
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" )
#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 cicatrización a 180 días",
ylim = c(0, max(pred_no, pred_si)),
col = "black")
lines(edades, pred_si,
lwd = 3,
lty = 2,
col = "black")
legend("bottomleft",
legend = c("Sin isquemia/isq leve", "Con isquemia moderada/grave"),
lwd = 3,
lty = c(1,2),
bty = "n")
#ANALISIS DE SENSIBILIDAD
#excluyendo del análisis a los pacientes con herida persistente seguidos 30 días o menos
sob1_menos30 <- sob1 %>%
filter(!(estado_unif == "persiste" & dias_out < 31))
tiempo_evento_s <- Surv(sob1_menos30$dias_out,sob1_menos30$cicat)
#modelo de cox
#Efecto de la isquemia sobre el tiempo a la cicatrización sin eventos competitivos ajustado por edad y edema
cox_model_s <- coxph(tiempo_evento_s ~ isq_factor+ edad + ede_factor , data = sob1_menos30)
summary(cox_model_s)## Call:
## coxph(formula = tiempo_evento_s ~ isq_factor + edad + ede_factor,
## data = sob1_menos30)
##
## n= 390, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71920 0.48714 0.16763 -4.29 0.000018 ***
## edad -0.01613 0.98400 0.00764 -2.11 0.035 *
## ede_factorUnilateral -0.20657 0.81337 0.15254 -1.35 0.176
## ede_factorBilateral -0.15199 0.85899 0.23499 -0.65 0.518
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.677
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.859 1.16 0.542 1.362
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007tab_model(cox_model_s,
show.ci = TRUE, # Mostrar intervalo de confianza
show.se = TRUE, # Mostrar error estándar
transform = "exp", # Para mostrar HR en lugar de log(HR)
dv.labels = "Modelo de Cox - Datos drugs")| Modelo de Cox - Datos drugs | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| isq factorSI | 0.49 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.035 |
| ede factor [Unilateral] | 0.81 | 0.12 | 0.00 – Inf | 0.176 |
| ede factorBilateral | 0.86 | 0.20 | 0.00 – Inf | 0.518 |
| Observations | 390 | |||
| R2 Nagelkerke | 0.099 | |||
confint(cox_model_s)## 2.5 % 97.5 %
## isq_factorSI -1.0478 -0.39065
## edad -0.0311 -0.00115
## ede_factorUnilateral -0.5055 0.09240
## ede_factorBilateral -0.6126 0.30859#Test de proporcionalidad de riesgos de Schoenfeld
test_ph_s <- cox.zph(cox_model_s)
test_ph_s## chisq df p
## isq_factor 1.7970 1 0.18
## edad 0.0335 1 0.85
## ede_factor 1.4034 2 0.50
## GLOBAL 3.5787 4 0.47#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph_s)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
#modelo de Fine and gray
modelo_fg_s <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor
, data = sob1_menos30,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg_s##
## Right-censored response of a competing.risks model
##
## No.Observations: 390
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 96
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor, data = sob1_menos30, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9733 0.378 0.16798 -5.79 0.0000000069
## edad -0.0173 0.983 0.00771 -2.24 0.0250000000
## ede_factorUnilateral -0.3477 0.706 0.14997 -2.32 0.0200000000
## ede_factorBilateral -0.1990 0.820 0.23122 -0.86 0.3900000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.378 2.65 0.272 0.525
## edad 0.983 1.02 0.968 0.998
## ede_factorUnilateral 0.706 1.42 0.526 0.948
## ede_factorBilateral 0.820 1.22 0.521 1.289
##
## Num. cases = 390
## Pseudo Log-likelihood = -1151
## Pseudo likelihood ratio test = 71.2 on 4 df,
##
## Convergence: TRUE# Calibracion global
ss <- Score(list("FGR" = modelo_fg_s),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1_menos30,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(ss, models = "FGR" )
#Considerando a los pacientes con herida persistente seguidos 30 días o menoscomo TODOS CICATRIZADOS sob1_sens
sob1_sens <- sob1
sob1_sens <- sob1 %>%
mutate(
estado_unif = ifelse(
estado_unif == "persiste" & dias_out < 31,
"cerro_herida",
estado_unif
)
)
table(sob1_sens$estado_unif)##
## amputacion cerro_herida fallecio persiste
## 57 239 14 96tiempo_evento_s1 <- Surv(sob1_sens$dias_out,sob1_sens$cicat)
#modelo de cox
#Efecto de la isquemia sobre el tiempo a la cicatrización sin eventos competitivos ajustado por edad y edema
cox_model_s1 <- coxph(tiempo_evento_s1 ~ isq_factor+ edad + ede_factor , data = sob1_sens)
summary(cox_model_s1)## Call:
## coxph(formula = tiempo_evento_s1 ~ isq_factor + edad + ede_factor,
## data = sob1_sens)
##
## n= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -0.71978 0.48686 0.16765 -4.29 0.000018 ***
## edad -0.01610 0.98402 0.00765 -2.11 0.035 *
## ede_factorUnilateral -0.20647 0.81345 0.15253 -1.35 0.176
## ede_factorBilateral -0.15043 0.86034 0.23500 -0.64 0.522
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.487 2.05 0.351 0.676
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.813 1.23 0.603 1.097
## ede_factorBilateral 0.860 1.16 0.543 1.364
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.5 on 4 df, p=0.00000003
## Wald test = 37.1 on 4 df, p=0.0000002
## Score (logrank) test = 38.9 on 4 df, p=0.00000007tab_model(cox_model_s1,
show.ci = TRUE, # Mostrar intervalo de confianza
show.se = TRUE, # Mostrar error estándar
transform = "exp", # Para mostrar HR en lugar de log(HR)
dv.labels = "Modelo de Cox - Datos drugs")| Modelo de Cox - Datos drugs | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| isq factorSI | 0.49 | 0.08 | 0.00 – Inf | <0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.035 |
| ede factor [Unilateral] | 0.81 | 0.12 | 0.00 – Inf | 0.176 |
| ede factorBilateral | 0.86 | 0.20 | 0.00 – Inf | 0.522 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.095 | |||
confint(cox_model_s1)## 2.5 % 97.5 %
## isq_factorSI -1.0484 -0.39118
## edad -0.0311 -0.00111
## ede_factorUnilateral -0.5054 0.09249
## ede_factorBilateral -0.6110 0.31016#Test de proporcionalidad de riesgos de Schoenfeld
test_ph_s1 <- cox.zph(cox_model_s1)
test_ph_s1## chisq df p
## isq_factor 1.8077 1 0.18
## edad 0.0327 1 0.86
## ede_factor 1.4240 2 0.49
## GLOBAL 3.6145 4 0.46#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph_s1)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
#evento cica: agrega como cicatrizados a los perdidos
sob1_sens <- sob1_sens %>%
mutate(
evento_cica = ifelse(
sob1$estado_unif == "persiste" & dias_out < 31,
1,
evento_cica
)
)
#modelo de Fine and gray
modelo_fg_s1 <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor
, data = sob1_sens,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg_s1##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 239 0
## 2 71 0
## unknown 0 96
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor, data = sob1_sens, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9534 0.385 0.16094 -5.924 0.0000000031
## edad -0.0168 0.983 0.00718 -2.339 0.0190000000
## ede_factorUnilateral -0.3064 0.736 0.14465 -2.118 0.0340000000
## ede_factorBilateral -0.2083 0.812 0.22283 -0.935 0.3500000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.385 2.59 0.281 0.528
## edad 0.983 1.02 0.970 0.997
## ede_factorUnilateral 0.736 1.36 0.554 0.977
## ede_factorBilateral 0.812 1.23 0.525 1.257
##
## Num. cases = 406
## Pseudo Log-likelihood = -1247
## Pseudo likelihood ratio test = 72.2 on 4 df,
##
## Convergence: TRUE# Calibracion global
s1 <- Score(list("FGR" = modelo_fg_s1),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1_sens,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(s1, models = "FGR" )
#Análisis exploratorio sobre el efecto de la isquemia sobre la AM tomando como competitivos muerte y cicatrizacion
#evento_ampu
sob1 <- sob1 %>%
mutate(
evento_ampu = case_when(
estado_unif %in% c("cerro_herida","fallecio") ~ 2,
estado_unif %in% c("amputacion") ~ 1,
estado_unif %in% c("perdido_seguimiento", "persiste") ~ 0,
TRUE ~ NA_real_
)
)
table(sob1$evento_ampu)##
## 0 1 2
## 112 57 237tiempo_evento_am <- Surv(sob1$dias_out,sob1$ampu)
#modelo de cox
#Efecto de la isquemia sobre el tiempo a la amputacion
cox_model_a <- coxph(tiempo_evento_am ~ isq_factor+ edad + ede_factor , data = sob1)
summary(cox_model_a)## Call:
## coxph(formula = tiempo_evento_am ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 406, number of events= 57
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 1.5974 4.9404 0.3586 4.45 0.0000084 ***
## edad 0.0232 1.0235 0.0138 1.68 0.093 .
## ede_factorUnilateral 0.6092 1.8389 0.2835 2.15 0.032 *
## ede_factorBilateral 0.1648 1.1792 0.4856 0.34 0.734
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 4.94 0.202 2.446 9.98
## edad 1.02 0.977 0.996 1.05
## ede_factorUnilateral 1.84 0.544 1.055 3.21
## ede_factorBilateral 1.18 0.848 0.455 3.05
##
## Concordance= 0.737 (se = 0.029 )
## Likelihood ratio test= 42.6 on 4 df, p=0.00000001
## Wald test = 33.9 on 4 df, p=0.0000008
## Score (logrank) test = 41.6 on 4 df, p=0.00000002tab_model(cox_model_a,
show.ci = TRUE, # Mostrar intervalo de confianza
show.se = TRUE, # Mostrar error estándar
transform = "exp", # Para mostrar HR en lugar de log(HR)
dv.labels = "Modelo de Cox - Datos drugs")| Modelo de Cox - Datos drugs | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| isq factorSI | 4.94 | 1.77 | 0.00 – Inf | <0.001 |
| edad | 1.02 | 0.01 | 0.00 – Inf | 0.093 |
| ede factor [Unilateral] | 1.84 | 0.52 | 0.00 – Inf | 0.032 |
| ede factorBilateral | 1.18 | 0.57 | 0.00 – Inf | 0.734 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.125 | |||
confint(cox_model_a)## 2.5 % 97.5 %
## isq_factorSI 0.89454 2.3003
## edad -0.00385 0.0503
## ede_factorUnilateral 0.05352 1.1648
## ede_factorBilateral -0.78693 1.1165#Test de proporcionalidad de riesgos de Schoenfeld
test_ph_a <- cox.zph(cox_model_a)
test_ph_a## chisq df p
## isq_factor 1.568891 1 0.21
## edad 0.000702 1 0.98
## ede_factor 3.858376 2 0.15
## GLOBAL 5.486547 4 0.24#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph_a)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
#modelo de Fine and gray
modelo_fg_a <- FGR(
formula = Hist(dias_out, evento_ampu) ~ isq_factor + edad+ede_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg_a##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 57 0
## 2 237 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_ampu) ~ isq_factor + edad +
## ede_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI 1.7075 5.52 0.3684 4.635 0.0000036
## edad 0.0213 1.02 0.0115 1.850 0.0640000
## ede_factorUnilateral 0.6072 1.84 0.2865 2.120 0.0340000
## ede_factorBilateral 0.1664 1.18 0.4689 0.355 0.7200000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 5.52 0.181 2.679 11.35
## edad 1.02 0.979 0.999 1.04
## ede_factorUnilateral 1.84 0.545 1.047 3.22
## ede_factorBilateral 1.18 0.847 0.471 2.96
##
## Num. cases = 406
## Pseudo Log-likelihood = -310
## Pseudo likelihood ratio test = 47.3 on 4 df,
##
## Convergence: TRUE# Calibracion global
sa <- Score(list("FGR" = modelo_fg_a),
formula = Hist(dias_out, evento_ampu) ~ 1,
data = sob1,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(sa, models = "FGR" )
table(sob1$evento_ampu)##
## 0 1 2
## 112 57 237#cif
#CIF
cifivhxa <- cuminc(Surv(dias_out, factor(evento_ampu)) ~ isq_factor, data=sob1)
cifivhxa## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 86 0.049 0.014 0.026, 0.082
## NO/leve 200 22 0.049 0.014 0.026, 0.082
## NO/leve 300 3 0.049 0.014 0.026, 0.082
## SI 100 71 0.244 0.034 0.181, 0.312
## SI 200 24 0.283 0.036 0.214, 0.356
## SI 300 8 0.296 0.038 0.224, 0.372## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 86 0.525 0.034 0.456, 0.590
## NO/leve 200 22 0.769 0.032 0.698, 0.825
## NO/leve 300 3 0.917 0.028 0.842, 0.957
## SI 100 71 0.221 0.034 0.159, 0.290
## SI 200 24 0.448 0.044 0.360, 0.532
## SI 300 8 0.549 0.051 0.445, 0.642## • Tests## outcome statistic df p.value
## 1 40.4 1.00 <0.001
## 2 50.8 1.00 <0.001ggcuminc(
cifivhxa,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO/leve" = "#56B4E9", "SI" = "#E69F00")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c(
"Evento de interes (amputacion)",
"Evento competitivo (muerte / cicatrizacion)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_step()`).
#altman
#altman
edades <- seq(30, 85, by = 1)
new_no_a <- data.frame(
isq_factor = factor(rep(levels(sob1$isq_factor)[1], length(edades)),
levels = levels(sob1$isq_factor)),
edad = edades,
ede_factor = factor(rep("Sin/local", length(edades)),
levels = levels(sob1$ede_factor))
)
new_si_a <- new_no_a
new_si_a$isq_factor <- factor(rep(levels(sob1$isq_factor)[2], length(edades)),
levels = levels(sob1$isq_factor))
pred_no_a <- predictRisk(modelo_fg_a,
newdata = new_no_a,
times = 180,
cause = 1)
pred_si_a <- predictRisk(modelo_fg_a,
newdata = new_si_a,
times = 180,
cause = 1)
plot(edades, pred_no_a,
type = "l",
lwd = 3,
xlab = "Edad (años)",
ylab = "Probabilidad de amputacion a 180 días",
ylim = c(0, max(pred_no_a, pred_si_a)),
col = "black")
lines(edades, pred_si_a,
lwd = 3,
lty = 2,
col = "black")
legend("bottomleft",
legend = c("Sin isquemia/isq leve", "Con isquemia moderada/grave"),
lwd = 3,
lty = c(1,2),
bty = "n")
#Gray test
cuminc_obja <- cmprsk::cuminc(
ftime = sob1$dias_out,
fstatus = sob1$evento_ampu,
group = sob1$isq_factor
)
cuminc_obja## Tests:
## stat pv df
## 1 40.4 0.00000000020845 1
## 2 50.8 0.00000000000105 1
## Estimates and Variances:
## $est
## 100 200 300
## NO/leve 1 0.0487 0.0487 0.0487
## SI 1 0.2440 0.2833 0.2961
## NO/leve 2 0.5253 0.7688 0.9167
## SI 2 0.2212 0.4480 0.5494
##
## $var
## 100 200 300
## NO/leve 1 0.000207 0.000207 0.000207
## SI 1 0.001141 0.001322 0.001441
## NO/leve 2 0.001165 0.001055 0.000769
## SI 2 0.001148 0.001970 0.002550infección como modificador del efecto de la isquemia sobre la cicatrización
#modelo de Fine and gray
table(sob1$infe_factor)##
## No Leve/Mod Grave
## 86 255 65table(sob1$evento_cica)##
## 0 1 2
## 112 223 71#no se puede hacer fine and gray. voy a hacer un modelo esratificado para la variable infección
sob1_grave <- subset(sob1, infe_factor == "Grave")
sob1_levemod <- subset(sob1, infe_factor == "Leve/Mod")
sob1_sin <- subset(sob1, infe_factor == "No")
#sin infeccion
modelo_fg_sin <- FGR( formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor , data = sob1_sin, cause = 1) # Aclara que Cicatrizacion es el evento de interés
modelo_fg_sin##
## Right-censored response of a competing.risks model
##
## No.Observations: 86
##
## Pattern:
##
## Cause event right.censored
## 1 48 0
## 2 11 0
## unknown 0 27
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor, data = sob1_sin, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -1.2354 0.291 0.4967 -2.487 0.013
## edad -0.0178 0.982 0.0155 -1.149 0.250
## ede_factorUnilateral 0.0648 1.067 0.4263 0.152 0.880
## ede_factorBilateral -0.3101 0.733 0.3725 -0.832 0.410
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.291 3.440 0.110 0.77
## edad 0.982 1.018 0.953 1.01
## ede_factorUnilateral 1.067 0.937 0.463 2.46
## ede_factorBilateral 0.733 1.364 0.353 1.52
##
## Num. cases = 86
## Pseudo Log-likelihood = -172
## Pseudo likelihood ratio test = 19.4 on 4 df,
##
## Convergence: TRUE#infeccion leve/mod
modelo_fg_levemod <- FGR( formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor , data = sob1_levemod, cause = 1) # Aclara que Cicatrizacion es el evento de interés
modelo_fg_levemod##
## Right-censored response of a competing.risks model
##
## No.Observations: 255
##
## Pattern:
##
## Cause event right.censored
## 1 147 0
## 2 36 0
## unknown 0 72
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor, data = sob1_levemod, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -0.9114 0.402 0.18345 -4.968 0.00000068
## edad -0.0145 0.986 0.00925 -1.570 0.12000000
## ede_factorUnilateral 0.0232 1.023 0.18241 0.127 0.90000000
## ede_factorBilateral -0.0382 0.962 0.28301 -0.135 0.89000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.402 2.488 0.281 0.576
## edad 0.986 1.015 0.968 1.004
## ede_factorUnilateral 1.023 0.977 0.716 1.463
## ede_factorBilateral 0.962 1.039 0.553 1.676
##
## Num. cases = 255
## Pseudo Log-likelihood = -689
## Pseudo likelihood ratio test = 38 on 4 df,
##
## Convergence: TRUE#infeccion grave
modelo_fg_grave <- FGR( formula = Hist(dias_out, evento_cica) ~ isq_factor + edad+ede_factor , data = sob1_grave, cause = 1) # Aclara que Cicatrizacion es el evento de interés
modelo_fg_grave##
## Right-censored response of a competing.risks model
##
## No.Observations: 65
##
## Pattern:
##
## Cause event right.censored
## 1 28 0
## 2 24 0
## unknown 0 13
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor + edad +
## ede_factor, data = sob1_grave, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI -1.9935 0.136 0.6891 -2.8929 0.00380
## edad -0.0115 0.989 0.0276 -0.4161 0.68000
## ede_factorUnilateral -1.4872 0.226 0.4127 -3.6038 0.00031
## ede_factorBilateral 0.0882 1.092 1.2836 0.0687 0.95000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 0.136 7.341 0.0353 0.526
## edad 0.989 1.012 0.9364 1.044
## ede_factorUnilateral 0.226 4.425 0.1007 0.507
## ede_factorBilateral 1.092 0.916 0.0882 13.517
##
## Num. cases = 65
## Pseudo Log-likelihood = -87.9
## Pseudo likelihood ratio test = 27.4 on 4 df,
##
## Convergence: TRUE#gráficos de cif para cada categoría de infección
#SIN INFECCIÓN
#cif
#CIF
cifivhxasin <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor, data=sob1_sin)
cifivhxasin## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 37 0.295 0.063 0.180, 0.420
## NO/leve 100 25 0.464 0.069 0.326, 0.592
## NO/leve 150 12 0.645 0.071 0.488, 0.765
## NO/leve 200 7 0.709 0.072 0.541, 0.824
## NO/leve 250 4 0.793 0.074 0.599, 0.901
## NO/leve 300 1 0.963 0.089 0.011, 1.00
## SI 50.0 19 0.097 0.054 0.024, 0.232
## SI 100 13 0.210 0.079 0.082, 0.377
## SI 150 10 0.210 0.079 0.082, 0.377
## SI 200 6 0.269 0.095 0.109, 0.460
## SI 250 3 0.269 0.095 0.109, 0.460
## SI 300 2 0.269 0.095 0.109, 0.460## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 37 0.019 0.019 0.001, 0.088
## NO/leve 100 25 0.037 0.026 0.007, 0.114
## NO/leve 150 12 0.037 0.026 0.007, 0.114
## NO/leve 200 7 0.037 0.026 0.007, 0.114
## NO/leve 250 4 0.037 0.026 0.007, 0.114
## NO/leve 300 1 0.037 0.026 0.007, 0.114
## SI 50.0 19 0.233 0.079 0.100, 0.398
## SI 100 13 0.268 0.083 0.124, 0.437
## SI 150 10 0.312 0.090 0.151, 0.488
## SI 200 6 0.312 0.090 0.151, 0.488
## SI 250 3 0.312 0.090 0.151, 0.488
## SI 300 2 0.312 0.090 0.151, 0.488## • Tests## outcome statistic df p.value
## 1 15.5 1.00 <0.001
## 2 12.1 1.00 <0.001ggcuminc(
cifivhxasin,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO/leve" = "#56B4E9", "SI" = "#E69F00")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c(
"Evento de interes (Cicatrizacion)",
"Evento competitivo (muerte / amputacion)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 4 rows containing missing values or values outside the scale range
## (`geom_step()`).
#INFECCIon leve/mod
#cif
#CIF
cifivhxasleve <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor, data=sob1_levemod)
cifivhxasleve##
## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 44 0.574 0.044 0.482, 0.655
## NO/leve 200 8 0.824 0.041 0.727, 0.890
## NO/leve 300 1 0.902 0.054 0.723, 0.968
## SI 100 50 0.210 0.041 0.136, 0.296
## SI 200 13 0.486 0.056 0.373, 0.590
## SI 300 5 0.579 0.060 0.452, 0.687## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 100 44 0.051 0.019 0.022, 0.096
## NO/leve 200 8 0.051 0.019 0.022, 0.096
## NO/leve 300 1 0.051 0.019 0.022, 0.096
## SI 100 50 0.231 0.042 0.154, 0.316
## SI 200 13 0.288 0.046 0.202, 0.380
## SI 300 5 0.288 0.046 0.202, 0.380## • Tests## outcome statistic df p.value
## 1 37.1 1.00 <0.001
## 2 22.2 1.00 <0.001ggcuminc(
cifivhxasleve,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO/leve" = "#56B4E9", "SI" = "#E69F00")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c(
"Evento de interes (amputacion)",
"Evento competitivo (muerte / cicatrizacion)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 10 rows containing missing values or values outside the scale range
## (`geom_step()`).
#INFECCIon grave
#cif
#CIF
cifivhxasgrave <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor, data=sob1_grave)
cifivhxasgrave##
## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 25 0.111 0.053 0.034, 0.239
## NO/leve 100 17 0.283 0.078 0.145, 0.438
## NO/leve 150 12 0.428 0.086 0.260, 0.586
## NO/leve 200 7 0.527 0.089 0.341, 0.683
## NO/leve 250 1 0.734 0.085 0.524, 0.862
## NO/leve 300 1 0.734 0.085 0.524, 0.862
## SI 50.0 16 0.000 0.000 NA, NA
## SI 100 8 0.047 0.048 0.003, 0.203
## SI 150 5 0.164 0.090 0.038, 0.369
## SI 200 5 0.164 0.090 0.038, 0.369
## SI 250 3 0.164 0.090 0.038, 0.369
## SI 300 1 0.320 0.126 0.106, 0.560## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 25 0.194 0.067 0.084, 0.338
## NO/leve 100 17 0.223 0.071 0.103, 0.372
## NO/leve 150 12 0.223 0.071 0.103, 0.372
## NO/leve 200 7 0.223 0.071 0.103, 0.372
## NO/leve 250 1 0.223 0.071 0.103, 0.372
## NO/leve 300 1 0.223 0.071 0.103, 0.372
## SI 50.0 16 0.448 0.095 0.261, 0.619
## SI 100 8 0.485 0.096 0.291, 0.654
## SI 150 5 0.544 0.104 0.324, 0.719
## SI 200 5 0.544 0.104 0.324, 0.719
## SI 250 3 0.602 0.109 0.360, 0.777
## SI 300 1 0.602 0.109 0.360, 0.777## • Tests## outcome statistic df p.value
## 1 9.28 1.00 0.002
## 2 8.13 1.00 0.004ggcuminc(
cifivhxasgrave,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO/leve" = "#56B4E9", "SI" = "#E69F00")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c(
"Evento de interes (cicatrizacion)",
"Evento competitivo (muerte / amputacion)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_step()`).
#Análisis secundario
table(sob1$isq_factor2)##
## NO LEVE MODERADA GRAVE
## 206 31 47 122#Efecto de la isquemia sobre el tiempo a la cicatrización sin eventos competitivos ajustado por edad y edema
cox_modelisq2 <- coxph(tiempo_evento ~ isq_factor2+ edad + ede_factor , data = sob1)
summary(cox_modelisq2)## Call:
## coxph(formula = tiempo_evento ~ isq_factor2 + edad + ede_factor,
## data = sob1)
##
## n= 406, number of events= 223
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factor2LEVE 0.0302 1.0306 0.2251 0.13 0.89339
## isq_factor2MODERADA -0.7480 0.4733 0.2394 -3.12 0.00178 **
## isq_factor2GRAVE -0.6937 0.4997 0.2008 -3.45 0.00055 ***
## edad -0.0163 0.9838 0.0077 -2.12 0.03421 *
## ede_factorUnilateral -0.2029 0.8163 0.1547 -1.31 0.18953
## ede_factorBilateral -0.1465 0.8638 0.2356 -0.62 0.53413
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factor2LEVE 1.031 0.97 0.663 1.602
## isq_factor2MODERADA 0.473 2.11 0.296 0.757
## isq_factor2GRAVE 0.500 2.00 0.337 0.741
## edad 0.984 1.02 0.969 0.999
## ede_factorUnilateral 0.816 1.22 0.603 1.105
## ede_factorBilateral 0.864 1.16 0.544 1.371
##
## Concordance= 0.63 (se = 0.021 )
## Likelihood ratio test= 40.6 on 6 df, p=0.0000003
## Wald test = 37.2 on 6 df, p=0.000002
## Score (logrank) test = 39 on 6 df, p=0.0000007tab_model(cox_modelisq2,
show.ci = TRUE, # Mostrar intervalo de confianza
show.se = TRUE, # Mostrar error estándar
transform = "exp", # Para mostrar HR en lugar de log(HR)
dv.labels = "Modelo de Cox - Datos drugs")| Modelo de Cox - Datos drugs | ||||
|---|---|---|---|---|
| Predictors | Estimates | std. Error | CI | p |
| isq factor2LEVE | 1.03 | 0.23 | 0.00 – Inf | 0.893 |
| isq factor2MODERADA | 0.47 | 0.11 | 0.00 – Inf | 0.002 |
| isq factor2GRAVE | 0.50 | 0.10 | 0.00 – Inf | 0.001 |
| edad | 0.98 | 0.01 | 0.00 – Inf | 0.034 |
| ede factor [Unilateral] | 0.82 | 0.13 | 0.00 – Inf | 0.190 |
| ede factorBilateral | 0.86 | 0.20 | 0.00 – Inf | 0.534 |
| Observations | 406 | |||
| R2 Nagelkerke | 0.096 | |||
confint(cox_modelisq2)## 2.5 % 97.5 %
## isq_factor2LEVE -0.4110 0.47135
## isq_factor2MODERADA -1.2172 -0.27882
## isq_factor2GRAVE -1.0872 -0.30013
## edad -0.0314 -0.00121
## ede_factorUnilateral -0.5061 0.10023
## ede_factorBilateral -0.6082 0.31525#
cifivhxas22 <- cuminc(Surv(dias_out, factor(evento_cica)) ~ isq_factor2, data=sob1)
cifivhxas22## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## NO 100 73 0.503 0.037 0.429, 0.572
## NO 200 20 0.726 0.036 0.647, 0.790
## NO 300 3 0.875 0.035 0.787, 0.928
## LEVE 100 13 0.478 0.096 0.284, 0.649
## LEVE 200 2 0.858 0.086 0.578, 0.958
## LEVE 300 0 0.968 0.064 0.190, 0.999
## MODERADA 100 22 0.240 0.068 0.122, 0.380
## MODERADA 200 12 0.463 0.083 0.297, 0.613
## MODERADA 300 3 0.636 0.106 0.394, 0.802
## GRAVE 100 49 0.162 0.035 0.100, 0.237
## GRAVE 200 12 0.376 0.053 0.273, 0.479
## GRAVE 300 5 0.442 0.058 0.327, 0.551## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## NO 100 73 0.082 0.020 0.049, 0.126
## NO 200 20 0.082 0.020 0.049, 0.126
## NO 300 3 0.082 0.020 0.049, 0.126
## LEVE 100 13 0.032 0.032 0.002, 0.144
## LEVE 200 2 0.032 0.032 0.002, 0.144
## LEVE 300 0 0.032 0.032 0.002, 0.144
## MODERADA 100 22 0.161 0.057 0.069, 0.287
## MODERADA 200 12 0.161 0.057 0.069, 0.287
## MODERADA 300 3 0.161 0.057 0.069, 0.287
## GRAVE 100 49 0.329 0.044 0.246, 0.415
## GRAVE 200 12 0.407 0.048 0.313, 0.498
## GRAVE 300 5 0.426 0.050 0.327, 0.522## • Tests## outcome statistic df p.value
## 1 57.7 3.00 <0.001
## 2 53.7 3.00 <0.001ggcuminc(
cifivhxas22,
outcome = c(1, 2),
linetype_aes = TRUE
) +
scale_color_manual(
name = "Isquemia",
values = c("NO" = "#56B4E9", "LEVE" = "#E69F00", "MODERADA"="#B58", "GRAVE"="#009E73")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c(
"Evento de interes (cicatrizacion)",
"Evento competitivo (muerte / amputacion)"
)
) +
labs(
x = "Tiempo (dias)",
y = "Incidencia acumulada"
) +
theme_minimal() +
theme(
legend.position = "bottom"
)## Warning: Removed 34 rows containing missing values or values outside the scale range
## (`geom_step()`).
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph_isq2 <- cox.zph(cox_modelisq2)
test_ph_isq2## chisq df p
## isq_factor2 3.8524 3 0.28
## edad 0.0327 1 0.86
## ede_factor 1.3943 2 0.50
## GLOBAL 5.6890 6 0.46#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph_isq2)## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.## Warning: ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
## ggtheme is not a valid theme.
## Please use `theme()` to construct themes.
#modelo de Fine and gray
modelo_fg_isq2 <- FGR(
formula = Hist(dias_out, evento_cica) ~ isq_factor2 + edad+ede_factor
, data = sob1,
cause = 1 # Aclara que Cicatrizacion es el evento de interés
)
modelo_fg_isq2##
## Right-censored response of a competing.risks model
##
## No.Observations: 406
##
## Pattern:
##
## Cause event right.censored
## 1 223 0
## 2 71 0
## unknown 0 112
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_cica) ~ isq_factor2 + edad +
## ede_factor, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factor2LEVE 0.0950 1.100 0.19799 0.480 0.630000000
## isq_factor2MODERADA -0.6968 0.498 0.22574 -3.087 0.002000000
## isq_factor2GRAVE -1.0946 0.335 0.20344 -5.380 0.000000074
## edad -0.0163 0.984 0.00782 -2.082 0.037000000
## ede_factorUnilateral -0.3405 0.711 0.15203 -2.240 0.025000000
## ede_factorBilateral -0.2122 0.809 0.23300 -0.911 0.360000000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factor2LEVE 1.100 0.909 0.746 1.621
## isq_factor2MODERADA 0.498 2.007 0.320 0.775
## isq_factor2GRAVE 0.335 2.988 0.225 0.499
## edad 0.984 1.016 0.969 0.999
## ede_factorUnilateral 0.711 1.406 0.528 0.958
## ede_factorBilateral 0.809 1.236 0.512 1.277
##
## Num. cases = 406
## Pseudo Log-likelihood = -1150
## Pseudo likelihood ratio test = 73.1 on 6 df,
##
## Convergence: TRUE# Calibracion global
s_isq2 <- Score(list("FGR" = modelo_fg_isq2),
formula = Hist(dias_out, evento_cica) ~ 1,
data = sob1,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(s_isq2, models = "FGR" )