Este documento incluye los clasificadores KNN, Naive Bayes, regresión logística y una breve implementación de arboles de decisión.

KNN - K Nearest Neighbors

Algoritmo de aprendizaje supervisado basado en distancias (euclidianas normalmente) para determinar, basado en un número “k” de vecinos, a quién se parece más la observación que tratamos de predecir.

Los pasos son los siguientes:

Este algoritmo debe contar con la etiqueta de la observación para definir a qué clase pertenece.

Ventajas

  • Muy simple
  • El entrenamiento también es facil
  • Trabaja con cualquier número de clases (variable independiente)
  • Fácil de agregar más datos

Desventajas

  • Alto costo de predicción (Mucho peor para grandes conjuntos de datos)
  • No es bueno cuando los datasets tienen demasidados features
  • No trabaja bien con datos categoricos
library("class")
library("dplyr")
library("naivebayes")
library("pROC")
library("rpart")
library("rpart.plot")
signs <- read.csv("knn_traffic_signs.csv")

train <- signs %>% filter(sample == "train") %>% select(-c(id, sample))
test  <- signs %>% filter(sample == "test") %>% select(-c(id, sample))
signs_pred <- knn(train = train[-1] , test[-1], cl = train$sign_type)

# Matriz de confusión
signs_actual <- test$sign_type
table(signs_pred, signs_actual)
##             signs_actual
## signs_pred   pedestrian speed stop
##   pedestrian         19     2    0
##   speed               0    17    0
##   stop                0     2   19
# Presición de los resultados
mean(signs_pred == signs_actual)
## [1] 0.9322034
  • El valor de k puede ser inicialmente elegido como la raiz cuadrada del número de observaciones en los datos de entrenamiento, pero este puede ser elegido con mayor exactitud para aumental el poder predictivo con herramientas de cross-validation o gridSearch.

  • Valores bajos de k general alta varianza y valores altos de k alto sesgo.

Podemos también extraer de las proabilidades

# Podemos usar el parametro prob pata extraer las probabilidades de la clase ganadora
sign_pred <- knn(train = train[-1],
                 test = test[-1], cl = train$sign_type, k =7, prob = TRUE)

# Este dato se extrae como atributo del modelo
sign_prob <- attr(sign_pred, "prob")

# Predicciones
head(sign_pred)
## [1] pedestrian pedestrian pedestrian stop       pedestrian pedestrian
## Levels: pedestrian speed stop
# Probabilidades (TAmbien nos servirían para ajustar el valor de K)
head(sign_prob)
## [1] 0.5714286 0.5714286 0.8571429 0.5714286 0.8571429 0.5714286

Al ser un método basado en distancias es necesario manejar las variables en la misma escala!!!

normalize <- function(x){
  return((x - min(x))/(max(x) - min(x)))
}

Naive Bayes (Bayesiano ingenuo)

Aplica métodos bayesianos (Teorema de bayes) para estimar la probabilidad condicional de un resultado.

en lugar de tratar el problema como la intersección de multiples eventos, el algoritmo hace un supuesto “ingenuo” sobre los datos. Asume que los eventos on independientes.

En otras palabras, el efecto de una característica particular es independiente de otras características.

Cuando los eventos son independientes, la probabilidad conjunta puede ser calculada multiplicando las probabilidades individuales. Bajo este supuesto, el algoritmo no necesita observar todas las posibles intersecciones, siendo computacionalmente más viable.

\[ P(h | D) = \frac{P(D | h)P(h)}{P(D)} \] Donde \(P(h)\) es la probabilidad de que la hipótesis sea cierta (independiente de los datos), también conocida como la probabilidad previa de \(h\)

\(P(D)\) es la probabilidad de los datos independientemente de la hipótesis. \(P(h|D)\) es la probabilidad de la hipótesis \(h\) dados los datos \(D\). Probabilidad posterior. \(P(D|h)\) es la probabilidad de los datos \(D\) dado que la hipótesis era cierta (Probabilidad posterior)

El problema de los datos infrecuentes

