Sobrevida
2025-08-09
##librarys
library(naniar)## Warning: package 'naniar' was built under R version 4.4.3library(dplyr)##
## Adjuntando el paquete: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(survival)
library(survminer)## Cargando paquete requerido: ggplot2## Cargando paquete requerido: ggpubr##
## 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:dplyr':
##
## selectlibrary(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(sjPlot)
options(scipen = 999, digits = 3, encoding = 'UTF-8')
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)
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, cumincR Markdown
library(readr)##
## Adjuntando el paquete: 'readr'## The following object is masked from 'package:scales':
##
## col_factordrugs <- read_csv("C:/Users/Administrador/Downloads/drugs.csv")## Rows: 628 Columns: 13## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (13): id, age, beck, hercoc, ivhx, ndrugtx, race, treat, sex, relapse, t...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.View(drugs)
#backup
drugs1=drugs##Datos faltantes
vis_miss(drugs)
##Pasar a factor
#drugs$sex
drugs$sex <- factor(drugs1$sex, levels = c(0, 1), labels = c( "Mujer", "Hombre"))
levels(drugs$sex)## [1] "Mujer" "Hombre"summary(drugs$sex)## Mujer Hombre
## 184 444table(drugs$sex)##
## Mujer Hombre
## 184 444#drugs$hercoc
drugs$hercoc <- factor(drugs$hercoc, levels = c(1, 2,3,4), labels = c( "Her-ox", "Her", "Ox", "Coc"))
levels(drugs$hercoc)## [1] "Her-ox" "Her" "Ox" "Coc"summary(drugs$hercoc)## Her-ox Her Ox Coc
## 123 118 182 205#ivhx
drugs$ivhx <- factor(drugs$ivhx, levels = c(1, 2,3), labels = c( "Nunca", "A veces", "Siempre"))
levels(drugs$ivhx)## [1] "Nunca" "A veces" "Siempre"summary(drugs$ivhx)## Nunca A veces Siempre
## 161 122 345#race
drugs$race<- factor(drugs1$race, levels = c(0, 1), labels = c( "Blanca", "Negra_hisp"))
levels(drugs$race)## [1] "Blanca" "Negra_hisp"summary(drugs$race)## Blanca Negra_hisp
## 189 439#treat
drugs$treat<- factor(drugs$treat, levels = c(0, 1), labels = c( "No", "Si"))
levels(drugs$treat)## [1] "No" "Si"summary(drugs$treat)## No Si
## 321 307#fentanyl
drugs$fentanyl<- factor(drugs$fentanyl, levels = c(0, 1), labels = c( "No", "Si"))
levels(drugs$fentanyl)## [1] "No" "Si"summary(drugs$fentanyl)## No Si
## 398 230#objeto de sobrevida
tiempo_evento <- Surv(drugs$time,drugs$relapse)##variables
summary(drugs)## id age beck hercoc ivhx
## Min. : 1 Min. :20.0 Min. : 0.0 Her-ox:123 Nunca :161
## 1st Qu.:158 1st Qu.:27.0 1st Qu.:10.0 Her :118 A veces:122
## Median :314 Median :32.0 Median :17.0 Ox :182 Siempre:345
## Mean :314 Mean :32.3 Mean :17.7 Coc :205
## 3rd Qu.:471 3rd Qu.:37.0 3rd Qu.:24.0
## Max. :628 Max. :56.0 Max. :54.0
## ndrugtx race treat sex relapse
## Min. : 0.0 Blanca :189 No:321 Mujer :184 Min. :0.000
## 1st Qu.: 1.0 Negra_hisp:439 Si:307 Hombre:444 1st Qu.:1.000
## Median : 3.0 Median :1.000
## Mean : 4.6 Mean :0.809
## 3rd Qu.: 6.0 3rd Qu.:1.000
## Max. :40.0 Max. :1.000
## time fentanyl death
## Min. : 2 No:398 Min. :0.000
## 1st Qu.: 74 Si:230 1st Qu.:0.000
## Median :159 Median :0.000
## Mean :224 Mean :0.073
## 3rd Qu.:340 3rd Qu.:0.000
## Max. :721 Max. :1.000##modelo full
cox_model_full <- coxph(tiempo_evento ~ age+beck+hercoc+ivhx+ndrugtx+race+treat+sex+fentanyl, data = drugs)
summary(cox_model_full)## Call:
## coxph(formula = tiempo_evento ~ age + beck + hercoc + ivhx +
## ndrugtx + race + treat + sex + fentanyl, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age -0.02890 0.97151 0.00781 -3.70 0.00022 ***
## beck 0.00828 1.00831 0.00464 1.78 0.07462 .
## hercocHer 0.12180 1.12953 0.14097 0.86 0.38755
## hercocOx 0.18460 1.20273 0.14980 1.23 0.21784
## hercocCoc 0.24164 1.27334 0.14951 1.62 0.10604
## ivhxA veces 0.58359 1.79247 0.15258 3.82 0.00013 ***
## ivhxSiempre 0.77523 2.17110 0.16189 4.79 0.0000017 ***
## ndrugtx 0.02767 1.02805 0.00804 3.44 0.00058 ***
## raceNegra_hisp -0.18685 0.82957 0.09563 -1.95 0.05072 .
## treatSi -0.27576 0.75900 0.09024 -3.06 0.00224 **
## sexHombre 0.08146 1.08487 0.10367 0.79 0.43199
## fentanylSi 0.16239 1.17632 0.11120 1.46 0.14419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 0.972 1.029 0.957 0.987
## beck 1.008 0.992 0.999 1.018
## hercocHer 1.130 0.885 0.857 1.489
## hercocOx 1.203 0.831 0.897 1.613
## hercocCoc 1.273 0.785 0.950 1.707
## ivhxA veces 1.792 0.558 1.329 2.417
## ivhxSiempre 2.171 0.461 1.581 2.982
## ndrugtx 1.028 0.973 1.012 1.044
## raceNegra_hisp 0.830 1.205 0.688 1.001
## treatSi 0.759 1.318 0.636 0.906
## sexHombre 1.085 0.922 0.885 1.329
## fentanylSi 1.176 0.850 0.946 1.463
##
## Concordance= 0.626 (se = 0.013 )
## Likelihood ratio test= 99.4 on 12 df, p=0.0000000000000007
## Wald test = 94.9 on 12 df, p=0.000000000000006
## Score (logrank) test = 98.5 on 12 df, p=0.000000000000001##modelo automático
exp(cbind(HR = coef(cox_model_full), confint(cox_model_full)))## HR 2.5 % 97.5 %
## age 0.972 0.957 0.987
## beck 1.008 0.999 1.018
## hercocHer 1.130 0.857 1.489
## hercocOx 1.203 0.897 1.613
## hercocCoc 1.273 0.950 1.707
## ivhxA veces 1.792 1.329 2.417
## ivhxSiempre 2.171 1.581 2.982
## ndrugtx 1.028 1.012 1.044
## raceNegra_hisp 0.830 0.688 1.001
## treatSi 0.759 0.636 0.906
## sexHombre 1.085 0.885 1.329
## fentanylSi 1.176 0.946 1.463#selección de modelo step
#Selección de vaariables para modelo final
cox_model <- stepAIC(cox_model_full, direction = "both")## Start: AIC=5810
## tiempo_evento ~ age + beck + hercoc + ivhx + ndrugtx + race +
## treat + sex + fentanyl
##
## Df AIC
## - hercoc 3 5807
## - sex 1 5809
## <none> 5810
## - fentanyl 1 5810
## - beck 1 5811
## - race 1 5812
## - treat 1 5817
## - ndrugtx 1 5818
## - age 1 5822
## - ivhx 2 5831
##
## Step: AIC=5807
## tiempo_evento ~ age + beck + ivhx + ndrugtx + race + treat +
## sex + fentanyl
##
## Df AIC
## - sex 1 5805
## <none> 5807
## - beck 1 5807
## - race 1 5808
## - fentanyl 1 5809
## + hercoc 3 5810
## - treat 1 5813
## - ndrugtx 1 5815
## - age 1 5820
## - ivhx 2 5827
##
## Step: AIC=5805
## tiempo_evento ~ age + beck + ivhx + ndrugtx + race + treat +
## fentanyl
##
## Df AIC
## <none> 5805
## - beck 1 5806
## - race 1 5806
## + sex 1 5807
## - fentanyl 1 5807
## + hercoc 3 5809
## - treat 1 5811
## - ndrugtx 1 5813
## - age 1 5818
## - ivhx 2 5828summary(cox_model)## Call:
## coxph(formula = tiempo_evento ~ age + beck + ivhx + ndrugtx +
## race + treat + fentanyl, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age -0.02966 0.97078 0.00770 -3.85 0.00012 ***
## beck 0.00770 1.00773 0.00460 1.67 0.09447 .
## ivhxA veces 0.55771 1.74667 0.15074 3.70 0.00022 ***
## ivhxSiempre 0.68468 1.98313 0.13608 5.03 0.00000049 ***
## ndrugtx 0.02702 1.02739 0.00807 3.35 0.00082 ***
## raceNegra_hisp -0.17905 0.83606 0.09557 -1.87 0.06099 .
## treatSi -0.25751 0.77297 0.08941 -2.88 0.00397 **
## fentanylSi 0.21233 1.23656 0.10239 2.07 0.03810 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 0.971 1.030 0.956 0.986
## beck 1.008 0.992 0.999 1.017
## ivhxA veces 1.747 0.573 1.300 2.347
## ivhxSiempre 1.983 0.504 1.519 2.589
## ndrugtx 1.027 0.973 1.011 1.044
## raceNegra_hisp 0.836 1.196 0.693 1.008
## treatSi 0.773 1.294 0.649 0.921
## fentanylSi 1.237 0.809 1.012 1.511
##
## Concordance= 0.624 (se = 0.013 )
## Likelihood ratio test= 96.2 on 8 df, p=<0.0000000000000002
## Wald test = 91.9 on 8 df, p=<0.0000000000000002
## Score (logrank) test = 95.4 on 8 df, p=<0.0000000000000002tab_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 |
| age | 0.97 | 0.01 | 0.00 – Inf | <0.001 |
| beck | 1.01 | 0.00 | 0.00 – Inf | 0.094 |
| ivhx [A veces] | 1.75 | 0.26 | 0.00 – Inf | <0.001 |
| ivhx [Siempre] | 1.98 | 0.27 | 0.00 – Inf | <0.001 |
| ndrugtx | 1.03 | 0.01 | 0.00 – Inf | 0.001 |
| race [Negra_hisp] | 0.84 | 0.08 | 0.00 – Inf | 0.061 |
| treat [Si] | 0.77 | 0.07 | 0.00 – Inf | 0.004 |
| fentanyl [Si] | 1.24 | 0.13 | 0.00 – Inf | 0.038 |
| Observations | 628 | |||
| R2 Nagelkerke | 0.142 | |||
confint(cox_model)## 2.5 % 97.5 %
## age -0.04475 -0.01456
## beck -0.00132 0.01672
## ivhxA veces 0.26227 0.85315
## ivhxSiempre 0.41797 0.95139
## ndrugtx 0.01119 0.04285
## raceNegra_hisp -0.36636 0.00826
## treatSi -0.43274 -0.08228
## fentanylSi 0.01165 0.41301AIC(cox_model) #5805 cox_model## [1] 5805##Que oasa si cambio endovenoso de categoría. Es decir, que nunca se drogó endovenoso
drugs$evno=drugs$ivhx
levels(drugs$evno)## [1] "Nunca" "A veces" "Siempre"table(drugs$evno)##
## Nunca A veces Siempre
## 161 122 345# Reordenar niveles para que "Siempre" sea la referencia
drugs$evno <- factor(drugs$evno,
levels = c("Siempre", "Nunca", "A veces"))
table(drugs$evno)##
## Siempre Nunca A veces
## 345 161 122levels(drugs$evno)## [1] "Siempre" "Nunca" "A veces"##nuevo modelo reordenado
cox_model_full2 <- coxph(tiempo_evento ~ age+beck+hercoc+evno+ndrugtx+race+treat+sex+fentanyl, data = drugs)
cox_model2 <- stepAIC(cox_model_full2, direction = "both")## Start: AIC=5810
## tiempo_evento ~ age + beck + hercoc + evno + ndrugtx + race +
## treat + sex + fentanyl
##
## Df AIC
## - hercoc 3 5807
## - sex 1 5809
## <none> 5810
## - fentanyl 1 5810
## - beck 1 5811
## - race 1 5812
## - treat 1 5817
## - ndrugtx 1 5818
## - age 1 5822
## - evno 2 5831
##
## Step: AIC=5807
## tiempo_evento ~ age + beck + evno + ndrugtx + race + treat +
## sex + fentanyl
##
## Df AIC
## - sex 1 5805
## <none> 5807
## - beck 1 5807
## - race 1 5808
## - fentanyl 1 5809
## + hercoc 3 5810
## - treat 1 5813
## - ndrugtx 1 5815
## - age 1 5820
## - evno 2 5827
##
## Step: AIC=5805
## tiempo_evento ~ age + beck + evno + ndrugtx + race + treat +
## fentanyl
##
## Df AIC
## <none> 5805
## - beck 1 5806
## - race 1 5806
## + sex 1 5807
## - fentanyl 1 5807
## + hercoc 3 5809
## - treat 1 5811
## - ndrugtx 1 5813
## - age 1 5818
## - evno 2 5828summary(cox_model2)## Call:
## coxph(formula = tiempo_evento ~ age + beck + evno + ndrugtx +
## race + treat + fentanyl, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age -0.02966 0.97078 0.00770 -3.85 0.00012 ***
## beck 0.00770 1.00773 0.00460 1.67 0.09447 .
## evnoNunca -0.68468 0.50425 0.13608 -5.03 0.00000049 ***
## evnoA veces -0.12697 0.88076 0.11801 -1.08 0.28199
## ndrugtx 0.02702 1.02739 0.00807 3.35 0.00082 ***
## raceNegra_hisp -0.17905 0.83606 0.09557 -1.87 0.06099 .
## treatSi -0.25751 0.77297 0.08941 -2.88 0.00397 **
## fentanylSi 0.21233 1.23656 0.10239 2.07 0.03810 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 0.971 1.030 0.956 0.986
## beck 1.008 0.992 0.999 1.017
## evnoNunca 0.504 1.983 0.386 0.658
## evnoA veces 0.881 1.135 0.699 1.110
## ndrugtx 1.027 0.973 1.011 1.044
## raceNegra_hisp 0.836 1.196 0.693 1.008
## treatSi 0.773 1.294 0.649 0.921
## fentanylSi 1.237 0.809 1.012 1.511
##
## Concordance= 0.624 (se = 0.013 )
## Likelihood ratio test= 96.2 on 8 df, p=<0.0000000000000002
## Wald test = 91.9 on 8 df, p=<0.0000000000000002
## Score (logrank) test = 95.4 on 8 df, p=<0.0000000000000002tab_model(cox_model2,
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 |
| age | 0.97 | 0.01 | 0.00 – Inf | <0.001 |
| beck | 1.01 | 0.00 | 0.00 – Inf | 0.094 |
| evno [Nunca] | 0.50 | 0.07 | 0.00 – Inf | <0.001 |
| evno [A veces] | 0.88 | 0.10 | 0.00 – Inf | 0.282 |
| ndrugtx | 1.03 | 0.01 | 0.00 – Inf | 0.001 |
| race [Negra_hisp] | 0.84 | 0.08 | 0.00 – Inf | 0.061 |
| treat [Si] | 0.77 | 0.07 | 0.00 – Inf | 0.004 |
| fentanyl [Si] | 1.24 | 0.13 | 0.00 – Inf | 0.038 |
| Observations | 628 | |||
| R2 Nagelkerke | 0.142 | |||
confint(cox_model2)## 2.5 % 97.5 %
## age -0.04475 -0.01456
## beck -0.00132 0.01672
## evnoNunca -0.95139 -0.41797
## evnoA veces -0.35827 0.10434
## ndrugtx 0.01119 0.04285
## raceNegra_hisp -0.36636 0.00826
## treatSi -0.43274 -0.08228
## fentanylSi 0.01165 0.41301AIC(cox_model2) #AIC 5805 modelo con cambio de levels## [1] 5805Las personas más jóvenes, quienes nunca usaron drogas endovenosas y aquellos que hicieron tratamiento son los que tienen menos riesgo de recaída
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph <- cox.zph(cox_model2)
test_ph## chisq df p
## age 0.177 1 0.674
## beck 5.571 1 0.018
## evno 0.209 2 0.901
## ndrugtx 0.794 1 0.373
## race 3.267 1 0.071
## treat 3.569 1 0.059
## fentanyl 0.752 1 0.386
## GLOBAL 14.552 8 0.068#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph)
# Gráficos log(-log) para la variable evno
ggsurvplot(survfit(tiempo_evento ~ evno, data = drugs), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
#Gráfico loglog variable age
# Gráficos log(-log) para la variable sex
ggsurvplot(survfit(tiempo_evento ~ race, data = drugs), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
#Gráfico loglog variable treat
ggsurvplot(survfit(tiempo_evento ~ treat, data = drugs), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
#Gráfico loglog variable fentanyl
ggsurvplot(survfit(tiempo_evento ~ fentanyl, data = drugs), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
#cómo se comporta la edad con el evento?
drugs<- drugs %>% mutate(age_q= ntile(age, 5))
drugs %>% dplyr::group_by(age_q) %>% dplyr::summarise(tar=mean(relapse)) %>% ggplot(aes(x=age_q, y=tar))+geom_point()+geom_smooth(method = "lm", col = "blue")## `geom_smooth()` using formula = 'y ~ x'
drugs %>%
dplyr::group_by(age_q) %>%
dplyr::summarise(tar = mean(relapse)) %>%
ggplot(aes(x = age_q, y = tar)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
labs(x = "Cuartiles de colesterol", y = "Proporción de eventos CHD") +
theme_minimal()
#dicotomizar edad en 4to quintilo
quantile(drugs$age, probs = 0.6, na.rm = TRUE)## 60%
## 34#como se comporta beck con el evento: beck prodría dicotomizarse en el quintilo 3
drugs<- drugs %>% mutate(beck_q= ntile(beck, 5))
drugs %>% dplyr::group_by(beck_q) %>% dplyr::summarise(tar=mean(relapse)) %>% ggplot(aes(x=beck_q, y=tar))+geom_point()+geom_smooth(method = "lm", col = "blue")## `geom_smooth()` using formula = 'y ~ x'
drugs %>%
dplyr::group_by(beck_q) %>%
dplyr::summarise(tar = mean(relapse)) %>%
ggplot(aes(x = beck_q, y = tar)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
labs(x = "Cuartiles de colesterol", y = "Proporción de eventos CHD") +
theme_minimal()
#cuál es el quintilo 3 de beck
# Calcular el 3er quintil (60%)
quantile(drugs$beck, probs = 0.6, na.rm = TRUE)## 60%
## 20#El tercer cuantil vale 20. Dicotomizo la variable
drugs$beck_cat <- as.factor(ifelse(drugs$beck>=20,"alto","bajo"))
table(drugs$beck_cat)##
## alto bajo
## 262 366#como se comporta ndrugtx con el evento
drugs<- drugs %>% mutate(ndrugtx_q= ntile(ndrugtx, 5))
drugs %>% dplyr::group_by(ndrugtx_q) %>% dplyr::summarise(tar=mean(relapse)) %>% ggplot(aes(x=ndrugtx_q, y=tar))+geom_point()+geom_smooth(method = "lm", col = "blue")## `geom_smooth()` using formula = 'y ~ x'
drugs %>%
dplyr::group_by(ndrugtx_q) %>%
dplyr::summarise(tar = mean(relapse)) %>%
ggplot(aes(x = ndrugtx_q, y = tar)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
labs(x = "Cuartiles de ndrugtx", y = "Proporción de eventos CHD") +
theme_minimal()
drugs <- drugs |>
mutate(ndrugtx_cat2 = factor(
if_else(ndrugtx >= 4, "ndrugtx_mas_4", "ndrugtx_menos_4"),
levels = c("ndrugtx_mas_4", "ndrugtx_menos_4")
))
#test schoenfeld 0.068 global. beck viola el supuesto es menor a 0.05##modelo 3 con beck dicotómica: tiene menor AIC y la p de schoenfeld en mayor a 0.05 para todas las variables. la p global de 0.23 en vez de 0.06
cox_model_full3 <- coxph(tiempo_evento ~ age+beck_cat+hercoc+evno+ndrugtx+race+treat+sex+fentanyl, data = drugs)
cox_model3 <- stepAIC(cox_model_full3, direction = "both")## Start: AIC=5808
## tiempo_evento ~ age + beck_cat + hercoc + evno + ndrugtx + race +
## treat + sex + fentanyl
##
## Df AIC
## - hercoc 3 5805
## - sex 1 5807
## <none> 5808
## - fentanyl 1 5808
## - race 1 5810
## - beck_cat 1 5811
## - treat 1 5815
## - ndrugtx 1 5817
## - age 1 5820
## - evno 2 5829
##
## Step: AIC=5805
## tiempo_evento ~ age + beck_cat + evno + ndrugtx + race + treat +
## sex + fentanyl
##
## Df AIC
## - sex 1 5804
## <none> 5805
## - race 1 5807
## - beck_cat 1 5807
## - fentanyl 1 5807
## + hercoc 3 5808
## - treat 1 5811
## - ndrugtx 1 5813
## - age 1 5818
## - evno 2 5826
##
## Step: AIC=5804
## tiempo_evento ~ age + beck_cat + evno + ndrugtx + race + treat +
## fentanyl
##
## Df AIC
## <none> 5804
## + sex 1 5805
## - race 1 5805
## - beck_cat 1 5806
## - fentanyl 1 5806
## + hercoc 3 5807
## - treat 1 5809
## - ndrugtx 1 5812
## - age 1 5816
## - evno 2 5826summary(cox_model3)## Call:
## coxph(formula = tiempo_evento ~ age + beck_cat + evno + ndrugtx +
## race + treat + fentanyl, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age -0.02915 0.97127 0.00769 -3.79 0.00015 ***
## beck_catbajo -0.18440 0.83161 0.09046 -2.04 0.04151 *
## evnoNunca -0.68836 0.50240 0.13595 -5.06 0.00000041 ***
## evnoA veces -0.12902 0.87895 0.11794 -1.09 0.27396
## ndrugtx 0.02720 1.02758 0.00812 3.35 0.00081 ***
## raceNegra_hisp -0.18098 0.83445 0.09559 -1.89 0.05832 .
## treatSi -0.24573 0.78213 0.08967 -2.74 0.00614 **
## fentanylSi 0.21028 1.23402 0.10230 2.06 0.03984 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 0.971 1.030 0.957 0.986
## beck_catbajo 0.832 1.202 0.696 0.993
## evnoNunca 0.502 1.990 0.385 0.656
## evnoA veces 0.879 1.138 0.698 1.108
## ndrugtx 1.028 0.973 1.011 1.044
## raceNegra_hisp 0.834 1.198 0.692 1.006
## treatSi 0.782 1.279 0.656 0.932
## fentanylSi 1.234 0.810 1.010 1.508
##
## Concordance= 0.624 (se = 0.013 )
## Likelihood ratio test= 97.5 on 8 df, p=<0.0000000000000002
## Wald test = 93.5 on 8 df, p=<0.0000000000000002
## Score (logrank) test = 97 on 8 df, p=<0.0000000000000002tab_model(cox_model3,
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 |
| age | 0.97 | 0.01 | 0.00 – Inf | <0.001 |
| beck cat [bajo] | 0.83 | 0.08 | 0.00 – Inf | 0.042 |
| evno [Nunca] | 0.50 | 0.07 | 0.00 – Inf | <0.001 |
| evno [A veces] | 0.88 | 0.10 | 0.00 – Inf | 0.274 |
| ndrugtx | 1.03 | 0.01 | 0.00 – Inf | 0.001 |
| race [Negra_hisp] | 0.83 | 0.08 | 0.00 – Inf | 0.058 |
| treat [Si] | 0.78 | 0.07 | 0.00 – Inf | 0.006 |
| fentanyl [Si] | 1.23 | 0.13 | 0.00 – Inf | 0.040 |
| Observations | 628 | |||
| R2 Nagelkerke | 0.144 | |||
confint(cox_model3)## 2.5 % 97.5 %
## age -0.04422 -0.01409
## beck_catbajo -0.36170 -0.00709
## evnoNunca -0.95481 -0.42191
## evnoA veces -0.36018 0.10213
## ndrugtx 0.01128 0.04312
## raceNegra_hisp -0.36833 0.00637
## treatSi -0.42148 -0.06998
## fentanylSi 0.00976 0.41079AIC(cox_model3) #5804 y beck se hce significativa cuando es menor de 20## [1] 5804En este caso, los pacientes con menos recaídas serían aquellos con beck menor de 20, que nunca se drogaron por vía endovenosa, que hicieron tratamiento internados y más jóvenes
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph3 <- cox.zph(cox_model3)
test_ph3## chisq df p
## age 0.163 1 0.686
## beck_cat 2.190 1 0.139
## evno 0.218 2 0.897
## ndrugtx 0.784 1 0.376
## race 3.244 1 0.072
## treat 3.518 1 0.061
## fentanyl 0.705 1 0.401
## GLOBAL 10.395 8 0.238#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph3)
##curvas de kaplan meyer para el modelo 3
# Graficar las curvas de Kaplan-Meier estratificadas por beck
km_global <- survfit(tiempo_evento ~ 1, data = drugs)#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 ~ beck_cat, data = drugs)
ggsurvplot(km_fit, data = drugs,
risk.table = TRUE, # ver tabla de pacientes en riesgo
pval = TRUE,# agrega p-valor de logrank
conf.int = TRUE,# Intervalo de confianza 95% por default
cumcensor = TRUE,# tabla de censurados
ggtheme = theme_minimal(),
palette = "simpsons")
#logrank
logrank <- survdiff(Surv(time, relapse) ~ beck_cat, data = drugs)
logrank## Call:
## survdiff(formula = Surv(time, relapse) ~ beck_cat, data = drugs)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## beck_cat=alto 262 221 190 4.97 8
## beck_cat=bajo 366 287 318 2.98 8
##
## Chisq= 8 on 1 degrees of freedom, p= 0.005# Graficar las curvas de Kaplan-Meier estratificadas por evno
km_fit1 <- survfit(tiempo_evento ~ evno, data = drugs)
ggsurvplot(km_fit1, data = drugs,
risk.table = TRUE, # ver tabla de pacientes en riesgo
pval = TRUE,# agrega p-valor de logrank
conf.int = TRUE,# Intervalo de confianza 95% por default
cumcensor = TRUE,# tabla de censurados
ggtheme = theme_minimal(),
palette = "simpsons")
#logrank
logrank <- survdiff(Surv(time, relapse) ~ evno, data = drugs)
logrank## Call:
## survdiff(formula = Surv(time, relapse) ~ evno, data = drugs)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## evno=Siempre 345 309 241.1 19.099 36.947
## evno=Nunca 161 100 171.2 29.629 45.701
## evno=A veces 122 99 95.6 0.118 0.146
##
## Chisq= 50 on 2 degrees of freedom, p= 0.00000000001# Graficar las curvas de Kaplan-Meier estratificadas por treat
km_fit2 <- survfit(tiempo_evento ~ treat, data = drugs)
ggsurvplot(km_fit2, data = drugs,
risk.table = TRUE, # ver tabla de pacientes en riesgo
pval = TRUE,# agrega p-valor de logrank
conf.int = TRUE,# Intervalo de confianza 95% por default
cumcensor = TRUE,# tabla de censurados
ggtheme = theme_minimal(),
palette = "simpsons")
#logrank
logrank <- survdiff(Surv(time, relapse) ~ treat, data = drugs)
logrank## Call:
## survdiff(formula = Surv(time, relapse) ~ treat, data = drugs)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## treat=No 321 266 234 4.30 8.05
## treat=Si 307 242 274 3.68 8.05
##
## Chisq= 8 on 1 degrees of freedom, p= 0.005##residuos de martingale y deviance (para ajuste del modelo a los datos reales)
# Martingale y Deviance
ggcoxdiagnostics(cox_model3, type = "martingale", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggcoxdiagnostics(cox_model3, type = "deviance", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
dicotomizar edad en 34 años
drugs <- drugs %>%
mutate(age_q = ntile(age, 10)) # 10 grupos iguales → decilos
# Resumir y graficar
drugs %>%
dplyr::group_by(age_q) %>%
dplyr::summarise(tar = mean(relapse)) %>%
ggplot(aes(x = age_q, y = tar)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
labs(x = "Cuartiles de colesterol", y = "Proporción de eventos CHD") +
theme_minimal()
quantile(drugs$age, probs = 0.9, na.rm = TRUE)## 90%
## 40drugs$age_cat <- as.factor(ifelse(drugs$age>=40,"alto","bajo"))
table(drugs$age_cat)##
## alto bajo
## 81 547drugs$age_cat <- factor(drugs$age_cat,
levels = c("bajo","alto"))
levels(drugs$age_cat)## [1] "bajo" "alto"##modelo edad cat
cox_model_full4 <- coxph(tiempo_evento ~ age_cat+beck_cat+hercoc+evno+ndrugtx+race+treat+sex+fentanyl, data = drugs)
cox_model4 <- stepAIC(cox_model_full4, direction = "both")## Start: AIC=5806
## tiempo_evento ~ age_cat + beck_cat + hercoc + evno + ndrugtx +
## race + treat + sex + fentanyl
##
## Df AIC
## - hercoc 3 5802
## - sex 1 5804
## <none> 5806
## - race 1 5807
## - fentanyl 1 5808
## - beck_cat 1 5810
## - treat 1 5812
## - ndrugtx 1 5813
## - age_cat 1 5820
## - evno 2 5823
##
## Step: AIC=5802
## tiempo_evento ~ age_cat + beck_cat + evno + ndrugtx + race +
## treat + sex + fentanyl
##
## Df AIC
## - sex 1 5800
## <none> 5802
## - race 1 5803
## - beck_cat 1 5805
## + hercoc 3 5806
## - fentanyl 1 5807
## - treat 1 5807
## - ndrugtx 1 5809
## - age_cat 1 5818
## - evno 2 5820
##
## Step: AIC=5800
## tiempo_evento ~ age_cat + beck_cat + evno + ndrugtx + race +
## treat + fentanyl
##
## Df AIC
## <none> 5800
## - race 1 5801
## + sex 1 5802
## - beck_cat 1 5803
## + hercoc 3 5804
## - fentanyl 1 5805
## - treat 1 5805
## - ndrugtx 1 5807
## - age_cat 1 5816
## - evno 2 5820summary(cox_model4)## Call:
## coxph(formula = tiempo_evento ~ age_cat + beck_cat + evno + ndrugtx +
## race + treat + fentanyl, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.57135 0.56476 0.14365 -3.98 0.0000696 ***
## beck_catbajo -0.20662 0.81333 0.09033 -2.29 0.0222 *
## evnoNunca -0.63853 0.52807 0.13319 -4.79 0.0000016 ***
## evnoA veces -0.13575 0.87306 0.11802 -1.15 0.2501
## ndrugtx 0.02496 1.02527 0.00793 3.15 0.0016 **
## raceNegra_hisp -0.16237 0.85013 0.09539 -1.70 0.0887 .
## treatSi -0.23877 0.78760 0.08979 -2.66 0.0078 **
## fentanylSi 0.25527 1.29082 0.10088 2.53 0.0114 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.565 1.771 0.426 0.748
## beck_catbajo 0.813 1.230 0.681 0.971
## evnoNunca 0.528 1.894 0.407 0.686
## evnoA veces 0.873 1.145 0.693 1.100
## ndrugtx 1.025 0.975 1.009 1.041
## raceNegra_hisp 0.850 1.176 0.705 1.025
## treatSi 0.788 1.270 0.661 0.939
## fentanylSi 1.291 0.775 1.059 1.573
##
## Concordance= 0.627 (se = 0.013 )
## Likelihood ratio test= 101 on 8 df, p=<0.0000000000000002
## Wald test = 96.2 on 8 df, p=<0.0000000000000002
## Score (logrank) test = 100 on 8 df, p=<0.0000000000000002tab_model(cox_model4,
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 |
| age cat [alto] | 0.56 | 0.08 | 0.00 – Inf | <0.001 |
| beck cat [bajo] | 0.81 | 0.07 | 0.00 – Inf | 0.022 |
| evno [Nunca] | 0.53 | 0.07 | 0.00 – Inf | <0.001 |
| evno [A veces] | 0.87 | 0.10 | 0.00 – Inf | 0.250 |
| ndrugtx | 1.03 | 0.01 | 0.00 – Inf | 0.002 |
| race [Negra_hisp] | 0.85 | 0.08 | 0.00 – Inf | 0.089 |
| treat [Si] | 0.79 | 0.07 | 0.00 – Inf | 0.008 |
| fentanyl [Si] | 1.29 | 0.13 | 0.00 – Inf | 0.011 |
| Observations | 628 | |||
| R2 Nagelkerke | 0.149 | |||
confint(cox_model4)## 2.5 % 97.5 %
## age_catalto -0.85290 -0.2898
## beck_catbajo -0.38365 -0.0296
## evnoNunca -0.89957 -0.3775
## evnoA veces -0.36707 0.0956
## ndrugtx 0.00942 0.0405
## raceNegra_hisp -0.34932 0.0246
## treatSi -0.41474 -0.0628
## fentanylSi 0.05755 0.4530AIC(cox_model4)## [1] 5800##modelo 4 age cat
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph4 <- cox.zph(cox_model4)
test_ph4## chisq df p
## age_cat 0.273 1 0.602
## beck_cat 2.068 1 0.150
## evno 0.158 2 0.924
## ndrugtx 0.499 1 0.480
## race 2.744 1 0.098
## treat 2.955 1 0.086
## fentanyl 0.774 1 0.379
## GLOBAL 9.069 8 0.336#Gráficos por variables de residuos de Schoenfeld en el tiempo
ggcoxzph(test_ph4, ylim = c(-10, 10))
##loglog edad
# Gráficos log(-log) para la variable evno
ggsurvplot(survfit(tiempo_evento ~ age_cat, data = drugs), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
# Martingale y Deviance
ggcoxdiagnostics(cox_model4, type = "martingale", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggcoxdiagnostics(cox_model4, type = "deviance", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
##discriminación
#C de harrell
summary(cox_model4)$concordance## C se(C)
## 0.6270 0.0129#Indice D de Somer (Dxy), Metricas Q D U, slope y R2 con paquete `rms`
dd <- datadist(lung); options(datadist = "dd")
cox_metrics <- cph(tiempo_evento ~ age_cat + beck_cat + evno + ndrugtx+race+treat+fentanyl,
data = drugs, 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.2539 0.2624 0.2466 0.0158 0.2381 200
## R2 0.1485 0.1600 0.1392 0.0207 0.1278 200
## Slope 1.0000 1.0000 0.9227 0.0773 0.9227 200
## D 0.0170 0.0185 0.0158 0.0027 0.0143 200
## U -0.0003 -0.0003 0.0003 -0.0006 0.0003 200
## Q 0.0173 0.0188 0.0155 0.0033 0.0140 200
## g 0.5418 0.5643 0.5172 0.0470 0.4948 200round(cox_rms[,5],2)## Dxy R2 Slope D U Q g
## 0.24 0.13 0.92 0.01 0.00 0.01 0.49##calibración
#9. 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 = "Gráfico de Calibración",
cex.main = 1.2,
cex.axis = 0.7,
cex.lab = 0.9)
##Cox model JP sin raza
cox_modelJP<- coxph(tiempo_evento ~ age_cat + beck_cat +
evno + ndrugtx_cat2 + treat + fentanyl, data = drugs)
cox_modelJP## Call:
## coxph(formula = tiempo_evento ~ age_cat + beck_cat + evno + ndrugtx_cat2 +
## treat + fentanyl, data = drugs)
##
## coef exp(coef) se(coef) z p
## age_catalto -0.57 0.57 0.14 -4 0.000081
## beck_catbajo -0.19 0.83 0.09 -2 0.038
## evnoNunca -0.62 0.54 0.13 -5 0.000005
## evnoA veces -0.13 0.88 0.12 -1 0.284
## ndrugtx_cat2ndrugtx_menos_4 -0.30 0.74 0.09 -3 0.001
## treatSi -0.24 0.79 0.09 -3 0.007
## fentanylSi 0.28 1.33 0.10 3 0.005
##
## Likelihood ratio test=100 on 7 df, p=<0.0000000000000002
## n= 628, number of events= 508AIC(cox_modelJP)## [1] 5800confint(cox_modelJP)## 2.5 % 97.5 %
## age_catalto -0.8465 -0.2843
## beck_catbajo -0.3631 -0.0100
## evnoNunca -0.8790 -0.3520
## evnoA veces -0.3588 0.1052
## ndrugtx_cat2ndrugtx_menos_4 -0.4867 -0.1199
## treatSi -0.4169 -0.0652
## fentanylSi 0.0849 0.4803tab_model(cox_modelJP,
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 |
| age cat [alto] | 0.57 | 0.08 | 0.00 – Inf | <0.001 |
| beck cat [bajo] | 0.83 | 0.07 | 0.00 – Inf | 0.038 |
| evno [Nunca] | 0.54 | 0.07 | 0.00 – Inf | <0.001 |
| evno [A veces] | 0.88 | 0.10 | 0.00 – Inf | 0.284 |
|
ndrugtx cat2 [ndrugtx_menos_4] |
0.74 | 0.07 | 0.00 – Inf | 0.001 |
| treat [Si] | 0.79 | 0.07 | 0.00 – Inf | 0.007 |
| fentanyl [Si] | 1.33 | 0.13 | 0.00 – Inf | 0.005 |
| Observations | 628 | |||
| R2 Nagelkerke | 0.147 | |||
test_ph <- cox.zph(cox_modelJP)
test_ph## chisq df p
## age_cat 0.233 1 0.629
## beck_cat 2.418 1 0.120
## evno 0.209 2 0.901
## ndrugtx_cat2 1.177 1 0.278
## treat 2.723 1 0.099
## fentanyl 0.720 1 0.396
## GLOBAL 7.795 7 0.351summary(cox_modelJP)$concordance## C se(C)
## 0.6266 0.0131#taller 3
#Creamos evento con 3 niveles Censura, Evento de Interes, Evento Competitivo
drugs <- drugs %>%
mutate(evento = case_when(
relapse == 1 ~ 1,
death == 1 & relapse == 0 ~ 2,
death == 0 & relapse == 0 ~ 0,
TRUE~ 0
))
table(drugs$evento)##
## 0 1 2
## 74 508 46##Tiempo
#Pasamos los dias a años dividiendo por 365.25 y asignamos el tiempo hasta el evento correspondiente
#Pasamos los dias a meses dividiendo por 365.25/12=30.44 y asignamos el tiempo hasta el evento correspondiente
drugs$tiempo_yr = drugs$time/365.25
drugs$tiempo_m = drugs$time/30.4375
#level fentanilo
drugs$fentanyl <- factor(drugs$fentanyl,levels = c("Si", "No"))
levels(drugs$fentanyl)## [1] "Si" "No"##crear modelo cox tiempo m evento
#Ajustamos un modelo de Cox para evento de interes
cox_STK <- coxph(formula = Surv(tiempo_m, evento == 1) ~ age_cat + beck_cat +
evno + ndrugtx_cat2 + treat + fentanyl, data = drugs)
summary(cox_STK)## Call:
## coxph(formula = Surv(tiempo_m, evento == 1) ~ age_cat + beck_cat +
## evno + ndrugtx_cat2 + treat + fentanyl, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.5654 0.5681 0.1434 -3.94 0.0000807 ***
## beck_catbajo -0.1866 0.8298 0.0901 -2.07 0.0383 *
## evnoNunca -0.6155 0.5404 0.1344 -4.58 0.0000047 ***
## evnoA veces -0.1268 0.8809 0.1183 -1.07 0.2840
## ndrugtx_cat2ndrugtx_menos_4 -0.3033 0.7384 0.0936 -3.24 0.0012 **
## treatSi -0.2411 0.7858 0.0897 -2.69 0.0072 **
## fentanylNo -0.2826 0.7538 0.1009 -2.80 0.0051 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.568 1.76 0.429 0.753
## beck_catbajo 0.830 1.21 0.696 0.990
## evnoNunca 0.540 1.85 0.415 0.703
## evnoA veces 0.881 1.14 0.699 1.111
## ndrugtx_cat2ndrugtx_menos_4 0.738 1.35 0.615 0.887
## treatSi 0.786 1.27 0.659 0.937
## fentanylNo 0.754 1.33 0.619 0.919
##
## Concordance= 0.627 (se = 0.013 )
## Likelihood ratio test= 99.5 on 7 df, p=<0.0000000000000002
## Wald test = 93.7 on 7 df, p=<0.0000000000000002
## Score (logrank) test = 97.4 on 7 df, p=<0.0000000000000002AIC(cox_STK)## [1] 5800##modelo para cada uno. causa específica
#Opcion con objeto creado
#si evento es 1, es porque tuvo recaída
#si evento es 2, es porque se murió
tiempo_recaida <- Surv(drugs$tiempo_yr, drugs$evento==1)
tiempo_muerte <- Surv(drugs$tiempo_yr, drugs$evento==2)
cox_recaida <- coxph(tiempo_recaida ~
age_cat + beck_cat +
evno + ndrugtx_cat2 + treat + fentanyl, data = drugs)
summary(cox_recaida)## Call:
## coxph(formula = tiempo_recaida ~ age_cat + beck_cat + evno +
## ndrugtx_cat2 + treat + fentanyl, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.5654 0.5681 0.1434 -3.94 0.0000807 ***
## beck_catbajo -0.1866 0.8298 0.0901 -2.07 0.0383 *
## evnoNunca -0.6155 0.5404 0.1344 -4.58 0.0000047 ***
## evnoA veces -0.1268 0.8809 0.1183 -1.07 0.2840
## ndrugtx_cat2ndrugtx_menos_4 -0.3033 0.7384 0.0936 -3.24 0.0012 **
## treatSi -0.2411 0.7858 0.0897 -2.69 0.0072 **
## fentanylNo -0.2826 0.7538 0.1009 -2.80 0.0051 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.568 1.76 0.429 0.753
## beck_catbajo 0.830 1.21 0.696 0.990
## evnoNunca 0.540 1.85 0.415 0.703
## evnoA veces 0.881 1.14 0.699 1.111
## ndrugtx_cat2ndrugtx_menos_4 0.738 1.35 0.615 0.887
## treatSi 0.786 1.27 0.659 0.937
## fentanylNo 0.754 1.33 0.619 0.919
##
## Concordance= 0.627 (se = 0.013 )
## Likelihood ratio test= 99.5 on 7 df, p=<0.0000000000000002
## Wald test = 93.7 on 7 df, p=<0.0000000000000002
## Score (logrank) test = 97.4 on 7 df, p=<0.0000000000000002cox_muerte <- coxph(tiempo_muerte ~
age_cat + beck_cat +
evno + ndrugtx_cat2 + treat + fentanyl, data = drugs)
summary(cox_muerte)## Call:
## coxph(formula = tiempo_muerte ~ age_cat + beck_cat + evno + ndrugtx_cat2 +
## treat + fentanyl, data = drugs)
##
## n= 628, number of events= 46
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto 0.1602 1.1737 0.4066 0.39 0.6936
## beck_catbajo 0.2433 1.2755 0.3376 0.72 0.4711
## evnoNunca 1.5554 4.7370 0.5750 2.71 0.0068 **
## evnoA veces 0.5958 1.8145 0.4889 1.22 0.2230
## ndrugtx_cat2ndrugtx_menos_4 -0.7697 0.4632 0.3559 -2.16 0.0306 *
## treatSi 0.0513 1.0527 0.3076 0.17 0.8674
## fentanylNo -2.1534 0.1161 0.5303 -4.06 0.000049 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 1.174 0.852 0.5291 2.604
## beck_catbajo 1.275 0.784 0.6581 2.472
## evnoNunca 4.737 0.211 1.5349 14.619
## evnoA veces 1.815 0.551 0.6959 4.731
## ndrugtx_cat2ndrugtx_menos_4 0.463 2.159 0.2306 0.930
## treatSi 1.053 0.950 0.5761 1.924
## fentanylNo 0.116 8.614 0.0411 0.328
##
## Concordance= 0.759 (se = 0.044 )
## Likelihood ratio test= 25.2 on 7 df, p=0.0007
## Wald test = 25 on 7 df, p=0.0008
## Score (logrank) test = 28.2 on 7 df, p=0.0002#CIF
cifivhx <- cuminc(Surv(tiempo_yr, factor(evento)) ~ evno, data=drugs)
cifivhx## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## Siempre 0.500 130 0.622 0.026 0.568, 0.671
## Siempre 1.00 54 0.842 0.020 0.798, 0.876
## Siempre 1.50 20 0.894 0.017 0.856, 0.923
## Nunca 0.500 97 0.373 0.039 0.297, 0.448
## Nunca 1.00 67 0.563 0.040 0.481, 0.636
## Nunca 1.50 39 0.626 0.039 0.545, 0.697
## A veces 0.500 52 0.544 0.046 0.450, 0.629
## A veces 1.00 24 0.785 0.038 0.698, 0.850
## A veces 1.50 12 0.837 0.035 0.755, 0.893## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## Siempre 0.500 130 0.000 0.000 NA, NA
## Siempre 1.00 54 0.000 0.000 NA, NA
## Siempre 1.50 20 0.021 0.009 0.009, 0.044
## Nunca 0.500 97 0.013 0.009 0.002, 0.041
## Nunca 1.00 67 0.013 0.009 0.002, 0.041
## Nunca 1.50 39 0.087 0.023 0.048, 0.139
## A veces 0.500 52 0.008 0.008 0.001, 0.041
## A veces 1.00 24 0.008 0.008 0.001, 0.041
## A veces 1.50 12 0.028 0.017 0.007, 0.075## • Tests## outcome statistic df p.value
## 1 48.3 2.00 <0.001
## 2 18.2 2.00 <0.001levels(drugs$ivhx)## [1] "Nunca" "A veces" "Siempre"table(drugs$ivhx)##
## Nunca A veces Siempre
## 161 122 345# kaplan meier usando cif
ggcuminc(cifivhx, outcome = c(1, 2), linetype_aes = TRUE) +
scale_color_manual(
name = "Drogadicción endovenosa",
values = c("#56b", "#56B4E9","#e69f00"),
labels = c("A veces", "Nunca","Siempre")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c("Evento de interés (CHD)", "Competitivo (muerte)")
) +
add_risktable(
tables.theme = theme(
axis.text.x = element_text(size = 5),
axis.text.y = element_text(size = 5))) +
theme(
legend.position = "bottom",
legend.title = element_text(size = 9),
legend.text = element_text(size = 8),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9))## Warning in ggplot2::geom_text(size = 3.5, hjust = 0.5, vjust = 0.5, angle = 0, : Ignoring unknown parameters: `tables.theme`
## Ignoring unknown parameters: `tables.theme`
##uso cig para treat pero de drugs1
drugs$treat_num=drugs1$treat
#CIF
cifivhx <- cuminc(Surv(tiempo_yr, factor(evento)) ~ treat_num, data=drugs)
cifivhx## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## 0 0.500 121 0.612 0.027 0.556, 0.664
## 0 1.00 62 0.800 0.023 0.751, 0.840
## 0 1.50 27 0.838 0.021 0.792, 0.875
## 1 0.500 158 0.472 0.029 0.415, 0.527
## 1 1.00 83 0.718 0.026 0.663, 0.765
## 1 1.50 44 0.791 0.024 0.740, 0.833## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## 0 0.500 121 0.003 0.003 0.000, 0.016
## 0 1.00 62 0.003 0.003 0.000, 0.016
## 0 1.50 27 0.040 0.012 0.021, 0.069
## 1 0.500 158 0.007 0.005 0.001, 0.022
## 1 1.00 83 0.007 0.005 0.001, 0.022
## 1 1.50 44 0.039 0.012 0.021, 0.067## • Tests## outcome statistic df p.value
## 1 8.23 1.00 0.004
## 2 0.619 1.00 0.43# kaplan meier usando cif
ggcuminc(cifivhx, outcome = c(1, 2), linetype_aes = TRUE) +
scale_color_manual(
name = "Tratamiento previo",
values = c("#56B4E9","#e69f00"),
labels = c("No", "Si")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c("Evento de interés (CHD)", "Competitivo (muerte)")
) +
add_risktable(
tables.theme = theme(
axis.text.x = element_text(size = 5),
axis.text.y = element_text(size = 5))) +
theme(
legend.position = "bottom",
legend.title = element_text(size = 9),
legend.text = element_text(size = 8),
axis.title = element_text(size = 10),
axis.text = element_text(size = 9))## Warning in ggplot2::geom_text(size = 3.5, hjust = 0.5, vjust = 0.5, angle = 0, : Ignoring unknown parameters: `tables.theme`
## Ignoring unknown parameters: `tables.theme`
##modelo de causa específica
# Ajustar modelo con CSC
modelo_CSC <- CSC(
formula = Hist(tiempo_m, evento) ~ age_cat + beck_cat +
evno + ndrugtx_cat2 + treat + fentanyl, data = drugs,
cause = 1 # Aclara que CHD es el evento de interés
)## Warning in FUN(X[[i]], ...): Variables named time in data will be ignored.
## Warning in FUN(X[[i]], ...): Variables named time in data will be ignored.# 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(tiempo_m, evento) ~ age_cat + beck_cat + evno +
## ndrugtx_cat2 + treat + fentanyl, data = drugs, cause = 1)
##
## Right-censored response of a competing.risks model
##
## No.Observations: 628
##
## Pattern:
##
## Cause event right.censored
## 1 508 0
## 2 46 0
## unknown 0 74
##
##
## ----------> Cause: 1
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ age_cat + beck_cat +
## evno + ndrugtx_cat2 + treat + fentanyl, x = TRUE, y = TRUE)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.5654 0.5681 0.1434 -3.94 0.0000807 ***
## beck_catbajo -0.1866 0.8298 0.0901 -2.07 0.0383 *
## evnoNunca -0.6155 0.5404 0.1344 -4.58 0.0000047 ***
## evnoA veces -0.1268 0.8809 0.1183 -1.07 0.2840
## ndrugtx_cat2ndrugtx_menos_4 -0.3033 0.7384 0.0936 -3.24 0.0012 **
## treatSi -0.2411 0.7858 0.0897 -2.69 0.0072 **
## fentanylNo -0.2826 0.7538 0.1009 -2.80 0.0051 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.568 1.76 0.429 0.753
## beck_catbajo 0.830 1.21 0.696 0.990
## evnoNunca 0.540 1.85 0.415 0.703
## evnoA veces 0.881 1.14 0.699 1.111
## ndrugtx_cat2ndrugtx_menos_4 0.738 1.35 0.615 0.887
## treatSi 0.786 1.27 0.659 0.937
## fentanylNo 0.754 1.33 0.619 0.919
##
## Concordance= 0.627 (se = 0.013 )
## Likelihood ratio test= 99.5 on 7 df, p=<0.0000000000000002
## Wald test = 93.7 on 7 df, p=<0.0000000000000002
## Score (logrank) test = 97.4 on 7 df, p=<0.0000000000000002
##
##
##
## ----------> Cause: 2
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ age_cat + beck_cat +
## evno + ndrugtx_cat2 + treat + fentanyl, x = TRUE, y = TRUE)
##
## n= 628, number of events= 46
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto 0.1602 1.1737 0.4066 0.39 0.6936
## beck_catbajo 0.2433 1.2755 0.3376 0.72 0.4711
## evnoNunca 1.5554 4.7370 0.5750 2.71 0.0068 **
## evnoA veces 0.5958 1.8145 0.4889 1.22 0.2230
## ndrugtx_cat2ndrugtx_menos_4 -0.7697 0.4632 0.3559 -2.16 0.0306 *
## treatSi 0.0513 1.0527 0.3076 0.17 0.8674
## fentanylNo -2.1534 0.1161 0.5303 -4.06 0.000049 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 1.174 0.852 0.5291 2.604
## beck_catbajo 1.275 0.784 0.6581 2.472
## evnoNunca 4.737 0.211 1.5349 14.619
## evnoA veces 1.815 0.551 0.6959 4.731
## ndrugtx_cat2ndrugtx_menos_4 0.463 2.159 0.2306 0.930
## treatSi 1.053 0.950 0.5761 1.924
## fentanylNo 0.116 8.614 0.0411 0.328
##
## Concordance= 0.759 (se = 0.044 )
## Likelihood ratio test= 25.2 on 7 df, p=0.0007
## Wald test = 25 on 7 df, p=0.0008
## Score (logrank) test = 28.2 on 7 df, p=0.0002##metricas. Score
# Evaluar performance predictiva con Score
score_CSC <- Score(
object = list("CSC model" = modelo_CSC),
formula = Hist(tiempo_m, evento) ~ 1,
data = drugs,
times = c(1:24),
metrics = "auc",
summary = "IPA",
cens.model = "km",
cause = 1 #Evaluamos solo el modelo causa Evento de Interes
)## Upper limit of followup is 23.6878850102669
## Results at times higher than 23.6878850102669 are not computed.score_CSC##
## Metric AUC:
##
## Results by model:
##
## model times AUC lower upper
## <fctr> <int> <char> <char> <char>
## 1: CSC model 1 69.1 62.7 75.6
## 2: CSC model 2 66.3 61.0 71.5
## 3: CSC model 3 66.9 62.3 71.4
## 4: CSC model 4 66.7 62.5 71.0
## 5: CSC model 5 67.3 63.1 71.5
## 6: CSC model 6 65.6 61.4 69.9
## 7: CSC model 7 66.4 62.0 70.7
## 8: CSC model 8 67.1 62.7 71.6
## 9: CSC model 9 67.2 62.7 71.7
## 10: CSC model 10 68.4 63.8 73.1
## 11: CSC model 11 68.9 64.2 73.7
## 12: CSC model 12 71.7 67.0 76.4
## 13: CSC model 13 72.3 67.5 77.1
## 14: CSC model 14 73.4 68.7 78.1
## 15: CSC model 15 74.4 69.7 79.0
## 16: CSC model 16 75.2 70.6 79.8
## 17: CSC model 17 76.1 71.5 80.7
## 18: CSC model 18 76.9 72.1 81.7
## 19: CSC model 19 77.6 71.8 83.4
## 20: CSC model 20 77.2 71.0 83.5
## 21: CSC model 21 77.3 70.2 84.4
## 22: CSC model 22 75.5 67.3 83.6
## 23: CSC model 23 75.5 67.3 83.6
## model times AUC lower upper##
## NOTE: Values are multiplied by 100 and given in %.## NOTE: The higher AUC the better.##
## Metric Brier:
##
## Results by model:
##
## model times Brier lower upper IPA
## <fctr> <int> <char> <char> <char> <char>
## 1: Null model 1 9.7 7.8 11.6 0.0
## 2: Null model 2 16.2 14.3 18.0 0.0
## 3: Null model 3 21.4 20.0 22.7 0.0
## 4: Null model 4 24.1 23.4 24.9 0.0
## 5: Null model 5 24.9 24.7 25.1 0.0
## 6: Null model 6 24.8 24.5 25.2 0.0
## 7: Null model 7 24.0 23.2 24.8 0.0
## 8: Null model 8 22.5 21.4 23.7 0.0
## 9: Null model 9 21.4 20.0 22.8 0.0
## 10: Null model 10 20.1 18.6 21.7 0.0
## 11: Null model 11 19.3 17.7 21.0 0.0
## 12: Null model 12 18.3 16.5 20.0 0.0
## 13: Null model 13 17.6 15.8 19.4 0.0
## 14: Null model 14 17.0 15.2 18.9 0.0
## 15: Null model 15 16.4 14.5 18.3 0.0
## 16: Null model 16 15.8 13.9 17.7 0.0
## 17: Null model 17 15.2 13.3 17.1 0.0
## 18: Null model 18 15.1 13.2 17.0 0.0
## 19: Null model 19 14.7 12.7 16.7 0.0
## 20: Null model 20 14.7 12.7 16.7 0.0
## 21: Null model 21 14.7 12.7 16.7 0.0
## 22: Null model 22 14.2 11.8 16.5 0.0
## 23: Null model 23 14.2 11.8 16.5 0.0
## 24: CSC model 1 9.2 7.4 11.0 4.3
## 25: CSC model 2 15.3 13.5 17.0 5.5
## 26: CSC model 3 19.7 18.3 21.1 7.8
## 27: CSC model 4 22.1 21.0 23.2 8.4
## 28: CSC model 5 22.6 21.5 23.7 9.2
## 29: CSC model 6 23.0 21.8 24.2 7.2
## 30: CSC model 7 22.2 20.9 23.5 7.6
## 31: CSC model 8 20.7 19.2 22.2 8.1
## 32: CSC model 9 19.9 18.3 21.4 7.2
## 33: CSC model 10 18.5 16.8 20.1 8.3
## 34: CSC model 11 17.6 16.0 19.3 8.7
## 35: CSC model 12 16.2 14.6 17.9 11.2
## 36: CSC model 13 15.5 13.9 17.2 11.7
## 37: CSC model 14 15.0 13.3 16.6 12.1
## 38: CSC model 15 14.3 12.6 16.0 12.7
## 39: CSC model 16 13.7 12.1 15.4 13.0
## 40: CSC model 17 13.2 11.5 14.8 13.4
## 41: CSC model 18 12.9 11.3 14.6 14.3
## 42: CSC model 19 12.4 10.7 14.1 15.8
## 43: CSC model 20 12.4 10.6 14.2 15.8
## 44: CSC model 21 12.3 10.4 14.1 16.8
## 45: CSC model 22 12.1 10.0 14.1 14.6
## 46: CSC model 23 12.1 10.0 14.1 14.6
## model times Brier lower upper IPA
##
## Results of model comparisons:
##
## times model reference delta.Brier lower upper p
## <int> <fctr> <fctr> <char> <char> <char> <num>
## 1: 1 CSC model Null model -0.4 -0.7 -0.2 0.0009763
## 2: 2 CSC model Null model -0.9 -1.4 -0.4 0.0009116
## 3: 3 CSC model Null model -1.7 -2.5 -0.9 0.0000405
## 4: 4 CSC model Null model -2.0 -3.0 -1.1 0.0000489
## 5: 5 CSC model Null model -2.3 -3.4 -1.2 0.0000250
## 6: 6 CSC model Null model -1.8 -2.9 -0.6 0.0023558
## 7: 7 CSC model Null model -1.8 -3.0 -0.7 0.0021518
## 8: 8 CSC model Null model -1.8 -3.0 -0.7 0.0019097
## 9: 9 CSC model Null model -1.6 -2.7 -0.4 0.0067539
## 10: 10 CSC model Null model -1.7 -2.8 -0.6 0.0026278
## 11: 11 CSC model Null model -1.7 -2.8 -0.6 0.0020580
## 12: 12 CSC model Null model -2.0 -3.1 -1.0 0.0001101
## 13: 13 CSC model Null model -2.0 -3.1 -1.0 0.0000759
## 14: 14 CSC model Null model -2.1 -3.0 -1.1 0.0000437
## 15: 15 CSC model Null model -2.1 -3.0 -1.1 0.0000219
## 16: 16 CSC model Null model -2.1 -3.0 -1.1 0.0000168
## 17: 17 CSC model Null model -2.0 -3.0 -1.1 0.0000129
## 18: 18 CSC model Null model -2.2 -3.1 -1.2 0.0000108
## 19: 19 CSC model Null model -2.3 -3.4 -1.2 0.0000385
## 20: 20 CSC model Null model -2.3 -3.5 -1.1 0.0001192
## 21: 21 CSC model Null model -2.5 -3.8 -1.1 0.0002679
## 22: 22 CSC model Null model -2.1 -3.6 -0.6 0.0069367
## 23: 23 CSC model Null model -2.1 -3.6 -0.6 0.0069367
## times model reference delta.Brier lower upper p##
## NOTE: Values are multiplied by 100 and given in %.## NOTE: The lower Brier the better, the higher IPA the better.score_CSC$AUC$score## model times AUC se lower upper
## <fctr> <int> <num> <num> <num> <num>
## 1: CSC model 1 0.691 0.0329 0.627 0.756
## 2: CSC model 2 0.663 0.0269 0.610 0.715
## 3: CSC model 3 0.669 0.0233 0.623 0.714
## 4: CSC model 4 0.667 0.0218 0.625 0.710
## 5: CSC model 5 0.673 0.0215 0.631 0.715
## 6: CSC model 6 0.656 0.0218 0.614 0.699
## 7: CSC model 7 0.664 0.0220 0.620 0.707
## 8: CSC model 8 0.671 0.0227 0.627 0.716
## 9: CSC model 9 0.672 0.0231 0.627 0.717
## 10: CSC model 10 0.684 0.0236 0.638 0.731
## 11: CSC model 11 0.689 0.0242 0.642 0.737
## 12: CSC model 12 0.717 0.0240 0.670 0.764
## 13: CSC model 13 0.723 0.0245 0.675 0.771
## 14: CSC model 14 0.734 0.0238 0.687 0.781
## 15: CSC model 15 0.744 0.0239 0.697 0.790
## 16: CSC model 16 0.752 0.0236 0.706 0.798
## 17: CSC model 17 0.761 0.0236 0.715 0.807
## 18: CSC model 18 0.769 0.0246 0.721 0.817
## 19: CSC model 19 0.776 0.0295 0.718 0.834
## 20: CSC model 20 0.772 0.0320 0.710 0.835
## 21: CSC model 21 0.773 0.0363 0.702 0.844
## 22: CSC model 22 0.755 0.0415 0.673 0.836
## 23: CSC model 23 0.755 0.0415 0.673 0.836
## model times AUC se lower upper#fine and gray: OJO da subhzard ratio. Es la probabilidad de acumular riesgo de evento a lo largo del tiempo. No mide C de harrell el f&G
# Ajustar modelo Fine & Gray para CHD
modelo_fg <- FGR(
formula = Hist(tiempo_m, evento) ~ age_cat + beck_cat +
evno + ndrugtx_cat2 + treat + fentanyl, data = drugs,
cause = 1 # Aclara que CHD es el evento de interés
)
modelo_fg##
## Right-censored response of a competing.risks model
##
## No.Observations: 628
##
## Pattern:
##
## Cause event right.censored
## 1 508 0
## 2 46 0
## unknown 0 74
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(tiempo_m, evento) ~ age_cat + beck_cat + evno +
## ndrugtx_cat2 + treat + fentanyl, data = drugs, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## age_catalto -0.563 0.569 0.1395 -4.04 0.000054
## beck_catbajo -0.178 0.837 0.0897 -1.98 0.048000
## evnoNunca -0.650 0.522 0.1367 -4.75 0.000002
## evnoA veces -0.169 0.844 0.1207 -1.40 0.160000
## ndrugtx_cat2ndrugtx_menos_4 -0.272 0.762 0.0930 -2.93 0.003400
## treatSi -0.231 0.793 0.0888 -2.60 0.009200
## fentanylNo -0.254 0.775 0.1010 -2.52 0.012000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## age_catalto 0.569 1.76 0.433 0.748
## beck_catbajo 0.837 1.19 0.702 0.998
## evnoNunca 0.522 1.92 0.399 0.682
## evnoA veces 0.844 1.18 0.666 1.070
## ndrugtx_cat2ndrugtx_menos_4 0.762 1.31 0.635 0.914
## treatSi 0.793 1.26 0.667 0.944
## fentanylNo 0.775 1.29 0.636 0.945
##
## Num. cases = 628
## Pseudo Log-likelihood = -2901
## Pseudo likelihood ratio test = 97.9 on 7 df,
##
## Convergence: TRUE##performance predicitiva de f&g
# Evaluar performance predictiva con Score
score_fg <- Score(
object = list("FG model" = modelo_fg),
formula = Hist(tiempo_m, evento) ~ 1,
data = drugs,
times = c(1:24),
metrics = "auc",
summary = "IPA",
cens.model = "km",
cause = 1
)## Upper limit of followup is 23.6878850102669
## Results at times higher than 23.6878850102669 are not computed.score_fg$AUC$score## model times AUC se lower upper
## <fctr> <int> <num> <num> <num> <num>
## 1: FG model 1 0.691 0.0328 0.627 0.755
## 2: FG model 2 0.664 0.0269 0.611 0.716
## 3: FG model 3 0.671 0.0232 0.625 0.716
## 4: FG model 4 0.669 0.0218 0.627 0.712
## 5: FG model 5 0.675 0.0215 0.633 0.717
## 6: FG model 6 0.658 0.0218 0.615 0.700
## 7: FG model 7 0.665 0.0220 0.622 0.708
## 8: FG model 8 0.673 0.0227 0.628 0.717
## 9: FG model 9 0.673 0.0230 0.628 0.718
## 10: FG model 10 0.686 0.0236 0.639 0.732
## 11: FG model 11 0.691 0.0242 0.644 0.738
## 12: FG model 12 0.719 0.0239 0.672 0.766
## 13: FG model 13 0.726 0.0243 0.678 0.773
## 14: FG model 14 0.736 0.0237 0.689 0.783
## 15: FG model 15 0.746 0.0237 0.700 0.792
## 16: FG model 16 0.755 0.0233 0.709 0.801
## 17: FG model 17 0.763 0.0232 0.718 0.809
## 18: FG model 18 0.771 0.0243 0.724 0.819
## 19: FG model 19 0.777 0.0292 0.720 0.835
## 20: FG model 20 0.774 0.0314 0.713 0.836
## 21: FG model 21 0.776 0.0350 0.707 0.844
## 22: FG model 22 0.759 0.0386 0.684 0.835
## 23: FG model 23 0.759 0.0386 0.684 0.835
## model times AUC se lower upper##SCORE COMPARADO CON IPA
score_comp <-Score(
list("FG model" = modelo_fg, "CSC model" = modelo_CSC),
formula = Hist(tiempo_m, evento) ~ 1,
data = drugs,
times = c(1:24),
cause = 1,
metrics ="auc",
summary = "IPA"
)## Upper limit of followup is 23.6878850102669
## Results at times higher than 23.6878850102669 are not computed.Visualizar gráfico
score_comp <-Score(
list("FG model" = modelo_fg, "CSC model" = modelo_CSC),
formula = Hist(tiempo_m, evento) ~ 1,
data = drugs,
times = c(1:24),
cause = 1,
metrics ="auc",
summary = "IPA"
)## Upper limit of followup is 23.6878850102669
## Results at times higher than 23.6878850102669 are not computed.Lo ponemos visual
#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.0000
## 2 Null model 2 0.0000
## 3 Null model 3 0.0000
## 4 Null model 4 0.0000
## 5 Null model 5 0.0000
## 6 Null model 6 0.0000
## 7 Null model 7 0.0000
## 8 Null model 8 0.0000
## 9 Null model 9 0.0000
## 10 Null model 10 0.0000
## 11 Null model 11 0.0000
## 12 Null model 12 0.0000
## 13 Null model 13 0.0000
## 14 Null model 14 0.0000
## 15 Null model 15 0.0000
## 16 Null model 16 0.0000
## 17 Null model 17 0.0000
## 18 Null model 18 0.0000
## 19 Null model 19 0.0000
## 20 Null model 20 0.0000
## 21 Null model 21 0.0000
## 22 Null model 22 0.0000
## 23 Null model 23 0.0000
## 24 FG model 1 0.0431
## 25 FG model 2 0.0541
## 26 FG model 3 0.0777
## 27 FG model 4 0.0846
## 28 FG model 5 0.0922
## 29 FG model 6 0.0724
## 30 FG model 7 0.0767
## 31 FG model 8 0.0799
## 32 FG model 9 0.0717
## 33 FG model 10 0.0827
## 34 FG model 11 0.0866
## 35 FG model 12 0.1108
## 36 FG model 13 0.1155
## 37 FG model 14 0.1196
## 38 FG model 15 0.1264
## 39 FG model 16 0.1294
## 40 FG model 17 0.1336
## 41 FG model 18 0.1429
## 42 FG model 19 0.1579
## 43 FG model 20 0.1576
## 44 FG model 21 0.1665
## 45 FG model 22 0.1427
## 46 FG model 23 0.1427
## 47 CSC model 1 0.0434
## 48 CSC model 2 0.0547
## 49 CSC model 3 0.0777
## 50 CSC model 4 0.0841
## 51 CSC model 5 0.0924
## 52 CSC model 6 0.0718
## 53 CSC model 7 0.0755
## 54 CSC model 8 0.0808
## 55 CSC model 9 0.0725
## 56 CSC model 10 0.0832
## 57 CSC model 11 0.0874
## 58 CSC model 12 0.1120
## 59 CSC model 13 0.1166
## 60 CSC model 14 0.1207
## 61 CSC model 15 0.1271
## 62 CSC model 16 0.1301
## 63 CSC model 17 0.1342
## 64 CSC model 18 0.1435
## 65 CSC model 19 0.1584
## 66 CSC model 20 0.1576
## 67 CSC model 21 0.1678
## 68 CSC model 22 0.1457
## 69 CSC model 23 0.1457# Gráfico
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 = "Índice de Precisión de Predicción (IPA) en el tiempo",
x = "Años 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, ##por auc
#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 = "Discriminación de Modelos en el Tiempo",
subtitle = "Línea punteada = Azar (0.5), Línea punteada = Adecuada (0.7)",
x = "Años de seguimiento",
y = "AUC (Área 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##calibración
# Calibracion global
s <- Score(list("FGR" = modelo_fg),
formula = Hist(tiempo_m, evento) ~ 1,
data = drugs,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(s, models = "FGR" )
#por tiempo. La AUC mejora luego de los 12 meses y pasa a ser mayor a
70
# Calibracion tiempo especifico
t <- Score(list("FGR" = modelo_fg),
formula = Hist(tiempo_m, evento) ~ 1,
data = drugs,
plots = "calibration",
summary = "IPA",
times = c(1:24),
cause = 1)## Upper limit of followup is 23.6878850102669
## Results at times higher than 23.6878850102669 are not computed.plotCalibration(t, models = "FGR", times = 2 )
plotCalibration(t, models = "FGR", times = 4 )
plotCalibration(t, models = "FGR", times = 6 )
plotCalibration(t, models = "FGR", times = 8 )
plotCalibration(t, models = "FGR", times = 12 )
plotCalibration(t, models = "FGR", times = 14 )
plotCalibration(t, models = "FGR", times = 16 )
plotCalibration(t, models = "FGR", times = 18 )
plotCalibration(t, models = "FGR", times = 20 )
plotCalibration(t, models = "FGR", times = 22 )
##