analise.Rmd

---
title: "Vinhos"
author: "FIAP-06IA - Carlos Martinelli, Jônatas Bertolazzo, Letticia Nicoli, Renato Ramos"
date: "6/18/2019"
output: html_document
---

## Dataset: Qualidade dos Vinhos
Dataset realicionado com amostras de vinhos tintos e brancos do norte de Portugal.  
O objetivo é estimar a qualidade do vinho com base em suas características físico-químicas.

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
library(dplyr)
library(corrgram)
library(ggplot2)
library(rpart) 
library(rpart.plot)
library(plotly)
library(caret)
library(qcc)
library(ggplot2)
library(ggcorrplot)
library(scorecard)
library(cluster)
library(tclust)
library(outliers)
library(kableExtra)
red_color = "#E98383"
white_color = "#83E9D2"
cores = c(red_color, white_color)
```

```{r input, include=FALSE}
Vinhos <- read.csv2("/Users/letticianicoli/Desktop/Projects/Fiap/vinhos/BaseWine_Red_e_White.csv", row.names=1)

attach(Vinhos)
```

# Ánalise Exploratória

### Amostra do dataset
```{r echo=FALSE}
kable(head(Vinhos)) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

### Verifica se há valores ausentes (NAs)
```{r echo=TRUE}
sapply(Vinhos, function(x)all(is.na(x)))
```
Observamos que não há nenhum valor ausente/faltante nesse dataset. Dessa forma, eliminamos a necessidade de tratar esses valores.

### Verifica se há preditores com variância próxima de zero
```{r echo=TRUE}
nearZeroVar(Vinhos)
```
Não diagnosticamos preditores com um valor exclusivo (ou seja, preditores de variação zero) ou preditores com poucos valores exclusivos em relação ao número de amostras.

###  Resumo do dataset
```{r}
str(Vinhos)
```
O dataset possui 6497 observações de 13 variáveis. Sendo elas numéricas, inteiras e fatoriais.  

Dentre essas observações, podemos observar que a quantidade de `vinhos brancos` é **maior** do que `vinhos tintos`, conforme o gráfico abaixo:
```{r echo=TRUE}
barplot(table(Vinho), col=c(red_color, white_color))
```

Executando o summary:  
```{r}
summary(Vinhos)
```
Observamos que as variáveis `residualsugar`, `chlorides`, `freesulfurdioxide` e `totalsulfurdioxide` possuem valores muito distantes entre mínimo e máximo. Isso pode indicar outliers e/ou desequilíbrio entre a quantidade de vinhos tintos e brancos, o que pode interferir nos resultados de classificação.  

Nessa tabela podemos observar a combinação entre o tipo do vinho e a nota de qualidade:  
```{r echo=FALSE}
table(Vinho, quality)
```
* Com esse output, notamos que apenas `vinhos brancos` possuem a nota máxima em relação a qualidade.  
* Tanto para vinhos tintos quanto brancos, as notas 5 e 6 concentram a maior parte dos vinhos.

###  Histograma das variáveis

* Simétrico  
Contêm a partir do centro do gráfico o maior número de dados. 
```{r}
p_fixedacidity <- plot_ly(x = fixedacidity, type="histogram", name = "Fixed Acidity")
p_pH <- plot_ly(x = pH, type="histogram", name = "PH")
subplot(p_fixedacidity, p_pH, nrows=1)
```

* Assimétrico à direita  
Indica a ocorrência de altos valores com baixa frequência. 
```{r}
p_residualsugar <- plot_ly(x = residualsugar, type="histogram", name = "Residual Sugar")
p_chlorides <- plot_ly(x = chlorides, type="histogram", name = "Chlorides")
p_sulphates <- plot_ly(x = sulphates, type="histogram", name = "Sulphates")
p_volatileacidity <- plot_ly(x = volatileacidity, type="histogram", name = "Volatile Acidity")
p_freesulfurdioxide <- plot_ly(x = freesulfurdioxide, type="histogram", name = "Free Sulfur Dioxide")
subplot(p_residualsugar, p_chlorides, p_sulphates,p_volatileacidity,p_freesulfurdioxide, nrows=2)
```

* Assimétrico à esquerda  
A frequência dos dados está concentrada nos altos valores. 
```{r}
p_alcohol <- plot_ly(x = alcohol, type="histogram", name = "Alcohol")
subplot(p_alcohol, nrows=1)
```

* Bimodal  
Há o aparecimento de dois picos.
```{r}
p_totalsulfurdioxide <- plot_ly(x = totalsulfurdioxide, type="histogram", name = "Total Sulfur Dioxide")
subplot(p_totalsulfurdioxide, nrows=1)
```

