Preparación

En primer lugar, cargamos las librerías necesarias

if (!require(readxl)) install.packages('readxl', dependencies = T)
## Loading required package: readxl
if (!require(haven)) install.packages('haven', dependencies = T)
## Loading required package: haven
if (!require(survival)) install.packages('survival', dependencies = T)
## Loading required package: survival
if (!require(ggfortify)) install.packages('ggfortify', dependencies = T)
## Loading required package: ggfortify
## Loading required package: ggplot2
library(readxl) 
library(haven) 
library(survival)
library(ggfortify)

A continuación cargamos la base de datos y observamos tanto los nombres como las características principales

CANCERLUNG <- read_excel("CANCERLUNG.xlsx")
## New names:
## • `` -> `...1`
names(CANCERLUNG)
##  [1] "...1"      "inst"      "time"      "status"    "age"       "sex"      
##  [7] "ph.ecog"   "ph.karno"  "pat.karno" "meal.cal"  "wt.loss"
summary(CANCERLUNG)
##       ...1            inst                time            status      
##  Min.   :  1.00   Length:228         Min.   :   5.0   Min.   :0.0000  
##  1st Qu.: 57.75   Class :character   1st Qu.: 166.8   1st Qu.:0.0000  
##  Median :114.50   Mode  :character   Median : 255.5   Median :1.0000  
##  Mean   :114.50                      Mean   : 305.2   Mean   :0.7237  
##  3rd Qu.:171.25                      3rd Qu.: 396.5   3rd Qu.:1.0000  
##  Max.   :228.00                      Max.   :1022.0   Max.   :1.0000  
##       age             sex           ph.ecog            ph.karno        
##  Min.   :39.00   Min.   :0.0000   Length:228         Length:228        
##  1st Qu.:56.00   1st Qu.:0.0000   Class :character   Class :character  
##  Median :63.00   Median :0.0000   Mode  :character   Mode  :character  
##  Mean   :62.45   Mean   :0.3947                                        
##  3rd Qu.:69.00   3rd Qu.:1.0000                                        
##  Max.   :82.00   Max.   :1.0000                                        
##   pat.karno           meal.cal           wt.loss         
##  Length:228         Length:228         Length:228        
##  Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character  
##                                                          
##                                                          
##

Observamos que hay un error en la primera columna, por lo que la eliminamos

CANCERLUNG <- CANCERLUNG [ , -1]

Para poder trabajar más comodamente, usamos las función attach

attach(CANCERLUNG)

Observación de las variables Sex y Status

Hacemos una nueva base de datos con dichas variables para poder realizar mejor la observación

Ej.a <- CANCERLUNG [ , c(3,5)]

Convertimos las variables

Ej.a$status <- as.factor(Ej.a$status)
levels(Ej.a$status) = c("Censurado","Muerto")
Ej.a$sex <- as.factor(Ej.a$sex)
levels(Ej.a$sex) = c("Hombre","Mujer")

Graficamos ambas variables

par(mfrow = c(1, 2))
plot(Ej.a$status, col = c ("Green", "Black"))
plot(Ej.a$sex, col = c ("Blue", "Pink"))

Observamos ambas variables y realizamos una tabla de contingencia con dichas variables. De esta forma comprobamos que existen 63 observaciones censuradas.

table(Ej.a$status)
## 
## Censurado    Muerto 
##        63       165
table(Ej.a$sex)
## 
## Hombre  Mujer 
##    138     90
table(Ej.a)
##            sex
## status      Hombre Mujer
##   Censurado     26    37
##   Muerto       112    53

Tal como podemos observar a continación, el porcentaje de mujeres que han censurado es mucho mayor

print(paste("Porcentaje de hombres que han censurado:",(26*100)/(26+112)))
## [1] "Porcentaje de hombres que han censurado: 18.8405797101449"
print(paste("Porcentaje de muejeres que han censurado:",(37*100)/(37+53)))
## [1] "Porcentaje de muejeres que han censurado: 41.1111111111111"

Graficamos también la tabla comparativa

plot(table(Ej.a), col = c("Blue", "Pink"), main = "Comparación entre las variables Sex y Status")

Curva de supervivencia de Kaplan-Meier

Recreamos la curva de supervivencia de Kaplan-Meier utilizando la función survfit() dentro de la biblioteca survival.

fit.surv <- survfit(Surv(time, status) ~ 1)
autoplot(fit.surv)

Como podemos observar, menos del 10% de los pacientes seguian vivos después del estudio.

Curva de Kaplan-Meier en función del sexo

A continación hacemos la diferenciación en función del sexo. A simpele vista ya podemos observar que los hombres de este estudio fallecieron antes.

sex.fit.surv <- survfit(Surv(time, status) ~ sex)
autoplot(sex.fit.surv)

Test log-rank

Para comprobar la diferencia entre ambos grupos relizaremos el test long-rank con la función survdiff()

logrank.test <- survdiff(Surv(time, status) ~ sex)
logrank.test
## Call:
## survdiff(formula = Surv(time, status) ~ sex)
## 
##         N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=0 138      112     91.6      4.55      10.3
## sex=1  90       53     73.4      5.68      10.3
## 
##  Chisq= 10.3  on 1 degrees of freedom, p= 0.001

El resultado del p-valor es de 0.001, lo que muestra que en efecto existen diferencias estadísticas

Modelo de riesgos proporcionales de Cox

Inicialmento conisderaremos un modelo que ultiliza la varaible sex como único predictor

