Trabajo 2 - Técnicas en aprendizaje estadístico
Camila Acosta Ramírez
Nicolás Amado Aristizabal
Sergio Andrés Arboleda Zuluaga
Esteban Palacio Londoño
Natalia Ramírez Ossa
Carlos Alejandro Torres Albarracín
01/12/2020
El presente trabajo muestra el desarrollo de algunos de los ejercicios prácticos propuestos en el libro “An Introduction to Statistical Learning with Applications in R”.
library(ggplot2)
library(pROC)
library(e1071)
library(ISLR)
library(ggplot2)
library(MLmetrics)
library(tree)
library(MASS)
library(class)
library(randomForest)
library(tidyverse)
library(glmnet)
library(gbm)Apartado 4.7
Punto 10
Introducción
En el siguiente estudio, se presentan resultados analíticos del set de datos Weekly del paquete ISLR, este set de datos contiene los datos del porcentaje de retorno del índice S&P 500 a lo largo de las semanas en 20 años, desde 1990 hasta 2010.
El set de datos contiene las variables año, el retorno porcentual de la semana, 5 variables del retorno porcentual de cada una de las semanas anteriores, una variable volumen (número promedio de acciones comerciadas diariamente, en billones), y finalmente la dirección, que indica si el retorno fue positivo o negativo dicha semana.
head(Weekly)## Year Lag1 Lag2 Lag3 Lag4 Lag5 Volume Today Direction
## 1 1990 0.816 1.572 -3.936 -0.229 -3.484 0.1549760 -0.270 Down
## 2 1990 -0.270 0.816 1.572 -3.936 -0.229 0.1485740 -2.576 Down
## 3 1990 -2.576 -0.270 0.816 1.572 -3.936 0.1598375 3.514 Up
## 4 1990 3.514 -2.576 -0.270 0.816 1.572 0.1616300 0.712 Up
## 5 1990 0.712 3.514 -2.576 -0.270 0.816 0.1537280 1.178 Up
## 6 1990 1.178 0.712 3.514 -2.576 -0.270 0.1544440 -1.372 DownLiteral a), Resumen de los datos
El resumen de los datos muestra un promedio de crecimiento porcentual positivo del día de \(0.14\%\), lo que indica que a lo largo del tiempo el porcentaje de retorno del S&P es bastante variable y tiende a mantener un equilibrio entre incremento y disminución del mismo.
La variable ‘Direction’ muestra \(605\) subidas contra \(484\) caídas, lo que indica que a la larga el S&P incrementa su valor en el tiempo.
summary(Weekly)## Year Lag1 Lag2 Lag3
## Min. :1990 Min. :-18.1950 Min. :-18.1950 Min. :-18.1950
## 1st Qu.:1995 1st Qu.: -1.1540 1st Qu.: -1.1540 1st Qu.: -1.1580
## Median :2000 Median : 0.2410 Median : 0.2410 Median : 0.2410
## Mean :2000 Mean : 0.1506 Mean : 0.1511 Mean : 0.1472
## 3rd Qu.:2005 3rd Qu.: 1.4050 3rd Qu.: 1.4090 3rd Qu.: 1.4090
## Max. :2010 Max. : 12.0260 Max. : 12.0260 Max. : 12.0260
## Lag4 Lag5 Volume Today
## Min. :-18.1950 Min. :-18.1950 Min. :0.08747 Min. :-18.1950
## 1st Qu.: -1.1580 1st Qu.: -1.1660 1st Qu.:0.33202 1st Qu.: -1.1540
## Median : 0.2380 Median : 0.2340 Median :1.00268 Median : 0.2410
## Mean : 0.1458 Mean : 0.1399 Mean :1.57462 Mean : 0.1499
## 3rd Qu.: 1.4090 3rd Qu.: 1.4050 3rd Qu.:2.05373 3rd Qu.: 1.4050
## Max. : 12.0260 Max. : 12.0260 Max. :9.32821 Max. : 12.0260
## Direction
## Down:484
## Up :605
##
##
##
## En el siguiente Diagrama de Caja se pueden identificar la varianza de los datos en cada año
Los años con una varianza menor se podrían relacionar con años con poca incertidumbre económica, mientras que los años con una varianza mayor podrían relacionarse con años en los que el S&P se comportó de una manera más volátil, un ejemplo de esto es el año 2008, donde ocurrió una crisis financiera y una crisis bursátil mundial
Weekly$YFactor <- as.factor(Weekly$Year) # Año como factor para producir un Box Plot
ggplot(data= Weekly, mapping = aes(x = YFactor, y= Today)) +
geom_boxplot() +
xlab("Año") +
ylab("Porcentaje de Retorno")
Gráfica histórica de los puntos
La gráfica histórica de los puntos muestra el constante incremento y disminución del porcentaje de retorno, sobresalen los puntos alrededor del año 2008 debido a la crisis financiera ya mencionada.
ggplot(data= Weekly, mapping = aes(x = c(1:length(Year)), y= Today)) +
geom_line() +
geom_point() +
xlab("Registro") +
ylab("Porcentaje de Retorno") +
labs(title="Serie histórica de Porcentaje de Retorno")
Se realiza una gráfica con menos puntos para poder visualizar de una mejor manera.
ggplot(data= Weekly[1:300,], mapping = aes(x = c(1:length(Year)), y= Today)) +
geom_line() +
geom_point() +
xlab("Registro") +
ylab("Porcentaje de Retorno") +
labs(title="Serie histórica de Porcentaje de Retorno")
Literal b), Regresión logística
Se realiza una regresión logística con todos los datos, tienendo como variable respuesta la dirección, las 5 variables de delay (Lag) más el volumen (Volume) como predictores.
weekly_glm <- glm(Direction ~ Lag1+Lag2+Lag3+Lag4+Lag5+Volume,
data= Weekly, family="binomial")
summary(weekly_glm)##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +
## Volume, family = "binomial", data = Weekly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6949 -1.2565 0.9913 1.0849 1.4579
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.26686 0.08593 3.106 0.0019 **
## Lag1 -0.04127 0.02641 -1.563 0.1181
## Lag2 0.05844 0.02686 2.175 0.0296 *
## Lag3 -0.01606 0.02666 -0.602 0.5469
## Lag4 -0.02779 0.02646 -1.050 0.2937
## Lag5 -0.01447 0.02638 -0.549 0.5833
## Volume -0.02274 0.03690 -0.616 0.5377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1496.2 on 1088 degrees of freedom
## Residual deviance: 1486.4 on 1082 degrees of freedom
## AIC: 1500.4
##
## Number of Fisher Scoring iterations: 4El único predictor que parece ser un estadísticamente significante es la variable Lag2, pero no demuestra mayor significancia.
Literal c), Matriz de confusión
contrasts(as.factor(Weekly$Direction)) # Nos indica que Down es un 0 y Up es un 1, para interpretar los resultados de la predicción## Up
## Down 0
## Up 1La matriz de confusión demuestra que el error más grande son los falsos positivos, en este caso el positivo se da cuando se predice un incremento en el porcentaje de retorno, y el falso es que realmente era un decrecimiento, la matriz de confusión también muestra que el modelo en la mayor parte de los datos predice un incremento, y tan solo alrededor de 100 veces predice una disminución.
predicted <- round(predict(weekly_glm, newdata=Weekly, type="response")) # Predicción
real <- weekly_glm$y # Reales
c_matrix <- ConfusionMatrix(predicted, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos es de %f%s", Accuracy(predicted,real)*100, "%"))## [1] "El porcentaje de aciertos es de 56.106520%"c_matrix## y_pred
## y_true 0 1
## 0 54 430
## 1 48 557Literal d), Regresión logística 2
Se realiza otra regresión logística, teniendo en cuenta esta vez solamente la variable Lag2, y usando los datos de entrenamiento desde el 1990 hasta el 2008.
weekly_tr <- subset(Weekly, Year <= 2008) # Datos de entrenamiento
weekly_vl <- subset(Weekly, Year > 2008) # Datos de validaciónweekly_glm2 <- glm(Direction ~ Lag2, data= weekly_tr, family="binomial")
summary(weekly_glm2)##
## Call:
## glm(formula = Direction ~ Lag2, family = "binomial", data = weekly_tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.536 -1.264 1.021 1.091 1.368
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.20326 0.06428 3.162 0.00157 **
## Lag2 0.05810 0.02870 2.024 0.04298 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.5 on 983 degrees of freedom
## AIC: 1354.5
##
## Number of Fisher Scoring iterations: 4La matriz de confusión muestra unos resultados similares a los anteriores, en los que el mayor error se presenta en los falsos positivos.
predicted <- round(predict(weekly_glm2, newdata=weekly_vl, type="response")) # Predicción
real <- ifelse(weekly_vl$Direction=="Up",1,0) # Reales
c_matrix2 <- ConfusionMatrix(predicted, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos es de %f%s", Accuracy(predicted,real)*100, "%"))## [1] "El porcentaje de aciertos es de 62.500000%"c_matrix2## y_pred
## y_true 0 1
## 0 9 34
## 1 5 56Literal e), LDA
weekly_lda <- lda(Direction ~ Lag2, data = weekly_tr)
weekly_lda## Call:
## lda(Direction ~ Lag2, data = weekly_tr)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
##
## Coefficients of linear discriminants:
## LD1
## Lag2 0.4414162Las predicciones realizadas por el LDA son identicas a las realizadas por el GLM.
predicted <- predict(weekly_lda, newdata=weekly_vl) # Predicción
real <- weekly_vl$Direction # Reales
c_matrix_lda <- ConfusionMatrix(predicted$class, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos es de %f%s", Accuracy(predicted$class,real)*100, "%"))## [1] "El porcentaje de aciertos es de 62.500000%"c_matrix_lda## y_pred
## y_true Down Up
## Down 9 34
## Up 5 56Literal f), QDA
weekly_qda <- qda(Direction ~ Lag2, data = weekly_tr)
weekly_qda## Call:
## qda(Direction ~ Lag2, data = weekly_tr)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581El QDA muestra los peores resultados en terminos de generalidad ya que no marca ninguna disminución.
predicted <- predict(weekly_qda, newdata=weekly_vl) # Predicción
real <- weekly_vl$Direction # Reales
c_matrix_lda <- ConfusionMatrix(predicted$class, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos es de %f%s", Accuracy(predicted$class,real)*100, "%"))## [1] "El porcentaje de aciertos es de 58.653846%"c_matrix_lda## y_pred
## y_true Down Up
## Down 0 43
## Up 0 61Literal g), KNN con k=1
# Preparación de los datos
x_tr <- cbind(weekly_tr$Lag2)
x_vl <- cbind(weekly_vl$Lag2)
y_tr <- weekly_tr$Direction
y_vl <- weekly_vl$DirectionLa aproximación por KNN muestra unos resultados más alejados a los otros modelos, y predice correctamente el \(50\%\) de los datos, con más variabilidad entre respuestas de incremento y disminución.
set.seed(2020)
predicted = knn(train = x_tr, test = x_vl, cl = y_tr, k = 1)
print(sprintf("El porcentaje de aciertos es de %f%s", Accuracy(predicted,y_vl)*100, "%"))## [1] "El porcentaje de aciertos es de 50.000000%"table(predicted, y_vl)## y_vl
## predicted Down Up
## Down 21 30
## Up 22 31Literal h), Resultados
Los modelos lineales (GLM y LDA) son los que demuestran el mejor porcentaje de aciertos, ambos con \(65\) aciertos, aunque con un claro desbalance a predecir incrementos en el porcentaje de retorno. El QDA es el modelo más desbalanceados ya que no predice ni una sola disminución. Por último el modelo de agrupación con KNN muestra los resultados con mejor generalidad ya que tiene una buena cantidad de aciertos \(52\) y estos se reparten en predicciones de incremento y disminución.
Literal i), Experimentación
Agregando la variable volumen como predictor
Agregar la variable volumen, aunque no aumente el porcentaje de aciertos, si le agrega cierto nivel de variación a las predicciones, ya que se predicen más disminuciones que en las aproximaciones anteriores. Se obtiene el mejor resultado al realizar KNN con \(k=1\) con un \(55\%\) de aciertos, aún así siendo un resultado bastante bajo para un problema de clasificación binaria.
GLM
weekly_glm <- glm(Direction ~ Lag2 + Volume, data = weekly_tr,
family = "binomial")
predicted <- round(predict(weekly_glm, newdata = weekly_vl, type="response"))
real <- ifelse(weekly_vl$Direction=="Up",1,0)
print(sprintf("El porcentaje de aciertos es de %f%s", Accuracy(predicted,real)*100, "%"))## [1] "El porcentaje de aciertos es de 53.846154%"ConfusionMatrix(predicted, real)## y_pred
## y_true 0 1
## 0 20 23
## 1 25 36LDA
weekly_lda <- lda(Direction ~ Lag2 + Volume, data= weekly_tr)
predicted <- predict(weekly_lda, newdata=weekly_vl) # Predicción
real <- weekly_vl$Direction # Reales
c_matrix_lda <- ConfusionMatrix(predicted$class, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos es de %f%s", Accuracy(predicted$class,real)*100, "%"))## [1] "El porcentaje de aciertos es de 53.846154%"c_matrix_lda## y_pred
## y_true Down Up
## Down 20 23
## Up 25 36QDA
weekly_qda <- qda(Direction ~ Lag2 + Volume, data= weekly_tr)
predicted <- predict(weekly_qda, newdata=weekly_vl) # Predicción
real <- weekly_vl$Direction # Reales
c_matrix_lda <- ConfusionMatrix(predicted$class, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos es de %f%s", Accuracy(predicted$class,real)*100, "%"))## [1] "El porcentaje de aciertos es de 47.115385%"c_matrix_lda## y_pred
## y_true Down Up
## Down 32 11
## Up 44 17KNN
# Preparación de los datos
x_tr <- cbind(weekly_tr$Lag2,weekly_tr$Volume)
x_vl <- cbind(weekly_vl$Lag2,weekly_vl$Volume)
y_tr <- weekly_tr$Direction
y_vl <- weekly_vl$Direction
for (k in c(1:10)){
set.seed(2020)
predicted = knn(train = x_tr, test = x_vl, cl = y_tr, k = k)
print(sprintf("El porcentaje de aciertos para k=%i es de %f%s",k,Accuracy(predicted,y_vl)*100, "%"))
}## [1] "El porcentaje de aciertos para k=1 es de 55.769231%"
## [1] "El porcentaje de aciertos para k=2 es de 54.807692%"
## [1] "El porcentaje de aciertos para k=3 es de 54.807692%"
## [1] "El porcentaje de aciertos para k=4 es de 48.076923%"
## [1] "El porcentaje de aciertos para k=5 es de 52.884615%"
## [1] "El porcentaje de aciertos para k=6 es de 48.076923%"
## [1] "El porcentaje de aciertos para k=7 es de 47.115385%"
## [1] "El porcentaje de aciertos para k=8 es de 50.961538%"
## [1] "El porcentaje de aciertos para k=9 es de 50.000000%"
## [1] "El porcentaje de aciertos para k=10 es de 52.884615%"KNN con \(k=1\) ha sido el método que ha demostrado el mejor comportamiento, y por ende se seguirá experimentando con el mismo.
KNN con k=1, usando todas las variables
El porcentaje de aciertos disminuye al usar todas las variables con KNN \(k=1\).
x_tr <- cbind(weekly_tr$Lag1,weekly_tr$Lag2, weekly_tr$Lag3, weekly_tr$Lag4,weekly_tr$Lag5,weekly_tr$Volume)
x_vl <- cbind(weekly_vl$Lag1,weekly_vl$Lag2, weekly_vl$Lag3, weekly_vl$Lag4,weekly_vl$Lag5,weekly_vl$Volume)
y_tr <- weekly_tr$Direction
y_vl <- weekly_vl$Direction
set.seed(2020)
predicted = knn(train = x_tr, test = x_vl, cl = y_tr, k = 1)
print(sprintf("El porcentaje de aciertos es de %f%s",Accuracy(predicted,y_vl)*100, "%"))## [1] "El porcentaje de aciertos es de 48.076923%"ConfusionMatrix(predicted, y_vl)## y_pred
## y_true Down Up
## Down 21 22
## Up 32 29KNN con k=1, variando los predictores de delay
Se realiza una comparación del comportamiento del modelo variando los predictores.
y_tr <- weekly_tr$Direction
y_vl <- weekly_vl$Direction
predictors <- c("Lag1", "Lag2", "Lag3", "Lag4", "Lag5")
for (i in predictors){
x_tr <- cbind(weekly_tr[[as.name(i)]], weekly_tr$Volume)
x_vl <- cbind(weekly_vl[[as.name(i)]], weekly_vl$Volume)
set.seed(2020)
predicted = knn(train = x_tr, test = x_vl, cl = y_tr, k = 1)
print(sprintf("El porcentaje de aciertos con el predictor %s es de %f%s",i,Accuracy(predicted,y_vl)*100, "%"))
}## [1] "El porcentaje de aciertos con el predictor Lag1 es de 52.884615%"
## [1] "El porcentaje de aciertos con el predictor Lag2 es de 55.769231%"
## [1] "El porcentaje de aciertos con el predictor Lag3 es de 46.153846%"
## [1] "El porcentaje de aciertos con el predictor Lag4 es de 45.192308%"
## [1] "El porcentaje de aciertos con el predictor Lag5 es de 51.923077%"Conclusión
Realizando diversos modelos variando los predictores de delay, se ha encontrado que el modelo más efectivo es el KNN con \(k=1\), utilizando Lag2 y Volume como predictores. Aún así ninguno de los modelos refleja buenos resultados, como máximo se alcanza un \(55\%\), porcentaje bastante bajo para un problema de clasificación binario. Es posible que se necesiten más datos, o sea necesario realizar limpieza en los datos para expulsar los datos anomalos.
Punto 11
Introducción
En este problema el objetivo es obtener un modelo para predecir si un automóvil obtiene un buen rendimiento de gasolina, basándose en el set de datos Auto.
El dataset cuenta con varios datos acerca de las especificaciones del automóvil, y las millas por galón como el dato que da a entender el rendimiento del carro.
summary(Auto)## mpg cylinders displacement horsepower weight
## Min. : 9.00 Min. :3.000 Min. : 68.0 Min. : 46.0 Min. :1613
## 1st Qu.:17.00 1st Qu.:4.000 1st Qu.:105.0 1st Qu.: 75.0 1st Qu.:2225
## Median :22.75 Median :4.000 Median :151.0 Median : 93.5 Median :2804
## Mean :23.45 Mean :5.472 Mean :194.4 Mean :104.5 Mean :2978
## 3rd Qu.:29.00 3rd Qu.:8.000 3rd Qu.:275.8 3rd Qu.:126.0 3rd Qu.:3615
## Max. :46.60 Max. :8.000 Max. :455.0 Max. :230.0 Max. :5140
##
## acceleration year origin name
## Min. : 8.00 Min. :70.00 Min. :1.000 amc matador : 5
## 1st Qu.:13.78 1st Qu.:73.00 1st Qu.:1.000 ford pinto : 5
## Median :15.50 Median :76.00 Median :1.000 toyota corolla : 5
## Mean :15.54 Mean :75.98 Mean :1.577 amc gremlin : 4
## 3rd Qu.:17.02 3rd Qu.:79.00 3rd Qu.:2.000 amc hornet : 4
## Max. :24.80 Max. :82.00 Max. :3.000 chevrolet chevette: 4
## (Other) :365Literal a), Creación de la variable binaria mpg01
Se procede a crear una variable binaria mpg01 la cuál es verdadera si el mpg es superior a la mediana del mpg a lo largo del dataset, la mediana en este caso es de \(22.75\).
Auto$mpg01 <- (Auto$mpg >= median(Auto$mpg))
median(Auto$mpg)## [1] 22.75Efectivamente el ejercicio anterior resulta en una división igualitaria de los autos.
summary(Auto$mpg01)## Mode FALSE TRUE
## logical 196 196FALSE indica un \(0\) y TRUE indica un \(1\).
contrasts(as.factor(Auto$mpg01))## TRUE
## FALSE 0
## TRUE 1Literal b), Análisis exploratorio
Se realizará un análisis de cada variable contra la variable mpg para decidir cuáles serán los predictores usados en los modelos.
Número de Cilindros
Se realiza una gráfica de puntos agregando la línea roja de la mediana, donde los valores encima o sobre esta tienen un valor de \(mpg01 = 1\) y las que estén por debajo un valor de \(mpg01 = 0\).
La gráfica muestra que para cada número de cilindros es posible tener un buen rendimiento de mpg, destaca que los autos con 4 cilindros tienen más ocurrencias en valores altos de mpg, y los de 8 cilindros en valores bajos, se alcanza a notar una pequeña tendencia a entre más cilindros, peor rendimiento de mpg, pero no es del todo clara.
# Cilindros
ggplot(data = Auto, mapping=aes(x = cylinders, y=mpg)) + geom_point() + geom_line(mapping= aes(x = cylinders,y = median(mpg), col="Mediana")) + labs(title="Comportamiento del No. de Cilindros",x = "Número de Cilindros", y = "Millas por galón", color="Leyenda") + scale_color_manual(values = c("Mediana" = "red")) 
Desplazamiento
En esta gráfica se nota una clara tendencia en la que los vehículos con un desplazamiento más bajo muestran el mejor rendimiento de mpg.
# Displacement
displacement_g <- ggplot(data = Auto, mapping=aes(x = displacement, y=mpg))
displacement_g + geom_point() + geom_line(mapping= aes(x = displacement,y = median(mpg), col="Mediana")) + labs(title="Comportamiento del Desplazamiento" ,x = "Desplazamiento en Pulgadas Cúbicas", y = "Millas por galón", color="Leyenda") + scale_color_manual(values = c("Mediana" = "red")) 
Caballos de fuerza
Al igual que con el desplazamiento, al tener menos caballos de fuerza mejor es el rendimiento en mpg.
horsepower_g <- ggplot(data = Auto, mapping=aes(x = horsepower, y=mpg))
horsepower_g + geom_point() + geom_line(mapping= aes(x = horsepower,y = median(mpg), col="Mediana")) + labs(title="Comportamiento de los Caballos de Fuerza",x = "Caballos de fuerza", y = "Millas por galón", color="Leyenda") + scale_color_manual(values = c("Mediana" = "red"))
Masa
Al igual que el desplazamiento y los caballos de fuerza, los carros más livianos presentan un mejor rendimiento en mpg. Este comportamiento de las 3 variables podría sugerir una correlación de las mismas.
weight_g <- ggplot(data = Auto, mapping=aes(x = weight, y=mpg))
weight_g + geom_point() + geom_line(mapping= aes(x = weight,y = median(mpg), col="Mediana")) + labs(title="Comportamiento de la Masa",x = "Masa del carro en libras", y = "Millas por galón", color="Leyenda") + scale_color_manual(values = c("Mediana" = "red"))
Tiempo de aceleración de 0 a 60 mph
La gráfica no sugiere una relación directa entre el tiempo de aceleración, al menos para los valores entre 11 y 23 segundos, aproximadamente, sin embargo se identifica que los autos que presentan una aceleración de 10 o menos segundos generalmente presentan un mal rendimiento en mpg.
accel_g <- ggplot(data = Auto, mapping=aes(x = acceleration, y=mpg))
accel_g + geom_point() + geom_line(mapping= aes(x = acceleration,y = median(mpg), col="Mediana")) + labs(titile="Comportamiento de la aceleración",x = "Tiempo en segundos para acelerar a 60 mph", y = "Millas por galón", color="Leyenda") + scale_color_manual(values = c("Mediana" = "red"))
País de origen
La diferenciación por país de origen demuestra una clara diferencia en el rendimiento en mpg, donde los autos americanos presentan el peor, y los europeos y japoneses el mejor.
# Factor de origen
Auto$forigin <- as.factor(Auto$origin)
origin_g <- ggplot(data = Auto, mapping=aes(x = forigin, y=mpg))
origin_g + geom_boxplot() + geom_line(mapping= aes(x = forigin,y = median(mpg), col="Mediana")) + labs(title="Comportamiento del Origen del Automóvil",x = "Origen del Automóvil", y = "Millas por galón", color="Leyenda") + scale_x_discrete(labels = c("1" = "Americano", "2" = "Europeo", "3" = "Japonés"))
Año del modelo
Se muestra que a medida que aumenta el año se reduce el mpg, sin embargo no hay una diferenciación absoluta, todos los años tienen al menos un carro con buen rendimiento en mpg y con mal rendimiento en mpg.
year_g <- ggplot(data = Auto, mapping=aes(x = year, y=mpg))
year_g + geom_point() + geom_line(mapping= aes(x = year,y = median(mpg), col="Mediana")) + labs(title="Comportamiento por Año" ,x = "Año", y = "Millas por galón", color="Leyenda") + scale_color_manual(values = c("Mediana" = "red")) + scale_x_continuous(breaks = c(70:82))
Nombre del carro
La variable de nombre del carro posee muchos niveles y una diferenciación por el mismo no sería del todo buena, ya que no se cuentan con muchos datos.
length(levels(Auto$name))## [1] 304Elección de predictores
Los predictores elegidos según el análisis serán: - Desplazamiento (displacement) - Caballos de fuerza (horsepower) - Masa del vehículo (weight) - Origen del auto (forigin)
formula <- mpg01 ~ displacement + horsepower + weight + foriginLiteral c), División del set de datos
Para la división del set de datos se toma una muestra aleatoria del \(80\%\) para entrenamiento y el \(20\%\) restante para validación.
# Se encuentra el tamaño del set de datos (tomando el 80% para entrenamiento)
sample_size = floor(0.8*nrow(Auto))
set.seed(2020)
# Dividir el set de datos aleatoriamente
picked = sample(seq_len(nrow(Auto)),size = sample_size)
auto_tr =Auto[picked,]
auto_vl =Auto[-picked,]Literal d), LDA
auto_lda <- lda(formula, data = auto_tr)
auto_lda## Call:
## lda(formula, data = auto_tr)
##
## Prior probabilities of groups:
## FALSE TRUE
## 0.5111821 0.4888179
##
## Group means:
## displacement horsepower weight forigin2 forigin3
## FALSE 273.8250 130.18125 3634.881 0.0687500 0.0500000
## TRUE 113.3889 78.35294 2312.418 0.2941176 0.3856209
##
## Coefficients of linear discriminants:
## LD1
## displacement -0.006321321
## horsepower 0.005530063
## weight -0.001196922
## forigin2 0.366420649
## forigin3 0.394776068predicted <- predict(auto_lda, newdata=auto_vl) # Predicción
real <- auto_vl$mpg01 # Reales
c_matrix_lda <- ConfusionMatrix(predicted$class, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos para es de %f%s", Accuracy(predicted$class,real)*100, "%"))## [1] "El porcentaje de aciertos para es de 87.341772%"c_matrix_lda## y_pred
## y_true FALSE TRUE
## FALSE 29 7
## TRUE 3 40Literal e), QDA
auto_qda <- qda(formula, data = auto_tr)
auto_qda## Call:
## qda(formula, data = auto_tr)
##
## Prior probabilities of groups:
## FALSE TRUE
## 0.5111821 0.4888179
##
## Group means:
## displacement horsepower weight forigin2 forigin3
## FALSE 273.8250 130.18125 3634.881 0.0687500 0.0500000
## TRUE 113.3889 78.35294 2312.418 0.2941176 0.3856209predicted <- predict(auto_qda, newdata=auto_vl) # Predicción
real <- auto_vl$mpg01 # Reales
c_matrix_qda <- ConfusionMatrix(predicted$class, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos para es de %f%s", Accuracy(predicted$class,real)*100, "%"))## [1] "El porcentaje de aciertos para es de 89.873418%"c_matrix_qda## y_pred
## y_true FALSE TRUE
## FALSE 30 6
## TRUE 2 41Literal f), Regresión logística
auto_glm <- glm(formula, data= auto_tr, family="binomial")
summary(auto_glm)##
## Call:
## glm(formula = formula, family = "binomial", data = auto_tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4910 -0.2100 -0.0057 0.3187 3.3680
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 12.1057415 1.8999057 6.372 1.87e-10 ***
## displacement -0.0101779 0.0081147 -1.254 0.2097
## horsepower -0.0323690 0.0158900 -2.037 0.0416 *
## weight -0.0026066 0.0009031 -2.886 0.0039 **
## forigin2 0.2320536 0.6745525 0.344 0.7308
## forigin3 0.1129592 0.6377234 0.177 0.8594
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 433.75 on 312 degrees of freedom
## Residual deviance: 156.08 on 307 degrees of freedom
## AIC: 168.08
##
## Number of Fisher Scoring iterations: 7predicted <- round(predict(auto_glm, newdata=auto_vl, type="response")) # Predicción
real <- auto_vl$mpg01 # Reales
c_matrix_glm <- ConfusionMatrix(predicted, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos para es de %f%s", Accuracy(predicted,real)*100, "%"))## [1] "El porcentaje de aciertos para es de 86.075949%"c_matrix_glm## y_pred
## y_true 0 1
## FALSE 29 7
## TRUE 4 39# Preparación de los datos
x_tr <- cbind(auto_tr$displacement,auto_tr$horsepower, auto_tr$weight, auto_tr$forigin)
x_vl <- cbind(auto_vl$displacement,auto_vl$horsepower, auto_vl$weight, auto_vl$forigin)
y_tr <- auto_tr$mpg01
y_vl <- auto_vl$mpg01La aproximación por KNN muestra unos resultados más alejados al resto en términos de efectividad, los \(k\) que presentaron el mejor desempeño fueron los \(k = 4\) y \(k = 7\)
set.seed(2020)
for (k in c(1:10)) {
predict = knn(train = x_tr, test = x_vl, cl = y_tr, k = k)
print(sprintf("El porcentaje de aciertos para el k=%i es de %f%s", k, Accuracy(predict,y_vl)*100, "%"))
}## [1] "El porcentaje de aciertos para el k=1 es de 84.810127%"
## [1] "El porcentaje de aciertos para el k=2 es de 82.278481%"
## [1] "El porcentaje de aciertos para el k=3 es de 82.278481%"
## [1] "El porcentaje de aciertos para el k=4 es de 86.075949%"
## [1] "El porcentaje de aciertos para el k=5 es de 83.544304%"
## [1] "El porcentaje de aciertos para el k=6 es de 84.810127%"
## [1] "El porcentaje de aciertos para el k=7 es de 86.075949%"
## [1] "El porcentaje de aciertos para el k=8 es de 84.810127%"
## [1] "El porcentaje de aciertos para el k=9 es de 83.544304%"
## [1] "El porcentaje de aciertos para el k=10 es de 82.278481%"Conclusión
El modelo que obtuvo un mejor porcentaje de aciertos fue el del QDA, seguido del LDA, por último la regresión logística y el mejor modelo de KNN obtuvieron el mismo porcentaje de aciertos.
Punto 12
Literal a), Funcion Power()
Power <- function(){
print(2^3)
}
Power()## [1] 8Literal b), Funcion Power2()
Power2 <- function(x,a){
print(x^a)
}
Power2(3,8)## [1] 6561Literal c), Uso de Power2()
Power2(10,3)## [1] 1000Power2(8,17)## [1] 2.2518e+15Power2(131,3)## [1] 2248091Literal d), Función Power3()
Power3 <- function(x,a){
return(x^a)
}Literal e), Gráfica de f(x)
plot(x = c(1:10), y = Power3(c(1:10), 2), type = "b", main = "f(x) = x^2", xlab = "x", ylab="y", las = 1)
plot(x = c(1:10), y = Power3(c(1:10), 2), type = "b", main = "f(x) = x^2", xlab = "x", ylab="y", las = 1, log="xy")
Literal f), Función PlotPower()
PlotPower <- function(x,a){
plot(x = x, y = Power3(x, a), type = "b", main = sprintf("f(x) = x^%i", a), xlab = "x", ylab="y", las = 1)
}
PlotPower(c(1:15), 3)
Punto 13
Introducción
En el siguiente punto el objetivo es predecir a partir de los datos si un suburbio tiene una tasa de criminalidad superior a la mediana. Se usará el set de datos Boston del paquete MASS.
Resumen y variables
El dataset cuenta con las siguientes variables:
- crim: Tasa de criminalidad per cápita.
- zn: Proporción de terreno residencial dividido en lotes, sobre 25000 pies cuadrados.
- indus: Proporción de areas de negocios no-minoristas por suburbio.
- chas: Variable dummy donde es 1 si el suburbio limita con el río Charles y 0 si no.
- nox: Concetración de óxido de nitrógeno en partes por 10 millones.
- rm: Número promedio de habitaciones por vivienda.
- age: Proporción de viviendas ocupadas por sus dueños, construidas antes de 1940.
- dis: Media ponderada de la distancia a 5 centros de empleo de Boston.
- rad: Índice de accesibilidad a carreteras radiales.
- ptratio: Radio de pupílo-maestro por suburbio.
- black: \(1000(Bk - 0.63)^2\) donde \(Bk\) es la proporción de negros en el suburbio.
- lstat: Status bajo de la población (porcentaje).
- medv: Mediana de viviendas ocupadas por sus dueños en 1000$
summary(Boston)## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08205 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00Variable crim01
Se procede a crear una variable binaria crim01 la cuál es verdadera si la tasa de criminalidad es superior a la mediana de la tasa a lo largo del dataset, la mediana en este caso es de \(0.25651\)
Boston$crim01 <- (Boston$crim >= median(Boston$crim))
median(Boston$crim)## [1] 0.25651FALSE indica un \(0\) y TRUE indica un \(1\)
contrasts(as.factor(Boston$crim01))## TRUE
## FALSE 0
## TRUE 1División del set de datos
Para la división del set de datos se toma una muestra aleatoria del \(80\%\) para entrenamiento y el \(20\%\) restante para validación.
# Se encuentra el tamaño del set de datos (tomando el 80% para entrenamiento)
sample_size = floor(0.8*nrow(Boston))
set.seed(2020)
# Dividir el set de datos aleatoriamente
picked = sample(seq_len(nrow(Boston)),size = sample_size)
boston_tr =Boston[picked,]
boston_vl =Boston[-picked,]Se realizarán diversos modelos con regresión logística, utilizando diferentes sets de predictores, los predictores se separarán según su tipo:
- Geográfico: zn, indus, chas, dis, rad
- Vivienda: rm, age, medv
- Poblacional: ptratio, black, lstat
Las variables nox y tax, se incluiran en los modelos finales ya que no se encasillan en ninguno de los grupos
De cada modelo luego se tomaran las variables que presentaron la mayor significancia y estas serán las que se usarán para cada uno de los modelos posteriores.
Regresión logística, predictores geográficos.
Las variables con mayor significancia fueron dis,rad, zn.
model_glm1 <- glm(crim01 ~ zn + indus + chas + dis + rad, data = boston_tr, family = "binomial")
summary(model_glm1)##
## Call:
## glm(formula = crim01 ~ zn + indus + chas + dis + rad, family = "binomial",
## data = boston_tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.72341 -0.51300 0.00348 0.00806 2.53603
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.02371 0.87044 -1.176 0.2396
## zn -0.03332 0.01523 -2.188 0.0287 *
## indus 0.03820 0.02808 1.360 0.1738
## chas 0.27502 0.55134 0.499 0.6179
## dis -0.60935 0.12642 -4.820 1.44e-06 ***
## rad 0.54076 0.11702 4.621 3.81e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 560.02 on 403 degrees of freedom
## Residual deviance: 255.65 on 398 degrees of freedom
## AIC: 267.65
##
## Number of Fisher Scoring iterations: 8predicted <- round(predict(model_glm1, newdata = boston_vl, type = "response"))
real <- boston_vl$crim01
acc <- Accuracy(predicted, real)* 100
print(sprintf("El porcentaje de aciertos fue del %f%s", acc, "%"))## [1] "El porcentaje de aciertos fue del 81.372549%"ConfusionMatrix(predicted, real) # Matriz de confusión## y_pred
## y_true 0 1
## FALSE 44 9
## TRUE 10 39Para las variables de vivienda, la de mayor significancia fue age.
model_glm2 <- glm(crim01 ~ rm + age + medv, data= boston_tr, family = "binomial")
summary(model_glm2)##
## Call:
## glm(formula = crim01 ~ rm + age + medv, family = "binomial",
## data = boston_tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9926 -0.6806 0.5119 0.6675 2.7922
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.807012 1.417786 -3.391 0.000698 ***
## rm 0.091781 0.252728 0.363 0.716486
## age 0.062754 0.006403 9.801 < 2e-16 ***
## medv -0.007696 0.020459 -0.376 0.706793
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 560.02 on 403 degrees of freedom
## Residual deviance: 378.67 on 400 degrees of freedom
## AIC: 386.67
##
## Number of Fisher Scoring iterations: 5predicted <- round(predict(model_glm2, newdata = boston_vl, type = "response"))
real <- boston_vl$crim01
acc <- Accuracy(predicted, real)* 100
print(sprintf("El porcentaje de aciertos fue del %f%s", acc, "%"))## [1] "El porcentaje de aciertos fue del 77.450980%"ConfusionMatrix(predicted, real) # Matriz de confusión## y_pred
## y_true 0 1
## FALSE 34 19
## TRUE 4 45Para las variables poblacionales, las de mayor significancia fueron las variables black y lstat.
model_glm3 <- glm(crim01 ~ ptratio + black + lstat, data = boston_tr, family = "binomial")
summary(model_glm3)##
## Call:
## glm(formula = crim01 ~ ptratio + black + lstat, family = "binomial",
## data = boston_tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -3.5860 -0.8707 0.0107 0.9099 1.9611
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.685163 2.485582 2.287 0.022181 *
## ptratio 0.034861 0.058975 0.591 0.554444
## black -0.021138 0.006211 -3.403 0.000666 ***
## lstat 0.142127 0.023182 6.131 8.74e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 560.02 on 403 degrees of freedom
## Residual deviance: 430.38 on 400 degrees of freedom
## AIC: 438.38
##
## Number of Fisher Scoring iterations: 7predicted <- round(predict(model_glm3, newdata = boston_vl, type = "response"))
real <- boston_vl$crim01
acc <- Accuracy(predicted, real)* 100
print(sprintf("El porcentaje de aciertos fue del %f%s", acc, "%"))## [1] "El porcentaje de aciertos fue del 74.509804%"ConfusionMatrix(predicted, real) # Matriz de confusión## y_pred
## y_true 0 1
## FALSE 40 13
## TRUE 13 36Predictores elegidos
Los predictores elegidos serán: - dis - rad - age - black - lstat - nox - tax
Regresión logística
Se realiza una regresión logística con los predictores escogidos.
model_glm <- glm(crim01 ~ dis + rad + age + black + lstat + nox + tax, data = boston_tr, family = "binomial")
summary(model_glm)##
## Call:
## glm(formula = crim01 ~ dis + rad + age + black + lstat + nox +
## tax, family = "binomial", data = boston_tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.90854 -0.25996 0.00008 0.00414 2.61963
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -19.450861 4.350661 -4.471 7.79e-06 ***
## dis 0.136719 0.168015 0.814 0.415800
## rad 0.666387 0.133441 4.994 5.92e-07 ***
## age 0.007796 0.010460 0.745 0.456059
## black -0.007433 0.005568 -1.335 0.181947
## lstat 0.004825 0.039166 0.123 0.901955
## nox 39.128187 6.635208 5.897 3.70e-09 ***
## tax -0.009910 0.002701 -3.669 0.000243 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 560.02 on 403 degrees of freedom
## Residual deviance: 190.33 on 396 degrees of freedom
## AIC: 206.33
##
## Number of Fisher Scoring iterations: 8El porcentaje de aciertos si aumenta respecto a los modelos anteriores, pero se sospecha que es por la inclusión de las variables nox y tax. Se realizará un modelo con todas las variables para ver el contraste.
predicted <- round(predict(model_glm, newdata = boston_vl, type = "response"))
real <- boston_vl$crim01
acc <- Accuracy(predicted, real) * 100
print(sprintf("El porcentaje de aciertos fue del %f%s", acc, "%"))## [1] "El porcentaje de aciertos fue del 83.333333%"ConfusionMatrix(predicted, real) # Matriz de confusión## y_pred
## y_true 0 1
## FALSE 44 9
## TRUE 8 41A continuación se realiza una regresión logística con todas las variables.
model_glm <- glm(crim01 ~ zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + black + lstat + medv, data = boston_tr, family = "binomial")
summary(model_glm)##
## Call:
## glm(formula = crim01 ~ zn + indus + chas + nox + rm + age + dis +
## rad + tax + ptratio + black + lstat + medv, family = "binomial",
## data = boston_tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.8878 -0.1486 0.0000 0.0022 3.2456
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -31.498628 6.802443 -4.630 3.65e-06 ***
## zn -0.071444 0.034115 -2.094 0.03624 *
## indus -0.004212 0.053722 -0.078 0.93750
## chas 0.478901 0.807239 0.593 0.55301
## nox 44.803595 8.223846 5.448 5.09e-08 ***
## rm -0.486381 0.778876 -0.624 0.53232
## age 0.016017 0.012749 1.256 0.20900
## dis 0.569745 0.228907 2.489 0.01281 *
## rad 0.725147 0.172624 4.201 2.66e-05 ***
## tax -0.008332 0.003035 -2.745 0.00605 **
## ptratio 0.271237 0.134621 2.015 0.04392 *
## black -0.008655 0.005567 -1.555 0.12001
## lstat 0.068000 0.057365 1.185 0.23586
## medv 0.174126 0.073687 2.363 0.01813 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 560.02 on 403 degrees of freedom
## Residual deviance: 171.87 on 390 degrees of freedom
## AIC: 199.87
##
## Number of Fisher Scoring iterations: 9El porcentaje de aciertos hasta ahora es el mejor de todos los modelos, sin embargo se realizará un modelo con las variables que parecen tener significancia del modelo.
predicted <- round(predict(model_glm, newdata = boston_vl, type = "response"))
real <- boston_vl$crim01
acc <- Accuracy(predicted, real)* 100
print(sprintf("El porcentaje de aciertos fue del %f%s", acc, "%"))## [1] "El porcentaje de aciertos fue del 87.254902%"ConfusionMatrix(predicted, real) # Matriz de confusión## y_pred
## y_true 0 1
## FALSE 45 8
## TRUE 5 44model_glm <- glm(crim01 ~ zn + nox + dis + rad + tax + ptratio + medv, data = boston_tr, family = "binomial")
summary(model_glm)##
## Call:
## glm(formula = crim01 ~ zn + nox + dis + rad + tax + ptratio +
## medv, family = "binomial", data = boston_tr)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.04249 -0.20026 0.00006 0.00343 3.06622
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -30.508465 5.637643 -5.412 6.25e-08 ***
## zn -0.062811 0.030345 -2.070 0.03846 *
## nox 43.648472 6.933642 6.295 3.07e-10 ***
## dis 0.363455 0.200774 1.810 0.07025 .
## rad 0.682218 0.149559 4.562 5.08e-06 ***
## tax -0.007849 0.002667 -2.943 0.00325 **
## ptratio 0.176106 0.113791 1.548 0.12171
## medv 0.077650 0.033965 2.286 0.02224 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 560.02 on 403 degrees of freedom
## Residual deviance: 181.24 on 396 degrees of freedom
## AIC: 197.24
##
## Number of Fisher Scoring iterations: 9El porcentaje de aciertos no mejora respecto al modelo con todas las variables, entonces se usarán todas las variables como predictores.
predicted <- round(predict(model_glm, newdata = boston_vl, type = "response"))
real <- boston_vl$crim01
acc <- Accuracy(predicted, real)* 100
print(sprintf("El porcentaje de aciertos fue del %f%s", acc, "%"))## [1] "El porcentaje de aciertos fue del 85.294118%"boston_fml <- crim01 ~ zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + black + lstat + medvboston_qda <- qda(boston_fml, data = boston_tr)
boston_qda## Call:
## qda(boston_fml, data = boston_tr)
##
## Prior probabilities of groups:
## FALSE TRUE
## 0.4950495 0.5049505
##
## Group means:
## zn indus chas nox rm age dis
## FALSE 21.770000 6.85325 0.06000000 0.4707995 6.396955 50.81350 5.135899
## TRUE 1.490196 15.20843 0.09313725 0.6369216 6.210858 85.17059 2.496911
## rad tax ptratio black lstat medv
## FALSE 4.1450 306.5650 17.90500 388.4158 9.24765 24.94350
## TRUE 14.7402 506.2353 18.82108 331.0183 15.59593 20.38284predicted <- predict(boston_qda, newdata=boston_vl) # Predicción
real <- boston_vl$crim01 # Reales
c_matrix_qda <- ConfusionMatrix(predicted$class, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos para es de %f%s", Accuracy(predicted$class,real)*100, "%"))## [1] "El porcentaje de aciertos para es de 89.215686%"c_matrix_qda## y_pred
## y_true FALSE TRUE
## FALSE 51 2
## TRUE 9 40boston_lda <- lda(boston_fml, data = boston_tr)
boston_lda## Call:
## lda(boston_fml, data = boston_tr)
##
## Prior probabilities of groups:
## FALSE TRUE
## 0.4950495 0.5049505
##
## Group means:
## zn indus chas nox rm age dis
## FALSE 21.770000 6.85325 0.06000000 0.4707995 6.396955 50.81350 5.135899
## TRUE 1.490196 15.20843 0.09313725 0.6369216 6.210858 85.17059 2.496911
## rad tax ptratio black lstat medv
## FALSE 4.1450 306.5650 17.90500 388.4158 9.24765 24.94350
## TRUE 14.7402 506.2353 18.82108 331.0183 15.59593 20.38284
##
## Coefficients of linear discriminants:
## LD1
## zn -0.0059989846
## indus 0.0354609587
## chas -0.1896438257
## nox 7.8437023277
## rm 0.0522461281
## age 0.0108681241
## dis 0.0339569971
## rad 0.0965538073
## tax -0.0022916846
## ptratio -0.0024566866
## black -0.0006313999
## lstat 0.0172151193
## medv 0.0418635460predicted <- predict(boston_lda, newdata=boston_vl) # Predicción
real <- boston_vl$crim01 # Reales
c_matrix_lda <- ConfusionMatrix(predicted$class, real) # Matriz de confusión
print(sprintf("El porcentaje de aciertos para es de %f%s", Accuracy(predicted$class,real)*100, "%"))## [1] "El porcentaje de aciertos para es de 83.333333%"c_matrix_lda## y_pred
## y_true FALSE TRUE
## FALSE 48 5
## TRUE 12 37# Preparación de los datos
x_tr <- cbind(boston_tr)
x_vl <- cbind(boston_vl)
y_tr <- boston_tr$crim01
y_vl <- boston_vl$crim01Aproximación por KNN
El k que genera el mejor porcentaje de aciertos es \(k=4\).
for (k in c(1:10)) {
set.seed(2020)
predict = knn(train = x_tr, test = x_vl, cl = y_tr, k = k)
print(sprintf("El porcentaje de aciertos para el k=%i es de %f%s", k, Accuracy(predict,y_vl)*100, "%"))
}## [1] "El porcentaje de aciertos para el k=1 es de 93.137255%"
## [1] "El porcentaje de aciertos para el k=2 es de 93.137255%"
## [1] "El porcentaje de aciertos para el k=3 es de 95.098039%"
## [1] "El porcentaje de aciertos para el k=4 es de 96.078431%"
## [1] "El porcentaje de aciertos para el k=5 es de 95.098039%"
## [1] "El porcentaje de aciertos para el k=6 es de 94.117647%"
## [1] "El porcentaje de aciertos para el k=7 es de 94.117647%"
## [1] "El porcentaje de aciertos para el k=8 es de 94.117647%"
## [1] "El porcentaje de aciertos para el k=9 es de 93.137255%"
## [1] "El porcentaje de aciertos para el k=10 es de 94.117647%"Los resultados del KNN con \(k=4\).
set.seed(2020)
predict = knn(train = x_tr, test = x_vl, cl = y_tr, k = 4)
print(sprintf("El porcentaje de aciertos para el k=%i es de %f%s", 4, Accuracy(predict,y_vl)*100,"%"))## [1] "El porcentaje de aciertos para el k=4 es de 96.078431%"table(predict, y_vl)## y_vl
## predict FALSE TRUE
## FALSE 51 2
## TRUE 2 47Conclusión
El modelo que presentó la mejor métrica de porcentaje de aciertos fue KNN con \(k=4\), con un \(96\%\). Las otras aproximaciones como el LDA, GLM, QDA presentan un rendimiento más bajo, yendo desde \(83\%\) hasta \(89\%\). Realizando el análisis de predictores con la regresión logística, se definió que todas las variables aportaban un poco de significancia a la predicción, presentando el mejor rendimiento al usar todas las variables en el modelo.
Apartado 8.4
Punto 7
En el laboratorio, aplicamos bosques aleatorios a los datos de Boston usando \(mtry = 6\) y usando \(ntree = 25\) y \(ntree = 500\). Cree un gráfico que muestre el error de prueba resultante de bosques aleatorios en este conjunto de datos para obtener un rango más completo de valores para \(mtry\) y \(ntree\). Puede modelar su diagrama según la Figura 8.10. Describe los resultados obtenidos.
Figura 8.10 (Libro)
Graficamente se puede observar el resultado de bosques aleatorios para el conjunto de datos Boston (Valores de la vivienda en los suburbios de Boston), donde el error de prueba \((MSE)\) se reduce en la medida que se van agregando más árboles y tiende a estabilizarse alredor de \(80\) árboles. El valor de \(m\) representa el número de predictores disponibles para dividir cada nodo del árbol y se llega a la conclusión que entre menor sea el valor de p, es decir \(m<p\) se tiene un menor \(MSE\) y por lo tanto hay una mejora en el modelo.
Punto 8
En el laboratorio, se aplicó un árbol de clasificación al conjunto de datos de Asientos para automóvil después de convertir Ventas en una variable de respuesta cualitativa. Ahora buscaremos predecir las ventas utilizando árboles de regresión y enfoques relacionados, tratando la respuesta como una variable cuantitativa.
Literal a), División del conjunto
Divida el conjunto de datos en un conjunto de entrenamiento y un conjunto de prueba.
set.seed(123456)
train <- sample(dim(Carseats)[1], dim(Carseats)[1]/2)
train_Carseats <- Carseats[train, ]
test_Carseats <- Carseats[-train, ]Literal b), Ajuste de un árbol de regresión
Ajuste un árbol de regresión al conjunto de entrenamiento. Trace el árbol e interprete los resultados. ¿Qué tasa de error de prueba obtiene?
##
## Regression tree:
## tree(formula = Sales ~ ., data = train_Carseats)
## Variables actually used in tree construction:
## [1] "ShelveLoc" "Price" "Age" "Income" "CompPrice"
## [6] "Education" "Population" "US"
## Number of terminal nodes: 16
## Residual mean deviance: 2.092 = 385 / 184
## Distribution of residuals:
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -3.55400 -1.03000 -0.01708 0.00000 0.94840 3.66700Del resumen del árbol de regresión ajustado se observa que de las 10 covariables solo \(8\) han sido utilizadas en la construcción del árbol, este árbol proporcioda una suma de cuadrados de los errores de \(184\).
El árbol de regresión contiene \(15\) nodos intermedios y \(16\) nodos terminales, se observa que el nivel Good del factor ShelveLoc no es de importancia y para la división solo se tienen en cuenta los niveles Bad y Medium, cuando este factor esta en Bad se encuentra otro nodo interno, si Price es menor a \(123.5\) se despliegan otros dos nodos intermedios, el de la izquierda relacionado con la variable Age y el de la derecha con CompPrice; si ShelveLoc es Medium se encuentra un nodo interno relacionado con la variable Price, si esta es menor a \(109.5\) se despliegan dos nodos intermedios que involucrans la variable CompPrice y US, cuando esta última es No la predicción para Sales es \(7.499\). La predicción para Sales cuando ShelveLoc es Bad, Price es menor a \(105.5\), Age es menor a \(50\), Price es menor a \(80.5\) es \(9.55\), de manera similar se interpretan las otras predicciones para la variable respuesta Sales.
set.seed(1234)
pred_Carseats <- predict(tree_Carseats, test_Carseats)
paste("El error de prueba es:", mean((test_Carseats$Sales - pred_Carseats)^2))## [1] "El error de prueba es: 5.0324577229447"Literar c), Determinar el nivel óptimo del árbol
Utilice la validación cruzada para determinar el nivel óptimo de complejidad del árbol. ¿La poda del árbol mejora la tasa de error de prueba?
set.seed(1234)
cv_Carseats <- cv.tree(tree_Carseats, FUN = prune.tree)
El nivel óptimo de complejidad del árbol es de \(8\) nodos terminales, a continuación, se presenta el ajuste del árbol de regresión con \(8\) nodos terminales.
p_Carseats <- prune.tree(tree_Carseats, best = 8)
ppruned <- predict(p_Carseats, test_Carseats)
paste("El error de prueba es:", mean((test_Carseats$Sales - ppruned)^2))## [1] "El error de prueba es: 4.87611548961474"esta tsa es menor a la obtenida con el árbol sin podar, es decir que la poda del árbol mejoro el ajuste.
Literal c), Importancia de variables
Utilice el método de \(bagging\) para analizar estos datos. Qué tasa de error de prueba que obtiene? Utilice la función de \(importance()\) para determinar qué variables son más importantes.
set.seed(1234)
embolsado_Carseats <- randomForest(Sales ~ ., data = train_Carseats,
mtry = 10, ntree = 500, importance = T)
embolsado_pred <- predict(embolsado_Carseats, test_Carseats)
paste("El error de prueba es:",
mean((test_Carseats$Sales - embolsado_pred)^2))## [1] "El error de prueba es: 2.98353481792896"esta tasa de error es mucho menor que el obtenido en los literales anteriores.
importance(embolsado_Carseats)## %IncMSE IncNodePurity
## CompPrice 24.7104427 181.982752
## Income 10.3732873 115.215011
## Advertising 16.3417929 109.833777
## Population -2.1288904 51.722072
## Price 48.1234357 356.453668
## ShelveLoc 60.2725206 457.515905
## Age 16.0661998 151.183964
## Education 0.7957303 36.296672
## Urban -0.3120038 7.093213
## US 7.1717362 16.766946Se observa que cuando la Variable ShelveLoc es retirada del modelo el \(MSE\) se incrementa en un \(60.27\%\) y cuando se retira Price dicho \(MSE\) aumenta en \(48.12\%\).
Los resultados indican que en todos los árboles considerados en el bosque aleatorio, el precio que cobra la empresa por los asientos de seguridad en cada sitio (Price) y el que indica la calidad de la ubicación de las estanterías para los asientos del automóvil en cada sitio (ShelveLoc) son las dos variables más importantes.
Literal e), Usar bosques aleatorios
Utilice bosques aleatorios para analizar estos datos. ¿Qué tasa de error de prueba obtiene? Utilice la función de \(importance()\) para determinar qué variables son más importantes. Describa el efecto de \(m\), el número de variables consideradas en cada división, sobre la tasa de error obtenida.
Para el ajuste se usa \(m=8\) ya que es el valor que disminuye el error de prueba.
bosques_Carseats <- randomForest(Sales ~ ., data = train_Carseats,
mtry = 8, ntree = 500, importance = T)
bosques_pred <- predict(bosques_Carseats, test_Carseats)
paste("El error de prueba es:",
mean((test_Carseats$Sales - bosques_pred)^2))## [1] "El error de prueba es: 3.01317456599399"un poco mayor a la obtenida cuando se usa el método de \(bagging\).
importance(bosques_Carseats)## %IncMSE IncNodePurity
## CompPrice 24.0654029 177.494697
## Income 9.3877917 117.014835
## Advertising 14.5452522 106.266632
## Population 0.2175140 56.304910
## Price 44.9379411 354.508881
## ShelveLoc 62.0919159 430.398167
## Age 17.5754922 152.541774
## Education -0.2325020 40.519237
## Urban 0.3138798 7.102657
## US 7.2092058 19.184242
Se observa que cuando la Variable ShelveLoc es retirada del modelo el \(MSE\) se incrementa en un \(57.76\%\) y cuando se retira Price dicho \(MSE\) aumenta en \(47.20\%\).
Punto 9
Este problema involucra al conjunto de datos de OJ que hace parte del paquete ISLR.
Literal a) División del conjunto de datos
Crear un conjunto de entrenamiento que contenga una muestra aleatoria de \(800\) observaciones, y un conjunto de pruebas que contiene las observaciones restantes.
set.seed(1234)
train <- sample(dim(OJ)[1], 800)
OJtrain <- OJ[train, ]
OJtest <- OJ[-train, ]Literal b), Ajustar un árbol
Ajustar un árbol a los datos de entrenamiento, con la respuesta “Purchase”. y las otras variables como predictores. Utilice la función \(summary()\) para producir estadísticas resumidas sobre el árbol, y describir la resultados obtenidos. ¿Cuál es la tasa de error de entrenamiento? ¿Cuántos nodos terminales que tiene el árbol?
set.seed(1234)
OJtree <- tree(Purchase ~ ., data = OJtrain)
summary(OJtree)##
## Classification tree:
## tree(formula = Purchase ~ ., data = OJtrain)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 7
## Residual mean deviance: 0.7482 = 593.3 / 793
## Misclassification error rate: 0.165 = 132 / 800Las estadísticas resumidas obtenidas por la función \(summary\) se tiene que el árbol usa \(2\) variables predictoras, las cuales son LoyalCH y PriceDiff. El árbol tiene \(7\) nodos terminales y su tasa de error es de \(0.165\).
Literal c)
Escriba el nombre del objeto de árbol para obtener una salida de texto detallada. Elija uno de los nodos terminales e interprete la información que se muestra.
Se escogió el nodo terminal “3)”. La variable dividida en este nodo es LoyalCH, el valor de division en este nodo es de \(0.48285\) y hay \(498\) puntos en el subarbol debajo de este nodo. La desviación para todos los puntos contenidos en la región por debajo de este nodo es de \(442.3\). La predicción de este nodo es \(Sales=CH\). Alrededor de \(83.735\%\) puntos en este nodo tienen a CH como valor de Sales, el otro \(16.265\%\) restante tienen a MM como valor de Sales.
Literal d)
Crear un gráfico del árbol e interpretar los resultados.
La variable mas importantes del árbol es LoyalCH, estando en los primeros tres nodos del árbol. Si la variable LoyalCH es menor que \(0.26\) el árbol predice MM, si LoyalCH es mayor a \(0.76\) el árbol predice CH.
Literal e)
Predecir la respuesta en los datos de la prueba, y producir una matriz de confusión que compara las etiquetas de prueba con las etiquetas de prueba predichas. ¿Cuál es la tasa de error de la prueba?
set.seed(1234)
OJpred <- predict(OJtree, OJtest, type = "class")
table(OJtest$Purchase, OJpred)## OJpred
## CH MM
## CH 126 44
## MM 17 83La tasa de error de prueba es de \(21.11%\).
Literal f)
Aplicar la función \(cv.tree()\) al conjunto de entrenamiento para determinar el tamaño óptimo del árbol.
OJcv = cv.tree(OJtree, FUN = prune.tree)Literal g)
Elaborar un gráfico con el tamaño del árbol en el eje x y validarlo de forma cruzada con la tasa de error de clasificación en el eje Y.
plot(OJcv$size, OJcv$dev, type = "b", xlab = "Tamaño del árbol", ylab = "Desviación",las=1)
Literal h)
¿Qué tamaño de árbol corresponde a la clasificación validada cruzada más baja tasa de error?
Un árbol con \(5\) nodos terminales arrojaría la menor tasa de error de validación cruzada.
Literal i)
Producir un árbol podado que corresponda al tamaño óptimo del árbol obtenido mediante validación cruzada. Si la validación cruzada no conduce a la selección de un árbol podado, luego crear un árbol podado con cinco nodos terminales.
set.seed(1234)
p_OJ <- prune.tree(OJtree, best = 5)
summary(p_OJ)##
## Classification tree:
## snip.tree(tree = OJtree, nodes = 5L)
## Variables actually used in tree construction:
## [1] "LoyalCH" "PriceDiff"
## Number of terminal nodes: 5
## Residual mean deviance: 0.7755 = 616.5 / 795
## Misclassification error rate: 0.165 = 132 / 800Las estadística resumidas obtenidas por la función \(summary\) se tiene que el árbol usa \(2\) variables predictoras, las cuales son LoyalCH y PriceDiff El árbol tiene \(5\) nodos terminales y su tasa de error es de \(0.1825\).
Literal j)
Comparar las tasas de error de entrenamiento entre los podados y los no podados árboles. ¿Cuál es más alto?
Al hacer el ajuste de los árboles se obtiene que la tasa de error de entrenamiento para el árbol sin podar es \(14.75%\) y para el árbol podado es \(18.25\) que es levemente mayor, por lo tanto la poda del árbol no mejor significativamente el error.
Literal k)
Comparar los índices de error de la prueba entre los podados y los no podados árboles. ¿Cuál es más alto?
## [1] "El MSE de prueba para el árbol no podado es 22.593"## [1] "El MSE de prueba para el árbol podado es 22.593"En cuanto a error de prueba es mejor usar el árbol podado para realizar las predicciones.
Punto 10
Para el punto número 10, se hará uso del método boosting, para predicción de salarios, en el set de datos Hitters.
Tratamiento de datos
Primero se eliminan las observaciones, para las que la información asociada a la variable salary es desconocida, y luego se aplica una transformación logarítmica.
Hitters <- na.omit(Hitters) #remoción de NA's.
Hitters$Salary <- log(Hitters$Salary) #Transformación logarítmicaSe generará un set de entrenamiento, con las primeras \(200\) observaciones y un set de prueba con las observaciones restantes.
set.seed(0)
filas <- nrow(Hitters)
variables <- ncol(Hitters) - 1 #Se omite la variable respuesta
train <- 1:200
test <- 201:filas
train_set <- Hitters[train, ]
test_set <- Hitters[test, ]Ahora, se generará una lista con un rango de parámetros lambda, con el fin de calcular el \(MSE\) asociado para cada lambda, tanto para los datos de entrenamiento, como de prueba. Se usará un modelo \(boosting\), con una distribución gaussiana, debido a que la variable respuesta es de tipo numérico.
set.seed(0)
lambda_set <- seq(1e-03, 0.05, by = 0.001) #Shrinkage lambdas.
#Listas de NA's con la longitud del lambda_set, para almacenar el MSE de training y testing
train_set_mse <- rep(NA, length(lambda_set))
test_set_mse <- rep(NA, length(lambda_set))
# Iteración y guardado de los MSE_ test y MSE_train
for (i in 1:length(lambda_set)) {
ld <- lambda_set[i]
boost.hitters <- gbm(Salary ~., data = train_set, distribution = "gaussian", n.trees =1000, interaction.depth = 4, shrinkage = ld)
y_hat <- predict(boost.hitters, newdata = train_set, n.trees = 1000)
#MSE_train
train_set_mse[i] <- mean((y_hat - train_set$Salary)^2)
y_hat <- predict(boost.hitters, newdata = test_set, n.trees = 1000)
#MSE_test
test_set_mse[i] <- mean((y_hat - test_set$Salary)^2)
}
#Gráfica MSE vs. Lambda
plot(lambda_set, train_set_mse, type = "b", pch = 19, col = "black", xlab = "Lambda Value", ylab = "train_MSE")
plot(lambda_set, test_set_mse, type = "b", pch = 20, col = "blue", xlab = "Lambda Value", ylab = "Test Set MSE")
grid()
Se realizará un modelo de regresión lineal con este set de datos, con el fin de realizar un comparativo con el \(MSE\) del modelo \(boosting\).
lr <- lm(Salary ~., data = train_set)
y_hat_tr <- predict(lr, newdata = train_set)
y_hat_t <- predict(lr, newdata = test_set)
MSE_train <- mean((y_hat_tr - train_set$Salary)^2)
MSE_test <- mean((y_hat_t - test_set$Salary)^2)
sprintf("MSE_train: %f, MSE_test: %f", MSE_train, MSE_test)## [1] "MSE_train: 0.320373, MSE_test: 0.491796"Se puede observar, que el \(MSE\), tanto en los datos de prueba, como en los datos de entrenamiento sigue siendo menor para el modelo \(boosting\), el cual posee un \(MSE_test\) de \(0.01086731\).
Por último se realiza un modelo de regresión tipo Lasso.
lss <- model.matrix(Salary ~., data = train_set)
lss2 <- model.matrix(Salary ~., data = test_set)
cv.out <- cv.glmnet(lss, train_set$Salary, alpha = 1)
bestld <- cv.out$lambda.1se
y_hat_tr <- predict(cv.out, s = bestld, newx = lss)
y_hat_t <- predict(cv.out, s = bestld,
newx = lss2)
MSE_train <- mean((y_hat_tr - train_set$Salary)^2)
MSE_test <- mean((y_hat_t - test_set$Salary)^2)
sprintf("lasso CV best value of lambda (one standard error): %f", bestld)## [1] "lasso CV best value of lambda (one standard error): 0.172520"sprintf("lasso regression MSE_train: %f", MSE_train)## [1] "lasso regression MSE_train: 0.431836"sprintf("lasso regression MSE_train: %f", MSE_test)## [1] "lasso regression MSE_train: 0.440617"En este caso, nuevamente sigue teniendo un \(MSE\) más bajo el modelo \(boosting\).
Ahora se identifican las variables más importantes para el modelo \(boosted\).
summary(boost.hitters)
## var rel.inf
## CAtBat CAtBat 27.0020209
## CWalks CWalks 10.6163551
## CHits CHits 7.6902661
## Walks Walks 6.8483932
## PutOuts PutOuts 6.1286755
## CRBI CRBI 6.0907884
## Years Years 4.4675957
## CRuns CRuns 4.4102640
## Assists Assists 3.7882081
## RBI RBI 3.7876623
## AtBat AtBat 3.7557979
## CHmRun CHmRun 3.5852874
## Runs Runs 2.8328845
## HmRun HmRun 2.7786882
## Hits Hits 2.7428142
## Errors Errors 2.1044565
## NewLeague NewLeague 0.7478656
## Division Division 0.4581761
## League League 0.1638003Es claro que las variables más importantes son CAtBat y CWalks.
Generación Modelo bagging.
bag.hitters <- randomForest(Salary ~., data = train_set, mtry = variables, ntree = 1000, importance = TRUE)
y_hat <- predict(bag.hitters, newdata = test_set)
mse <- mean((test_set$Salary - y_hat)^2)
sprintf("Bagging MSE_test: %f", MSE_test)## [1] "Bagging MSE_test: 0.440617"En este caso, el argumento \(mtyr=variables\), nos indica que todos los predictores deberían ser considerados en cada una de las divisiones del arbol, en otras palabras, que se debe hacer uso del método \(bagging\).
Punto 11
Para la realización de este punto, se hace uso del set de datos Caravan
En principio se creará un set de datos, donde las primeras \(1000\) observaciones serán consideradas como el set de entrenamiento y las restantes como el set de prueba.
set.seed(0)
df <- read.csv("data.csv")
#Se omiten los NA's
df <- na.omit(df)
filas <- nrow(df)
variables <- ncol(df) - 1 #Se omite la variable respuesta
train <- 1:1000
test <- 1001:filas
#Estas variables no presentan variación en el modelo, por lo que son eliminadas ("variable 50: PVRAAUT has no variation.variable 71: AVRAAUT has no variation").
df$PVRAAUT <- NULL
df$AVRAAUT <- NULL
#La variable 'Purchase' se transformará en una variable binaria (0, 1), para hacer uso del boosting model.
Purchase <- ifelse(df$Purchase == "Yes", 1, 0)
df$Purchase <- Purchase
train_set <- df[train, ]
test_set <- df[test, ]Para abordar este problema, se realiza un modelo tipo \(boosting\) haciendo uso del set de entrenamineto, la variable respuesta será \(Purchase\), se hará uso de \(1000\) árboles y un valor \(\lambda=0.01\). En este caso se hará uso de la distribución Bernoulli, debido a que la variable respuesta es binaria.
#Entrenamiento del modelo.
ld <- 0.01 #Shrinkage lambda
boost.caravan <- gbm(Purchase ~., data = train_set, distribution = "bernoulli", n.trees = 1000, interaction.depth = 2, shrinkage = ld)
summary(boost.caravan)
## var rel.inf
## PPERSAUT PPERSAUT 9.83047713
## MKOOPKLA MKOOPKLA 5.67101078
## MOPLHOOG MOPLHOOG 5.57015086
## PBRAND PBRAND 5.53090647
## MGODGE MGODGE 5.50146339
## MBERMIDD MBERMIDD 4.66053217
## MOSTYPE MOSTYPE 3.71528842
## MINK3045 MINK3045 3.29167392
## MGODPR MGODPR 3.15189659
## MAUT2 MAUT2 2.75860820
## MBERARBG MBERARBG 2.41585129
## ABRAND ABRAND 2.34750399
## MSKC MSKC 2.07742330
## MSKA MSKA 2.03539655
## MAUT1 MAUT1 2.00290825
## MRELGE MRELGE 1.91903070
## MSKB1 MSKB1 1.88713704
## PWAPART PWAPART 1.88354710
## MFWEKIND MFWEKIND 1.86568208
## MGODOV MGODOV 1.75362888
## MINK7512 MINK7512 1.68185238
## MFGEKIND MFGEKIND 1.62053883
## MBERHOOG MBERHOOG 1.47386158
## MBERARBO MBERARBO 1.45496977
## MINKM30 MINKM30 1.45332376
## MHKOOP MHKOOP 1.44894919
## MGODRK MGODRK 1.42859782
## MRELOV MRELOV 1.35136073
## MINKGEM MINKGEM 1.33514219
## MAUT0 MAUT0 1.20352717
## MZFONDS MZFONDS 1.17292916
## MINK4575 MINK4575 1.16903944
## MOSHOOFD MOSHOOFD 1.13985697
## MHHUUR MHHUUR 1.04951077
## MGEMLEEF MGEMLEEF 1.03005616
## APERSAUT APERSAUT 0.99864198
## MSKB2 MSKB2 0.97672510
## PBYSTAND PBYSTAND 0.87858256
## MOPLMIDD MOPLMIDD 0.87431940
## PMOTSCO PMOTSCO 0.83289524
## MFALLEEN MFALLEEN 0.82850043
## MZPART MZPART 0.82588303
## PLEVEN PLEVEN 0.72568141
## MGEMOMV MGEMOMV 0.70719044
## MSKD MSKD 0.69605715
## MBERBOER MBERBOER 0.63402319
## MBERZELF MBERZELF 0.38036763
## MRELSA MRELSA 0.32946406
## MINK123M MINK123M 0.21200829
## MOPLLAAG MOPLLAAG 0.19163868
## MAANTHUI MAANTHUI 0.02438834
## PWABEDR PWABEDR 0.00000000
## PWALAND PWALAND 0.00000000
## PBESAUT PBESAUT 0.00000000
## PAANHANG PAANHANG 0.00000000
## PTRACTOR PTRACTOR 0.00000000
## PWERKT PWERKT 0.00000000
## PBROM PBROM 0.00000000
## PPERSONG PPERSONG 0.00000000
## PGEZONG PGEZONG 0.00000000
## PWAOREG PWAOREG 0.00000000
## PZEILPL PZEILPL 0.00000000
## PPLEZIER PPLEZIER 0.00000000
## PFIETS PFIETS 0.00000000
## PINBOED PINBOED 0.00000000
## AWAPART AWAPART 0.00000000
## AWABEDR AWABEDR 0.00000000
## AWALAND AWALAND 0.00000000
## ABESAUT ABESAUT 0.00000000
## AMOTSCO AMOTSCO 0.00000000
## AAANHANG AAANHANG 0.00000000
## ATRACTOR ATRACTOR 0.00000000
## AWERKT AWERKT 0.00000000
## ABROM ABROM 0.00000000
## ALEVEN ALEVEN 0.00000000
## APERSONG APERSONG 0.00000000
## AGEZONG AGEZONG 0.00000000
## AWAOREG AWAOREG 0.00000000
## AZEILPL AZEILPL 0.00000000
## APLEZIER APLEZIER 0.00000000
## AFIETS AFIETS 0.00000000
## AINBOED AINBOED 0.00000000
## ABYSTAND ABYSTAND 0.00000000Es claro que para este modelo, las variables más importantes son PPERSAUT y MKOOPKLA.
# Predicción del error de prueba
y_hat <- predict(boost.caravan, newdata = test_set, n.trees = 1000)
p_hat <- exp(y_hat)/(1 + exp(y_hat))
# Conversión de los odd logarítmicos en probabilidades
si_compra <- rep(0, length(test))
si_compra[p_hat > 0.2]<-1
# Creación de matriz de confusión.
table(si_compra, test_set$Purchase)##
## si_compra 0 1
## 0 4357 253
## 1 176 36Ahora se genera un modelo de regresión lógistica, con el fin de generar un comparativo.
lg <- glm(Purchase ~., data = train_set, family = "binomial")
#Predicciones
y_hat <- predict(lg, newdata = test_set)
#Conversión de probabiliades
p_hat <- exp(y_hat)/(1 + exp(y_hat))
si_compra <- rep(0, length(test))
si_compra[p_hat > 0.2] <- 1
#Matriz de confusión
table(si_compra, test_set$Purchase)##
## si_compra 0 1
## 0 4183 231
## 1 350 58El modelo de regresión lógistica presenta una mayor cantidad de observaciones mal clasificadas.
Punto 12
Aplicar el potenciamiento, el embolsamiento y los bosques aleatorios a un conjunto de datos de su elección. Asegúrate de que los modelos encajen en datos de entrenamiento y de que evalúen su rendimiento en datos de prueba. ¿Cómo de precisos son los resultados comparados a métodos simples como la regresión lineal o logística? ¿Cuál de estos que los enfoques de la investigación dan el mejor resultado?
Por comodidad se usará la base de datos Caravan.
train <- 1:1000
library(tidyverse)
Caravan <- Caravan%>%
mutate(Purchase=case_when(Purchase == "Yes" ~ 1, TRUE~0))
Ctrain <- Caravan[train, ]
Ctest <- Caravan[-train, ]#Usando regresión logistica
logist <- glm(Purchase ~ ., data = Ctrain, family = binomial)
logistp <- predict(logist, Ctest, type = "response")
prediccion <- ifelse(logistp > 0.2, 1, 0)
table(Ctest$Purchase, prediccion)## prediccion
## 0 1
## 0 4183 350
## 1 231 58El error de prueba es de \(12.04\%\).
#Usando potenciamiento
set.seed(1234)
library(gbm)
cpotencido <- gbm(Purchase ~ ., data = Ctrain, n.trees = 1000,
shrinkage = 0.01,
distribution = "bernoulli")
potenciado <- predict(cpotencido, Ctest, n.trees = 1000, type = "response")
potenciadop <- ifelse(potenciado > 0.2, 1, 0) #mayor a 20%
table(Ctest$Purchase, potenciadop)## potenciadop
## 0 1
## 0 4406 127
## 1 258 31Test de error es de \(7.98\%\).
#Usando embolsamiento
set.seed(1234)
cembolsado <- randomForest(Purchase ~ ., data = Ctrain, mtry = 6)
cpredembolsado <- predict(cembolsado, newdata = Ctest)
embolsadop <- ifelse(cpredembolsado > 0.2, 1, 0) #mayor a 20%
table(embolsadop, Ctest$Purchase)##
## embolsadop 0 1
## 0 4310 249
## 1 223 40Test de error es de \(9.78\%\).
#Usando bosques aleatorios
set.seed(1234)
cbosques <- randomForest(Purchase ~ ., data = Ctrain, mtry = 2)
cbosquesp <- predict(cbosques, newdata = Ctest)
bp <- ifelse(cbosquesp > 0.2, 1, 0) #mayor a 20%
table(bp, Ctest$Purchase)##
## bp 0 1
## 0 4447 270
## 1 86 19Test de error es de \(7.38\%\)
El error de test más bajo fue el de bosques aleatorios.
Apartado 9.7
Punto 4
El objetivo de este punto es crear dos conjuntos de puntos (clases) no lineales para comparar el rendimiento de las máquinas de soporte vectorial con un kernel diferente.
Se crean dos conjuntos de puntos y se grafican para identificar si son no lineales.
set.seed(69)
x <- rnorm(100)
y <- 4 * x^2 + 0 + rnorm(100)
set <- sample(100, 50)
y[set] <- y[set] + 3
y[-set] <- y[-set] - 3
plot(x[set], y[set], col = "red", xlab = "X", ylab = "Y", ylim = c(-5, 25))
points(x[-set], y[-set], col = "blue")
Se crean los datasets de entrenamiento y de prueba.
z <- rep(0, 100)
z[set] <- 1
data <- data.frame(x = x, y = y, z = as.factor(z))
train <- sample(100, 50)
data.train <- data[train, ]
data.test <- data[-train, ]Entrenamos la maquina de soporte vectorial con un kernel linear
svm.linear <- svm(z ~ ., data = data.train, kernel = "linear", cost = 10)
plot(svm.linear, data.train)
Se visualizan los errores
table(predict = predict(svm.linear, data.train), truth = data.train$z)## truth
## predict 0 1
## 0 19 1
## 1 5 25Se entrena a la máquina de soporte vectorial con un kernel radial.
svm.radial <- svm(z ~ ., data = data.train, kernel = "radial", gamma = 1, cost = 10)
plot(svm.radial, data.train)
Se visualizan los errores
table(predict = predict(svm.radial, data.train), truth = data.train$z)## truth
## predict 0 1
## 0 24 0
## 1 0 26El entrenamiento con un kernel radial tiene menos errores que el entrenamiento con kernel lineal por lo tanto tiene mejor rendimiento.
Punto 5
El objetivo de este punto es aprender cómo obtener límites de decisión para clasificación a partir de una regresión logística aplicando transformaciones no lineales
Literal a), Creacion del dataset
Se crean los puntos (como nos indica el libro)
x1 = runif (10000) -0.5
x2 = runif (10000) -0.5
y =1*( x1 ^2 - x2 ^2 > 0)Literal b), Graficar el dataset
Se realiza el gráfico del dataset creado anteriormente usando diferentes colores para cada clase.
plot(x1[y == 0], x2[y == 0], col = "red", xlab = "X1", ylab = "X2", pch = "+")
points(x1[y == 1], x2[y == 1], col = "blue", pch = "x")
Literal c), Regresión logística
Se realiza el ajuste de una regresión logística al dataset.
logit.fit <- glm(y ~ x1 + x2, family = "binomial")
summary(logit.fit)##
## Call:
## glm(formula = y ~ x1 + x2, family = "binomial")
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.180 -1.168 -1.157 1.187 1.199
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.02245 0.02000 -1.122 0.262
## x1 -0.02213 0.06978 -0.317 0.751
## x2 0.03735 0.06896 0.542 0.588
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 13862 on 9999 degrees of freedom
## Residual deviance: 13861 on 9997 degrees of freedom
## AIC: 13867
##
## Number of Fisher Scoring iterations: 3Las variables \(X1\) y \(X2\) no son siginificativas para el modelo ajustado.
Literal d), Aplicación del modelo
Se aplica el modelo y se realiza un gráfico con el fin de visualizar de una manera mas clara.
data = data.frame(x1 = x1, x2 = x2, y = y)
lm.prob = predict(logit.fit, data, type = "response")
lm.pred = ifelse(lm.prob > 0.5, 1, 0)
data.pos = data[lm.pred == 1, ]
data.neg = data[lm.pred == 0, ]
plot(data.pos$x1, data.pos$x2, col = "blue", xlab = "X1", ylab = "X2", pch = "+",xlim = c(-0.5, 0.5), ylim = c(-0.5, 0.5))
points(data.neg$x1, data.neg$x2, col = "red", pch = "x")
Como se esperaba el límite de decisión es lineal y no se ajusta correctamente a los datos.
Literal e), Regresión logística con transformaciones no lineales
Ahora se ajusta una regresión logística, pero usando funciones no lineales como “predictors”.
logitnl.fit <- glm(y ~ poly(x1, 2) + poly(x2, 2) + I(x1 * x2), data = data, family = "binomial")
summary(logitnl.fit)##
## Call:
## glm(formula = y ~ poly(x1, 2) + poly(x2, 2) + I(x1 * x2), family = "binomial",
## data = data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.4171 0.0000 0.0000 0.0000 0.3745
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1569 1173 -1.337 0.181
## poly(x1, 2)1 -72032 53319 -1.351 0.177
## poly(x1, 2)2 6015847 4300776 1.399 0.162
## poly(x2, 2)1 38052 29447 1.292 0.196
## poly(x2, 2)2 -6040513 4322410 -1.397 0.162
## I(x1 * x2) -235 1465 -0.160 0.873
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1.3862e+04 on 9999 degrees of freedom
## Residual deviance: 4.5154e-01 on 9994 degrees of freedom
## AIC: 12.452
##
## Number of Fisher Scoring iterations: 25Ninguna de las fuunciones no lineales usadas logro agregar significancia al modelo.
Literal f), Aplicación del nuevo modelo con trasformaciones no lineales
Se aplica el nuevo modelo y se realiza un gráfico con el fin de visualizar de una manera mas clara.
lm.prob = predict(logitnl.fit, data, type = "response")
lm.pred = ifelse(lm.prob > 0.5, 1, 0)
data.pos = data[lm.pred == 1, ]
data.neg = data[lm.pred == 0, ]
plot(data.pos$x1, data.pos$x2, col = "blue", xlab = "X1", ylab = "X2", pch = "+",xlim = c(-0.5, 0.5), ylim = c(-0.5, 0.5))
points(data.neg$x1, data.neg$x2, col = "red", pch = "x")
Ahora el límite de decisión es no lineal y se ajusta mejor a los datos
Literal g), Ajuste de un clasificador vectorial
Para este caso se utilizará una máquina de soporte vectorial con kernel linear y se grafican los resultados.
svm.fit = svm(as.factor(y) ~ x1 + x2, data, kernel = "linear", cost = 10)
svm.pred = predict(svm.fit, data)
data.pos = data[svm.pred == 1, ]
data.neg = data[svm.pred == 0, ]
plot(data.pos$x1, data.pos$x2, col = "blue", xlab = "X1", ylab = "X2", pch = "+",xlim = c(-0.5, 0.5), ylim = c(-0.5, 0.5))
points(data.neg$x1, data.neg$x2, col = "red", pch = "x")
De igual forma que con la regresión logística el límite de decisión es lineal y no se ajusta correctamente a los datos.
Literal h), Máquina de soporte vectorial con kernel no-linear
Se ajusta un máquina de soporte vectorial con kernel radial.
svm.fit = svm(as.factor(y) ~ x1 + x2, data, kernel="radial", gamma = 1, cost=10)
svm.pred = predict(svm.fit, data)
data.pos = data[svm.pred == 1, ]
data.neg = data[svm.pred == 0, ]
plot(data.pos$x1, data.pos$x2, col = "blue", xlab = "X1", ylab = "X2", pch = "+",xlim = c(-0.5, 0.5), ylim = c(-0.5, 0.5))
points(data.neg$x1, data.neg$x2, col = "red", pch = 4)
Literal i), Resultados
En este punto se identifica que si se usa una máquina de soporte vectorial con kernel no-lineal o una regresión logística con transformaciones puede ser muy útil para encontrar límites de decisión no lineales, en cambio si se usa un kernel lineal o no se usan transformaciones los resultados están muy alejados de la realidad para este tipo de problemas.
Punto 6
Al final de la Sección \(9.6.1\), se afirma que en el caso de datos que son apenas separables linealmente, un clasificador de vectores de soporte con un valor pequeño de \(cost\) que clasifica erróneamente un par de observaciones de entrenamiento puede funcionar mejor en datos de prueba que uno con un valor de \(cost\) enorme que no clasifica erróneamente ninguna observación de entrenamiento. Ahora investigará este reclamo.
Literal a), simulación de datos
Genere datos de dos clases con \(p = 2\) de tal manera que las clases sean apenas separables linealmente.
Se inicia fijando una semilla, con el fin de obtener resultados iguales si el ejericio es simulado por otra persona, se simularan dos variables de una \(N(0,1)\), se crea la variable \(y\) que contiene las verdaderas clases a las que pertenecen las observaciones, se hace un cambio en la media para garantizar que las clases sean realmente separables.
set.seed (123456)
x <- matrix (rnorm (1000*2) , ncol =2)
y <- c(rep (-1,50) , rep (1 ,50) )
x[y==1 ,] <- x[y==1,] + 3
Observe que las clases son apenas linealmente separables.
Literal b), Clasificador de soporte vectorial.
Calcule las tasas de error de validación cruzada para clasificadores de vectores de soporte con un rango de valores de costo. ¿Cuántos valores de entrenamiento se clasifican erróneamente para cada valor de costo considerado y cómo se relaciona esto con los errores de validación cruzada obtenidos?
Es necesario crear un data frame con la variable respuesta codificada como factor y generar la división del conjunto en entrenamiento y prueba, a demás se cambiaran los nombres de las columnas para que sean más claros y finalmente se procede con el ajuste del clasificador de soporte vectorial con diferentes valores del parámetro \(Cost\), para esto se hace uso de la función \(tune()\) de la librería \(e1071\).
#Creación y corrección de nombres de columnas
data <- data.frame(x, as.factor(y))
colnames(data) <- c("X1", "X2", "Y")
#División del conjunto de datos
set.seed(123456)
index <- sample(nrow(data), 700)
train <- data[index, ]
test <- data[-index, ]
#Ajuste del clasificador lineal
tune.out <- tune(svm, Y ~ ., data = train,
kernel = "linear",
ranges = list(cost = c(0.001, 0.01, 0.1, 1, 5, 10, 15, 20)),
scale = TRUE)
Como es necesario identificar cuantas observaciones del conjunto de entrenamiento están mal clasificadas, se debe ajustar un calsificador de soporte vectorial para cada valor considerado del parámetro \(cost\), realizar las predicciones y calcular su tasa de error de mala clasificación.
\(Cost=0.001\)
svmfit0.001 <- svm(Y~., data = train, kernel ="linear", cost = 0.001,
scale = T)
predict0.001 <- predict(svmfit0.001, train)
table(Predict = predict0.001 , Truth = train$Y)## Truth
## Predict -1 1
## -1 343 9
## 1 5 343## [1] "Error de entrenamiento: 2 %"La matriz anterior permite concluir que con valor del parámetro \(cost=0.001\), \(686\) de las observaciones se han clasificado de forma correcta, las restantes \(14\) de las \(700\) observaciones del conjunto de entrenamiento quedaron mal clasificadas y por lo tanto el \(MSE\) de entrenamiento es de \(2\%\).
\(Cost=0.01\)
svmfit0.01 <- svm(Y~., data = train, kernel ="linear", cost = 0.01,
scale = T)
predict0.01 <- predict(svmfit0.01, train)
table(Predict = predict0.01 , Truth= train$Y)## Truth
## Predict -1 1
## -1 342 7
## 1 6 345## [1] "Error de entrenamiento: 1.857 %"La matriz anterior permite concluir que con valor el parámetro \(cost=0.01\), \(687\) de las observaciones se han clasificado de forma correcta, solo \(13\) de las \(700\) observaciones del conjunto de entrenamiento quedaron mal clasificadas y por lo tanto el \(MSE\) de entrenamiento es del \(1.857\%\).
\(Cost=0.1\)
svmfit0.1 <- svm(Y~., data = train, kernel ="linear", cost = 0.1,
scale = T)
predict0.1 <- predict(svmfit0.1, train)
table(Predict = predict0.1 , Truth= train$Y)## Truth
## Predict -1 1
## -1 341 5
## 1 7 347## [1] "Error de entrenamiento: 1.714 %"\(Cost=1\)
svmfit1 <- svm(Y~., data = train, kernel ="linear", cost = 1,
scale = T)
predict1 <- predict(svmfit1, train)
table(Predict = predict1 , Truth= train$Y)## Truth
## Predict -1 1
## -1 341 5
## 1 7 347## [1] "Error de entrenamiento: 1.714 %"Observe que para los valores de \(cost=0.1\) y \(cost=1\) las tasas de mala clasificación son iguales, ambas son pequeñas.
\(Cost=5\)
svmfit5 <- svm(Y~., data = train, kernel ="linear", cost = 5,
scale = T)
predict5 <- predict(svmfit5, train)
table(Predict = predict5 , Truth= train$Y)## Truth
## Predict -1 1
## -1 343 5
## 1 5 347## [1] "Error de entrenamiento: 1.429 %"Cuando el parámetro \(cost\) es fijado en \(0.01\), \(690\) de las observaciones se han clasificado de forma correcta, solo \(10\) de las \(700\) observaciones del conjunto de entrenamiento quedaron mal clasificadas y por lo tanto el \(MSE\) de entrenamiento es del \(1.429\%\).
\(Cost=10\)
svmfit10 <- svm(Y~., data = train, kernel ="linear", cost = 10,
scale = T)
predict10 <- predict(svmfit10, train)
table(Predict = predict10 , Truth= train$Y)## Truth
## Predict -1 1
## -1 344 5
## 1 4 347## [1] "Error de entrenamiento: 1.286 %"\(Cost=15\)
svmfit15 <- svm(Y~., data = train, kernel ="linear", cost = 15,
scale = T)
predict15 <- predict(svmfit15, train)
table(Predict = predict15 , Truth= train$Y)## Truth
## Predict -1 1
## -1 344 5
## 1 4 347## [1] "Error de entrenamiento: 1.286 %"\(Cost=20\)
svmfit20 <- svm(Y~., data = train, kernel ="linear", cost = 20,
scale = T)
predict20 <- predict(svmfit20, train)
table(Predict = predict20 , Truth= train$Y)## Truth
## Predict -1 1
## -1 344 5
## 1 4 347## [1] "Error de entrenamiento: 1.286 %"A partir del valor de \(cost=10\), los valores de \(MSE\) se estabilizan, es decir que no cambian cuando se incrementa el valor del parámetro.
El resultado de la aplicación de la función \(tune()\) se guardo en un objeto llamado \(tune.out\), este objeto contiene un dataframe con información sobre los errores producidos con los diferentes valores de \(Cost\), dicho dataframe se presenta a continuación.
tune.out$performances## cost error dispersion
## 1 0.001 0.02142857 0.01214052
## 2 0.010 0.01714286 0.01475422
## 3 0.100 0.02142857 0.01543033
## 4 1.000 0.01714286 0.01475422
## 5 5.000 0.01428571 0.01649572
## 6 10.000 0.01428571 0.01649572
## 7 15.000 0.01571429 0.01572150
## 8 20.000 0.01571429 0.01572150La columna error contiene el error de validación para cada uno de los valores del parámetro \(cost\), observe que estos valores son un poco diferentes a los obtenidos calculando el \(MSE\) a partir de la matriz de confusión, esto puede deberse a que internamente la función \(tune()\) realiza una validación cruzada con \(k-folds\).
Literal c), Validación de resultados
Generar un conjunto de datos de prueba apropiado y calcular los errores de prueba correspondientes a cada uno de los valores de costo considerados. ¿Qué valor de costo conduce a la menor cantidad de errores de prueba y cómo se compara esto con los valores de costo que producen la menor cantidad de errores de entrenamiento y la menor cantidad de errores de validación cruzada?
No es necesario generar un nuevo conjunto de datos de prueba ya que los simulados en el literal (a) fueron divididos en entrenamiento y validación.
\(Cost=0.001\)
svmfit0.001 <- svm(Y~., data = train, kernel ="linear", cost = 0.001,
scale = T)
predict0.001 <- predict(svmfit0.001, test)
table(Predict = predict0.001 , Truth = test$Y)## Truth
## Predict -1 1
## -1 148 1
## 1 4 147## [1] "Error de entrenamiento: 1.667 %"La matriz anterior permite concluir que con valor del parámetro \(cost=0.001\), \(295\) de las observaciones se han clasificado de forma correcta, las restantes \(5\) de las \(300\) observaciones del conjunto de prueba quedaron mal clasificadas y por lo tanto el \(MSE\) de prueba es de \(1.667\%\).
\(Cost=0.01\)
svmfit0.01 <- svm(Y~., data = train, kernel ="linear", cost = 0.01,
scale = T)
predict0.01 <- predict(svmfit0.01, test)
table(Predict = predict0.01 , Truth= test$Y)## Truth
## Predict -1 1
## -1 147 1
## 1 5 147## [1] "Error de entrenamiento: 2 %"\(Cost=0.1\)
svmfit0.1 <- svm(Y~., data = train, kernel ="linear", cost = 0.1,
scale = T)
predict0.1 <- predict(svmfit0.1, test)
table(Predict = predict0.1 , Truth= test$Y)## Truth
## Predict -1 1
## -1 147 1
## 1 5 147## [1] "Error de entrenamiento: 2 %"\(Cost=1\)
svmfit1 <- svm(Y~., data = train, kernel ="linear", cost = 1,
scale = T)
predict1 <- predict(svmfit1, test)
table(Predict = predict1 , Truth= test$Y)## Truth
## Predict -1 1
## -1 147 1
## 1 5 147## [1] "Error de entrenamiento: 2 %"Para valores de \(cost\) igual a \((0.01, 0.1, 1)\) las tasas de mala clasificación en el conjunto de prueba son iguales.
\(Cost=5\)
svmfit5 <- svm(Y~., data = train, kernel ="linear", cost = 5,
scale = T)
predict5 <- predict(svmfit5, test)
table(Predict = predict5 , Truth= test$Y)## Truth
## Predict -1 1
## -1 145 2
## 1 7 146## [1] "Error de entrenamiento: 3 %"Cuando el parámetro \(cost\) es fijado en \(5\), \(291\) de las observaciones se han clasificado de forma correcta, solo \(9\) de las \(300\) observaciones del conjunto de entrenamiento quedaron mal clasificadas y por lo tanto el \(MSE\) de prueba es del \(3\%\).
\(Cost=10\)
svmfit10 <- svm(Y~., data = train, kernel ="linear", cost = 10,
scale = T)
predict10 <- predict(svmfit10, test)
table(Predict = predict10 , Truth= test$Y)## Truth
## Predict -1 1
## -1 148 2
## 1 4 146## [1] "Error de entrenamiento: 2 %"\(Cost=15\)
svmfit15 <- svm(Y~., data = train, kernel ="linear", cost = 15,
scale = T)
predict15 <- predict(svmfit15, test)
table(Predict = predict15 , Truth= test$Y)## Truth
## Predict -1 1
## -1 148 2
## 1 4 146## [1] "Error de entrenamiento: 2 %"\(Cost=20\)
svmfit20 <- svm(Y~., data = train, kernel ="linear", cost = 20,
scale = T)
predict20 <- predict(svmfit20, test)
table(Predict = predict20 , Truth= test$Y)## Truth
## Predict -1 1
## -1 148 2
## 1 4 146## [1] "Error de entrenamiento: 2 %"A partir del valor de \(cost=5\), los valores de \(MSE\) se estabilizan, es decir que no cambian cuando se incrementa el valor del parámetro.
Literal d), discusión de resultados
En el conjunto de entrenamiento se evidencia que el valor óptimo para el parámetro \(cost\) se encuentra entre los sigueintes valores \((10, 15, 20)\), mientras que para el conjunto de validación el valor de \(cost\) que produce la tasa de mala clasificación más baja es \(0.001\), cuando se usa este valor en el entrenamiento se obtiene que la tasa de mala clasificación es del \(2\%\), por lo tanto, es preferible un valor del parámetro que produzca mayores observaciones mal clasificadas en el conjunto de entrenamiento pero que tenga un desempeño mejor en el conjunto de validación.
Punto 7
En este problema, utilizará enfoques de vectores de soporte para predecir si un automóvil determinado obtiene un consumo alto o bajo de combustible según el conjunto de datos de Auto.
A continuación se presenta las primeras seis filas del conjunto de datos Auto, este conjunto contiene 9 variables.
library(ISLR)
head(Auto)## mpg cylinders displacement horsepower weight acceleration year origin
## 1 18 8 307 130 3504 12.0 70 1
## 2 15 8 350 165 3693 11.5 70 1
## 3 18 8 318 150 3436 11.0 70 1
## 4 16 8 304 150 3433 12.0 70 1
## 5 17 8 302 140 3449 10.5 70 1
## 6 15 8 429 198 4341 10.0 70 1
## name mpg01 forigin
## 1 chevrolet chevelle malibu FALSE 1
## 2 buick skylark 320 FALSE 1
## 3 plymouth satellite FALSE 1
## 4 amc rebel sst FALSE 1
## 5 ford torino FALSE 1
## 6 ford galaxie 500 FALSE 1-Mpg: Millas por galón.
-Cylinders: Número de cilindros entre 4 y 8.
-Displacement: Desplazamiento del motor (pulgadas cúbicos).
-Horsepower: Potencia del motor.
-Weight: Peso del vehículo (libras).
-Acceleration: Tiempo para acelerar de 0 60 mph (segundos).
-Year: Año del modelo (módulo 100).
-Origin: Origen del automóvil (1: Americano, 2: Europeo, 3: Japonés).
-Name: Nombre del vehículo.
Literal a), Creación de una variable categórica
Cree una variable binaria que tome el valor de \(1\) para los automóviles con rendimiento de gasolina por encima de la mediana y el valor de \(0\) para automóviles con rendimiento de combustible por debajo de la mediana.
medianamillaje <- median(Auto$mpg)
Auto$High <- ifelse(Auto$mpg < medianamillaje, "0", "1")
Auto$High <- as.factor(Auto$High)La mediana de millas por galón es de \(22.75\) y la proporción de las categorías de la nueva variable llamada High es de \(50\%\).
Literal b), Ajuste de un clasificador de soporte vectorial
Ajuste un clasificador de vectores de soporte a los datos con varios valores de \(cost\), para predecir si un automóvil posee un alto o bajo consumo de combustible. Informe los errores de validación cruzada asociados con diferentes valores de este parámetro. Comente sus resultados.
Inicialmente el conjunto de datos será dividido, el \(70\%\) de las observaciones se asignan a la base de entrenamiento y el restante a la base de prueba.
Luego de esto se usa la función tune() perteneciente a la librería e1071 que realiza una validación cruzada de diez folds en un conjunto de modelos de interés. El parámetro \(Cost\) tomara valores en el vector \((0.001, 0.01, 0.1, 1, 5, 10, 100)\), como se muestra a continuación,
library(e1071)
set.seed(123456)
tune.out.l <- tune(svm, High~., data=train, kernel ="linear",
ranges =list(cost=c(0.001, 0.01, 0.1, 1, 5, 10, 100)),
scale = T)En la salida \(tune.out.l\) se guarda un conjunto de datos que contiene los valores del parámetro \(cost\), sus errores y dispersión, en base a esta información se construye el siguiente gráfico.
Se evidencia que el mínimo error se produce cuando el parámetro \(cost=0.01\), dicho error es igual a \(0\), a continuación se presenta el mismo gráfico pero con los errores para el conjunto de prueba.
Cuando el parámetro \(cost\) es \(0.01\) se obtiene un valor de \(Error=0.025\), es decir que hay muy pocas observaciones mal clasificadas.
Litera c), Usar kernell radial
Ahora repita (b), esta vez usando SVMs con núcleos de base radial y polinomial, con diferentes valores de \(\gamma\) y \(degree\) y \(cost\). Comente sus resultados.
Kernell Radial
Nuevamente se hace uso de la función tune(), variando los parámetros de \(Cost\) y \(\gamma\).
set.seed(1233456)
tune.out.r <- tune(svm, High~., data = train, kernel = "radial",
ranges =list(cost=c(0.001, 0.01, 0.1, 1, 5, 10, 100),
gamma=c(0.1, 0.5, 1, 2, 3, 4, 5)),
scale = T)
Se identifica que la combinación óptima de los parámetros \(Cost\) y \(\gamma\) son \(1\) y \(0.1\) resectivamente.
Se procede a ajustar la SMV de Kernell radial, con la combinación óptima de los parámetros \(Cost\) y \(\gamma\).
##
## Call:
## best.tune(method = svm, train.x = High ~ ., data = train, ranges = list(cost = c(0.001,
## 0.01, 0.1, 1, 5, 10, 100), gamma = c(0.1, 0.5, 1, 2, 3, 4, 5)),
## kernel = "radial", scale = T)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 1
##
## Number of Support Vectors: 67
##
## ( 34 33 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1Después de realizar el mejor modelo se obtiene que hay \(67\) vectores de soporte, \(34\) de una clase y \(33\) de la otra. Ahora se predicen las etiquetas de clase de las observaciones de prueba creadas anteriomente, para el mejor modelo usando validación cruzada.
## Truht
## Predict 0 1
## 0 59 0
## 1 0 59Finalmente se puede concluir que con valor asignado a \(cost=1\) y \(\gamma=0.1\), todas las observaciones del conjunto de prueba se han clasificado de manera y por lo tanto el \(MSE\) de prueba es del \(0\%\).
Kernell Polinomial
Nuevamente se hace uso de la función tune(), variando los parámetros de \(Cost\) y \(degree\), ya que el kernell es polinomial.
set.seed(123456)
tune.out.p <- tune(svm, High~., data = train, kernel = "polynomial",
ranges =list(cost = c(0.001,0.01, 0.1, 1, 10,100),
degree = c(2, 3, 4, 5, 6)),
scale = TRUE)
Este kernell produce errores bastante grandes, aún así se identifica que la combinación óptima de los parámetros \(Cost\) y \(degree\) son \(100\) y \(2\) resectivamente.
##
## Call:
## best.tune(method = svm, train.x = High ~ ., data = train, ranges = list(cost = c(0.001,
## 0.01, 0.1, 1, 10, 100), degree = c(2, 3, 4, 5, 6)), kernel = "polynomial",
## scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 100
## degree: 2
## coef.0: 0
##
## Number of Support Vectors: 224
##
## ( 111 113 )
##
##
## Number of Classes: 2
##
## Levels:
## 0 1Después de realizar el ajuste del mejor modelo se obtiene que hay \(224\) vectores de soporte, \(111\) de una clase y \(113\) de la otra. Ahora se predicen las etiquetas de clase de las observaciones de prueba creadas anteriomente, para el mejor modelo usando validación cruzada.
## Truht
## Predict 0 1
## 0 50 1
## 1 9 58Finalmente se puede concluir que con valor asignado a \(cost=100\) y \(degree = 2\), \(108\) de las observaciones se han clasificado de forma correcta y \(10\) de ellas quedaron mal clasificadas, por lo tanto el \(MSE\) de prueba es de \(8.47\%\).
Literal d), Gráficos
Haga algunas gráficas para respaldar sus afirmaciones en (b) y (c).
Clasificador de soporte vectorial
Es necesario verificar si las dos clases son completamente separables.
En los gráficos presentados anteriormente se puede concluir que las observaciones no son completamente separables, se ajustó el clasificador de soporte vectorial con un kernel lineal, donde la región del espacio de atributos que se asignó fue a la clase \(0\) que se muestra en amarillo claro. Recordar que se obtuvo que el valor óptimo de \(cost=5\), con \(47\) vectores de soporte. Finalmente por medio de las curvas ROC a pesar de que los datos no son completamente separables, el clasificador de soporte se ajusta bien a los datos.
SVM Kernell Radial
Nuevamente necesario verificar si las dos clases son completamente separables bajo este modelo.
Los gráficos presentados anteriormente permiten concluir que las observaciones no son completamente separables, se ajusto una máquina de soporte vectorial en la cual se obtuvo que el valor óptimo de \(cost=1\) siendo grande, \(gamma=0.1\) y con \(67\) vectores de soporte, las curvas ROC muestran que la máquina de soporte vectorial con kernel radial es un buen clasificador.
SVM Kernell Polinomial
Los gráficos previamente mostrados permiten concluir que las observaciones no son completamente separables, se ajustó el clasificador de soporte vectorial con un kernel polinomial, donde la región del espacio de atributos que se asignó fue a la clase \(1\) que se muestra en purpura y en algunos casos se evidencio una región de color amarillo claro que pertenece a la clase \(0\). Recordar que se obtuvo que el valor óptimo de \(cost=100\) siendo grande, con \(224\). Finalmente a partir de las curvas ROC se observa que se produce mayor error de test, lo que indica que no es tan buen clasificador.
Punto 8
Este problema involucra el conjunto de datos OJ que forma parte del paquete ISLR.
Se presenta el encabezado del conjunto de datos a trabajar, este contiene \(18\) variables.
## Purchase WeekofPurchase StoreID PriceCH PriceMM DiscCH DiscMM SpecialCH
## 1 CH 237 1 1.75 1.99 0.00 0.0 0
## 2 CH 239 1 1.75 1.99 0.00 0.3 0
## 3 CH 245 1 1.86 2.09 0.17 0.0 0
## 4 MM 227 1 1.69 1.69 0.00 0.0 0
## 5 CH 228 7 1.69 1.69 0.00 0.0 0
## 6 CH 230 7 1.69 1.99 0.00 0.0 0
## SpecialMM LoyalCH SalePriceMM SalePriceCH PriceDiff Store7 PctDiscMM
## 1 0 0.500000 1.99 1.75 0.24 No 0.000000
## 2 1 0.600000 1.69 1.75 -0.06 No 0.150754
## 3 0 0.680000 2.09 1.69 0.40 No 0.000000
## 4 0 0.400000 1.69 1.69 0.00 No 0.000000
## 5 0 0.956535 1.69 1.69 0.00 Yes 0.000000
## 6 1 0.965228 1.99 1.69 0.30 Yes 0.000000
## PctDiscCH ListPriceDiff STORE
## 1 0.000000 0.24 1
## 2 0.000000 0.24 1
## 3 0.091398 0.23 1
## 4 0.000000 0.00 1
## 5 0.000000 0.00 0
## 6 0.000000 0.30 0También es de importancia identificar la distribución de la variable respuesta, la mayoría pertenece a la categoría CH, en una prorporción de \(61\%\) y el restante (\(39\%\)) para MM.
Literal a), División del conjunto
Cree un conjunto de entrenamiento que contenga una muestra aleatoria de 800 observaciones y un conjunto de prueba que contenga las observaciones restantes.
set.seed(123456)
index <- sample(1:nrow(OJ), 800)
train <- OJ[index,]
test <- OJ[-index,]Literal b), Clasificador de soporte vectorial
Ajuste un clasificador de vector de soporte a los datos de entrenamiento usando \(cost = 0.01\), con Purchase como respuesta y las otras variables como predictores. Utilice la función \(summary()\) para producir estadísticas de resumen y describir los resultados obtenidos.
svmfit <- svm(Purchase~., data = train, kernel ="linear", cost = 0.01,
scale = TRUE)
summary(svmfit)##
## Call:
## svm(formula = Purchase ~ ., data = train, kernel = "linear", cost = 0.01,
## scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: linear
## cost: 0.01
##
## Number of Support Vectors: 443
##
## ( 221 222 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMLa variable respuesta Purchase es un factor que posee dos clases, una de ellas es CH y la otra MM, se uso un kernell lineal para el ajuste de clasificador de soporte vectorial, este clasificador tiene \(443\) vectores de soporte, tomo \(221\) de una clase y los restantes \(222\) de la otra clase, al asignar un valor tan pequeño al parámetro \(cost\) la margen es grande y habrán mucho vectores de sporte.
Literal c), Tasas de error
¿Cuáles son las tasas de error de entrenamiento y prueba?
predict_train <- predict(svmfit, train)
table(Predict = predict_train, Truth = train$Purchase)## Truth
## Predict CH MM
## CH 423 75
## MM 59 243Se obtiene que \(134\) de las \(800\) observaciones de entrenamiento estan mal clasificadas, es decir, que se obtiene una tasa de mala clasificación del \(16.75\%\).
predict_test <- predict(svmfit, test)
table(Predict = predict_test, Truth = test$Purchase)## Truth
## Predict CH MM
## CH 152 23
## MM 19 76En el conjunto de prueba \(42\) de las \(270\) observaciones estan mal clasificadas, el casificador de soporte vectorial tiene una tasa de mala clasificación del \(15.55\%\) en el conjunto de prueba.
Literal d), Validación cruzada para elección de parámetros
Utilice la función tune () para seleccionar un costo óptimo. Considere valores en el rango de 0.01 a 10.
set.seed(123456)
tune.out <- tune(svm ,Purchase~., data = train, kernel ="linear",
ranges =list(cost=c(seq(0.01, 10, 0.1))), scale = T)
El gráfico anterior permite concluir que el valor óptimo para el parámetro \(Cost\) es \(0.01\).
Literal e), Tasas de error
Calcule las tasas de error de entrenamiento y prueba utilizando este nuevo valor de costo.
svmoptimo <- svm(Purchase~., data = train, kernel ="linear", cost = 0.01,
scale = TRUE)predict_trainoptimo <- predict(svmoptimo, train)
table(Predict = predict_trainoptimo, Truth = train$Purchase)## Truth
## Predict CH MM
## CH 423 75
## MM 59 243Se obtiene que \(134\) de las \(800\) observaciones de entrenamiento estan mal clasificadas, es decir, que se ontiene una tasa de mala clasificación del \(16.75\%\).
predict_testoptimo <- predict(svmoptimo, test)
table(Predict = predict_testoptimo, Truth = test$Purchase)## Truth
## Predict CH MM
## CH 152 23
## MM 19 76En el conjunto de prueba \(42\) de las \(270\) observaciones estan mal clasificadas, el casificador de soporte vectorial tiene una tasa de mala clasificación del \(15.55\%\) que igual a las tasas de error calculadas en el literal (c).
Literal f), Máquina de soporte vectorial con kernel radial
Repita las partes (b) a (e) usando una máquina de vectores de soporte con un núcleo radial. Utilice el valor predeterminado para gamma.
Ajuste de la máquina de soporte vectorial kernel radial y \(cost=0.01\)
svmfit <- svm(Purchase~., data = train, kernel ="radial", cost = 0.01,
scale = TRUE)
summary(svmfit)##
## Call:
## svm(formula = Purchase ~ ., data = train, kernel = "radial", cost = 0.01,
## scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: radial
## cost: 0.01
##
## Number of Support Vectors: 639
##
## ( 318 321 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMSe uso un kernel radial para el ajuste de la máquina de soporte vectorial, este clasificador tiene \(639\) vectores de soporte, tomo \(318\) de una clase y los restantes \(321\) de la otra clase.
Ahora se calculan las tasas de error producidas por esta máquina.
predict_train <- predict(svmfit, train)
table(Predict = predict_train, Truth = train$Purchase)## Truth
## Predict CH MM
## CH 482 318
## MM 0 0Se obtiene que \(318\) de las \(800\) observaciones de entrenamiento estan mal clasificadas, es decir, que se obtiene una tasa de mala clasificación del \(39.75\%\).
predict_test <- predict(svmfit, test)
table(Predict = predict_test, Truth = test$Purchase)## Truth
## Predict CH MM
## CH 171 99
## MM 0 0En el conjunto de prueba \(99\) de las \(270\) observaciones estan mal clasificadas, la máquina de soporte vectorial con kernel radial tiene una tasa de mala clasificación del \(36.66\%\) en el conjunto de prueba.
Uso de la función tune() para encontrar el valor óptimo de \(cost\)
set.seed(123456)
tune.out <- tune(svm, Purchase~., data = train, kernel ="radial",
ranges =list(cost=c(seq(0.01, 10, 0.1))), scale = TRUE)
El gráfico anterior perite concluir que el valor óptimo para el parámetro \(Cost\) es \(0.51\).
svmfit <- svm(Purchase~., data = train, kernel ="radial", cost = 0.51,
scale = TRUE)Ahora se calculan las tasas de error producidas por la máquina con el parámetro de \(cost\) óptimo.
## Truth
## Predict CH MM
## CH 438 79
## MM 44 239Se obtiene que \(123\) de las \(800\) observaciones de entrenamiento estan mal clasificadas, es decir, que se obtiene una tasa de mala clasificación del \(15.375\%\).
predict_test <- predict(svmfit, test)
table(Predict = predict_test, Truth = test$Purchase)## Truth
## Predict CH MM
## CH 156 22
## MM 15 77En el conjunto de prueba \(37\) de las \(270\) observaciones estan mal clasificadas, la máquina de sorporte vectorial con kernel radial tiene una tasa de mala clasificación del \(13.70%\) en el conjunto de prueba. Ambos errores son menores a los obtenidos cuando se uso el parámetro \(cost=0.01\).
Literal g), Máquina de soporte vectorial con kernel polinomial
Repita las partes (b) a (e) usando una máquina de vectores de soporte con un núcleo polinomial. Establecer \(degree = 2\).
Ajuste de la máquina de soporte vectorial con kernel polinomial, \(cost=0.01\) y \(degree = 2\)
svmfit <- svm(Purchase~., data = train, kernel ="polynomial", cost = 0.01,
scale = TRUE, degree = 2)
summary(svmfit)##
## Call:
## svm(formula = Purchase ~ ., data = train, kernel = "polynomial",
## cost = 0.01, degree = 2, scale = TRUE)
##
##
## Parameters:
## SVM-Type: C-classification
## SVM-Kernel: polynomial
## cost: 0.01
## degree: 2
## coef.0: 0
##
## Number of Support Vectors: 639
##
## ( 318 321 )
##
##
## Number of Classes: 2
##
## Levels:
## CH MMSe uso un kernel polinomial de \(2\) grados para el ajuste de la máquina de soporte vectorial, este clasificador tiene \(639\) vectores de soporte, tomo \(318\) de una clase y los restantes \(321\) de la otra clase.
Ahora se calculan las tasas de error producidas por esta máquina.
## Truth
## Predict CH MM
## CH 482 312
## MM 0 6Se obtiene que \(312\) de las \(800\) observaciones de entrenamiento estan mal clasificadas, es decir, que se obtiene una tasa de mala clasificación del \(39\%\).
## Truth
## Predict CH MM
## CH 171 98
## MM 0 1En el conjunto de prueba \(98\) de las \(270\) observaciones estan mal clasificadas, la máquina de soporte vectorial con kernel polinomial tiene una tasa de mala clasificación del \(36.2\%\) en el conjunto de prueba.
Uso de la función tune() para encontrar el valor óptimo de \(cost\)
set.seed(123456)
tune.out <- tune(svm, Purchase~., data = train, kernel ="polynomial",
degree = 2, ranges =list(cost=c(seq(0.01, 10, 0.1))),
scale = TRUE)
El gráfico anterior permite concluir que el valor óptimo para el parámetro \(Cost\) es \(3.41\).
svmfit <- svm(Purchase~., data = train, kernel ="polynomial", cost = 3.41,
scale = TRUE, degree = 2)Ahora se calculan las tasas de error producidas por la máquina con el parámetro de \(cost\) óptimo.
## Truth
## Predict CH MM
## CH 448 93
## MM 34 225Se obtiene que \(127\) de las \(800\) observaciones de entrenamiento estan mal clasificadas, es decir, que se obtiene una tasa de mala clasificación del \(15.875\%\).
## Truth
## Predict CH MM
## CH 157 26
## MM 14 73En el conjunto de prueba \(40\) de las \(270\) observaciones estan mal clasificadas, la máquina de soporte vectorial con kernell polinomial tiene una tasa de mala clasificación del \(14.8148\%\) en el conjunto de prueba.
Ambos errores son menores a los obtenidos cuando se uso el parámetro \(cost=0.01\).
Literal h), Resultados
En general, ¿qué enfoque parece dar los mejores resultados con estos datos?
El método que parece proporcionar los mejores resultados es la maquina de soporte vectorial con un kernell radial con \(\gamma = 0.1\) y parámetro \(cost=0.51\).
Finalmente se adjunta el enlace para acceder al ensayo, este fue elaborado sobre el tema de inseguridad en Medellín.