* Multimodal  
Há o aparecimento de vários picos. 
```{r}
p_citricacid <- plot_ly(x = citricacid, type="histogram", name = "Citric Acid")
p_density <- plot_ly(x = density, type="histogram", name = "Density")
subplot(p_citricacid,p_density, nrows=1)
```

### Encontrando as variáveis mais relevantes 
```{r echo=TRUE}
library(randomForest)
rfModel = randomForest( Vinhos$quality ~ ., data=Vinhos, ntree=500 ) 
varImpPlot(rfModel) 
```

###  Ánalise de outliers

Abaixo temos o boxplot de todas as variáveis do dataset:  

```{r echo=FALSE}
par (mfrow=c(1,3))

boxplot(fixedacidity ~ Vinho, main='fixedacidity',col=cores)
boxplot(volatileacidity ~ Vinho , main='volatileacidity',col=cores)
boxplot(citricacid ~ Vinho, main='citricacid',col=cores)
```

```{r echo=FALSE}
par (mfrow=c(1,3))
boxplot(residualsugar ~ Vinho, main='residualsugar',col=cores)
boxplot(chlorides ~ Vinho, main='chlorides',col=cores)
boxplot(freesulfurdioxide ~ Vinho, main='freesulfurdioxide' ,col=cores)

```

```{r echo=FALSE}
par (mfrow=c(1,3))
boxplot(totalsulfurdioxide ~ Vinho, main='totalsulfurdioxide',col=cores)
boxplot(density ~ Vinho, main='density',col=cores)
boxplot(pH ~ Vinho, main='pH',col=cores)
```

```{r echo=FALSE}
par (mfrow=c(1,3))
boxplot(sulphates ~ Vinho, main='sulphates',col=cores)
boxplot(alcohol ~ Vinho, main='alcohol',col=cores)
boxplot(quality ~ Vinho, main='quality', col=cores)
```

Com os gráficos acima foi observado que todas as variáveis possuem candidatos a outliers.  
No dataset existem variáveis como `residualsugar` que possuem o valor máximo muito acima do terceiro quartil, isso pode gerar distorções nos algoritmos que serão executado a seguir. Além disso há uma concentração de outliers nessas variáriveis, o que traz a necessidade de removê-los para evitar distorção nos passos seguintes.  

Nivelamento  
```{r echo=TRUE}
VinhosOut <- Vinhos[Vinhos$quality > quantile(Vinhos$quality, .25) - 1.5*IQR(Vinhos$quality) & Vinhos$quality < quantile(Vinhos$quality, .75) + 1.5*IQR(Vinhos$quality), ]
```

###  Correlação
Para realizar a correlação transformamos o campo `Vinho` (fator) para tipo numérico. Além de *normalizar* os dados para evitar algum tipo de distorção.

```{r echo=TRUE}
VinhosOut$Vinho <- as.numeric(VinhosOut$Vinho)
norm_vinhos <- VinhosOut %>% mutate_at(c(1,2,3,4,5,6,7,8,9,10,11,13), list( ~ c(scale(.))))
```
Obs.: Dataset sem outliers.

```{r echo=FALSE}

matcor = cor(norm_vinhos)
ggcorrplot(matcor, hc.order = TRUE, 
           type = "lower", 
           lab = TRUE, 
           lab_size = 3, 
           method="circle", 
           colors = c("tomato2", "white", "springgreen3"), 
           title="Correlação dos Vinhos", 
           ggtheme=theme_dark)
```

Outra forma de visualização da correlação das variáveis:
```{r echo=FALSE}
corrgram(matcor, type = "cor", lower.panel = panel.shade, upper.panel = panel.pie)
```

Com o gráfico de correlação podemos observar alguns pontos:  
* O `alcool` tem uma anti-correlação alta com `density`  
* O `freesulfurdioxide` tem uma correlação possitiva alta com`totalsulfurdioxide`  
* O `volatileacidity` tem uma anti-correlação alta com o tipo do `vinho`  
* O `totalsulfurdioxide` tem uma correlação possitiva alta com o tipo do `vinho`  
* As outras correlações não são tão significativas no dataset, nesse momento.  

## Explicando a variável Quality

### Preparação para separação de dados: treinamento e teste  

Iremos separar a base de dados em 25% para testes e 75% para treino. Conforme output abaixo:
```{r echo=TRUE}
norm_vinhos = data.frame(norm_vinhos)
dt_list = split_df(norm_vinhos, ratio = 0.75, seed = 66)
train = dt_list$train
test = dt_list$test
```

