---
title: "Métodos de Agregación Homogéneos"
author: "Mathias Bourel"
date: "17/11/2018"
output:
  pdf_document: default
  html_document: default
---
```{r}
rm(list=ls()) 
```


#1 - Paquetes que vamos/podemos usar
```{r}
library(rpart) #arboles
library(ipred) #bagging
library(randomForest)#random forests
#library(caret)#paquete con varias cosas de machine learning
```

#2 - Regresión
```{r}
data=read.csv("housing_data.csv",sep=";", dec=",", header=T)
summary(data)
names(data)
data$CHAS=as.factor(data$CHAS)
summary(data$CHAS)

```


Construir un árbol de regresión considerando las variables CRIM, ZN, CHAS, NOX, RM, DIS, TAX y PTRATIO.
Hallar el arbol optimo dado por el algoritmo de la poda con el minimo error de validación cruzada.

Hallar el arbol optimo dado por el algoritmo de la poda con la regla 1-SE.

Hallar el arbol optimo con paramétro de costo complejidad igual a cp=0.01 (es el que devuelve la función rpart)
```{r}
tree=rpart(MEDV~CRIM+ZN+CHAS+NOX+RM+DIS+TAX+PTRATIO,cp=0.01,data)
cp.opt = tree$cptable[which.min(tree$cptable[,"xerror"]),"CP"]
tree2=rpart(MEDV~CRIM+ZN+CHAS+NOX+RM+DIS+TAX+PTRATIO,cp=cp.opt,data)
tsse=min(tree$cptable[,"xerror"])+tree$cptable[which.min(tree$cptable[,"xerror"]),"xstd"]
cpopt1SE=tree$cptable[min(tree$cptable[tree$cptable[,"xerror"]<tsse,"nsplit"])+1,"CP"]
tree3=rpart(MEDV~CRIM+ZN+CHAS+NOX+RM+DIS+TAX+PTRATIO,cp=cpopt1SE,data)
plot(tree,margin=0.1)
text(tree)

```



Calcule el error del árbol sobre la muestra de entrenamietno

```{r}
pred=predict(tree) 
error.tree.lrn= sqrt(mean((data[,1]-pred)^2)) 
error.tree.lrn
MEDV=data$MEDV
error.tree.lrn/mean(MEDV)
```

Realize el error sobre la muestra de prueba (media y SD de 20 iteraciones y partiendo la muestra en 2/3 - 1/3)

```{r}
K=20
	error.tree=list(matrix(NA, K))
	n = nrow(data)
	for(k in 1:K)	{
		smp=sample(n,round(n/3))
		learn=data[-smp,]
		test=data[smp,]
			tree.lrn=rpart(MEDV~CRIM+ZN+CHAS+NOX+RM+DIS+TAX+PTRATIO,cp=0.01,learn)
			pred=predict(tree.lrn,test)
			error.tree[k] = sqrt(mean((test[,1]-pred)^2)) 
	
	}				
	
mean.error.tree=mean(unlist(error.tree))
sd.error.tree=sd(unlist(error.tree))
mean.error.tree
sd.error.tree

	
```


	


#2-1 BAGGING
Construya un modelo Bagging con la función bagging del paquete ipred. 
Pida que el estimador se construya con 25 remuestras bootstrap y 
que la salida le brinde el error de las observaciones ``out of bag''. 
Indique que para la construcción de cada árbol se tome en cuenta un cp=0.01. 
Realice el gráfico del error cuadrático medio en función de la cantidad de arboles. 
Obtenga el error por muestra de prueba y su desviación estándar de 
la misma forma que se hizo en CART.


```{r}
bag=ipredbagg(data[,1], data[,c(2,3,5,6,7,9,11,12)], control=rpart.control(cp = 0.01),nbagg=25, coob=F)
bag
summary(bag)
print(bag)

bag2=bagging(MEDV ~ CRIM + ZN + CHAS + NOX + RM + DIS  + TAX + PTRATIO,control=rpart.control(cp=0.01),nbagg=25, coob=T,data=data)

bag3=ipredbagg(data[,1], data[,c(2,3,5,6,7,9,11,12)], control=rpart.control(cp = 0.01),nbagg=100, coob=TRUE)
bag3

K=20
	error.bag=list(matrix(NA, K))
	n = nrow(data)
	for(k in 1:K)	{
		smp=sample(n,round(n/3))
		learn=data[-smp,]
		test=data[smp,]
			bag.lrn=ipredbagg(learn[,1], learn[,c(2,3,5,6,7,9,11,12)], control=rpart.control(cp = 0.01),nbagg=25)
			pred=predict(tree.lrn,test)
			error.bag[k] = sqrt(mean((test[,1]-pred)^2)) 
	
		}				
	
mean.error.bag=mean(unlist(error.bag))
sd.error.bag=sd(unlist(error.bag))
mean.error.bag
sd.error.bag

```


