Análisis Discriminante
Pamela E. Pairo
Es una técnica de clasificación en donde la variable respuesta es cuantitativa.
El objetivo es:
Separar los grupos, identificando las variables que permiten esa separación
Clasificar a nuevos individuos
La idea es contruir funciones discriminantes (FD) que maximicen la variabilidad entre grupos con el objetivo de discriminarlos mejor. Esto se logra al maximizar la varianza entre grupos en relación con la varianza total
Base de datos
Vamos a trabajar con la base de datos de Iris
data(iris)
glimpse(iris)## Rows: 150
## Columns: 5
## $ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9, 5.4, 4.~
## $ Sepal.Width <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1, 3.7, 3.~
## $ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5, 1.5, 1.~
## $ Petal.Width <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1, 0.2, 0.~
## $ Species <fct> setosa, setosa, setosa, setosa, setosa, setosa, setosa, s~Realicemos algunos scatterplots
library(GGally)
ggpairs(iris, legend = 1, columns = 1:4, aes(color = Species, alpha = 0.5),
upper = list(continuous = "points"))+
theme(legend.position = "bottom")
Al conocerse la especie de pertenencia de cada individuo, se puede comprobar la efectividad del método de clasificación observando el porcentaje de casos bien clasificados. Para ello, vamos a separar el dataset en train y test
set.seed(123)#setear la semilla
# Create data split for train and test
df_split <- initial_split(iris,
prop = 0.9,
strata = Species)#para conservar la proporción de las clases
df_train <- df_split %>%
training()
df_test <- df_split %>%
testing()
# Número de datos en test y train
paste0("Total del dataset de entrenamiento: ", nrow(df_train))## [1] "Total del dataset de entrenamiento: 135"paste0("Total del dataset de testeo: ", nrow(df_test))## [1] "Total del dataset de testeo: 15"Creamos tres subsets de datos para cada especie
setosa<- subset(df_train[,1:4], df_train$Species == "setosa")
versicolor<- subset(df_train[,1:4], df_train$Species == "versicolor")
virginica <- subset(df_train[,1:4], df_train$Species=="virginica")Análisis Discriminante Lineal (LDA)
Este tipo de análisis es válido solo si se satisfacen los siguientes supuestos:
1- Normalidad multivariada
2- Independencia de las observaciones
3- Homocedasticidad.
Verefiquemos supuestos
1- Normalidad multivariada
¿Cuál es la hipótesis estadística?
mvShapiro.Test(as.matrix(setosa))##
## Generalized Shapiro-Wilk test for Multivariate Normality by
## Villasenor-Alva and Gonzalez-Estrada
##
## data: as.matrix(setosa)
## MVW = 0.96115, p-value = 0.0371mvShapiro.Test(as.matrix(versicolor)) ##
## Generalized Shapiro-Wilk test for Multivariate Normality by
## Villasenor-Alva and Gonzalez-Estrada
##
## data: as.matrix(versicolor)
## MVW = 0.97362, p-value = 0.4165mvShapiro.Test(as.matrix(virginica))##
## Generalized Shapiro-Wilk test for Multivariate Normality by
## Villasenor-Alva and Gonzalez-Estrada
##
## data: as.matrix(virginica)
## MVW = 0.98463, p-value = 0.97642- Independencia de las observaciones
Viene dada por el diseño
3- Homocedasticidad
Ho las matrices de varianzas-covarianzas de los grupos son iguales.
Analizamos igualdad de matrices de varianzas y covarianzas:
boxM(df_train[,-5],df_train[,5])##
## Box's M-test for Homogeneity of Covariance Matrices
##
## data: df_train[, -5]
## Chi-Sq (approx.) = 136.76, df = 20, p-value < 2.2e-16Entonces ¿qué concluimos?
Vamos a proseguir como si se hubiese cumplido todos los supuestos. Tener en cuenta que los resultados del LDA no van a ser confiables en este caso.
LDA
model_lda <- lda(Species~., data =df_train)
model_lda## Call:
## lda(Species ~ ., data = df_train)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 4.991111 3.393333 1.471111 0.2377778
## versicolor 5.940000 2.757778 4.262222 1.3266667
## virginica 6.566667 2.977778 5.506667 2.0133333
##
## Coefficients of linear discriminants:
## LD1 LD2
## Sepal.Length 0.7633602 -0.2360469
## Sepal.Width 1.8553244 2.4183855
## Petal.Length -2.2547164 -0.7663631
## Petal.Width -2.9175526 2.6999226
##
## Proportion of trace:
## LD1 LD2
## 0.9913 0.0087El objeto lda (en este ejemplo, model_lda) contiene los siguientes componentes:
prior: las probabilidades previas utilizadas.
means: la media de cada clase.
svd: son los valores singulares, informan el cociente entre desvíos estándar entre y dentro de grupos para cada FD; si se elevan al cuadrado se obtiene el autovalor de cada FD
counts: El número de observaciones por clase.
Coefficients of linear discriminants: Muestra la combinación lineal de variables predictoras que se utilizan para formar la regla de decisión LDA.
Por ejemplo:
\(\ LD1= 0.76*Sepal.Length+ 1.85*Sepal.Width - 2.25*Petal.Length -2.91*Petal.Width\)
LD1 explica el 99% de la proporción de varianza entre clases.
La 2da función discriminante es independiente de la primera (ortogonal) y es la que mejor separa los grupos usando la variación remanente o residual, después que la 1ra función discriminante ha sido determinada
prop = model_lda$svd^2/sum(model_lda$svd^2)
prop #varianza entre grupos explicada por cada FD## [1] 0.991319021 0.008680979lda.data <- cbind(df_train, predict(model_lda)$x)
ggplot(lda.data, aes(LD1, LD2)) +
geom_point(aes(color = Species))
Veamos que predice el modelo en df_test
predictions <- model_lda %>% predict(df_test)
predictions## $class
## [1] setosa setosa setosa setosa setosa versicolor
## [7] versicolor versicolor versicolor versicolor virginica virginica
## [13] virginica virginica virginica
## Levels: setosa versicolor virginica
##
## $posterior
## setosa versicolor virginica
## 2 1.000000e+00 1.249522e-19 2.400232e-39
## 16 1.000000e+00 1.695972e-31 2.933345e-53
## 23 1.000000e+00 4.323480e-28 7.878366e-50
## 34 1.000000e+00 2.371974e-32 4.218843e-55
## 44 1.000000e+00 7.284354e-18 2.344684e-35
## 51 3.482024e-20 9.999108e-01 8.921915e-05
## 60 2.108587e-22 9.994953e-01 5.047370e-04
## 70 1.717244e-19 9.999965e-01 3.547080e-06
## 89 2.753852e-19 9.999581e-01 4.187694e-05
## 92 7.431814e-24 9.981819e-01 1.818098e-03
## 101 1.698805e-55 4.638431e-09 1.000000e+00
## 119 4.077561e-65 4.227345e-10 1.000000e+00
## 125 1.930807e-42 7.142563e-05 9.999286e-01
## 131 2.082064e-46 7.590232e-05 9.999241e-01
## 143 1.949097e-41 6.623572e-04 9.993376e-01
##
## $x
## LD1 LD2
## 2 7.395395 -0.76528994
## 16 9.794555 2.89496048
## 23 9.181468 1.06310065
## 34 10.079800 1.99514454
## 44 6.781428 1.34699455
## 51 -1.572111 -0.06640228
## 60 -2.070049 -0.23762022
## 70 -1.260504 -1.62569284
## 89 -1.367295 -0.02978819
## 92 -2.404729 -0.26100090
## 101 -8.061370 2.31431193
## 119 -9.737127 -0.93873481
## 125 -5.912590 1.36983307
## 131 -6.624276 -0.85112223
## 143 -5.776469 0.05107727Veamos la matriz de confusión en df_train
table(predict(model_lda,type="class")$class,df_train$Species)##
## setosa versicolor virginica
## setosa 45 0 0
## versicolor 0 43 1
## virginica 0 2 44partimat muestra una matriz de gráficos para cada
combinación de dos variables. Cada gráfico muestra una vista diferente
de los mismos datos. Las regiones coloreadas delimitan cada área de
clasificación. Se predice que cualquier observación que se encuentre
dentro de una región pertenece a una clase específica. Cada gráfico
también incluye la tasa de error aparente para esa vista de los
datos.
library(klaR)
partimat (Species~. , data=df_train , method="lda")
Veamos la performance en df_test
lda.test <- predict(model_lda,df_test)
df_test$lda <- lda.test$class
table(df_test$lda,df_test$Species)##
## setosa versicolor virginica
## setosa 5 0 0
## versicolor 0 5 0
## virginica 0 0 5Biplot
#install.packages("devtools")
#library(devtools)
#install_github("fawda123/ggord")
library(ggord)
ggord(model_lda, df_train$Species, xlim = c(-10, 11))
Análisis discriminante cuadrático (QDA)
El discriminante se dice cuadrático porque el término de segundo orden no se cancela como en el caso del discriminante lineal. QDA no asume la igualdad en la matriz de varianzas/covarianzas. En otras palabras, para QDA la matriz de covarianza puede ser diferente para cada clase.
Vamos a aplicar QDA, a pesar de que no se satisface el supuesto de normalidad multivariada.
model_qda <- qda(Species ~ ., df_train)
model_qda## Call:
## qda(Species ~ ., data = df_train)
##
## Prior probabilities of groups:
## setosa versicolor virginica
## 0.3333333 0.3333333 0.3333333
##
## Group means:
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 4.991111 3.393333 1.471111 0.2377778
## versicolor 5.940000 2.757778 4.262222 1.3266667
## virginica 6.566667 2.977778 5.506667 2.0133333partimat(Species ~ ., data=df_train, method="qda")
table(predict(model_qda,type="class")$class,df_train$Species)##
## setosa versicolor virginica
## setosa 45 0 0
## versicolor 0 43 1
## virginica 0 2 44lda.test_qda <- predict(model_qda,df_test)
df_test$qda <- lda.test_qda$class
table(df_test$qda,df_test$Species)##
## setosa versicolor virginica
## setosa 5 0 0
## versicolor 0 5 0
## virginica 0 0 5Algunos comentarios
El ánalisis discriminante es una técnica sensible a la presencia de outliers
Si el supuesto de normalidad multivariada se sostiene pero hay outliers es recomendable recurrir a la versión robusta de QDA.
Si se rechaza la normalidad multivariada, el modelo robusto no es adecuado.
¿Qué hacer si no se cumple el supuesto de normalidad multivariada?
Hay diversas opiniones al respecto. Algunos mencionan que LDA es robusto a la falta de normalidad multivariada si los datos son lo suficientemente grandes.1 Sin embargo, Lachenbruch et al. (1973)2 han demostrado que los resultados de aplicar LDA si se ven muy afectados por la falta de normalidad multivariada luego de usar tres tipos de distribuciones no normales. Por otra parte, también se ha encontrado que QDA si es robusto a la falta de normalidad si las distribuciones de los datos se alejan ligeramente de distribución teórica3
Actualmente, se han propuesto alternativas no paramétricas para los casos en que el supuesto de normalidad no se cumple.4 haciendo uso de matrices de dispersión, por ejemplo. Además, Sharipah Soaad Syed Yahaya et al. (2016)5 proponen el uso de dos estimadores robusto llamados one-step M-estimator (MOM) y winsorized modified one-step M-estimator (WMOM), que permiten trabajar con aquellos datos que no solo no siguen una distribución normal sino que tambien son heterocedásticos.
Finalmente, se propone también el uso de algoritmos mas complejos que permiten proseguir cuando no se cumplen los supuestos y no hay necesidad de transformar los datos para ello. Uno de los modelos propuestos es el uso de Árboles de decisión para clasificación6
Discriminant Function Analysis Introductory Overview - Assumptions↩︎
Robustness of the linear and quadratic discriminant function to certain types of non-normality↩︎
How non-normality affects the quadratic discriminant function↩︎
Classification Trees as an Alternative to Linear Discriminant Analysis↩︎