fit.cox <- coxph(Surv(time, status) ~ sex)
summary(fit.cox)
## Call:
## coxph(formula = Surv(time, status) ~ sex)
## 
##   n= 228, number of events= 165 
## 
##        coef exp(coef) se(coef)      z Pr(>|z|)   
## sex -0.5310    0.5880   0.1672 -3.176  0.00149 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##     exp(coef) exp(-coef) lower .95 upper .95
## sex     0.588      1.701    0.4237     0.816
## 
## Concordance= 0.579  (se = 0.021 )
## Likelihood ratio test= 10.63  on 1 df,   p=0.001
## Wald test            = 10.09  on 1 df,   p=0.001
## Score (logrank) test = 10.33  on 1 df,   p=0.001

Como los valores del test se han redondeado, vamos a obtener la puntuación excta

summary(fit.cox)$logtest 
##        test          df      pvalue 
## 10.63360473  1.00000000  0.00111051
summary(fit.cox)$waldtest 
##         test           df       pvalue 
## 10.090000000  1.000000000  0.001491229
summary(fit.cox)$sctest 
##         test           df       pvalue 
## 10.325148496  1.000000000  0.001312297

Tal como hemos podido observar, tanto en el test de razón de verosimilitud, como en el de _Wald: y el de score, el p-valor es significativo

Modelo con todos los posibles predictores

A contiunación, ajustamos un modelo con todos los posibles predictores de la base de datos. Al realizar una primera prospectiva, observamos que muchas variables tienes forma de character lo cual satura el sistema al crear demasiadas interacciones, por lo que, en primer lugar, adecuamos la base de datos

CANCERLUNG3 <- CANCERLUNG
CANCERLUNG3$inst <- as.numeric(CANCERLUNG3$inst)
## Warning: NAs introduced by coercion
CANCERLUNG3$ph.ecog <- as.numeric(CANCERLUNG3$ph.ecog)
## Warning: NAs introduced by coercion
CANCERLUNG3$ph.karno <- as.numeric(CANCERLUNG3$ph.karno)
## Warning: NAs introduced by coercion
CANCERLUNG3$pat.karno <- as.numeric(CANCERLUNG3$pat.karno)
## Warning: NAs introduced by coercion
CANCERLUNG3$meal.cal <- as.numeric(CANCERLUNG3$meal.cal)
## Warning: NAs introduced by coercion
CANCERLUNG3$wt.loss <- as.numeric(CANCERLUNG3$wt.loss)
## Warning: NAs introduced by coercion

A continuación ya podemos realizar el estudio completo

fit.all <- coxph(Surv(time, status) ~ . , data = CANCERLUNG3)
fit.all
## Call:
## coxph(formula = Surv(time, status) ~ ., data = CANCERLUNG3)
## 
##                 coef  exp(coef)   se(coef)      z        p
## inst      -3.037e-02  9.701e-01  1.312e-02 -2.315 0.020619
## age        1.281e-02  1.013e+00  1.194e-02  1.073 0.283403
## sex       -5.666e-01  5.674e-01  2.014e-01 -2.814 0.004890
## ph.ecog    9.074e-01  2.478e+00  2.386e-01  3.803 0.000143
## ph.karno   2.658e-02  1.027e+00  1.163e-02  2.286 0.022231
## pat.karno -1.091e-02  9.891e-01  8.141e-03 -1.340 0.180160
## meal.cal   2.602e-06  1.000e+00  2.677e-04  0.010 0.992244
## wt.loss   -1.671e-02  9.834e-01  7.911e-03 -2.112 0.034647
## 
## Likelihood ratio test=33.7  on 8 df, p=4.609e-05
## n= 167, number of events= 120 
##    (61 observations deleted due to missingness)

Así pues, tal como se puede observar, las variables que influyen en el modelo son el sexo y la escala de ECOG, seguido de la institución y la perdidia de peso en los últimos seis meses. El análisis de las variables nos muestra que la escala ECOG es casi cinco veces más indicadora que el sexo.

Modelo alternativo

Como alternativa a lo realizado, usaremos el modelo de regresión aditiva de Aalen para daros censurados. En este caso nos muestra que los únicos datos que muestran un p-valor significativo sonel sexo y la escala de ECOG.

aa_fit <-aareg(Surv(time, status) ~ . ,
               data = CANCERLUNG3)

summary(aa_fit)
##               slope      coef se(coef)       z       p
## Intercept  7.82e-03  8.37e-03 1.71e-02  0.4880 0.62500
## inst      -6.25e-05 -1.75e-04 1.25e-04 -1.4000 0.16100
## age       -1.06e-05  5.13e-06 1.16e-04  0.0442 0.96500
## sex       -2.58e-03 -5.99e-03 1.90e-03 -3.1600 0.00160
## ph.ecog    2.67e-03  7.22e-03 2.61e-03  2.7700 0.00563
## ph.karno   2.51e-05  1.94e-04 1.23e-04  1.5800 0.11500
## pat.karno -5.74e-05 -1.12e-04 9.59e-05 -1.1700 0.24400
## meal.cal  -1.45e-06 -3.17e-06 2.33e-06 -1.3600 0.17300
## wt.loss   -6.73e-05 -1.48e-04 8.16e-05 -1.8200 0.06920
## 
## Chisq=22.99 on 8 df, p=0.00337; test weights=aalen

Para finalizar, graficaos dicho modelo para observar como las covariables cambian con el tiempo

autoplot(aa_fit)