Con el fin de mostrar una aplicaciĂ³n en R de la metĂ³dologia de anĂ¡lisis discriminante lineal, flexible y penalizado, comparando cada metodo, se plantean dos ejemplos clasicos dentro la teoria, el primero es sobre la base de datos Iris con 3 clases para la variable de Respuesta; el segundo, seria sobre la base de datos Vowel (vocal).

Ejemplo con IRIS

La base de datos IRIS es un dataset que trae por defecto R, la cual se compone de 150 observaciones de flores de la planta Iris, a las flores se les midieron 5 variables su especie, el ancho del petalo, el largo del petalo, el ancho del sefalo y el largo del sepalo. Se presenta entonces las siguientes descriptivas considerando que la variables especies es categorica y cuentas con 3 clases Virginica, Setosa y Versicolor.

set.seed(123)
# Carga de los datos
data("iris")
# Descriptivas basicas de cada una de las variables
summary(iris)
##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500  
##        Species  
##  setosa    :50  
##  versicolor:50  
##  virginica :50  
##                 
##                 
## 
ds<-sqrt(diag(var(iris[,1:4])))
ds
## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##    0.8280661    0.4358663    1.7652982    0.7622377

De estas descriptivas se resalta el hecho de que las variables continuas presenta baja variabilidad y que hay la misma cantidad de datos para cada especie de flor.

Entonces, para el desarrollo del ejercio se va a considerar 120 obaservaciones como datos de entrenamiento (traing data) y 30 como datos de prueba (test data). Dado que las escalas o unidades de medida en las que fueron medidas las variables predictoras X pueden afectar las estimaciones que se realizan en el anĂ­lisis discriminante por lo que se procede a su estandarizaciĂ³n de estas variables antes de iniciar con los anĂ¡lisis discriminante.

# Librerias implementadas
library(tidyverse)
library(caret)
library(dplyr)

# DefiniciĂ³n de los datos de entrenamiento y los de prueba.

obs<-createDataPartition(iris$Species,p = 0.8, list = FALSE)

train.data <- iris[obs, ]
test.data <- iris[-obs, ]
summary(train.data$Species)
##     setosa versicolor  virginica 
##         40         40         40
summary(test.data$Species)
##     setosa versicolor  virginica 
##         10         10         10
# EstadarizaciĂ³n de las variables de forma conjunta.

preproc.param <- train.data %>% preProcess(method = c("center", "scale"))

train.est <- preproc.param %>% predict(train.data)
test.est <- preproc.param %>% predict(test.data)

AnĂ¡lisis Discriminante Lineal.

Ahora bien, como se sabe para la aplicaciĂ³n de este anĂ¡lisis es necesario que las varaibles predictoras (X-Input) sea normales, entonces para comprobar tal supuesto se procede a mirar como es el grafico QQplot de estas variables.

# Graficos QQ-plot para probar normalidad

par(mfrow=c(2,2))
for(i in 1:4){
qqnorm(iris[,i], pch = 1, frame = FALSE)
qqline(iris[,i], col = "steelblue", lwd = 2)
}

De estos Graficos podemos apreciar que las mediciones de ancho y largo del sepalo presenta un comportamiento normal, mientras, que las mediciones de ancho y largo del petalo no estarian presentado dicho comportamiento tan marco como en el otro caso, aĂºn asi para efectos practicos e ilustrativos se asumen normales con la misma matriz de covarianzas. AdĂ©mas de comprobar la normalidad se debe garantizar que entre la categorias de la variable de respuesta categorica G las variables X tengan diferentes medias pero igual matriz de covarianzas, para este caso tambien vamos a asumir este factor. Cabe mencionar, que tambien se debe considerar la precencia de puntos atipicos y de haberlos se debe eliminar.

Entonces, se tiene que un anĂ¡lisis discrimininante lineal (LDA) en R se hace mediante la funciĂ³n lda() de la libreria MASS. Se ajusta el modelo de regresiĂ³n del LDA con los datos de entrenamiento estandarizados (train.est)

library(MASS)

# Ajuste del Modelo

modelo<- lda(Species~.,data = train.est)

modelo
## Call:
## lda(Species ~ ., data = train.est)
## 
## 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
# Predicciones de las Clases

prediccion<- modelo %>% predict(test.est)

names(prediccion)
## [1] "class"     "posterior" "x"
# PrecisiĂ³n del modelo

mean(prediccion$class==test.est$Species)
## [1] 0.9666667
#Grafico de la clasificaciĂ³n

lda.data <- cbind(train.est, predict(modelo)$x)

