I. Se tiene la base de datos denominada Academico , que podemos descargar desde la plataforma de Moodle. Esta base de datos fue creada para discriminar la admisión de estudiantes a la Escuela Graduada de Negocios. Se hicieron tres pruebas (Economía, Finanzas y Matemáticas), y se cuenta con la puntuación de la universidad de procedencia y su GPA. La codificación de los grupos es (1) admitido y (2) no admitido. Realice un análisis de discriminante para evaluar si la función discriminante que se esta aplicando es adecuada.

Primero descargamos y modificamos la base de datos:

library(DT)
library(readxl)
Academico <- read.csv("Academico.csv")
Data <- Academico[,-1]
Data_mod <- Data[,-4:-5]
str(Data)
## 'data.frame':    54 obs. of  6 variables:
##  $ Química     : num  52.2 50.9 52.6 52.6 50.5 ...
##  $ Física      : num  51.2 51.5 48.9 52.5 51.9 ...
##  $ Matemáticas : num  49.9 52.7 47.9 48.6 50.1 ...
##  $ Clas_colegio: int  6 7 5 5 5 6 7 7 6 5 ...
##  $ GPA         : num  3.04 3.04 3.44 3.44 2.88 3.28 3.76 2.8 3.28 3.36 ...
##  $ Grupo       : int  1 1 1 1 1 1 1 1 1 1 ...
class(Data)
## [1] "data.frame"
Data1 <- as.data.frame(Data_mod)
class(Data1)
## [1] "data.frame"
datatable(Data1)

Decidimos eliminar las variables 4 “Clas_colegio” y 5 “GPA” de la base de datos debido a que generaban gráficas con magnitudes significativamente diferentes que dificultaban la interpretación de los datos.

Notamos que hay algunas variables categóricas que se clasifican como numéricas, por lo que debemos hacer ajustes:

Data$Grupo <- factor(Data$Grupo,levels=c(1,2),labels=c("admitido","no admitido"))
Data$Clas_colegio <- factor(Data$Clas_colegio,levels = c(2,3,4,5,6,7))
str(Data)
## 'data.frame':    54 obs. of  6 variables:
##  $ Química     : num  52.2 50.9 52.6 52.6 50.5 ...
##  $ Física      : num  51.2 51.5 48.9 52.5 51.9 ...
##  $ Matemáticas : num  49.9 52.7 47.9 48.6 50.1 ...
##  $ Clas_colegio: Factor w/ 6 levels "2","3","4","5",..: 5 6 4 4 4 5 6 6 5 4 ...
##  $ GPA         : num  3.04 3.04 3.44 3.44 2.88 3.28 3.76 2.8 3.28 3.36 ...
##  $ Grupo       : Factor w/ 2 levels "admitido","no admitido": 1 1 1 1 1 1 1 1 1 1 ...

Determinamos la probabilidad (estimación de las probabilidades a priori):

prop.table(table(Data$Grupo))
## 
##    admitido no admitido 
##   0.6296296   0.3703704

Ahora analizamos las matrices de covarianzas para las dos poblaciones, para poder determinar si se pueden suponer iguales o diferentes, para esto utilizamos la función boxM :

library(biotools)
boxM(Data1[,1:3], grouping = Data1[,4])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Data1[, 1:3]
## Chi-Sq (approx.) = 2.0783, df = 6, p-value = 0.9124

Con ayuda del p-value, podemos suponer que las matrices de covarianza son iguales, por lo que es conveniente utilizar una regla de discriminación lineal (LDA).

Ahora, para aplicar LDA en R, se puede utilizar la función lda() de la librería MASS:

library(MASS)
modelo1 <- lda(Grupo ~ .,Data1)
modelo1
## Call:
## lda(Grupo ~ ., data = Data1)
## 
## Prior probabilities of groups:
##         1         2 
## 0.6296296 0.3703704 
## 
## Group means:
##    Química   Física Matemáticas
## 1 51.64941 50.07824    49.98618
## 2 46.41950 45.41100    45.42450
## 
## Coefficients of linear discriminants:
##                     LD1
## Química     -0.11629608
## Física      -0.07169761
## Matemáticas -0.01478072

Analizamos el error de predicción para este modelo:

Predicciones <- predict(modelo1)
table(Data1$Grupo,Predicciones$class, dnn = c("Clase real", "Clase predicha"))
##           Clase predicha
## Clase real  1  2
##          1 28  6
##          2 10 10
tasa_error <- mean(Data1$Grupo != Predicciones$class)
round(tasa_error*100,1)
## [1] 29.6

Podemos observar que hubieron 16 errores, en cuanto al grupo 1 “admitidos” el margen de error es de 6/34, mientras que para el grupo 2 “no admitidos” el margen de error es de 10/20. En total, hubieron 16 errores de las 54 observaciones. Es decir, que el modelo creado tiene un margen de error del 29.6%, lo que puede no ser ideal para discriminar.

Para poder ver gráficamente estos resultados del modelo creado con “LDA” creamos el siguiente gráfico:

par(mfrow = c(1, 1), mar = c(2, 2, 2, 2))
library(klaR)
Data1$Grupo <- as.factor(Data1$Grupo)
partimat(Grupo ~ ., Data1, method = "lda", image.colors = c("lightgreen", "skyblue2"), col.mean = "red", plot.matrix =T)

En conclusión, tras evaluar los resultados arrojados por el modelo creado, determinamos que la función discriminante “LDA” no es la más adecuada para hacer este análisis discriminante.

II. Aplique la técnica de análisis de discriminante a la base de datos iris de R y compruebe si es sencillo discriminar cada especie.

Primero descargamos nuestra base de datos:

Data2 <- iris
str(Data2)
## 'data.frame':    150 obs. of  5 variables:
##  $ Sepal.Length: num  5.1 4.9 4.7 4.6 5 5.4 4.6 5 4.4 4.9 ...
##  $ Sepal.Width : num  3.5 3 3.2 3.1 3.6 3.9 3.4 3.4 2.9 3.1 ...
##  $ Petal.Length: num  1.4 1.4 1.3 1.5 1.4 1.7 1.4 1.5 1.4 1.5 ...
##  $ Petal.Width : num  0.2 0.2 0.2 0.2 0.2 0.4 0.3 0.2 0.2 0.1 ...
##  $ Species     : Factor w/ 3 levels "setosa","versicolor",..: 1 1 1 1 1 1 1 1 1 1 ...
datatable(Data2)

No es necesario hacer modificaciones a la base de datos.

Ahora, podemos explorar gráficamente los datos. Primero, creamos una visualización de cada par de variables:

pairs(Data2[,-5], col = c("blue", "green3","red")[Data2[,5]], pch = 19)

Estas visualizaciones nos permiten entender que las variables con información sobre el largo y ancho de los pétalos respecto a las otras dos variables, ayuda a discriminar bastante bien los tres grupos.

Decidimos hacer dos gráficas en 3D, la primera con la variable del ancho del sépalo respecto a las dos variables sobre pétalos, y la segunda sobre la variable del largo del sépalo respecto a las mismas dos variables:

library(plotly)

fig <- plot_ly(Data2, x = ~Sepal.Length, y = ~Petal.Length, z = ~Petal.Width, color = ~Species, colors = c("blue", "green3","red"))
fig <- fig %>% add_markers()
fig
fig2 <- plot_ly(Data2, x = ~Sepal.Width, y = ~Petal.Length, z = ~Petal.Width, color = ~Species, colors = c("blue", "green3","red"))
fig_2 <- fig2 %>% add_markers()
fig_2

Pudimos notar que la especie “setosa” tiene carasterísticas muy particulares, por lo que es sencillo discriminarla.

Ahora analizamos las matrices de covarianzas para las tres poblaciones, para poder determinar si se pueden suponer iguales o diferentes:

library(biotools)
boxM(Data2[,-5], grouping = Data2[, 5])
## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  Data2[, -5]
## Chi-Sq (approx.) = 140.94, df = 20, p-value < 2.2e-16

Debido a que el p-value es de 2.2e-16, podemos suponer que las matrices de covarianzas no son iguales, por lo que deberíamos aplicar clasificación “QDA”.

Antes de crear el modelo con clasificación QDA, analizamos la normalidad de las variables utilizando las librería nortest y MVN para saber si los resultados se verán afectados por la falta de este supuesto:

# Lilli
library(nortest)

L1 <- lillie.test(iris$Sepal.Length) 
L2 <- lillie.test(iris$Sepal.Width)   
L3 <- lillie.test(iris$Petal.Length)   
L4 <- lillie.test(iris$Petal.Width)  

L1
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  iris$Sepal.Length
## D = 0.088654, p-value = 0.005788
L2
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  iris$Sepal.Width
## D = 0.10566, p-value = 0.0003142
L3
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  iris$Petal.Length
## D = 0.19815, p-value = 7.901e-16
L4
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  iris$Petal.Width
## D = 0.17283, p-value = 7.33e-12
# Anderson-Darling

library(nortest)
AD1 <- ad.test(iris$Sepal.Length)
AD2 <- ad.test(iris$Sepal.Width)
AD3 <- ad.test(iris$Petal.Length)
AD4 <- ad.test(iris$Petal.Width)

AD1
## 
##  Anderson-Darling normality test
## 
## data:  iris$Sepal.Length
## A = 0.8892, p-value = 0.02251
AD2
## 
##  Anderson-Darling normality test
## 
## data:  iris$Sepal.Width
## A = 0.90796, p-value = 0.02023
AD3
## 
##  Anderson-Darling normality test
## 
## data:  iris$Petal.Length
## A = 7.6785, p-value < 2.2e-16
AD4
## 
##  Anderson-Darling normality test
## 
## data:  iris$Petal.Width
## A = 5.1057, p-value = 1.125e-12
# Cramer-Von Mises

library(nortest)
CVM1 <- cvm.test(iris$Sepal.Length)    
CVM2 <- cvm.test(iris$Sepal.Width)
CVM3 <- cvm.test(iris$Petal.Length)
CVM4 <- cvm.test(iris$Petal.Width)

