Descriminante Lineal

El Análisis Discriminante Lineal o Linear Discrimiant Analysis (LDA) es un método de clasificación de variables cualitativas en el que dos o más grupos son conocidos a priori y nuevas observaciones se clasifican en uno de ellos en función de sus características. Haciendo uso del teorema de Bayes, LDA estima la probabilidad de que una observación, dado un determinado valor de los predictores, pertenezca a cada una de las clases de la variable cualitativa, P(Y=k|X=x). Finalmente se asigna la observación a la clase k para la que la probabilidad predicha es mayor.

Es una alternativa a la regresión logística cuando la variable cualitativa tiene más de dos niveles. Si bien existen extensiones de la regresión logística para múltiples clases, el LDA presenta una serie de ventajas:

Si las clases están bien separadas, los parámetros estimados en el modelo de regresión logística son inestables. El método de LDA no sufre este problema. Si el número de observaciones es bajo y la distribución de los predictores es aproximadamente normal en cada una de las clases, LDA es más estable que la regresión logística.

Cuando se trata de un problema de clasificación con solo dos niveles, ambos métodos suelen llegar a resultados similares.

Condiciones necesarias para el LDA

Las condiciones que se deben cumplir para que un Análisis Discriminante Lineal sea válido son:

Precisión del LDA

Una vez que las normas de clasificación se han establecido, se tiene que evaluar como de buena es la clasificación resultante. En otras palabras, evaluar el porcentaje de aciertos en las clasificaciones.

Las matrices de confusión son una de las mejores formas de evaluar la capacidad de acierto que tiene un modelo LDA. Muestran el número de verdaderos positivos, verdaderos negativos, falsos positivos y falsos negativos. El método LDA busca los límites de decisión que más se aproximan al clasificador de Bayes, que por definición, tiene el menor ratio de error total de entre todos los clasificadores (si se cumple la condición de normalidad). Por lo tanto, el LDA intenta conseguir el menor número de clasificaciones erróneas posibles, pero no diferencia entre falsos positivos o falsos negativos. Si se quiere intentar reducir el número de errores de clasificación en una dirección determinada (por ejemplo, menos falsos negativos) se puede modificar el límite de decisión, aunque como consecuencia aumentará el número de falsos positivos.

Cuando para evaluar el error de clasificación se emplean las mismas observaciones con las que se ha creado el modelo, se obtiene lo que se denomina el training error. Si bien esta es una forma sencilla de estimar la precisión en la clasificación, tiende a ser excesivamente optimista. Es más adecuado evaluar el modelo empleando observaciones nuevas que el modelo no ha visto, obteniendo así el test error. En el capítulo Validación de modelos de regresión se describen diferentes estrategias para estimar el test error.

Ejemplo Se desea utilizar el conjunto de datos iris para realizar una analisis LDA

normalidad<-function(var){
  shapiro.test(var)$p
}

# Verificamos normalidad para cada grupo 
iris %>% 
  group_by(Species) %>% 
  summarise_all(normalidad)->tb1
kable(tb1)
Species Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa 0.4595132 0.2715264 0.0548115 0.0000009
versicolor 0.4647370 0.3379951 0.1584778 0.0272778
virginica 0.2583147 0.1808960 0.1097754 0.0869542 
####Normalidad Multivariada
library(MVN)

mvn(iris[-5],mvnTest = "mardia")
## $multivariateNormality
##              Test          Statistic              p value Result
## 1 Mardia Skewness    67.430508778062 4.75799820400869e-07     NO
## 2 Mardia Kurtosis -0.230112114481001    0.818004651478012    YES
## 3             MVN               <NA>                 <NA>     NO
## 
## $univariateNormality
##           Test     Variable Statistic   p value Normality
## 1 Shapiro-Wilk Sepal.Length    0.9761  0.0102      NO    
## 2 Shapiro-Wilk Sepal.Width     0.9849  0.1012      YES   
## 3 Shapiro-Wilk Petal.Length    0.8763  <0.001      NO    
## 4 Shapiro-Wilk Petal.Width     0.9018  <0.001      NO    
## 
## $Descriptives
##                n     Mean   Std.Dev Median Min Max 25th 75th       Skew
## Sepal.Length 150 5.843333 0.8280661   5.80 4.3 7.9  5.1  6.4  0.3086407
## Sepal.Width  150 3.057333 0.4358663   3.00 2.0 4.4  2.8  3.3  0.3126147
## Petal.Length 150 3.758000 1.7652982   4.35 1.0 6.9  1.6  5.1 -0.2694109
## Petal.Width  150 1.199333 0.7622377   1.30 0.1 2.5  0.3  1.8 -0.1009166
##                Kurtosis
## Sepal.Length -0.6058125
## Sepal.Width   0.1387047
## Petal.Length -1.4168574
## Petal.Width  -1.3581792
# prueba de box

