1 Definición de LDA

El análisis discriminante (LDA acrónimo de Linear Discriminant Analysis) se utiliza para predecir la probabilidad de pertenecer a una clase (o categoría) determinada en función de una o varias variables predictoras. Trabaja con variables predictoras continuas y/o categóricas.

Análisis discriminante lineal (LDA): utiliza combinaciones lineales de predictores para predecir la clase de una observación determinada. Supone que las variables predictoras (p) están distribuidas normalmente y que las clases tienen varianzas idénticas (para el análisis univariante, \(p = 1\)) o matrices de covarianza idénticas (para el análisis multivariante, \(p > 1\)).

Bibliografía Recomendada:

[Libro: Machine Learning Essentials. Practical Guide in R] (https://www.datanovia.com/en/product/machine-learning-essentials-practical-guide-in-r/?url=/5-bookadvisor/54-machine-learning-essentials/)

Carga de paquetes de R requeridos:
tidyverse para facilitar la manipulación y visualización de datos caret para facilitar el flujo de trabajo de aprendizaje automático

library(tidyverse)            # para facilitar la manipulación y visualización de datos
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.2     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.0
## ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
## ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
## ✔ purrr     1.0.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(caret)                # Para aplicar modelos
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
theme_set(theme_classic())    # Diseño

1.1 Ejemplo de la base de datos

Recogemos 150 registros de 3 especies de flores. Para cada registro se recogen 4 variables numéricas:
- Sepal.length (largo de Sépalo) - Sepal.Width (ancho de Sépalo) - Petal.Length (largo de Pétalo) - Petal.Width (ancho de Pétalo)

Las especies: - Setosa - Versicolor - Virginica

# Cargamos la base de datos
data("iris")
# Es necesario crar una parte para "entrenar" al modelo y otra para "testear"
set.seed(123)                                   # Para que la patición aleatoria sea fija y se siga el ejemplo            
entreno  <- iris$Species %>%
  createDataPartition(p = 0.8, list = FALSE)    # el 80% de los registros, se envían a "entreno"
train.data <- iris[entreno, ]                   # train.data es la parte de datos de entreno
test.data <- iris[-entreno, ]                   # test.data es la parte de datos para testear el modelo

Estudiemos la base de datos mediante histogramas por especies.

library(ggplot2)
library(ggpubr)

plot1 <- ggplot(data = iris, aes(x = Sepal.Length)) +
         geom_density(aes(colour = Species)) + theme_bw()
plot2 <- ggplot(data = iris, aes(x = Sepal.Width)) +
         geom_density(aes(colour = Species)) + theme_bw()
plot3 <- ggplot(data = iris, aes(x = Petal.Length)) +
         geom_density(aes(colour = Species)) + theme_bw()
plot4 <- ggplot(data = iris, aes(x = Petal.Width)) +
         geom_density(aes(colour = Species)) + theme_bw()
# la función grid.arrange del paquete grid.extra permite ordenar
# graficos de ggplot2
ggarrange(plot1, plot2, plot3, plot4, common.legend = TRUE, legend = "bottom")

Un análisis por pares

library(psych)
## Warning: package 'psych' was built under R version 4.3.1
## 
## Attaching package: 'psych'
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
pairs.panels(iris[1:4],
             gap = 0,
             bg = c("red", "green", "blue")[iris$Species],
             pch = 21)

Evaluamos la normalidad de las muestras por cada variable numérica y por especie.

4 x 3 = 12 análisis.

Para P-value < 0.05 la muestra no se distribuye bajo la normalidad estadística

#Contraste de normalidad Shapiro-Wilk para cada variable en cada especie
library(reshape2)
## 
## Attaching package: 'reshape2'
## The following object is masked from 'package:tidyr':
## 
##     smiths
library(knitr)
library(dplyr)
datos_tidy <- melt(iris, value.name = "valor")
## Using Species as id variables
kable(datos_tidy %>% group_by(Species, variable) %>% summarise(p_value_Shapiro.test = round(shapiro.test(valor)$p.value,5)))
## `summarise()` has grouped output by 'Species'. You can override using the
## `.groups` argument.
Species variable p_value_Shapiro.test
setosa Sepal.Length 0.45951
setosa Sepal.Width 0.27153
setosa Petal.Length 0.05481
setosa Petal.Width 0.00000
versicolor Sepal.Length 0.46474
versicolor Sepal.Width 0.33800
versicolor Petal.Length 0.15848
versicolor Petal.Width 0.02728
virginica Sepal.Length 0.25831
virginica Sepal.Width 0.18090
virginica Petal.Length 0.10978
virginica Petal.Width 0.08695

Hemos de tener en cuenta que la variable de “ancho de pétalo” de la especie Versicolor no se distribuye de forma normal, con una p-value 0.02 y Setosa “ancho de pétalo” con una p-value = 0

Escalamos los datos de manera “normalizada”.

preproc.param <- train.data %>% 
  preProcess(method = c("center", "scale"))
# Trabajamos el modelo con los datos normalizados.
train.transformed <- preproc.param %>% predict(train.data)
test.transformed <- preproc.param %>% predict(test.data)
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
# Creamos el modelo
set.seed(123)   
model <- lda(Species~., data = train.transformed)
# Generamos las predicciones
predictions <- model %>% predict(test.transformed)
# Evaluamos el modelo mediante la métrica de precisión
mean(predictions$class==test.transformed$Species)
## [1] 0.9666667
# Creamos una matriz de confusión
p1 <- predict(model, test.transformed)$class
tab1 <- table(Predicted = p1, Actual = test.transformed$Species)
tab1
##             Actual
## Predicted    setosa versicolor virginica
##   setosa         10          0         0
##   versicolor      0         10         1
##   virginica       0          0         9

96.6% de acierto del modelo.

model
## Call:
## lda(Species ~ ., data = train.transformed)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa       -1.0112835  0.78048647   -1.2900001  -1.2453195
## versicolor    0.1014181 -0.68674658    0.2566029   0.1472614
## virginica     0.9098654 -0.09373989    1.0333972   1.0980581
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Length  0.6794973  0.04463786
## Sepal.Width   0.6565085 -1.00330120
## Petal.Length -3.8365047  1.44176147
## Petal.Width  -2.2722313 -1.96516251
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9902 0.0098

El LDA determina las medias de los grupos y calcula, para cada individuo, la probabilidad de pertenecer a los diferentes grupos. A continuación, el individuo se ve afectado al grupo con la puntuación de probabilidad más alta.

Las salidas contienen los siguientes elementos:lda()

  • Probabilidades previas de los grupos: la proporción de observaciones de entrenamiento en cada grupo. Por ejemplo, hay un 31% de las observaciones de entrenamiento en el grupo setosa
  • Grupo significa: grupo centro de gravedad. Muestra la media de cada variable en cada grupo.
  • Coeficientes de discriminantes lineales: muestra la combinación lineal de variables predictoras que se utilizan para formar la regla de decisión LDA. por ejemplo. Semejantemente.LD1 = 0.6794973Sepal.Length + 0.6565085 Sepal.Width - 3.8365047Petal.Length - 2.2722313Petal.WidthLD2 = 0.04463786Sepal.Length - 1.00330120Sepal.Width + 1.44176147Petal.Length - 1.96516251Petal.Width

También tenemos la Proportion of trace , el porcentaje de separaciones archivadas por la primera función discriminante LD1 es 99,02% .

Dibujamos las predicciones y la situación de los registros en función de las 2 variables LDA1 y LDA2 para ver si es buena la discriminación.

lda.data <- cbind(train.transformed, predict(model)$x)
ggplot(lda.data, aes(LD1, LD2)) +
  geom_point(aes(color = Species))

Para visualizar las particiones, y evaluar qué variable categórica es la mejor.

library(klaR)
## Warning: package 'klaR' was built under R version 4.3.3
partimat(Species ~ Sepal.Width + Sepal.Length + Petal.Length + Petal.Width,
         data = iris, method = "lda", prec = 200,
         image.colors = c("darkgoldenrod1", "snow2", "skyblue2"),
         col.mean = "firebrick")

Ahora podemos crear un histograma apilado de valores de función discriminante.

p <- predict(model, train.transformed)
ldahist(data = p$x[,1], g = train.transformed$Species)

En el gráfico anterior, tenemos el histograma de LD1, y podemos ver que la separación entre setosa y las otras dos especies es bastante grande y no se superpone. Por el contrario, existe cierta superposición entre versicolor y virginica . Ya dijimos que el porcentaje de separación archivado por LD1 es 99,32% , es decir, podemos ver una separación muy clara en el histograma de arriba. Ahora podemos intentar hacer lo mismo con LD2

p <- predict(model, train.transformed)
ldahist(data = p$x[,2], g = train.transformed$Species)

La variable LDA2 no es buena discriminante, hay mucha superposición de valores.

Webgrafia

http://www.sthda.com/english/articles/36-classification-methods-essentials/146-discriminant-analysis-essentials-in-r/

https://cienciadedatos.net/documentos/28_linear_discriminant_analysis_lda_y_quadratic_discriminant_analysis_qda#Ejemplo_con_Iris_data

https://www.andreaperlato.com/mlpost/linear-discriminant-analysis/