Deliverable5.Rmd

---
title: "Forecasting modeling of categorical target Assignment"
author: "Katerina Dimitrova, Jose Romero, Sergi Munoz"
date: "5 June, 2018"
output: pdf_document
editor_options: 
  chunk_output_type: console
---
  
  ```{r setup, include=FALSE}
#knitr::opts_chunk$set(echo = TRUE)
```

```{r,include=FALSE}
setwd("C://Users/Sergi/Desktop/Sergi/ADEI") #Change 
#setwd("D:/ADEI/ADEI.git/trunk") #Change

```

#Previous work

## Load requiered packages

```{r, echo=FALSE, include=FALSE}
rm(list=ls())
# Load Required Packages: to be increased over the course

requiredPackages <- c("effects","FactoMineR","car", "factoextra","RColorBrewer","ggplot2","missMDA","mvoutlier","dplyr","caret", "ggmap","ggthemes","knitr","MVA", "ROCR")
missingPackages <- requiredPackages[!(requiredPackages %in% installed.packages()[,"Package"])]

if(length(missingPackages)) install.packages(missingPackages)
lapply(requiredPackages, require, character.only = TRUE)

# Useful function
calcQ <- function(x) {
  s.x <- summary(x)
  iqr<-s.x[5]-s.x[2]
  list(souti=s.x[2]-3*iqr, mouti=s.x[2]-1.5*iqr, min=s.x[1], q1=s.x[2], q2=s.x[3],        q3=s.x[5], max=s.x[6], mouts=s.x[5]+1.5*iqr, souts=s.x[5]+3*iqr ) }

```

# Statistical Modelling

## Load your sample after data cleaning and validation


```{r}
#load("Taxi5000_raw_DataDefinitivev1.RData")
load("C:/Users/Sergi/Desktop/Sergi/ADEI/Taxi5000_raw_DataDefinitivev1.RData")

names(df)
#summary(df)
#vars_con
#vars_dis
#vars_res


```
## Split data 70-30

In order to build the model, and test it afterwords, we will generate work and test samples (consisting on a 70-30 split).

```{r}
llwork<-sample(1:nrow(df),0.70*nrow(df),replace=FALSE)
llwork<-sort(llwork);length(llwork)
dfwork<-df[llwork,]
dftest<-df[-llwork,]
```

## Selecting the best model

We start modelling and executing some basic tests of Bayesian information criterion (deviance) and colinearity to achieve the better suited model for our binary target.
In first place, we take a look at catdes and cor outputs to grasp those numerical and discrete coefficients which explain better our targey Any Tip. In vars_cexp we keep the continuous explicative variables to look at their correlation table so we can decide wether to try to combine them or avoid it in order to don't carry colinearity in our model.

```{r}

names(dfwork)
catdes(df,num.var=22)
vars_cexp<-vars_con[c(1:4,5,6,7,8,11,13)]
dfX <- dfwork[,c("AnyTip",vars_cexp)]
cor(dfwork[,c(vars_cexp)])

m50<-glm(AnyTip~.,family="binomial",data=dfX)
summary(m50)
Anova(m50,test="Wald")
vif(m50)

m51<-step(m50,k=log(nrow(dfwork)))
summary(m51)
anova (m51,m50, test = "Chisq")
Anova(m51,test="Wald")
vif(m51)

```

So far we took general model m50 and applied step on it. The resulting one is the addition of Pickup_latitude + Dropoff_longitude + Trip_distance.
We want to consider adding factors to our model. For that purpose we get those appearing in catdes as associated to anytip. 

```{r}
#We try claHP 
m52<-glm(AnyTip~ (Pickup_latitude + Dropoff_longitude + Trip_distance)*claHP ,family="binomial",data=dfwork)
summary(m52)
#exaggerated colinearity with cluster factor
vif(m52)

#As this last model has quite big complexity, we apply step on it and it gets rid of claHP and any possible interaction. So the same m51 model with only numerical coefficients.
m62<-step(m52,k=log(nrow(dfwork)))
summary(m62)

#We try again without trip distance 
m53<-glm(AnyTip~ (Pickup_latitude + Dropoff_longitude)*claHP ,family="binomial",data=dfwork)
summary(m53)
vif(m53)

#And we apply again step on it, resulting Dropoff_longitude + claHP
m63<-step(m53,k=log(nrow(dfwork)))
summary(m63)
vif(m63)

#This two models are significally differents (p.value under threshold).
anova(m53,m63, test = "Chisq")

#Now we try to add a new factor to the model, even though we know it adds too much complexity and we will need to simpify it through step method.
m54<-glm(AnyTip~ (Pickup_latitude + Dropoff_longitude)+(claHP*f.distance) ,family="binomial",data=dfwork)
summary(m54)
#Too much colinearity claHP
vif(m54)

m64<-step(m54,k=log(nrow(dfwork)))
#Dropoff_longitude + claHP + f.distance
summary(m64)
#we would choose m64

#but colinearity test points that the model is not good enough (claHP over 4 points)
vif(m64)
```

