Teoría

La regresión logística es un modelo estadístico de clasificación binaria, que estima la probabilidad de que ocurra un evento (valor 1) frente a que no ocurra (valor 0), en función de variables independientes. # Instalar paquetes y llamar librerías

#install.packages("titanic")
library(titanic)
#install.packages("caret")
library(caret)
## Loading required package: ggplot2
## Loading required package: lattice
#install.packages("tidyverse")
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.0     ✔ readr     2.1.6
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ lubridate 1.9.5     ✔ tibble    3.3.1
## ✔ purrr     1.2.1     ✔ tidyr     1.3.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ purrr::lift()   masks caret::lift()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Crear la base de datos

df <- read.csv("/Users/giuliana/Downloads/heart.csv")

Entender la base de datos

summary(df)
##       age             sex               cp            trestbps    
##  Min.   :29.00   Min.   :0.0000   Min.   :0.0000   Min.   : 94.0  
##  1st Qu.:48.00   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:120.0  
##  Median :56.00   Median :1.0000   Median :1.0000   Median :130.0  
##  Mean   :54.43   Mean   :0.6956   Mean   :0.9424   Mean   :131.6  
##  3rd Qu.:61.00   3rd Qu.:1.0000   3rd Qu.:2.0000   3rd Qu.:140.0  
##  Max.   :77.00   Max.   :1.0000   Max.   :3.0000   Max.   :200.0  
##       chol          fbs            restecg          thalach     
##  Min.   :126   Min.   :0.0000   Min.   :0.0000   Min.   : 71.0  
##  1st Qu.:211   1st Qu.:0.0000   1st Qu.:0.0000   1st Qu.:132.0  
##  Median :240   Median :0.0000   Median :1.0000   Median :152.0  
##  Mean   :246   Mean   :0.1493   Mean   :0.5298   Mean   :149.1  
##  3rd Qu.:275   3rd Qu.:0.0000   3rd Qu.:1.0000   3rd Qu.:166.0  
##  Max.   :564   Max.   :1.0000   Max.   :2.0000   Max.   :202.0  
##      exang           oldpeak          slope             ca        
##  Min.   :0.0000   Min.   :0.000   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:0.0000   1st Qu.:0.000   1st Qu.:1.000   1st Qu.:0.0000  
##  Median :0.0000   Median :0.800   Median :1.000   Median :0.0000  
##  Mean   :0.3366   Mean   :1.072   Mean   :1.385   Mean   :0.7541  
##  3rd Qu.:1.0000   3rd Qu.:1.800   3rd Qu.:2.000   3rd Qu.:1.0000  
##  Max.   :1.0000   Max.   :6.200   Max.   :2.000   Max.   :4.0000  
##       thal           target      
##  Min.   :0.000   Min.   :0.0000  
##  1st Qu.:2.000   1st Qu.:0.0000  
##  Median :2.000   Median :1.0000  
##  Mean   :2.324   Mean   :0.5132  
##  3rd Qu.:3.000   3rd Qu.:1.0000  
##  Max.   :3.000   Max.   :1.0000
str(df)
## 'data.frame':    1025 obs. of  14 variables:
##  $ age     : int  52 53 70 61 62 58 58 55 46 54 ...
##  $ sex     : int  1 1 1 1 0 0 1 1 1 1 ...
##  $ cp      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ trestbps: int  125 140 145 148 138 100 114 160 120 122 ...
##  $ chol    : int  212 203 174 203 294 248 318 289 249 286 ...
##  $ fbs     : int  0 1 0 0 1 0 0 0 0 0 ...
##  $ restecg : int  1 0 1 1 1 0 2 0 0 0 ...
##  $ thalach : int  168 155 125 161 106 122 140 145 144 116 ...
##  $ exang   : int  0 1 1 0 0 0 0 1 0 1 ...
##  $ oldpeak : num  1 3.1 2.6 0 1.9 1 4.4 0.8 0.8 3.2 ...
##  $ slope   : int  2 0 0 2 1 1 0 1 2 1 ...
##  $ ca      : int  2 0 0 1 3 0 3 1 0 2 ...
##  $ thal    : int  3 3 3 3 2 2 1 3 3 2 ...
##  $ target  : int  0 0 0 0 0 1 0 0 0 0 ...
df <- df[, c("target", "age","sex","cp",
"oldpeak","slope","ca","thal")]
df <- na.omit(df)
df$target <- as.factor(df$target)
df$sex <- as.factor(df$sex)
df$cp <- as.factor(df$cp)

