Carga de librerias

library(neuralnet) # Para el algoritmo redes neuronales
library(randomForest) # Para el algoritmo Bosque Aleatorio
## randomForest 4.7-1.1
## Type rfNews() to see new features/changes/bug fixes.
library(readxl) 
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:randomForest':
## 
##     combine
## The following object is masked from 'package:neuralnet':
## 
##     compute
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggplot2)
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:randomForest':
## 
##     margin
library(corrplot)
## corrplot 0.92 loaded
library(ggcorrplot)

Lectura de datos

BreastFinal <- read_excel("E:/U-MAQUINAS DE APRENDIZAJE/EXAMEN FINAL/BreastFinal.xlsx")
BreastFinal$Class<-as.factor(BreastFinal$Class)

Variable target: clases de tejido
- car: carcinoma
- fad: fibroadenoma
- mas: Mastopatía
- gla: Glandular
- con: Conectivo
- adi: Adiposo

Variables predictoras: - Impedancia (ohm) a frecuencia cero (I0)
- Ángulo de fase a 500 KHz (PA500)
- Pendiente de alta frecuencia del ángulo de fase (HFS)
- Distancia de impedancia entre extremos espectrales (DA)
- Área bajo el espectro (AREA)
- Área normalizada por DA (A.DA)
- Máximo del espectro (Max.IP)
- Distancia entre I0 y parte real del punto de máxima frecuencia (DR)
- Longitud de la curva espectral (P)

Renombramos las siguientes 2 columnas para no usar caracteres especiales:

BreastFinal<-BreastFinal %>% rename (A.DA = `A/DA`, Max.IP = `Max IP`)
head(BreastFinal)
dim(BreastFinal)
## [1] 106  11
summary(BreastFinal)
##      Case #       Class          I0             PA500        
##  Min.   :  1.00   adi:22   Min.   : 103.0   Min.   :0.01239  
##  1st Qu.: 27.25   car:21   1st Qu.: 250.0   1st Qu.:0.06741  
##  Median : 53.50   con:14   Median : 384.9   Median :0.10542  
##  Mean   : 53.50   fad:15   Mean   : 784.3   Mean   :0.12013  
##  3rd Qu.: 79.75   gla:16   3rd Qu.:1488.0   3rd Qu.:0.16960  
##  Max.   :106.00   mas:18   Max.   :2800.0   Max.   :0.35832  
##       HFS                 DA               Area                A.DA        
##  Min.   :-0.06632   Min.   :  19.65   Min.   :    70.43   Min.   :  1.596  
##  1st Qu.: 0.04398   1st Qu.:  53.85   1st Qu.:   409.65   1st Qu.:  8.180  
##  Median : 0.08657   Median : 120.78   Median :  2219.58   Median : 16.134  
##  Mean   : 0.11469   Mean   : 190.57   Mean   :  7335.16   Mean   : 23.474  
##  3rd Qu.: 0.16650   3rd Qu.: 255.33   3rd Qu.:  7615.20   3rd Qu.: 30.953  
##  Max.   : 0.46775   Max.   :1063.44   Max.   :174480.48   Max.   :164.072  
##      Max.IP              DR                P         
##  Min.   :  7.969   Min.   : -9.258   Min.   : 125.0  
##  1st Qu.: 26.894   1st Qu.: 41.781   1st Qu.: 270.2  
##  Median : 44.216   Median : 97.833   Median : 454.1  
##  Mean   : 75.381   Mean   :166.711   Mean   : 810.6  
##  3rd Qu.: 83.672   3rd Qu.:232.990   3rd Qu.:1301.6  
##  Max.   :436.100   Max.   :977.552   Max.   :2896.6

