library(tidyverse)
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.3.6 ✔ purrr 0.3.4
## ✔ tibble 3.1.7 ✔ dplyr 1.0.9
## ✔ tidyr 1.2.0 ✔ stringr 1.4.0
## ✔ readr 2.1.2 ✔ forcats 0.5.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(gtsummary)
library(ggpubr)
library(survival)
library(survminer)
##
## Attaching package: 'survminer'
##
## The following object is masked from 'package:survival':
##
## myeloma
library(janitor)
##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.test
library(readxl)
dbmax <- read_excel("C:/Users/fidel/OneDrive - CINVESTAV/Tesis Max/dbmax.xlsx")
dbmax
## # A tibble: 87 × 23
## edad ECOG sintomas metas…¹ contr…² ECOG2 sinto…³ ester…⁴ contr…⁵ ECOG3
## <dbl> <dbl> <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <dbl>
## 1 39 1 sin sintomas no si 1 mas de… si si 1
## 2 NA 1 sin sintomas si si 1 cefalea no si 1
## 3 NA 1 sin sintomas si si 1 mas de… si si 1
## 4 59 1 deficit motor no no 1 sin si… no no 1
## 5 NA 1 sin sintomas no no 1 sin si… no no 1
## 6 57 1 sin sintomas no no 4 mas de… si no 4
## 7 45 1 otros no no 1 sin si… no no 1
## 8 NA 1 sin sintomas no no 4 defici… si no 4
## 9 NA 1 sin sintomas no no 1 sin si… no si 0
## 10 NA 1 sin sintomas no no 1 sin si… no si 1
## # … with 77 more rows, 13 more variables: sintomas3 <chr>, esteroides3 <chr>,
## # rescate <chr>, seg6meses <chr>, seguimiento6m <dbl>, seg12meses <chr>,
## # seguimiento12m <dbl>, res6meses <chr>, respuesta <dbl>, meses <dbl>,
## # PTV <dbl>, GTV <dbl>, dosis <dbl>, and abbreviated variable names
## # ¹`metastasis extra`, ²`control sistemico`, ³sintomas2, ⁴esteroides,
## # ⁵controlsistemico2
Recodificación de las variables
dbmax <- dbmax %>% mutate(respuesta6m = case_when(respuesta == 1 ~ 1,
respuesta == 2 ~ 1,
respuesta == 3 ~ 0,
respuesta == 4 ~ 0,
respuesta == 5 ~ 0,
respuesta == 0 ~ 1))
dbmax %>% select(respuesta,respuesta6m)
## # A tibble: 87 × 2
## respuesta respuesta6m
## <dbl> <dbl>
## 1 3 0
## 2 3 0
## 3 3 0
## 4 2 1
## 5 2 1
## 6 2 1
## 7 2 1
## 8 2 1
## 9 2 1
## 10 2 1
## # … with 77 more rows
dbmax <- dbmax %>% mutate(r6m = case_when(respuesta6m == 1 ~ "Progresión",
respuesta6m == 0 ~ "Respuesta"))
dbmax %>% select(respuesta6m, r6m)
## # A tibble: 87 × 2
## respuesta6m r6m
## <dbl> <chr>
## 1 0 Respuesta
## 2 0 Respuesta
## 3 0 Respuesta
## 4 1 Progresión
## 5 1 Progresión
## 6 1 Progresión
## 7 1 Progresión
## 8 1 Progresión
## 9 1 Progresión
## 10 1 Progresión
## # … with 77 more rows
Tabla EDA
dbmax %>% select(edad,seg6meses,r6m, PTV,GTV, dosis) %>% tbl_summary(by=r6m) %>% add_p()
## Warning for variable 'edad':
## simpleWarning in wilcox.test.default(x = DATA[[1L]], y = DATA[[2L]], ...): cannot compute exact p-value with ties
## Warning: The `fmt_missing()` function is deprecated and will soon be removed
## * Use the `sub_missing()` function instead
| Characteristic | Progresión, N = 721 | Respuesta, N = 151 | p-value2 |
|---|---|---|---|
| edad | 50 (44, 58) | 56 (44, 57) | >0.9 |
| Unknown | 45 | 6 | |
| seg6meses | 0.3 | ||
| mas de 6 meses | 52 (72%) | 13 (87%) | |
| menor o igual a 6 meses | 20 (28%) | 2 (13%) | |
| PTV | 2,418 (472, 8,728) | 6,368 (2,290, 9,179) | 0.075 |
| GTV | 1,741 (422, 5,267) | 3,777 (1,824, 8,497) | 0.080 |
| dosis | 0.2 | ||
| 15 | 9 (12%) | 1 (6.7%) | |
| 18 | 13 (18%) | 0 (0%) | |
| 20 | 1 (1.4%) | 2 (13%) | |
| 21 | 2 (2.8%) | 0 (0%) | |
| 22 | 27 (38%) | 7 (47%) | |
| 24 | 8 (11%) | 3 (20%) | |
| 25 | 1 (1.4%) | 0 (0%) | |
| 27 | 3 (4.2%) | 1 (6.7%) | |
| 30 | 8 (11%) | 1 (6.7%) | |
| 1 Median (IQR); n (%) | |||
| 2 Wilcoxon rank sum test; Fisher's exact test | |||
survfit(Surv(meses, respuesta6m) ~ 1, data = dbmax)
## Call: survfit(formula = Surv(meses, respuesta6m) ~ 1, data = dbmax)
##
## n events median 0.95LCL 0.95UCL
## [1,] 87 72 14 14 21
KM plot
ggsurvplot(
survfit(Surv(meses, respuesta6m) ~ 1, data = dbmax),
xlab = "Meses",
ylab = "Probabilidad de supervivencia total",
surv.median.line = "hv")
Buscar puntos de corte para PTV y GTV
library(cutpointr)
dbmax$PTV = as.double(dbmax$PTV)
opcpptv <- cutpointr(dbmax, PTV, r6m, direction = "<=", pos_class = "Respuesta",
neg_class = "Progresión", method = maximize_metric, metric = youden)
opcpptv
## # A tibble: 1 × 16
## direction optimal_cutpoint method youden acc sensitivity
## <chr> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 <= 12890.6 maximize_metric 0.100000 0.298851 0.933333
## specificity AUC pos_class neg_class prevalence outcome predictor
## <dbl> <dbl> <chr> <chr> <dbl> <chr> <chr>
## 1 0.166667 0.352778 Respuesta Progresión 0.172414 r6m PTV
## data roc_curve boot
## <list> <list> <lgl>
## 1 <tibble [87 × 2]> <roc_cutpointr [72 × 10]> NA
summary(opcpptv)
## Method: maximize_metric
## Predictor: PTV
## Outcome: r6m
## Direction: <=
##
## AUC n n_pos n_neg
## 0.3528 87 15 72
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 12890.6 0.1 0.2989 0.9333 0.1667 14 1 60 12
##
## Predictor summary:
## Data Min. 5% 1st Qu. Median Mean 3rd Qu. 95% Max.
## Overall 51.86 134.1910 608.150 4016.50 6310.481 8990.89 21224.96 56650.02
## Progresión 51.86 125.2585 471.575 2418.40 6229.256 8728.16 22210.63 56650.02
## Respuesta 705.53 854.8820 2290.010 6368.15 6700.361 9179.35 14196.59 17243.91
## SD NAs
## 8679.335 0
## 9291.674 0
## 4970.623 0
pvptv <- cutpointr(dbmax, PTV, r6m, direction = "<=", pos_class = "Respuesta",
neg_class = "Progresión", method = maximize_metric, metric = p_chisquared)
pvptv
## # A tibble: 1 × 16
## direction optimal_cutpoint method p_chisquared acc sensitivity
## <chr> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 <= 9861.56 maximize_metric 0.960628 0.298851 0.8
## specificity AUC pos_class neg_class prevalence outcome predictor
## <dbl> <dbl> <chr> <chr> <dbl> <chr> <chr>
## 1 0.194444 0.352778 Respuesta Progresión 0.172414 r6m PTV
## data roc_curve boot
## <list> <list> <lgl>
## 1 <tibble [87 × 2]> <roc_cutpointr [72 × 10]> NA
#punto de corte optimo = 12890.6
PUNTO DE CORTE OPTIMO PARA PTV = 12890.6 (y=0.1 Sen = 93% spe = 16% AUC=0.35)
dbmax$GTV = as.double(dbmax$GTV)
opcpgtv <- cutpointr(dbmax, GTV, r6m, direction = "<=", pos_class = "Respuesta",
neg_class = "Progresión", method = maximize_metric, metric = youden)
opcpgtv
## # A tibble: 1 × 16
## direction optimal_cutpoint method youden acc sensitivity
## <chr> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 <= 12916.7 maximize_metric 0.125 0.275862 1
## specificity AUC pos_class neg_class prevalence outcome predictor
## <dbl> <dbl> <chr> <chr> <dbl> <chr> <chr>
## 1 0.125 0.355556 Respuesta Progresión 0.172414 r6m GTV
## data roc_curve boot
## <list> <list> <lgl>
## 1 <tibble [87 × 2]> <roc_cutpointr [72 × 10]> NA
summary(opcpgtv)
## Method: maximize_metric
## Predictor: GTV
## Outcome: r6m
## Direction: <=
##
## AUC n n_pos n_neg
## 0.3556 87 15 72
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 12916.7 0.125 0.2759 1 0.125 15 0 63 9
##
## Predictor summary:
## Data Min. 5% 1st Qu. Median Mean 3rd Qu. 95%
## Overall 17.25 115.2470 509.9350 2115.44 10597.684 6026.595 17921.40
## Progresión 17.25 113.3735 421.5875 1741.13 11725.415 5267.185 18036.83
## Respuesta 333.12 593.8070 1823.8850 3776.60 5184.577 8497.140 12195.69
## Max. SD NAs
## 515908.0 55241.928 0
## 515908.0 60708.911 0
## 12916.7 4114.264 0
PUNTO DE CORTE OPTIMO PARA GTV = 12916.7 (y=0.1 Sen = 93% spe = 16% AUC=0.35)
Recodificación puntos de corte
#RECODIFICAMOS y creamos columna nueva para los puntos de corte
dbmax <- dbmax %>% mutate(volPTV = case_when(PTV <= 12890.6 ~ "<=12890.6",
PTV > 12890.6 ~ ">12890.6"))
dbmax <- dbmax %>% mutate(volGTV = case_when(GTV <= 12916.7 ~ "<=12916.7",
PTV > 12916.7 ~ ">12916.7"))
dbmax %>% select(volPTV, PTV, volGTV, GTV,respuesta6m ,meses)
## # A tibble: 87 × 6
## volPTV PTV volGTV GTV respuesta6m meses
## <chr> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 <=12890.6 8497. <=12916.7 8497. 0 10
## 2 <=12890.6 8497. <=12916.7 8497. 0 10
## 3 <=12890.6 8497. <=12916.7 8497. 0 10
## 4 <=12890.6 8476 <=12916.7 542. 1 10
## 5 <=12890.6 8476 <=12916.7 542. 1 10
## 6 <=12890.6 598. <=12916.7 177. 1 36
## 7 <=12890.6 4100. <=12916.7 4100. 1 21
## 8 <=12890.6 4100. <=12916.7 4100. 1 21
## 9 <=12890.6 4100. <=12916.7 4100. 1 21
## 10 <=12890.6 4100. <=12916.7 4100. 1 21
## # … with 77 more rows
CURVAS KM para los puntos de corte
#Modelos PTV y GTV
dbmax
## # A tibble: 87 × 27
## edad ECOG sintomas metas…¹ contr…² ECOG2 sinto…³ ester…⁴ contr…⁵ ECOG3
## <dbl> <dbl> <chr> <chr> <chr> <dbl> <chr> <chr> <chr> <dbl>
## 1 39 1 sin sintomas no si 1 mas de… si si 1
## 2 NA 1 sin sintomas si si 1 cefalea no si 1
## 3 NA 1 sin sintomas si si 1 mas de… si si 1
## 4 59 1 deficit motor no no 1 sin si… no no 1
## 5 NA 1 sin sintomas no no 1 sin si… no no 1
## 6 57 1 sin sintomas no no 4 mas de… si no 4
## 7 45 1 otros no no 1 sin si… no no 1
## 8 NA 1 sin sintomas no no 4 defici… si no 4
## 9 NA 1 sin sintomas no no 1 sin si… no si 0
## 10 NA 1 sin sintomas no no 1 sin si… no si 1
## # … with 77 more rows, 17 more variables: sintomas3 <chr>, esteroides3 <chr>,
## # rescate <chr>, seg6meses <chr>, seguimiento6m <dbl>, seg12meses <chr>,
## # seguimiento12m <dbl>, res6meses <chr>, respuesta <dbl>, meses <dbl>,
## # PTV <dbl>, GTV <dbl>, dosis <dbl>, respuesta6m <dbl>, r6m <chr>,
## # volPTV <chr>, volGTV <chr>, and abbreviated variable names
## # ¹`metastasis extra`, ²`control sistemico`, ³sintomas2, ⁴esteroides,
## # ⁵controlsistemico2
PTV <- survfit(Surv(meses, respuesta6m) ~ volPTV, data=dbmax)
GTV <- survfit(Surv(meses, respuesta6m) ~ volGTV, data = dbmax)
Grafica Kaplan Meier dosis PTV
PTV
## Call: survfit(formula = Surv(meses, respuesta6m) ~ volPTV, data = dbmax)
##
## n events median 0.95LCL 0.95UCL
## volPTV=<=12890.6 74 60 14 14 21
## volPTV=>12890.6 13 12 17 7 NA
ggsurvplot(PTV,
xlab = "Meses",
ylab = "Probabilidad de supervivencia total",
pval= TRUE,
pval.method=TRUE)
Grafica Kaplan Meier dosis GTV
GTV
## Call: survfit(formula = Surv(meses, respuesta6m) ~ volGTV, data = dbmax)
##
## n events median 0.95LCL 0.95UCL
## volGTV=<=12916.7 78 63 14 14 21
## volGTV=>12916.7 9 9 29 10 NA
ggsurvplot(GTV,
xlab = "Meses",
ylab = "Probabilidad de supervivencia total",
pval= TRUE,
pval.method=TRUE)
Punto de corte para dosis = 21
dbmax <- dbmax %>% mutate(dosisrt = case_when(dosis <= 21 ~ "<=21",
dosis > 21 ~ ">21"))
Modelo para dosis punto de corte maxito
dosis <- survfit(Surv(meses, respuesta6m) ~ dosisrt, data=dbmax)
dosis
## Call: survfit(formula = Surv(meses, respuesta6m) ~ dosisrt, data = dbmax)
##
## n events median 0.95LCL 0.95UCL
## dosisrt=<=21 28 25 8 6 21
## dosisrt=>21 59 47 14 14 21
Curva y analisis KM
ggsurvplot(dosis,
xlab = "Meses",
ylab = "Probabilidad de supervivencia total",
pval= TRUE,
pval.method=TRUE)
Punto de corte optimo para la dosis
optdosis <- cutpointr(dbmax, dosis, r6m, direction = "<=", pos_class = "Respuesta",
neg_class = "Progresión", method = maximize_metric, metric = youden)
optdosis
## # A tibble: 1 × 16
## direction optimal_cutpoint method youden acc sensitivity
## <chr> <dbl> <chr> <dbl> <dbl> <dbl>
## 1 <= 27 maximize_metric 0.0444444 0.252874 0.933333
## specificity AUC pos_class neg_class prevalence outcome predictor
## <dbl> <dbl> <chr> <chr> <dbl> <chr> <chr>
## 1 0.111111 0.428241 Respuesta Progresión 0.172414 r6m dosis
## data roc_curve boot
## <list> <list> <lgl>
## 1 <tibble [87 × 2]> <roc_cutpointr [10 × 10]> NA
summary(optdosis)
## Method: maximize_metric
## Predictor: dosis
## Outcome: r6m
## Direction: <=
##
## AUC n n_pos n_neg
## 0.4282 87 15 72
##
## optimal_cutpoint youden acc sensitivity specificity tp fn fp tn
## 27 0.0444 0.2529 0.9333 0.1111 14 1 64 8
##
## Predictor summary:
## Data Min. 5% 1st Qu. Median Mean 3rd Qu. 95% Max. SD NAs
## Overall 15 15.0 18 22 21.85057 24 30.0 30 4.090627 0
## Progresión 15 15.0 18 22 21.70833 24 30.0 30 4.240773 0
## Respuesta 15 18.5 22 22 22.53333 24 27.9 30 3.313752 0
Punto de corte con youden = 27
dbmax <- dbmax %>% mutate(dosispco = case_when(dosis <= 27 ~ "<=27",
dosis > 27 ~ ">27"))
dbmax %>% select(dosis,dosispco)
## # A tibble: 87 × 2
## dosis dosispco
## <dbl> <chr>
## 1 24 <=27
## 2 24 <=27
## 3 24 <=27
## 4 24 <=27
## 5 24 <=27
## 6 21 <=27
## 7 18 <=27
## 8 18 <=27
## 9 18 <=27
## 10 18 <=27
## # … with 77 more rows
Analisis supervivencia dosis pc 27
dosis27 <- survfit(Surv(meses, respuesta6m) ~ dosispco, data=dbmax)
dosis27
## Call: survfit(formula = Surv(meses, respuesta6m) ~ dosispco, data = dbmax)
##
## n events median 0.95LCL 0.95UCL
## dosispco=<=27 78 64 14 14 21
## dosispco=>27 9 8 12 7 NA
ggsurvplot(dosis27,
xlab = "Meses",
ylab = "Probabilidad de supervivencia total",
pval= TRUE,
pval.method=TRUE)