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)
Diego Jacobi, Argentina
18 octubre, 2022
Abstract
El objetivo de este pagina es un resumen personal de pasos a seguir en el proceso de clasificaciĆ³n. No tiene por objetivo explicar la matemĆ”tica interna de cada operaciĆ³n.<br> Incluye conceptos aprendidos en la asignatura āMĆ©todos predictivosā, dictada por el Dr.Ā NoĆ© Amir, de la Universidad TecnolĆ³gica En Linea (UTEL) en el programa de āMaestrĆa en Ciencia de Datos para Negociosā.
ā 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()
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
# 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)' )## 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
Los datos que los clasificadores no pudieron procesar adecuadamente son los siguientes.
# erors: data that is not correctly predicted
errors <-
rbind(
cbind( model="AD", DF.test[ which( DF.test[,Cidx] != pred_ad_class ), ]),
cbind( model="LDA", DF.test[ which( DF.test[,Cidx] != pred_lda$class ), ]),
cbind( model="QDA", DF.test[ which( DF.test[,Cidx] != pred_qda$class ), ]),
cbind( model="NB", DF.test[ which( DF.test[,Cidx] != pred_nb_class ), ]),
cbind( model="kNN(3)", DF.test[ which( DF.test[,Cidx] != pred_knn[[1]] ), ]),
cbind( model="kNN(5)", DF.test[ which( DF.test[,Cidx] != pred_knn[[2]] ), ]),
cbind( model="kNN(7)", DF.test[ which( DF.test[,Cidx] != pred_knn[[3]] ), ]),
cbind( model="kNN(8)", DF.test[ which( DF.test[,Cidx] != pred_knn[[4]] ), ]),
cbind( model="NN", DF.test[ which( DF.test[,Cidx] != pred_NN_class ), ])
) ;errorserrors[ duplicated( errors[-1] ), ]Por 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
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) )
) %>% as.data.frame() %>% setNames( .[1,] ) %>% .[-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.