CVM1
## 
##  Cramer-von Mises normality test
## 
## data:  iris$Sepal.Length
## W = 0.1274, p-value = 0.04706
CVM2
## 
##  Cramer-von Mises normality test
## 
## data:  iris$Sepal.Width
## W = 0.18065, p-value = 0.009336
CVM3
## 
##  Cramer-von Mises normality test
## 
## data:  iris$Petal.Length
## W = 1.2223, p-value = 7.37e-10
CVM4
## 
##  Cramer-von Mises normality test
## 
## data:  iris$Petal.Width
## W = 0.72156, p-value = 4.338e-08

Para las pruebas de normalidad multivariadas debemos eliminar la variable “Species” para que solo se trabaje con las variables númericas, por lo que modificamos nuestra base de datos:

Data_num <- Data2[,-5]
head(Data_num)
##   Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1          5.1         3.5          1.4         0.2
## 2          4.9         3.0          1.4         0.2
## 3          4.7         3.2          1.3         0.2
## 4          4.6         3.1          1.5         0.2
## 5          5.0         3.6          1.4         0.2
## 6          5.4         3.9          1.7         0.4
# Royston
library(MVN)
Royston <- mvn(Data_num, mvnTest = "royston")
Royston$multivariateNormality
##      Test        H      p value MVN
## 1 Royston 50.39667 3.098229e-11  NO
# Henze-Zirkler
library(MVN)
HZ <- mvn(Data_num, mvnTest = "hz")
HZ$multivariateNormality
##            Test       HZ p value MVN
## 1 Henze-Zirkler 2.336394       0  NO

Pruebas gráficas de normalidad:

# Q-Q plot
par(mfrow = c(2, 2))
qqnorm(Data2[,1], pch = 1, main="Q-Q plot x1")
qqline(Data2[,1], col = "blue", lwd = 2)
qqnorm(Data2[,2], pch = 1, main="Q-Q plot x2")
qqline(Data2[,2], col = "blue", lwd = 2)
qqnorm(Data2[,3], pch = 1, main="Q-Q plot x3")
qqline(Data2[,3], col = "blue", lwd = 2)
qqnorm(Data2[,4], pch = 1, main="Q-Q plot x4")
qqline(Data2[,4], col = "blue", lwd = 2)

# Histogramas
par(mfrow = c(2, 2))
x<- seq(-4, 4, length = 100)
hist(Data2[,1], prob = TRUE, xlab=" ", main="Histograma x1")
lines(x,dnorm(x),col="red",lwd=2)
hist(Data2[,2], prob = TRUE,  xlab=" ", main="Histograma x2")
lines(x,dnorm(x),col="red",lwd=2)
hist(Data2[,3], prob = TRUE,  xlab=" ", main="Histograma x3")
lines(x,dnorm(x),col="red",lwd=2)
hist(Data2[,4], prob = TRUE,  xlab=" ", main="Histograma x4")
lines(x,dnorm(x),col="red",lwd=2)

Tras analizar las distintas pruebas de normalidad univariadas, multivariadas y gráficas realizadas, podemos determinar que ninguna de las variables cumple con el supuesto de normalidad, por lo que los resultados que nos arroje el modelo con clasificación QDA puede verse afectado.

La prueba de homogeneidad de varianzas sugiere que las varianzas entre los grupos no son igual, lo que indica que usar el método LDA podría no ser adecuado y el QDA sería mejor. Sin embargo, esta prueba es muy sensible a la falta de normalidad en los datos, lo que puede hacer que parezca que la diferencia está en la matriz de covarianza cuando en realidad es por la falta de normalidad. Por lo tanto, asumiremos que la matriz de covarianza es constante y que LDA puede funcionar bien para clasificar los datos. Vamos a evaluar qué tan precisa es esta suposición al probar el modelo y en las conclusiones explicaremos esta decisión.

Ahora si creamos el modelo con clasificación LDA:

library(MASS)
modelo2 <- lda(Species ~ ., data= Data2)
modelo2
## Call:
## lda(Species ~ ., data = Data2)
## 
## 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

Hacemos la evaluación del error de clasificación:

predicciones <- predict(modelo2)
table(Data2$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

Podemos observar que en cuanto a la especie “versicolor” el margen de error es de 2/50, mientras que para la especie “virginica” el margen de error es de 1/50. Es decir, que solo hubieron 3 errores de las 150 observaciones. Es decir, que el modelo creado solo tiene un margen de error del 2%, lo que es muy bueno.

Finalmente podemos visualizar como se realizó la clasificación para cada par de variables:

library(klaR)
partimat(Species ~ .,Data2, method = "qda", image.colors = c("lightgreen", "skyblue2", "orange"), col.mean = "red")

En conclusión, pudimos determinar que con el modelo creado con la clasificacion “LDA”, si es fácil discriminar cada especie de las flores en la base de datos (especialmente la especie “setosa”), con un 98% de confianza.