library(biotools)
## ---
## biotools version 3.1
# Ho: matriz de covarianza es constante
boxM(data = iris[, -5], grouping = iris[, 5])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  iris[, -5]
## Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16
# ser rechaza  la Ho

# LDA
modelo<-lda(Species~.,data=iris)
modelo
## Call:
## lda(Species ~ ., data = iris)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa            5.006       3.428        1.462       0.246
## versicolor        5.936       2.770        4.260       1.326
## virginica         6.588       2.974        5.552       2.026
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Length  0.8293776  0.02410215
## Sepal.Width   1.5344731  2.16452123
## Petal.Length -2.2012117 -0.93192121
## Petal.Width  -2.8104603  2.83918785
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9912 0.0088
prediccion<-predict(modelo,iris[-5])
m_confusion<-table(iris$Species,prediccion$class,
                   dnn=c("Real","Predicho"))
m_confusion
##             Predicho
## Real         setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         2
##   virginica       0          1        49
mosaicplot(m_confusion,col=2:4)

precision=mean(iris$Species==prediccion$class)
precision
## [1] 0.98
error=1-precision
error
## [1] 0.02

Ejemplo QDA

library(MASS)
modelo2<-qda(Species~.,data=iris)
prediccion2<-predict(modelo2,iris[-5])
m_confusion2<- table(iris$Species,prediccion2$class,
                     dnn=c("Real","predicho"))
m_confusion2
##             predicho
## Real         setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         2
##   virginica       0          1        49
precision2=mean(iris$Species==prediccion2$class)
error2=1-precision2
precision2
## [1] 0.98
error2
## [1] 0.02

Gráficos de la clasificación

library(klaR)
partimat(Species~.,
         data=iris,
         method="lda",
         image.colors=c("orange","yellow","salmon"),
         col.mean="blue")

partimat(Species~.,
         data=iris,
         method="qda",
         image.colors=c("orange","yellow","salmon"),
         col.mean="blue")

Regresión Logística

La regresión logística forma parte de los llamados modelos lineales generalizados y trata de estimar la probabilidad de ocurrencia de un evento binario (éxito/fracaso, cara/cruz) en función de una serie de variables predictoras.

Es importante tener en cuenta que, aunque la regresión logística permite clasificar, se trata de un modelo de regresión que modela el logaritmo de la probabilidad de pertenecer a cada grupo. La asignación final se hace en función de las probabilidades predichas.

Condiciones

  • Independencia: las observaciones tienen que ser independientes unas de otras.

  • Relación lineal entre el logaritmo natural de odds y la variable continua: patrones en forma de U son una clara violación de esta condición.

  • La regresión logística no precisa de una distribución normal de la variable continua independiente.

  • Número de observaciones: no existe una norma establecida al respecto, pero se recomienda entre 50 a 100 observaciones.

Ejemplo de un modelo

EL siguiente conjunto de datos se quiere modelar la variable donated (donado) en funcion de las variables bad_address(variable que indica si la dirección de correo esta incorrecta) interest_religion(interes por la religión), interest_veterans(interes por los veteranos)

#modelo 1  
data<-read.csv("donacion.csv")
library(dplyr)
glimpse(data)
## Observations: 93,462
## Variables: 14
## $ X                 <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 1...
## $ donated           <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ veteran           <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ bad_address       <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ age               <int> 60, 46, NA, 70, 78, NA, 38, NA, NA, 65, NA, ...
## $ has_children      <int> 0, 1, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0,...
## $ wealth_rating     <int> 0, 3, 1, 2, 1, 0, 2, 3, 1, 0, 1, 2, 1, 0, 2,...
## $ interest_veterans <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ interest_religion <int> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ pet_owner         <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ catalog_shopper   <int> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,...
## $ recency           <fct> CURRENT, CURRENT, CURRENT, CURRENT, CURRENT,...
## $ frequency         <fct> FREQUENT, FREQUENT, FREQUENT, FREQUENT, FREQ...
## $ money             <fct> MEDIUM, HIGH, MEDIUM, MEDIUM, MEDIUM, MEDIUM...
# explorando variable respuesta
prop.table(table(data$donated))
## 
##          0          1 
## 0.94959449 0.05040551
# fabricando el modelo
mod1 <-glm(donated~interest_religion+
             interest_veterans+
             bad_address,
           data=data,
           family = binomial)

#calculando las  variable respuesta 
pred_p<-predict(mod1,type="response")

# convertimos la probabilidad de clasificación
pred<-ifelse(pred_p>0.05,1,0)
sum(pred)
## [1] 16258
# evaluamos la precisión del modelo

mean(data$donated==pred)
## [1] 0.794815
table(data$donated,pred,
      dnn=c("actual","predicho")) -> conf