Crear el modelo

modelo <- glm(target ~ ., data=df, family=binomial)
summary(modelo)
## 
## Call:
## glm(formula = target ~ ., family = binomial, data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.29415    0.79117   5.428 5.71e-08 ***
## age         -0.03969    0.01069  -3.712 0.000205 ***
## sex1        -1.55890    0.22591  -6.901 5.18e-12 ***
## cp1          1.69405    0.27080   6.256 3.96e-10 ***
## cp2          2.40989    0.22975  10.489  < 2e-16 ***
## cp3          2.41829    0.32258   7.497 6.54e-14 ***
## oldpeak     -0.61905    0.11051  -5.602 2.12e-08 ***
## slope        0.74670    0.17606   4.241 2.22e-05 ***
## ca          -0.73568    0.09967  -7.381 1.57e-13 ***
## thal        -0.86820    0.14287  -6.077 1.23e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1420.24  on 1024  degrees of freedom
## Residual deviance:  772.88  on 1015  degrees of freedom
## AIC: 792.88
## 
## Number of Fisher Scoring iterations: 6

Probar el modelo

prueba <- df[c(5, 20), ]

probabilidad <- predict(modelo, newdata = prueba, type = "response")
cbind(prueba, Probabilidad_Enfermedad = probabilidad)
##    target age sex cp oldpeak slope ca thal Probabilidad_Enfermedad
## 5       0  62   0  0     1.9     1  3    2              0.07313319
## 20      1  58   1  2     0.0     2  0    2              0.93087047

En los dos casos que se hacen en la prueba (filas 5 y 20), el modelo predice una probabilidad de aproximadamente 98.4% de que sean setosa. Esto significa que el modelo está muy seguro de que estas observaciones pertenecen a esa especie. Además, coincide con la variable real Species, que también es setosa en ambos casos, por lo que el modelo está clasificando correctamente estos ejemplos. En general, para estas observaciones el modelo muestra un desempeño muy bueno y alta confianza en su predicción.