Al ser una multiplicación de probabilidades, si los datos históricos no tienen una combinación de eventos dado que son infrecuentes (no ocurrieron en el set de datos de prueba), la multiplicación total será cero, y llevará a que toda la predicción indique que hay cero por ciento de probabilidad de que el evento ocurra.

Para solucionar esto se puede usar la correción de laplace, que consiste en sumar un valor (1% de probabilidad normalmente) para evitar que la predicción llegue a a ser cero en caso que un evento sea improbable

locations <- read.csv("locations.csv")

nb <- naive_bayes(location ~ hourtype, data = locations)

future_location <- predict(nb, future_condition)

print(locmodel)

# podemos predecir la probabilidad de estar en cierto lugar un sabado a la 9 ma

predict(locmodel, saturday9am , type = "prob")

Con corrección de laplace…

# En este caso nos da cero
predict(locmodel, weekend_afternoon, type = "prob")

# Usamos la corrección de laplace para datos infrecuentes
locmodel2 <- naive_bayes(location ~ daytype + hourtype , data= locations, laplace = 1)

# En la corrección la probabilidad ahora es más precisa
predict(locmodel2, weekend_afternoon, type = "prob")

Naive Bayes trabaja muy bien con variables predictoras categoricas, por ellos es una buena idea usar “binning” (Apilar datos en rangos, categorías que tengan sentido, percentiles, etc) sobre datos numericos

Regresión logística

La regresión logística es un modelo de probabilidad no lineal estimado con máxima verosimilitud, que busca estimar los parametros \(\beta\) que responden a una probabilidad entre 0 y 1. La función que usamos para predecir:

\[ \rho(Z) = \frac{e^( \beta_{0} + \beta_{1}X_{1} + ... + \beta_{p}X_{p} )}{1 + e^( \beta_{0} + \beta_{1}X_{1} + ... + \beta_{p}X_{p} )} \]

Odds ratio

Podemos demás usar los coeficientes para hacer comparaciones en los efectos de las variables sobre la probabilidad:

  • Si el coeficiente \(\beta\) es positivo: \(e^\beta\) (es \(odds_ratio\) veces más probable que ocurra el evento si aumenta en 1 o es de x categoria la variable)

  • Si el coficiente \(\beta\) es negativo: \(\frac{1}{e^\beta}\) (es \(odds_ratio\) veces menos probable que ocurra el evento si aumenta en 1 o es de x categoria la variable)

Podemos examinar los resultados de donaciones de las personas:

donors <- read.csv("donors.csv")
str(donors)
## 'data.frame':    93462 obs. of  13 variables:
##  $ donated          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ veteran          : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ bad_address      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ age              : int  60 46 NA 70 78 NA 38 NA NA 65 ...
##  $ has_children     : int  0 1 0 0 1 0 1 0 0 0 ...
##  $ wealth_rating    : int  0 3 1 2 1 0 2 3 1 0 ...
##  $ interest_veterans: int  0 0 0 0 0 0 0 0 0 0 ...
##  $ interest_religion: int  0 0 0 0 1 0 0 0 0 0 ...
##  $ pet_owner        : int  0 0 0 0 0 0 1 0 0 0 ...
##  $ catalog_shopper  : int  0 0 0 0 1 0 0 0 0 0 ...
##  $ recency          : Factor w/ 2 levels "CURRENT","LAPSED": 1 1 1 1 1 1 1 1 1 1 ...
##  $ frequency        : Factor w/ 2 levels "FREQUENT","INFREQUENT": 1 1 1 1 1 2 2 1 2 2 ...
##  $ money            : Factor w/ 2 levels "HIGH","MEDIUM": 2 1 2 2 2 2 2 2 2 2 ...
table(donors$donated)
## 
##     0     1 
## 88751  4711
# Modelo de regresión logística
donation_model <- glm(donated ~ bad_address + interest_religion +interest_veterans , 
                      data = donors, family = "binomial")