Validando a consistência de qualidade entre as bases de treino e teste
```{r}
prop.table(table(train$quality))
```

```{r}
prop.table(table(test$quality))
```
As proporções estão, consideravelemente, bem distribuídas entre a qualidade. Dessa forma, conseguimos realizar um bom treinamento para o modelo.

### Regressão Linear

No primeiro modelo mantemos todas as varáveis:

```{r}

modelo1 <- lm(train$quality ~ train$fixedacidity+train$volatileacidity+
                train$citricacid+train$residualsugar+
                train$chlorides+train$freesulfurdioxide+train$totalsulfurdioxide+
                train$density+train$pH+train$sulphates+train$alcohol+
                train$Vinho)
summary(modelo1)

```

Nesse primeiro modelo a variável `citricacid` não possui muita relevância pois o seu p-value está próximo a 1.  
Para ter certeza sobre essa informação, utilizamos o método `stepwise` para identificar quais são as variáveis relavantes para o modelo.

```{r}
stepwise<-step(modelo1)
```

Confiança desse modelo é:
```{r}
confint(stepwise)
```

Modelo sem a variável de `citricidade`:
```{r}
modeloSemCitricidade <- lm(train$quality ~ train$fixedacidity + train$volatileacidity + train$residualsugar + train$chlorides + train$freesulfurdioxide + train$totalsulfurdioxide+ train$density + train$pH + train$sulphates + train$alcohol + train$Vinho)

summary(modeloSemCitricidade)  
```
O R-quadrado é de mais ou menos 30, o que significa que a regressão linear não descreve o modelo com tanta precisão.

Abaixo o modelo para base testes:

```{r echo=TRUE}
modeloTestSemCitricidade <- lm(test$quality ~ test$fixedacidity + test$volatileacidity + test$residualsugar + test$chlorides + test$freesulfurdioxide + test$totalsulfurdioxide+ test$density + test$pH + test$sulphates + test$alcohol + test$Vinho)

summary(modeloTestSemCitricidade) 
```

O modelo criado através da técnica de regressão linear não descreve muito bem a nota de qualidade dos vinhos, com uma acertividade de aproximadamente 30%. Não será necessário fazer o modelo de predição, devido ao baixo índice de acertividade.

### Análise de Resíduos

```{r echo=FALSE}
par(mfrow=c(2,2))
plot(modeloTestSemCitricidade)
```
Esse gráfico de resíduos apresenta a distância entre o valor estimado x valor real. Então, quanto mais próximo de zero o ponto estiver melhor é a assertividade do modelo.

Executa a predição do modelo:
```{r}
erro_usando_media <- mean((test$quality - mean(test$quality))^2)
sqrt(1-erro_usando_media)
```
O percentual de acerto do modelo é de 62,43%

Executando o teste de shapiro no modelo:

```{r}
shapiro.test(residuals(modeloTestSemCitricidade))
```
Como o p-value resultou em um valor menor que 0,05, não podemos assumir a normalidade.  
Pode-se dizer que temos um modelo de regressão com pouca assertividade.  
Isso corrobora com o valor d R-Quadrado demonstrando que **modelo não é assertivo**.

### Árvore de Regressão

A seguir será executado a árvore de regressão para comparar com a regressão linear.  
A árvore de regressão é um método de **aprendizado supervisionado** utilizado para classificação e regressão.  
A variável target é a quality e as variáveis que utilizaremos para prever o seu valor são: `fixedacidity`,  `volatileacidity`, `citricacid`, `residualsugar`, `chlorides`, `freesulfurdioxide`, `totalsulfurdioxide`, `density`, `pH`,`sulphates`,`alcohol`.


```{r include=FALSE}
attach(train)
```


```{r echo=TRUE}
arvore_regressao = rpart (quality ~ fixedacidity + volatileacidity + citricacid + residualsugar + chlorides + freesulfurdioxide + totalsulfurdioxide + density + pH + sulphates + alcohol,data=train, cp = 0.007, minsplit = 15, maxdepth=30)
```

```{r echo=TRUE}
rpart.plot(arvore_regressao, type=4, extra=1, under=FALSE, clip.right.labs=TRUE,
           fallen.leaves=FALSE,   digits=2, varlen=-10, faclen=20,
           cex=0.4, tweak=1.7,
           compress=TRUE, 
           snip=FALSE)
```
Para esse plot, conseguimos visualizr variáveis representativas como: `alcohol`, `volatileacidity`, `residualsugar`.  

### Testando a árvore de regressão

Erro utilizando o modelo de árvore de regressão:

```{r echo=TRUE}
Val_pred_tree = predict(arvore_regressao,interval = "prediction", level = 0.95) 

mse_tree = mean((quality - Val_pred_tree)^2)
sqrt(mse_tree)
```

Erro utilizando média:

```{r echo=TRUE}
erro_usando_media = mean((train$quality - mean(train$quality))^2)
sqrt(erro_usando_media)
```

Pode-se dizer que o modelo de árvore de regressão é mais acertivo que o modelo de regressão linear. E pode-se dizer que a árvore tem uma acertividade melhor do que informando apenas a média como explicação da qualidade.

E para **árvore de regressão** pode-se dizer que a quantide de `alcohol` é fundamental para a qualidade do vinho seguindo de `volatileacidity`.

### Classificando Vinhos em Bom ou Ruim considerando a variável Quality

Adicionando a coluna de classificação dos vinhos e já separando base em teste/treino:

```{r echo=TRUE}
vinhos_com_classificacao = VinhosOut
vinhos_com_classificacao$classificacao = ifelse(vinhos_com_classificacao$quality >= 6,  T, F)
dt_list_log = split_df(vinhos_com_classificacao, ratio = 0.75, seed = 66)
train_log = dt_list_log$train
test_log = dt_list_log$test
attach(train_log)
```

Executando o modelo com todos as variáveis:

```{r echo=TRUE}
modelo_logistico <- rpart (as.factor(classificacao) ~ fixedacidity+volatileacidity+citricacid+residualsugar+chlorides+freesulfurdioxide+totalsulfurdioxide+density+pH+sulphates+alcohol, maxdepth=20, train_log)
```

Resultado do modelo:

```{r}
rpart.plot(modelo_logistico, type=5, extra=104, under=FALSE, clip.right.labs=TRUE,
           fallen.leaves=FALSE,
           digits=2, varlen=-8, faclen=10,
           cex=0.6, tweak=1,
           compress=TRUE,
           snip=FALSE)
```

### Testando o modelo

#### Matriz de confusão

```{r echo=TRUE}
previsto.com.modelo<-predict(modelo_logistico, train_log, type='class')

matriz.de.confusao<-table(train_log$classificacao, previsto.com.modelo)
matriz.de.confusao
```

Calculando a diagonal da matriz

```{undefined echo=TRUE}
diagonal <- diag(matriz.de.confusao)
Acc <-  sum(diagonal)/sum(matriz.de.confusao)
Acc

```

### Executando o algoritmo para base de testes

```{r echo=TRUE}
previsto.valid<-predict(modelo_logistico, test_log , type='class')

test$previsto=previsto.valid
test$classificacao <- ifelse(test$quality >= 6, T, F)
test$errou = ifelse(test$previsto != test$classificacao, 1, 0)

```

Matriz de confusão para a base teste

```{r echo=TRUE}
previsto.com.modelo<-predict(modelo_logistico, test_log, type='class')

matriz.de.confusao<-table(test_log$classificacao, previsto.com.modelo)
matriz.de.confusao
```

Calculando a diagonal da matriz

```{undefined echo=TRUE}
diagonal <- diag(matriz.de.confusao)
Acc <-  sum(diagonal)/sum(matriz.de.confusao)
Acc

```

### Clusterizando os Vinhos

Como **técnica não supervisionada**, vamos testar se o algoritmo de clusterização será adequado para agrupar dois conjunto de vinhos, categorizando-os como vinhos bons e vinhos ruins. 

A variável `quality`, que identifica a nota do vinho, será a variável utilizada para correlacionar com as demais variáveis para identificar se existe algum agrupamento entre os vinhos.

De acordo com Luis Costa de Oliveira, Sara Oliveira, Maria Eugenia em seu artigo 'Avaliação das características físico-químicas e colorimétricas de
vinhos finos', a cor não é uma característica físico-química. Portanto foi removida a coluna tipo de vinho.

```{r}
vinhos_noColor <- norm_vinhos[1:12]
kable(head(vinhos_noColor)) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

### Método hierárquico

```{r echo=TRUE}
hier_cluster<-hclust(dist(vinhos_noColor),method='ward.D2')
d <- dist(norm_vinhos, method = "euclidean") 
plot(hier_cluster, ylab='distancia', cex=0.6)