I0: El valor medio de esta variable es mayor que su mediana. Esto indica una distribución asimétrica a la derecha. La distribución de la variable tiene una diferencia bastante grande entre sus valores mínimo y máximo.
PA500: La mediana de esta variable es inferior a su valor medio. Esto indica una distribución asimétrica a la izquierda. La distribución de la variable tiene una diferencia bastante grande entre sus valores mínimo y máximo.
HFS: La mediana de esta variable es ligeramente inferior a su valor medio. Esto indica una distribución ligeramente sesgada hacia la izquierda. La distribución de la variable tiene una menor diferencia entre sus valores mínimo y máximo.
DA: El valor medio de esta variable es mayor que su mediana. Esto indica una distribución asimétrica a la derecha. La distribución de la variable tiene una diferencia bastante grande entre sus valores mínimo y máximo.
Área: El valor medio de esta variable es mucho mayor que su mediana. Esto indica una distribución asimétrica a la derecha. La distribución de la variable tiene una diferencia bastante grande entre sus valores mínimo y máximo.
A.DA: El valor medio de esta variable es mayor que su mediana. Esto indica una distribución asimétrica a la derecha. La distribución de la variable tiene una diferencia bastante grande entre sus valores mínimo y máximo.
Max.IP: La mediana de esta variable es ligeramente inferior a su valor medio. Esto indica una distribución ligeramente sesgada hacia la izquierda. La distribución de la variable tiene una diferencia bastante grande entre sus valores mínimo y máximo.
DR: La mediana de esta variable es menor que su valor medio. Esto indica una distribución asimétrica a la izquierda. La distribución de la variable tiene una diferencia bastante grande entre sus valores mínimo y máximo.
P: El valor medio de esta variable es mucho mayor que su mediana. Esto indica una distribución asimétrica a la derecha. La distribución de la variable tiene una diferencia bastante grande entre sus valores mínimo y máximo.

Valores perdidos

sum(is.na(BreastFinal))
## [1] 0

No existen valores perdidos en la data.

Analizando mediante visualizacion

Analisis univariado

Grafico de barras
Tabla_class <- BreastFinal%>% group_by(Class) %>% summarise(Total=n()) %>% mutate(Porcentaje = round(Total/sum(Total)*100, 1))

ggplot(Tabla_class, aes(x = Class, y=Total,fill=Class) ) +          
  geom_bar(width = 0.9, stat="identity", position = position_dodge())+
  ylim(c(0,100))+
  labs(x="Clase de tejido", y= "Frecuencia \n (Porcentajes)") + 
  labs(fill = "")+
  geom_text(aes(label=paste0(Total," ", "", "(", Porcentaje, "%",")")), vjust=-0.9, color="black", hjust=0.5,position = position_dodge(0.9),  angle=0,  size=4.0) + 
  scale_fill_discrete(name = "Clase", labels = c("adi", "car", "con", "fat", "gla", "mas")) + 
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  facet_wrap(~"Variable Class")

Se observa que la data esta balanceada, ya que los porcentajes de cada clase son superiores al 10%.

Analisis Bivariado

Graficos de dispersion
ggplot(BreastFinal, aes(x = DA, y = I0)) +
  geom_point(aes(color = factor(Class))) +
  labs(x = "Distancia de impedancia entre extremos espectrales \n DA",
       y = "Impedancia (ohm) a frecuencia cero \n I0",
       title = "Grafico de dispersion")

Se observa que en la dispersión que las variables DA y I0, están separando bien a las clases diferenciados por los colores.

ggplot(BreastFinal, aes(x=PA500, y=HFS)) +
  geom_point(aes(color = Class)) +
  geom_smooth(method = "lm", col = "#C42126", se = FALSE, size = 1)+
  labs(x = "Ángulo de fase a 500 KHz \n PA500",
       y = "Pendiente de alta frecuencia del ángulo de fase \n HFS",
       title = "Grafico de dispersion") 
## `geom_smooth()` using formula = 'y ~ x'

Se observa que en la dispersión que las variables PA500 y HFS se pueden ajustar a un modelo de regresión lineal.

ggplot(BreastFinal, aes(x = DA, y = A.DA)) +
  geom_point(aes(color = Class)) +
  geom_smooth(method = "lm", col = "#C42126", se = FALSE, size = 1) +
  labs(x = "Distancia de impedancia entre extremos espectrales",
       y = "Área normalizada por DA ",
       title = "Grafico de dispersion")
## `geom_smooth()` using formula = 'y ~ x'

Se observa que en la dispersión que las variables DA y A.DA, están separando bien a las clases diferenciados por los colores.

