Análisis de Supervivencia con algoritmos de Machine Learning
library(data.table)
library(readr)
library(ggplot2)
library(survival)
library(survminer)
library(caret)
library(gdata)
library(auxfun)
library(mice)
library(ggmice)
library(randomForestSRC)
library(ggRandomForests)
library(Hmisc)
library(scales)
library(dplyr)
library(e1071)
library(survivalsvm)
library(kableExtra)
library(compareC)
# Carga de datos
lake <- read_delim("~/Máster Bioinformática/TFM/DATOS/lake.csv",
delim = ";", escape_double = FALSE, locale = locale(decimal_mark = ","),
trim_ws = TRUE)
lake3 <- readxl::read_excel("~/Máster Bioinformática/TFM/DATOS/lake2.xlsx")
names(lake3)[32] <- "tiempo"Limpieza de datos
El primer paso antes de realizar el análisis es comprobar cuántos valores perdidos tenemos y si es necesario llevar a cabo alguna transformación.
Comenzamos visualizando el patrón general de NAs:
plot_pattern(lake, rotate = TRUE)
Calculamos el % en la base de datos:
# % NAs en toda la base de datos
porc_NA <- round((sum(is.na(lake)) / (ncol(lake) * nrow(lake))) *100,2)
porc_NA## [1] 42.21Tenemos un 42.21% de valores perdidos.
A continuación eliminamos algunas variables no relevantes. Este paso dependerá de los datos con los que trabajemos.
En este caso, eliminamos la procedencia de los datos, el nombre del paciente, el número de paciente, el número de visita, la fecha de nacimiento (ya tenemos una variabe edad) y los factores de riesgo individuales, que ya que se encuentran en la variable factor_riesgo_total. También la fecha inicio Lake, las fechas correspondientes a las semanas de seguimiento y la fecha VIH porque ya tenemos el tiempo en meses.
lake <- lake[,-c(1:5,7,10:15,17,18,19,23:24,60:61,97:98,137:138,174:175, 174:175)]Solucionamos las incongruencias con los tiempos de infección, que muestra valores negativos:
summary(lake$tpo_vih_meses)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## -183.667 1.533 4.900 20.960 21.133 285.300 7lake$tpo_vih_meses <- ifelse(lake$tpo_vih_meses <=0, NA, lake$tpo_vih_meses)Una vez solucionado nos queda así:
summary(lake$tpo_vih_meses)## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 0.2333 2.1583 6.1667 28.1088 23.5000 285.3000 18Calculamos el porcentaje de NAs por variable y eliminamos aquellas que presentan un 20% o más.
porc_missing <- round(unlist(lapply(lake, function(x) sum(is.na(x))))/
nrow(lake)*100,2)
porc_missing## sexo especificar a19
## 5.17 91.38 100.00
## a28 Estado edad
## 100.00 100.00 0.86
## Grupo CargaViral_0 CD4A_0
## 0.00 20.69 13.79
## CD4P_0 CD8A_0 CD8P_0
## 13.79 13.79 13.79
## Hematocrito_0 Hemoglobina_0 Plaquetas_0
## 11.21 11.21 11.21
## Leucocitos_0 LinfosTotales_0 Glucosa_mg_0
## 11.21 11.21 13.79
## Urea_mg_0 Creatinina_mumol_0 Sodio_0
## 17.24 14.66 31.03
## Potasio_0 Cloro_0 Calcio_0
## 31.03 67.24 46.55
## Bilirrubina_mumol_0 GPT_0 GOT_0
## 15.52 13.79 14.66
## GGT_0 ProteinasTotales_0 Albumina_0
## 16.38 27.59 37.93
## Colesterol_mg_0 LDL_mg_0 HDL_mg_0
## 14.66 40.52 37.07
## Trigliceridos_mg_0 Amilasa_0 pH_0
## 14.66 53.45 77.59
## Bicarbonato_0 AcidoLactico_0 AcidoPiruvico_0
## 73.28 75.86 99.14
## VHC_0 VHB_0 Embarazo_0
## 7.76 10.34 78.45
## CargaViral_12 CD4A_12 CD4P_12
## 21.55 19.83 19.83
## CD8A_12 CD8P_12 Hematocrito_12
## 20.69 19.83 20.69
## Hemoglobina_12 Plaquetas_12 Leucocitos_12
## 20.69 20.69 20.69
## LinfosTotales_12 Glucosa_mg_12 Urea_mg_12
## 19.83 18.97 18.97
## Creatinina_mumol_12 Sodio_12 Potasio_12
## 18.97 30.17 30.17
## Cloro_12 Calcio_12 Bilirrubina_mumol_12
## 58.62 43.10 18.97
## GPT_12 GOT_12 GGT_12
## 18.97 18.97 20.69
## ProteinasTotales_12 Albumina_12 Colesterol_mg_12
## 25.00 34.48 18.97
## LDL_mg_12 HDL_mg_12 Trigliceridos_mg_12
## 26.72 24.14 18.97
## Amilasa_12 pH_12 Bicarbonato_12
## 40.52 46.55 47.41
## AcidoLactico_12 AcidoPiruvico_12 VHC_12
## 41.38 86.21 100.00
## VHB_12 Embarazo_12 CargaViral_24
## 100.00 100.00 37.93
## CD4A_24 CD4P_24 CD8A_24
## 34.48 35.34 34.48
## CD8P_24 Hematocrito_24 Hemoglobina_24
## 35.34 34.48 34.48
## Plaquetas_24 Leucocitos_24 LinfosTotales_24
## 34.48 34.48 36.21
## Glucosa_mg_24 Urea_mg_24 Creatinina_mumol_24
## 34.48 33.62 33.62
## Sodio_24 Potasio_24 Cloro_24
## 39.66 39.66 62.07
## Calcio_24 Bilirrubina_mumol_24 GPT_24
## 46.55 34.48 33.62
## GOT_24 GGT_24 ProteinasTotales_24
## 33.62 33.62 38.79
## Albumina_24 Colesterol_mg_24 LDL_mg_24
## 44.83 33.62 42.24
## HDL_mg_24 Trigliceridos_mg_24 Amilasa_24
## 37.07 33.62 54.31
## pH_24 Bicarbonato_24 AcidoLactico_24
## 60.34 54.31 51.72
## AcidoPiruvico_24 VHC_24 VHB_24
## 87.93 100.00 100.00
## Embarazo_24 cv50_0 cv50_12
## 97.41 20.69 21.55
## cv50_24 CargaViral_36 CD4A_36
## 37.93 50.86 50.86
## CD4P_36 CD8A_36 CD8P_36
## 50.86 50.86 50.86
## Hematocrito_36 Hemoglobina_36 Plaquetas_36
## 51.72 51.72 51.72
## Leucocitos_36 LinfosTotales_36 Glucosa_mg_36
## 51.72 52.59 50.86
## Urea_mg_36 Creatinina_mumol_36 Sodio_36
## 50.86 50.86 53.45
## Potasio_36 Cloro_36 Calcio_36
## 53.45 62.93 62.07
## Bilirrubina_mumol_36 GPT_36 GOT_36
## 50.86 50.86 50.86
## GGT_36 ProteinasTotales_36 Albumina_36
## 51.72 56.03 55.17
## Colesterol_mg_36 LDL_mg_36 HDL_mg_36
## 50.86 59.48 52.59
## Trigliceridos_mg_36 Amilasa_36 pH_36
## 50.86 57.76 69.83
## Bicarbonato_36 AcidoLactico_36 AcidoPiruvico_36
## 62.93 62.93 87.07
## VHC_36 VHB_36 Embarazo_36
## 43.97 43.97 43.97
## CargaViral_48 CD4A_48 CD4P_48
## 60.34 60.34 60.34
## CD8A_48 CD8P_48 Hematocrito_48
## 60.34 60.34 59.48
## Hemoglobina_48 Plaquetas_48 Leucocitos_48
## 59.48 59.48 60.34
## LinfosTotales_48 Glucosa_mg_48 Urea_mg_48
## 60.34 59.48 59.48
## Creatinina_mumol_48 Sodio_48 Potasio_48
## 59.48 61.21 61.21
## Cloro_48 Calcio_48 Bilirrubina_mumol_48
## 68.10 65.52 59.48
## GPT_48 GOT_48 GGT_48
## 59.48 59.48 59.48
## ProteinasTotales_48 Albumina_48 Colesterol_mg_48
## 62.07 62.93 59.48
## LDL_mg_48 HDL_mg_48 Trigliceridos_mg_48
## 63.79 60.34 59.48
## Amilasa_48 pH_48 Bicarbonato_48
## 66.38 70.69 68.97
## AcidoLactico_48 AcidoPiruvico_48 VHC_48
## 67.24 88.79 56.90
## VHB_48 Embarazo_48 cv50_36
## 56.90 56.90 50.86
## cv50_48 tpo_vih_meses factor_riesgo_total
## 60.34 15.52 6.03
## diff_cd4_48_0 diff_cd4p_48_0 diff_col_48_0
## 66.38 65.52 67.24
## diff_HDL_48_0 diff_LDL_48_0
## 75.00 79.31dd_limpio <- remove.vars(lake,names(porc_missing[porc_missing >=20]))##
## Changing in lake
## Dropping variables: especificar, a19, a28, Estado, CargaViral_0, Sodio_0, Potasio_0, Cloro_0, Calcio_0, ProteinasTotales_0, Albumina_0, LDL_mg_0, HDL_mg_0, Amilasa_0, pH_0, Bicarbonato_0, AcidoLactico_0, AcidoPiruvico_0, Embarazo_0, CargaViral_12, CD8A_12, Hematocrito_12, Hemoglobina_12, Plaquetas_12, Leucocitos_12, Sodio_12, Potasio_12, Cloro_12, Calcio_12, GGT_12, ProteinasTotales_12, Albumina_12, LDL_mg_12, HDL_mg_12, Amilasa_12, pH_12, Bicarbonato_12, AcidoLactico_12, AcidoPiruvico_12, VHC_12, VHB_12, Embarazo_12, CargaViral_24, CD4A_24, CD4P_24, CD8A_24, CD8P_24, Hematocrito_24, Hemoglobina_24, Plaquetas_24, Leucocitos_24, LinfosTotales_24, Glucosa_mg_24, Urea_mg_24, Creatinina_mumol_24, Sodio_24, Potasio_24, Cloro_24, Calcio_24, Bilirrubina_mumol_24, GPT_24, GOT_24, GGT_24, ProteinasTotales_24, Albumina_24, Colesterol_mg_24, LDL_mg_24, HDL_mg_24, Trigliceridos_mg_24, Amilasa_24, pH_24, Bicarbonato_24, AcidoLactico_24, AcidoPiruvico_24, VHC_24, VHB_24, Embarazo_24, cv50_0, cv50_12, cv50_24, CargaViral_36, CD4A_36, CD4P_36, CD8A_36, CD8P_36, Hematocrito_36, Hemoglobina_36, Plaquetas_36, Leucocitos_36, LinfosTotales_36, Glucosa_mg_36, Urea_mg_36, Creatinina_mumol_36, Sodio_36, Potasio_36, Cloro_36, Calcio_36, Bilirrubina_mumol_36, GPT_36, GOT_36, GGT_36, ProteinasTotales_36, Albumina_36, Colesterol_mg_36, LDL_mg_36, HDL_mg_36, Trigliceridos_mg_36, Amilasa_36, pH_36, Bicarbonato_36, AcidoLactico_36, AcidoPiruvico_36, VHC_36, VHB_36, Embarazo_36, CargaViral_48, CD4A_48, CD4P_48, CD8A_48, CD8P_48, Hematocrito_48, Hemoglobina_48, Plaquetas_48, Leucocitos_48, LinfosTotales_48, Glucosa_mg_48, Urea_mg_48, Creatinina_mumol_48, Sodio_48, Potasio_48, Cloro_48, Calcio_48, Bilirrubina_mumol_48, GPT_48, GOT_48, GGT_48, ProteinasTotales_48, Albumina_48, Colesterol_mg_48, LDL_mg_48, HDL_mg_48, Trigliceridos_mg_48, Amilasa_48, pH_48, Bicarbonato_48, AcidoLactico_48, AcidoPiruvico_48, VHC_48, VHB_48, Embarazo_48, cv50_36, cv50_48, diff_cd4_48_0, diff_cd4p_48_0, diff_col_48_0, diff_HDL_48_0, diff_LDL_48_0porc_NA <- round((sum(is.na(dd_limpio)) / (ncol(dd_limpio) * nrow(dd_limpio))) *100,2)
porc_NA## [1] 14.14Ahora tenemos un 14.14% de NAs general y el patrón gráfico nos queda como sigue:
plot_pattern(dd_limpio, rotate = TRUE)
gg <- plot_pattern(dd_limpio, rotate = TRUE)
ggtoppt(gg)En este caso tuvimos que calcular el tiempo de seguimiento y determinar la variable éxito/fracaso aparte de la base de datos, por lo que las incoportamos al data.frame:
#Incorporamos el tiempo de seguimiento y el status
dd_limpio$tiempo <- lake3$tiempo #contamos el tiempo de seguimiento total
dd_limpio$status <- lake3$status_sin0
dd_limpio <- dd_limpio[!is.na(dd_limpio$tiempo),]
dd_limpio$status <- ifelse(dd_limpio$status == 1, 0, 1)
dd_limpio$status <- ifelse(is.na(dd_limpio$status) == TRUE, 1, dd_limpio$status)
dd_limpio$edad <- ifelse(dd_limpio$edad == 0, NA, dd_limpio$edad)
table(dd_limpio$tiempo)##
## 0 12 24 36 48
## 13 20 21 11 47table(dd_limpio$status)##
## 0 1
## 78 34Análisis descriptivo
Tranformación de los datos
Comenzamos modificando etiquetas para que las tablas y los gráficos sean más comprensibles.
lake2_2 <- dd_limpio
lake2_2$factor_riesgo_total <- factor(lake2_2$factor_riesgo_total, levels = 1:5,
labels = c("ADVP", "Heterosexual","Homosexual", "Hemofilia", "Otros"))
lake2_2$Grupo <- factor(lake2_2$Grupo, levels = c(-1,0), labels = c("EFV + Kivexa", "Kaletra + Kivexa"))
lake2_2$sexo <- factor(lake2_2$sexo, levels = c(1,2), labels = c("Hombre", "Mujer"))
lake2_2$tiempo <- as.factor(lake2_2$tiempo)
lake2_2$status <- as.factor(lake2_2$status)
lake2_2$edad <- ifelse(lake2_2$edad == 0, NA, lake2_2$edad)Generamos un par de tablas con los datos más relevantes:
table1::table1( ~ sexo + edad + tpo_vih_meses | Grupo, lake2_2, overall = "Total")| EFV + Kivexa (N=56) |
Kaletra + Kivexa (N=56) |
Total (N=112) |
|
|---|---|---|---|
| sexo | |||
| Hombre | 46 (82.1%) | 45 (80.4%) | 91 (81.3%) |
| Mujer | 8 (14.3%) | 7 (12.5%) | 15 (13.4%) |
| Missing | 2 (3.6%) | 4 (7.1%) | 6 (5.4%) |
| edad | |||
| Mean (SD) | 38.5 (8.55) | 37.6 (8.13) | 38.1 (8.32) |
| Median [Min, Max] | 38.0 [21.0, 59.0] | 37.0 [20.0, 58.0] | 37.0 [20.0, 59.0] |
| Missing | 1 (1.8%) | 1 (1.8%) | 2 (1.8%) |
| tpo_vih_meses | |||
| Mean (SD) | 36.5 (62.1) | 20.5 (44.3) | 28.0 (53.7) |
| Median [Min, Max] | 13.9 [0.233, 285] | 4.32 [0.467, 214] | 6.05 [0.233, 285] |
| Missing | 12 (21.4%) | 6 (10.7%) | 18 (16.1%) |
table1::table1( ~ tiempo + status | Grupo, lake2_2, overall = "Total")| EFV + Kivexa (N=56) |
Kaletra + Kivexa (N=56) |
Total (N=112) |
|
|---|---|---|---|
| tiempo | |||
| 0 | 8 (14.3%) | 5 (8.9%) | 13 (11.6%) |
| 12 | 11 (19.6%) | 9 (16.1%) | 20 (17.9%) |
| 24 | 9 (16.1%) | 12 (21.4%) | 21 (18.8%) |
| 36 | 4 (7.1%) | 7 (12.5%) | 11 (9.8%) |
| 48 | 24 (42.9%) | 23 (41.1%) | 47 (42.0%) |
| status | |||
| 0 | 41 (73.2%) | 37 (66.1%) | 78 (69.6%) |
| 1 | 15 (26.8%) | 19 (33.9%) | 34 (30.4%) |
Y un gráfico con la evolución de la carga viral:
lake_cv50 <- lake[, c("cv50_0", "cv50_12", "cv50_24", "cv50_36", "cv50_48")]
lake_cv50 <- lapply(lake_cv50, factor)
lake_cv50$id <- 1:116
long_cv <- melt(setDT(lake_cv50), id.vars = c("id"), variable.name = "Carga_viral")
long_cv$value <- factor(long_cv$value, levels = 0:1,
labels = c("Superior al límite", "Indetectable"))
long_cv2 <- na.omit(long_cv) %>%
group_by(Carga_viral, value) %>%
summarise(count=n()) %>%
mutate(perc=count/sum(count))
ggplot(long_cv2, aes(x = Carga_viral, y = perc*100, fill = value)) +
geom_bar(stat = "identity", position = 'dodge') +
ylim(0,100) +
labs(y = "%", x = "Semana", fill = "Carga viral")
gg <- ggplot(long_cv2, aes(x = Carga_viral, y = perc*100, fill = value)) +
geom_bar(stat = "identity", position = 'dodge') +
ylim(0,100) +
labs(y = "%", x = "Semana", fill = "Carga viral")
# ggtoppt(gg)
#
# ggtoppt(export = TRUE, "graficos.pptx")Generamos unas curvas de Kaplan Meier para ver cómo funcionan los modelos tradicionales.
Hacemos una con todas las variables:
sur1 <- survfit(Surv(tiempo, status) ~ 1, data = dd_limpio)
ggsurvplot(sur1,
risk.table = T,
surv.median.line = "hv",
linetype = "strata",
break.time.by = 12,
xlab = "Semanas",
conf.int = TRUE,
title = "Todos")
Y otra por grupos:
lake_surv <- dd_limpio[, c("sexo", "edad", "Grupo", "tpo_vih_meses", "tiempo", "status")]
lake_surv$Grupo <- factor(lake_surv$Grupo, levels = c(-1,0), labels = c("EFV + Kivexa", "Kaletra + Kivexa"))
sur2 <- survfit(Surv(tiempo, status) ~ Grupo, data = lake_surv)
ggsurvplot(sur2,
risk.table = T,
surv.median.line = "hv",
pval = TRUE,
break.time.by = 12,
conf.int = TRUE,
xlab = "Semanas",
title = "")
Algoritmos de Machine Learning
Preparación de los datos
Comenzamos dividiendo la base de datos en dos partes: entrenamiento (67%) y test (33%).
dd_limpio_ml <- dd_limpio
dd_limpio_ml$tiempo <- ifelse(dd_limpio_ml$tiempo == 0 , 0.0, dd_limpio_ml$tiempo)
dd_limpio_ml$Grupo <- ifelse(dd_limpio_ml$Grupo == -1, 1, 0)
#Implantamos la semilla
set.seed(123)
#Separamos las dos partes creando un objeto que contenga el 67% de los datos.
objeto1 <- sample(nrow(dd_limpio), nrow(dd_limpio)*0.67)
#Al subset train le añadimos el 67% de los datos
train <- dd_limpio[objeto1, ]
#Al subset test le añadimos el resto (33%)
test <- dd_limpio[-objeto1, ]Árbol de supervivencia
Primero realizamos un modelo con 100 árboles.
train_forest <- train
test_forest <- test
rf_1 <- rfsrc(Surv(tiempo, status) ~ .,
data = train_forest,
ntree = 100,
na.action = "na.impute",
tree.err = TRUE,importance = TRUE)
rf_1## Sample size: 75
## Number of deaths: 24
## Was data imputed: yes
## Number of trees: 100
## Forest terminal node size: 15
## Average no. of terminal nodes: 4.71
## No. of variables tried at each split: 7
## Total no. of variables: 37
## Resampling used to grow trees: swor
## Resample size used to grow trees: 47
## Analysis: RSF
## Family: surv
## Splitting rule: logrank *random*
## Number of random split points: 10
## (OOB) CRPS: 0.19838121
## (OOB) Requested performance error: 0.60320934Representamos la estimación del error de predicción en función del número de árboles en el bosque:
plot(gg_error(rf_1)) 
Hacemos las predicciones y volvemos a representar gráficamente el error:
rf_test <- predict(rf_1, newdata = test_forest, na.action = "na.impute", importance = TRUE)
rf_test## Sample size of test (predict) data: 37
## Was test data imputed: yes
## Number of grow trees: 100
## Average no. of grow terminal nodes: 4.71
## Total no. of grow variables: 37
## Resampling used to grow trees: swor
## Resample size used to grow trees: 47
## Analysis: RSF
## Family: surv
## CRPS: 0.19586458
## Requested performance error: 0.36050725plot(gg_error(rf_test)) 
c_100rf <- compareC::estC(test_forest$tiempo, test_forest$status, rf_test$predicted)También podemos representar qué variables son las que más influyen en los resultados: en azul aparecen las variables que reducen el error y en rojo las que lo aumentan.
plot(gg_vimp(rf_test))
El índice C de Harrell para este modelo es: 0.3558
Repetimos el modelo pero aumentando el número de árboles:
rf_2 <- rfsrc(Surv(tiempo, status) ~ .,
data = train_forest,
ntree = 300,
na.action = "na.impute",
tree.err = TRUE,importance = TRUE)
rf_2## Sample size: 75
## Number of deaths: 24
## Was data imputed: yes
## Number of trees: 300
## Forest terminal node size: 15
## Average no. of terminal nodes: 4.6467
## No. of variables tried at each split: 7
## Total no. of variables: 37
## Resampling used to grow trees: swor
## Resample size used to grow trees: 47
## Analysis: RSF
## Family: surv
## Splitting rule: logrank *random*
## Number of random split points: 10
## (OOB) CRPS: 0.20195885
## (OOB) Requested performance error: 0.71699489plot(gg_error(rf_2))
rf_test2 <- predict(rf_2, newdata = test_forest, na.action = "na.impute", importance = TRUE)
rf_test2## Sample size of test (predict) data: 37
## Was test data imputed: yes
## Number of grow trees: 300
## Average no. of grow terminal nodes: 4.6467
## Total no. of grow variables: 37
## Resampling used to grow trees: swor
## Resample size used to grow trees: 47
## Analysis: RSF
## Family: surv
## CRPS: 0.19411631
## Requested performance error: 0.35326087plot(gg_error(rf_test2))
plot(gg_vimp(rf_test2))
c_300rf <- compareC::estC(test_forest$tiempo, test_forest$status, rf_test2$predicted)Ahora el índice C es 0.3483.
Naive Bayes
De este algoritmo también vamos a hacer dos versiones: una sin la corrección de Laplace y otra con ella.
Laplace = 0
bayes_lap0_model<-naiveBayes(Surv(tiempo,status) ~ ., data=train, laplace=0, na.action = na.pass)
summary(bayes_lap0_model)## Length Class Mode
## apriori 9 table numeric
## tables 37 -none- list
## levels 9 -none- character
## isnumeric 37 -none- logical
## call 4 -none- callbayes_lap0_pred <- predict(bayes_lap0_model, test)
train_labels <- Surv(train$tiempo, train$status)
test_labels <- as.factor(Surv(test$tiempo,test$status))
c_nb0 <- compareC::estC(test$tiempo, test$status, bayes_lap0_pred)Obtenemos un índice C de 0.6067.
Laplace 1
bayes_lap1_model <- naiveBayes(Surv(tiempo,status) ~ ., data=train, laplace=1, na.action = na.pass)
summary(bayes_lap1_model)## Length Class Mode
## apriori 9 table numeric
## tables 37 -none- list
## levels 9 -none- character
## isnumeric 37 -none- logical
## call 4 -none- callbayes_lap1_pred <- predict(bayes_lap1_model, test)
c_nb1 <- compareC::estC(test$tiempo, test$status, bayes_lap1_pred)Obtenemos un índice C de 0.6067.
Support Vector Machine
En el caso del SVM, debemos eliminar los NAs para que no den problemas con los modelos.
#Implantamos la semilla
set.seed(123)
dd_svm <- na.omit(dd_limpio_ml)
objeto_svm <- sample(nrow(dd_svm), nrow(dd_svm)*0.67)
train_svm <- dd_svm[objeto_svm,]
test_svm <- dd_svm[-objeto_svm,]Vamos a probar tres variantes:
Kernel lineal
surv_lin <- survivalsvm(Surv(tiempo, status) ~.,
data = train_svm,
type = "regression", gamma.mu = 1,
kernel = "lin_kernel")
print(surv_lin)##
## survivalsvm result
##
## Call:
##
## survivalsvm(Surv(tiempo, status) ~ ., data = train_svm, type = "regression", gamma.mu = 1, kernel = "lin_kernel")
##
## Survival svm approach : regression
## Type of Kernel : lin_kernel
## Optimization solver used : quadprog
## Number of support vectors retained : 3
## survivalsvm version : 0.0.5pred_surv_lin <- predict(surv_lin, newdata = test_svm)
print(pred_surv_lin)##
##
## survivalsvm prediction
##
## Type of survivalsvm : regression
## Type of kernel : lin_kernel
## Optimization solver used in model : quadprog
## predicted risk ranks : 16.69 -9.29 -11.73 19.65 34.58 ...
## survivalsvm version : 0.0.5c_svmlin <- compareC::estC(test_svm$tiempo, test_svm$status, pred_surv_lin$predicted)Obtenemos una C de 0.6176
Kernel radial
surv_rbf <- survivalsvm(Surv(tiempo, status) ~.,
data = train_svm,
type = "regression", gamma.mu = 1,
opt.meth = "quadprog", kernel = "rbf_kernel")
print(surv_rbf)##
## survivalsvm result
##
## Call:
##
## survivalsvm(Surv(tiempo, status) ~ ., data = train_svm, type = "regression", gamma.mu = 1, opt.meth = "quadprog", kernel = "rbf_kernel")
##
## Survival svm approach : regression
## Type of Kernel : rbf_kernel
## Optimization solver used : quadprog
## Number of support vectors retained : 25
## survivalsvm version : 0.0.5pred_surv_rbf <- predict(surv_rbf, test_svm)
print(pred_surv_rbf)##
##
## survivalsvm prediction
##
## Type of survivalsvm : regression
## Type of kernel : rbf_kernel
## Optimization solver used in model : quadprog
## predicted risk ranks : 40.8 40.8 40.8 40.8 40.8 ...
## survivalsvm version : 0.0.5c_svmrbf <- compareC::estC(test_svm$tiempo, test_svm$status, pred_surv_rbf$predicted)Obtenemos una C de 0.5
Kernel aditivo
surv_add <- survivalsvm(Surv(tiempo, status) ~.,
data = train_svm[,2:39],
type = "regression", gamma.mu = 1,
opt.meth = "quadprog", kernel = "add_kernel")
print(surv_add)##
## survivalsvm result
##
## Call:
##
## survivalsvm(Surv(tiempo, status) ~ ., data = train_svm[, 2:39], type = "regression", gamma.mu = 1, opt.meth = "quadprog", kernel = "add_kernel")
##
## Survival svm approach : regression
## Type of Kernel : add_kernel
## Optimization solver used : quadprog
## Number of support vectors retained : 8
## survivalsvm version : 0.0.5pred_surv_add <- predict(surv_add, newdata = test_svm)
print(pred_surv_add)##
##
## survivalsvm prediction
##
## Type of survivalsvm : regression
## Type of kernel : add_kernel
## Optimization solver used in model : quadprog
## predicted risk ranks : 36.92 34.33 38.96 35.75 35.53 ...
## survivalsvm version : 0.0.5conindex(pred_surv_add, test_svm$tiempo)## C Index
## 0.4095238c_svmadd <- compareC::estC(test_svm$tiempo, test_svm$status, pred_surv_add$predicted)Obtenemos una C de 0.1765
Resumen de los resultados
En la siguiente tabla se muestran los valores de C para los diferentes algoritmos:
c_values <- data.frame(names = c("Árbol de Supervivencia n = 100", "Árbol de supervivencia n = 300",
"Naive Bayes L = 0", "Naive Bayes L = 1",
"SVM Lineal", "SVM RBF", "SVM Aditivo"),
values = c(c_100rf, c_300rf, c_nb0, c_nb1, c_svmlin, c_svmrbf, c_svmadd))
kable(c_values, caption = "Resumen de los valores de C para todos los algoritmos",
col.names = c("Algoritmo", "C-Index"),
digits = 4,
format = "html", table.attr = "style='width:100%;'", align = "c") %>%
kable_styling(bootstrap_options = c("striped", "hover","condensed"))| Algoritmo | C-Index |
|---|---|
| Árbol de Supervivencia n = 100 | 0.3558 |
| Árbol de supervivencia n = 300 | 0.3483 |
| Naive Bayes L = 0 | 0.6067 |
| Naive Bayes L = 1 | 0.6067 |
| SVM Lineal | 0.6176 |
| SVM RBF | 0.5000 |
| SVM Aditivo | 0.1765 |
max_c <- max(c_values$values)El algoritmo con mayor valor de C es SVM Lineal, con un índice C de 0.6176.