groups <- cutree(hier_cluster, k=4)
rect.hclust(hier_cluster, k=4, border="red") 
```
* Divisão por 4 clusters.

### Método não hierárquico

Utilizando o K-means para descobrir a quantidade ideal de clusters dentro de 10 iterações.

```{r echo=TRUE}
set.seed(45)
wss = 0
for (i in 1:10) {
  wine_cluster <- kmeans(vinhos_noColor, centers = i)
  wss[i] <- wine_cluster$tot.withinss
}
plot(1:10, wss, type = "b",  main="Elbow method",
     xlab = "Number of Clusters", 
     ylab = "Within groups sum of squares",pch=8, col="red")
```

Com base na plotagem, podemos determinar que, após o cluster 3, não vemos uma grande queda na soma das distâncias quadradas dentro de cada cluster, portanto, podemos considerar o valor de K como 3 e prosseguir com o agrupamento.

### K-Means Segmentation
```{r}
set.seed(45)
k=3
w_cluster = kmeans(vinhos_noColor,  centers = k, nstart=20)
w_cluster
```

```{r}
segmento<-w_cluster$cluster
table (segmento)
```

Quantidade de interações até chegar nos clusters:
```{r}
Qte_iter<-w_cluster$iter
Qte_iter
```

```{r}
kable(head(aggregate(vinhos_noColor,by=list(segmento),FUN=mean))) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
```

```{r}
clus_teste <- tkmeans(vinhos_noColor , k = 3, alpha = 0.01)
plot(clus_teste)
```


```{r}
clusplot(vinhos_noColor, clus_teste$cluster, color=TRUE, shade=TRUE,
         labels=2, lines=0 , cex=0.75)
```

# Componentes Principais (PCA)

Com a dataset sem o tipo do vinho foi testado o método dos Componentes Principais:

```{r echo=TRUE}
acpcor <- prcomp(vinhos_noColor, scale = TRUE) 
summary(acpcor)
```

```{r}
plot(1:ncol(vinhos_noColor), acpcor$sdev^2, type = "b", xlab = "Componente",
     ylab = "Variância", pch = 20, cex.axis = 0.8, cex.lab = 0.8)
```

Pode-se observar que o número de componentes prinpais que explicam a maior parte dos componentes é 3, por tanto, decidiu-se por utilizar 3 componentes.

### Exibição de cada componente
```{r}

par (mfrow=c(2,3))
escore1 <- acpcor$x[, 1]
hist(escore1)

escore2 <- acpcor$x[, 2]
hist(escore2)

escore3 <- acpcor$x[, 3]
hist(escore3)

escore4 <- acpcor$x[, 4]
hist(escore4)

escore5 <- acpcor$x[, 5]
hist(escore5)

par (mfrow=c(1,1))
```

Criando um novo dataset com os componentes princpais
```{r echo=TRUE}
vinhos_cpa <-cbind(escore1,escore2, escore3,escore4 ,escore5)
```


### Utilizando o método hierárquico com PCA

```{r}
hier_cluster<-hclust(dist(vinhos_cpa),method='ward.D2')
d <- dist(vinhos_cpa, method = "euclidean") # distance matrix
plot(hier_cluster, ylab='distancia', cex=0.6)
```


* Observando o dendograma podemos cortar em 3 grandes grupos:

```{r}
hier_cluster<-hclust(dist(vinhos_cpa),method='ward.D2')
d <- dist(vinhos_cpa, method = "euclidean") # distance matrix
plot(hier_cluster, ylab='distancia', cex=0.6)

groups <- cutree(hier_cluster, k=3) # cut tree into 5 clusters
# draw dendogram with red borders around the 5 clusters
rect.hclust(hier_cluster, k=3, border="blue") 
```

### Utilizando o método não hierárquico com PCA

Determinando a quantidade de clusters

```{r}
set.seed(45)
wss <- (nrow(vinhos_cpa )-1)*sum(apply(vinhos_cpa ,2,var))
for (i in 2:15) wss[i] <- sum(kmeans(vinhos_cpa ,iter.max=100,
                                     centers=i)$withinss)
plot(1:15, wss, type="b", xlab="Number of Clusters",
     ylab="Within groups sum of squares") 
```

Com base na plotagem, podemos determinar que, após o cluster 3, não vemos uma grande queda na soma das distâncias quadradas dentro de cada cluster, portanto, podemos considerar o valor de K como 3 e prosseguir com o agrupamento.

```{r echo=TRUE}
set.seed(333)
output_cluster<-kmeans(vinhos_cpa,3,iter=100)

clus_vinhos_cpa<-output_cluster$cluster
table (clus_vinhos_cpa)
```

Assim os vinhos ficam agrupados da seguinte maneira

```{r}
kable(head(aggregate(norm_vinhos,by=list(clus_vinhos_cpa),FUN=mean))) %>%
  kable_styling(bootstrap_options = c("striped", "hover"))

```