summary(donation_model)
## 
## Call:
## glm(formula = donated ~ bad_address + interest_religion + interest_veterans, 
##     family = "binomial", data = donors)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.3480  -0.3192  -0.3192  -0.3192   2.5678  
## 
## Coefficients:
##                   Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)       -2.95139    0.01652 -178.664   <2e-16 ***
## bad_address       -0.30780    0.14348   -2.145   0.0319 *  
## interest_religion  0.06724    0.05069    1.327   0.1847    
## interest_veterans  0.11009    0.04676    2.354   0.0186 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37330  on 93461  degrees of freedom
## Residual deviance: 37316  on 93458  degrees of freedom
## AIC: 37324
## 
## Number of Fisher Scoring iterations: 5
# Estimamos las probabilidades
donors$donation_prob <- predict(donation_model, type = "response")

#Buscamos la probabilidad promedio del prospecto donador para estableces un umbral
mean(donors$donated)
## [1] 0.05040551
# Clasificamos a las personas basados en el umbral que fijamos (proporción de donantes de los datos históricos)
donors$donation_pred <- ifelse(donors$donation_prob  > 0.0504, 1, 0)

# Calculamos la presición del modelo
mean(donors$donation_pred == donors$donated)
## [1] 0.794815

Un gran reto que tienen los modelos de clasificación es que en algunos casos el evento de interés (en el ejemplo anterior es si la persona donó o no) son muy raros, por lo que es natural que las predicciones sean muy altas, ya que se predicen bien los resultados negativos (no donó), pero no se está revisando con precisión si lo hace igual o mejor con resultados positivos

La curva ROC (Receive operating curve) es excelente para distinguir entre predicciones positivas (1) o negativas (0) del resultado esperado

  • El modelo de linea base que no es mejor que uno aleatorio tiene un AUC (Area under de curve) de 0.5

  • Un modelo perfeceto tiene un AUC de 1

Importante tener en cuenta que curvas ROC pueden tener el mismo AUC pero diferente forma, por lo que hay que prestar atención cómo están siendo las predicciones a lo largo de un rango de la curva (más verdaderos positivos, más verdaderos negativos)

