Aplicación de algoritmos clasificadores al dataset de Iris.
AD (Arbol de Decision),
LDA (Lineal Discriminant Analisys),
QDA (Quadratic Discriminant Analisys),
NB(Naive-Bayes o Discriminante Bayesiano),
k-NN (k-Nearest Neighbors)
Pilar Rocío Sayán Mejía
19 Noviembre, 2024
✔ General
INICIO
Se intenta mantener el codigo en secciones bien definidas.
- PREPARACION
- Cargar el Dataset
- Dividirlo en conjunto de Entrenamiento y Prueba
- Crear un subconjunto de valores normalizados
- PRUEBA de ATRIBUTOS
- Analisis univariante
- Analisis multivariante
- ALGORITMOs de Clasificacion
- Los mencionados
- METRICAS de RENDIMIENTO y CALIDAD
- Tabla de confusión
- Visualizaciones de confusión
- Métricas del paquete caret para cada clasificador
Las funciones de ayuda son pequeñas funciones utilizadas en la secuencia de análisis con el objetivo de reducir partes del proceso y hacerlo mas legible. Es posible que no sean aplicables al 100% de los escenarios, Y que no se hayan pulido al 100%.
Los algoritmos:
- AD (Arboles de Decisión) consiste en dividir iterativamente el dataset usando una función d ponderación o coste para encontrar las variables que aportan mayor información.
- LDA () produce límites de decisión lineales a partir de una combinación lineal de las variables, lo que se traduce en menor flexibilidad y por lo tanto menor problema de varianza, pero si la separación de los grupos no es lineal, tendrá una varianza grande.
- QDA () produce límites cuadráticos y por lo tanto curvos, lo que
aporta mayor flexibilidad en la linea de decisión permitiendo ajustarse
mejor a los datos.
- LDA tiende a conseguir mejores clasificaciones que QDA cuando hay pocas observaciones con las que entrenar al modelo
- QDA es mas apropiado si no es asumible que existe una matriz de covarianza común entre clases y si se dispone de una mayor cantidad de datos.
- NB (Naive-Bayes) …
- k-NN (k-Nearest Neighbors) es un algoritmo que se ejecuta ante cada predicción calculando las distancias de los vecinos y eligiendo los K mas cercanos.
- NN (Neural Networks), redes de elementos simples (neuronas artificiales) interconectadas masivamente en paralelo, adaptativos y con organización jerárquica similar a la de un sistema neuronal biológico.
Funciones de ayuda
# NORMALIDAD. Chequear con Histograma por clase
histo <- function( variable,...){
h = hist( variable, probability=TRUE, col=grey(0.8), ... )
u = variable %>% mean %>% round(4)
s = variable %>% sd %>% round(4)
x = variable %>% {seq( min(.), max(.), length.out=50 )}
lines( x, dnorm(x,u,s), col='red', lwd=2)
lines( density(variable), col='green' )
abline( v = u , lty=2, col='red' )
#abline( v = u+s , lty=2, col='red' )
#abline( v = u-s , lty=2, col='red' )
abline( v = median(variable) , lty=2, col='green' )
return(h)
}
histo_multi = function( datos, classes, mfcol=NULL ){
oldpar <- par(no.readonly=TRUE)
on.exit( par(oldpar) )
lev = levels(classes)
mfcol = if( is.null(mfcol) ) c( length(lev) , ncol(datos) ) else mfcol
par(mfcol = mfcol )
for ( k in 1:ncol(datos) ) {
varname <- names(datos)[k]
for ( i in 1:length(lev) ) {
x <- datos[ classes==lev[i], varname]
histo(x, main=paste( lev[i], " &\n", varname ), xlab="", ylab="" )
}
}
}
#histo_multi( iris[,-c(4,5)], iris[,5] ,c(3,3) )
##### NORMALIDAD. Chequear con Quantile Plot por clase
qqploto = function(x, title="varname" ){
qqnorm(x, main=title , xlab="", ylab="",
pch = 19 ) #, col = i + 1)
qqline(x, col="red")
abline( h=mean(x, na.rm=TRUE), col="blue")
abline( v=0, col="blue")
}
qqplot_multi = function( datos, classes ){
oldpar <- par(no.readonly=TRUE) # store 'par' setup
on.exit( par(oldpar) ) # recover 'par' setup
datos <- select( datos, where(is.numeric) )
lev = levels(classes)
par(mfcol = c( length(lev), ncol(datos)) )
for ( k in 1:ncol(datos) ) {
varname <- names(datos)[k]
x0 <- seq( min(datos[, k], na.rm=TRUE),
max(datos[, k], na.rm=TRUE),
length.out=50)
for ( i in 1:length(lev) ) {
x <- datos[ classes==lev[i], varname]
qqploto( x, paste( lev[i]," VS ", varname ) )
}
}
}
#qqplot_multi( iris[,-5], iris[,5] )
panel.cor <- function(x, y, digits = 3, prefix = "", cex.cor, ...)
{
usr <- par("usr")
par(usr = c(0, 1, 0, 1))
corrrr = cor(x, y)
if ( abs(corrrr) > 0.6 ) digits = digits-1
if ( abs(corrrr) < 0.1 ) digits = digits+1
if ( abs(corrrr) < 0.01 ) digits = digits+1
r <- max(abs(corrrr), 0.2)
txt <- format(c(corrrr, 0.123456789), digits = digits)[1]
txt <- paste0(prefix, txt)
if(missing(cex.cor)) cex.cor <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex.cor * r)
}
panel.hist <- function(x, ...)
{
usr <- par("usr")
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot=FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
# grafica sobre una grafica preexistente
rect(breaks[-nB], 0, breaks[-1], y, col='gray')
u = mean(x) %>% round(4)
s = sd(x) %>% round(4)
xn = x %>% {seq( min(.), max(.), length.out=50 )}
d = dnorm(xn,u,s)
lines( xn, d, col='red', lwd=2)
lines( density(x), col='green' )
# colors <- list(...)$col
# classes <- split(x, colors )
# colors <- unique(colors)
#for( i in 1:length(colors) )
# lines( density(classes[[i]]), col=colors[i] )
abline( v = u , lty=2, col='red' )
#abline( v = u+s , lty=2, col='red' )
#abline( v = u-s , lty=2, col='red' )
abline( v = median(x) , lty=2, col='green' )
}
ploto_pairs = function( datos, grupos, ...) {
oldpar <- par(no.readonly=TRUE) # store 'par' setup
on.exit( par(oldpar) ) # recover 'par' setup
pairs( datos,
col=rainbow( nlevels( factor(grupos) ) )[ factor(grupos) ],
lower.panel=panel.cor, diag.panel = panel.hist, upper.panel = panel.smooth, ... )
}
#ploto_pairs( iris[,-5], iris[,5] )
my.math.formula <- function( coef, vars=NULL )
paste0( c('+','-')[ as.numeric(coef<0)+1 ],
abs( round( coef ,3) ),
my.if.null( vars, names(coef) ), collapse=' ' )
my.math.normalize <- function(x) (x - min(x, na.rm=TRUE))/(max(x, na.rm=TRUE)-min(x, na.rm=TRUE))
# coef = c( x=4, y=-5, z=-6)
# my.math.formula( coef )
my.download_btn <- function( DF, oname, ...) downloadthis::download_this(DF,
output_name = oname, output_extension = ".csv",
button_label = "Download data as csv", button_type = "default",
has_icon = TRUE, icon = "fa fa-save" ) ✔ iris (con R)
El dataset Iris contiene mediciones de 150 flores de 3 especies diferentes de flores iris (virginica, setosa y versicolor) con 50 observaciones cada una y 4 variables: Largo y ancho del sépalo (sepal.length, sepal.width) y largo y ancho del pétalo (petal.length y petal.width), todas en centímetros. leer mas
Creada en 1936 por R.A. Fisher, es un dataset para clasificación multivariada que se puede encontrar en UCI Irvine Machine Learning Repository asi como tambien en los multiples lenguajes de programación para ciencia de datos.
Ilustracion de Flores de IRIS
• CARGAR
• RESAMPLE
• NORMALICE
El primer paso es cargar el dataset y dividirlo en ‘Entrenamiento’ y ‘Prueba’.
Y hacer los ajustes necesarios como posibles reducciones de variables y remoción de outliers si los hubiere.
En el siguiente código, se exponen 2 métodos para ello.
- Metodo 1: 1º Mezclar y luego segmentar
- Metodo 2: NO mezclar y segmentar aleatoriamente
Debido al pequeño conjunto de datos del dataset iris (150 flores), el
fraccionamiento produce un conjunto de prueba pequeño, de menos de 50
valores ( menos de 30 por cada flor ).
Debido a esto las diferencias entre los diferentes algoritmos no queda
inmediatamente evidente.
Para ser capaces de graficar las diferencias logradas en cada algoritmo se utiliza el conjunto completo de datos tanto para entrenamiento como para prueba.
set.seed(123)
# ==============================================================================
# ==== PREPARACION
# ==============================================================================
# ✔ CARGAR DATASET --------------------------------------------------------------
DF = iris[,]
Cidx = 5 # Class column index
Cname = names(DF)[ Cidx ] # Class column name
fmla = as.formula( paste(Cname,"~ .") ) # formula expresion
fmla2 = as.formula( paste("as.numeric(",Cname,")~ .") ) # formula expresion
# ✔ RESAMPLE -----------------------------------------------------------------
set.seed(123)
# METODO 1: primero mezclar y luego particionar
# DF <- DF[ sample(nrow(DF)) ,] # shuffle
# DF.test <- DF[ 1:(nrow(DF)*0.5) ,]
# DF.train <- DF[ (nrow(DF)*0.5+1):nrow(DF) ,]
# METODO 2: mantener orden y crear un vector booleano aleatorio
# DF.frac = sample( c(TRUE, FALSE), nrow(DF), replace=TRUE, prob=c( 0.7, 0.3 ) )
# DF.test = DF[ !DF.frac, ]
# DF.train = DF[ DF.frac, ]
# METODO 3: usar el mismo dataset para test y train
DF.test <- DF
DF.train <- DF
# VERIFICAR: si hay desbalance de clases:
ptable = rbind( table( DF.test[, Cidx] ) / nrow(DF.test),
table( DF.train[,Cidx] ) / nrow(DF.train) )
barplot( ptable, col=c(2,3), border=0, main="Check for Class unbalance", beside=TRUE,
legend=c( paste0('test(',nrow(DF.test),')'),paste0('train(',nrow(DF.train),')') ),
args.legend=list( x="center" ) )
# ✔ NORMALICE -----------------------------------------------------------------
DFn.test <- DF.test
DFn.train <- DF.train
DFn.test[, -Cidx] <- sapply( DFn.test[, -Cidx] , my.math.normalize )
DFn.train[,-Cidx] <- sapply( DFn.train[,-Cidx] , my.math.normalize )
DF %>% head(3)
# ✔ REDUCCION DE VARIABLES -------------------------------------------------------
# ...saltearemos esta sección hasta encontrar evidencia de posibilidad de reducir variables...
# ✔ OUTLIERS --------------------------------------------------------------
# ...saltearemos esta sección hasta encontrar evidencia de valores outliers...
• INSPECCION
De la siguientes gráficas se observa sencillamente que las variables dimensionales del Pétalo de la flor constituyen medidas significativas de diferenciación, especialmente para la clase Setosa.
par(mfcol = c(2,2) )
barplot( Sepal.Length~Species , aggregate( Sepal.Length~Species, iris, sum), col=c(6,7,8))
barplot( Sepal.Width ~Species , aggregate( Sepal.Width ~Species, iris, sum), col=c(6,7,8))
barplot( Petal.Length~Species , aggregate( Petal.Length~Species, iris, sum), col=c(6,7,8))
barplot( Petal.Width ~Species , aggregate( Petal.Width ~Species, iris, sum), col=c(6,7,8))
par(mfcol = c(1,1) )
• CORRELACIONES
• HISTOGRAMAS
• QQ-plot
• DENSIDAD x GRUPOS
El siguiente paso es realizar pruebas sobre los atributos de los datos, para verificar si se confirman los supuestos de las técnicas de clasificación.
- CORRELACIONES que se incluyen en la matriz de gráficas de dispersión de los pares de variables. Además también, en la misma se incluyen histogramas para cada variable (sin sub-clasificacion).
- HISTOGRAMAS para cada variable seccionadas por las clases de predicción, y en la misma se incluye, curva de densidad y curva normal. En la mayoría de los algoritmos se presupone distribución normal.
- QQ-PLOTs para cada variable seccionadas por las clases de predicción, para observar el ajuste con los intervalos de la curva normal. Es decir, que se confirmariía normal si la recta de tendencia atraviesa diagonalmente los cuadrantes y la media se ubica en el centro.
Observaciones:
- Los histogramas de cada variable muestran que algunos subconjuntos se concentran mientras que otros se dispersan.
- En las variables cuyos subconjuntos se separan habrá mayor capacidad de clasificación.
- En las variables cuyos subconjuntos se separan se obtendrá una prueba de normalidad univariante insatisfactoria, entonces habrá que segmentarlas por clase para constatar que cada clase se distribuye normal.
- Los histogramas segmentados por clases se observan normales, sin embargo, esto no es así para la variable Petal.Width, lo cual se confirma en el gráfico QQ-plot.
- Hay una fuerte correlación con la variable de clasificación pero también se observa fuerte correlación entre las variables individuales indicando que no son independientes.
- Finalmente las gráficas de densidad por grupos permiten observar la superposición de los datos.
# ==============================================================================
# ==== PRUEBAS de ATRIBUTOS (propiedades intrínsecas de las variables)
# ==============================================================================
# ✔ CORRELACIONES --------------------------------------------------------------
ploto_pairs( DF, DF[,Cidx] )
# ✔ HISTOGRAMAS --------------------------------------------------------------
histo_multi( DF[,-Cidx], DF[,Cidx] )
# ✔ QQ-plot --------------------------------------------------------------
qqplot_multi( DF[,-Cidx], DF[,Cidx] )
# ✔ DENSIDAD x GRUPOS --------------------------------------------------------------
plot1 <- ggplot( DF, aes_string( x=names(DF)[1], colour=Cname ) ) + geom_density() + theme_bw()## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.plot2 <- ggplot( DF, aes_string( x=names(DF)[2], colour=Cname ) ) + geom_density() + theme_bw()
plot3 <- ggplot( DF, aes_string( x=names(DF)[3], colour=Cname ) ) + geom_density() + theme_bw()
plot4 <- ggplot( DF, aes_string( x=names(DF)[4], colour=Cname ) ) + geom_density() + theme_bw()
ggpubr::ggarrange( plot1, plot2, plot3, plot4, common.legend = TRUE, legend = "bottom")
• ANOVA
En el ANOVA (ANalysis Of VAriance) o AOV, el estadístico F, sigue la distribución F bajo la hipótesis nula sobre el modelo, de que ninguna de las variables explicativas influye sobre la observada Y.
El estadístico \(p\) para cada variable explicativa evalúa la hipótesis nula de que dicha variable sea NO Explicativa y su coeficiente sea igual a cero (no hay efecto). Un valor \(p < 0,05\) (bajo) indica que se puede rechazar la hipótesis nula y confirmar que la variable Es explicativa porque los cambios en dicha variable se relacionan con los cambios de la variable observada.
A continuación puede observarse que todos los valores \(p\) se encuentran por debajo de 0.05, lo cual confirma que todas son variables explicativas, sin embargo, se observa puntualmente Sepal.Length como una variable de Mean Squared Error considerablemente mas elevado que los demás. Esto podría indicar menor grado de Explicabilidad o que la variable se comporta en relación no-lineal.
aov( as.integer(Species) ~ . , data=iris )
.Last.value %>% summary## Call:
## aov(formula = as.integer(Species) ~ ., data = iris)
##
## Terms:
## Sepal.Length Sepal.Width Petal.Length Petal.Width Residuals
## Sum of Squares 61.24021 11.35617 18.44591 1.99710 6.96061
## Deg. of Freedom 1 1 1 1 145
##
## Residual standard error: 0.2190986
## Estimated effects may be unbalanced
## Length Class Mode
## help_type 0 -none- NULL• NORMALIDAD univariante
A continuación ejecutas las pruebas de normalidad univariante.
Se observan resultados pobres en las variables mas dispersas como fue previsto anteriormente.
# ✔ NORMALIDAD univariante --------------------------------------------------------------
# ••••••••••• Shapiro-Wilks
sapply( DF[,-Cidx] , function(x) round( stats::shapiro.test(x)$p.value ,4) )
# ••••••••••• Anderson-Darling
sapply( DF[,-Cidx] , function(x) round( nortest::ad.test(x)$p.value ,4) )
# ••••••••••• Lilliefors (Kolmogorov-Smirnov)
sapply( DF[,-Cidx] , function(x) round( nortest::lillie.test(x)$p.value ,4) )## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.0102 0.1012 0.0000 0.0000
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.0225 0.0202 0.0000 0.0000
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## 0.0058 0.0003 0.0000 0.0000• NORMALIDAD univariante X clase
Continuamos con las pruebas univariantes fraccionadas por la variable clasificadora.
Se observan mejores resultados de los p-valores. Sin embargo, confirmamos lo predicho para la variable Petal.Width. Su distribución no es normal, o no hay prueba evidente de ello.
# ✔ NORMALIDAD univariante X grupo -------------------------------------------------------
# ••••••••••• Shapiro-Wilks
mapply( function(col) tapply( col, DF[,Cidx],
function(x) round( stats::shapiro.test(x)$p.value ,4) ) , DF[,-Cidx] )
# ••••••••••• Anderson-Darling
aggregate( DF[,-Cidx], by=list( DF[,Cidx] ),
function(x) round( nortest::ad.test(x)$p.value ,4) ) %>% print.data.frame()
# ••••••••••• Lilliefors (Kolmogorov-Smirnov)
aggregate( DF[,-Cidx], by=list( DF[,Cidx] ),
function(x) round( nortest::lillie.test(x)$p.value ,4) ) %>% print.data.frame()
# ••••••••••• ...otros...
aggregate( DF[,-Cidx], by=list( DF[,Cidx] ),
function(x) round( nortest::lillie.test(x)$p.value ,4) ) %>% print.data.frame()
aggregate( DF[,-Cidx], by=list( DF[,Cidx] ),
function(x) round( nortest::cvm.test(x)$p.value ,4) ) %>% print.data.frame()
aggregate( DF[,-Cidx], by=list( DF[,Cidx] ),
function(x) round( nortest::pearson.test(x)$p.value,4) ) %>% print.data.frame()
aggregate( DF[,-Cidx], by=list( DF[,Cidx] ),
function(x) round( nortest::sf.test(x)$p.value ,4) ) %>% print.data.frame()## Sepal.Length Sepal.Width Petal.Length Petal.Width
## setosa 0.4595 0.2715 0.0548 0.0000
## versicolor 0.4647 0.3380 0.1585 0.0273
## virginica 0.2583 0.1809 0.1098 0.0870
## Group.1 Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 setosa 0.3352 0.2102 0.0108 0.0000
## 2 versicolor 0.4333 0.1406 0.1446 0.0144
## 3 virginica 0.1475 0.1018 0.1074 0.0508
## Group.1 Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 setosa 0.0969 0.1863 0.0049 0.0000
## 2 versicolor 0.2942 0.0663 0.0839 0.0081
## 3 virginica 0.0959 0.0400 0.1100 0.0658
## Group.1 Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 setosa 0.0969 0.1863 0.0049 0.0000
## 2 versicolor 0.2942 0.0663 0.0839 0.0081
## 3 virginica 0.0959 0.0400 0.1100 0.0658
## Group.1 Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 setosa 0.2597 0.2324 0.0069 0.0000
## 2 versicolor 0.4039 0.0980 0.1498 0.0213
## 3 virginica 0.1522 0.0850 0.1674 0.0618
## Group.1 Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 setosa 0.2386 0.0512 0.000 0.0000
## 2 versicolor 0.3692 0.1670 0.029 0.0000
## 3 virginica 0.1006 0.0065 0.450 0.0088
## Group.1 Sepal.Length Sepal.Width Petal.Length Petal.Width
## 1 setosa 0.5357 0.1181 0.0324 0.0000
## 2 versicolor 0.6153 0.3462 0.1643 0.0448
## 3 virginica 0.2210 0.1205 0.1301 0.1354• NORMALIDAD multivariante
Seguimos con las pruebas multivariantes de normalidad, las cuales fracasan obteniendo bajos niveles de p-valores en todos los casos (excepto en la clase Sepal.Width).
# ✔ NORMALIDAD multivariante --------------------------------------------------------------
# ••••••••••• Henze-Zirkler
mvnormalTest::mhz( DF[,-Cidx] )
# ••••••••••• Royston
mvnTest::R.test( DF[,-Cidx], qqplot=FALSE )
# ••••••••••• Mardia (Skewness and Kurtosis)
mvnormalTest::mardia( DF[,-Cidx] )$mv.test %>% print.data.frame## $mv.test
## Statistic p-value Result
## 2.3364 0 NO
##
## $uv.shapiro
## W p-value UV.Normality
## Sepal.Length 0.9761 0.0102 No
## Sepal.Width 0.9849 0.1012 Yes
## Petal.Length 0.8763 0 No
## Petal.Width 0.9018 0 No
##
## Royston test for Multivariate Normality
##
## data : DF[, -Cidx]
##
## R : 50.3738
## p-value : 3.133482e-11
##
## Result : Data are not multivariate normal (sig.level = 0.05)
## Test Statistic p-value Result
## 1 Skewness 67.4305 0 NO
## 2 Kurtosis 0.0509 0.9594 YES
## 3 MV Normality <NA> <NA> NO• HOMOGENEIDAD
Las pruebas de homogeneidad de varianzas obtienen p-valores menores que 0.05, lo cual rechaza la hipotesis nula de que las varianzas de los grupos generados por la variable de clasificación son iguales. Efectivamente, observamos en las graficas anteriores que la variabilidad es diferente.
Algo a tener en cuenta es que el modelo QDA (Quadratic discriminant) no asume que exista homogeneidad en la matriz de covarianza.
# ✔ HOMOGENEIDAD --------------------------------------------------------------
# ••••••••••• Bartlett
bartlett.test( DF[,-Cidx], g=DF[,Cidx] )$p.value
# ••••••••••• Box's M-test
DF[,-Cidx] %>% cov # matriz de covarianza comun
lapply( split( DF[,-Cidx], DF[,Cidx] ) , cov ) # matriz de covarianza por clases
biotools::boxM( DF[,-Cidx], DF[,Cidx])$p.value## [1] 1.795553e-63
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.6856935 -0.0424340 1.2743154 0.5162707
## Sepal.Width -0.0424340 0.1899794 -0.3296564 -0.1216394
## Petal.Length 1.2743154 -0.3296564 3.1162779 1.2956094
## Petal.Width 0.5162707 -0.1216394 1.2956094 0.5810063
## $setosa
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.12424898 0.099216327 0.016355102 0.010330612
## Sepal.Width 0.09921633 0.143689796 0.011697959 0.009297959
## Petal.Length 0.01635510 0.011697959 0.030159184 0.006069388
## Petal.Width 0.01033061 0.009297959 0.006069388 0.011106122
##
## $versicolor
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.26643265 0.08518367 0.18289796 0.05577959
## Sepal.Width 0.08518367 0.09846939 0.08265306 0.04120408
## Petal.Length 0.18289796 0.08265306 0.22081633 0.07310204
## Petal.Width 0.05577959 0.04120408 0.07310204 0.03910612
##
## $virginica
## Sepal.Length Sepal.Width Petal.Length Petal.Width
## Sepal.Length 0.40434286 0.09376327 0.30328980 0.04909388
## Sepal.Width 0.09376327 0.10400408 0.07137959 0.04762857
## Petal.Length 0.30328980 0.07137959 0.30458776 0.04882449
## Petal.Width 0.04909388 0.04762857 0.04882449 0.07543265
##
## [1] 3.352034e-20• ALGORITMO y PREDICCION
Finalmente, aplicamos los algoritmos y realizamos las predicciones sobre el conjunto reservado para pruebas.
El arbol de decisión indica que es capaz de clasificar el 100% de los datos de setosa con la variable Petal.Length<2.5, mientras que no le resulta tan determinista para los otros dos conjuntos de clases.
# ==============================================================================
# ==== ALGORITMOs
# ==============================================================================
# ✔ CLASIFICACION --------------------------------------------------------------
# ••••••••••• AD (Arbol de Decision)
model_ad <- rpart::rpart( fmla , data=DF, parms=list(split="information") )
model_ad %>% rpart.plot::rpart.plot( type=1, extra=100, cex = 1.5 )
# ••••••••••• LDA (Discriminante Lineal)
model_lda <- MASS::lda( fmla, data=DF )
# ••••••••••• NBayesiano (Naive-Bayes)
model_nb <- naivebayes::naive_bayes( fmla , data=DF)
# ••••••••••• QDA (Quadratic discriminant)
model_qda <- MASS::qda( fmla , data=DF)
# ••••••••••• k-NN (k-nearest neighbors) --> https://www.youtube.com/watch?v=ENSHwuJU5sU
# ... k-NN no genera un modelo ...
# ••••••••••• NN (Neural Network)
#set.seed(123);
model_NN <- neuralnet::neuralnet( Species ~ . , DF,
hidden=c(8,6,6), rep=3, stepmax=10000,
#algorithm='backprop', learningrate=0.02,
err.fct='sse' , threshold=0.01,
act.fct='logistic' , linear.output=FALSE )
model_NN %>% plot(rep="best")
# ✔ PREDICCIONES --------------------------------------------------------------
# predict() llama a predict.xyz() y ejecuta el modelo sobre un conjunto de datos
pred_ad_class <- predict( model_ad , DF.test[-Cidx], type="class" )
pred_lda <- predict( model_lda , DF.test[-Cidx])
pred_nb_class <- predict( model_nb , DF.test[-Cidx])
pred_qda <- predict( model_qda, DF.test[-Cidx])
pred_knn <- list()
pred_knn[[1]] <- class::knn( DF[,-Cidx], DF.test[,-Cidx], DF[,Cidx], l = 0, k = 3)
pred_knn[[2]] <- class::knn( DF[,-Cidx], DF.test[,-Cidx], DF[,Cidx], l = 0, k = 5)
pred_knn[[3]] <- class::knn( DF[,-Cidx], DF.test[,-Cidx], DF[,Cidx], l = 0, k = 7)
pred_knn[[4]] <- class::knn( DF[,-Cidx], DF.test[,-Cidx], DF[,Cidx], l = 0, k = 9)
pred_NN <- predict( model_NN , DF.test[,-Cidx] )
pred_NN_class <- factor(apply(pred_NN, 1, which.max ), labels=model_NN$model.list$response )• CONFUSION
Las confusiones pueden representarse en tabla, o en un gráfico de dispersión.
En cada ejecución, el algoritmo de particionamiento del dataset puede generar distintos subconjuntos dependiendo del valor que obtenga la semilla aleatoria. Sin embargo, el algoritmo que produce mejores resultados consistentemente es el LDA.
# ==============================================================================
# ==== METRICAS de RENDIMIENTO y CALIDAD
# ==============================================================================
# ✔ CONFUSION --------------------------------------------------------------
table( DF.test[,Cidx], `AD` =pred_ad_class ) # Confusion: 6 items para test=train
table( DF.test[,Cidx], `LDA` =pred_lda$class ) # Confusion: 3 items
table( DF.test[,Cidx], `NB` =pred_nb_class ) # Confusion: 6 items
table( DF.test[,Cidx], `QDA` =pred_qda$class ) # Confusion: 3 items
table( DF.test[,Cidx], `k-NN(3)`=pred_knn[[1]] ) # Confusion: 6 items
table( DF.test[,Cidx], `k-NN(5)`=pred_knn[[2]] ) # Confusion: 5 items
table( DF.test[,Cidx], `k-NN(7)`=pred_knn[[3]] ) # Confusion: 4 items
table( DF.test[,Cidx], `k-NN(9)`=pred_knn[[4]] ) # Confusion: 3 items
table( DF.test[,Cidx], `NN` =pred_NN_class ) # Confusion: 1 items## AD
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 49 1
## virginica 0 5 45
## LDA
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 48 2
## virginica 0 1 49
## NB
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 47 3
## virginica 0 3 47
## QDA
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 48 2
## virginica 0 1 49
## k-NN(3)
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 47 3
## virginica 0 3 47
## k-NN(5)
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 47 3
## virginica 0 2 48
## k-NN(7)
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 47 3
## virginica 0 1 49
## k-NN(9)
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 48 2
## virginica 0 1 49
## NN
## setosa versicolor virginica
## setosa 50 0 0
## versicolor 0 49 1
## virginica 0 0 50# Visualizacion de Confusiones
conf_count <-
c( AD = sum( DF.test[,Cidx] != pred_ad_class ) ,
LDA = sum( DF.test[,Cidx] != pred_lda$class ),
NB = sum( DF.test[,Cidx] != pred_nb_class ) ,
QDA = sum( DF.test[,Cidx] != pred_qda$class ),
kNN3 = sum( DF.test[,Cidx] != pred_knn[[1]] ),
kNN5 = sum( DF.test[,Cidx] != pred_knn[[2]] ),
kNN7 = sum( DF.test[,Cidx] != pred_knn[[3]] ),
kNN9 = sum( DF.test[,Cidx] != pred_knn[[4]] ),
NN = sum( DF.test[,Cidx] != pred_NN_class )
)
barplot( conf_count, legend=row.names(conf_count), col=terrain.colors(length(conf_count))[conf_count] )
# Visualizacion de Confusion
plot1 <- ggplot2::ggplot( DF.test, aes_string( Cname, pred_ad_class , color=Cname) ) +
geom_jitter(width=0.2, height=0.2) + labs( title='Arbol de Decision' )
plot2 <- ggplot2::ggplot( DF.test, aes_string( Cname, pred_lda$class, color=Cname) ) +
geom_jitter(width=0.2, height=0.2) + labs( title='LDA (An. Discriminante Lineal)' )
plot3 <- ggplot2::ggplot( DF.test, aes_string( Cname, pred_nb_class , color=Cname) ) +
geom_jitter(width=0.2, height=0.2) + labs( title='Naive-Bayes' )
plot4 <- ggplot2::ggplot( DF.test, aes_string( Cname, pred_qda$class , color=Cname) ) +
geom_jitter(width=0.2, height=0.2) + labs( title='QDA (Quadratic discriminant)' )
ggpubr::ggarrange( plot1, plot2, plot3, plot4, common.legend = TRUE, legend = "bottom")
plot1 <- ggplot2::ggplot( DF.test, aes_string( Cname, pred_knn[[1]] , color=Cname) ) +
geom_jitter(width=0.2, height=0.2) + labs( title='k-NN (k-Nearest Neighbors) K=3' )
plot2 <- ggplot2::ggplot( DF.test, aes_string( Cname, pred_knn[[2]], color=Cname) ) +
geom_jitter(width=0.2, height=0.2) + labs( title='k-NN (k-Nearest Neighbors) K=5' )
plot3 <- ggplot2::ggplot( DF.test, aes_string( Cname, pred_knn[[3]] , color=Cname) ) +
geom_jitter(width=0.2, height=0.2) + labs( title='k-NN (k-Nearest Neighbors) K=7' )
plot4 <- ggplot2::ggplot( DF.test, aes_string( Cname, pred_knn[[4]] , color=Cname) ) +
geom_jitter(width=0.2, height=0.2) + labs( title='k-NN (k-Nearest Neighbors) K=8' )
ggpubr::ggarrange( plot1, plot2, plot3, plot4, common.legend = TRUE, legend = "bottom")
ggplot2::ggplot( DF.test, aes_string( Cname, pred_NN_class , color=Cname) ) +
geom_jitter(width=0.2, height=0.2) + labs( title='NN (Neural Networks)' )
Los datos que los clasificadores no pudieron procesar adecuadamente son los siguientes.
# Crear una lista con los modelos y sus respectivas predicciones
models <- list(
"AD" = pred_ad_class,
"LDA" = pred_lda$class,
"QDA" = pred_qda$class,
"NB" = pred_nb_class,
"kNN(3)" = pred_knn[[1]],
"kNN(5)" = pred_knn[[2]],
"kNN(7)" = pred_knn[[3]],
"kNN(9)" = pred_knn[[4]],
"NN" = pred_NN_class
)
# Inicializar un data frame vacío para almacenar los errores
errors <- data.frame()
# Iterar sobre cada modelo y calcular los errores
for (model_name in names(models)) {
pred <- models[[model_name]]
incorrect_rows <- DF.test[which(DF.test[, Cidx] != pred), ]
# Si hay errores, añadirlos; si no, añadir una fila vacía
if (nrow(incorrect_rows) > 0) {
model_errors <- cbind(model = model_name, incorrect_rows)
errors <- rbind(errors, model_errors)
}
}
# Remover duplicados si es necesario
errors <- errors[!duplicated(errors), ]
# Resultado
errorsPor ultimo calculamos un conjunto de métricas en los cuales el LDA obtiene las calificaciones mas altas para todas ellas. Las métricas se calculan con la función confusionMatrix del paquete caret para el resultado final, y separadas por clase. Sin embargo, es posible calcularlas sin convocar módulos externos.
# ✔ CONFUSION metrics -------------------------------------------------
cm_ad = caret::confusionMatrix( data=pred_ad_class , reference=DF.test[,Cidx])
cm_lda = caret::confusionMatrix( data=pred_lda$class, reference=DF.test[,Cidx])
cm_nb = caret::confusionMatrix( data=pred_nb_class , reference=DF.test[,Cidx])
cm_qda = caret::confusionMatrix( data=pred_qda$class, reference=DF.test[,Cidx])
cm_knn3 = caret::confusionMatrix( data=pred_knn[[1]] , reference=DF.test[,Cidx])
cm_knn5 = caret::confusionMatrix( data=pred_knn[[2]] , reference=DF.test[,Cidx])
cm_knn7 = caret::confusionMatrix( data=pred_knn[[3]] , reference=DF.test[,Cidx])
cm_knn9 = caret::confusionMatrix( data=pred_knn[[4]] , reference=DF.test[,Cidx])
cm_NN = caret::confusionMatrix( data=pred_NN_class , reference=DF.test[,Cidx])# ••••••••••• overall
result <- cbind(
c( model='AD' , cm_ad$overall %>% round(4) ),
c( model='LDA', cm_lda$overall %>% round(4) ),
c( model='NB' , cm_nb$overall %>% round(4) ),
c( model='QDA', cm_qda$overall %>% round(4) ),
c( model='kNN(3)', cm_knn3$overall %>% round(4) ),
c( model='kNN(5)', cm_knn5$overall %>% round(4) ),
c( model='kNN(7)', cm_knn7$overall %>% round(4) ),
c( model='kNN(9)', cm_knn9$overall %>% round(4) ),
c( model='NN', cm_NN$overall %>% round(4) )
)
result_df <- as.data.frame(result)
names(result_df) <- result_df[1,]
result_df <- result_df[-1,]
# ••••••••••• byClass
confMatrixByClass <-
lapply( 1:nrow(cm_ad$byClass), function(k) { # iterate over classes results
cbind(
c( cm_ad$byClass[k,] %>% round(4) ),
c( cm_lda$byClass[k,] %>% round(4) ),
c( cm_nb$byClass[k,] %>% round(4)),
c( cm_qda$byClass[k,] %>% round(4) ),
c( cm_knn3$byClass[k,] %>% round(4) ),
c( cm_knn5$byClass[k,] %>% round(4) ),
c( cm_knn7$byClass[k,] %>% round(4) ),
c( cm_knn9$byClass[k,] %>% round(4) ),
c( cm_NN$byClass[k,] %>% round(4) )
) %>% as.data.frame() %>%
setNames( c( 'AD','LDA','NB','QDA',
'kNN(3)','kNN(5)','kNN(7)','kNN(9)',
'NN' ) )
})
# get class names in the correct order
classnames = str_split( row.names(cm_ad$byClass), ": ", simplify=TRUE)[,2]
# prepare to merge 3 dataframes
for( i in 1:length(confMatrixByClass) ){
confMatrixByClass[[i]] <- add_column( confMatrixByClass[[i]],
metric=row.names( confMatrixByClass[[i]] ) , .before = 1 )
confMatrixByClass[[i]] <- add_column( confMatrixByClass[[i]],
class=classnames[i] , .before = 1 )
row.names(confMatrixByClass[[i]]) <- NULL
}
confMatrixByClass_merged = merge( confMatrixByClass[[1]], confMatrixByClass[[2]], all=TRUE )
confMatrixByClass_merged = merge( confMatrixByClass_merged, confMatrixByClass[[3]], all=TRUE )
# print merged dataframe of 'caret' metrics
confMatrixByClass_merged %>% arrange( metric, class )
accuracies <- confMatrixByClass_merged[ confMatrixByClass_merged$metric=="Balanced Accuracy",] %>% select(-c(2))
row.names(accuracies) <- accuracies[,1]
accuracies <- accuracies[,2:ncol(accuracies)] %>% as.matrix
accuracies %>% barplot( beside=TRUE, main="Accuracy",
col=c(2,3,4), border=0, legend=row.names(accuracies),
args.legend=list( x=ncol(accuracies)*4, y=0.8 ) )
• CONCLUSIONES
El algoritmo de modelado que produjo menos confusiones fue el NN, sin embargo, LDA, QDA y kNN(9) resultaron muy competitivos.
Se debe considerar que el LDA y el QDA, produjeron resultados con mayor facilidad y menor uso del cómputo, mientras que la red neural NN podría utilizarse en un escenario de procesamiento contínuo.
Las confusiones de LDA y QDA pueden deberse particularmente a la no-lineal observada en Petal.Width y la alta correlación observada entre ella y Petal.Length que resulta evidente ya que refieren al crecimiento de un mismo aspecto, el Pétalo.
Sin embargo, las pruebas realizadas indican que los clasificadores empeoran su resultado al remover una de esas variables.
En todo los modelos aplicados, teniendo en cuenta las 150 muestras, se obtuvo eficiencias de clasificación superiores al 95%, obteniendo 97.5% para el LDA.
Aquellos errores o confusiones encontrados se observan particularmente en las zonas de unión entre grupos debido a muy poca evidencia clasificadora en este punto. Por ese motivo, sería interesante observar el origen de los datos e identificar si no pertenecen a un rango de valores outliers.
Como propuesta de mejora, se podría estudiar la posibilidad de utilizar relaciones entre dichas variables como relación de aspecto Largo/Ancho del Petalo y del Tallo, o de Largo de Pétalo vs Largo de Tallo. Es necesario identificar una relación que sea representativa de la edad de la flor, ya que todas las medidas están relacionadas con ello.
✔ Algoritmos explicados con animaciones
k-NN (k-Nearest Neighbors)
La idea general es hallar los K vecinos más cercanos y calcular la clase más frecuente para asignarla al nuevo elemento.
