Análisis para la tesis de mi amiguito Max :)

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

Exploratorio

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

ANALISIS DE SUPERVIVENCIA

graficas de kaplan-Meier

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)