Introdución a la regresión logística con R

Navegar al directorio donde tenemos el dataset guardado

getwd()
## [1] "/home/diego/Proyecto Curso Big Data/R/Regresion-logistica"
setwd("/home/diego/Proyecto Curso Big Data/R/Regresion-logistica/")

Read the dataset framingham

ds=read.csv("framingham.csv")

Get the structure of the dataset

str(ds)
## 'data.frame':    4240 obs. of  16 variables:
##  $ male           : int  1 0 1 0 0 0 0 0 1 1 ...
##  $ age            : int  39 46 48 61 46 43 63 45 52 43 ...
##  $ education      : int  4 2 1 3 3 2 1 2 1 1 ...
##  $ currentSmoker  : int  0 0 1 1 1 0 0 1 0 1 ...
##  $ cigsPerDay     : int  0 0 20 30 23 0 0 20 0 30 ...
##  $ BPMeds         : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ prevalentStroke: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ prevalentHyp   : int  0 0 0 1 0 1 0 0 1 1 ...
##  $ diabetes       : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ totChol        : int  195 250 245 225 285 228 205 313 260 225 ...
##  $ sysBP          : num  106 121 128 150 130 ...
##  $ diaBP          : num  70 81 80 95 84 110 71 71 89 107 ...
##  $ BMI            : num  27 28.7 25.3 28.6 23.1 ...
##  $ heartRate      : int  80 95 75 65 85 77 60 79 76 93 ...
##  $ glucose        : int  77 76 70 103 85 99 85 78 79 88 ...
##  $ TenYearCHD     : int  0 0 0 1 0 0 1 0 0 0 ...

Dividir el dataset en dos: datos de entrenamiento y datos de test. Necesitamos la librería caTools, si no la tienes descargada, descargala:

#install.packages("caTools")

Cargar libreria

library(caTools)

Establecer la semilla

set.seed(1000)

Probabilidades con un ratio del 65% para train y 35% para test

split = sample.split(ds$TenYearCHD, SplitRatio = 0.65)

Vamos a dividir el dataset en dos subconjuntos. Uno para el entrenamiento del modelo al que llamaremos train, y otro para pruebarlo que llamaremos test.

train = subset(ds, split == TRUE)
test = subset(ds, split == FALSE)

La función de regresión logística es: glm() El primer parámetro indica cual es la variable dependiente. El resto, tras el símbolo ~ son las variables independientes. En este caso el ‘.’ indica que se utilicen todas las variables restantes del dataset

framinghamLog = glm(TenYearCHD ~ ., data = train, family = binomial)
summary(framinghamLog)
## 
## Call:
## glm(formula = TenYearCHD ~ ., family = binomial, data = train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8487  -0.6007  -0.4257  -0.2842   2.8369  
## 
## Coefficients:
##                  Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     -7.886574   0.890729  -8.854  < 2e-16 ***
## male             0.528457   0.135443   3.902 9.55e-05 ***
## age              0.062055   0.008343   7.438 1.02e-13 ***
## education       -0.058923   0.062430  -0.944  0.34525    
## currentSmoker    0.093240   0.194008   0.481  0.63080    
## cigsPerDay       0.015008   0.007826   1.918  0.05514 .  
## BPMeds           0.311221   0.287408   1.083  0.27887    
## prevalentStroke  1.165794   0.571215   2.041  0.04126 *  
## prevalentHyp     0.315818   0.171765   1.839  0.06596 .  
## diabetes        -0.421494   0.407990  -1.033  0.30156    
## totChol          0.003835   0.001377   2.786  0.00533 ** 
## sysBP            0.011344   0.004566   2.485  0.01297 *  
## diaBP           -0.004740   0.008001  -0.592  0.55353    
## BMI              0.010723   0.016157   0.664  0.50689    
## heartRate       -0.008099   0.005313  -1.524  0.12739    
## glucose          0.008935   0.002836   3.150  0.00163 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 2020.7  on 2384  degrees of freedom
## Residual deviance: 1792.3  on 2369  degrees of freedom
##   (371 observations deleted due to missingness)
## AIC: 1824.3
## 
## Number of Fisher Scoring iterations: 5

Vamos a crear el primer conjunto de predicciones con el modelo creado anteriormente

predictTest = predict(framinghamLog, type = "response", newdata = test)

Vamos a ver los resultados de la predicción.

table(test$TenYearCHD, predictTest > 0.2)
##    
##     FALSE TRUE
##   0   865  210
##   1    88  110

Precisión: casos positivos / casos totales. Calcular para ejemplo anterior teniendo en cuenta el resultado de la tabla.

Para realizar las predicciones primero debemos cargar el paquete ROCR. Y luego, usaremos la función de predicción del paquete ROCR para hacer nuestras predicciones.

#install.packages("ROCR")
library(ROCR)
## Loading required package: gplots
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess

Predicciones vs reales

predictionsROCR = prediction(predictTest, test$TenYearCHD)

ROC Curve

ROCRPerf = performance(predictionsROCR, "tpr", "fpr")
plot(ROCRPerf, colorize = TRUE, print.cutoffs.at=seq(0,1,0.1),text.adj=c(-0.2,1.7))

Calcular AUC fuera de la muestra.

auc = as.numeric(performance(predictionsROCR,"auc")@y.values)
auc
## [1] 0.7421095

Esto nos dará el valor AUC en nuestro conjunto de pruebas. Así que tenemos un AUC de alrededor del 72% en nuestro conjunto de pruebas,lo que significa que el modelo puede diferenciar entre riesgo bajo pacientes y pacientes de alto riesgo bastante bien.