Diagrama de Cajas
par(mfrow= c(2,5))
for  ( i in  3 : ncol ( BreastFinal ) )  {
   boxplot ( BreastFinal [ , i ] , main =  names(BreastFinal) [ i ] ,
           col =  rainbow(ncol(BreastFinal)-2),
           notch = T) 
}

par(mfrow= c(3,3))

boxplot(formula = I0 ~ Class, data =  BreastFinal)
boxplot(formula = PA500 ~ Class, data =  BreastFinal)
boxplot(formula = HFS ~ Class, data =  BreastFinal)
boxplot(formula = DA ~ Class, data =  BreastFinal)
boxplot(formula = Area ~ Class, data =  BreastFinal)
boxplot(formula = A.DA ~ Class, data =  BreastFinal)
boxplot(formula = Max.IP ~ Class, data =  BreastFinal)
boxplot(formula = DR ~ Class, data =  BreastFinal)
boxplot(formula = P ~ Class, data =  BreastFinal)

Cuando se analiza el box plot de cada una de las variables, se observa que existe una diferencia significativa entre las escalas. Por esta razón, el conjunto de datos debe estandarizarse. También por el elevado número de observaciones atípicas que se observan en algunas variables.

Grafico de Area
ggplot(BreastFinal, aes(x=DR, y=P, fill=Class)) +geom_area() +
  labs(x="Distancia entre I0 y parte real del punto de máxima frecuencia \n DR",
       y="Longitud de la curva espectral \n P", 
       title= "Grafico de Area")

Se observa en el gráfico de áreas que las variables DR y P, están separando bien a las clases diferenciandolos por los colores.

Analisis Multivariado

Correlation matrix

corr=cor_pmat(BreastFinal[3:11])
head(corr)
##                 I0        PA500          HFS           DA         Area
## I0    0.000000e+00 2.984613e-05 7.721654e-01 6.593602e-27 4.305070e-10
## PA500 2.984613e-05 0.000000e+00 2.524054e-08 3.598703e-01 3.945145e-01
## HFS   7.721654e-01 2.524054e-08 0.000000e+00 2.750666e-01 3.407799e-02
## DA    6.593602e-27 3.598703e-01 2.750666e-01 0.000000e+00 5.676282e-19
## Area  4.305070e-10 3.945145e-01 3.407799e-02 5.676282e-19 0.000000e+00
## A.DA  3.148793e-12 1.778436e-02 1.799832e-04 5.757673e-14 3.814249e-28
##               A.DA       Max.IP           DR            P
## I0    3.148793e-12 2.254868e-27 3.996944e-19 1.542477e-87
## PA500 1.778436e-02 6.078945e-01 4.324086e-01 2.837045e-04
## HFS   1.799832e-04 9.111752e-05 9.061169e-01 2.964251e-01
## DA    5.757673e-14 1.234475e-20 4.634323e-69 2.278755e-22
## Area  3.814249e-28 2.859278e-19 1.901359e-15 1.247825e-10
## A.DA  0.000000e+00 3.739071e-26 2.190473e-09 1.190821e-15

Corellogram

ggcorrplot(corr) + labs(title = "CORELLOGRAM")

En esta matriz el mayor coeficiente de correlación se observa entre DR y HSF con 0.9061. Esto indica una fuerte relación lineal positiva entre DR y HSF. El segundo coeficiente de correlación más alto se observa entre entre I0 y HSF con 0.772. Esto indica una fuerte relación lineal positiva entre I0 y HSF.

Los coeficiente de correlación cercanos a cero se observan en los cruces de la variables donde los cuadros no están coloreados.

Eleccion de las 6 variables

Como la variable DR y HFS, I0 y HFS, DR y PA500, DA y PA500, tienen buena correlación entonces eliminamos a la variable HFS, DR y PA500, ya que las 6 variables que queden podrian representar un mejor modelo, ademas de los análisis graficos tambien podemos corroborar ello.

BreastFinal<-BreastFinal[,c(-1,-4,-5,-10)]

Data de entrenamiento y evaluación

70% Data de entrenamiento y 30% de evaluación

set.seed(888)

