Optimización campaña sanitaria
Data Science: Fernando Cabello Rodriguez
24/11/2020
#CONTEXTO
- Se necesita saber que pacientes son los que presentan mayor probabilidad de muerte, conforme a, variables obtenidas de realizar determinados análisis de diagnósticos clínico. De esta forma, se podrá saber a que pacientes dirigirse primero, para realizar una campaña de prevención de aparición de eventos que pueden causar la muerte.
DETALLE DEL PROYECTO REALIZADO
Se ha trabajado sobre un histórico de parámetros o variables de laboratorio y se ha obtenido un modelo predictivo de alta calidad.
El modelo es estable y consigue separar satisfactoriamente a los clientes por su mayor probabilidad de muerte
- Preparamos el entorno. 1.1 Descargamos las librerías con las que vamos atrabajr
paquetes <- c('data.table',#para leer y escribir datos de forma rapida
'dplyr',#para manipulación de datos
'tidyr',#para manipulación de datos
'ggplot2',#para gráficos
'randomForest',#para crear los modelos
'ROCR',#para evaluar modelos
'purrr',#para usar la función map que aplica la misma funciona a varios componentes de un dataframe
'smbinning',#para calcular la para importancia de las variables
'rpart',#para crear arboles de decisión
'rpart.plot'#para el gráfico del árbol
)
#Crea un vector lógico con si están instalados o no
instalados <- paquetes %in% installed.packages()
#Si hay al menos uno no instalado los instala
if(sum(instalados == FALSE) > 0) {
install.packages(paquetes[!instalados])
}
lapply(paquetes,require,character.only = TRUE)## Loading required package: data.table## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:data.table':
##
## between, first, last## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union## Loading required package: tidyr## Loading required package: ggplot2## Loading required package: randomForest## randomForest 4.6-14## Type rfNews() to see new features/changes/bug fixes.##
## Attaching package: 'randomForest'## The following object is masked from 'package:ggplot2':
##
## margin## The following object is masked from 'package:dplyr':
##
## combine## Loading required package: ROCR## Loading required package: purrr##
## Attaching package: 'purrr'## The following object is masked from 'package:data.table':
##
## transpose## Loading required package: smbinning## Loading required package: sqldf## Loading required package: gsubfn## Loading required package: proto## Loading required package: RSQLite## Loading required package: partykit## Loading required package: grid## Loading required package: libcoin## Loading required package: mvtnorm## Loading required package: Formula## Loading required package: rpart## Loading required package: rpart.plot## [[1]]
## [1] TRUE
##
## [[2]]
## [1] TRUE
##
## [[3]]
## [1] TRUE
##
## [[4]]
## [1] TRUE
##
## [[5]]
## [1] TRUE
##
## [[6]]
## [1] TRUE
##
## [[7]]
## [1] TRUE
##
## [[8]]
## [1] TRUE
##
## [[9]]
## [1] TRUE
##
## [[10]]
## [1] TRUE1.2 Cargamos el dataset
df <- fread('heart_failure_clinical_records_dataset.csv')- Análisis exploratorio 2.1 Análisis exploratorio general y tipo de datos
glimpse(df)## Rows: 299
## Columns: 13
## $ age <dbl> 75, 55, 65, 50, 65, 90, 75, 60, 65, 80, 75...
## $ anaemia <int> 0, 0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, ...
## $ creatinine_phosphokinase <int> 582, 7861, 146, 111, 160, 47, 246, 315, 15...
## $ diabetes <int> 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, ...
## $ ejection_fraction <int> 20, 38, 20, 20, 20, 40, 15, 60, 65, 35, 38...
## $ high_blood_pressure <int> 1, 0, 0, 0, 0, 1, 0, 0, 0, 1, 1, 1, 0, 1, ...
## $ platelets <dbl> 265000, 263358, 162000, 210000, 327000, 20...
## $ serum_creatinine <dbl> 1.90, 1.10, 1.30, 1.90, 2.70, 2.10, 1.20, ...
## $ serum_sodium <int> 130, 136, 129, 137, 116, 132, 137, 131, 13...
## $ sex <int> 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, ...
## $ smoking <int> 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 1, 0, 0, ...
## $ time <int> 4, 6, 7, 7, 8, 8, 10, 10, 10, 10, 10, 10, ...
## $ DEATH_EVENT <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ...Conclusiones: -Tenemos un dataframe con 13 variables y 299 registros relacionados con la presentcia de infartos cardiovasculares. - Tenemos variables que se deben pasar a factores, anemia, diabetes, high_blood_pressure, sex, smoking, DEATH_EVENT. - Se debe predecir la probabilidad de una futura presencia de DEATH_EVENT en los pacientes que no han sufrido muerte por infarto, de acuerdo a estos datos de laboratorio.
Reservamos las variables a pasar a factores
a_factores <- c('anaemia', 'diabetes', 'high_blood_pressure', 'sex', 'smoking', 'DEATH_EVENT')2.2 Calidad de datos: Estadísticos básicos Hacemos un summary, con lapply que sale en formato de lista y se lee mejor
lapply(df, summary)## $age
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 40.00 51.00 60.00 60.83 70.00 95.00
##
## $anaemia
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4314 1.0000 1.0000
##
## $creatinine_phosphokinase
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 23.0 116.5 250.0 581.8 582.0 7861.0
##
## $diabetes
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.4181 1.0000 1.0000
##
## $ejection_fraction
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 14.00 30.00 38.00 38.08 45.00 80.00
##
## $high_blood_pressure
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.3512 1.0000 1.0000
##
## $platelets
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 25100 212500 262000 263358 303500 850000
##
## $serum_creatinine
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.500 0.900 1.100 1.394 1.400 9.400
##
## $serum_sodium
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 113.0 134.0 137.0 136.6 140.0 148.0
##
## $sex
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 1.0000 0.6488 1.0000 1.0000
##
## $smoking
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.3211 1.0000 1.0000
##
## $time
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 4.0 73.0 115.0 130.3 203.0 285.0
##
## $DEATH_EVENT
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.0000 0.0000 0.0000 0.3211 1.0000 1.0000Time -> tiene un máximo de 285, lo apuntamos para su posterior revisión.
2.3 Calidad de datos: Análisi de nulos
colSums(is.na(df))## age anaemia creatinine_phosphokinase
## 0 0 0
## diabetes ejection_fraction high_blood_pressure
## 0 0 0
## platelets serum_creatinine serum_sodium
## 0 0 0
## sex smoking time
## 0 0 0
## DEATH_EVENT
## 0Conclusión: No hay nulos, pero sospechamos que pueden estar reconvertidos a ceros. 2.4 - Calidad de datos: Análisis de ceros
No es algo que se haga siempre, pero en el análisis general superior habíamos detectado muchos ceros. Vamos a construir una función concreta para analizar esto
contar_ceros <- function(variable){
temp <- transmute(df, if_else(variable==0,1,0))
sum(temp)
}
num_ceros <- sapply(df, contar_ceros)
num_ceros <- data.frame(VARIABLE=names(num_ceros), CEROS=as.numeric(num_ceros), stringsAsFactors = F)
num_ceros <- num_ceros%>%
arrange(desc(CEROS))%>%
mutate(PORCENTAJE = CEROS / nrow(df)*100)Visualizamos la variable creada
num_ceros## VARIABLE CEROS PORCENTAJE
## 1 smoking 203 67.89298
## 2 DEATH_EVENT 203 67.89298
## 3 high_blood_pressure 194 64.88294
## 4 diabetes 174 58.19398
## 5 anaemia 170 56.85619
## 6 sex 105 35.11706
## 7 age 0 0.00000
## 8 creatinine_phosphokinase 0 0.00000
## 9 ejection_fraction 0 0.00000
## 10 platelets 0 0.00000
## 11 serum_creatinine 0 0.00000
## 12 serum_sodium 0 0.00000
## 13 time 0 0.00000Todas la variables tiene información, asi que el número de ceros no influye en nuestra calidad de los datos, debido a que hay un gran orcentaje de información en cada una de nuestras variables.
2.5 - Calidad de datos: Análisis de atípicos
2.5.1 - Analizamos las que son de tipo numérico
out <- function(variable){
t(t(head(sort(variable,decreasing = T),20)))
}
lapply(df,function(x){
if(is.double(x)) out(x)
})## $age
## [,1]
## [1,] 95
## [2,] 95
## [3,] 94
## [4,] 90
## [5,] 90
## [6,] 90
## [7,] 87
## [8,] 86
## [9,] 85
## [10,] 85
## [11,] 85
## [12,] 85
## [13,] 85
## [14,] 85
## [15,] 82
## [16,] 82
## [17,] 82
## [18,] 81
## [19,] 80
## [20,] 80
##
## $anaemia
## NULL
##
## $creatinine_phosphokinase
## NULL
##
## $diabetes
## NULL
##
## $ejection_fraction
## NULL
##
## $high_blood_pressure
## NULL
##
## $platelets
## [,1]
## [1,] 850000
## [2,] 742000
## [3,] 621000
## [4,] 543000
## [5,] 533000
## [6,] 507000
## [7,] 504000
## [8,] 497000
## [9,] 481000
## [10,] 461000
## [11,] 454000
## [12,] 451000
## [13,] 451000
## [14,] 448000
## [15,] 427000
## [16,] 422000
## [17,] 418000
## [18,] 406000
## [19,] 406000
## [20,] 404000
##
## $serum_creatinine
## [,1]
## [1,] 9.4
## [2,] 9.0
## [3,] 6.8
## [4,] 6.1
## [5,] 5.8
## [6,] 5.0
## [7,] 4.4
## [8,] 4.0
## [9,] 3.8
## [10,] 3.7
## [11,] 3.5
## [12,] 3.5
## [13,] 3.4
## [14,] 3.2
## [15,] 3.0
## [16,] 3.0
## [17,] 2.9
## [18,] 2.7
## [19,] 2.7
## [20,] 2.7
##
## $serum_sodium
## NULL
##
## $sex
## NULL
##
## $smoking
## NULL
##
## $time
## NULL
##
## $DEATH_EVENT
## NULL2.5.2 - Analizamos las que son de tipo integer
out <- function(variable){
t(t(table(variable))) #la doble traspuesta es un truco para que se visualice la salida, si no lo que crearìa es una colección de dataframes que no se ven bien
}
lapply(df,function(x){
if(is.integer(x)) out(x)
})## $age
## NULL
##
## $anaemia
##
## variable [,1]
## 0 170
## 1 129
##
## $creatinine_phosphokinase
##
## variable [,1]
## 23 1
## 30 1
## 47 3
## 52 1
## 53 1
## 54 1
## 55 1
## 56 2
## 57 1
## 58 1
## 59 3
## 60 3
## 61 2
## 62 1
## 63 1
## 64 3
## 66 4
## 68 3
## 69 3
## 70 1
## 72 1
## 75 1
## 76 1
## 78 1
## 80 2
## 81 2
## 84 3
## 86 1
## 88 1
## 90 1
## 91 1
## 92 1
## 93 1
## 94 1
## 95 1
## 96 2
## 97 1
## 99 1
## 101 1
## 102 2
## 103 1
## 104 1
## 109 2
## 110 1
## 111 1
## 112 1
## 113 2
## 115 3
## 118 1
## 119 1
## 121 1
## 122 2
## 123 1
## 124 1
## 125 1
## 127 1
## 128 1
## 129 4
## 130 1
## 131 1
## 132 2
## 133 1
## 135 2
## 143 2
## 144 1
## 145 1
## 146 1
## 148 2
## 149 1
## 151 1
## 154 1
## 156 1
## 157 2
## 159 1
## 160 1
## 161 1
## 166 1
## 167 2
## 168 2
## 170 1
## 171 1
## 176 1
## 180 1
## 185 1
## 190 1
## 191 1
## 193 1
## 196 2
## 198 1
## 200 1
## 203 1
## 207 1
## 211 1
## 212 2
## 213 1
## 220 1
## 224 2
## 231 3
## 232 1
## 233 1
## 235 1
## 244 1
## 245 1
## 246 1
## 248 1
## 249 1
## 250 2
## 253 1
## 257 1
## 258 1
## 260 1
## 270 1
## 280 1
## 281 1
## 291 1
## 292 1
## 298 1
## 305 1
## 308 1
## 315 1
## 318 1
## 320 1
## 326 1
## 328 1
## 335 1
## 336 1
## 337 1
## 358 1
## 364 1
## 369 1
## 371 1
## 379 1
## 395 1
## 400 1
## 418 1
## 427 1
## 446 1
## 478 1
## 482 1
## 514 1
## 553 1
## 571 1
## 572 1
## 577 1
## 582 47
## 588 1
## 607 1
## 615 1
## 618 1
## 624 1
## 646 1
## 655 1
## 675 1
## 707 1
## 719 1
## 720 1
## 737 1
## 748 1
## 754 1
## 776 1
## 789 1
## 805 1
## 835 2
## 855 1
## 892 1
## 897 1
## 898 1
## 910 1
## 936 1
## 943 1
## 972 1
## 981 1
## 1021 1
## 1051 1
## 1082 1
## 1185 1
## 1199 1
## 1202 1
## 1211 1
## 1380 1
## 1419 1
## 1548 1
## 1610 1
## 1688 1
## 1767 1
## 1808 1
## 1820 1
## 1846 1
## 1876 1
## 1896 1
## 2017 1
## 2060 1
## 2261 1
## 2281 1
## 2334 1
## 2413 1
## 2442 1
## 2522 1
## 2656 1
## 2695 1
## 2794 1
## 3964 1
## 3966 1
## 4540 1
## 5209 1
## 5882 1
## 7702 1
## 7861 1
##
## $diabetes
##
## variable [,1]
## 0 174
## 1 125
##
## $ejection_fraction
##
## variable [,1]
## 14 1
## 15 2
## 17 2
## 20 18
## 25 36
## 30 34
## 35 49
## 38 40
## 40 37
## 45 20
## 50 21
## 55 3
## 60 31
## 62 2
## 65 1
## 70 1
## 80 1
##
## $high_blood_pressure
##
## variable [,1]
## 0 194
## 1 105
##
## $platelets
## NULL
##
## $serum_creatinine
## NULL
##
## $serum_sodium
##
## variable [,1]
## 113 1
## 116 1
## 121 1
## 124 1
## 125 1
## 126 1
## 127 3
## 128 2
## 129 2
## 130 9
## 131 5
## 132 14
## 133 10
## 134 32
## 135 16
## 136 40
## 137 38
## 138 23
## 139 22
## 140 35
## 141 12
## 142 11
## 143 3
## 144 5
## 145 9
## 146 1
## 148 1
##
## $sex
##
## variable [,1]
## 0 105
## 1 194
##
## $smoking
##
## variable [,1]
## 0 203
## 1 96
##
## $time
##
## variable [,1]
## 4 1
## 6 1
## 7 2
## 8 2
## 10 6
## 11 2
## 12 1
## 13 1
## 14 2
## 15 2
## 16 1
## 20 2
## 22 1
## 23 2
## 24 1
## 26 3
## 27 1
## 28 2
## 29 2
## 30 5
## 31 1
## 32 1
## 33 3
## 35 1
## 38 1
## 40 1
## 41 1
## 42 1
## 43 3
## 44 1
## 45 1
## 50 1
## 54 2
## 55 1
## 59 1
## 60 3
## 61 1
## 63 1
## 64 1
## 65 2
## 66 1
## 67 1
## 68 1
## 71 1
## 72 2
## 73 2
## 74 4
## 75 1
## 76 1
## 77 1
## 78 2
## 79 5
## 80 2
## 82 2
## 83 3
## 85 2
## 86 1
## 87 5
## 88 5
## 90 4
## 91 2
## 94 3
## 95 5
## 96 1
## 97 1
## 100 1
## 104 2
## 105 1
## 106 1
## 107 6
## 108 3
## 109 3
## 110 1
## 111 1
## 112 2
## 113 2
## 115 2
## 117 1
## 118 1
## 119 1
## 120 4
## 121 4
## 123 1
## 126 1
## 129 1
## 130 1
## 134 1
## 135 1
## 140 1
## 145 2
## 146 5
## 147 4
## 148 1
## 150 1
## 154 1
## 162 1
## 170 1
## 171 1
## 172 3
## 174 3
## 175 1
## 180 3
## 185 1
## 186 6
## 187 7
## 188 1
## 192 2
## 193 1
## 194 1
## 195 1
## 196 2
## 197 2
## 198 1
## 200 1
## 201 2
## 205 3
## 206 1
## 207 3
## 208 1
## 209 5
## 210 2
## 211 1
## 212 3
## 213 3
## 214 5
## 215 4
## 216 1
## 220 1
## 230 2
## 231 1
## 233 2
## 235 1
## 237 2
## 240 1
## 241 1
## 244 5
## 245 5
## 246 3
## 247 1
## 250 7
## 256 2
## 257 1
## 258 2
## 270 2
## 271 1
## 278 1
## 280 1
## 285 1
##
## $DEATH_EVENT
##
## variable [,1]
## 0 203
## 1 96Conclusiones: - Se puede observar 47 pacientes con un valor de 587 u/l de creatinine_phosphokinase. - Los valores de eyección se concentran entre 20-50, algo reducido, con respecto a los valores normales de esta variable.
2.6 - Análisis longitudinal
longi <- df%>%
summarise_all(mean)%>% # Calular la media de cada variable
t()%>% #trasponerlo para tenerlo en una sola columna y leerlo mejor
as.data.frame()#reconvertirlo a dataframe porque t() lo pasa a matriz
data.frame(variable=rownames(longi), media= longi$V1)%>%
arrange(desc(variable)) #ordenar por el nombre para tener la visión longitudinal## variable media
## 1 time 1.302609e+02
## 2 smoking 3.210702e-01
## 3 sex 6.488294e-01
## 4 serum_sodium 1.366254e+02
## 5 serum_creatinine 1.393880e+00
## 6 platelets 2.633580e+05
## 7 high_blood_pressure 3.511706e-01
## 8 ejection_fraction 3.808361e+01
## 9 diabetes 4.180602e-01
## 10 DEATH_EVENT 3.210702e-01
## 11 creatinine_phosphokinase 5.818395e+02
## 12 anaemia 4.314381e-01
## 13 age 6.083389e+01longi ## V1
## age 6.083389e+01
## anaemia 4.314381e-01
## creatinine_phosphokinase 5.818395e+02
## diabetes 4.180602e-01
## ejection_fraction 3.808361e+01
## high_blood_pressure 3.511706e-01
## platelets 2.633580e+05
## serum_creatinine 1.393880e+00
## serum_sodium 1.366254e+02
## sex 6.488294e-01
## smoking 3.210702e-01
## time 1.302609e+02
## DEATH_EVENT 3.210702e-01Conclusiones: todos los datos parecen normales. Vamos a comprobar los 47 pacientes con un valor de 587 u/l de creatinine_phosphokinase.
hist(df$creatinine_phosphokinase,breaks = 200)
Vemos una frecuencia desproporcionada de clientes con valores entre 50 y 400, lo cual es normal porque los valores normales oscilan aproximadamente entre 30-400. Por otro lado los valores altos indican mayor porbabilidad de infarto de miocardio.
hist(df$ejection_fraction, breaks = 20)
Se puede observar que hay una frecuencia de aproximadamente un 80% de los valores totales de fracción de eyección con valores comprendidos entre 35-40, presentando una fracción de eyección algo reducida. 2.9 - Acciones resultado del analisis de calidad de datos y exploratorio
df <-df %>%
mutate_at(a_factores,as.factor)Grafico la cantidad de fallecidos
ggplot(data= df, aes(x = DEATH_EVENT))+ geom_bar()+
xlab('Fallecimientos')+
ylab('Pacientes')+
ggtitle('Histograma de fallecimientos')
Se puede visualizar la cantidad de pacientes fallecidos (1) y no fallecidos(0) 3- Creación de variables sintéticas. Al ser pocas variables no se hace preselección de las variables más predictivas, pasando directamente a la creación de variables sintéticas. En concreto solo se va a discretizar las variables age, creatinine_phosphokinase, ejection_fraction, platelets, serum_creatinine , serum_sodium. 3.1 Discretizar. Se crea una función para discretizar las variables selecionadas.
discretizar <- function(vi, target){
temp_df <- data.frame(vi=vi, target=target)
temp_df$target <- as.numeric(as.character(temp_df$target))
disc <- smbinning(temp_df, y = 'target', x = 'vi')
return(disc)
}Discretizamos age
disc_temp_EDAD <- discretizar(df$age, df$DEATH_EVENT)
df_temp <- select(df, age, DEATH_EVENT)
df_temp <- as.data.frame(df_temp)
df_temp <- smbinning.gen(df_temp, disc_temp_EDAD, chrname = 'EDAD_DISC')
df <- cbind(df,df_temp[3]) %>% select(-age)Discretizamos de forma manual creatinine_phosphokinase
ggplot(df,aes(creatinine_phosphokinase)) + geom_density() + scale_x_continuous(limits = c(0, 8000))
Parece que la gran mayoría esta por debajo de 2000
ggplot(df,aes(creatinine_phosphokinase)) + geom_density()+scale_x_continuous(limits = c(0, 2000))## Warning: Removed 18 rows containing non-finite values (stat_density).
Realizamos un histograma para aproximar mejor la forma que queremos conseguir
ggplot(df, aes(creatinine_phosphokinase))+geom_histogram(bins = 10)+
scale_x_continuous(limits = c(0, 4000)) ## Warning: Removed 5 rows containing non-finite values (stat_bin).## Warning: Removed 2 rows containing missing values (geom_bar).
Ya sabemos que queremos una forma decreciente ahora veamos como se comporta la variable target para ver si podremos generar algo monotónico
ggplot(df,aes(creatinine_phosphokinase, fill=DEATH_EVENT))+ geom_histogram(bins = 30, position = 'fill') + scale_x_continuous(limits = c(0, 5000))## Warning: Removed 4 rows containing non-finite values (stat_bin).## Warning: Removed 24 rows containing missing values (geom_bar).
Sabiendo ambas cosas vamos a apoyarnos en los deciles para intuir donde podemos hacer buenos cortes
as.data.frame(quantile(df$creatinine_phosphokinase,prob = seq(0, 1, length = 11)))## quantile(df$creatinine_phosphokinase, prob = seq(0, 1, length = 11))
## 0% 23.0
## 10% 67.6
## 20% 100.2
## 30% 130.4
## 40% 176.8
## 50% 250.0
## 60% 425.2
## 70% 582.0
## 80% 620.4
## 90% 1203.8
## 100% 7861.0df <- df %>% mutate(creatinine_phosphokinase_DISC = as.factor(case_when(
creatinine_phosphokinase <= 100 ~ '01_MENOR_100',
creatinine_phosphokinase >= 100 & creatinine_phosphokinase <= 200 ~'02_DE_100_A_200',
creatinine_phosphokinase >= 200 & creatinine_phosphokinase <= 600 ~'03_DE_200_A_600',
creatinine_phosphokinase >= 600 & creatinine_phosphokinase <= 2000 ~'04_DE_600_A_2000',
creatinine_phosphokinase >=2000 ~ '05_MAYOR_2000',
TRUE ~'00_ERROR'))
)Borramos la variable original
df <- select(df,-creatinine_phosphokinase)Discretizamos manualmente ejection_fraction
ggplot(df, aes(ejection_fraction))+geom_density()+scale_x_continuous(limits = c(0, 100))
ggplot(df, aes(ejection_fraction))+geom_histogram(bins = 10)+scale_x_continuous(limits = c(0, 80))## Warning: Removed 2 rows containing missing values (geom_bar).
Esta variable no es monotónica.
Vamos como se comporta la variable target para ver si podremos generar algo monotónico
ggplot(df,aes(ejection_fraction, fill=DEATH_EVENT))+ geom_histogram(bins = 10, position = 'fill') + scale_x_continuous(limits = c(0, 80))## Warning: Removed 6 rows containing missing values (geom_bar).
Vamos a apoyarnos en deciles para intuir donde podemos hacer buenos cortes
as.data.frame(quantile(df$ejection_fraction, prob = seq(0,1, length=11)))## quantile(df$ejection_fraction, prob = seq(0, 1, length = 11))
## 0% 14
## 10% 25
## 20% 30
## 30% 30
## 40% 35
## 50% 38
## 60% 38
## 70% 40
## 80% 47
## 90% 60
## 100% 80df <- df %>% mutate(ejection_fraction_DISC= as.factor(case_when(
ejection_fraction <= 30 ~ '01_MENOR_30',
ejection_fraction >= 30 & ejection_fraction <= 60 ~'02_DE_30_A_60',
ejection_fraction >= 60 ~'03_MAYOR_DE_60',
TRUE ~ '00_ERROR'))
)Borramos la variable
df <- select(df, -ejection_fraction)Probamos a discretizar de forma automática la variable serum_creatinine, pero primero vamos a visualizar como se comporta la variable con respecto a la target
ggplot(df, aes(serum_creatinine, color = DEATH_EVENT))+geom_density()
Discretizamos de manera manual
ggplot(df, aes(serum_creatinine))+geom_density()
ggplot(df, aes(serum_creatinine))+geom_density()+scale_x_continuous(limits = c(0, 4))## Warning: Removed 7 rows containing non-finite values (stat_density).
Vamos a determinar más en profundidad si tiene o no forma monotónica
ggplot(df, aes(serum_creatinine)) + geom_histogram(bins=10)+scale_x_continuous(limits = c(0,9))## Warning: Removed 1 rows containing non-finite values (stat_bin).## Warning: Removed 2 rows containing missing values (geom_bar).
Vemos que tiene forma decreciente.
El siguiente paso es determinar como se comporta la variable target para ver si podremos generar algo monotónico
ggplot(df, aes(serum_creatinine, fill= DEATH_EVENT))+geom_histogram(bins = 10, position = 'fill')+ scale_x_continuous(limits = c(0, 5))## Warning: Removed 5 rows containing non-finite values (stat_bin).## Warning: Removed 4 rows containing missing values (geom_bar).
Vemos que se comporta de forma monotónica. Ahora vamos a a apoyarnos en los deciles para intuir donde podemos hacer buenos cortes y discretizamos
as.data.frame(quantile(df$serum_creatinine, probs = seq(0,1, length=11)))## quantile(df$serum_creatinine, probs = seq(0, 1, length = 11))
## 0% 0.5
## 10% 0.8
## 20% 0.9
## 30% 1.0
## 40% 1.0
## 50% 1.1
## 60% 1.2
## 70% 1.3
## 80% 1.7
## 90% 2.1
## 100% 9.4df <- df%>% mutate(serum_creatinine_DISC=as.factor(case_when(
serum_creatinine <=1.0~'01_MENOR_1.0',
serum_creatinine > 1.0 & serum_creatinine <= 2.0~'02_DE_1.0_A_2.0',
serum_creatinine > 2.0 ~'03_MAYOR_DE_02',
TRUE~'00_ERROR'))
)Borramos la variable
df <- select(df, -serum_creatinine)Ahora vamos a discretizar platelets.
Para elo vamos a ver como se comporta la variable
ggplot(df, aes(platelets))+geom_density()
Acotamos el valor de la variable
ggplot(df, aes(platelets))+geom_density()+scale_x_continuous(limits= c(0,500000))## Warning: Removed 7 rows containing non-finite values (stat_density).
Hacemos un histograma para ver si se compota de manera descencdente
ggplot(df, aes(platelets))+geom_histogram(bins = 20)
Vamos a comprobar como se comporta con respecto a la target (DEATH_EVENT)
ggplot(df, aes(platelets, color = DEATH_EVENT))+geom_density()
Parece que la variable apenas no discrimina asi que no vamos a usar esta variable para para la predición. Borramos la variable
df <- select(df, -platelets)Vamos adiscretizar serum_sodium, pero previamente vamos a comprobar como se comporta con respecto a la target
ggplot(df, aes(serum_sodium, color= DEATH_EVENT))+ geom_density()
Vemos que discrimina (paraece haber una distribución entre difente entre los que mueren y no), por lo que vamos a discretizar
ggplot(df, aes(serum_sodium))+geom_density()
Vamos a limitar el eje x
ggplot(df, aes(serum_sodium))+geom_density()+ scale_x_continuous(limits=c(120, 150))## Warning: Removed 2 rows containing non-finite values (stat_density).
ggplot(df, aes(serum_sodium))+geom_histogram(bins = 10)+ scale_x_continuous(limits=c(120, 150))## Warning: Removed 2 rows containing non-finite values (stat_bin).## Warning: Removed 2 rows containing missing values (geom_bar).
as.data.frame(quantile(df$serum_sodium, probs = seq(0,1,length=11)))## quantile(df$serum_sodium, probs = seq(0, 1, length = 11))
## 0% 113.0
## 10% 132.0
## 20% 134.0
## 30% 135.0
## 40% 136.0
## 50% 137.0
## 60% 138.0
## 70% 139.0
## 80% 140.0
## 90% 141.2
## 100% 148.0Discretizamos serum_sodium
df <- df%>% mutate(serum_sodium_DISC=as.factor(case_when(
serum_sodium <=115~'01_MENOR_1.0',
serum_sodium > 115 & serum_sodium <= 140~'02_DE_1.0_A_2.0',
serum_sodium >140 ~'03_MAYOR_DE_02',
TRUE~'00_ERROR'))
)Borramos la variable origina
df <- select(df, -serum_sodium)Discretizamos de manera automática la variable time.
disc_temp_time <- discretizar(df$time, df$DEATH_EVENT)
df_temp <- select(df, time, DEATH_EVENT)
df_temp <- as.data.frame(df_temp)
df_temp <- smbinning.gen(df_temp, disc_temp_time, chrname = 'time_DISC')
df <- cbind(df,df_temp[3]) %>% select(-time)Vamos a hacer una inspeccion visual de todas las variables a ver si han salido bien
df %>%
select_if(is.factor)%>%
gather()%>%
ggplot(aes(value)) + geom_bar() +
facet_wrap(~key, scales = 'free') + theme(axis.text = element_text(size=4))## Warning: attributes are not identical across measure variables;
## they will be dropped
Ahora vamos a analizar la penetración de la target en cada categoría para ver si las variables han salido monotónicas
a <- function(var1, var2){
df_temp <- data.frame(var1=df[[var1]], var2 = df[[var2]])
df_temp %>%
group_by(var1) %>%
summarise(Conteo = n (), Porc = mean(as.numeric(as.character(var2))))%>%
ggplot(aes(var1, Porc))+geom_bar(stat = 'identity')+xlab(var1)
}
df2_nombres <- df%>%select_if(is.factor) %>% names()
lapply(df2_nombres, function(x){a(x, 'DEATH_EVENT')}) ## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)
## `summarise()` ungrouping output (override with `.groups` argument)## [[1]]
##
## [[2]]
##
## [[3]]
##
## [[4]]
##
## [[5]]
##
## [[6]]
##
## [[7]]
##
## [[8]]
##
## [[9]]
##
## [[10]]
##
## [[11]]
##
## [[12]]
Ahora incluiremos un id al dataframe y ordenaremos los nombres de las variables
df<-df%>%mutate(id = row_number())
centrales <- setdiff(names(df), c('id', 'centrales'))
df <- df %>%select(
id, one_of(centrales), DEATH_EVENT)Limpiamos el entorno, dejando solo el dataframe
a_borrar <- setdiff(ls(), 'df')
rm(list = c(a_borrar, 'a_borrar'))Guardamos en otro cahe temporal
saveRDS(df, 'cacheT2.rds')Cargamos el cache
df <- readRDS('cacheT2.rds')- Modelización
4.1 - Preparar las funciones que vamos a necesitar
Función para crear una matriz de confusión
confusion <- function(real, scoring, umbral){
conf<- table(real, scoring >=umbral)
if(ncol(conf)==2)return(conf) else return(NULL)
}Funcion para calcular las métricas de los modelos: acierto, precisión, cobertura y F1
metricas<-function(matriz_conf){
acierto <- (matriz_conf[1,1] + matriz_conf[2,2]) / sum(matriz_conf) *100
precision <- matriz_conf[2,2] / (matriz_conf[2,2] + matriz_conf[1,2]) *100
cobertura <- matriz_conf[2,2] / (matriz_conf[2,2] + matriz_conf[2,1]) *100
F1 <- 2*precision*cobertura/(precision+cobertura)
salida<-c(acierto,precision,cobertura,F1)
return(salida)
}Función para probar distintos umbrales y ver el efecto sobre precisión y cobertura?
umbrales<-function(real,scoring){
umbrales<-data.frame(umbral=rep(0,times=19),acierto=rep(0,times=19),precision=rep(0,times=19),cobertura=rep(0,times=19),F1=rep(0,times=19))
cont <- 1
for (cada in seq(0.05,0.95,by = 0.05)){
datos<-metricas(confusion(real,scoring,cada))
registro<-c(cada,datos)
umbrales[cont,]<-registro
cont <- cont + 1
}
return(umbrales)
}Funciones que calculan la curva ROC y el AUC
roc<-function(prediction){
r<-performance(prediction,'tpr','fpr')
plot(r)
}
auc<-function(prediction){
a<-performance(prediction,'auc')
return(a@y.values[[1]])
}Establecemos una semilla para que nos salgan los mismos resultados
set.seed(12345)Generamos una variable aleatoria con una distribución 70-30
df$random <- sample(0:1,size = nrow(df), replace = T, prob = c(0.3, 0.7))Creamos los dos dataframes
train <- filter(df, random==1)
test <- filter(df, random==0)
df$random <- NULL4.2 - Creación del modelo de regresión Nota: Vamos a probar dos algoritmos diferentes para ver cual funciona mejor y aprender como se comparan
independientes <- setdiff(names(df), c('id', 'DEATH_EVENT'))
target <- 'DEATH_EVENT'Creamos la formula para usar en el modelo
formula <- reformulate(independientes, target)4.3.1 - Modelizamos con regresion logistica
formula_rl <- formula
rl <- glm(formula_rl, train, family = binomial(link = 'logit'))
summary(rl)##
## Call:
## glm(formula = formula_rl, family = binomial(link = "logit"),
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.7735 -0.4574 -0.2517 0.3581 2.6980
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) -11.70019 1455.39828 -0.008
## anaemia1 -0.36734 0.53425 -0.688
## diabetes1 0.67734 0.53843 1.258
## high_blood_pressure1 -0.15296 0.51489 -0.297
## sex1 -0.22324 0.58630 -0.381
## smoking1 -0.43554 0.60485 -0.720
## EDAD_DISC02 > 70 2.39680 0.67055 3.574
## creatinine_phosphokinase_DISC02_DE_100_A_200 0.75324 0.76042 0.991
## creatinine_phosphokinase_DISC03_DE_200_A_600 0.08841 0.75555 0.117
## creatinine_phosphokinase_DISC04_DE_600_A_2000 0.48823 0.85438 0.571
## creatinine_phosphokinase_DISC05_MAYOR_2000 -0.80774 1.52275 -0.530
## ejection_fraction_DISC02_DE_30_A_60 -1.63141 0.52458 -3.110
## ejection_fraction_DISC03_MAYOR_DE_60 2.01759 1.74533 1.156
## serum_creatinine_DISC02_DE_1.0_A_2.0 1.04515 0.56577 1.847
## serum_creatinine_DISC03_MAYOR_DE_02 1.30003 0.76119 1.708
## serum_sodium_DISC02_DE_1.0_A_2.0 13.41647 1455.39777 0.009
## serum_sodium_DISC03_MAYOR_DE_02 13.26975 1455.39797 0.009
## time_DISC02 > 73 -3.78640 0.60013 -6.309
## Pr(>|z|)
## (Intercept) 0.993586
## anaemia1 0.491720
## diabetes1 0.208397
## high_blood_pressure1 0.766420
## sex1 0.703382
## smoking1 0.471477
## EDAD_DISC02 > 70 0.000351 ***
## creatinine_phosphokinase_DISC02_DE_100_A_200 0.321903
## creatinine_phosphokinase_DISC03_DE_200_A_600 0.906846
## creatinine_phosphokinase_DISC04_DE_600_A_2000 0.567700
## creatinine_phosphokinase_DISC05_MAYOR_2000 0.595802
## ejection_fraction_DISC02_DE_30_A_60 0.001871 **
## ejection_fraction_DISC03_MAYOR_DE_60 0.247682
## serum_creatinine_DISC02_DE_1.0_A_2.0 0.064704 .
## serum_creatinine_DISC03_MAYOR_DE_02 0.087659 .
## serum_sodium_DISC02_DE_1.0_A_2.0 0.992645
## serum_sodium_DISC03_MAYOR_DE_02 0.992725
## time_DISC02 > 73 0.00000000028 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 240.60 on 188 degrees of freedom
## Residual deviance: 123.18 on 171 degrees of freedom
## AIC: 159.18
##
## Number of Fisher Scoring iterations: 14Seleccionamos las variables más predictoras según rl. Estas son las que tienen tres estrellas, dos estrellas y un punto.
a_mantener <- c('EDAD_DISC', 'ejection_fraction_DISC', 'serum_creatinine_DISC','time_DISC')formula_rl <- reformulate(a_mantener, target)
rl <- glm(formula_rl, train, family = binomial(link = 'logit'))
summary(rl)##
## Call:
## glm(formula = formula_rl, family = binomial(link = "logit"),
## data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4315 -0.4176 -0.2361 0.3573 2.6809
##
## Coefficients:
## Estimate Std. Error z value
## (Intercept) 1.5493 0.5863 2.642
## EDAD_DISC02 > 70 1.9962 0.5732 3.482
## ejection_fraction_DISC02_DE_30_A_60 -1.5006 0.4999 -3.002
## ejection_fraction_DISC03_MAYOR_DE_60 1.4178 1.5444 0.918
## serum_creatinine_DISC02_DE_1.0_A_2.0 1.1701 0.5411 2.163
## serum_creatinine_DISC03_MAYOR_DE_02 1.3534 0.6830 1.981
## time_DISC02 > 73 -3.6144 0.5485 -6.590
## Pr(>|z|)
## (Intercept) 0.008231 **
## EDAD_DISC02 > 70 0.000497 ***
## ejection_fraction_DISC02_DE_30_A_60 0.002686 **
## ejection_fraction_DISC03_MAYOR_DE_60 0.358595
## serum_creatinine_DISC02_DE_1.0_A_2.0 0.030578 *
## serum_creatinine_DISC03_MAYOR_DE_02 0.047537 *
## time_DISC02 > 73 0.0000000000441 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 240.60 on 188 degrees of freedom
## Residual deviance: 127.95 on 182 degrees of freedom
## AIC: 141.95
##
## Number of Fisher Scoring iterations: 5Y calculamos el pseudo R cuadrado:
pr2_rl <- 1 - (rl$deviance/rl$null.deviance)
pr2_rl## [1] 0.4682044Aplicamos el modelo al conjunto de test, generando un vector con las probabilidades
rl_predict <- predict(rl, test, type = 'response')Vemos que pinta tiene
plot(rl_predict~test$DEATH_EVENT)
Ahora tenemos que transformar la probabilidad de muete del paciente Con la funcion umbrales probamos diferentes cortes
umb_rl <- umbrales(test$DEATH_EVENT, rl_predict)
umb_rl## umbral acierto precision cobertura F1
## 1 0.05 55.45455 40.00000 96.96970 56.63717
## 2 0.10 74.54545 54.38596 93.93939 68.88889
## 3 0.15 80.90909 62.00000 93.93939 74.69880
## 4 0.20 81.81818 63.26531 93.93939 75.60976
## 5 0.25 81.81818 63.26531 93.93939 75.60976
## 6 0.30 82.72727 68.42105 78.78788 73.23944
## 7 0.35 80.90909 66.66667 72.72727 69.56522
## 8 0.40 80.90909 66.66667 72.72727 69.56522
## 9 0.45 85.45455 79.31034 69.69697 74.19355
## 10 0.50 85.45455 79.31034 69.69697 74.19355
## 11 0.55 84.54545 80.76923 63.63636 71.18644
## 12 0.60 84.54545 80.76923 63.63636 71.18644
## 13 0.65 84.54545 80.76923 63.63636 71.18644
## 14 0.70 84.54545 80.76923 63.63636 71.18644
## 15 0.75 84.54545 80.76923 63.63636 71.18644
## 16 0.80 81.81818 84.21053 48.48485 61.53846
## 17 0.85 80.00000 86.66667 39.39394 54.16667
## 18 0.90 79.09091 91.66667 33.33333 48.88889
## 19 0.95 72.72727 80.00000 12.12121 21.05263Seleccionamos el umbral que maximiza la F1
umbral_final_rl<-umb_rl[which.max(umb_rl$F1),1]
umbral_final_rl## [1] 0.2Evaluamos la matriz de confusión y las métricas con el umbral optimizado
confusion(test$DEATH_EVENT,rl_predict,umbral_final_rl)##
## real FALSE TRUE
## 0 59 18
## 1 2 31rl_metricas<-filter(umb_rl,umbral==umbral_final_rl)
rl_metricas## umbral acierto precision cobertura F1
## 1 0.2 81.81818 63.26531 93.93939 75.60976Evaluamos la ROC
#creamos el objeto prediction
rl_prediction<-prediction(rl_predict,test$DEATH_EVENT)
#visualizamos la ROC
roc(rl_prediction)
Sacamos las métricas definitivas incluyendo el AUC
rl_metricas<-cbind(rl_metricas,AUC=round(auc(rl_prediction),2)*100)
print(t(rl_metricas))## [,1]
## umbral 0.20000
## acierto 81.81818
## precision 63.26531
## cobertura 93.93939
## F1 75.60976
## AUC 90.000004.3 Modelizamos con Random Forest
Creamos el modelo
formula_rf <- formula
rf <- randomForest(formula_rf, train, importance = T)
rf##
## Call:
## randomForest(formula = formula_rf, data = train, importance = T)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 3
##
## OOB estimate of error rate: 19.05%
## Confusion matrix:
## 0 1 class.error
## 0 113 13 0.1031746
## 1 23 40 0.3650794importance(rf)[,3:4]## MeanDecreaseAccuracy MeanDecreaseGini
## anaemia -1.7925556 2.421144
## diabetes -1.8244133 2.855078
## high_blood_pressure 2.6158302 3.052748
## sex -1.4122765 2.449528
## smoking -0.7211306 2.289087
## EDAD_DISC 14.4912948 6.610334
## creatinine_phosphokinase_DISC -0.9126329 7.947236
## ejection_fraction_DISC 15.7227233 7.942199
## serum_creatinine_DISC 10.3887072 7.918194
## serum_sodium_DISC 0.6863976 1.814491
## time_DISC 49.4740931 27.026928Visualizamos las variables mas importantes
varImpPlot(rf)
Como hay dos criterios vamos a crear una única variable agregada y visualizarla para tener una mejor idea de la importancia de cada variable
importancia <- importance(rf)[,3:4]
importancia_norm <- as.data.frame(scale(importancia))
df <- as.data.frame(df)
#normalizamos para poner las dos variables en la misma escala. los valores negativos son las que menos predicen y los positivos las que mas.
importancia_norm <- importancia_norm%>% mutate(Variable = rownames(importancia_norm),
Imp_tot= MeanDecreaseAccuracy + MeanDecreaseGini)%>%
mutate(Imp_tot= Imp_tot + abs(min(Imp_tot))) %>%
arrange(desc(Imp_tot))%>%
select(Variable, Imp_tot, MeanDecreaseAccuracy, MeanDecreaseGini)Hacemos un gráfico para ver la curva de caída de importancia
ggplot(importancia_norm, aes(reorder(Variable, -Imp_tot), Imp_tot))+geom_bar(stat='identity')+theme(axis.text.x=element_text(angle=90, size=7))
importancia_norm## Variable Imp_tot MeanDecreaseAccuracy
## 1 time_DISC 6.74148986 2.7151609
## 2 ejection_fraction_DISC 1.90514077 0.5117803
## 3 EDAD_DISC 1.64100303 0.4313893
## 4 serum_creatinine_DISC 1.55361007 0.1635613
## 5 creatinine_phosphokinase_DISC 0.81983476 -0.5742207
## 6 high_blood_pressure 0.37492874 -0.3438731
## 7 serum_sodium_DISC 0.07813769 -0.4698317
## 8 diabetes 0.05778658 -0.6337442
## 9 smoking 0.05172666 -0.5617189
## 10 sex 0.02874150 -0.6068388
## 11 anaemia 0.00000000 -0.6316644
## MeanDecreaseGini
## 1 2.821565337
## 2 0.188596867
## 3 0.004850106
## 4 0.185285156
## 5 0.189291867
## 6 -0.485961760
## 7 -0.656794242
## 8 -0.513232829
## 9 -0.591318020
## 10 -0.569183315
## 11 -0.573099168Vemos una clara diferencia de la variable máspredictora, time_DISC, con el respecto al resto de variables.Podemos coger por ejemplo hasta high_blood_presure incluido, que tiene una importancia total de lrededor de 0.2. También se podrían hacer diferentes pruebas, pero en este caso lo dejaremos asi
a_mantener <- importancia_norm%>%
filter(Imp_tot > 0.2)%>%
select(Variable)
#Extraemos los nombres como un vector
a_mantener <- as.character((a_mantener$Variable))Creamos de nuevo el modelo con las nuevas variables
formula_rf <- reformulate(a_mantener, target)
rf<-randomForest(formula_rf, train, importance = T)
rf##
## Call:
## randomForest(formula = formula_rf, data = train, importance = T)
## Type of random forest: classification
## Number of trees: 500
## No. of variables tried at each split: 2
##
## OOB estimate of error rate: 16.93%
## Confusion matrix:
## 0 1 class.error
## 0 112 14 0.1111111
## 1 18 45 0.2857143Aplicamos el modelo al conjunto de test, generando un vector con las probabilidades
rf_predict<-predict(rf,test,type = 'prob')[,2]plot(rf_predict~test$DEATH_EVENT)
umb_rf<-umbrales(test$DEATH_EVENT,rf_predict)
umb_rf## umbral acierto precision cobertura F1
## 1 0.05 75.45455 55.17241 96.969697 70.329670
## 2 0.10 75.45455 55.17241 96.969697 70.329670
## 3 0.15 76.36364 56.14035 96.969697 71.111111
## 4 0.20 76.36364 56.60377 90.909091 69.767442
## 5 0.25 77.27273 58.00000 87.878788 69.879518
## 6 0.30 80.00000 62.22222 84.848485 71.794872
## 7 0.35 83.63636 69.23077 81.818182 75.000000
## 8 0.40 81.81818 68.57143 72.727273 70.588235
## 9 0.45 83.63636 72.72727 72.727273 72.727273
## 10 0.50 83.63636 72.72727 72.727273 72.727273
## 11 0.55 81.81818 70.96774 66.666667 68.750000
## 12 0.60 80.00000 68.96552 60.606061 64.516129
## 13 0.65 81.81818 78.26087 54.545455 64.285714
## 14 0.70 81.81818 80.95238 51.515152 62.962963
## 15 0.75 82.72727 85.00000 51.515152 64.150943
## 16 0.80 77.27273 78.57143 33.333333 46.808511
## 17 0.85 74.54545 72.72727 24.242424 36.363636
## 18 0.90 70.90909 66.66667 6.060606 11.111111
## 19 0.95 70.90909 100.00000 3.030303 5.882353umbral_final_rf<-umb_rf[which.max(umb_rf$F1),1]
umbral_final_rf## [1] 0.35confusion(test$DEATH_EVENT,rf_predict,umbral_final_rf)##
## real FALSE TRUE
## 0 65 12
## 1 6 27rf_metricas<-filter(umb_rf,umbral==umbral_final_rf)
rf_metricas## umbral acierto precision cobertura F1
## 1 0.35 83.63636 69.23077 81.81818 75rf_prediction <- prediction (rf_predict, test$DEATH_EVENT)
roc(rf_prediction)
rf_metricas<- cbind(rf_metricas, AUC=round(auc(rf_prediction),2)*100)
print(t(rf_metricas))## [,1]
## umbral 0.35000
## acierto 83.63636
## precision 69.23077
## cobertura 81.81818
## F1 75.00000
## AUC 89.00000Comparamos los dos modelos
comparativa <- rbind(rl_metricas, rf_metricas)rownames(comparativa)<- c('Regresión Logistica', 'Random Forest')
t(comparativa)## Regresión Logistica Random Forest
## umbral 0.20000 0.35000
## acierto 81.81818 83.63636
## precision 63.26531 69.23077
## cobertura 93.93939 81.81818
## F1 75.60976 75.00000
## AUC 90.00000 89.00000Parece que el modelo de Regresión Logistica presenta un AUC mayor asi que nos quedamos con este modelo. A partir de ahora escribimos el scoring definitivo final en el dataset y guardamos el modelo.
df$SCORING_EVENTO <- predict(rl, df, type='response')Hacemos un head y un select de dos variables del df, dandome únicamente los 20 primeros clientes.
head(select(df, id, SCORING_EVENTO),100)## id SCORING_EVENTO
## 1 1 0.99112495
## 2 2 0.77184863
## 3 3 0.93816186
## 4 4 0.93816186
## 5 5 0.94797761
## 6 6 0.96764910
## 7 7 0.99112495
## 8 8 0.77184863
## 9 9 0.98428376
## 10 10 0.96764910
## 11 11 0.96764910
## 12 12 0.82480700
## 13 13 0.93816186
## 14 14 0.77184863
## 15 15 0.82480700
## 16 16 0.96139387
## 17 17 0.88542440
## 18 18 0.82480700
## 19 19 0.82480700
## 20 20 0.77184863
## 21 21 0.93816186
## 22 22 0.93816186
## 23 23 0.51215716
## 24 24 0.51215716
## 25 25 0.99112495
## 26 26 0.96139387
## 27 27 0.88542440
## 28 28 0.77184863
## 29 29 0.80250641
## 30 30 0.99112495
## 31 31 0.96139387
## 32 32 0.96764910
## 33 33 0.51215716
## 34 34 0.93816186
## 35 35 0.51215716
## 36 36 0.80250641
## 37 37 0.88542440
## 38 38 0.88542440
## 39 39 0.94797761
## 40 40 0.80250641
## 41 41 0.93816186
## 42 42 0.93816186
## 43 43 0.77184863
## 44 44 0.88542440
## 45 45 0.77184863
## 46 46 0.77184863
## 47 47 0.82480700
## 48 48 0.51215716
## 49 49 0.99260002
## 50 50 0.82480700
## 51 51 0.82480700
## 52 52 0.93816186
## 53 53 0.98688082
## 54 54 0.51215716
## 55 55 0.80250641
## 56 56 0.99112495
## 57 57 0.80250641
## 58 58 0.51215716
## 59 59 0.93816186
## 60 60 0.99112495
## 61 61 0.82480700
## 62 62 0.80250641
## 63 63 0.77184863
## 64 64 0.51215716
## 65 65 0.98428376
## 66 66 0.94797761
## 67 67 0.93816186
## 68 68 0.97195380
## 69 69 0.93816186
## 70 70 0.93816186
## 71 71 0.51215716
## 72 72 0.51215716
## 73 73 0.88542440
## 74 74 0.77184863
## 75 75 0.93816186
## 76 76 0.82480700
## 77 77 0.02749749
## 78 78 0.08350599
## 79 79 0.17227569
## 80 80 0.02749749
## 81 81 0.08350599
## 82 82 0.08350599
## 83 83 0.32921015
## 84 84 0.40144613
## 85 85 0.11252998
## 86 86 0.02749749
## 87 87 0.08350599
## 88 88 0.02749749
## 89 89 0.02749749
## 90 90 0.29007615
## 91 91 0.02749749
## 92 92 0.02749749
## 93 93 0.08350599
## 94 94 0.29007615
## 95 95 0.02749749
## 96 96 0.02749749
## 97 97 0.29007615
## 98 98 0.08350599
## 99 99 0.29007615
## 100 100 0.08350599Accedemos a la probabilidad de que aparezca evento del paciente con un id=44
df[df$id=='44', 'SCORING_EVENTO']## [1] 0.8854244Filtramos DEATH_EVENT = 0
df <- df%>% filter(DEATH_EVENT==0)Para este caso será necesario realizar una campaña de prevención de muerte, para aquello que tiene mayor probabilidad
tamaño_campaña = 100 ## Tamaño campaña = a registros en los que DEAT_EVENT=0
dimension_prevencion <- df%>%
arrange(desc(SCORING_EVENTO))%>%
slice(1:tamaño_campaña)%>%
select(id,SCORING_EVENTO)Previsualizamos la salida
head(dimension_prevencion, 50)## id SCORING_EVENTO
## 1 65 0.9842838
## 2 39 0.9479776
## 3 21 0.9381619
## 4 34 0.9381619
## 5 44 0.8854244
## 6 15 0.8248070
## 7 57 0.8025064
## 8 63 0.7718486
## 9 74 0.7718486
## 10 103 0.7504811
## 11 194 0.7504811
## 12 230 0.7504811
## 13 24 0.5121572
## 14 58 0.5121572
## 15 71 0.5121572
## 16 72 0.5121572
## 17 118 0.4461624
## 18 191 0.4461624
## 19 84 0.4014461
## 20 102 0.4014461
## 21 135 0.4014461
## 22 136 0.4014461
## 23 159 0.4014461
## 24 213 0.4014461
## 25 216 0.4014461
## 26 226 0.4014461
## 27 236 0.4014461
## 28 237 0.4014461
## 29 212 0.3435897
## 30 138 0.3292102
## 31 204 0.3292102
## 32 229 0.3292102
## 33 248 0.3292102
## 34 283 0.3292102
## 35 90 0.2900762
## 36 97 0.2900762
## 37 99 0.2900762
## 38 101 0.2900762
## 39 113 0.2900762
## 40 156 0.2900762
## 41 227 0.2900762
## 42 238 0.2900762
## 43 242 0.2900762
## 44 271 0.2900762
## 45 79 0.1722757
## 46 205 0.1722757
## 47 208 0.1722757
## 48 244 0.1722757
## 49 290 0.1722757
## 50 104 0.1125300Vamos a ver gráficamente si de esta forma estamos aprovechando el potencial de nuestro modelo
penetracion_target <- mean(as.numeric(as.character(df$DEATH_EVENT)))
df %>%
arrange(desc(SCORING_EVENTO)) %>%
ggplot(aes(y = SCORING_EVENTO, x=seq_along(SCORING_EVENTO))) + ##seq_along para establecer una secuencia de números naturales de SCORING_EVENTO en el eje x,
geom_line() +
geom_vline(xintercept = tamaño_campaña, col = 'orange') +
geom_hline(yintercept = penetracion_target, col='blue') +
labs(x='PACIENTES ORDENADOS POR SCORING')
Se puede ver en la gráfica que con la campaña que hemos establecido, estamos cogiendo todos los pacientes con una probabiidad alta.