mosaicplot(conf,col=2:4)

library(pROC)
roc(data$donated,pred) -> ROC
plot(ROC)

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

##Ejemplo 2 (modelo sin valores perdidos) y stepwise

data2<-na.omit(data)
modelo_nulo<-glm(donated~1,data=data2,
                family = binomial)
modelo_completo<-glm(donated~.,
                     data=data2,
                     family=binomial)
step(modelo_nulo, scope = list(lower = modelo_nulo, upper = modelo_completo), direction = "forward")->modelo_paso
## Start:  AIC=28715.77
## donated ~ 1
## 
##                     Df Deviance   AIC
## + frequency          1    28502 28506
## + money              1    28621 28625
## + wealth_rating      1    28705 28709
## + has_children       1    28705 28709
## + age                1    28707 28711
## + interest_veterans  1    28709 28713
## + catalog_shopper    1    28710 28714
## + pet_owner          1    28711 28715
## <none>                    28714 28716
## + interest_religion  1    28712 28716
## + recency            1    28713 28717
## + bad_address        1    28714 28718
## + X                  1    28714 28718
## + veteran            1    28714 28718
## 
## Step:  AIC=28505.88
## donated ~ frequency
## 
##                     Df Deviance   AIC
## + money              1    28441 28447
## + wealth_rating      1    28493 28499
## + has_children       1    28494 28500
## + interest_veterans  1    28498 28504
## + catalog_shopper    1    28499 28505
## + age                1    28499 28505
## + pet_owner          1    28499 28505
## <none>                    28502 28506
## + interest_religion  1    28501 28507
## + recency            1    28501 28507
## + bad_address        1    28502 28508
## + X                  1    28502 28508
## + veteran            1    28502 28508
## 
## Step:  AIC=28446.96
## donated ~ frequency + money
## 
##                     Df Deviance   AIC
## + wealth_rating      1    28431 28439
## + has_children       1    28432 28440
## + interest_veterans  1    28438 28446
## + catalog_shopper    1    28438 28446
## + age                1    28439 28447
## + pet_owner          1    28439 28447
## <none>                    28441 28447
## + interest_religion  1    28440 28448
## + recency            1    28441 28449
## + bad_address        1    28441 28449
## + veteran            1    28441 28449
## + X                  1    28441 28449
## 
## Step:  AIC=28439.43
## donated ~ frequency + money + wealth_rating
## 
##                     Df Deviance   AIC
## + has_children       1    28421 28431
## + interest_veterans  1    28429 28439
## + catalog_shopper    1    28429 28439
## + age                1    28429 28439
## <none>                    28431 28439
## + pet_owner          1    28430 28440
## + interest_religion  1    28431 28441
## + recency            1    28431 28441
## + bad_address        1    28431 28441
## + veteran            1    28431 28441
## + X                  1    28431 28441
## 
## Step:  AIC=28431.06
## donated ~ frequency + money + wealth_rating + has_children
## 
##                     Df Deviance   AIC
## + pet_owner          1    28418 28430
## + catalog_shopper    1    28418 28430
## + interest_veterans  1    28418 28430
## <none>                    28421 28431
## + interest_religion  1    28420 28432
## + recency            1    28421 28433
## + age                1    28421 28433
## + bad_address        1    28421 28433
## + X                  1    28421 28433
## + veteran            1    28421 28433
## 
## Step:  AIC=28429.73
## donated ~ frequency + money + wealth_rating + has_children + 
##     pet_owner
## 
##                     Df Deviance   AIC
## <none>                    28418 28430
## + interest_veterans  1    28416 28430
## + catalog_shopper    1    28416 28430
## + age                1    28417 28431
## + recency            1    28417 28431
## + interest_religion  1    28417 28431
## + bad_address        1    28418 28432
## + X                  1    28418 28432
## + veteran            1    28418 28432
predict(modelo_paso,type="response")->prob

roc(data2$donated,prob)->ROC


auc(ROC)
## Area under the curve: 0.5824
plot(ROC)

Ejercicio:

Resolver el ejemplo anterior, haciendo las siguientes modificaciones.

  • hacer inputación de la media, para los valores perdidos de la variable age

  • construir una nueva variable indicadora de que se efectuó inputación a de la media en los valores perdidos, 1=se efectuó imputación, 0= no se efectuó imputación

  • Después de hacer esto, eliminar la variable age

  • Codificar la variable wealt_rating, para que sea tipo factor de la siguiente manera: 0= Unknown, 1=Low, 2=Medium, 3= High, guarde el resultado en data$wealth_levels

  • codificar la variable anterior para que se un factor ordenado, Medium, Unknown, Low, High, utilice relevel(data$wealth_levels, ref = “Medium”)

Bibliografía