Minería de datos

Análisis discriminante lineal (LDA)

El Análisis Discriminante Lineal o Linear Discrimiant Analysis (LDA) es un método de clasificación supervisado 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 posterior 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, 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.

El proceso de un análisis discriminante puede resumirse en 6 pasos:

Disponer de un conjunto de datos de entrenamiento (training data) en el que se conoce a que grupo pertenece cada observación.
Calcular las probabilidades a priori (prior probabilities): la proporción esperada de observaciones que pertenecen a cada grupo.
Determinar si la varianza o matriz de covarianzas es homogénea en todos los grupos. De esto dependerá que se emplee LDA o QDA.
Estimar los parámetros necesarios para las funciones de probabilidad condicional, verificando que se cumplen las condiciones para hacerlo.
Calcular el resultado de la función discriminante. El resultado de esta determina a qué grupo se asigna cada observación.
Utilizar validación cruzada (cross-validation) para estimar las probabilidades de clasificaciones erróneas.

Teorema de Bayes para clasificación

Considérense dos eventos $A$ y $B$, sabemos que el teorema de Bayes establece que: \[ P(B|A)=\dfrac{P(AB)}{P(A)}\] Supóngase que se desea clasificar una nueva observación en una de las $K$ clases de una variable cualitativa $Y$, siendo $K≥2$, a partir de un solo predictor $X$. Se dispone de las siguientes definiciones:

Se define como probabilidad a priori $(π_k)$ la probabilidad de que una observación aleatoria pertenezca a la clase $k$.
Se define $f_k(X)≡P(X=x|Y=k)$ como la función de densidad de probabilidad condicional de $X$ para una observación que pertenece a la clase $k$. Cuanto mayor sea $f_k(X)$ mayor la probabilidad de que una observación de la clase $k$ adquiera una valor de $X≈x$.
Se define como probabilidad posterior $P(Y=k|X=x)$ la probabilidad de que una observación pertenezca a la clase $k$ siendo $x$ el valor del predictor.

Aplicando del teorema de Bayes se pueden conocer la probabilidad posterior para cada clase:

\[P(Y=k|X=x)=\dfrac{π_k \, f_k(x)}{\displaystyle\sum_{j=1}^k π_jf_j(x)}\] Para que la clasificación basada en Bayes sea posible, se necesita conocer la probabilidad poblacional de que una observación cualquiera pertenezca a cada clase ($π_k$) y la probabilidad poblacional de que una observación que pertenece a la clase $k$ adquiera el valor $x$ en el predictor ($f_k(X)≡P(X=x|Y=k)$). En la práctica, raramente se dispone de esta información, por lo que los parámetros tienen que ser estimados a partir de la muestra. Como consecuencia, el clasificador LDA obtenido se aproxima al clasificador de Bayes.

Estimación de $π_k$ y $f_k(X)$

La capacidad del LDA para clasificar correctamente las observaciones depende de cómo de buenas sean las estimaciones de $π_k$ y $f_k(X)$. Cuanto más cercanas al valor real, más se aproximará el clasificador LDA al clasificador de Bayes. En el caso de la prior probability ($π_k$) la estimación suele ser sencilla, la probabilidad de que una observación cualquiera pertenezca a la clase $k$ es igual al número de observaciones de esa clase entre el número total de observaciones $\hat{π}_k=\frac{n_k}{N}$.

La estimación de $f_k(X)$ no es tan directa y para conseguirla se requiere de ciertas asunciones. Si se considera que $f_k(X)$ se distribuye de forma normal en las $K$ clases, entonces se puede estimar su valor a partir de la densidad normal.

Si además se asume que la varianza es constante en todos los grupos $σ^2_1=σ_2^2=\dots=σ^2_k=σ^2$, entonces, el sumatorio $∑^K_{j=1} π_jf_j(x)$ se simplifica en gran medida permitiendo calcular la probabilidad posterior según la ecuación:

\[P(Y=k|X=x)=\dfrac{π_k\frac{1}{2π}\exp\left(−\frac{1}{2σ^2}(x−μ_k)^2\right)}{\displaystyle\sum_{j=1}^k π_j\frac{1}{2πσ}\exp \left(−\frac{1}{2σ^2}(x−μ_j)^2\right)}\] Esta ecuación se simplifica aun más mediante una transformación logarítmica de sus dos términos:

\[\log(P(Y=k|X=x))=x \frac{μ_k}{σ^2}−\frac{μ^2_k}{σ^2}+\log(π_k)\] El término lineal en el nombre Análisis discriminante lineal se debe al hecho de que la función discriminatoria es lineal respecto de $X$.

En la práctica, a pesar de tener una certeza considerable de que $X$ se distribuye de forma normal dentro de cada clase, los valores $μ_1,\dots,μ_k$, $π_1,\dots,π_k$ y $σ^2$ se desconocen, por lo que tienen que ser estimados a partir de las observaciones. Las estimaciones empleadas en LDA son:

\[\begin{align*} \hat{μ}_k = & \frac{1}{n_k} \displaystyle\sum_{i:y_1} x_i \\ \hat{σ}_k = & \frac{1}{N−K} \displaystyle\sum_{k=1}^K \displaystyle\sum_{i:y_1} (x_i−\hat{μ}_k)^2 \\ \hat{\pi}_k = & \frac{n_k}{N} \end{align*}\]

$\hat{μ}_k$ es la media de las observaciones del grupo $k$, $\hat{\sigma}_k$ es la media ponderada de las varianzas muestrales de las $K$ clases y $\hat{\pi}_k$ la proporción de observaciones de la clase $k$ respecto al tamaño total de la muestra.

La clasificación de Bayes consiste en asignar cada observación $X = x$ a aquella clase para la que $P(Y=k|X=x)$ sea mayor. En el caso particular de una variable cualitativa Y con solo dos niveles, se puede expresar la regla de clasificación como un cociente entre las dos probabilidades posteriores. Se asignará la observación a la clase 1 si $\frac{P(Y=1|X=x)}{P(Y=2|X=x)}>1$, y a la clase 2 si es menor. En este caso particular el límite de decisión de Bayes viene dado por $x=\frac{μ1+μ2}{2}$.

La siguiente imagen muestra dos grupos distribuidos de forma normal con medias $μ_1= -1.25$, $μ_2= 1.25$ y varianzas $σ^2_1=σ^2_2= 1$. Dado que se conoce el valor real de las medias y varianzas poblacionales (esto en la realidad no suele ocurrir), se puede calcular el límite de decisión de Bayes $x== 0 (linea discontinua).

library(ggplot2)

## Warning in as.POSIXlt.POSIXct(Sys.time()): unable to identify current timezone 'H':
## please set environment variable 'TZ'

ggplot(data.frame(x = c(-4, 4)), aes(x)) +
stat_function(fun = dnorm, args = list(mean = -1.25, sd = 1),
              color = "firebrick") + 
stat_function(fun = dnorm, args = list(mean = 1.25, sd = 1), color = "green3") +
geom_vline(xintercept = 0, linetype = "longdash") +
theme_bw()

Si en lugar de conocer la verdadera distribución poblacional de cada grupo solo se dispone de muestras, escenario que suele ocurrir en los casos reales, el límite de decisión LDA se aproxima al verdadero límite de decisión de Bayes pero no es exacto. Cuanto más representativas sean las muestras mejor la aproximación.

set.seed(1981)
library(ggplot2)
n <-30
grupo_a <- rnorm(n, mean = -1.25, sd = 1)
grupo_b <- rnorm(n, mean = 1.25, sd = 1)
datos <- data.frame(valor = c(grupo_a, grupo_b),
                    grupo = rep(c("A","B"), each = n))

ggplot(data = datos, aes(x = valor, fill = grupo)) +
geom_histogram(alpha = 0.5, position = "identity") +
geom_vline(xintercept = 0, linetype = "longdash") +
geom_vline(xintercept = (mean(grupo_a) + mean(grupo_b))/2)  +
annotate(geom = "text", x = 1.5, y = 9, label = "Límite decisión LDA") +
annotate(geom = "text", x = -1.5, y = 10, label = "Límite decisión Bayes") +
theme_bw() + 
theme(legend.position = "top")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Extensión del LDA para múltiples predictores

Los conceptos anteriormente descritos empleando un único predictor pueden generalizarse para introducir múltiples predictores en el modelo. La diferencia reside en que $X$, en lugar de ser un único valor, es un vector formado por el valor de $p$ predictores $X=(X_1,X_2,\dots,X_n)$ y que, en lugar de proceder de una distribución normal, procede de una distribución normal multivariante.