ggplot(lda.data, aes(LD1, LD2)) +
  geom_point(aes(color = Species))

Entonces se tiene que el modelo de clasificaciĂ³n predice correctamente el 96.6666667% de las veces.

AnĂ¡lisis Discriminante Flexible.

A diferencia del LDA este anĂ¡lisis como se base en la una regrciĂ³n no parametrica no necesita que las variables X-input sean normales pero si que esten estandarizadas para evitar problemas de estiamciĂ³n del modelo clasificador. R para la realizaciĂ³n de este anĂ¡lisis implementa la funciĂ³n fda() de la libreria mda la cual cuenta con 4 opciones para la funciĂ³n de regreciĂ³n a considerar polyreg, mars, bruto y penalizada (gen.ridge), para este caso se va aconsiderar regreciĂ³n bruto. Entonces se tiene:

library(mda)

# Ajuste del modelo FDA

modelo.fda<- fda(Species~.,data = train.est,method = bruto)

modelo.fda
## Call:
## fda(formula = Species ~ ., data = train.est, method = bruto)
## 
## Dimension: 2 
## 
## Percent Between-Group Variance Explained:
##  v1  v2 
##  98 100 
## 
## Degrees of Freedom (per dimension): 12.85484 
## 
## Training Misclassification Error: 0.03333 ( N = 120 )
# PredicciĂ³n de las Clases

predicted.classes <- modelo.fda %>% predict(test.est)

# PrecisiĂ³n del modelo

mean(predicted.classes == test.est$Species)
## [1] 0.9333333
# Grafico de la clasificaciĂ³n
plot(modelo.fda,train.est[,1:4],group = "pred", col=c("red","blue","green"))

Aqui la precisiĂ³n del modelo es del 93.3333333%, es decir, el modelo clasificara correctamente dicho porcentaje de veces.

Ejemplo Vowel

La base de datos vowel consiste en un conjunto de datos de pronunciaciĂ³n de 11 vocales foneticas diferentes, dichas vocales fueron codificadas como 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10. La base de datos tambien cuenta con una variable que indica quien de los 14 participantes pronuncio la vocal junto con el sexo de ellos, ademĂ¡s midieron 10 variables continuas correspondientes a unos atributos medidos por los investigadores (no sabe a que hacen referencia pues los autores en la pagina de descarga de la base de datos no expefican nada al respecto). La base de datos cuenta con 990 observaciĂ³nes de las cuales 572 se van a considerar como training data y 418 como test data, por parte de la variables, se considera al tipo de vocal como variable de respuesta categorica G y las variables de atributos con las variables predictoras X. Se presenta entonces las siguientes descriptivas del total de las observaciones:

# Semilla
set.seed(123)

# Carga de los datos

vowel<-read.table("C:/Users/cesar/OneDrive/Escritorio/Seminario/Segunda Ronda/vowel.dat",sep=",")
vowel<-vowel[,-(1:3)]
colnames(vowel)<-c("F0","F1", "F2", "F3", "F4", "F5", "F6", "F7", "F8", "F9","Vocal")

vowel$Vocal<-as.factor(vowel$Vocal)
# Estadisticas descriptivas basicas

summary(vowel)
##        F0               F1               F2                 F3         
##  Min.   :-5.211   Min.   :-1.274   Min.   :-2.48700   Min.   :-1.4090  
##  1st Qu.:-3.888   1st Qu.: 1.052   1st Qu.:-0.97575   1st Qu.:-0.0655  
##  Median :-3.146   Median : 1.877   Median :-0.57250   Median : 0.4335  
##  Mean   :-3.204   Mean   : 1.882   Mean   :-0.50777   Mean   : 0.5155  
##  3rd Qu.:-2.603   3rd Qu.: 2.738   3rd Qu.:-0.06875   3rd Qu.: 1.0960  
##  Max.   :-0.941   Max.   : 5.074   Max.   : 1.43100   Max.   : 2.3770  
##                                                                        
##        F4                F5                F6                  F7          
##  Min.   :-2.1270   Min.   :-0.8360   Min.   :-1.537000   Min.   :-1.29300  
##  1st Qu.:-0.7690   1st Qu.: 0.1960   1st Qu.:-0.307000   1st Qu.:-0.09575  
##  Median :-0.2990   Median : 0.5520   Median : 0.022000   Median : 0.32800  
##  Mean   :-0.3057   Mean   : 0.6302   Mean   :-0.004365   Mean   : 0.33655  
##  3rd Qu.: 0.1695   3rd Qu.: 1.0285   3rd Qu.: 0.296500   3rd Qu.: 0.77000  
##  Max.   : 1.8310   Max.   : 2.3270   Max.   : 1.403000   Max.   : 2.03900  
##                                                                            
##        F8                 F9               Vocal    
##  Min.   :-1.61300   Min.   :-1.68000   0      : 90  
##  1st Qu.:-0.70400   1st Qu.:-0.54800   1      : 90  
##  Median :-0.30250   Median :-0.15650   2      : 90  
##  Mean   :-0.30298   Mean   :-0.07134   3      : 90  
##  3rd Qu.: 0.09375   3rd Qu.: 0.37100   4      : 90  
##  Max.   : 1.30900   Max.   : 1.39600   5      : 90  
##                                        (Other):450
ds<-sqrt(diag(var(vowel[,1:10])))
ds
##        F0        F1        F2        F3        F4        F5        F6        F7 
## 0.8689872 1.1752720 0.7119483 0.7592613 0.6646023 0.6038711 0.4619268 0.5733020 
##        F8        F9 
## 0.5701616 0.6039855