Observaciones:
1) Se puede encontrar más información en la página: http://www.inside-r.org/packages/cran/ipred/docs/bagging

2) Realizar las mismas simulaciones con el paquete caret, el paquete randomForest, el paquete adabag (mirar los helps)

#2-2 RANDOM FORESTS
Misma consigna que en 2.1, ahora con la función randomForest. 
Construya un modelo Random Forests con 100 árboles
A partir de la salida del modelo indique cuántas variables se 
utilizaron en cada partición y cuál es el porcentaje de varianza explicada por el modelo. 
Realice el gráfico de número de árboles vs el error cuadrático medio. 
Entiende que considerar 100 árboles para construir el bosque es adecuado? 
Obtenga el gráfico de importancia de variables.
Obtenga el error por muestra de prueba de la misma forma que se hizo antes.

```{r}
#?randomForest
rf=randomForest(MEDV~CRIM+ZN+CHAS+NOX+RM+DIS+TAX+PTRATIO,ntree=100,na.action=na.omit,data)
rf
plot(rf)
varImpPlot(rf)

#consideramos ahora 500 arboles
rf=randomForest(MEDV~CRIM+ZN+CHAS+NOX+RM+DIS+TAX+PTRATIO,ntree=500,na.action=na.omit,data)
rf
plot(rf)
varImpPlot(rf)

K=20
	error.rf=list(matrix(NA, K))
	n = nrow(data)
	for(k in 1:K)	{
		smp=sample(n,round(n/3))
		learn=data[-smp,]
		test=data[smp,]
			rf.lrn=randomForest(MEDV~CRIM+ZN+CHAS+NOX+RM+DIS+TAX+PTRATIO,ntree=500,na.action=na.omit,learn)
			pred=predict(rf.lrn,test)
			error.rf[k] = sqrt(mean((test[,1]-pred)^2)) 
	
		}				
	
mean.error.rf=mean(unlist(error.rf))
sd.error.rf=sd(unlist(error.rf))
mean.error.rf
sd.error.rf

errores=c(mean.error.tree,mean.error.bag,mean.error.rf)
sd=c(sd.error.tree,sd.error.bag,sd.error.rf)
tabla=cbind(errores,sd)
rownames(tabla) = c("CART","BAGGING","RF")
tabla


```




#3- CLASIFICACIÓN

#3-1 Boosting
```{r}
library(adabag)

data(iris)
iris.adaboost <- boosting(Species~., data=iris, boos=TRUE, 
                          mfinal=3)
importanceplot(iris.adaboost)
#Ver como se calcula esta importancia de variables.

```

#3-2 Comparación.
Implementar CART, RF y BOOSTING y comparar sus errores de clasificación
considerando un error por muestra de prueba con 20 iteraciones y los demás
valores por defecto que usamos para regresión.

```{r}
#Vamos a recodificar la variable MEDV

data$MEDVcat <- ifelse(data$MEDV > 20, c("high"), c("low")) 
summary(data)

MEDVcat=as.factor(data$MEDVcat)
data=data[,-1]
data=cbind(data,MEDVcat)
summary(data)

```

# 3-3 Curvas ROC

```{r}
library(rpart)  
library(rattle) # para data set
library(ROCR)

datos         <- weather
datos         <- within(datos, rm("Date","Location","RISK_MM")) #borra columnas dummy
set.seed(42)  #  fija la secuencia de numeros aleatorios
sampleTrain   <- sample(nrow(datos),(nrow(datos)*.6)) 
Train         <- datos[sampleTrain,]
Test          <- datos[-sampleTrain,]
 
modelo.rpart  <- rpart(RainTomorrow ~ .,Train, method="class")
 
 
# PREDICCION
#------------------------------------------------------------------------------
predict.rpart <- predict(modelo.rpart,Test,type = "prob")[,2] #prob. clase=yes
predict.rocr  <- prediction (predict.rpart,Test$RainTomorrow)
perf.rocr     <- performance(predict.rocr,"tpr","fpr") #True y False postivie.rate
 
 
# GRAFICO CURVA ROC
#------------------------------------------------------------------------------
auc <- as.numeric(performance(predict.rocr ,"auc")@y.values)
plot(perf.rocr,type='o', main = paste('Area Bajo la Curva =',round(auc,2)))  
#abline(a=0, b= 1)
 
```