BreastFinal.train.idx <- sample(x = nrow(BreastFinal), size = nrow(BreastFinal)*0.70)
BreastFinal.train <- BreastFinal[BreastFinal.train.idx,]
BreastFinal.test <- BreastFinal[-BreastFinal.train.idx,]

Para la red neuroral realizaremos escalamiento dado que se corroboro en el diagrama de cajas ya que existe una diferencia significativa entre las escalas. Para ello creamos la función normalize que realiza un escalado de la variables en un rango [0,1]

normalize <- function(x) {
  return ((x - min(x)) / (max(x) - min(x)))
}
BreastFinal1<-BreastFinal
BreastFinal1.train<-BreastFinal.train
BreastFinal1.test<-BreastFinal.test

BreastFinal1$I0<-normalize(BreastFinal1$I0)
BreastFinal1$DA<-normalize(BreastFinal1$DA)
BreastFinal1$Area<-normalize(BreastFinal1$Area)
BreastFinal1$A.DA<-normalize(BreastFinal1$A.DA)
BreastFinal1$Max.IP<-normalize(BreastFinal1$Max.IP)
BreastFinal1$P<-normalize(BreastFinal1$P)

BreastFinal1.train$I0<-normalize(BreastFinal1.train$I0)
BreastFinal1.train$DA<-normalize(BreastFinal1.train$DA)
BreastFinal1.train$Area<-normalize(BreastFinal1.train$Area)
BreastFinal1.train$A.DA<-normalize(BreastFinal1.train$A.DA)
BreastFinal1.train$Max.IP<-normalize(BreastFinal1.train$Max.IP)
BreastFinal1.train$P<-normalize(BreastFinal1.train$P)

BreastFinal1.test$I0<-normalize(BreastFinal1.test$I0)
BreastFinal1.test$DA<-normalize(BreastFinal1.test$DA)
BreastFinal1.test$Area<-normalize(BreastFinal1.test$Area)
BreastFinal1.test$A.DA<-normalize(BreastFinal1.test$A.DA)
BreastFinal1.test$Max.IP<-normalize(BreastFinal1.test$Max.IP)
BreastFinal1.test$P<-normalize(BreastFinal1.test$P)

Agregamos las clases

BreastFinal1$car <- BreastFinal1$Class=="car"
BreastFinal1$fad <- BreastFinal1$Class=="fad"
BreastFinal1$mas <- BreastFinal1$Class=="mas"
BreastFinal1$gla <- BreastFinal1$Class=="gla"
BreastFinal1$con <- BreastFinal1$Class=="con"
BreastFinal1$adi <- BreastFinal1$Class=="adi"

BreastFinal1.train$car <- BreastFinal1.train$Class=="car"
BreastFinal1.train$fad <- BreastFinal1.train$Class=="fad"
BreastFinal1.train$mas <- BreastFinal1.train$Class=="mas"
BreastFinal1.train$gla <- BreastFinal1.train$Class=="gla"
BreastFinal1.train$con <- BreastFinal1.train$Class=="con"
BreastFinal1.train$adi <- BreastFinal1.train$Class=="adi"

BreastFinal1.test$car <- BreastFinal1.test$Class=="car"
BreastFinal1.test$fad <- BreastFinal1.test$Class=="fad"
BreastFinal1.test$mas <- BreastFinal1.test$Class=="mas"
BreastFinal1.test$gla <- BreastFinal1.test$Class=="gla"
BreastFinal1.test$con <- BreastFinal1.test$Class=="con"
BreastFinal1.test$adi <- BreastFinal1.test$Class=="adi"

Para el random forest es un modelo basado en árboles y, por lo tanto no requiere escalado.

BreastFinal2<-BreastFinal
BreastFinal2.train<-BreastFinal.train
BreastFinal2.test<-BreastFinal.test

Red Neuronal

Modelo

set.seed(111)

BreastFinal1.net <- neuralnet(car+fad+mas+gla+con+adi ~
                                I0+DA+Area+A.DA+Max.IP+P,
                              data=BreastFinal1.train, 
                              hidden=c(10,10,10),
                              act.fct = 'tanh',
                              threshold = 0.001)
par(mfrow = c(1,1))
plot(BreastFinal1.net, rep="best")