Un vector sigue una distribución $k$-normal multivariante si cada uno de los elementos individuales que lo forman sigue una distribución normal y lo mismo para toda combinación lineal de sus $k$ elementos. Las siguientes imágenes muestran representaciones gráficas de distribuciones normales multivariante de 2 elementos (distribución normal bivariante).

mu1 <- 0 # set mean x1
mu2 <- 0 # set mean x2
s11 <- 10 # set variance x1
s22 <- 10 # set variance x2
s12 <- 15 # set covariance x1 and x2
rho <- 0.5 # set correlation coefficient  x1 and x2
x1 <- seq(-10,10,length = 41) # generate vector  x1
x2 <- x1 # copy x1 to x2

f <- function(x1,x2) # multivariate function
{
  term1 <- 1/(2 * pi * sqrt(s11*s22*(1 - rho^2)))
  term2 <- -1/(2 * (1 - rho^2))
  term3 <- (x1 - mu1)^2/s11
  term4 <- (x2 - mu2)^2/s22
  term5 <- -2*rho*((x1 - mu1)*(x2 - mu2))/(sqrt(s11)*sqrt(s22))
  term1*exp(term2*(term3 + term4 - term5))
} 

z <- outer(x1,x2,f) # calculate density values

persp(x1, x2, z, # 3-D plot
  main = "Distribución multivariante con dos predictores",
  col = "lightgreen",
  theta = 30, phi = 20,
  r = 50,
  d = 0.1,
  expand = 0.5,
  ltheta = 90, lphi = 180,
  shade = 0.75,
  ticktype = "simple",
  nticks = 5)

# Otra forma de representar una distribución bivariante
library(mvtnorm)
library(scatterplot3d)

ss <- matrix(c(1,0,0,1), ncol = 2) 
x1000 <- rmvnorm(n = 1000, mean = c(0,0), sigma = ss)
scatterplot3d(x1000[,1], x1000[,2],
              dmvnorm(x1000, mean = c(0,0), sigma = ss),
              highlight = TRUE, xlab = "x", ylab = "y", zlab = "z")

Condiciones de LDA

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

Cada predictor que forma parte del modelo se distribuye de forma normal en cada una de las clases de la variable respuesta. En el caso de múltiples predictores, las observaciones siguen una distribución normal multivariante en todas las clases.
La varianza del predictor es igual en todas las clases de la variable respuesta. En el caso de múltiples predictores, la matriz de covarianza es igual en todas las clases. Si esto no se cumple se recurre a Análisis Discriminante Cuadrático (QDA).
Cuando la condición de normalidad no se cumple, LDA pierde precisión pero aun así puede llegar a clasificaciones relativamente buenas.

Ejemplo datos insectos

Un equipo de biólogos quiere generar un modelo estadístico que permita identificar a que especie $(a \text{ o } b)$ pertenece un determinado insecto. Para ello se han medido tres variables (longitud de las patas, diámetro del abdomen y diámetro del órgano sexual) en 10 individuos de cada una de las dos especies.

input <- ("
especie pata abdomen organo_sexual 
a 191 131 53
a 185 134 50
a 200 137 52
a 173 127 50
a 171 128 49
a 160 118 47
a 188 134 54
a 186 129 51
a 174 131 52
a 163 115 47
b 186 107 49
b 211 122 49
b 201 144 47
b 242 131 54
b 184 108 43
b 211 118 51
b 217 122 49
b 223 127 51
b 208 125 50
b 199 124 46
")
datos <- read.table(textConnection(input), header = TRUE)

Exploración gráfica de los datos

library(ggplot2)
library(ggpubr)

## Loading required package: magrittr

p1 <- ggplot(data = datos, aes(x = pata, fill = especie)) +
      geom_histogram(position = "identity", alpha = 0.5)
p2 <- ggplot(data = datos, aes(x = abdomen, fill = especie)) +
      geom_histogram(position = "identity", alpha = 0.5)
p3 <- ggplot(data = datos, aes(x = organo_sexual, fill = especie)) +
      geom_histogram(position = "identity", alpha = 0.5)
ggarrange(p1, p2, p3, nrow = 3, common.legend = TRUE)

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

A nivel individual, la longitud de la pata parece ser la variable que más se diferencia entre especies (menor solapamiento entre poblaciones).

pairs(x = datos[, c("pata","abdomen","organo_sexual")],
      col = c("blue", "green4")[datos$especie], pch = 19)

El par de variables abdomen-pata y el par pata-organo_sexual parecen separar bien las dos especies.

library(scatterplot3d)
scatterplot3d(datos$pata, datos$abdomen, datos$organo_sexual,
              color = c("firebrick", "green3")[datos$especie], pch = 19,
              grid = TRUE, xlab = "pata", ylab = "abdomen",
              zlab = "organo sexual", angle = 65, cex.axis = 0.6)
legend("topleft",
       bty = "n", cex = 0.8,
       title = "Especie",
       c("a", "b"), fill = c("firebrick", "green3"))

La representación de las tres variables de forma simultánea parece indicar que las dos especies sí están bastante separadas en el espacio 3D generado.

Probabilidades a priori

Como no se dispone de información sobre la abundancia relativa de las especies a nivel poblacional, se considera como probabilidad a pirori de cada especie el número de observaciones de la especie entre el número de observaciones totales.

\[\hat{π}_a=\hat{π}_b=\dfrac{10}{20}=0.5\]

Homogeneidad de Varianza

De entre los diferentes test que contrastan la homogeneidad de varianza, el más recomendable cuando solo hay un predictor, dado que se asume que se distribuye de forma normal, es el test de Barttlet. Cuando se emplean múltiples predictores, se tiene que contrastar que la matriz de covarianzas $∑$ es constante en todos los grupos, siendo recomendable comprobar también la homogeneidad de varianza para cada predictor a nivel individual.

El test Box fue desarrollado por el matemático Box (1949) como una extensión del test de Barttlet para escenarios multivariante y permite contrastar la igualdad de matrices entre grupos. El test Box es muy sensible a violaciones de la normalidad multivariante, por lo que esta debe ser contrastada con anterioridad.

# Representación mediante Histograma de cada variable para cada especie 
par(mfcol = c(2, 3))
for (k in 2:4) {
  j0 <- names(datos)[k]
  #br0 <- seq(min(datos[, k]), max(datos[, k]), le = 11)
  x0 <- seq(min(datos[, k]), max(datos[, k]), le = 50)
  for (i in 1:2) {
    i0 <- levels(datos$especie)[i]
    x <- datos[datos$especie == i0, j0]
    hist(x, proba = T, col = grey(0.8), main = paste("especie", i0),
    xlab = j0)
    lines(x0, dnorm(x0, mean(x), sd(x)), col = "red", lwd = 2)
  }
}

# Representación de cuantiles normales de cada variable para cada especie 
par(mfcol = c(2, 3))
for (k in 2:4) {
  j0 <- names(datos)[k]
  x0 <- seq(min(datos[, k]), max(datos[, k]), le = 50)
  for (i in 1:2) {
    i0 <- levels(datos$especie)[i]
    x <- datos[datos$especie == i0, j0]
    qqnorm(x, main = paste("especie", i0, j0), pch = 19, col = i + 1)
    qqline(x)
  }
}

# Contraste de normalidad Shapiro-Wilk para cada variable en cada especie
library(reshape2)
library(knitr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

datos_tidy <- melt(datos, value.name = "valor")

## Using especie as id variables

kable(datos_tidy %>% group_by(especie, variable) %>% summarise(p_value_Shapiro.test = shapiro.test(valor)$p.value))

especie	variable	p_value_Shapiro.test
a	pata	0.7763034
a	abdomen	0.1845349
a	organo_sexual	0.6430844
b	pata	0.7985711
b	abdomen	0.5538213
b	organo_sexual	0.8217855

# Misma operación con aggregate
aggregate(formula = valor ~ especie + variable, data = datos_tidy,
          FUN = function(x){shapiro.test(x)$p.value})

##   especie      variable     valor
## 1       a          pata 0.7763034
## 2       b          pata 0.7985711
## 3       a       abdomen 0.1845349
## 4       b       abdomen 0.5538213
## 5       a organo_sexual 0.6430844
## 6       b organo_sexual 0.8217855

No hay evidencias de falta de normalidad univariante en ninguna de las variables empleadas como predictores en ninguno de los grupos.

El paquete MVN contiene funciones que permiten realizar los tres test de hipótesis comúnmente empleados para evaluar la normalidad multivariante (Mardia, Henze-Zirkler y Royston) y también funciones para identificar outliers que puedan influenciar en el contraste. Para información detallada de cada uno consultar https://cran.r-project.org/web/packages/MVN/vignettes/MVN.pdf.

library(MVN)

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

## sROC 0.1-2 loaded

outliers <- mvn(data = datos[,-1], mvnTest = "hz", multivariateOutlierMethod = "quan")

royston_test <- mvn(data = datos[,-1], mvnTest = "royston", multivariatePlot = "qq")

royston_test$multivariateNormality

##      Test         H   p value MVN
## 1 Royston 0.4636176 0.9299447 YES

hz_test <- mvn(data = datos[,-1], mvnTest = "hz")
hz_test$multivariateNormality

##            Test        HZ    p value MVN
## 1 Henze-Zirkler 0.7870498 0.07666139 YES

A pesar de los 5 outliers detectados, ninguno de los dos test encuentran evidencias significativas $(α=0.05)$ de falta de normalidad multivariante.

Finalmente, mediante la función boxM() del paquete biotools se realiza el contraste de matrices de covarianza.

library(biotools)

## Loading required package: rpanel

## Loading required package: tcltk

## Package `rpanel', version 1.1-4: type help(rpanel) for summary information

## Loading required package: tkrplot

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## Loading required package: lattice

## Loading required package: SpatialEpi

## Loading required package: sp

## ---
## biotools version 3.1

##

boxM(data = datos[, 2:4], grouping = datos[, 1])

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  datos[, 2:4]
## Chi-Sq (approx.) = 9.831, df = 6, p-value = 0.132

Se puede aceptar que la matriz de covarianza es igual en todos los grupos.

Estimación de los parámetros de la función de densidad $(\hat{μ}(X),∑)$ y cálculo de la función discriminante.

Estos dos pasos se realizan mediante la función lda del paquete MASS. lda realiza la clasificación mediante la aproximación de Fisher.

modelo_lda <- lda(formula = especie ~ pata + abdomen + organo_sexual,
                  data = datos)

Una vez obtenidas las funciones discriminantes, se puede clasificar un nuevo insecto en función de sus medidas. Por ejemplo, un nuevo espécimen cuyas medidas sean: pata = 194, abdomen = 124, organo_sexual = 49.

nuevas_observaciones <- data.frame(pata = 194, abdomen = 124,
                                   organo_sexual = 49)
predict(object = modelo_lda, newdata = nuevas_observaciones)

## $class
## [1] b
## Levels: a b
## 
## $posterior
##            a         b
## 1 0.05823333 0.9417667
## 
## $x
##         LD1
## 1 0.5419421

El resultado muestra que, según la función discriminante, la probabilidad posterior de que el espécimen pertenezca a la especie $b$ es del $94.2\%$ frente al $5.8\%$ de que pertenezca a la especie $a$.

Evaluación de los errores de clasificación

predicciones <- predict(object = modelo_lda, newdata = datos[, -1],
                        method = "predictive")
table(datos$especie, predicciones$class,
      dnn = c("Clase real", "Clase predicha"))

##           Clase predicha
## Clase real  a  b
##          a 10  0
##          b  0 10

trainig_error <- mean(datos$especie != predicciones$class) * 100
paste("trainig_error=", trainig_error, "%")

## [1] "trainig_error= 0 %"

Empleando las mismas observaciones con las que se ha generado el modelo discriminante (trainig data), la precisión de clasificación es del $100\%$. Evaluar un modelo con los mismos datos con los que se ha creado suele resultar en estimaciones de la precisión demasiado optimistas (error muy bajo). Por lo tanto, la estimación del test error mediante validación cruzada es más adecuada para obtener una evaluación realista del modelo.

La siguiente imagen muestra la representación de las observaciones, coloreadas por la verdadera especie a la que pertenecen y acompañadas por una etiqueta con la especie que ha predicho el LDA.

with(datos, {
  s3d <- scatterplot3d(pata, abdomen, organo_sexual,
                       color = c("firebrick", "green3")[datos$especie],
                       pch = 19, grid = TRUE, xlab = "pata", ylab = "abdomen",
                       zlab = "organo sexual", angle = 65, cex.axis = 0.6)
  
  s3d.coords <- s3d$xyz.convert(pata, abdomen, organo_sexual)
  # convierte coordenadas 3D en proyecciones 2D
    
  text(s3d.coords$x, s3d.coords$y, # cordenadas x, y
       labels = datos$especie,     # texto
       cex = .8, pos = 4)   
    
  legend("topleft", 
         bty = "n", cex = 0.8,
         title = "Especie",
         c("a", "b"), fill = c("firebrick", "green3"))
})

Ejemplo con Iris data

data("iris")
kable(head(iris, n = 10))

Sepal.Length	Sepal.Width	Petal.Length	Petal.Width	Species
5.1	3.5	1.4	0.2	setosa
4.9	3.0	1.4	0.2	setosa
4.7	3.2	1.3	0.2	setosa
4.6	3.1	1.5	0.2	setosa
5.0	3.6	1.4	0.2	setosa
5.4	3.9	1.7	0.4	setosa
4.6	3.4	1.4	0.3	setosa
5.0	3.4	1.5	0.2	setosa
4.4	2.9	1.4	0.2	setosa
4.9	3.1	1.5	0.1	setosa

probabilidades a priori

Como no se dispone de información sobre la abundancia relativa de las especies a nivel poblacional, se considera como probabilidad previa de cada especie el número de observaciones de la especie entre el número de observaciones totales.

\[\hat{π}_{\text{setosa}}=\hat{π}_{\text{versicolor}}=\hat{π}_{\text{virginica}}=\dfrac{50}{150}=\dfrac{1}{3}\]

Normalidad univariante, normalidad multivariante y homogeneidad de varianza

#representación mediante histograma de cada variable para cada especie 
par(mfcol = c(3, 4))
for (k in 1:4) {
  j0 <- names(iris)[k]
  x0 <- seq(min(iris[, k]), max(iris[, k]), le = 50)
  for (i in 1:3) {
    i0 <- levels(iris$Species)[i]
    x <- iris[iris$Species == i0, j0]
    hist(x, proba = T, col = grey(0.8), main = paste("especie", i0),
    xlab = j0)
    lines(x0, dnorm(x0, mean(x), sd(x)), col = "red", lwd = 2)
  }
}

#representación de cuantiles normales de cada variable para cada especie 
for (k in 1:4) {
  j0 <- names(iris)[k]
  x0 <- seq(min(iris[, k]), max(iris[, k]), le = 50)
  for (i in 1:3) {
    i0 <- levels(iris$Species)[i]
    x <- iris[iris$Species == i0, j0]
    qqnorm(x, main = paste(i0, j0), pch = 19, col = i + 1) 
    qqline(x)
  }
}

#Contraste de normalidad Shapiro-Wilk para cada variable en cada especie
library(reshape2)
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)))

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

Observamos que la variable petal.width no parece tener una distribución normal en los grupos setosa y versicolor.

Normalidad multivariante

library(MVN)
outliers <- mvn(data = iris[,-5], mvnTest = "hz", multivariateOutlierMethod = "quan")

royston_test <- mvn(data = iris[,-5], mvnTest = "royston", multivariatePlot = "qq")

royston_test$multivariateNormality

##      Test        H      p value MVN
## 1 Royston 50.39667 3.098229e-11  NO

hz_test <- mvn(data = iris[,-5], mvnTest = "hz")
hz_test$multivariateNormality

##            Test       HZ p value MVN
## 1 Henze-Zirkler 2.336394       0  NO

Ambos test muestran evidencias significativas de falta de normalidad multivariante. LDA tiene cierta robustez frente a la falta de normalidad multivariante, pero es importante tenerlo en cuenta en la conclusión del análisis.

library(biotools)
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

El test Box muestra evidencias de que la matriz de covarianza no es constante en todos los grupos, lo que a priori descartaría el método LDA en favor del QDA. Sin embargo, como el test Box’s M es muy sensible a la falta de normalidad multivariante, con frecuencia resulta significativo no porque la matriz de covarianza no sea constante sino por la falta de normalidad, cosa que ocurre para los datos de Iris. Por esta razón se va a asumir que la matriz de covarianza sí es constante y que LDA puede alcanzar una buena precisión en la clasificación. En la evaluación del modelo se verá como de buena es esta aproximación. Además, en las conclusiones se debe explicar la asunción hecha.

Función discriminante

library(MASS)
modelo_lda <- lda(Species ~ Sepal.Width + Sepal.Length + Petal.Length +
                  Petal.Width, data = iris)
modelo_lda

## Call:
## lda(Species ~ Sepal.Width + Sepal.Length + Petal.Length + Petal.Width, 
##     data = iris)
## 
## Prior probabilities of groups:
##     setosa versicolor  virginica 
##  0.3333333  0.3333333  0.3333333 
## 
## Group means:
##            Sepal.Width Sepal.Length Petal.Length Petal.Width
## setosa           3.428        5.006        1.462       0.246
## versicolor       2.770        5.936        4.260       1.326
## virginica        2.974        6.588        5.552       2.026
## 
## Coefficients of linear discriminants:
##                     LD1         LD2
## Sepal.Width   1.5344731  2.16452123
## Sepal.Length  0.8293776  0.02410215
## Petal.Length -2.2012117 -0.93192121
## Petal.Width  -2.8104603  2.83918785
## 
## Proportion of trace:
##    LD1    LD2 
## 0.9912 0.0088

Evaluación de los errores de clasificación

predicciones <- predict(object = modelo_lda, newdata = iris[, -5])
table(iris$Species, predicciones$class, dnn = c("Clase real", "Clase predicha"))

##             Clase predicha
## Clase real   setosa versicolor virginica
##   setosa         50          0         0
##   versicolor      0         48         2
##   virginica       0          1        49

trainig_error <- mean(iris$Species != predicciones$class) * 100
paste("trainig_error =", trainig_error, "%")

## [1] "trainig_error = 2 %"

Solo 3 de las 150 predicciones que ha realizado el modelo han sido erróneas. El error es muy bajo ($2\%$), lo que apunta a que el modelo es bueno. Sin embargo, para validarlo es necesario un nuevo set de datos con el que calcular el test error o recurrir a validación cruzada.

Visualización de las clasificaciones

La función partimat() del paquete klar permite representar los límites de clasificación de un modelo discriminante lineal o cuadrático para cada par de predictores. Cada color representa una región de clasificación acorde al modelo, se muestra el centroide de cada región y el valor real de las observaciones.

library(klaR)
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")

Análisis Discriminante Cuadrático

El clasificador cuadrático o Quadratic Discriminat Analysis (QDA) se asemeja en gran medida al LDA, con la única diferencia de que el QDA considera que cada clase $k$ tiene su propia matriz de covarianza ($∑_k$) y, como consecuencia, la función discriminante toma forma cuadrática:

\[\log(P(Y=k|X=x))=−\dfrac{1}{2}\log|Σ_k|−\dfrac{1}{2}(x−μ_k)^T Σ^{−1}_k(x−μ_k)+\log(π_k)\]

Para poder calcular la probabilidad posteriro a partir de esta ecuación discriminante es necesario estimar, para cada clase, $∑_k$, $μ_k$ y $π_k$ a partir de la muestra.

QDA genera límites de decisión curvos por lo que puede aplicarse a situaciones en las que la separación entre grupos no es lineal.

Ejemplo QDA con 2 predictores

Se dispone de los siguientes datos simulados.

set.seed(8558)
grupoA_x <- seq(from = -3, to = 4, length.out = 100)
grupoA_y <- 6 + 0.15 * grupoA_x - 0.3 * grupoA_x^2 + rnorm(100, sd = 1)
grupoA <- data.frame(variable_z = grupoA_x, variable_w = grupoA_y, grupo = "A")

grupoB_x <- rnorm(n = 100, mean = 0.5, sd = 0.8)
grupoB_y <- rnorm(n = 100, mean = 2, sd = 0.8)
grupoB <- data.frame(variable_z = grupoB_x, variable_w = grupoB_y, grupo = "B")

datos <- rbind(grupoA, grupoB)
plot(datos[, 1:2], col = datos$grupo, pch = 19)

La separación entre los grupos no es de tipo lineal, sino que muestra cierta curvatura. En este tipo de escenarios el método QDA es más adecuado que el LDA.

library(ggplot2)
library(ggpubr)
p1 <- ggplot(data = datos, aes(x = variable_z, fill = grupo)) +
      geom_histogram(position = "identity", alpha = 0.5)
p2 <- ggplot(data = datos, aes(x = variable_w, fill = grupo)) +
      geom_histogram(position = "identity", alpha = 0.5)
ggarrange(p1, p2, nrow = 2, common.legend = TRUE, legend = "bottom")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

La variable $W$ permite discriminar entre grupos mejor que la variable $z$.

Normalidad univariante, normalidad multivariante y homogeneidad de varianza

Distribución de los predictores de forma individual:

# Representación mediante histograma de cada variable para cada grupo 
par(mfcol = c(2, 2))
for (k in 1:2) {
  j0 <- names(datos)[k]
  x0 <- seq(min(datos[, k]), max(datos[, k]), le = 50)
  for (i in 1:2) {
    i0 <- levels(datos$grupo)[i]
    x <- datos[datos$grupo == i0, j0]
    hist(x, proba = T, col = grey(0.8), main = paste("grupo", i0),
    xlab = j0)
    lines(x0, dnorm(x0, mean(x), sd(x)), col = "red", lwd = 2)
  }
}

#representación de cuantiles normales de cada variable para cada grupo 
for (k in 1:2) {
  j0 <- names(datos)[k]
  x0 <- seq(min(datos[, k]), max(datos[, k]), le = 50)
  for (i in 1:2) {
    i0 <- levels(datos$grupo)[i]
    x <- datos[datos$grupo == i0, j0]
    qqnorm(x, main = paste(i0, j0), pch = 19, col = i + 1)
    qqline(x)
  }
}

#Contraste de normalidad Shapiro-Wilk para cada variable en cada grupo
library(reshape2)
datos_tidy <- melt(datos, value.name = "valor")

## Using grupo as id variables

library(dplyr)
datos_tidy %>%
  group_by(grupo, variable) %>% 
  summarise(p_value_Shapiro.test = round(shapiro.test(valor)$p.value,5))

## # A tibble: 4 x 3
## # Groups:   grupo [2]
##   grupo variable   p_value_Shapiro.test
##   <fct> <fct>                     <dbl>
## 1 A     variable_z              0.00172
## 2 A     variable_w              0.0932 
## 3 B     variable_z              0.625  
## 4 B     variable_w              0.810

La variable $Z$ no se distribuye de forma normal en el grupo $A$.

library(MVN)
outliers <- mvn(data = datos[,-3], mvnTest = "hz", multivariateOutlierMethod = "quan")

royston_test <- mvn(data = datos[,-3], mvnTest = "royston", multivariatePlot = "qq")

royston_test$multivariateNormality

##      Test        H     p value MVN
## 1 Royston 29.11024 4.77274e-07  NO

hz_test <- mvn(data = datos[,-3], mvnTest = "hz")
hz_test$multivariateNormality

##            Test       HZ      p value MVN
## 1 Henze-Zirkler 6.739874 5.317968e-14  NO

Ambos test muestran evidencias significativas de falta de normalidad multivariante. QDA tiene cierta robustez frente a la falta de normalidad multivariante, pero es importante tenerlo en cuenta en la conclusión del análisis.

library(MASS)
modelo_qda <- qda(grupo ~ variable_z + variable_w, data = datos)
modelo_qda

## Call:
## qda(grupo ~ variable_z + variable_w, data = datos)
## 
## Prior probabilities of groups:
##   A   B 
## 0.5 0.5 
## 
## Group means:
##   variable_z variable_w
## A  0.5000000   4.615307
## B  0.4864889   1.992911

predicciones <- predict(object = modelo_qda, newdata = datos)
table(datos$grupo, predicciones$class,
      dnn = c("Clase real", "Clase predicha"))

##           Clase predicha
## Clase real  A  B
##          A 97  3
##          B  7 93

trainig_error <- mean(datos$grupo != predicciones$class) * 100
paste("trainig_error=",trainig_error,"%")

## [1] "trainig_error= 5 %"

library(klaR)
partimat(formula = grupo ~ variable_z + variable_w, data = datos,
         method = "qda", prec = 400,
         image.colors = c("darkgoldenrod1", "skyblue2"),
         col.mean = "firebrick")

Minería de datos

Análisis discriminante lineal y Análisis discriminante cuadrático

Jairo Ayala

27/01/2020

Análisis discriminante lineal (LDA)

Teorema de Bayes para clasificación

Estimación de \(π_k\) y \(f_k(X)\)

Extensión del LDA para múltiples predictores

Condiciones de LDA

Ejemplo datos insectos

Exploración gráfica de los datos

Probabilidades a priori

Homogeneidad de Varianza

Estimación de los parámetros de la función de densidad \((\hat{μ}(X),∑)\) y cálculo de la función discriminante.

Evaluación de los errores de clasificación

Ejemplo con Iris data

probabilidades a priori

Normalidad univariante, normalidad multivariante y homogeneidad de varianza

Normalidad multivariante

Función discriminante

Evaluación de los errores de clasificación

Visualización de las clasificaciones

Análisis Discriminante Cuadrático

Ejemplo QDA con 2 predictores

Normalidad univariante, normalidad multivariante y homogeneidad de varianza