library(readr)
library(data.table)
## Warning: package 'data.table' was built under R version 4.4.3
library(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 errors
library(broom)
library(vcd)
## Warning: package 'vcd' was built under R version 4.4.3
## Cargando paquete requerido: grid
library(sjPlot)
##
## Adjuntando el paquete: 'sjPlot'
##
## The following object is masked from 'package:ggplot2':
##
## set_theme
library(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 analyses
library(gtsummary)
## Warning: package 'gtsummary' was built under R version 4.4.3
library(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':
##
## summarise
library(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, units
library(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, var
library(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-10
library(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':
##
## lift
library(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_intercepts
options(scipen = 999, digits = 3, encoding = 'UTF-8')
library(naniar)
## Warning: package 'naniar' was built under R version 4.4.3
library(dplyr)
library(survival)
library(survminer)
##
## Adjuntando el paquete: 'survminer'
##
## The following object is masked from 'package:survival':
##
## myeloma
library(MASS)
##
## Adjuntando el paquete: 'MASS'
##
## The following object is masked from 'package:gtsummary':
##
## select
##
## The following object is masked from 'package:dplyr':
##
## select
library(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] TRUE
library(ggsurvfit)
## Warning: package 'ggsurvfit' was built under R version 4.4.3
library(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.20
library(rms)
library(cmprsk)
## Warning: package 'cmprsk' was built under R version 4.4.3
library(prodlim)
## Warning: package 'prodlim' was built under R version 4.4.3
library(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_factor
library(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':
##
## trial
library(readxl)
sob <- read_excel("C:/Users/Administrador/Downloads/sob.xlsx")
## New names:
## • `` -> `...32`
## • `` -> `...33`
View(sob)
sob1 <- sob %>%
filter(!if_any(7, is.na))
sob1 <- sob1 %>%
filter(!if_any(19:28, is.na))
sob1 <- sob1 %>%
filter(!if_any(5, is.na))
sob1 <- sob1 %>%
filter(!if_any(6, is.na))
table(sob1$FALLO)
##
## 0 1
## 135 112
sob1 <- sob1 %>%
mutate(across(13:29, as.factor))
sob1$cicat_factor <- factor(sob1$cicat, levels = c(0, 1), labels = c( "No Cicatrizo", "Cicatrizo"))
levels(sob1$cicat_factor)
## [1] "No Cicatrizo" "Cicatrizo"
summary(sob1$cicat_factor)
## No Cicatrizo Cicatrizo
## 112 135
names(sob1)<- tolower(names (sob1))
#ISQuEMIA
sob1 <- sob1 %>%
mutate(isq_factor = factor(
if_else(isq %in% c(2, 3), "SI", "NO/leve")
))
table(sob1$isq)
##
## 0 1 2 3
## 128 13 29 77
table(sob1$isq_factor)
##
## NO/leve SI
## 141 106
#localizacion
sob1 <- sob1 %>%
mutate(loca_factor = factor(if_else(local %in% c(2, 3), "Met/Tar", "Fal")),
loca_factor = factor(loca_factor, levels = rev(levels(loca_factor))))
table(sob1$loca)
## Warning: Unknown or uninitialised column: `loca`.
## < table of extent 0 >
table(sob1$loca_factor)
##
## Met/Tar Fal
## 136 111
#topog
sob1<- sob1 %>%
mutate(topo_factor = factor(if_else(topo %in% c(1,2), "dor/plan/lat", "Mas de dos")))
sob1 <- sob1 %>%
mutate(topo_factor = factor(if_else(topo %in% c(1, 2), "dor/plan/lat", "Mas de dos")),
topo_factor = factor(topo_factor, levels = rev(levels(topo_factor))))
table(sob1$topo)
##
## 1 2 3
## 98 48 101
table(sob1$topo_factor)
##
## Mas de dos dor/plan/lat
## 101 146
# zon
sob1 <- sob1 %>%
mutate(zon_factor = factor(if_else(zon %in% c(2, 3), "Dos o mas", "Una")))
table(sob1$zon)
##
## 0 1 2 3
## 1 144 51 51
table(sob1$zon_factor)
##
## Dos o mas Una
## 102 145
#infe
sob1 <- sob1 %>%
mutate(infe_factor = case_when(
inf == 0 ~ "No",
inf %in% c(1, 2) ~ "Leve/Mod",
inf == 3 ~ "Grave",
TRUE ~ NA_character_ # conserva los NA originales
)) %>%
mutate(infe_factor = factor(infe_factor, levels = c("No", "Leve/Mod", "Grave")))
sob1$infe_factor <- factor (sob1$infe_factor, levels= c ("Grave","Leve/Mod","No"))
table(sob1$infe)
## Warning: Unknown or uninitialised column: `infe`.
## < table of extent 0 >
table(sob1$infe_factor)
##
## Grave Leve/Mod No
## 51 143 52
#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
## 90 58 79 20
table(sob1$ede_factor)
##
## Bilateral Unilateral Sin/local
## 20 79 148
#neur
sob1 <- sob1 %>%
mutate(neur_factor = factor(
if_else(neu %in% c(1, 2, 3), "Leve/mod/charcot", "No"),
levels = c("No", "Leve/mod/charcot")
))
table(sob1$neur)
## Warning: Unknown or uninitialised column: `neur`.
## < table of extent 0 >
table(sob1$neur_factor)
##
## No Leve/mod/charcot
## 3 244
#tisu (prof)
sob1 <- sob1 %>%
mutate(tisu_factor = factor(if_else(prof %in% c(2, 3), "parcial/total", "Piel")))
table(sob1$tisu)
## Warning: Unknown or uninitialised column: `tisu`.
## < table of extent 0 >
table(sob1$tisu_factor)
##
## parcial/total Piel
## 231 16
#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
## 98 103 46
table(sob1$area_factor)
##
## Menor a 10 10 a 40 Mas de 40
## 98 103 46
#cica(fase)
sob1 <- sob1 %>%
mutate(cica_factor = factor(if_else(fase %in% c(2, 3), "Granu/infla", "Epit")))
sob1$cica_factor <- factor (sob1$cica_factor, levels= c ("Granu/infla","Epit"))
table(sob1$cica)
## Warning: Unknown or uninitialised column: `cica`.
## < table of extent 0 >
table(sob1$cica_factor)
##
## Granu/infla Epit
## 237 10
tiempo_evento <- Surv(sob1$dias_out,sob1$fallo)
cox_model <- coxph(tiempo_evento ~ isq_factor+ ede_factor , data = sob1)
summary(cox_model)
## Call:
## coxph(formula = tiempo_evento ~ isq_factor + ede_factor, data = sob1)
##
## n= 247, number of events= 112
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 0.5831 1.7917 0.1946 3.00 0.0027 **
## ede_factorUnilateral -0.0204 0.9799 0.3691 -0.06 0.9560
## ede_factorSin/local -0.4101 0.6636 0.3622 -1.13 0.2575
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 1.792 0.558 1.223 2.62
## ede_factorUnilateral 0.980 1.021 0.475 2.02
## ede_factorSin/local 0.664 1.507 0.326 1.35
##
## Concordance= 0.611 (se = 0.029 )
## Likelihood ratio test= 11.6 on 3 df, p=0.009
## Wald test = 11.6 on 3 df, p=0.009
## Score (logrank) test = 11.7 on 3 df, p=0.008
tab_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 factor [SI]
|
1.79
|
0.35
|
0.00 – Inf
|
0.003
|
|
ede factor [Unilateral]
|
0.98
|
0.36
|
0.00 – Inf
|
0.956
|
|
ede factor [Sin/local]
|
0.66
|
0.24
|
0.00 – Inf
|
0.258
|
|
Observations
|
247
|
|
R2 Nagelkerke
|
0.047
|
confint(cox_model)
## 2.5 % 97.5 %
## isq_factorSI 0.202 0.965
## ede_factorUnilateral -0.744 0.703
## ede_factorSin/local -1.120 0.300
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph <- cox.zph(cox_model)
test_ph
## chisq df p
## isq_factor 0.347 1 0.556
## ede_factor 5.520 2 0.063
## GLOBAL 6.176 3 0.103
#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.
#fallo_factor
sob1$fallo_factor <- factor(sob1$fallo, levels = c(0, 1), labels = c( "No fallo", "Fallo"))
#ver eventos
sob1 <- sob1 %>%
mutate(isq_factor1 = factor(if_else(prof %in% c("NO/leve"), 0, 1)))
# 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 (días)",
ylab = "Probabilidad acumulada del evento"
)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
## Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## 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 141 47 61.4 3.38 7.58
## isq_factor=SI 106 65 50.6 4.10 7.58
##
## Chisq= 7.6 on 1 degrees of freedom, p= 0.006
# 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'
cox_model_ede <- coxph(tiempo_evento ~ isq_factor + edad , data = sob1)
summary(cox_model_ede)
## Call:
## coxph(formula = tiempo_evento ~ isq_factor + edad, data = sob1)
##
## n= 246, number of events= 112
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 0.20838 1.23169 0.22173 0.94 0.3473
## edad 0.02945 1.02988 0.00999 2.95 0.0032 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 1.23 0.812 0.798 1.90
## edad 1.03 0.971 1.010 1.05
##
## Concordance= 0.628 (se = 0.031 )
## Likelihood ratio test= 16 on 2 df, p=0.0003
## Wald test = 16.5 on 2 df, p=0.0003
## Score (logrank) test = 16.5 on 2 df, p=0.0003
tab_model(cox_model_ede,
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 factor [SI]
|
1.23
|
0.27
|
0.00 – Inf
|
0.347
|
|
edad
|
1.03
|
0.01
|
0.00 – Inf
|
0.003
|
|
Observations
|
246
|
|
R2 Nagelkerke
|
0.064
|
confint(cox_model_ede)
## 2.5 % 97.5 %
## isq_factorSI -0.22620 0.643
## edad 0.00986 0.049
# Martingale y Deviance
ggcoxdiagnostics(cox_model_ede, type = "martingale", linear.predictions = FALSE)
## `geom_smooth()` using formula = 'y ~ x'
ggcoxdiagnostics(cox_model_ede, type = "deviance", linear.predictions = FALSE)
## `geom_smooth()` using formula = 'y ~ x'
#Test de proporcionalidad de_ riesgos de Schoenfeld
test_ph <- cox.zph(cox_model_ede)
test_ph
## chisq df p
## isq_factor 0.514 1 0.47
## edad 0.199 1 0.66
## GLOBAL 0.521 2 0.77
#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.
wilcox.test(sob1$edad~sob1$fallo)
##
## Wilcoxon rank sum test with continuity correction
##
## data: sob1$edad by sob1$fallo
## W = 4768, p-value = 0.0000008
## alternative hypothesis: true location shift is not equal to 0
library(psych)
##
## 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':
##
## %+%, alpha
describeBy(
x = sob1$edad,
group = sob1$fallo_factor,
mat = TRUE
)
## item group1 vars n mean sd median trimmed mad min max range skew
## X11 1 No fallo 1 134 54.5 10.8 54.5 54.5 9.64 25 89 64 0.0145
## X12 2 Fallo 1 112 61.7 11.1 61.5 61.5 9.64 30 89 59 0.1052
## kurtosis se
## X11 0.492 0.933
## X12 0.342 1.053
wilcox.test(sob1$edad~sob1$isq_factor)
##
## Wilcoxon rank sum test with continuity correction
##
## data: sob1$edad by sob1$isq_factor
## W = 3448, p-value = 0.0000000000006
## alternative hypothesis: true location shift is not equal to 0
library(psych)
describeBy(
x = sob1$edad,
group = sob1$isq_factor,
mat = TRUE
)
## item group1 vars n mean sd median trimmed mad min max range skew
## X11 1 NO/leve 1 140 53.2 9.82 54 53.6 8.90 25 77 52 -0.342
## X12 2 SI 1 106 63.8 10.82 63 63.5 9.64 35 89 54 0.203
## kurtosis se
## X11 0.3012 0.83
## X12 -0.0138 1.05
table(sob1$estado)
##
## ampu AMPU cerro CERRO FALLECIO murio MURIO PERDIDO
## 1 42 5 130 4 3 3 4
## persiste PERSISTE
## 4 51
sob1 <- sob1 %>%
mutate(
estado = tolower(trimws(estado))
)
sob1 <- sob1 %>%
mutate(
estado_unif = case_when(
estado %in% c("ampu") ~ "amputacion",
estado %in% c("cerro") ~ "cerro_herida",
estado %in% c("murio", "fallecio") ~ "fallecio",
estado %in% c("persiste") ~ "persiste",
estado %in% c("perdido") ~ "perdido_seguimiento",
TRUE ~ NA_character_
)
)
table(sob1$estado_unif)
##
## amputacion cerro_herida fallecio perdido_seguimiento
## 43 135 10 4
## persiste
## 55
kruskal.test(sob1$edad~sob1$estado_unif)
##
## Kruskal-Wallis rank sum test
##
## data: sob1$edad by sob1$estado_unif
## Kruskal-Wallis chi-squared = 29, df = 4, p-value = 0.00001
describeBy(
x = sob1$edad,
group = sob1$estado_unif,
mat = TRUE
)
## item group1 vars n mean sd median trimmed mad min max
## X11 1 amputacion 1 43 63.6 8.19 63.0 63.0 7.41 47 88
## X12 2 cerro_herida 1 134 54.5 10.79 54.5 54.5 9.64 25 89
## X13 3 fallecio 1 10 63.6 22.84 72.5 64.8 22.98 30 88
## X14 4 perdido_seguimiento 1 4 62.0 10.98 64.5 62.0 8.15 47 72
## X15 5 persiste 1 55 59.8 10.11 58.0 59.3 8.90 41 89
## range skew kurtosis se
## X11 41 0.6770 0.5553 1.250
## X12 64 0.0145 0.4918 0.933
## X13 58 -0.3333 -1.7265 7.222
## X14 25 -0.4074 -1.9368 5.492
## X15 48 0.5176 -0.0857 1.363
cox_model_edad <- coxph(tiempo_evento ~ edad , data = sob1)
summary(cox_model_edad)
## Call:
## coxph(formula = tiempo_evento ~ edad, data = sob1)
##
## n= 246, number of events= 112
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## edad 0.03425 1.03484 0.00865 3.96 0.000075 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## edad 1.03 0.966 1.02 1.05
##
## Concordance= 0.623 (se = 0.031 )
## Likelihood ratio test= 15.1 on 1 df, p=0.0001
## Wald test = 15.7 on 1 df, p=0.00008
## Score (logrank) test = 15.6 on 1 df, p=0.00008
tab_model(cox_model_edad,
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
|
|
edad
|
1.03
|
0.01
|
0.00 – Inf
|
<0.001
|
|
Observations
|
246
|
|
R2 Nagelkerke
|
0.061
|
confint(cox_model_edad)
## 2.5 % 97.5 %
## edad 0.0173 0.0512
#C de harrell
summary(cox_model)$concordance
## C se(C)
## 0.6114 0.0286
#C de harrell
summary(cox_model_ede)$concordance
## C se(C)
## 0.6282 0.0311
#C de harrell
#Creamos evento con 3 niveles Censura, Evento de Interes, Evento Competitivo
sob1 <- sob1 %>%
mutate(
evento_fallo = case_when(
estado_unif %in% c("amputacion", "persiste") ~ 1,
estado_unif %in% c("fallecio", "perdido_seguimiento") ~ 2,
estado_unif == "cerro_herida" ~ 0,
TRUE ~ NA_real_
)
)
table(sob1$evento_fallo)
##
## 0 1 2
## 135 98 14
#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_fallo <- Surv(sob1$dias_out, sob1$evento_fallo==1)
tiempo_muerte <- Surv(sob1$dias_out, sob1$evento_fallo==2)
cox_fallo <- coxph(tiempo_fallo ~
isq_factor + edad+
ede_factor , data = sob1)
summary(cox_fallo)
## Call:
## coxph(formula = tiempo_fallo ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 246, number of events= 98
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 0.4608 1.5853 0.2392 1.93 0.054 .
## edad 0.0258 1.0262 0.0111 2.33 0.020 *
## ede_factorUnilateral -0.1380 0.8711 0.3935 -0.35 0.726
## ede_factorSin/local -0.5398 0.5829 0.3877 -1.39 0.164
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 1.585 0.631 0.992 2.53
## edad 1.026 0.975 1.004 1.05
## ede_factorUnilateral 0.871 1.148 0.403 1.88
## ede_factorSin/local 0.583 1.716 0.273 1.25
##
## Concordance= 0.668 (se = 0.029 )
## Likelihood ratio test= 18.8 on 4 df, p=0.0008
## Wald test = 18.1 on 4 df, p=0.001
## Score (logrank) test = 18.5 on 4 df, p=0.001
cox_muerte <- coxph(tiempo_muerte ~
isq_factor + edad+
ede_factor , data = sob1)
summary(cox_muerte)
## Call:
## coxph(formula = tiempo_muerte ~ isq_factor + edad + ede_factor,
## data = sob1)
##
## n= 246, number of events= 14
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -1.0078 0.3650 0.6601 -1.53 0.127
## edad 0.0682 1.0706 0.0285 2.40 0.017 *
## ede_factorUnilateral 0.3000 1.3498 1.0822 0.28 0.782
## ede_factorSin/local -0.2036 0.8158 1.0805 -0.19 0.851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.365 2.739 0.1001 1.33
## edad 1.071 0.934 1.0125 1.13
## ede_factorUnilateral 1.350 0.741 0.1618 11.26
## ede_factorSin/local 0.816 1.226 0.0982 6.78
##
## Concordance= 0.676 (se = 0.095 )
## Likelihood ratio test= 6.99 on 4 df, p=0.1
## Wald test = 6.6 on 4 df, p=0.2
## Score (logrank) test = 6.69 on 4 df, p=0.2
#CIF
cifivhx <- cuminc(Surv(dias_out, factor(evento_fallo)) ~ isq_factor, data=sob1)
cifivhx
##
## ── cuminc() ────────────────────────────────────────────────────────────────────
## • Failure type "1"
## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 95 0.064 0.021 0.032, 0.114
## NO/leve 100 59 0.145 0.035 0.085, 0.222
## NO/leve 150 33 0.274 0.052 0.178, 0.379
## NO/leve 200 18 0.420 0.068 0.286, 0.548
## NO/leve 250 3 0.696 0.084 0.498, 0.828
## NO/leve 300 2 0.696 0.084 0.498, 0.828
## SI 50.0 75 0.247 0.042 0.169, 0.333
## SI 100 52 0.333 0.047 0.243, 0.426
## SI 150 28 0.476 0.055 0.365, 0.580
## SI 200 15 0.633 0.061 0.501, 0.740
## SI 250 6 0.796 0.057 0.655, 0.885
## SI 300 4 0.796 0.057 0.655, 0.885
## • Failure type "2"
## strata time n.risk estimate std.error 95% CI
## NO/leve 50.0 95 0.065 0.021 0.032, 0.115
## NO/leve 100 59 0.065 0.021 0.032, 0.115
## NO/leve 150 33 0.065 0.021 0.032, 0.115
## NO/leve 200 18 0.065 0.021 0.032, 0.115
## NO/leve 250 3 0.065 0.021 0.032, 0.115
## NO/leve 300 2 0.065 0.021 0.032, 0.115
## SI 50.0 75 0.028 0.016 0.008, 0.074
## SI 100 52 0.028 0.016 0.008, 0.074
## SI 150 28 0.055 0.025 0.020, 0.118
## SI 200 15 0.055 0.025 0.020, 0.118
## SI 250 6 0.055 0.025 0.020, 0.118
## SI 300 4 0.055 0.025 0.020, 0.118
## • Tests
## outcome statistic df p.value
## 1 10.7 1.00 0.001
## 2 0.512 1.00 0.47
levels(sob1$isq_factor)
## [1] "NO/leve" "SI"
table(sob1$isq_factor)
##
## NO/leve SI
## 141 106
# kaplan meier usando 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 interés (amputación / persiste)",
"Evento competitivo (muerte / pérdida)"
)
) +
labs(
x = "Tiempo (días)",
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_CSC <- CSC(
formula = Hist(dias_out, evento_fallo) ~ 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_fallo) ~ isq_factor + edad +
## ede_factor, data = sob1, cause = 1)
##
## Right-censored response of a competing.risks model
##
## No.Observations: 247
##
## Pattern:
##
## Cause event right.censored
## 1 98 0
## 2 14 0
## unknown 0 135
##
##
## ----------> Cause: 1
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 246, number of events= 98
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI 0.4608 1.5853 0.2392 1.93 0.054 .
## edad 0.0258 1.0262 0.0111 2.33 0.020 *
## ede_factorUnilateral -0.1380 0.8711 0.3935 -0.35 0.726
## ede_factorSin/local -0.5398 0.5829 0.3877 -1.39 0.164
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 1.585 0.631 0.992 2.53
## edad 1.026 0.975 1.004 1.05
## ede_factorUnilateral 0.871 1.148 0.403 1.88
## ede_factorSin/local 0.583 1.716 0.273 1.25
##
## Concordance= 0.668 (se = 0.029 )
## Likelihood ratio test= 18.8 on 4 df, p=0.0008
## Wald test = 18.1 on 4 df, p=0.001
## Score (logrank) test = 18.5 on 4 df, p=0.001
##
##
##
## ----------> Cause: 2
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ isq_factor + edad +
## ede_factor, x = TRUE, y = TRUE)
##
## n= 246, number of events= 14
## (1 observation deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## isq_factorSI -1.0078 0.3650 0.6601 -1.53 0.127
## edad 0.0682 1.0706 0.0285 2.40 0.017 *
## ede_factorUnilateral 0.3000 1.3498 1.0822 0.28 0.782
## ede_factorSin/local -0.2036 0.8158 1.0805 -0.19 0.851
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## isq_factorSI 0.365 2.739 0.1001 1.33
## edad 1.071 0.934 1.0125 1.13
## ede_factorUnilateral 1.350 0.741 0.1618 11.26
## ede_factorSin/local 0.816 1.226 0.0982 6.78
##
## Concordance= 0.676 (se = 0.095 )
## Likelihood ratio test= 6.99 on 4 df, p=0.1
## Wald test = 6.6 on 4 df, p=0.2
## Score (logrank) test = 6.69 on 4 df, p=0.2
sob1$edad[104] <- 55 # reemplazá 65 por la edad correcta
modelo_fg <- FGR(
formula = Hist(dias_out, evento_fallo) ~ isq_factor + ede_factor + edad
, data = sob1,
cause = 1 # Aclara que CHD es el evento de interés
)
modelo_fg
##
## Right-censored response of a competing.risks model
##
## No.Observations: 247
##
## Pattern:
##
## Cause event right.censored
## 1 98 0
## 2 14 0
## unknown 0 135
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(dias_out, evento_fallo) ~ isq_factor + ede_factor +
## edad, data = sob1, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## isq_factorSI 0.6052 1.832 0.2311 2.618 0.0088
## ede_factorUnilateral -0.1674 0.846 0.3269 -0.512 0.6100
## ede_factorSin/local -0.5491 0.577 0.3100 -1.771 0.0770
## edad 0.0144 1.015 0.0091 1.583 0.1100
##
## exp(coef) exp(-coef) 2.5% 97.5%
## isq_factorSI 1.832 0.546 1.164 2.88
## ede_factorUnilateral 0.846 1.182 0.446 1.61
## ede_factorSin/local 0.577 1.732 0.315 1.06
## edad 1.015 0.986 0.997 1.03
##
## Num. cases = 247
## Pseudo Log-likelihood = -443
## Pseudo likelihood ratio test = 16.6 on 4 df,
##
## Convergence: TRUE
# Evaluar performance predictiva con Score
score_fg <- Score(
object = list("FG model" = modelo_fg),
formula = Hist(dias_out, evento_fallo) ~ 1,
data = sob1,
times = c(1:24),
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 1 0.822 0.0831 0.659 0.985
## 2: FG model 2 0.822 0.0831 0.659 0.985
## 3: FG model 3 0.822 0.0831 0.659 0.985
## 4: FG model 4 0.822 0.0831 0.659 0.985
## 5: FG model 5 0.822 0.0831 0.659 0.985
## 6: FG model 6 0.822 0.0831 0.659 0.985
## 7: FG model 7 0.711 0.1076 0.500 0.922
## 8: FG model 8 0.711 0.1076 0.500 0.922
## 9: FG model 9 0.711 0.1076 0.500 0.922
## 10: FG model 10 0.711 0.0819 0.550 0.871
## 11: FG model 11 0.694 0.0761 0.545 0.843
## 12: FG model 12 0.694 0.0761 0.545 0.843
## 13: FG model 13 0.740 0.0766 0.590 0.890
## 14: FG model 14 0.704 0.0748 0.557 0.850
## 15: FG model 15 0.692 0.0684 0.558 0.826
## 16: FG model 16 0.692 0.0684 0.558 0.826
## 17: FG model 17 0.697 0.0642 0.571 0.823
## 18: FG model 18 0.722 0.0596 0.605 0.838
## 19: FG model 19 0.743 0.0544 0.636 0.849
## 20: FG model 20 0.743 0.0544 0.636 0.849
## 21: FG model 21 0.751 0.0525 0.649 0.854
## 22: FG model 22 0.731 0.0541 0.625 0.837
## 23: FG model 23 0.720 0.0533 0.615 0.824
## 24: FG model 24 0.717 0.0517 0.615 0.818
## model times AUC se lower upper
#score comparado con ipa
score_comp <-Score(
list("FG model" = modelo_fg, "CSC model" = modelo_CSC),
formula = Hist(dias_out, evento_fallo) ~ 1,
data = sob1,
times = c(1:24),
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 1 0.00000
## 2 Null model 2 0.00000
## 3 Null model 3 0.00000
## 4 Null model 4 0.00000
## 5 Null model 5 0.00000
## 6 Null model 6 0.00000
## 7 Null model 7 0.00000
## 8 Null model 8 0.00000
## 9 Null model 9 0.00000
## 10 Null model 10 0.00000
## 11 Null model 11 0.00000
## 12 Null model 12 0.00000
## 13 Null model 13 0.00000
## 14 Null model 14 0.00000
## 15 Null model 15 0.00000
## 16 Null model 16 0.00000
## 17 Null model 17 0.00000
## 18 Null model 18 0.00000
## 19 Null model 19 0.00000
## 20 Null model 20 0.00000
## 21 Null model 21 0.00000
## 22 Null model 22 0.00000
## 23 Null model 23 0.00000
## 24 Null model 24 0.00000
## 25 FG model 1 0.00608
## 26 FG model 2 0.00608
## 27 FG model 3 0.00608
## 28 FG model 4 0.00608
## 29 FG model 5 0.00608
## 30 FG model 6 0.00608
## 31 FG model 7 0.00462
## 32 FG model 8 0.00462
## 33 FG model 9 0.00462
## 34 FG model 10 0.00529
## 35 FG model 11 0.00727
## 36 FG model 12 0.00727
## 37 FG model 13 0.01835
## 38 FG model 14 0.02225
## 39 FG model 15 0.02490
## 40 FG model 16 0.02490
## 41 FG model 17 0.02627
## 42 FG model 18 0.03883
## 43 FG model 19 0.05404
## 44 FG model 20 0.05404
## 45 FG model 21 0.05852
## 46 FG model 22 0.05474
## 47 FG model 23 0.05286
## 48 FG model 24 0.05280
## 49 CSC model 1 0.00530
## 50 CSC model 2 0.00530
## 51 CSC model 3 0.00530
## 52 CSC model 4 0.00530
## 53 CSC model 5 0.00530
## 54 CSC model 6 0.00530
## 55 CSC model 7 0.00435
## 56 CSC model 8 0.00435
## 57 CSC model 9 0.00435
## 58 CSC model 10 0.00453
## 59 CSC model 11 0.00548
## 60 CSC model 12 0.00548
## 61 CSC model 13 0.02051
## 62 CSC model 14 0.02331
## 63 CSC model 15 0.02229
## 64 CSC model 16 0.02229
## 65 CSC model 17 0.02362
## 66 CSC model 18 0.03795
## 67 CSC model 19 0.05578
## 68 CSC model 20 0.05578
## 69 CSC model 21 0.05928
## 70 CSC model 22 0.05571
## 71 CSC model 23 0.05382
## 72 CSC model 24 0.05311
# 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")
# Interpretamos que a medida que pasa el tiempo el modelo logra separar mejor,
#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.5) +
geom_hline(yintercept = 0.7, linetype = "dotted", alpha = 0.7) +
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_fallo) ~ 1,
data = sob1,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(s, models = "FGR" )
