TH3arreglado
2025-08-12
##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)## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!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, cuminclibrary(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"))
drugs$fentanyl <- factor(drugs$fentanyl,
levels = c("Si", "No"))
levels(drugs$fentanyl)## [1] "Si" "No"summary(drugs$fentanyl)## Si No
## 230 398#objeto de sobrevida
tiempo_evento <- Surv(drugs$time,drugs$relapse)#categorización devariablesy ordenamiento
###ivhx
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"###EDAD
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 edad", y = "Proporción de recaídas") +
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"drugs<- drugs %>% mutate(age_q= ntile(age, 5))
#como se comporta beck con el evento: beck prodría dicotomizarse en el quintilo 3
drugs<- drugs %>% mutate(beck_q= ntile(beck, 5))
drugs<- drugs %>% mutate(beck_q2= ntile(beck, 10))
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_q2) %>%
dplyr::summarise(tar = mean(relapse)) %>%
ggplot(aes(x = beck_q2, y = tar)) +
geom_col(fill = "steelblue") +
geom_text(aes(label = round(tar, 2)), vjust = -0.3) +
labs(x = "Cuartiles de beck", y = "Proporción de recaídas") +
theme_minimal()
#cuál es el quintilo 3 de beck
# Calcular el 3er quintil (40%)
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###NDRUGTX
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 recaídas") +
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")
))##Summary drugs
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 evno age_q age_cat
## Min. : 2 Si:230 Min. :0.000 Siempre:345 Min. :1 bajo:547
## 1st Qu.: 74 No:398 1st Qu.:0.000 Nunca :161 1st Qu.:2 alto: 81
## Median :159 Median :0.000 A veces:122 Median :3
## Mean :224 Mean :0.073 Mean :3
## 3rd Qu.:340 3rd Qu.:0.000 3rd Qu.:4
## Max. :721 Max. :1.000 Max. :5
## beck_q beck_q2 beck_cat ndrugtx_q ndrugtx_cat2
## Min. :1 Min. : 1.00 alto:262 Min. :1 ndrugtx_mas_4 :279
## 1st Qu.:2 1st Qu.: 3.00 bajo:366 1st Qu.:2 ndrugtx_menos_4:349
## Median :3 Median : 5.00 Median :3
## Mean :3 Mean : 5.49 Mean :3
## 3rd Qu.:4 3rd Qu.: 8.00 3rd Qu.:4
## Max. :5 Max. :10.00 Max. :5###modelo automático
cox_model_full4 <- coxph(tiempo_evento ~ age_cat+beck_cat+hercoc+evno+ndrugtx_cat2+race+treat+sex+fentanyl, data = drugs)
cox_model4 <- stepAIC(cox_model_full4, direction = "both")## Start: AIC=5804
## tiempo_evento ~ age_cat + beck_cat + hercoc + evno + ndrugtx_cat2 +
## race + treat + sex + fentanyl
##
## Df AIC
## - hercoc 3 5801
## - sex 1 5802
## <none> 5804
## - race 1 5805
## - fentanyl 1 5806
## - beck_cat 1 5807
## - treat 1 5810
## - ndrugtx_cat2 1 5813
## - age_cat 1 5817
## - evno 2 5821
##
## Step: AIC=5801
## tiempo_evento ~ age_cat + beck_cat + evno + ndrugtx_cat2 + race +
## treat + sex + fentanyl
##
## Df AIC
## - sex 1 5799
## <none> 5801
## - race 1 5801
## - beck_cat 1 5803
## + hercoc 3 5804
## - treat 1 5806
## - fentanyl 1 5807
## - ndrugtx_cat2 1 5809
## - age_cat 1 5816
## - evno 2 5818
##
## Step: AIC=5799
## tiempo_evento ~ age_cat + beck_cat + evno + ndrugtx_cat2 + race +
## treat + fentanyl
##
## Df AIC
## <none> 5799
## - race 1 5800
## + sex 1 5801
## - beck_cat 1 5802
## + hercoc 3 5802
## - treat 1 5804
## - fentanyl 1 5805
## - ndrugtx_cat2 1 5807
## - age_cat 1 5815
## - evno 2 5817summary(cox_model4)## Call:
## coxph(formula = tiempo_evento ~ age_cat + beck_cat + evno + ndrugtx_cat2 +
## race + treat + fentanyl, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.5633 0.5693 0.1433 -3.93 0.0000848 ***
## beck_catbajo -0.1902 0.8268 0.0901 -2.11 0.0349 *
## evnoNunca -0.6167 0.5397 0.1343 -4.59 0.0000044 ***
## evnoA veces -0.1236 0.8837 0.1182 -1.05 0.2958
## ndrugtx_cat2ndrugtx_menos_4 -0.2950 0.7446 0.0936 -3.15 0.0016 **
## raceNegra_hisp -0.1519 0.8591 0.0955 -1.59 0.1120
## treatSi -0.2401 0.7866 0.0898 -2.67 0.0075 **
## fentanylNo -0.2780 0.7573 0.1009 -2.75 0.0059 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.569 1.76 0.430 0.754
## beck_catbajo 0.827 1.21 0.693 0.987
## evnoNunca 0.540 1.85 0.415 0.702
## evnoA veces 0.884 1.13 0.701 1.114
## ndrugtx_cat2ndrugtx_menos_4 0.745 1.34 0.620 0.895
## raceNegra_hisp 0.859 1.16 0.712 1.036
## treatSi 0.787 1.27 0.660 0.938
## fentanylNo 0.757 1.32 0.621 0.923
##
## Concordance= 0.627 (se = 0.013 )
## Likelihood ratio test= 102 on 8 df, p=<0.0000000000000002
## Wald test = 95.8 on 8 df, p=<0.0000000000000002
## Score (logrank) test = 99.6 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.57 | 0.08 | 0.00 – Inf | <0.001 |
| beck cat [bajo] | 0.83 | 0.07 | 0.00 – Inf | 0.035 |
| evno [Nunca] | 0.54 | 0.07 | 0.00 – Inf | <0.001 |
| evno [A veces] | 0.88 | 0.10 | 0.00 – Inf | 0.296 |
|
ndrugtx cat2 [ndrugtx_menos_4] |
0.74 | 0.07 | 0.00 – Inf | 0.002 |
| race [Negra_hisp] | 0.86 | 0.08 | 0.00 – Inf | 0.112 |
| treat [Si] | 0.79 | 0.07 | 0.00 – Inf | 0.007 |
| fentanyl [No] | 0.76 | 0.08 | 0.00 – Inf | 0.006 |
| Observations | 628 | |||
| R2 Nagelkerke | 0.150 | |||
confint(cox_model4)## 2.5 % 97.5 %
## age_catalto -0.844 -0.2824
## beck_catbajo -0.367 -0.0135
## evnoNunca -0.880 -0.3535
## evnoA veces -0.355 0.1081
## ndrugtx_cat2ndrugtx_menos_4 -0.478 -0.1114
## raceNegra_hisp -0.339 0.0354
## treatSi -0.416 -0.0641
## fentanylNo -0.476 -0.0802AIC(cox_model4)## [1] 5799##Test de proporcionalidad de Schoenfeld
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph <- cox.zph(cox_model4)
test_ph## chisq df p
## age_cat 0.262 1 0.609
## beck_cat 2.301 1 0.129
## evno 0.209 2 0.901
## ndrugtx_cat2 1.254 1 0.263
## race 2.555 1 0.110
## treat 2.805 1 0.094
## fentanyl 0.826 1 0.363
## GLOBAL 10.117 8 0.257#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 de proporcionalidad para ndrugtx_cat2
ggsurvplot(survfit(tiempo_evento ~ ndrugtx_cat2, data = drugs), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
#Gráfico loglog beck
ggsurvplot(survfit(tiempo_evento ~ beck_cat, data = drugs), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
#Gráfico loglog variable age_cat
ggsurvplot(survfit(tiempo_evento ~ age_cat, 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")
#Gráfico loglog variable race
ggsurvplot(survfit(tiempo_evento ~ race, data = drugs), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
Kaplan meier
# 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# Graficar las curvas de Kaplan-Meier estratificadas por age_cat
km_fit2 <- survfit(tiempo_evento ~ age_cat, 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) ~ age_cat, data = drugs)
logrank## Call:
## survdiff(formula = Surv(time, relapse) ~ age_cat, data = drugs)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## age_cat=bajo 547 451 430.6 0.968 6.39
## age_cat=alto 81 57 77.4 5.384 6.39
##
## Chisq= 6.4 on 1 degrees of freedom, p= 0.01# Graficar las curvas de Kaplan-Meier estratificadas por ndrugtx_cat2
km_fit2 <- survfit(tiempo_evento ~ ndrugtx_cat2, 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) ~ ndrugtx_cat2, data = drugs)
logrank## Call:
## survdiff(formula = Surv(time, relapse) ~ ndrugtx_cat2, data = drugs)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## ndrugtx_cat2=ndrugtx_mas_4 279 244 193 13.32 21.8
## ndrugtx_cat2=ndrugtx_menos_4 349 264 315 8.18 21.8
##
## Chisq= 21.8 on 1 degrees of freedom, p= 0.000003# Graficar las curvas de Kaplan-Meier estratificadas por race
km_fit2 <- survfit(tiempo_evento ~ race, 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) ~ race, data = drugs)
logrank## Call:
## survdiff(formula = Surv(time, relapse) ~ race, data = drugs)
##
## N Observed Expected (O-E)^2/E (O-E)^2/V
## race=Blanca 189 164 143 3.15 4.41
## race=Negra_hisp 439 344 365 1.23 4.41
##
## Chisq= 4.4 on 1 degrees of freedom, p= 0.04##residuos de martingale y deviance (para ajuste del modelo a los datos reales)
# 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.627 0.013##índice de Sommer
#Indice D de Somer (Dxy), Metricas Q D U, slope y R2 con paquete `rms`
dd <- datadist(drugs); options(datadist = "dd")
cox_metrics <- cph(tiempo_evento ~ age_cat + beck_cat + evno + ndrugtx_cat2+treat+fentanyl+race,
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.2536 0.2633 0.2460 0.0174 0.2362 200
## R2 0.1499 0.1614 0.1410 0.0204 0.1295 200
## Slope 1.0000 1.0000 0.9315 0.0685 0.9315 200
## D 0.0172 0.0187 0.0160 0.0027 0.0145 200
## U -0.0003 -0.0003 0.0003 -0.0006 0.0003 200
## Q 0.0175 0.0191 0.0158 0.0033 0.0142 200
## g 0.5470 0.5693 0.5246 0.0447 0.5023 200round(cox_rms[,5],2)## Dxy R2 Slope D U Q g
## 0.24 0.13 0.93 0.01 0.00 0.01 0.50##Calibración ##calibración
# Calibración. Evaluación Gráfica
invisible(capture.output({cal <- calibrate(cox_metrics, method = "boot", B = 200, u = 180)}))
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)
#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##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+race, data = drugs)
summary(cox_STK)## Call:
## coxph(formula = Surv(tiempo_m, evento == 1) ~ age_cat + beck_cat +
## evno + ndrugtx_cat2 + treat + fentanyl + race, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.5633 0.5693 0.1433 -3.93 0.0000848 ***
## beck_catbajo -0.1902 0.8268 0.0901 -2.11 0.0349 *
## evnoNunca -0.6167 0.5397 0.1343 -4.59 0.0000044 ***
## evnoA veces -0.1236 0.8837 0.1182 -1.05 0.2958
## ndrugtx_cat2ndrugtx_menos_4 -0.2950 0.7446 0.0936 -3.15 0.0016 **
## treatSi -0.2401 0.7866 0.0898 -2.67 0.0075 **
## fentanylNo -0.2780 0.7573 0.1009 -2.75 0.0059 **
## raceNegra_hisp -0.1519 0.8591 0.0955 -1.59 0.1120
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.569 1.76 0.430 0.754
## beck_catbajo 0.827 1.21 0.693 0.987
## evnoNunca 0.540 1.85 0.415 0.702
## evnoA veces 0.884 1.13 0.701 1.114
## ndrugtx_cat2ndrugtx_menos_4 0.745 1.34 0.620 0.895
## treatSi 0.787 1.27 0.660 0.938
## fentanylNo 0.757 1.32 0.621 0.923
## raceNegra_hisp 0.859 1.16 0.712 1.036
##
## Concordance= 0.627 (se = 0.013 )
## Likelihood ratio test= 102 on 8 df, p=<0.0000000000000002
## Wald test = 95.8 on 8 df, p=<0.0000000000000002
## Score (logrank) test = 99.6 on 8 df, p=<0.0000000000000002AIC(cox_STK)## [1] 5799##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_m, drugs$evento==1)
tiempo_muerte <- Surv(drugs$tiempo_m, drugs$evento==2)
cox_recaida <- coxph(tiempo_recaida ~
age_cat + beck_cat +
evno + ndrugtx_cat2 + treat + race + fentanyl, data = drugs)
summary(cox_recaida)## Call:
## coxph(formula = tiempo_recaida ~ age_cat + beck_cat + evno +
## ndrugtx_cat2 + treat + race + fentanyl, data = drugs)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.5633 0.5693 0.1433 -3.93 0.0000848 ***
## beck_catbajo -0.1902 0.8268 0.0901 -2.11 0.0349 *
## evnoNunca -0.6167 0.5397 0.1343 -4.59 0.0000044 ***
## evnoA veces -0.1236 0.8837 0.1182 -1.05 0.2958
## ndrugtx_cat2ndrugtx_menos_4 -0.2950 0.7446 0.0936 -3.15 0.0016 **
## treatSi -0.2401 0.7866 0.0898 -2.67 0.0075 **
## raceNegra_hisp -0.1519 0.8591 0.0955 -1.59 0.1120
## fentanylNo -0.2780 0.7573 0.1009 -2.75 0.0059 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.569 1.76 0.430 0.754
## beck_catbajo 0.827 1.21 0.693 0.987
## evnoNunca 0.540 1.85 0.415 0.702
## evnoA veces 0.884 1.13 0.701 1.114
## ndrugtx_cat2ndrugtx_menos_4 0.745 1.34 0.620 0.895
## treatSi 0.787 1.27 0.660 0.938
## raceNegra_hisp 0.859 1.16 0.712 1.036
## fentanylNo 0.757 1.32 0.621 0.923
##
## Concordance= 0.627 (se = 0.013 )
## Likelihood ratio test= 102 on 8 df, p=<0.0000000000000002
## Wald test = 95.8 on 8 df, p=<0.0000000000000002
## Score (logrank) test = 99.6 on 8 df, p=<0.0000000000000002cox_muerte <- coxph(tiempo_muerte ~
age_cat + beck_cat +
evno + ndrugtx_cat2 + treat + race +fentanyl, data = drugs)
summary(cox_muerte)## Call:
## coxph(formula = tiempo_muerte ~ age_cat + beck_cat + evno + ndrugtx_cat2 +
## treat + race + fentanyl, data = drugs)
##
## n= 628, number of events= 46
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto 0.1610 1.1747 0.4064 0.40 0.6920
## beck_catbajo 0.2453 1.2780 0.3386 0.72 0.4688
## evnoNunca 1.5500 4.7116 0.5788 2.68 0.0074 **
## evnoA veces 0.5931 1.8095 0.4898 1.21 0.2260
## ndrugtx_cat2ndrugtx_menos_4 -0.7633 0.4661 0.3637 -2.10 0.0359 *
## treatSi 0.0535 1.0549 0.3088 0.17 0.8625
## raceNegra_hisp -0.0342 0.9663 0.4168 -0.08 0.9345
## fentanylNo -2.1579 0.1156 0.5334 -4.05 0.000052 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 1.175 0.851 0.5297 2.605
## beck_catbajo 1.278 0.782 0.6582 2.481
## evnoNunca 4.712 0.212 1.5153 14.650
## evnoA veces 1.810 0.553 0.6928 4.726
## ndrugtx_cat2ndrugtx_menos_4 0.466 2.145 0.2285 0.951
## treatSi 1.055 0.948 0.5760 1.932
## raceNegra_hisp 0.966 1.035 0.4270 2.187
## fentanylNo 0.116 8.653 0.0406 0.329
##
## Concordance= 0.759 (se = 0.044 )
## Likelihood ratio test= 25.2 on 8 df, p=0.001
## Wald test = 24.9 on 8 df, p=0.002
## Score (logrank) test = 28.2 on 8 df, p=0.0004#CIF
cifivhx <- cuminc(Surv(tiempo_m, factor(evento)) ~ treat, data=drugs)
cifivhx## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## No 5.00 138 0.558 0.028 0.502, 0.611
## No 10.0 70 0.774 0.024 0.724, 0.817
## No 15.0 56 0.819 0.022 0.771, 0.857
## No 20.0 8 0.846 0.021 0.798, 0.882
## Si 5.00 184 0.386 0.028 0.332, 0.441
## Si 10.0 99 0.664 0.027 0.608, 0.715
## Si 15.0 68 0.767 0.024 0.715, 0.811
## Si 20.0 9 0.795 0.023 0.745, 0.837## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## No 5.00 138 0.003 0.003 0.000, 0.016
## No 10.0 70 0.003 0.003 0.000, 0.016
## No 15.0 56 0.003 0.003 0.000, 0.016
## No 20.0 8 0.087 0.020 0.054, 0.131
## Si 5.00 184 0.007 0.005 0.001, 0.022
## Si 10.0 99 0.007 0.005 0.001, 0.022
## Si 15.0 68 0.007 0.005 0.001, 0.022
## Si 20.0 9 0.097 0.020 0.062, 0.140## • Tests## outcome statistic df p.value
## 1 8.23 1.00 0.004
## 2 0.619 1.00 0.43Gráfico
levels(drugs$treat)## [1] "No" "Si"table(drugs$treat)##
## No Si
## 321 307# kaplan meier usando cif
ggcuminc(cifivhx, outcome = c(1, 2), linetype_aes = TRUE) +
scale_color_manual(
name = "Tratamiento",
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`
table(drugs$treat)##
## No Si
## 321 307##modelo de causa específica
# Ajustar modelo con CSC
modelo_CSC <- CSC(
formula = Hist(tiempo_m, evento) ~ age_cat + beck_cat +
evno + ndrugtx_cat2 + race +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 + race + 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 + race + treat + fentanyl, x = TRUE,
## y = TRUE)
##
## n= 628, number of events= 508
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.5633 0.5693 0.1433 -3.93 0.0000848 ***
## beck_catbajo -0.1902 0.8268 0.0901 -2.11 0.0349 *
## evnoNunca -0.6167 0.5397 0.1343 -4.59 0.0000044 ***
## evnoA veces -0.1236 0.8837 0.1182 -1.05 0.2958
## ndrugtx_cat2ndrugtx_menos_4 -0.2950 0.7446 0.0936 -3.15 0.0016 **
## raceNegra_hisp -0.1519 0.8591 0.0955 -1.59 0.1120
## treatSi -0.2401 0.7866 0.0898 -2.67 0.0075 **
## fentanylNo -0.2780 0.7573 0.1009 -2.75 0.0059 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.569 1.76 0.430 0.754
## beck_catbajo 0.827 1.21 0.693 0.987
## evnoNunca 0.540 1.85 0.415 0.702
## evnoA veces 0.884 1.13 0.701 1.114
## ndrugtx_cat2ndrugtx_menos_4 0.745 1.34 0.620 0.895
## raceNegra_hisp 0.859 1.16 0.712 1.036
## treatSi 0.787 1.27 0.660 0.938
## fentanylNo 0.757 1.32 0.621 0.923
##
## Concordance= 0.627 (se = 0.013 )
## Likelihood ratio test= 102 on 8 df, p=<0.0000000000000002
## Wald test = 95.8 on 8 df, p=<0.0000000000000002
## Score (logrank) test = 99.6 on 8 df, p=<0.0000000000000002
##
##
##
## ----------> Cause: 2
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ age_cat + beck_cat +
## evno + ndrugtx_cat2 + race + treat + fentanyl, x = TRUE,
## y = TRUE)
##
## n= 628, number of events= 46
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto 0.1610 1.1747 0.4064 0.40 0.6920
## beck_catbajo 0.2453 1.2780 0.3386 0.72 0.4688
## evnoNunca 1.5500 4.7116 0.5788 2.68 0.0074 **
## evnoA veces 0.5931 1.8095 0.4898 1.21 0.2260
## ndrugtx_cat2ndrugtx_menos_4 -0.7633 0.4661 0.3637 -2.10 0.0359 *
## raceNegra_hisp -0.0342 0.9663 0.4168 -0.08 0.9345
## treatSi 0.0535 1.0549 0.3088 0.17 0.8625
## fentanylNo -2.1579 0.1156 0.5334 -4.05 0.000052 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 1.175 0.851 0.5297 2.605
## beck_catbajo 1.278 0.782 0.6582 2.481
## evnoNunca 4.712 0.212 1.5153 14.650
## evnoA veces 1.810 0.553 0.6928 4.726
## ndrugtx_cat2ndrugtx_menos_4 0.466 2.145 0.2285 0.951
## raceNegra_hisp 0.966 1.035 0.4270 2.187
## treatSi 1.055 0.948 0.5760 1.932
## fentanylNo 0.116 8.653 0.0406 0.329
##
## Concordance= 0.759 (se = 0.044 )
## Likelihood ratio test= 25.2 on 8 df, p=0.001
## Wald test = 24.9 on 8 df, p=0.002
## Score (logrank) test = 28.2 on 8 df, p=0.0004##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.7 63.3 76.0
## 2: CSC model 2 66.5 61.3 71.8
## 3: CSC model 3 66.8 62.2 71.3
## 4: CSC model 4 66.3 62.0 70.6
## 5: CSC model 5 67.1 62.9 71.3
## 6: CSC model 6 65.7 61.4 70.0
## 7: CSC model 7 66.6 62.3 70.9
## 8: CSC model 8 67.4 63.0 71.9
## 9: CSC model 9 67.3 62.8 71.8
## 10: CSC model 10 68.9 64.3 73.4
## 11: CSC model 11 69.6 64.9 74.2
## 12: CSC model 12 72.4 67.8 77.0
## 13: CSC model 13 73.0 68.3 77.6
## 14: CSC model 14 74.0 69.4 78.5
## 15: CSC model 15 74.9 70.3 79.5
## 16: CSC model 16 75.8 71.2 80.4
## 17: CSC model 17 76.8 72.3 81.3
## 18: CSC model 18 77.8 73.1 82.5
## 19: CSC model 19 78.5 72.9 84.2
## 20: CSC model 20 78.1 72.0 84.2
## 21: CSC model 21 78.8 72.2 85.4
## 22: CSC model 22 77.0 69.5 84.5
## 23: CSC model 23 77.0 69.5 84.5
## 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.4
## 25: CSC model 2 15.2 13.5 17.0 5.7
## 26: CSC model 3 19.7 18.3 21.2 7.6
## 27: CSC model 4 22.3 21.1 23.4 7.8
## 28: CSC model 5 22.7 21.6 23.8 8.8
## 29: CSC model 6 23.1 21.9 24.3 7.0
## 30: CSC model 7 22.2 20.8 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.4 16.8 20.0 8.6
## 34: CSC model 11 17.6 15.9 19.2 9.2
## 35: CSC model 12 16.1 14.5 17.8 11.7
## 36: CSC model 13 15.4 13.8 17.1 12.1
## 37: CSC model 14 14.9 13.3 16.6 12.5
## 38: CSC model 15 14.2 12.6 15.9 13.2
## 39: CSC model 16 13.6 12.0 15.3 13.7
## 40: CSC model 17 13.0 11.4 14.7 14.2
## 41: CSC model 18 12.8 11.2 14.4 15.2
## 42: CSC model 19 12.3 10.5 14.0 16.7
## 43: CSC model 20 12.3 10.5 14.1 16.5
## 44: CSC model 21 12.1 10.3 13.9 17.8
## 45: CSC model 22 12.0 10.0 14.0 15.4
## 46: CSC model 23 12.0 10.0 14.0 15.4
## 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.00082292
## 2: 2 CSC model Null model -0.9 -1.5 -0.4 0.00062828
## 3: 3 CSC model Null model -1.6 -2.4 -0.8 0.00007010
## 4: 4 CSC model Null model -1.9 -2.9 -0.9 0.00018389
## 5: 5 CSC model Null model -2.2 -3.3 -1.1 0.00007656
## 6: 6 CSC model Null model -1.7 -2.9 -0.6 0.00324463
## 7: 7 CSC model Null model -1.8 -3.0 -0.7 0.00214657
## 8: 8 CSC model Null model -1.8 -3.0 -0.7 0.00204783
## 9: 9 CSC model Null model -1.6 -2.7 -0.4 0.00735528
## 10: 10 CSC model Null model -1.7 -2.8 -0.6 0.00215395
## 11: 11 CSC model Null model -1.8 -2.8 -0.7 0.00132079
## 12: 12 CSC model Null model -2.1 -3.2 -1.1 0.00005941
## 13: 13 CSC model Null model -2.1 -3.2 -1.1 0.00004179
## 14: 14 CSC model Null model -2.1 -3.1 -1.1 0.00002590
## 15: 15 CSC model Null model -2.2 -3.1 -1.2 0.00001269
## 16: 16 CSC model Null model -2.2 -3.1 -1.2 0.00000774
## 17: 17 CSC model Null model -2.2 -3.1 -1.2 0.00000532
## 18: 18 CSC model Null model -2.3 -3.3 -1.3 0.00000395
## 19: 19 CSC model Null model -2.5 -3.6 -1.3 0.00002045
## 20: 20 CSC model Null model -2.4 -3.6 -1.2 0.00008578
## 21: 21 CSC model Null model -2.6 -4.0 -1.3 0.00011218
## 22: 22 CSC model Null model -2.2 -3.7 -0.7 0.00468228
## 23: 23 CSC model Null model -2.2 -3.7 -0.7 0.00468228
## 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.697 0.0323 0.633 0.760
## 2: CSC model 2 0.665 0.0267 0.613 0.718
## 3: CSC model 3 0.668 0.0232 0.622 0.713
## 4: CSC model 4 0.663 0.0219 0.620 0.706
## 5: CSC model 5 0.671 0.0215 0.629 0.713
## 6: CSC model 6 0.657 0.0218 0.614 0.700
## 7: CSC model 7 0.666 0.0219 0.623 0.709
## 8: CSC model 8 0.674 0.0226 0.630 0.719
## 9: CSC model 9 0.673 0.0230 0.628 0.718
## 10: CSC model 10 0.689 0.0234 0.643 0.734
## 11: CSC model 11 0.696 0.0238 0.649 0.742
## 12: CSC model 12 0.724 0.0234 0.678 0.770
## 13: CSC model 13 0.730 0.0239 0.683 0.776
## 14: CSC model 14 0.740 0.0233 0.694 0.785
## 15: CSC model 15 0.749 0.0234 0.703 0.795
## 16: CSC model 16 0.758 0.0232 0.712 0.804
## 17: CSC model 17 0.768 0.0231 0.723 0.813
## 18: CSC model 18 0.778 0.0238 0.731 0.825
## 19: CSC model 19 0.785 0.0287 0.729 0.842
## 20: CSC model 20 0.781 0.0311 0.720 0.842
## 21: CSC model 21 0.788 0.0338 0.722 0.854
## 22: CSC model 22 0.770 0.0385 0.695 0.845
## 23: CSC model 23 0.770 0.0385 0.695 0.845
## 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+race, 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 + race, data = drugs, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## age_catalto -0.563 0.570 0.1376 -4.09 0.0000430
## beck_catbajo -0.182 0.833 0.0891 -2.05 0.0410000
## evnoNunca -0.650 0.522 0.1372 -4.74 0.0000021
## evnoA veces -0.164 0.849 0.1195 -1.37 0.1700000
## ndrugtx_cat2ndrugtx_menos_4 -0.265 0.768 0.0931 -2.84 0.0045000
## treatSi -0.231 0.794 0.0884 -2.61 0.0090000
## fentanylNo -0.251 0.778 0.1005 -2.49 0.0130000
## raceNegra_hisp -0.159 0.853 0.0907 -1.76 0.0790000
##
## exp(coef) exp(-coef) 2.5% 97.5%
## age_catalto 0.570 1.76 0.435 0.746
## beck_catbajo 0.833 1.20 0.700 0.992
## evnoNunca 0.522 1.92 0.399 0.683
## evnoA veces 0.849 1.18 0.672 1.073
## ndrugtx_cat2ndrugtx_menos_4 0.768 1.30 0.640 0.921
## treatSi 0.794 1.26 0.668 0.944
## fentanylNo 0.778 1.28 0.639 0.948
## raceNegra_hisp 0.853 1.17 0.714 1.019
##
## Num. cases = 628
## Pseudo Log-likelihood = -2900
## Pseudo likelihood ratio test = 101 on 8 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.699 0.0320 0.636 0.761
## 2: FG model 2 0.667 0.0266 0.614 0.719
## 3: FG model 3 0.669 0.0231 0.624 0.714
## 4: FG model 4 0.664 0.0218 0.622 0.707
## 5: FG model 5 0.671 0.0215 0.629 0.713
## 6: FG model 6 0.657 0.0218 0.614 0.700
## 7: FG model 7 0.666 0.0220 0.623 0.709
## 8: FG model 8 0.673 0.0227 0.629 0.718
## 9: FG model 9 0.673 0.0230 0.628 0.718
## 10: FG model 10 0.688 0.0234 0.642 0.734
## 11: FG model 11 0.696 0.0238 0.649 0.743
## 12: FG model 12 0.724 0.0234 0.678 0.770
## 13: FG model 13 0.731 0.0239 0.684 0.777
## 14: FG model 14 0.741 0.0233 0.695 0.786
## 15: FG model 15 0.750 0.0233 0.705 0.796
## 16: FG model 16 0.760 0.0231 0.715 0.805
## 17: FG model 17 0.769 0.0229 0.725 0.814
## 18: FG model 18 0.780 0.0236 0.733 0.826
## 19: FG model 19 0.786 0.0284 0.731 0.842
## 20: FG model 20 0.783 0.0306 0.723 0.842
## 21: FG model 21 0.789 0.0330 0.724 0.854
## 22: FG model 22 0.771 0.0365 0.700 0.843
## 23: FG model 23 0.771 0.0365 0.700 0.843
## 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.0441
## 25 FG model 2 0.0566
## 26 FG model 3 0.0758
## 27 FG model 4 0.0785
## 28 FG model 5 0.0874
## 29 FG model 6 0.0708
## 30 FG model 7 0.0775
## 31 FG model 8 0.0804
## 32 FG model 9 0.0718
## 33 FG model 10 0.0852
## 34 FG model 11 0.0911
## 35 FG model 12 0.1159
## 36 FG model 13 0.1207
## 37 FG model 14 0.1243
## 38 FG model 15 0.1313
## 39 FG model 16 0.1364
## 40 FG model 17 0.1414
## 41 FG model 18 0.1521
## 42 FG model 19 0.1671
## 43 FG model 20 0.1648
## 44 FG model 21 0.1779
## 45 FG model 22 0.1511
## 46 FG model 23 0.1511
## 47 CSC model 1 0.0444
## 48 CSC model 2 0.0572
## 49 CSC model 3 0.0759
## 50 CSC model 4 0.0785
## 51 CSC model 5 0.0879
## 52 CSC model 6 0.0703
## 53 CSC model 7 0.0763
## 54 CSC model 8 0.0812
## 55 CSC model 9 0.0725
## 56 CSC model 10 0.0856
## 57 CSC model 11 0.0916
## 58 CSC model 12 0.1167
## 59 CSC model 13 0.1214
## 60 CSC model 14 0.1249
## 61 CSC model 15 0.1316
## 62 CSC model 16 0.1366
## 63 CSC model 17 0.1415
## 64 CSC model 18 0.1521
## 65 CSC model 19 0.1671
## 66 CSC model 20 0.1646
## 67 CSC model 21 0.1784
## 68 CSC model 22 0.1537
## 69 CSC model 23 0.1537# 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 )
Hago subgrupo de los pacientes tratados
tratados <- drugs[drugs$treat == "Si", ]
tiempo_evtratados <- Surv(tratados$time,tratados$relapse)
#modelo automático
cox_model_tratados <- coxph(tiempo_evtratados ~ age_cat+beck_cat+hercoc+evno+ndrugtx_cat2+race+sex+fentanyl, data = tratados)
cox_model_tratados <- stepAIC(cox_model_tratados, direction = "both")## Start: AIC=2445
## tiempo_evtratados ~ age_cat + beck_cat + hercoc + evno + ndrugtx_cat2 +
## race + sex + fentanyl
##
## Df AIC
## - hercoc 3 2439
## - sex 1 2443
## - ndrugtx_cat2 1 2443
## <none> 2445
## - fentanyl 1 2445
## - race 1 2445
## - beck_cat 1 2448
## - evno 2 2450
## - age_cat 1 2450
##
## Step: AIC=2439
## tiempo_evtratados ~ age_cat + beck_cat + evno + ndrugtx_cat2 +
## race + sex + fentanyl
##
## Df AIC
## - sex 1 2437
## - ndrugtx_cat2 1 2438
## <none> 2439
## - race 1 2440
## - fentanyl 1 2441
## - beck_cat 1 2443
## + hercoc 3 2445
## - evno 2 2446
## - age_cat 1 2446
##
## Step: AIC=2437
## tiempo_evtratados ~ age_cat + beck_cat + evno + ndrugtx_cat2 +
## race + fentanyl
##
## Df AIC
## - ndrugtx_cat2 1 2436
## <none> 2437
## - race 1 2439
## - fentanyl 1 2439
## + sex 1 2439
## - beck_cat 1 2441
## + hercoc 3 2443
## - evno 2 2444
## - age_cat 1 2444
##
## Step: AIC=2436
## tiempo_evtratados ~ age_cat + beck_cat + evno + race + fentanyl
##
## Df AIC
## <none> 2436
## - race 1 2437
## + ndrugtx_cat2 1 2437
## - fentanyl 1 2438
## + sex 1 2438
## - beck_cat 1 2440
## + hercoc 3 2442
## - age_cat 1 2442
## - evno 2 2445summary(cox_model_tratados)## Call:
## coxph(formula = tiempo_evtratados ~ age_cat + beck_cat + evno +
## race + fentanyl, data = tratados)
##
## n= 307, number of events= 242
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.560 0.571 0.208 -2.69 0.00718 **
## beck_catbajo -0.315 0.730 0.132 -2.38 0.01728 *
## evnoNunca -0.661 0.516 0.188 -3.52 0.00044 ***
## evnoA veces -0.165 0.848 0.165 -1.00 0.31697
## raceNegra_hisp -0.261 0.770 0.141 -1.85 0.06399 .
## fentanylNo -0.287 0.751 0.151 -1.89 0.05838 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.571 1.75 0.380 0.859
## beck_catbajo 0.730 1.37 0.564 0.946
## evnoNunca 0.516 1.94 0.357 0.746
## evnoA veces 0.848 1.18 0.614 1.171
## raceNegra_hisp 0.770 1.30 0.584 1.015
## fentanylNo 0.751 1.33 0.558 1.010
##
## Concordance= 0.619 (se = 0.018 )
## Likelihood ratio test= 44.5 on 6 df, p=0.00000006
## Wald test = 42.8 on 6 df, p=0.0000001
## Score (logrank) test = 44.4 on 6 df, p=0.00000006tab_model(cox_model_tratados,
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.12 | 0.00 – Inf | 0.007 |
| beck cat [bajo] | 0.73 | 0.10 | 0.00 – Inf | 0.017 |
| evno [Nunca] | 0.52 | 0.10 | 0.00 – Inf | <0.001 |
| evno [A veces] | 0.85 | 0.14 | 0.00 – Inf | 0.317 |
| race [Negra_hisp] | 0.77 | 0.11 | 0.00 – Inf | 0.064 |
| fentanyl [No] | 0.75 | 0.11 | 0.00 – Inf | 0.058 |
| Observations | 307 | |||
| R2 Nagelkerke | 0.135 | |||
confint(cox_model_tratados)## 2.5 % 97.5 %
## age_catalto -0.968 -0.1517
## beck_catbajo -0.574 -0.0556
## evnoNunca -1.030 -0.2927
## evnoA veces -0.488 0.1582
## raceNegra_hisp -0.538 0.0152
## fentanylNo -0.584 0.0102AIC(cox_model_tratados)## [1] 2436cox_model_tratados1 <- coxph(tiempo_evtratados ~ age_cat+beck_cat+evno+fentanyl, data = tratados)
cox_model_tratados1## Call:
## coxph(formula = tiempo_evtratados ~ age_cat + beck_cat + evno +
## fentanyl, data = tratados)
##
## coef exp(coef) se(coef) z p
## age_catalto -0.6 0.6 0.2 -2.7 0.007
## beck_catbajo -0.3 0.7 0.1 -2.2 0.026
## evnoNunca -0.7 0.5 0.2 -3.5 0.0005
## evnoA veces -0.2 0.9 0.2 -0.9 0.344
## fentanylNo -0.3 0.7 0.1 -2.3 0.023
##
## Likelihood ratio test=41 on 5 df, p=0.00000009
## n= 307, number of events= 242summary(cox_model_tratados1)## Call:
## coxph(formula = tiempo_evtratados ~ age_cat + beck_cat + evno +
## fentanyl, data = tratados)
##
## n= 307, number of events= 242
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.560 0.571 0.208 -2.69 0.00713 **
## beck_catbajo -0.292 0.747 0.132 -2.22 0.02634 *
## evnoNunca -0.656 0.519 0.187 -3.50 0.00047 ***
## evnoA veces -0.156 0.856 0.165 -0.95 0.34444
## fentanylNo -0.337 0.714 0.148 -2.27 0.02324 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.571 1.75 0.380 0.859
## beck_catbajo 0.747 1.34 0.577 0.966
## evnoNunca 0.519 1.93 0.359 0.750
## evnoA veces 0.856 1.17 0.620 1.182
## fentanylNo 0.714 1.40 0.534 0.955
##
## Concordance= 0.613 (se = 0.018 )
## Likelihood ratio test= 41.2 on 5 df, p=0.00000009
## Wald test = 39.2 on 5 df, p=0.0000002
## Score (logrank) test = 40.9 on 5 df, p=0.0000001tab_model(cox_model_tratados1,
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.12 | 0.00 – Inf | 0.007 |
| beck cat [bajo] | 0.75 | 0.10 | 0.00 – Inf | 0.026 |
| evno [Nunca] | 0.52 | 0.10 | 0.00 – Inf | <0.001 |
| evno [A veces] | 0.86 | 0.14 | 0.00 – Inf | 0.344 |
| fentanyl [No] | 0.71 | 0.11 | 0.00 – Inf | 0.023 |
| Observations | 307 | |||
| R2 Nagelkerke | 0.126 | |||
confint(cox_model_tratados1)## 2.5 % 97.5 %
## age_catalto -0.968 -0.1521
## beck_catbajo -0.550 -0.0344
## evnoNunca -1.023 -0.2881
## evnoA veces -0.478 0.1669
## fentanylNo -0.628 -0.0460AIC(cox_model_tratados1)## [1] 2437##Test de proporcionalidad de Schoenfeld
#Test de proporcionalidad de riesgos de Schoenfeld
test_ph <- cox.zph(cox_model_tratados1)
test_ph## chisq df p
## age_cat 0.576 1 0.45
## beck_cat 0.341 1 0.56
## evno 0.313 2 0.86
## fentanyl 1.479 1 0.22
## GLOBAL 3.310 5 0.65#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_evtratados ~ evno, data = tratados), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
# Gráficos log(-log) para la variable beck
ggsurvplot(survfit(tiempo_evtratados ~ beck_cat, data = tratados), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
#Gráfico loglog variable age_cat
ggsurvplot(survfit(tiempo_evtratados ~ age_cat, data = tratados), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
#Gráfico loglog variable fentanyl
ggsurvplot(survfit(tiempo_evtratados ~ fentanyl, data = tratados), fun = "cloglog",
xlab = "Tiempo (dias)",
ylab = "Log(-Log(S(t)))",
palette = "simpsons")
# Martingale y Deviance
ggcoxdiagnostics(cox_model_tratados1, type = "martingale", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
ggcoxdiagnostics(cox_model_tratados1, type = "deviance", linear.predictions = FALSE)## `geom_smooth()` using formula = 'y ~ x'
#C de harrell
summary(cox_model_tratados1)$concordance## C se(C)
## 0.6128 0.0184#Indice D de Somer (Dxy), Metricas Q D U, slope y R2 con paquete `rms`
dd <- datadist(tratados); options(datadist = "dd")
cox_metrics <- cph(tiempo_evtratados ~ age_cat + beck_cat + evno + fentanyl,
data = tratados, 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.2257 0.2373 0.2210 0.0163 0.2094 200
## R2 0.1257 0.1379 0.1145 0.0234 0.1022 200
## Slope 1.0000 1.0000 0.9176 0.0824 0.9176 200
## D 0.0163 0.0181 0.0147 0.0034 0.0129 200
## U -0.0008 -0.0008 0.0006 -0.0014 0.0006 200
## Q 0.0171 0.0189 0.0142 0.0048 0.0123 200
## g 0.5116 0.5301 0.4781 0.0520 0.4596 200round(cox_rms[,5],2)## Dxy R2 Slope D U Q g
## 0.21 0.10 0.92 0.01 0.00 0.01 0.46# Calibración. Evaluación Gráfica
invisible(capture.output({cal <- calibrate(cox_metrics, method = "boot", B = 200, u = 180)}))
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)
#Creamos evento con 3 niveles Censura, Evento de Interes, Evento Competitivo
tratados <- tratados %>%
mutate(eventotrat = case_when(
relapse == 1 ~ 1,
death == 1 & relapse == 0 ~ 2,
death == 0 & relapse == 0 ~ 0,
TRUE~ 0
))#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
tratados$tiempo_yr = tratados$time/365.25
tratados$tiempo_m = tratados$time/30.4375#Ajustamos un modelo de Cox para evento de interes
cox_STK1 <- coxph(formula = Surv(tiempo_m, eventotrat == 1) ~ age_cat + beck_cat +
evno + fentanyl, data = tratados)
summary(cox_STK1)## Call:
## coxph(formula = Surv(tiempo_m, eventotrat == 1) ~ age_cat + beck_cat +
## evno + fentanyl, data = tratados)
##
## n= 307, number of events= 242
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.560 0.571 0.208 -2.69 0.00713 **
## beck_catbajo -0.292 0.747 0.132 -2.22 0.02634 *
## evnoNunca -0.656 0.519 0.187 -3.50 0.00047 ***
## evnoA veces -0.156 0.856 0.165 -0.95 0.34444
## fentanylNo -0.337 0.714 0.148 -2.27 0.02324 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.571 1.75 0.380 0.859
## beck_catbajo 0.747 1.34 0.577 0.966
## evnoNunca 0.519 1.93 0.359 0.750
## evnoA veces 0.856 1.17 0.620 1.182
## fentanylNo 0.714 1.40 0.534 0.955
##
## Concordance= 0.613 (se = 0.018 )
## Likelihood ratio test= 41.2 on 5 df, p=0.00000009
## Wald test = 39.2 on 5 df, p=0.0000002
## Score (logrank) test = 40.9 on 5 df, p=0.0000001AIC(cox_STK1)## [1] 2437#Opcion con objeto creado
#si evento es 1, es porque tuvo recaída
#si evento es 2, es porque se murió
tiempo_recaida1 <- Surv(tratados$tiempo_m, tratados$eventotrat==1)
tiempo_muerte1 <- Surv(tratados$tiempo_m, tratados$eventotrat==2)
cox_recaida1 <- coxph(tiempo_recaida1 ~
age_cat + beck_cat +
evno + fentanyl, data = tratados)
summary(cox_recaida1)## Call:
## coxph(formula = tiempo_recaida1 ~ age_cat + beck_cat + evno +
## fentanyl, data = tratados)
##
## n= 307, number of events= 242
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.560 0.571 0.208 -2.69 0.00713 **
## beck_catbajo -0.292 0.747 0.132 -2.22 0.02634 *
## evnoNunca -0.656 0.519 0.187 -3.50 0.00047 ***
## evnoA veces -0.156 0.856 0.165 -0.95 0.34444
## fentanylNo -0.337 0.714 0.148 -2.27 0.02324 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.571 1.75 0.380 0.859
## beck_catbajo 0.747 1.34 0.577 0.966
## evnoNunca 0.519 1.93 0.359 0.750
## evnoA veces 0.856 1.17 0.620 1.182
## fentanylNo 0.714 1.40 0.534 0.955
##
## Concordance= 0.613 (se = 0.018 )
## Likelihood ratio test= 41.2 on 5 df, p=0.00000009
## Wald test = 39.2 on 5 df, p=0.0000002
## Score (logrank) test = 40.9 on 5 df, p=0.0000001cox_muerte1 <- coxph(tiempo_muerte1 ~
age_cat + beck_cat +
evno + fentanyl, data =tratados)
summary(cox_muerte1)## Call:
## coxph(formula = tiempo_muerte1 ~ age_cat + beck_cat + evno +
## fentanyl, data = tratados)
##
## n= 307, number of events= 26
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto 0.755 2.127 0.468 1.61 0.1067
## beck_catbajo 0.657 1.928 0.531 1.24 0.2159
## evnoNunca 1.323 3.755 0.628 2.11 0.0351 *
## evnoA veces 0.513 1.670 0.609 0.84 0.4001
## fentanylNo -2.044 0.130 0.645 -3.17 0.0015 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 2.13 0.470 0.8502 5.320
## beck_catbajo 1.93 0.519 0.6816 5.456
## evnoNunca 3.76 0.266 1.0971 12.855
## evnoA veces 1.67 0.599 0.5059 5.512
## fentanylNo 0.13 7.721 0.0366 0.459
##
## Concordance= 0.716 (se = 0.065 )
## Likelihood ratio test= 15.4 on 5 df, p=0.009
## Wald test = 15.1 on 5 df, p=0.01
## Score (logrank) test = 17.4 on 5 df, p=0.004Edad CIF
#CIF
cifivhx1 <- cuminc(Surv(tiempo_m, factor(eventotrat)) ~ age_cat, data=tratados)
cifivhx1## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## bajo 5.00 157 0.400 0.030 0.340, 0.458
## bajo 10.0 82 0.682 0.029 0.622, 0.735
## bajo 15.0 55 0.786 0.025 0.730, 0.831
## bajo 20.0 7 0.814 0.024 0.760, 0.856
## alto 5.00 27 0.297 0.073 0.164, 0.442
## alto 10.0 17 0.548 0.080 0.380, 0.689
## alto 15.0 13 0.649 0.077 0.476, 0.777
## alto 20.0 2 0.674 0.076 0.501, 0.798## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## bajo 5.00 157 0.004 0.004 0.000, 0.020
## bajo 10.0 82 0.004 0.004 0.000, 0.020
## bajo 15.0 55 0.004 0.004 0.000, 0.020
## bajo 20.0 7 0.077 0.018 0.047, 0.118
## alto 5.00 27 0.024 0.024 0.002, 0.112
## alto 10.0 17 0.024 0.024 0.002, 0.112
## alto 15.0 13 0.024 0.024 0.002, 0.112
## alto 20.0 2 0.208 0.085 0.073, 0.390## • Tests## outcome statistic df p.value
## 1 3.46 1.00 0.063
## 2 5.78 1.00 0.016levels(tratados$age_cat)## [1] "bajo" "alto"table(tratados$age_cat)##
## bajo alto
## 266 41# kaplan meier usando cif
ggcuminc(cifivhx1, outcome = c(1, 2), linetype_aes = TRUE) +
scale_color_manual(
name = "Edad menor a 40",
values = c("#56B4E9","#e69f00"),
labels = c("Bajo", "Alto")
) +
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`
Fentanyl CIF
#CIF
cifivhx1 <- cuminc(Surv(tiempo_m, factor(eventotrat)) ~ fentanyl, data=tratados)
cifivhx1## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## Si 5.00 57 0.486 0.048 0.390, 0.576
## Si 10.0 22 0.797 0.039 0.707, 0.861
## Si 15.0 10 0.908 0.029 0.833, 0.950
## Si 20.0 0 0.926 0.026 0.853, 0.964
## No 5.00 127 0.329 0.034 0.264, 0.396
## No 10.0 77 0.589 0.036 0.516, 0.655
## No 15.0 58 0.688 0.034 0.617, 0.749
## No 20.0 9 0.719 0.033 0.650, 0.777## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## Si 5.00 57 0.000 0.000 NA, NA
## Si 10.0 22 0.000 0.000 NA, NA
## Si 15.0 10 0.000 0.000 NA, NA
## Si 20.0 0 0.074 0.026 0.033, 0.137
## No 5.00 127 0.010 0.007 0.002, 0.034
## No 10.0 77 0.010 0.007 0.002, 0.034
## No 15.0 58 0.010 0.007 0.002, 0.034
## No 20.0 9 0.103 0.027 0.058, 0.164## • Tests## outcome statistic df p.value
## 1 22.4 1.00 <0.001
## 2 0.049 1.00 0.82# kaplan meier usando cif
ggcuminc(cifivhx1, outcome = c(1, 2), linetype_aes = TRUE) +
scale_color_manual(
name = "Fentanilo",
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`
Beck
#CIF
cifivhx1 <- cuminc(Surv(tiempo_m, factor(eventotrat)) ~ beck_cat, data=tratados)
cifivhx1## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## alto 5.00 63 0.470 0.046 0.378, 0.557
## alto 10.0 34 0.708 0.042 0.616, 0.781
## alto 15.0 21 0.819 0.036 0.736, 0.879
## alto 20.0 3 0.845 0.034 0.764, 0.900
## bajo 5.00 121 0.333 0.035 0.266, 0.401
## bajo 10.0 65 0.637 0.036 0.562, 0.702
## bajo 15.0 47 0.734 0.033 0.664, 0.792
## bajo 20.0 6 0.764 0.032 0.694, 0.819## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## alto 5.00 63 0.000 0.000 NA, NA
## alto 10.0 34 0.000 0.000 NA, NA
## alto 15.0 21 0.000 0.000 NA, NA
## alto 20.0 3 0.040 0.021 0.013, 0.096
## bajo 5.00 121 0.011 0.008 0.002, 0.035
## bajo 10.0 65 0.011 0.008 0.002, 0.035
## bajo 15.0 47 0.011 0.008 0.002, 0.035
## bajo 20.0 6 0.129 0.028 0.080, 0.190## • Tests## outcome statistic df p.value
## 1 5.26 1.00 0.022
## 2 3.21 1.00 0.073table(tratados$beck_cat)##
## alto bajo
## 120 187# kaplan meier usando cif
ggcuminc(cifivhx1, outcome = c(1, 2), linetype_aes = TRUE) +
scale_color_manual(
name = "Beck",
values = c("#56B4E9","#e69f00"),
labels = c("Mayor a 20", "Menor a 20")
) +
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`
EVNO
#CIF
cifivhx1 <- cuminc(Surv(tiempo_m, factor(eventotrat)) ~ evno, data=tratados)
cifivhx1## ## ── cuminc() ────────────────────────────────────────────────────────────────────## • Failure type "1"## strata time n.risk estimate std.error 95% CI
## Siempre 5.00 85 0.459 0.040 0.379, 0.535
## Siempre 10.0 38 0.754 0.035 0.678, 0.815
## Siempre 15.0 21 0.864 0.028 0.798, 0.910
## Siempre 20.0 2 0.880 0.027 0.815, 0.924
## Nunca 5.00 62 0.228 0.046 0.144, 0.323
## Nunca 10.0 39 0.505 0.055 0.393, 0.608
## Nunca 15.0 34 0.566 0.055 0.451, 0.665
## Nunca 20.0 5 0.614 0.054 0.499, 0.710
## A veces 5.00 37 0.415 0.062 0.294, 0.532
## A veces 10.0 22 0.652 0.060 0.520, 0.756
## A veces 15.0 13 0.795 0.052 0.670, 0.876
## A veces 20.0 2 0.826 0.049 0.705, 0.901## • Failure type "2"## strata time n.risk estimate std.error 95% CI
## Siempre 5.00 85 0.000 0.000 NA, NA
## Siempre 10.0 38 0.000 0.000 NA, NA
## Siempre 15.0 21 0.000 0.000 NA, NA
## Siempre 20.0 2 0.041 0.018 0.015, 0.088
## Nunca 5.00 62 0.024 0.017 0.005, 0.076
## Nunca 10.0 39 0.024 0.017 0.005, 0.076
## Nunca 15.0 34 0.024 0.017 0.005, 0.076
## Nunca 20.0 5 0.197 0.051 0.108, 0.305
## A veces 5.00 37 0.000 0.000 NA, NA
## A veces 10.0 22 0.000 0.000 NA, NA
## A veces 15.0 13 0.000 0.000 NA, NA
## A veces 20.0 2 0.090 0.040 0.032, 0.187## • Tests## outcome statistic df p.value
## 1 22.3 2.00 <0.001
## 2 9.71 2.00 0.008table(tratados$evno)##
## Siempre Nunca A veces
## 157 85 65# kaplan meier usando cif
ggcuminc(cifivhx1, outcome = c(1, 2), linetype_aes = TRUE) +
scale_color_manual(
name = "Drogadicción endovenosa",
values = c("#56B4E9","#e69f00","#56b"),
labels = c("Siempre", "Nunca","A veces")
) +
scale_linetype_manual(
name = "Tipo de evento",
values = c("solid", "dashed"),
labels = c("Evento de interés (Recaída)", "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 438
# Ajustar modelo con CSC
modelo_CSC1 <- CSC(
formula = Hist(tiempo_m, eventotrat) ~ age_cat + beck_cat +
evno + fentanyl, data = tratados,
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_CSC1## CSC(formula = Hist(tiempo_m, eventotrat) ~ age_cat + beck_cat +
## evno + fentanyl, data = tratados, cause = 1)
##
## Right-censored response of a competing.risks model
##
## No.Observations: 307
##
## Pattern:
##
## Cause event right.censored
## 1 242 0
## 2 26 0
## unknown 0 39
##
##
## ----------> Cause: 1
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ age_cat + beck_cat +
## evno + fentanyl, x = TRUE, y = TRUE)
##
## n= 307, number of events= 242
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto -0.560 0.571 0.208 -2.69 0.00713 **
## beck_catbajo -0.292 0.747 0.132 -2.22 0.02634 *
## evnoNunca -0.656 0.519 0.187 -3.50 0.00047 ***
## evnoA veces -0.156 0.856 0.165 -0.95 0.34444
## fentanylNo -0.337 0.714 0.148 -2.27 0.02324 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 0.571 1.75 0.380 0.859
## beck_catbajo 0.747 1.34 0.577 0.966
## evnoNunca 0.519 1.93 0.359 0.750
## evnoA veces 0.856 1.17 0.620 1.182
## fentanylNo 0.714 1.40 0.534 0.955
##
## Concordance= 0.613 (se = 0.018 )
## Likelihood ratio test= 41.2 on 5 df, p=0.00000009
## Wald test = 39.2 on 5 df, p=0.0000002
## Score (logrank) test = 40.9 on 5 df, p=0.0000001
##
##
##
## ----------> Cause: 2
##
## Call:
## coxph(formula = survival::Surv(time, status) ~ age_cat + beck_cat +
## evno + fentanyl, x = TRUE, y = TRUE)
##
## n= 307, number of events= 26
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_catalto 0.755 2.127 0.468 1.61 0.1067
## beck_catbajo 0.657 1.928 0.531 1.24 0.2159
## evnoNunca 1.323 3.755 0.628 2.11 0.0351 *
## evnoA veces 0.513 1.670 0.609 0.84 0.4001
## fentanylNo -2.044 0.130 0.645 -3.17 0.0015 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_catalto 2.13 0.470 0.8502 5.320
## beck_catbajo 1.93 0.519 0.6816 5.456
## evnoNunca 3.76 0.266 1.0971 12.855
## evnoA veces 1.67 0.599 0.5059 5.512
## fentanylNo 0.13 7.721 0.0366 0.459
##
## Concordance= 0.716 (se = 0.065 )
## Likelihood ratio test= 15.4 on 5 df, p=0.009
## Wald test = 15.1 on 5 df, p=0.01
## Score (logrank) test = 17.4 on 5 df, p=0.004# Evaluar performance predictiva con Score
score_CSC <- Score(
object = list("CSC model" = modelo_CSC1),
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 68.7 62.5 74.9
## 2: CSC model 2 64.0 58.9 69.2
## 3: CSC model 3 66.1 61.6 70.7
## 4: CSC model 4 64.1 59.8 68.5
## 5: CSC model 5 64.5 60.1 68.8
## 6: CSC model 6 64.0 59.7 68.3
## 7: CSC model 7 65.5 61.2 69.8
## 8: CSC model 8 65.2 60.8 69.7
## 9: CSC model 9 65.7 61.2 70.2
## 10: CSC model 10 66.7 62.1 71.3
## 11: CSC model 11 66.8 62.0 71.5
## 12: CSC model 12 68.7 63.9 73.5
## 13: CSC model 13 69.7 64.8 74.5
## 14: CSC model 14 70.7 65.9 75.5
## 15: CSC model 15 72.2 67.5 76.9
## 16: CSC model 16 72.9 68.1 77.6
## 17: CSC model 17 73.6 68.8 78.4
## 18: CSC model 18 74.0 68.8 79.3
## 19: CSC model 19 75.5 69.2 81.8
## 20: CSC model 20 75.6 68.8 82.5
## 21: CSC model 21 77.1 69.7 84.4
## 22: CSC model 22 75.5 66.2 84.9
## 23: CSC model 23 75.5 66.2 84.9
## 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.3 7.4 11.2 3.7
## 25: CSC model 2 15.5 13.6 17.4 3.8
## 26: CSC model 3 20.0 18.4 21.6 6.3
## 27: CSC model 4 23.1 21.7 24.5 4.4
## 28: CSC model 5 23.8 22.6 25.1 4.4
## 29: CSC model 6 23.7 22.6 24.8 4.6
## 30: CSC model 7 22.5 21.4 23.7 6.1
## 31: CSC model 8 21.4 20.2 22.7 5.0
## 32: CSC model 9 20.3 18.9 21.7 5.2
## 33: CSC model 10 19.1 17.7 20.5 5.2
## 34: CSC model 11 18.3 16.8 19.8 5.5
## 35: CSC model 12 16.9 15.4 18.5 7.4
## 36: CSC model 13 16.2 14.6 17.7 8.0
## 37: CSC model 14 15.5 14.0 17.1 8.8
## 38: CSC model 15 14.7 13.1 16.3 10.3
## 39: CSC model 16 14.1 12.5 15.7 10.6
## 40: CSC model 17 13.5 11.9 15.1 11.1
## 41: CSC model 18 13.3 11.7 14.9 12.0
## 42: CSC model 19 12.7 11.0 14.4 13.7
## 43: CSC model 20 12.6 10.9 14.4 14.1
## 44: CSC model 21 12.4 10.6 14.1 16.0
## 45: CSC model 22 12.0 9.9 14.1 15.5
## 46: CSC model 23 12.0 9.9 14.1 15.5
## 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.5 -0.2 0.0001226
## 2: 2 CSC model Null model -0.6 -1.0 -0.2 0.0036177
## 3: 3 CSC model Null model -1.4 -2.0 -0.7 0.0000973
## 4: 4 CSC model Null model -1.1 -2.0 -0.1 0.0273646
## 5: 5 CSC model Null model -1.1 -2.2 0.0 0.0546931
## 6: 6 CSC model Null model -1.1 -2.3 0.1 0.0622190
## 7: 7 CSC model Null model -1.5 -2.7 -0.3 0.0163992
## 8: 8 CSC model Null model -1.1 -2.4 0.1 0.0737826
## 9: 9 CSC model Null model -1.1 -2.3 0.1 0.0629019
## 10: 10 CSC model Null model -1.0 -2.2 0.1 0.0832020
## 11: 11 CSC model Null model -1.1 -2.2 0.1 0.0705656
## 12: 12 CSC model Null model -1.4 -2.5 -0.2 0.0166992
## 13: 13 CSC model Null model -1.4 -2.5 -0.3 0.0111214
## 14: 14 CSC model Null model -1.5 -2.5 -0.4 0.0054932
## 15: 15 CSC model Null model -1.7 -2.7 -0.7 0.0007534
## 16: 16 CSC model Null model -1.7 -2.6 -0.7 0.0006374
## 17: 17 CSC model Null model -1.7 -2.6 -0.7 0.0004861
## 18: 18 CSC model Null model -1.8 -2.8 -0.8 0.0003616
## 19: 19 CSC model Null model -2.0 -3.2 -0.9 0.0005861
## 20: 20 CSC model Null model -2.1 -3.3 -0.9 0.0008705
## 21: 21 CSC model Null model -2.4 -3.7 -1.0 0.0005198
## 22: 22 CSC model Null model -2.2 -3.7 -0.7 0.0047164
## 23: 23 CSC model Null model -2.2 -3.7 -0.7 0.0047164
## 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.687 0.0316 0.625 0.749
## 2: CSC model 2 0.640 0.0264 0.589 0.692
## 3: CSC model 3 0.661 0.0233 0.616 0.707
## 4: CSC model 4 0.641 0.0222 0.598 0.685
## 5: CSC model 5 0.645 0.0220 0.601 0.688
## 6: CSC model 6 0.640 0.0220 0.597 0.683
## 7: CSC model 7 0.655 0.0221 0.612 0.698
## 8: CSC model 8 0.652 0.0228 0.608 0.697
## 9: CSC model 9 0.657 0.0230 0.612 0.702
## 10: CSC model 10 0.667 0.0235 0.621 0.713
## 11: CSC model 11 0.668 0.0242 0.620 0.715
## 12: CSC model 12 0.687 0.0246 0.639 0.735
## 13: CSC model 13 0.697 0.0248 0.648 0.745
## 14: CSC model 14 0.707 0.0245 0.659 0.755
## 15: CSC model 15 0.722 0.0241 0.675 0.769
## 16: CSC model 16 0.729 0.0242 0.681 0.776
## 17: CSC model 17 0.736 0.0246 0.688 0.784
## 18: CSC model 18 0.740 0.0269 0.688 0.793
## 19: CSC model 19 0.755 0.0322 0.692 0.818
## 20: CSC model 20 0.756 0.0349 0.688 0.825
## 21: CSC model 21 0.771 0.0374 0.697 0.844
## 22: CSC model 22 0.755 0.0475 0.662 0.849
## 23: CSC model 23 0.755 0.0475 0.662 0.849
## model times AUC se lower upper##fyg 462
# Ajustar modelo Fine & Gray para CHD
modelo_fg1 <- FGR(
formula = Hist(tiempo_m, eventotrat) ~ age_cat + beck_cat +
evno + fentanyl, data = tratados,
cause = 1 # Aclara que CHD es el evento de interés
)
modelo_fg1##
## Right-censored response of a competing.risks model
##
## No.Observations: 307
##
## Pattern:
##
## Cause event right.censored
## 1 242 0
## 2 26 0
## unknown 0 39
##
##
## Fine-Gray model: analysis of cause 1
##
## Competing Risks Regression
##
## Call:
## FGR(formula = Hist(tiempo_m, eventotrat) ~ age_cat + beck_cat +
## evno + fentanyl, data = tratados, cause = 1)
##
## coef exp(coef) se(coef) z p-value
## age_catalto -0.575 0.563 0.214 -2.691 0.00710
## beck_catbajo -0.298 0.742 0.130 -2.298 0.02200
## evnoNunca -0.680 0.507 0.189 -3.594 0.00033
## evnoA veces -0.161 0.851 0.169 -0.956 0.34000
## fentanylNo -0.329 0.720 0.147 -2.237 0.02500
##
## exp(coef) exp(-coef) 2.5% 97.5%
## age_catalto 0.563 1.78 0.370 0.855
## beck_catbajo 0.742 1.35 0.575 0.957
## evnoNunca 0.507 1.97 0.350 0.734
## evnoA veces 0.851 1.18 0.611 1.185
## fentanylNo 0.720 1.39 0.540 0.960
##
## Num. cases = 307
## Pseudo Log-likelihood = -1217
## Pseudo likelihood ratio test = 43.7 on 5 df,
##
## Convergence: TRUEPredicción Fine and Gray
# Evaluar performance predictiva con Score
score_fg1 <- Score(
object = list("FG model" = modelo_fg1),
formula = Hist(tiempo_m, eventotrat) ~ 1,
data = tratados,
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_fg1$AUC$score## model times AUC se lower upper
## <fctr> <int> <num> <num> <num> <num>
## 1: FG model 1 0.675 0.0514 0.574 0.776
## 2: FG model 2 0.614 0.0420 0.532 0.696
## 3: FG model 3 0.644 0.0359 0.574 0.715
## 4: FG model 4 0.651 0.0328 0.587 0.715
## 5: FG model 5 0.647 0.0323 0.584 0.711
## 6: FG model 6 0.639 0.0315 0.577 0.701
## 7: FG model 7 0.666 0.0308 0.606 0.727
## 8: FG model 8 0.663 0.0310 0.603 0.724
## 9: FG model 9 0.663 0.0312 0.602 0.724
## 10: FG model 10 0.670 0.0310 0.609 0.731
## 11: FG model 11 0.673 0.0316 0.611 0.735
## 12: FG model 12 0.695 0.0320 0.632 0.758
## 13: FG model 13 0.704 0.0319 0.642 0.767
## 14: FG model 14 0.709 0.0315 0.648 0.771
## 15: FG model 15 0.739 0.0298 0.680 0.797
## 16: FG model 16 0.740 0.0297 0.682 0.799
## 17: FG model 17 0.737 0.0305 0.678 0.797
## 18: FG model 18 0.733 0.0325 0.670 0.797
## 19: FG model 19 0.763 0.0367 0.692 0.835
## 20: FG model 20 0.766 0.0406 0.686 0.845
## 21: FG model 21 0.798 0.0440 0.711 0.884
## 22: FG model 22 0.774 0.0509 0.674 0.874
## 23: FG model 23 0.774 0.0509 0.674 0.874
## model times AUC se lower upperScore comparado con IPA 500
score_comp <-Score(
list("FG model" = modelo_fg1, "CSC model" = modelo_CSC1),
formula = Hist(tiempo_m, eventotrat) ~ 1,
data = tratados,
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.#Extraigo las IPAs por tiempo y modelo. Grafico
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.0318
## 25 FG model 2 0.0261
## 26 FG model 3 0.0533
## 27 FG model 4 0.0666
## 28 FG model 5 0.0675
## 29 FG model 6 0.0592
## 30 FG model 7 0.0863
## 31 FG model 8 0.0805
## 32 FG model 9 0.0700
## 33 FG model 10 0.0747
## 34 FG model 11 0.0757
## 35 FG model 12 0.0980
## 36 FG model 13 0.1051
## 37 FG model 14 0.1043
## 38 FG model 15 0.1261
## 39 FG model 16 0.1229
## 40 FG model 17 0.1178
## 41 FG model 18 0.1128
## 42 FG model 19 0.1511
## 43 FG model 20 0.1566
## 44 FG model 21 0.1993
## 45 FG model 22 0.1632
## 46 FG model 23 0.1632
## 47 CSC model 1 0.0316
## 48 CSC model 2 0.0258
## 49 CSC model 3 0.0532
## 50 CSC model 4 0.0664
## 51 CSC model 5 0.0676
## 52 CSC model 6 0.0597
## 53 CSC model 7 0.0866
## 54 CSC model 8 0.0812
## 55 CSC model 9 0.0711
## 56 CSC model 10 0.0757
## 57 CSC model 11 0.0766
## 58 CSC model 12 0.0984
## 59 CSC model 13 0.1056
## 60 CSC model 14 0.1047
## 61 CSC model 15 0.1261
## 62 CSC model 16 0.1230
## 63 CSC model 17 0.1183
## 64 CSC model 18 0.1139
## 65 CSC model 19 0.1524
## 66 CSC model 20 0.1578
## 67 CSC model 21 0.2004
## 68 CSC model 22 0.1723
## 69 CSC model 23 0.1723# 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 = "Meses de seguimiento", y = "IPA", color = "Modelo"
) +
theme_minimal() +
theme(legend.position = "bottom", plot.title.position = "plot")
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 = "Meses 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# Calibracion global
s <- Score(list("FGR" = modelo_fg1),
formula = Hist(tiempo_m, eventotrat) ~ 1,
data = tratados,
plots = "calibration",
summary = "IPA",
cause = 1)
plotCalibration(s, models = "FGR" )
# Calibracion tiempo especifico
t <- Score(list("FGR" = modelo_fg1),
formula = Hist(tiempo_m, eventotrat) ~ 1,
data = tratados,
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 )
##tabla 1 para toda la cohorte
library(tableone)
drugs$Recaida <- factor(drugs$relapse, levels = c(0, 1), labels = c( "No", "Si"))
table(drugs$Recaida)##
## No Si
## 120 508#Las variables van directamente sin el nombre del dataframe
catvars = c( "age_cat","beck_cat","ndrugtx_cat2","evno","sex","race", "Recaida")
#edad y enojo son continuas
vars = c("age_cat","beck_cat","ndrugtx_cat2","evno","sex","race","age","beck","ndrugtx","time","fentanyl","treat","hercoc")
tabla_1 <- CreateTableOne(vars = vars, strata = "Recaida", factorVars = catvars, data = drugs)
print(tabla_1)## Stratified by Recaida
## No Si p
## n 120 508
## age_cat = alto (%) 24 (20.0) 57 (11.2) 0.015
## beck_cat = bajo (%) 79 (65.8) 287 (56.5) 0.078
## ndrugtx_cat2 = ndrugtx_menos_4 (%) 85 (70.8) 264 (52.0) <0.001
## evno (%) <0.001
## Siempre 36 (30.0) 309 (60.8)
## Nunca 61 (50.8) 100 (19.7)
## A veces 23 (19.2) 99 (19.5)
## sex = Hombre (%) 80 (66.7) 364 (71.7) 0.333
## race = Negra_hisp (%) 95 (79.2) 344 (67.7) 0.019
## age (mean (SD)) 33.10 (6.72) 32.12 (5.97) 0.117
## beck (mean (SD)) 16.49 (9.94) 17.96 (9.34) 0.126
## ndrugtx (mean (SD)) 3.01 (3.10) 4.92 (5.78) <0.001
## time (mean (SD)) 519.62 (157.21) 153.92 (121.91) <0.001
## fentanyl = No (%) 104 (86.7) 294 (57.9) <0.001
## treat = Si (%) 65 (54.2) 242 (47.6) 0.236
## hercoc (%) 0.136
## Her-ox 23 (19.2) 100 (19.7)
## Her 15 (12.5) 103 (20.3)
## Ox 43 (35.8) 139 (27.4)
## Coc 39 (32.5) 166 (32.7)
## Stratified by Recaida
## test
## n
## age_cat = alto (%)
## beck_cat = bajo (%)
## ndrugtx_cat2 = ndrugtx_menos_4 (%)
## evno (%)
## Siempre
## Nunca
## A veces
## sex = Hombre (%)
## race = Negra_hisp (%)
## age (mean (SD))
## beck (mean (SD))
## ndrugtx (mean (SD))
## time (mean (SD))
## fentanyl = No (%)
## treat = Si (%)
## hercoc (%)
## Her-ox
## Her
## Ox
## Cockableone(tabla_1) #se visualiza mejor| No | Si | p | test | |
|---|---|---|---|---|
| n | 120 | 508 | ||
| age_cat = alto (%) | 24 (20.0) | 57 (11.2) | 0.015 | |
| beck_cat = bajo (%) | 79 (65.8) | 287 (56.5) | 0.078 | |
| ndrugtx_cat2 = ndrugtx_menos_4 (%) | 85 (70.8) | 264 (52.0) | <0.001 | |
| evno (%) | <0.001 | |||
| Siempre | 36 (30.0) | 309 (60.8) | ||
| Nunca | 61 (50.8) | 100 (19.7) | ||
| A veces | 23 (19.2) | 99 (19.5) | ||
| sex = Hombre (%) | 80 (66.7) | 364 (71.7) | 0.333 | |
| race = Negra_hisp (%) | 95 (79.2) | 344 (67.7) | 0.019 | |
| age (mean (SD)) | 33.10 (6.72) | 32.12 (5.97) | 0.117 | |
| beck (mean (SD)) | 16.49 (9.94) | 17.96 (9.34) | 0.126 | |
| ndrugtx (mean (SD)) | 3.01 (3.10) | 4.92 (5.78) | <0.001 | |
| time (mean (SD)) | 519.62 (157.21) | 153.92 (121.91) | <0.001 | |
| fentanyl = No (%) | 104 (86.7) | 294 (57.9) | <0.001 | |
| treat = Si (%) | 65 (54.2) | 242 (47.6) | 0.236 | |
| hercoc (%) | 0.136 | |||
| Her-ox | 23 (19.2) | 100 (19.7) | ||
| Her | 15 (12.5) | 103 (20.3) | ||
| Ox | 43 (35.8) | 139 (27.4) | ||
| Coc | 39 (32.5) | 166 (32.7) |
#Predictivo puntual Recaída
# Crear paciente "nuevo" sin con tto
newdata1 <- data.frame(
age_cat = factor("alto", levels = levels(drugs$age_cat)),
beck_cat = factor("bajo", levels = levels(drugs$beck_cat)),
treat = factor("Si", levels = levels(drugs$treat)),
fentanyl = factor("No", levels = levels(drugs$fentanyl)),
ndrugtx_cat2=factor("ndrugtx_menos_4", levels=levels(drugs$ndrugtx_cat2)),
evno=factor("Nunca", levels=levels(drugs$evno)),
race=factor("Negra_hisp", levels=levels(drugs$race))
)
# Predecir probabilidad acumulada de evento a 12 y 18 meses
predictRisk(modelo_fg, newdata = newdata1, times = c(12, 18))## [,1] [,2]
## [1,] 0.305 0.359# Crear paciente "nuevo" sin tto
newdata <- data.frame(
age_cat = factor("alto", levels = levels(drugs$age_cat)),
beck_cat = factor("bajo", levels = levels(drugs$beck_cat)),
treat = factor("No", levels = levels(drugs$treat)),
fentanyl = factor("No", levels = levels(drugs$fentanyl)),
ndrugtx_cat2=factor("ndrugtx_menos_4", levels=levels(drugs$ndrugtx_cat2)),
evno=factor("Nunca", levels=levels(drugs$evno)),
race=factor("Negra_hisp", levels=levels(drugs$race))
)
# Predecir probabilidad acumulada de evento a 12 y 22 meses
predictRisk(modelo_fg, newdata = newdata, times = c(12, 18))## [,1] [,2]
## [1,] 0.368 0.429Predicción puntualen Pacientes que recibieron tratamiento
newdata3 <- data.frame(
age_cat = factor("alto", levels = levels(tratados$age_cat)),
beck_cat = factor("bajo", levels = levels(tratados$beck_cat)),
fentanyl = factor("No", levels = levels(tratados$fentanyl)),
evno=factor("Nunca", levels=levels(tratados$evno))
)
# Predecir probabilidad acumulada de evento a 12 y 22 meses
predictRisk(modelo_fg1, newdata = newdata3, times = c(12, 18))## [,1] [,2]
## [1,] 0.322 0.392