En este primer caso de estudio, el conjunto de datos de vino, tenemos 13 concentraciones químicas que describen muestras de vino de tres cultivares.
library(car)
## Loading required package: carData
library(rattle)
## Loading required package: tibble
## Loading required package: bitops
## Rattle: A free graphical interface for data science with R.
## Versión 5.4.0 Copyright (c) 2006-2020 Togaware Pty Ltd.
## Escriba 'rattle()' para agitar, sacudir y rotar sus datos.
#install.packages('rattle')
data(wine, package='rattle')
attach(wine)
head(wine)
scatterplotMatrix(wine[2:6])
El propósito del análisis discriminante lineal (LDA) en este ejemplo es encontrar las combinaciones lineales de las variables originales (las 13 concentraciones químicas aquí) que dan la mejor separación posible entre los grupos (cultivares de vino aquí) en nuestro conjunto de datos. El análisis discriminante lineal también se conoce como “análisis discriminante canónico”, o simplemente “análisis discriminante”.
Si queremos separar los vinos por cultivar, los vinos provienen de tres cultivares diferentes, por lo que el número de grupos G = 3 y el número de variables es 13 (13 concentraciones de sustancias químicas; p = 13). El número máximo de funciones discriminantes útiles que pueden separar los vinos por cultivar es el mínimo de G- 1 yp, por lo que en este caso es el mínimo de 2 y 13, que es 2. Así, podemos encontrar como máximo 2 útiles funciones discriminantes para separar los vinos por cultivar, utilizando las 13 variables de concentración química.
Puede realizar un análisis discriminante lineal utilizando la función “lda()” del paquete R MASS. Para usar esta función, primero necesitamos instalar el paquete MASS de R.
# install.packages('MASS')
library(MASS)
wine.lda <- lda(Type ~ ., data=wine)
Para obtener los valores de las cargas(loadings) de las funciones discriminantes para los datos del vino(winw data), podemos escribir:
wine.lda
## Call:
## lda(Type ~ ., data = wine)
##
## Prior probabilities of groups:
## 1 2 3
## 0.3314607 0.3988764 0.2696629
##
## Group means:
## Alcohol Malic Ash Alcalinity Magnesium Phenols Flavanoids
## 1 13.74475 2.010678 2.455593 17.03729 106.3390 2.840169 2.9823729
## 2 12.27873 1.932676 2.244789 20.23803 94.5493 2.258873 2.0808451
## 3 13.15375 3.333750 2.437083 21.41667 99.3125 1.678750 0.7814583
## Nonflavanoids Proanthocyanins Color Hue Dilution Proline
## 1 0.290000 1.899322 5.528305 1.0620339 3.157797 1115.7119
## 2 0.363662 1.630282 3.086620 1.0562817 2.785352 519.5070
## 3 0.447500 1.153542 7.396250 0.6827083 1.683542 629.8958
##
## Coefficients of linear discriminants:
## LD1 LD2
## Alcohol -0.403399781 0.8717930699
## Malic 0.165254596 0.3053797325
## Ash -0.369075256 2.3458497486
## Alcalinity 0.154797889 -0.1463807654
## Magnesium -0.002163496 -0.0004627565
## Phenols 0.618052068 -0.0322128171
## Flavanoids -1.661191235 -0.4919980543
## Nonflavanoids -1.495818440 -1.6309537953
## Proanthocyanins 0.134092628 -0.3070875776
## Color 0.355055710 0.2532306865
## Hue -0.818036073 -1.5156344987
## Dilution -1.157559376 0.0511839665
## Proline -0.002691206 0.0028529846
##
## Proportion of trace:
## LD1 LD2
## 0.6875 0.3125
Esto significa que la primera función discriminante es una combinación lineal de las variables: -0,403* Alcohol + 0,165* Málico+… -0,003*Prolina. Por conveniencia, el valor de cada función discriminante (por ejemplo, la primera función discriminante) se escala de modo que su valor medio sea cero y su varianza sea uno.
La “proporción de traza” que se imprime cuando escribe “wine.lda” (la variable devuelta por la función lda()) es la separación porcentual lograda por cada función discriminante. Por ejemplo, para los datos del vino obtenemos los mismos valores que acabamos de calcular (68,75% y 31,25%).
Un histograma apilado de los valores LDA
Una buena forma de mostrar los resultados de un análisis discriminante lineal (LDA) es hacer un histograma apilado de los valores de la función discriminante para las muestras de diferentes grupos (diferentes cultivares de vino en nuestro ejemplo).
Podemos hacer esto usando la función “ldahist ()” en R. Por ejemplo, para hacer un histograma apilado de los valores de la primera función discriminante para muestras de vino de los tres cultivares de vino diferentes, escribimos:
wine.lda.values <- predict(wine.lda)
ldahist(data = wine.lda.values$x[,1], g=Type)
La segunda función discriminante separa esos cultivares, haciendo un histograma apilado de los valores de la segunda función discriminante:
ldahist(data = wine.lda.values$x[,2], g=Type)
Diagramas de dispersión de las funciones discriminantes
Podemos obtener una gráfica de dispersión de las dos mejores funciones discriminantes, con los puntos de datos etiquetados por cultivar, escribiendo:
plot(wine.lda.values$x[,1],wine.lda.values$x[,2]) # make a scatterplot
text(wine.lda.values$x[,1],wine.lda.values$x[,2],Type,cex=0.7,pos=4,col="red") # add labels
De el scatterplot de la primeras dos funciones discriminantes funciones discriminantes, podemos ver que los vinos de los tres cultivares están bien separados en el diagrama de dispersión. La primera función discriminante (eje x) separa muy bien los cultivares 1 y 3, pero no separa perfectamente los cultivares 1 y 3, o los cultivares 2 y 3.
La segunda función discriminante (eje y) logra una separación bastante buena de los cultivares 1 y 3, y los cultivares 2 y 3, aunque no es totalmente perfecta.
Para lograr una muy buena separación de los tres cultivares, sería mejor usar la primera y la segunda función discriminante juntas, ya que la primera función discriminante puede separar muy bien los cultivares 1 y 3, y la segunda función discriminante puede separar los cultivares 1 y 2, y los cultivares 2 y 3, razonablemente bien.
El conjunto de datos proporciona datos de admisión para los solicitantes de escuelas de posgrado en negocios. El objetivo es utilizar los puntajes GPA y GMAT para predecir la probabilidad de admisión (admitir, no admitir y límite).
url <- 'http://www.biz.uiowa.edu/faculty/jledolter/DataMining/admission.csv'
admit <- read.csv(url)
head(admit)
Podemos hacer un diagrama de dispersión con los datos disponibles:
str(admit)
## 'data.frame': 85 obs. of 3 variables:
## $ GPA : num 2.96 3.14 3.22 3.29 3.69 3.46 3.03 3.19 3.63 3.59 ...
## $ GMAT: int 596 473 482 527 505 693 626 663 447 588 ...
## $ De : chr "admit" "admit" "admit" "admit" ...
admit$De=as.factor(admit$De)
adm=as.data.frame(admit)
head(adm)
attach(adm)
plot(GPA,GMAT,col=De)
Podemos hacer un diagrama de dispersión con los datos a mano: Comencemos haciendo el LDA, observe que en este caso tenemos 3 clases:
library(MASS)
## linear discriminant analysis
m1=lda(De~.,adm)
m1
## Call:
## lda(De ~ ., data = adm)
##
## Prior probabilities of groups:
## admit border notadmit
## 0.3647059 0.3058824 0.3294118
##
## Group means:
## GPA GMAT
## admit 3.403871 561.2258
## border 2.992692 446.2308
## notadmit 2.482500 447.0714
##
## Coefficients of linear discriminants:
## LD1 LD2
## GPA 5.008766354 1.87668220
## GMAT 0.008568593 -0.01445106
##
## Proportion of trace:
## LD1 LD2
## 0.9673 0.0327
predict(m1,newdata=data.frame(GPA=3.21,GMAT=497),type="class")
## $class
## [1] admit
## Levels: admit border notadmit
##
## $posterior
## admit border notadmit
## 1 0.5180421 0.4816015 0.0003563717
##
## $x
## LD1 LD2
## 1 1.252409 0.318194
Discriminante Cuadratico
m2=qda(De~.,adm)
m2
## Call:
## qda(De ~ ., data = adm)
##
## Prior probabilities of groups:
## admit border notadmit
## 0.3647059 0.3058824 0.3294118
##
## Group means:
## GPA GMAT
## admit 3.403871 561.2258
## border 2.992692 446.2308
## notadmit 2.482500 447.0714
¿Qué modelo es mejor? Para responder a esta pregunta, evaluamos el análisis discriminante lineal seleccionando aleatoriamente a 60 de 85 estudiantes, estimando los parámetros de los datos de entrenamiento y clasificando a los 25 estudiantes restantes de la muestra reservada. Repetimos esto 100 veces.
n=85
nt=60
neval=n-nt
rep=100
### LDA
set.seed(123456789)
errlin=dim(rep)
for (k in 1:rep) {
train=sample(1:n,nt)
## linear discriminant analysis
m1=lda(De~.,adm[train,])
predict(m1,adm[-train,])$class
tablin=table(adm$De[-train],predict(m1,adm[-train,])$class)
errlin[k]=(neval-sum(diag(tablin)))/neval
}
merrlin=mean(errlin)
merrlin
## [1] 0.0916
### QDA
set.seed(123456789)
errqda=dim(rep)
for (k in 1:rep) {
train=sample(1:n,nt)
## quadratic discriminant analysis
m1=qda(De~.,adm[train,])
predict(m1,adm[-train,])$class
tablin=table(adm$De[-train],predict(m1,adm[-train,])$class)
errqda[k]=(neval-sum(diag(tablin)))/neval
}
merrqda=mean(errlin)
merrqda
## [1] 0.0916