LS0tCnRpdGxlOiAiSGVhcnQgZGlzZWFzZSBSZWdyZXNpw7NuIExvZ8Otc3RpY2EiCmF1dGhvcjogIkdpdWxpYW5hIE1hbmNlcmEgRmxvcmVzIgpkYXRlOiAiYHIgU3lzLkRhdGUoKWAiCm91dHB1dDogCiAgaHRtbF9kb2N1bWVudDoKICAgIHRvYzogVFJVRQogICAgdG9jX2Zsb2F0OiBUUlVFCiAgICBjb2RlX2Rvd25sb2FkOiBUUlVFCiAgICB0aGVtZTogY29zbW8KLS0tCjxjZW50ZXI+CiFbXShodHRwczovL2ltYWdlcy5jdGZhc3NldHMubmV0L3V0N3J6djh5ZWhwZi8xRGhDM3VYM0VlS25qVTAyTFd5VFhILzljODJlNmFlODI2NjJlZDU5MDNlYWZiNDBkODg4ZDkwLzhfTWFpbl9UeXBlc19vZl9IZWFydF9EaXNlYXNlLmpwZz93PTE4MDAmaD05MDAmZmw9cHJvZ3Jlc3NpdmUmcT01MCZmbT1qcGcpCjwvY2VudGVyPgoKIyA8c3BhbiBzdHlsZT0iY29sb3I6Ymx1ZSI+IFRlb3LDrWEgPC9zcGFuPgpMYSAqKnJlZ3Jlc2nDs24gbG9nw61zdGljYSoqIGVzIHVuIG1vZGVsbyBlc3RhZMOtc3RpY28gZGUgY2xhc2lmaWNhY2nDs24KYmluYXJpYSwgcXVlIGVzdGltYSBsYSBwcm9iYWJpbGlkYWQgZGUgcXVlIG9jdXJyYSB1biBldmVudG8gKHZhbG9yIDEpCmZyZW50ZSBhIHF1ZSBubyBvY3VycmEgKHZhbG9yIDApLCBlbiBmdW5jacOzbiBkZSB2YXJpYWJsZXMgaW5kZXBlbmRpZW50ZXMuCiMgPHNwYW4gc3R5bGU9ImNvbG9yOmJsdWUiPiBJbnN0YWxhciBwYXF1ZXRlcyB5IGxsYW1hciBsaWJyZXLDrWFzIDwvc3Bhbj4KYGBge3J9CiNpbnN0YWxsLnBhY2thZ2VzKCJ0aXRhbmljIikKbGlicmFyeSh0aXRhbmljKQojaW5zdGFsbC5wYWNrYWdlcygiY2FyZXQiKQpsaWJyYXJ5KGNhcmV0KQojaW5zdGFsbC5wYWNrYWdlcygidGlkeXZlcnNlIikKbGlicmFyeSh0aWR5dmVyc2UpCmBgYAojIDxzcGFuIHN0eWxlPSJjb2xvcjpibHVlIj4gQ3JlYXIgbGEgYmFzZSBkZSBkYXRvcyA8L3NwYW4+CmBgYHtyfQpkZiA8LSByZWFkLmNzdigiL1VzZXJzL2dpdWxpYW5hL0Rvd25sb2Fkcy9oZWFydC5jc3YiKQpgYGAKIyA8c3BhbiBzdHlsZT0iY29sb3I6Ymx1ZSI+IEVudGVuZGVyIGxhIGJhc2UgZGUgZGF0b3MgPC9zcGFuPgpgYGB7cn0Kc3VtbWFyeShkZikKc3RyKGRmKQpkZiA8LSBkZlssIGMoInRhcmdldCIsICJhZ2UiLCJzZXgiLCJjcCIsCiJvbGRwZWFrIiwic2xvcGUiLCJjYSIsInRoYWwiKV0KZGYgPC0gbmEub21pdChkZikKZGYkdGFyZ2V0IDwtIGFzLmZhY3RvcihkZiR0YXJnZXQpCmRmJHNleCA8LSBhcy5mYWN0b3IoZGYkc2V4KQpkZiRjcCA8LSBhcy5mYWN0b3IoZGYkY3ApCmBgYAojIDxzcGFuIHN0eWxlPSJjb2xvcjpibHVlIj4gQ3JlYXIgZWwgbW9kZWxvIDwvc3Bhbj4KYGBge3J9Cm1vZGVsbyA8LSBnbG0odGFyZ2V0IH4gLiwgZGF0YT1kZiwgZmFtaWx5PWJpbm9taWFsKQpzdW1tYXJ5KG1vZGVsbykKYGBgCiMgPHNwYW4gc3R5bGU9ImNvbG9yOmJsdWUiPiBQcm9iYXIgZWwgbW9kZWxvIDwvc3Bhbj4KYGBge3J9CnBydWViYSA8LSBkZltjKDUsIDIwKSwgXQoKcHJvYmFiaWxpZGFkIDwtIHByZWRpY3QobW9kZWxvLCBuZXdkYXRhID0gcHJ1ZWJhLCB0eXBlID0gInJlc3BvbnNlIikKY2JpbmQocHJ1ZWJhLCBQcm9iYWJpbGlkYWRfRW5mZXJtZWRhZCA9IHByb2JhYmlsaWRhZCkKYGBgCgpFbiBsb3MgZG9zIGNhc29zIHF1ZSBzZSBoYWNlbiBlbiBsYSBwcnVlYmEgKGZpbGFzIDUgeSAyMCksIGVsIG1vZGVsbyBwcmVkaWNlIHVuYSBwcm9iYWJpbGlkYWQgZGUgYXByb3hpbWFkYW1lbnRlIDk4LjQlIGRlIHF1ZSBzZWFuIHNldG9zYS4gRXN0byBzaWduaWZpY2EgcXVlIGVsIG1vZGVsbyBlc3TDoSBtdXkgc2VndXJvIGRlIHF1ZSBlc3RhcyBvYnNlcnZhY2lvbmVzIHBlcnRlbmVjZW4gYSBlc2EgZXNwZWNpZS4gQWRlbcOhcywgY29pbmNpZGUgY29uIGxhIHZhcmlhYmxlIHJlYWwgU3BlY2llcywgcXVlIHRhbWJpw6luIGVzIHNldG9zYSBlbiBhbWJvcyBjYXNvcywgcG9yIGxvIHF1ZSBlbCBtb2RlbG8gZXN0w6EgY2xhc2lmaWNhbmRvIGNvcnJlY3RhbWVudGUgZXN0b3MgZWplbXBsb3MuIEVuIGdlbmVyYWwsIHBhcmEgZXN0YXMgb2JzZXJ2YWNpb25lcyBlbCBtb2RlbG8gbXVlc3RyYSB1biBkZXNlbXBlw7FvIG11eSBidWVubyB5IGFsdGEgY29uZmlhbnphIGVuIHN1IHByZWRpY2Npw7NuLgoK