# Creamos la curva ROC del modelo (objecto)
ROC <- roc(donors$donated, donors$donation_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(ROC, col = "blue")

# AUC
auc(ROC)
## Area under the curve: 0.5102

Por defecto la función glm omite las observaciones con NA en alguno de los predictores, por lo que hay que tener cuidado la selección tanto de observaciones como de variables.

Aunque algunas veces los NA puede representar incluso una variable en si misma!

  • **Efector interacción*+ Multiplicación entre variables predictoras que buscaría amplificar, disminuir o anular los efectos de las variables
donors$wealth_levels <- factor(donors$wealth_rating, levels = c(0, 1, 2, 3), labels =c("Unknown", "Low", "Medium", "High"))


donors$wealth_levels <- relevel(donors$wealth_levels, ref = "Medium")
summary(glm(donated ~ wealth_levels, family = 'binomial', data = donors))
## 
## Call:
## glm(formula = donated ~ wealth_levels, family = "binomial", data = donors)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.3320  -0.3243  -0.3175  -0.3175   2.4582  
## 
## Coefficients:
##                      Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -2.91894    0.03614 -80.772   <2e-16 ***
## wealth_levelsUnknown -0.04373    0.04243  -1.031    0.303    
## wealth_levelsLow     -0.05245    0.05332  -0.984    0.325    
## wealth_levelsHigh     0.04804    0.04768   1.008    0.314    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37330  on 93461  degrees of freedom
## Residual deviance: 37323  on 93458  degrees of freedom
## AIC: 37331
## 
## Number of Fisher Scoring iterations: 5
summary(donors$age)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##    1.00   48.00   62.00   61.65   75.00   98.00   22546
donors$imputed_age <- ifelse(is.na(donors$age), round(mean(donors$age, na.rm = T),2), donors$age)

donors$missing_age <- ifelse(is.na(donors$age),1, 0)
rfm_model <- glm(donated ~ money + recency * frequency, data = donors, family = 'binomial')

summary(rfm_model)
## 
## Call:
## glm(formula = donated ~ money + recency * frequency, family = "binomial", 
##     data = donors)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -0.3696  -0.3696  -0.2895  -0.2895   2.7924  
## 
## Coefficients:
##                                   Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                       -3.01142    0.04279 -70.375   <2e-16 ***
## moneyMEDIUM                        0.36186    0.04300   8.415   <2e-16 ***
## recencyLAPSED                     -0.86677    0.41434  -2.092   0.0364 *  
## frequencyINFREQUENT               -0.50148    0.03107 -16.143   <2e-16 ***
## recencyLAPSED:frequencyINFREQUENT  1.01787    0.51713   1.968   0.0490 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 37330  on 93461  degrees of freedom
## Residual deviance: 36938  on 93457  degrees of freedom
## AIC: 36948
## 
## Number of Fisher Scoring iterations: 6
rfm_prob <- predict(rfm_model, type = "response")

ROC <- roc(donors$donated, rfm_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(ROC, col = "red")

auc(ROC)
## Area under the curve: 0.5785

El modelo mejoró levemente dada la experticie e intución implementada en el modelo.

Regresión Stepwise

Podemos construir el modelo con menos experticie y más iteraciones y decisiones conducidas por los datos. Este tipo de modelos hacen un análisis hacia atras y hacia adelante con todas las variables hasta encontrar la mejor combinación de predictores que general el modelo más potente.

# Modelo nulo sin predictores
null_model <- glm(donated ~ 1, data = donors, family = "binomial")

# Modelo con todos los predictores
full_model <- glm(donated ~ ., data = donors, family = "binomial")

# Usamos la función step() hacia adelante para empezar a incluir variables en el modelo y mirar como se comporta
step_model <- step(null_model, scope = list(lower = null_model, upper = full_model), direction = "forward")
## Start:  AIC=37332.13
## donated ~ 1
## Warning in add1.glm(fit, scope$add, scale = scale, trace = trace, k = k, : using
## the 70916/93462 rows from a combined fit
##                     Df Deviance   AIC
## + frequency          1    28502 37122
## + money              1    28621 37241
## + wealth_rating      1    28705 37326
## + has_children       1    28705 37326
## + age                1    28707 37328
## + imputed_age        1    28707 37328
## + wealth_levels      3    28704 37328
## + interest_veterans  1    28709 37330
## + donation_prob      1    28710 37330
## + donation_pred      1    28710 37330
## + catalog_shopper    1    28710 37330
## + pet_owner          1    28711 37331
## <none>                    28714 37332
## + interest_religion  1    28712 37333
## + recency            1    28713 37333
## + bad_address        1    28714 37334
## + veteran            1    28714 37334
## 
## Step:  AIC=37024.77
## donated ~ frequency
## Warning in add1.glm(fit, scope$add, scale = scale, trace = trace, k = k, : using
## the 70916/93462 rows from a combined fit
##                     Df Deviance   AIC
## + money              1    28441 36966
## + wealth_rating      1    28493 37018
## + wealth_levels      3    28490 37019
## + has_children       1    28494 37019
## + donation_prob      1    28498 37023
## + interest_veterans  1    28498 37023
## + catalog_shopper    1    28499 37024
## + donation_pred      1    28499 37024
## + age                1    28499 37024
## + imputed_age        1    28499 37024
## + pet_owner          1    28499 37024
## <none>                    28502 37025
## + interest_religion  1    28501 37026
## + recency            1    28501 37026
## + bad_address        1    28502 37026
## + veteran            1    28502 37027
## 
## Step:  AIC=36949.71
## donated ~ frequency + money
## Warning in add1.glm(fit, scope$add, scale = scale, trace = trace, k = k, : using
## the 70916/93462 rows from a combined fit
##                     Df Deviance   AIC
## + wealth_levels      3    28427 36942
## + wealth_rating      1    28431 36942
## + has_children       1    28432 36943
## + interest_veterans  1    28438 36948
## + donation_prob      1    28438 36949
## + catalog_shopper    1    28438 36949
## + donation_pred      1    28439 36949
## + age                1    28439 36949
## + imputed_age        1    28439 36949
## + pet_owner          1    28439 36949
## <none>                    28441 36950
## + interest_religion  1    28440 36951
## + recency            1    28441 36951
## + bad_address        1    28441 36951
## + veteran            1    28441 36952
## 
## Step:  AIC=36945.48
## donated ~ frequency + money + wealth_levels
## Warning in add1.glm(fit, scope$add, scale = scale, trace = trace, k = k, : using
## the 70916/93462 rows from a combined fit
##                     Df Deviance   AIC
## + has_children       1    28416 36937
## + age                1    28424 36944
## + imputed_age        1    28424 36944
## + interest_veterans  1    28424 36945
## + donation_prob      1    28424 36945
## + catalog_shopper    1    28425 36945
## + donation_pred      1    28425 36945
## <none>                    28427 36945
## + pet_owner          1    28425 36946
## + interest_religion  1    28426 36947
## + recency            1    28427 36947
## + bad_address        1    28427 36947
## + veteran            1    28427 36947
## 
## Step:  AIC=36938.4
## donated ~ frequency + money + wealth_levels + has_children
## Warning in add1.glm(fit, scope$add, scale = scale, trace = trace, k = k, : using
## the 70916/93462 rows from a combined fit
##                     Df Deviance   AIC
## + pet_owner          1    28413 36937
## + donation_prob      1    28413 36937
## + catalog_shopper    1    28413 36937
## + interest_veterans  1    28413 36937
## + donation_pred      1    28414 36938
## <none>                    28416 36938
## + interest_religion  1    28415 36939
## + age                1    28416 36940
## + imputed_age        1    28416 36940
## + recency            1    28416 36940
## + bad_address        1    28416 36940
## + veteran            1    28416 36940
## 
## Step:  AIC=36932.25
## donated ~ frequency + money + wealth_levels + has_children + 
##     pet_owner
## Warning in add1.glm(fit, scope$add, scale = scale, trace = trace, k = k, : using
## the 70916/93462 rows from a combined fit
##                     Df Deviance   AIC
## <none>                    28413 36932
## + donation_prob      1    28411 36932
## + interest_veterans  1    28411 36932
## + catalog_shopper    1    28412 36933
## + donation_pred      1    28412 36933
## + age                1    28412 36933
## + imputed_age        1    28412 36933
## + recency            1    28413 36934
## + interest_religion  1    28413 36934
## + bad_address        1    28413 36934
## + veteran            1    28413 36934
# PRedecimos de la misma forma que con glm()
step_prob <- predict(step_model, type = 'response')

ROC <- roc(donors$donated, step_prob)
## Setting levels: control = 0, case = 1
## Setting direction: controls < cases
plot(ROC, col = "red")

auc(ROC)
## Area under the curve: 0.5849

Arboles de clasificación

Los arboles de decisión son modelos bastante simples que realizan particiones para llegar a una decisión sobre qué valor o clasificación tendrán nuevas predicciones.

El primero split será el más “puro” o el que cree las particiones más homogeneas, y asi de manera sucesiva hasta que se hagan los splits necesarios o deseados.

loans <- read.csv("loans.csv")

loan_model <- rpart(keep ~ ., data = loans, method = "class", control = rpart.control(cp = 0))


predictions <- predict(loan_model, loans[, -1], type = "class")


summary(loan_model)
## Call:
## rpart(formula = keep ~ ., data = loans, method = "class", control = rpart.control(cp = 0))
##   n= 39732 
## 
##          CP nsplit rel error       xerror         xstd
## 1 0.5835396      0 1.0000000 1.0000000000 7.951922e-03
## 2 0.4164604      1 0.4164604 0.4165487977 5.697064e-03
## 3 0.0000000      2 0.0000000 0.0000884017 8.840058e-05
## 
## Variable importance
##    rand default 
##      50      50 
## 
## Node number 1: 39732 observations,    complexity param=0.5835396
##   predicted class=0  expected loss=0.2847075  P(node) =1
##     class counts: 28420 11312
##    probabilities: 0.715 0.285 
##   left son=2 (33131 obs) right son=3 (6601 obs)
##   Primary splits:
##       rand           < 0.1666702 to the right, improve=8100.52000, (0 missing)
##       default        < 0.5       to the left,  improve=6745.58200, (0 missing)
##       credit_score   splits as  RLR, improve=  59.89707, (0 missing)
##       loan_purpose   splits as  LLRRLRLRRRLRRL, improve=  44.57217, (0 missing)
##       recent_inquiry splits as  LR, improve=  41.13932, (0 missing)
## 
## Node number 2: 33131 observations,    complexity param=0.4164604
##   predicted class=0  expected loss=0.1421931  P(node) =0.8338619
##     class counts: 28420  4711
##    probabilities: 0.858 0.142 
##   left son=4 (28420 obs) right son=5 (4711 obs)
##   Primary splits:
##       default            < 0.5       to the left,  improve=8082.25600, (0 missing)
##       credit_score       splits as  RLR, improve=  71.93522, (0 missing)
##       loan_purpose       splits as  LLLLLLLLLLLRLL, improve=  48.46663, (0 missing)
##       credit_utilization splits as  RLL, improve=  41.99130, (0 missing)
##       recent_inquiry     splits as  LR, improve=  31.79880, (0 missing)
## 
## Node number 3: 6601 observations
##   predicted class=1  expected loss=0  P(node) =0.1661381
##     class counts:     0  6601
##    probabilities: 0.000 1.000 
## 
## Node number 4: 28420 observations
##   predicted class=0  expected loss=0  P(node) =0.7152925
##     class counts: 28420     0
##    probabilities: 1.000 0.000 
## 
## Node number 5: 4711 observations
##   predicted class=1  expected loss=0  P(node) =0.1185694
##     class counts:     0  4711
##    probabilities: 0.000 1.000
rpart.plot(loan_model)

rpart.plot(loan_model, type = 3, box.palette = c("red", "green"), fallen.leaves = TRUE)

Lo verdaderamente correcto en cualqueir tarea de regresión o clasificación es dividir los datos en entrenamiento y prueba, de manera que tengamos un set de datos que el algoritmo entrenado no ha visto y que podemos evaluar de manera fiel los resultados de las predicción con alguna medida de erroe deseada.

# Partimos los datos en un 75-25%
nrow(loans) * 0.75
## [1] 29799
# Creamos una partición
sample_rows <- sample(11312, 8484)

#usamos los indices alatorios en los datos
loans_train <- loans[sample_rows,]

# Cremoa el datset de prueba
loans_test <- loans[-sample_rows,]


# Modelo sencillo de arbol 
loan_model <- rpart(keep ~ ., data = loans_train, method = "class", control = rpart.control(cp = 0))

# predecimos las clasificaiones
loans_test$pred <- predict(loan_model, loans_test, type = 'class')

# REvismos la matriz de confusión
table(loans_test$keep, loans_test$pred)
##    
##         0     1
##   0 22406    19
##   1     0  8823
# Calculamos la presición promedio
mean(loans_test$keep == loans_test$pred)
## [1] 0.999392

Tambien podemos modificar los parametros de manera que consigamos optimizar el modelo

loans_test <- loans[-sample_rows,]
loan_model <- rpart(keep ~ ., data = loans_train,control =  rpart.control(cp = 0, maxdepth = 6, minsplit = 500))

loans_test$pred <- predict(loan_model, loans_test, type = "class")
mean(loans_test$pred == loans_test$keep)
loan_model <- rpart(keep ~ ., data  = loans_train, control = rpart.control(cp = 0), method = "class")


plotcp(loan_model)

loan_model_pruned <- prune(loan_model, cp = 0.0014)


loans_test$pred <- predict(loan_model_pruned, loans_test, type = 'class')
mean(loans_test$pred == loans_test$keep, na.rm = T)
## [1] 0.999392