Analisis de Supervivencia
Ricardo Rodriguez Ovilla
4/2/2021
Librerias
library(KMsurv); library(LifeTables); library(mclust); library(foreign)
library(splines); library(survival); library(nlme); library(muhaz)
library(TH.data); library(imputeTS); library(dplyr); library(tidyverse);library(ggplot2)
library(highcharter);library(survminer)Este es una propuesta para la modelacion de los tiempos de supervivencia.
Las variables numericas son las siguientes: Age, meal.cal, wt.loss
Datos Faltantes.
Antes de realizar o decidir que variables sufriran una codificacion procedemos ha hacer analisis exploratorio de datos.
En este caso el set de datos cancer posee 6 variables con datos faltantes las cuales seran imputadas a partir de de algun metodo , aun decidire cual.
map_dbl(cancer, .f = function(x){sum(is.na(x))})## inst time status age sex ph.ecog ph.karno pat.karno
## 1 0 0 0 0 1 1 3
## meal.cal wt.loss
## 47 14Para una mejor comprension de los datos faltantes presentes en las variables podemos ver que la variable meal.cal contiene un 20% de los datos …….continuar
datos_long <- cancer %>% gather(key = "variable", value = "valor")
datos_long %>%
group_by(variable) %>%
summarize(porcentaje_NA = 100 * sum(is.na(valor)) / length(valor)) %>%
ggplot(aes(x = reorder(variable, desc(porcentaje_NA)), y = porcentaje_NA)) +
geom_col() +
labs(title = "Porcentaje valores ausentes por variable en Cancer Data",
x = "Variable", y = "Porcentaje NAs") +
theme_bw()
Imputacion
datoscancer <- na_interpolation(cancer)Censoring and Event Plots
o <- datoscancer %>%
arrange(sex) %>%
mutate(index = 1:n())
## En este caso las leyendas las debemos mejorar lo hago despues. agregar titulo y otro tem mas bonito.
ggplot(o,
aes(xend = 0,
y = index,
x = time,
yend = index,
colour = factor(sex),
shape = factor(status))) +
geom_segment() +
geom_point() 
o_2 <- datoscancer %>%
arrange(ph.ecog) %>%
mutate(index = 1:n())
## En este caso las leyendas las debemos mejorar lo hago despues. En este casp hay NA´s
ggplot(o_2,
aes(xend = 0,
y = index,
x = time,
yend = index,
colour = factor(ph.ecog),
shape = factor(status))) +
geom_segment() +
geom_point() 
o_3 <- datoscancer %>%
arrange(ph.karno) %>%
mutate(index = 1:n())
## En este caso las leyendas las debemos mejorar lo hago despues. En este casp hay NA´s
ggplot(o_3,
aes(xend = 0,
y = index,
x = time,
yend = index,
colour = factor(ph.karno),
shape = factor(status))) +
geom_segment() +
geom_point() 
o_4 <- datoscancer %>%
arrange(pat.karno) %>%
mutate(index = 1:n())
## En este caso las leyendas las debemos mejorar lo hago despues. En este casp hay NA´s
ggplot(o_4,
aes(xend = 0,
y = index,
x = time,
yend = index,
colour = factor(pat.karno),
shape = factor(status))) +
geom_segment() +
geom_point() 
Histogramas
Tiempo del evento
hchart(density(cancer$time), type = "area", color = "#B71C1C", name = "Time") %>% hc_add_theme(hc_theme_elementary())Riesgos Proporcionales.
Los resultados de la regresión de Cox se pueden interpretar de la siguiente manera:
Significancia estadística: La columna marcada con “pvalue” da el valor estadístico de Wald. Corresponde a la razón de cada coeficiente de regresión a su error estándar (z = coef / se (coef)). La estadística de Wald evalúa si el coeficiente \(\beta\) de una variable dada es estadísticamente significativamente diferente de 0.
Los coeficientes de regresión: La segunda característica a tener en cuenta en los resultados del modelo de Cox es el signo de los coeficientes de regresión (“beta”). Un signo positivo significa que el peligro (riesgo de muerte) es mayor y, por tanto, el pronóstico peor, para los sujetos con valores más altos de esa variable.
Intervalos de confianza de las razones de riesgo. El resultado resumido también proporciona intervalos de confianza superior e inferior del 95% para el índice de riesgo (exp (coef)). Los coeficientes exponenciados (exp (coef)), también conocidos como cocientes de riesgo, dan el tamaño del efecto de las covariables.
covariates <- c("sex","wt.loss","meal.cal","age","ph.karno","pat.karno","ph.ecog")
univ_formulas <- sapply(covariates,
function(x) as.formula(paste('Surv(time, status)~', x)))
univ_models <- lapply( univ_formulas, function(x){coxph(x, data = cancer)})
# datos
univ_results <- lapply(univ_models,
function(x){
x <- summary(x)
p.value<-signif(x$wald["pvalue"], digits=2)
wald.test<-signif(x$wald["test"], digits=2)
beta<-signif(x$coef[1], digits=2);#coeficiente beta
HR <-signif(x$coef[2], digits=2);#exp(beta)
HR.confint.lower <- signif(x$conf.int[,"lower .95"], 2)
HR.confint.upper <- signif(x$conf.int[,"upper .95"],2)
HR <- paste0(HR, " (",
HR.confint.lower, "-", HR.confint.upper, ")")
res<-c(beta, HR, wald.test, p.value)
names(res)<-c("beta", "HR (95% CI for HR)", "wald.test",
"p.value")
return(res)
#regresa (exp(cbind(coef(x),confint(x))))
})
res <- t(as.data.frame(univ_results, check.names = F))
as.data.frame(res)## beta HR (95% CI for HR) wald.test p.value
## sex -0.53 0.59 (0.42-0.82) 10 0.0015
## wt.loss 0.0013 1 (0.99-1) 0.05 0.83
## meal.cal -0.00012 1 (1-1) 0.29 0.59
## age 0.019 1 (1-1) 4.1 0.042
## ph.karno -0.016 0.98 (0.97-1) 7.9 0.005
## pat.karno -0.02 0.98 (0.97-0.99) 13 0.00028
## ph.ecog 0.48 1.6 (1.3-2) 18 2.7e-05El resultado anterior muestra los coeficientes beta de regresión, los tamaños del efecto (expresados como cocientes de riesgo) y la significación estadística de cada una de las variables en relación con la supervivencia general. Cada factor se evalúa mediante regresiones de Cox univariadas independientes.
De arriba:
La variable sexo, edad, ph.karno y pat.karno son significantes para las regresiones univariadas del modelo.
Los coeficientes de las variables numericas que tienen coeficientes positivos nos indican que a medida de que esta aumenta el riesgo del evento se vuelve mas grande para el set de datos cancer wt.loss y age estan en esta clasificacion. Mientras que disminuir las calorias, en este caso una buena dieta, reduce el risgo del evento.
La variable sex tiene un coeficiente negativo por lo que se interpreta com; ser mujer (sexo = 2) se asocia con una mejor supervivencia.
Multivariante
Ahora, queremos describir cómo los factores impactan conjuntamente en la supervivencia. Para responder a esta pregunta, realizaremos un análisis de regresión de Cox multivariante. Solo incluiremos las varibales que parecer ser significativas en las regresiones univariantes.
# Data preparation and computing cox model
cox_mult <- coxph(Surv(time, status) ~ age + sex + ph.karno+pat.karno, data = lung)
summary(cox_mult )## Call:
## coxph(formula = Surv(time, status) ~ age + sex + ph.karno + pat.karno,
## data = lung)
##
## n= 224, number of events= 161
## (4 observations deleted due to missingness)
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age 0.010658 1.010715 0.009530 1.118 0.26345
## sex -0.508150 0.601608 0.169298 -3.002 0.00269 **
## ph.karno -0.004552 0.995459 0.006972 -0.653 0.51382
## pat.karno -0.016793 0.983347 0.006509 -2.580 0.00988 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age 1.0107 0.9894 0.9920 1.0298
## sex 0.6016 1.6622 0.4317 0.8383
## ph.karno 0.9955 1.0046 0.9819 1.0092
## pat.karno 0.9833 1.0169 0.9709 0.9960
##
## Concordance= 0.644 (se = 0.025 )
## Likelihood ratio test= 24.18 on 4 df, p=7e-05
## Wald test = 23.76 on 4 df, p=9e-05
## Score (logrank) test = 24.35 on 4 df, p=7e-05En este caso notamos que la variable age y ph.karno conjuntamente no parecer ser significativa por lo que decidimos hacer un tratamiento buscando agrupaciones que nos ayuden a pasar la prueba de significancia, no se tuvo exito en este tratamiento ya que a pesar de hacer agrupaciones no podiamos conseguir la mejor agrupacion.
datoscancer <- na_interpolation(cancer)
datos = datoscancer
datos <- datos %>%
mutate(Age_grupo = case_when(age >= 39 & age <= 43 ~ "adulto",
age > 43 ~ "anciano"))
datos$Age_grupo <- as.factor(datos$Age_grupo)
# Modificamos la variable ph
datos <- datos %>%
mutate(ph.karno_G = case_when(ph.karno >= 50 & age <= 60 ~ "bueno",
age > 60 ~ "malo"))
datos$ph.karno_G <- as.factor(datos$ph.karno_G)
## Datos Modificados
# Data preparation and computing cox model
cox_mult <- coxph(Surv(time, status) ~ Age_grupo+ph.karno_G + sex +pat.karno, data = datos)
summary(cox_mult )## Call:
## coxph(formula = Surv(time, status) ~ Age_grupo + ph.karno_G +
## sex + pat.karno, data = datos)
##
## n= 228, number of events= 165
##
## coef exp(coef) se(coef) z Pr(>|z|)
## Age_grupoanciano 0.823065 2.277469 1.010493 0.815 0.415348
## ph.karno_Gmalo 0.102561 1.108005 0.162538 0.631 0.528041
## sex -0.524985 0.591564 0.167335 -3.137 0.001705 **
## pat.karno -0.019710 0.980483 0.005604 -3.517 0.000436 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## Age_grupoanciano 2.2775 0.4391 0.3143 16.5039
## ph.karno_Gmalo 1.1080 0.9025 0.8057 1.5237
## sex 0.5916 1.6904 0.4262 0.8212
## pat.karno 0.9805 1.0199 0.9698 0.9913
##
## Concordance= 0.639 (se = 0.024 )
## Likelihood ratio test= 25.23 on 4 df, p=5e-05
## Wald test = 24.31 on 4 df, p=7e-05
## Score (logrank) test = 24.95 on 4 df, p=5e-05Nos quedaremos con el siguiente modelo.
# Data preparation and computing cox model
cox_final <- coxph(Surv(time, status) ~ sex+pat.karno+ph.ecog, data =datos)
summary(cox_final)## Call:
## coxph(formula = Surv(time, status) ~ sex + pat.karno + ph.ecog,
## data = datos)
##
## n= 228, number of events= 165
##
## coef exp(coef) se(coef) z Pr(>|z|)
## sex -0.548909 0.577579 0.167431 -3.278 0.00104 **
## pat.karno -0.009922 0.990127 0.006752 -1.470 0.14169
## ph.ecog 0.372220 1.450952 0.135080 2.756 0.00586 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## sex 0.5776 1.7314 0.4160 0.8019
## pat.karno 0.9901 1.0100 0.9771 1.0033
## ph.ecog 1.4510 0.6892 1.1135 1.8907
##
## Concordance= 0.658 (se = 0.025 )
## Likelihood ratio test= 31.5 on 3 df, p=7e-07
## Wald test = 31.71 on 3 df, p=6e-07
## Score (logrank) test = 32.53 on 3 df, p=4e-07El valor \(p_{value}\) para las tres pruebas generales (probabilidad, Wald y puntuación) es significativo, lo que indica que el modelo es significativo. Estas pruebas evalúan la hipótesis nula de que todas las betas \(\beta\) son 0. En el ejemplo anterior, las estadísticas de la prueba están muy de acuerdo y la hipótesis nula ómnibus se rechaza rotundamente.
En el análisis de Cox multivariado, las covariables sexo y pat.karno siguen siendo significativas (\(p <0.05\)).
El \(p_{value}\) para el sexo es 0.001545, con un índice de riesgo HR = exp (coef) = 0,58, lo que indica una fuerte relación entre el sexo de los pacientes y un menor riesgo de muerte. Las razones de riesgo de las covariables se pueden interpretar como efectos multiplicativos sobre el riesgo. Por ejemplo, si se mantienen constantes las otras covariables, ser mujer (sexo = 2) reduce el riesgo en un factor de 0.58 o 42%. Concluimos que ser mujer se asocia con un buen pronóstico.
Validacion de Supuestos.
Prueba de la suposición de riesgos proporcionales
El supuesto de riesgos proporcionales (PH) se puede verificar mediante pruebas estadísticas y diagnósticos gráficos basados en los residuos de Schoenfeld escalados.
En principio, los residuos de Schoenfeld son independientes del tiempo que muestra un patrón no aleatorio contra el tiempo es evidencia de una violación de la suposición de PH.
El supuesto de riesgo proporcional está respaldado por una relación no significativa entre los residuos y el tiempo, y refutado por una relación significativa.
test.ph <- cox.zph(cox_final)
test.ph## chisq df p
## sex 2.96 1 0.086
## pat.karno 5.17 1 0.023
## ph.ecog 1.86 1 0.172
## GLOBAL 7.42 3 0.060De la salida anterior, la prueba no es estadísticamente significativa para cada una de las covariables, y la prueba global tampoco es estadísticamente significativa. Por lo tanto, podemos asumir los riesgos proporcionales.
ggcoxzph(test.ph)
En la figura anterior, la línea continua es un ajuste de spline suavizado al gráfico, con las líneas discontinuas que representan una banda de error estándar de +/- 2 alrededor del ajuste.
Cabe resaltar que las desviaciones sistemáticas de una línea horizontal son indicativas de peligros no proporcionales, ya que los peligros proporcionales suponen que las estimaciones \(\beta_1\) , \(\beta_2\) y \(\beta_3\) no varían mucho con el tiempo.
De la inspección gráfica, no hay patrón con el tiempo. El supuesto de riesgos proporcionales parece estar respaldado por las covariables sexo (que es, recordemos, un factor de dos niveles, que representa las dos bandas en el gráfico), pérdida de peso y edad.
Prueba de observaciones influyentes
Para probar observaciones influyentes o valores atípicos, podemos visualizar:
los residuales de desviación o los valores de dfbeta La función ggcoxdiagnostics () proporciona una solución conveniente para verificar las observaciones influyentes. El formato simplificado es el siguiente:
ggcoxdiagnostics(cox_final, type = "dfbeta",
linear.predictions = FALSE, ggtheme = theme_bw())
Los gráficos de índice anteriores muestran que la comparación de las magnitudes de los valores dfbeta más grandes con los coeficientes de regresión sugiere que ninguna de las observaciones es muy influyente individualmente, a pesar de que algunos de los valores dfbeta para la edad y la pérdida de peso son grandes en comparación con los demás.
También es posible comprobar los valores atípicos visualizando los residuos de desviación. El residuo de desviación es una transformación normalizada del residuo de martingala. Estos residuos deben distribuirse aproximadamente simétricamente alrededor de cero con una desviación estándar de 1.
ggcoxdiagnostics(cox_final, type = "martingale",
linear.predictions = FALSE, ggtheme = theme_bw())
Prueba de no linealidad
A menudo, asumimos que las covariables continuas tienen una forma lineal. Sin embargo, esta suposición debe verificarse.
Trazar los residuos de Martingala contra covariables continuas es un enfoque común que se usa para detectar la no linealidad o, en otras palabras, para evaluar la forma funcional de una covariable. Para una covariable continua dada, los patrones en la gráfica pueden sugerir que la variable no se ajusta correctamente.
La no linealidad no es un problema para las variables categóricas, por lo que solo examinamos los gráficos de los residuos de martingala y los residuos parciales frente a una variable continua.
Los residuos de martingala pueden presentar cualquier valor en el rango (-INF, +1):
Un valor de los residuos martinguale cercano a 1 representa a los individuos que “murieron demasiado pronto”, y los valores negativos grandes corresponden a los individuos que “vivieron demasiado tiempo”.
ggcoxfunctional(Surv(time, status) ~ pat.karno+ log(pat.karno) + sqrt(pat.karno), data = datos)
ggcoxfunctional(Surv(time, status) ~ age + log(age) + sqrt(age), data = datos)