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
# 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