From this moment, we try to simplify our interactions with factors to achieve some balanced model. 
```{r}

m55<-glm(AnyTip~ Pickup_latitude + claHP + f.distance + Pickup_latitude:claHP,family="binomial",data=dfwork)
summary(m55)
vif(m55)

#Here, with this step, it seems to get some decent model taking into consideration just and addition of factors clusters and f.distance. Although, it contains coefficient (groups of the factors) which are irrelevant. 
m65<-step(m55,k=log(nrow(dfwork)))
summary(m65)
vif(m65)

# We change to two numerical again and let's see which results we have from this interaction.
m56<-glm(AnyTip~ Pickup_latitude + Dropoff_longitude * f.distance, family="binomial",data=dfwork)
summary(m56)
vif(m56)

#Since m56 it's not fitting well, we try a new step. m66 seems interesting, all their coefficients are significative somehow and it has no colinearity.
m66<-step(m56,k=log(nrow(dfwork)))
summary(m66)
vif(m66)

# With all models built so far we got the best estimator on m51 and m62 (which are actually the same).
BIC(m51,m52,m53,m54,m55,m56,m62,m63,m64,m65,m66)


m58<-glm(AnyTip~Dropoff_longitude:claHP + Dropoff_longitude:f.distance, family="binomial",data=dfwork)
summary(m58)
vif(m58)
BIC(m51,m58)
```

We decide to incorporate f.total as a coefficient
```{r}

m59<-glm(AnyTip~ Dropoff_longitude + Pickup_latitude + f.total + f.distance + claHP, family="binomial",data=dfwork)
summary(m59)
vif(m59)
BIC(m51,m59)

#We apply step in order to simplify
m69<-step(m59,k=log(nrow(dfwork)))
summary(m69)
#colinearity
vif(m69)

m1<-glm(AnyTip~ (Dropoff_longitude + Trip_distance)*f.total, family="binomial",data=dfwork)
summary(m1)
vif(m1)

m11<-step(m1,k=log(nrow(dfwork)))
summary(m11)
#still colinearity present
vif(m11)

#This seems to be also a good model
m2<-glm(AnyTip~ Dropoff_longitude + f.espeed + Trip_distance+ f.total, family="binomial",data=dfwork)
summary(m2)
vif(m2)
Anova(m2)
BIC(m2)

BIC(m51,m2)


```
We finally choose m2.


## Diagnostics

### Influent Data

In the plots below one big element can be clearly distinguished over the rest: 759257 (in the boxplot 1844 but is the same element).

```{r}
# Potentially influent data (we don't see any from hatvalues analysis)
quantile(hatvalues(m2),seq(0,1,0.1))
mean(hatvalues(m2))
hh<-5*mean(hatvalues(m2))
llhat<-which(hatvalues(m2)>hh);length(llhat)
hh
Boxplot(hatvalues(m2),labels=rownames(dfwork))
abline(h=hh,lwd=2,col="red",lty=2)
dfwork[llhat,]

ll1 <- influencePlot(m2); ll1
Boxplot(cooks.distance(m2))
Boxplot(rstudent(m2),id.n=15)

```

So we took a look at this very differentiate element. We can observe that the most influential element (759257) was outlier in the second dimension of PCA too.

```{r}
llcoo<-which(cooks.distance(m2)>0.01);length(llcoo)
dfwork[llcoo[1],]
```

We decide to remove the outlier.

```{r}

df1<-dfwork[-llcoo[2],]

```

### Predicting

Finally, we execute the predictions and generate a confusion matrix with the results of it.
To conclude, the diagonal of this confusion matrix shows a performance of 70% hit rate in our prediction model against our testing data set, while the model 0 without any coefficient would give us only a 41% of hit rate.

Curve ROC of our choosen method is displayed at the end.

```{r}
tfit2<-predict(m2,type="response",newdata=dftest)
fit.AnyTip<-factor(ifelse(predict(m2,type="response")<0.5,0,1),labels=c("fit.No","fit.Yes"))
tt<-table(fit.AnyTip,dfwork$AnyTip);tt
sum(tt)
100*sum(diag(tt))/sum(tt)

m0<-glm(AnyTip~1, family="binomial", data=dfwork)
fit0<-predict(m0,type="response")
fit.AnyTip0<-factor(ifelse(fit0<0.5,0,1),labels=c("fit.Yes"))
tt0<-table(fit.AnyTip0,dfwork$AnyTip);tt0;sum(tt0)
100*sum(tt0[1,2])/sum(tt0)

# Plot ROC curve
library("ROCR")
dadesroc<-prediction(predict(m2,type="response"),dfwork$AnyTip)
par(mfrow=c(1,2))
plot(performance(dadesroc,"err"))
plot(performance(dadesroc,"tpr","fpr"))
abline(0,1,lty=2)

```