Tarea3
Gonzalo Estebaranz
25 de noviembre de 2017
Introducción
Para esta tarea se ha utilizado un dataset correspondiente a información escolar de una muestra de alumnos. El objetivo principal va a ser predecir la nota G3 de los alumnos. Dicho dataset está disponible en: https://archive.ics.uci.edu/ml/machine-learning-databases/00320/student.zip
#Fijamos directorio de trabajo
setwd("C:/Users/Gonzalo/Documents/datos/tarea3")
rm(list=ls())
ls()
fileURL <- "http://archive.ics.uci.edu/ml/machine-learning-databases/00320/student.zip"
download.file(fileURL,"student.zip", mode="wb")
unzip("student.zip") #Almacenamos la fecha de la descarga
fechaDescarga <- date()
fechaDescarga## [1] "Sat Nov 25 19:06:52 2017"#Librerías necesarias#
library(knitr)###CARGA DE LOS DATOS Y ANÁLISIS DESCRIPTIVO###
#Lectura de los datos#
setwd("C:/Users/Gonzalo/Documents/datos/tarea3")
mat <- file("student-mat.csv","r")
studentMat <- read.csv2(mat)
#Realizamos un primer anális descriptivo de las variable con summary#
summary(studentMat)## school sex age address famsize Pstatus Medu
## GP:349 F:208 Min. :15.0 R: 88 GT3:281 A: 41 Min. :0.000
## MS: 46 M:187 1st Qu.:16.0 U:307 LE3:114 T:354 1st Qu.:2.000
## Median :17.0 Median :3.000
## Mean :16.7 Mean :2.749
## 3rd Qu.:18.0 3rd Qu.:4.000
## Max. :22.0 Max. :4.000
## Fedu Mjob Fjob reason
## Min. :0.000 at_home : 59 at_home : 20 course :145
## 1st Qu.:2.000 health : 34 health : 18 home :109
## Median :2.000 other :141 other :217 other : 36
## Mean :2.522 services:103 services:111 reputation:105
## 3rd Qu.:3.000 teacher : 58 teacher : 29
## Max. :4.000
## guardian traveltime studytime failures schoolsup
## father: 90 Min. :1.000 Min. :1.000 Min. :0.0000 no :344
## mother:273 1st Qu.:1.000 1st Qu.:1.000 1st Qu.:0.0000 yes: 51
## other : 32 Median :1.000 Median :2.000 Median :0.0000
## Mean :1.448 Mean :2.035 Mean :0.3342
## 3rd Qu.:2.000 3rd Qu.:2.000 3rd Qu.:0.0000
## Max. :4.000 Max. :4.000 Max. :3.0000
## famsup paid activities nursery higher internet romantic
## no :153 no :214 no :194 no : 81 no : 20 no : 66 no :263
## yes:242 yes:181 yes:201 yes:314 yes:375 yes:329 yes:132
##
##
##
##
## famrel freetime goout Dalc
## Min. :1.000 Min. :1.000 Min. :1.000 Min. :1.000
## 1st Qu.:4.000 1st Qu.:3.000 1st Qu.:2.000 1st Qu.:1.000
## Median :4.000 Median :3.000 Median :3.000 Median :1.000
## Mean :3.944 Mean :3.235 Mean :3.109 Mean :1.481
## 3rd Qu.:5.000 3rd Qu.:4.000 3rd Qu.:4.000 3rd Qu.:2.000
## Max. :5.000 Max. :5.000 Max. :5.000 Max. :5.000
## Walc health absences G1
## Min. :1.000 Min. :1.000 Min. : 0.000 Min. : 3.00
## 1st Qu.:1.000 1st Qu.:3.000 1st Qu.: 0.000 1st Qu.: 8.00
## Median :2.000 Median :4.000 Median : 4.000 Median :11.00
## Mean :2.291 Mean :3.554 Mean : 5.709 Mean :10.91
## 3rd Qu.:3.000 3rd Qu.:5.000 3rd Qu.: 8.000 3rd Qu.:13.00
## Max. :5.000 Max. :5.000 Max. :75.000 Max. :19.00
## G2 G3
## Min. : 0.00 Min. : 0.00
## 1st Qu.: 9.00 1st Qu.: 8.00
## Median :11.00 Median :11.00
## Mean :10.71 Mean :10.42
## 3rd Qu.:13.00 3rd Qu.:14.00
## Max. :19.00 Max. :20.00## Cambiamos nombres de las variables
#Modificamos los nombres de las variables para que vayan en minúsculas#
names(studentMat) <- tolower(names(studentMat))
attach(studentMat)
#Eliminamos school, mjob, fjob, reason y guardian que son nominales para ver la matriz de correlaciones#
student<-studentMat[,-c(1,9,10,11,12)]
#Convertimos variables binarias no/yes a 0/1 y todas a numericas
student$schoolsup <- as.numeric(studentMat$schoolsup) - 1
student$famsup <- as.numeric(studentMat$famsup) - 1
student$paid <- as.numeric(studentMat$paid) - 1
student$activities <- as.numeric(studentMat$activities) - 1
student$nursery <- as.numeric(studentMat$nursery) - 1
student$higher <- as.numeric(studentMat$higher) - 1
student$internet <- as.numeric(studentMat$internet) - 1
student$romantic <- as.numeric(studentMat$romantic) - 1
#Matriz de correlación solo aplicable a variables numericas#
student<-sapply(student,as.numeric)
matCor <- cor(student)
#Nos quedamos con las correlaciones respecto a G3 (variable relativa a la nota final)
matCor <- as.matrix(matCor[,28])
colnames(matCor)<-"G3"
matCor## G3
## sex 0.10345565
## age -0.16157944
## address 0.10575606
## famsize 0.08140711
## pstatus -0.05800898
## medu 0.21714750
## fedu 0.15245694
## traveltime -0.11714205
## studytime 0.09781969
## failures -0.36041494
## schoolsup -0.08278821
## famsup -0.03915715
## paid 0.10199624
## activities 0.01609970
## nursery 0.05156790
## higher 0.18246462
## internet 0.09848337
## romantic -0.12996995
## famrel 0.05136343
## freetime 0.01130724
## goout -0.13279147
## dalc -0.05466004
## walc -0.05193932
## health -0.06133460
## absences 0.03424732
## g1 0.80146793
## g2 0.90486799
## g3 1.00000000#Podemos ver que las relaciones son bastante leves#
#G1 y G2 no nos interesan como correladas con G3 (la nota final es resultado de las dos notas anteriores)#
#Aun así, los alumnos con más failures tienen una nota más baja.
#Aquellos que quieren ir a la universidad tienen una nota más alta (higher)
#Buscamos dicha relación graficamente#
student2<-data.frame(student)
if(! "gridExtra" %in% installed.packages()) install.packages("gridExtra", depend = TRUE)
library(gridExtra)
library(ggplot2)
miGgplot1 <- qplot(failures, g3, data = student2)
miGgplot1
miGgplot2 <- boxplot(g3~higher, data = student2)
Análisis exploratorio apoyado en algún método NO supervisado
Vamos a realizar un clustering o agrupamiento de nuestros alumnos en aprobados/suspensos en función del resto de variables. Dado que nuestras variables son de diferentes tipos utilizamos la distancia de Gower para nuestro cluster jerárquico
##Librerias usadas##
if(! "cluster" %in% installed.packages()) install.packages("cluster", depend = TRUE)
library(cluster)
if(! "dplyr" %in% installed.packages()) install.packages("dplyr", depend = TRUE)
library(dplyr)
library(Hmisc)#Nos quedamos las variables más influyentes en G3#
student3<-student2[,c(6,7,10,16,18,21,28)]
#Creamos una variable binaria que diga 1 si el alumno prueba y 0 si no#
notag3 <-cut2(student3$g3, c(10,20))
levels(notag3) <- c(0,1)
student4<-cbind(student3,notag3)
#Utilizamos la distancia de Gower,útil con diferentes tipos de variables (binarias, categóricas,...)#
student.mod <- student4 %>% select(-c(g3,notag3))
idsamp <- sample(1:dim(student.mod)[1], 360)
student.samp <- student.mod[idsamp, ]
head(student.samp)## medu fedu failures higher romantic goout
## 146 1 1 0 1 0 2
## 26 2 2 2 1 0 2
## 370 4 4 0 1 1 2
## 19 3 2 3 1 0 5
## 264 3 3 0 1 0 3
## 41 2 2 1 1 1 3D<-daisy(student.samp,"gower")Aplicamos clustering jerarquico
H.fit <- hclust(D, method="ward.D")
plot(H.fit, hang = -1, cex=0.8, labels=student4$notag3[idsamp])
groups <- cutree(H.fit, k=2)
rect.hclust(H.fit, k=2, border="red")
##Solo se explica algo más del 50%. El modelo no parece muy adecuado##
clusplot(student.samp, groups, color=TRUE, shade=TRUE,
labels=2, lines=0, main= 'Customer segments')
Analisis exploratorio apoyado en algun metodo supervisado
En primer lugar vamos a realizar una regresión logística que intente predecir si un alumno aprueba o suspende. En este modelo no se incluyen variables como G1,G2 o G3 ya que no sería util si ya incluimos sus notas previas. Tampoco se incluyen otras como Pstatus, trabajo de los padres, reason o guardian.
#Eliminamos algunas variables que no queremos incluir en el estudio por ser poco útiles#
studentMat<-cbind(studentMat,notag3)
studentMat<-studentMat[,-c(6,9,10,11,12,31,32,33)]
#Realizamos la muestra de test y de training#
set.seed(12345)
train.sample <- sample(1:395,size=(350),replace=F)
train.student <- studentMat[train.sample,]
test.student <- studentMat[-train.sample,]
dim(train.student)## [1] 350 26#Realizamos el modelo de regresión logística#
logistic.reg.model <- glm(notag3 ~.,family=binomial,data=train.student)Sin embargo, este modelo de regresión tiene muchas variables. En estos casos, lo adecuado es buscar un equilibrio de tal forma que podamos explicar lo máximo posible con el menor número de variables. Por ello, realizamos stepwize regression.
logistic.reg.model<-step(logistic.reg.model)Parece que failures y goout son las variables más significativas a la hora de aprobar o suspender.
summary(logistic.reg.model)##
## Call:
## glm(formula = notag3 ~ age + failures + schoolsup + famsup +
## goout + walc + absences, family = binomial, data = train.student)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.2130 -1.0221 0.6129 0.8160 2.1581
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.91198 1.89330 3.123 0.00179 **
## age -0.20883 0.10961 -1.905 0.05675 .
## failures -0.90823 0.18470 -4.917 8.78e-07 ***
## schoolsupyes -0.71355 0.37199 -1.918 0.05508 .
## famsupyes -0.49386 0.26984 -1.830 0.06722 .
## goout -0.39783 0.12773 -3.115 0.00184 **
## walc 0.19438 0.11184 1.738 0.08221 .
## absences -0.02393 0.01450 -1.650 0.09892 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 443.21 on 349 degrees of freedom
## Residual deviance: 383.11 on 342 degrees of freedom
## AIC: 399.11
##
## Number of Fisher Scoring iterations: 4Tratamos ahora de realizar un arbol de decisión para nuestros datos que, de igual manera, nos permita predecir si un alumno suspende o aprueba la asignatura.
#Librerias usadas#
if(! "rpart" %in% installed.packages()) install.packages("rpart", depend = TRUE)
library(rpart)
library(rpart.plot)
library(RColorBrewer)
library(rattle)## Rattle: A free graphical interface for data science with R.
## Versión 5.1.0 Copyright (c) 2006-2017 Togaware Pty Ltd.
## Escriba 'rattle()' para agitar, sacudir y rotar sus datos.#Arbol de decision#
mytree <- rpart(notag3 ~ ., data = train.student, method = "class")
#Representamos el arbol#
fancyRpartPlot(mytree)
#Podemos ver que hay mayores posibilidades de aprobar al dedicar más tiempo de estudio, faltar menos a clase o cuando la educación de la madre es superior#EVALUACIÓN DEL MODELO 1 (REGRESIÓN)
Utilizamos Curvas ROC, el área bajo la curva y la precisión para determinar la calidad de nuestro modelo y cuál es mejor.
#Librerias usadas#
if(! "ROCR" %in% installed.packages()) install.packages("ROCR", depend = TRUE)
library(ROCR)#Aplicamos el modelo de regresión al conjunto de test#
test.nota.output1 <- predict(logistic.reg.model,newdata=test.student,type='response')
head(test.nota.output1)## 10 23 36 41 52 56
## 0.8891681 0.9375093 0.8891681 0.4417075 0.7753486 0.7275851#Obtenemos un vector con el porcentaje de que la nota sea OK
#Si mayor que 0.5 -> OK. Si no, suspenso.
test.nota.output1 <- ifelse(test.nota.output1 > 0.5,1,0)
head(test.nota.output1)## 10 23 36 41 52 56
## 1 1 1 0 1 1clasif.error1 <- mean(test.nota.output1 != test.student$notag3)
clasif.error1## [1] 0.3111111##Vemos la buena precision (73%) que tiene nuestro modelo##
precision1 <- 1 - clasif.error1
precision1## [1] 0.6888889##Calculamos la curva ROC##
pr1 <- prediction(test.nota.output1, test.student$notag3)
prf1 <- performance(pr1, measure = "tpr", x.measure = "fpr")
plot(prf1)
auc1 <- performance(pr1, measure = "auc")
auc1 <- auc1@y.values[[1]]
auc1## [1] 0.6##Obtenemos un valor de 0.65. Bastante aceptable para nuestro modelo##EVALUACIÓN DEL MODELO 2 (ARBOL DE DECISIÓN)
#Calculamos su error#
test.nota.output2 <- predict(mytree, newdata = test.student,type="vector")
test.nota.output2 <- ifelse(test.nota.output2 > 0.5,1,0)
head(test.nota.output2)## 10 23 36 41 52 56
## 1 1 1 1 1 1clasif.error2 <- mean(test.nota.output2 != test.student$notag3)
clasif.error2## [1] 0.3333333##Vemos la precision (66%) que tiene nuestro modelo##
precision2 <- 1 - clasif.error2
precision2## [1] 0.6666667pr2 <- prediction(test.nota.output2, test.student$notag3)
prf2 <- performance(pr2, measure = "tpr", x.measure = "fpr")
plot(prf2)
auc2 <- performance(pr2, measure = "auc")
auc2<- auc2@y.values[[1]]
auc2## [1] 0.5##Obtenemos un area bajo la curva de 0'5##Parece con todo ello , que el modelo óptimo sería el de la regresión ya que tiene mayor precisión (0.68 frente a 0.66), y un área bajo la curva superior (0.6 frente a 0.5)