Entonces mis datos de entrena miento y de prueba se especifican y estandarizan a continuaciĂ³n, considerando para los anĂ¡lisis los estandarizados para evitar proablemas de estimaciĂ³n debido a la influencia de las unidades de medidas de las variables predictoras.

# DefiniciĂ³n de los datos de entrenamiento

ods <- vowel$Vocal %>% createDataPartition(p = 562/990, list = FALSE)

train.data <-vowel [ods, ]
test.data <-vowel[-ods, ]

# EstarizaciĂ³n conjunta de las variables

preproc.param <- train.data %>% preProcess(method = c("center", "scale"))

train.ets <- preproc.param %>% predict(train.data)
test.ets <- preproc.param %>% predict(test.data)

ComparaciĂ³n de los AnĂ¡lisis Discriminantes

Ahora, con este ejemplo lo que se quiere es hacer una comparaciĂ³n mas notoria entre los metodos de anĂ¡lisis discriminantes considerados (Lineal, Flexible y Penalizado) por lo que para el caso del lineal se van asumir que las variables predictoras son normales con igual matriz de covarianza, para el caso del flexible se considera regreciĂ³n bruto y mars.

# LDA

modelo.lda <- lda(Vocal~., data = train.ets)

# Predicciones

pdclda <- modelo.lda %>% predict(test.ets)

# FDA MARS

modelo.fda1 <- fda(Vocal~., data = train.ets,method=mars)

# Predicciones

pdcfda1 <- modelo.fda1 %>% predict(test.ets)

# FDA bruto

modelo.fda2 <- fda(Vocal~., data = train.ets,method=bruto)

# Predicciones

pdcfda2 <- modelo.fda2 %>% predict(test.ets)

# PDA

modelo.pda <- fda(Vocal~., data = train.ets,method=gen.ridge)

# Predicciones

pdcpda <- modelo.pda %>% predict(test.ets)

Entonces podemos apreciar que las preciones de cada uno de los modelos de clasificaciĂ³n serian 60.2870813% para el LDA, 66.2679426% para el FDA MARS, 75.3588517% para el FDA Bruto y 60.2870813% para el PDA.

Ahora bien, graficamente podemos apreciar que el FDA con regresiĂ³n Bruto es que mejor discrimina las categorias y mejor clasifica a los datos.

par(mfrow=c(1,2))
#LDA
plot(pdclda$x[,1],pdclda$x[,2],col=c("blue4","brown3","darkgoldenrod1","darkorchid4","forestgreen","darkturquoise","lightyellow4","slategray4","seagreen3","yellow","tan3")[vowel$Vocal],xlab = "LDA1",ylab = "LDA2",main = "LDA")

#FDA1
plot(modelo.fda1,data = train.ets[,1:10],group = "pred",col=c("blue4","brown3","darkgoldenrod1","darkorchid4","forestgreen","darkturquoise","lightyellow4","slategray4","seagreen3","yellow","tan3")[vowel$Vocal])

par(mfrow=c(1,2))
#FDA2
plot(modelo.fda2,data = train.ets[,1:10],group = "pred",col=c("blue4","brown3","darkgoldenrod1","darkorchid4","forestgreen","darkturquoise","lightyellow4","slategray4","seagreen3","yellow","tan3")[vowel$Vocal])

#PDA
plot(modelo.pda,data = train.ets[,1:10],group = "pred",col=c("blue4","brown3","darkgoldenrod1","darkorchid4","forestgreen","darkturquoise","lightyellow4","slategray4","seagreen3","yellow","tan3")[vowel$Vocal])