El modelo lo generamos con 3 capas ocultas cada uno con 10, 10, 10 neuronas respectivamente usamos la función de activacion “tanh” y la red neuronal continuará iterando hasta que la mejora en el error de entrenamiento sea menor o igual a 0.001.

Prediccion

BreastFinal1.prediccion <- predict(BreastFinal1.net, BreastFinal1.test[,-c(8:13)])

idx <- apply(BreastFinal1.prediccion, 1, which.max)
predicted <- c('car', 'fad', 'mas','gla', 'con', 'adi')[idx]
tab<-table(predicted, BreastFinal1.test$Class)
tab
##          
## predicted adi car con fad gla mas
##       adi  11   0   2   0   0   0
##       car   0   3   0   0   0   0
##       con   0   0   4   0   0   0
##       fad   0   0   0   2   0   0
##       gla   0   0   0   3   1   3
##       mas   0   0   0   1   2   0
sum(diag(tab)/sum(tab))
## [1] 0.65625

Se tiene una tasa de correcta clasificación del 65.625%.

Random Forest

Bosque aleatorio con 500 arboles y las 6 variables

set.seed(888)
modelo_bosque<-randomForest(Class~.,data=BreastFinal2.train,ntree=500,mtry=6)
modelo_bosque
## 
## Call:
##  randomForest(formula = Class ~ ., data = BreastFinal2.train,      ntree = 500, mtry = 6) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 6
## 
##         OOB estimate of  error rate: 28.38%
## Confusion matrix:
##     adi car con fad gla mas class.error
## adi  10   0   1   0   0   0  0.09090909
## car   0  15   0   0   0   3  0.16666667
## con   1   0   6   0   0   1  0.25000000
## fad   0   0   0   5   2   2  0.44444444
## gla   0   0   0   2  11   0  0.15384615
## mas   0   4   0   2   3   6  0.60000000
plot(modelo_bosque,main="Tasa de error Bosque Aleatorio")
legend(350,0.5,legend=colnames(modelo_bosque$err.rate),fill=1:9)

Se puede obervar como el error se mantiene constante a partir de los 500 arboles. Por lo que el uso de 500 árboles son suficientes para una óptima clasificación.

Número óptimo de variables

obb.values<-vector(length=(ncol(BreastFinal2)-1))
mtry<-vector(length=(ncol(BreastFinal2)-1))
for(i in 1:(ncol(BreastFinal2)-1)){
  variablesmod<-randomForest(Class~.,data=BreastFinal2.train,ntree=500,mtry=i)
  obb.values[i]<-variablesmod$err.rate[nrow(variablesmod$err.rate),1]
  mtry[i]<-i
}
data.frame(mtry,obb.values)

Se observa que con 4 variables se obtiene el mínimo error 0.2837838 Por lo que usar 3 variables es óptimo para el modelo.

Importancia de las variables

importance(modelo_bosque)
##        MeanDecreaseGini
## I0            12.978606
## DA            12.028147
## Area           6.342615
## A.DA           9.177859
## Max.IP         9.984883
## P              9.421946
varImpPlot(modelo_bosque,col="violetred4",n.var=(ncol(BreastFinal2)-1),main="Variables Importantes")

Se observa que las variables mas importantes según el coeficiente de Gini son I0, DA y Max.IP .Mientras que las otras 3 variables no serian tan importantes.

Validación del modelo con la data test

table(predict(modelo_bosque,BreastFinal2.test),BreastFinal2.test$Class)
##      
##       adi car con fad gla mas
##   adi  10   0   1   0   0   0
##   car   0   2   0   0   1   0
##   con   1   0   5   0   0   0
##   fad   0   0   0   3   0   0
##   gla   0   0   0   2   1   2
##   mas   0   1   0   1   1   1
sum(predict(modelo_bosque,BreastFinal2.test)==BreastFinal2.test$Class)/ length(BreastFinal2.test$Class)*100
## [1] 68.75

De los 32 datos de la data de evaluación, fueron predichos correctamente 22 datos. El porcentaje de aciertos que obtuvo el modelo es de 68.75%, el cual es un porcentaje muy bueno.

Conclusiones: