Aula ibpad

Aula_4.1_modelagem.Rmd

---
title: "Aula 4.1 - Modelagem"
output:
  html_notebook:
    css: custom.css
    number_sections: yes
    toc: yes
    toc_float: true
---

```{r}
library(tidyverse)
library(ggExtra)
library(ggalt)
library(broom)
library(GGally)
library(rpart)
library(rpart.plot)
```


**AVISO**

Esta aula não substitui um bom curso de modelagem ou de Machine Learning. Para isso, sugiro dois recursos:

![](https://images-na.ssl-images-amazon.com/images/I/41sEbo-RgjL._SY346_.jpg)
![](https://images-na.ssl-images-amazon.com/images/I/41jy3mwLy9L.jpg)


A parte de modelagem deste curso se destina a apresentar uma introdução de conceitos de aprendizagem supervisionada (regressão e classificação) e não-supervisionada (clusterização).


# Inteligência Artificial e Machine Learning

A Inteligência Artificial é uma área de estudos da computação que se interessa pelo estudo e criação de sistemas que possam exibir um comportamento inteligente realizar tarefas complexas com um nível de competência que é equivalente ou superior ao de um especialista humano.

Machine Learning é um subcampo da Inteligência Artificial que permite dar aos computadores a habilidade de aprender sem que sejam explicitamente programados para isso. É um método de análise de dados que automatiza o desenvolvimento de modelos analíticos. 


Existem dois tipos de problemas de Machine Learning:

* Aprendizado supervisionado: Caracterizado por um conjunto de variáveis explanatórias ($X$) e uma variável resposta ($y$). Um algoritmo é então aplicado para aprender uma função `f` tal que $Y = f(x)$ seja uma aproximação aceitável de acordo com uma métrica de erro ou acurácia. Se $y$ for uma variável numérica, tem-se um modelo de regressão, caso contrário se trata de um problema de classificação.  

* Aprendizado não-supervisionado: existem apenas as variáveis de input $X$, não há um output desejado. O objetivo é modelar a estrutura dos dados para encontrar padrões ou grupos. Isso é chamado de não-supervisionado porque não há uma resposta correta. A interpretação dos resultados é mais complexa e depende do feeling humano. Uma aplicação comum são os algoritmos de clusterização, como **k-means**, que divide um conjunto de dados em **k** grupos ou clusteres, que são homogêneos ou semelhantes entre si.

## Aprendizado supervisionado

### Regressão

A técnica chamada de regressão é usada para predizer o valor de uma variável Y (chamada de variável resposta ou dependente) baseado em uma ou mais variáveis X (variável explanatória ou independente). Se a regressão utiliza apenas uma variável explanatória, é chamada de regressão simples. O objetivo da regressão é representar a relação entre as variáveis resposta e explanatória por meio de uma equação matemática linear do tipo:

$Y = \beta_1 + \beta_2X + \epsilon$

onde $\beta_1$ é a interceptação da reta com o eixo vertical e $\beta2$ o coeficiente de inclinação associado à variável explanatória. Tais elementos são chamados coeficientes da regressão. O termo $\epsilon$ representa o termo do erro, que é a parte de Y que a regressão é incapaz de explicar (por existir outras variáveis que explicariam Y mas que não foram incorporadas ao modelo).

Neste capítulo, usaremos como exemplo o dataset `x15` extraído [deste site](http://people.sc.fsu.edu/~jburkardt/datasets/regression/), que traz dados que tentam explicar o consumo de petróleo baseado em outros inputs. Os dados não vêm já formatados, o que vai acabar servindo para mostrar um exemplo de limpeza de dados em R.



```{r regressao-ch02}

url <- "http://people.sc.fsu.edu/~jburkardt/datasets/regression/x15.txt"
readLines(url)
# Os dados de verdade estão presentes entre as linhas 42 a 89. 
# Os nomes das colunas estão presentes entre as linhas 36 a 41
dados <- readLines(url)[42:89]
# como salvamos o dataset como se fosse um string em uma variavel...
# (ao inves de importar de um arquivo de texto), precisamos indicar que a variavel...
# dados corresponde a um arquivo de texto que contem um dataset
con <- textConnection(dados)
df <- read.table(con, header = FALSE)
close(con)
# adicionando o nome das colunas
names(df) <- c("linha", "imposto", "renda_media", "km_asfaltado", "n_carteiras",
               "consumo")
# retirar a primeira coluna já que nao serve para nada
df <- df[,-1]

head(df)

```

O objetivo aqui é descobrir qual das 4 variáveis explica a variável `consumo`.

#### Correlação

Correlação é um indicador estatístico que mede o nível de dependência linear entre duas variáveis. Está definida no intervalo [-1, +1]. Se a correlação é negativa, indica que as variáveis são inversamente proporcinais: quando uma aumenta, a outra diminui. Se é positiva, indica que as variáveis são diretamente proporcionais.

Medir a correlação no R é muito simples:

```{r}
# Usando a função cor
cor(df$consumo, df$n_carteiras)
cor(df$consumo, df$renda_media)


```
Para analisar todas as correlações entre todas as variáveis de uma só vez, podemos tanto calcular uma matriz de correlação como plotar todas as correlações com o auxílio do pacote `GGally`:

```{r}
# matriz de correlacao:
cor(df)
# usando o pacote GGally
library(GGally)
ggpairs(df)

```

Percebe-se pela matriz de correlação (e principalmente pelo gráfico) que só valeria a pena usar como variáveis explanatórias do nosso objeto de estudo o número de carteiras de habilitação e talvez o imposto sobre o consumo de petróleo.

#### Análise explatória e gráfica

No gráfico acima, mais precisamente na parte referente ao gráfico de dispersão entre as variáveis consumo e número de carteiras, é notório que existem quatro pontos que se destacam dos demais: três para cima e um para baixo. Como eles estão distantes verticalmente do "bolo" de pontos, o que é um indício de que são outliers em consumo. Tratamento de outliers é um tópico fundamental na construção de modelos.

Primeiramente, vamos analisar a distribuição da variável de consumo por meio de um histograma:

```{r}
# pelo ggplot2, o grafico fica com intervalos meio feios:
ggplot(df, aes(x = consumo)) +
  geom_histogram() +
  labs(x = "Consumo", y = "Quantidade")

```

O histograma confirma que os valores acima de 800 peças e abaixo de 400 são anomalias dada a distribuição normal que a variável `consumo` aparenta ter. 

Visto que é possível descrever `consumo` como uma distribuição normal, pode-se assumir que a probabilidade de que um valor esteja fora do intervalo $\mu(consumo) \pm 3 \times \sigma(consumo)$ é de apenas 0,27%. Vamos destacar esses outliers com o auxílio dos pacotes `ggalt` e `ggExtra`, que são extensões do `ggplot2`.

```{r}
# calcular limites superior e inferior
lim_sup <- mean(df$consumo) + 2 * sd(df$consumo)
lim_inf <- mean(df$consumo) - 2 * sd(df$consumo)

# por curiosidade, como ficaria com o pipe:
lim_sup <- df$consumo %>% {mean(.) + 3 * sd(.)}
lim_inf <- df$consumo %>% {mean(.) - 3 * sd(.)}
# criar dataframe sem outliers
df_sem_out <- df %>% filter(lim_inf < consumo & consumo < lim_sup)


p <- ggplot(df, aes(x = n_carteiras, y = consumo)) +
  geom_point() + 
  # adicionar curva da regressao
  geom_smooth(method = "lm") +
  # plotar circulo que deixe de fora os outliers
  geom_encircle(data = df_sem_out,
                aes(x = n_carteiras, y = consumo), color = 'red') + 
  labs(x = "Número de carteiras de motorista", y = "Consumo de petróleo")

# plotar histogramas nas margens do grafico
p_marg <- ggMarginal(p, type = "histogram")
p_marg
```

```{r, echo = FALSE}
grid::grid.newpage()
grid::grid.draw(p_marg)

```


Agora o outlier onde `consumo > 800` e `n_carteiras > 0,65` ficou mais visível. Portanto, faz sentido removê-lo da análise.

Como ficam as correlacões sem o outlier?

```{r}
# comparar as matrizes de correlacao
cor(df)
cor(df_sem_out)

```

#### Modelagem por regressão simples

No R, é bem simples ajustar um modelo de regressão. Usando a variável `n_carteiras` como explanatória e `consumo` como resposta, um modelo é construído da seguinte maneira:

```{r}
modelo.simples <- lm(consumo ~ n_carteiras, data = df_sem_out)
summary(modelo.simples)


```

Com o modelo criado, é possível descrever a relação entre `consumo` e `n_carteiras` matematicamente por meio da seguinte equação:  

$consumo = -123.5 + 1217.8 \times n\_carteiras$

Vamos deixar para analisar os diagnósticos da regressão no próximo item:

#### Regressão multivariada

Suponha também que você deseja incorporar a variável `imposto` ao modelo. Antes de fazer isso, vamos ver como as duas variáveis explanatórias se relacionam com a resposta em um gráfico tridimensional:

```{r, fig.height=5}

ggplot(df, aes(x = n_carteiras, y = consumo, color = imposto)) +
  geom_point() +
  scale_color_continuous(low = "green",  high = "red")

```

Valores altos de impostos aparentam estar associados com valores baixos de consumo.

Para adicionar uma nova variável ao modelo, fazemos:

```{r}
modelo.mult <- lm(consumo ~ n_carteiras + imposto, data = df_sem_out)
summary(modelo.mult)
# Usando o pacote broom para formatar o output dos modelos de regressao
# concatenando os dois modelos em um dataframe so

# metricas dos regressores
modelo.simples %>% tidy
modelo.mult %>% tidy

# metricas do modelo
rbind(modelo.simples %>% glance,
     modelo.mult %>% glance)



```

Agora vamos à análise dos indicadores da regressão:

##### Hipótese nula da regressão

A presença de um valor-p indica que existe uma hipótese nula sendo testada. Na regressão linear, a hipótese nula é a de que os coeficientes das variáveis explanatórias são iguais a zero. A hipótese alternativa é a de que os coeficientes não são iguais a zero, ou seja, existe uma relação matemático entre as variáveis do modelo.

##### valor-p

Nós podemos considerar um modelo linear estatisticamente significante apenas se os valores-p, tanto dos coeficientes como do modelo, são menores que um nível de significância pré-determinado, que idealmente é 0,05.

##### valor-t ou estatística

Um valor grande do valor-t indica que é pouco provável que os coeficientes não sejam iguais a zero puramente por coincidência. Assim, quanto maior o valor-t, melhor.

##### R-quadrado e R-quadrado ajustado

R-quadrado é a proporção da variação da variável resposta que é explicada pelo modelo. Quanto maior, melhor o modelo, supostamente.

Se continuarmos adicionando variáveis ao modelo de regressão, o R-quadrado apenas tende a crescer, intuitivamente. Isso acontecerá mesmo que a variável explanatória adicionada não seja significante. Para evitar esse problema que tornaria a comparação entre modelos praticamente inviável, o R-quadrado ajustado "penaliza" o valor do R-quadrado pelo número de variáveis adicionadas. Semelhantemente ao R-quadrado, quanto maior, melhor.


##### Análise dos resíduos

Um indicador visual da qualidade de um modelo é a distribuição dos modelos: um bom modelo apresentará resídos que seguem uma distribuição normal com média 0.

Um modelo de regressão pressupõe que seus resíduos (subtração entre o valor real e o ajustado) seguem uma distribuição normal e não possuem nenhum tipo de relação matemática com os regressores do modelo (ou mesmo com variáveis independentes não usadas no modelo).

Para analisar o primeiro pressuposto do modelo múltiplo, fazemos um histograma dos resíduos:
```{r}
modelo.mult$residuals %>% hist
```

Para analisar o segundo pressuposto, plotamos os resíduos contra todas as variáveis do dataset:

```{r}
par(mfrow=c(2,2))
plot(df_sem_out$imposto, residuals(modelo.mult), xlab = "imposto")
plot(df_sem_out$renda_media, residuals(modelo.mult), xlab = "renda media")
plot(df_sem_out$km_asfaltado, residuals(modelo.mult), xlab = "km asfaltado")
plot(df_sem_out$n_carteiras, residuals(modelo.mult), xlab = "n_carteiras")


```

Os resíduos do modelo múltiplo aparentam ter alguma relação com a variável `renda_media`, o que indica que ela pode ser incorporada ao modelo.


#### Regressão como modelo preditivo {#modelo-preditivo}

Um dos objetivos da regressão, além de descrever matematicamente a relação entre duas ou mais variáveis, é prever o valor da variável dependente baseado em novos valores da variável independente. Não é possível afirmar que um modelo apresentará um bom desempenho preditivo analisando apenas as métricas da regressão do tópico anterior. É necessário testar o modelo em dados que ele nunca viu.

A prática comum em Data Science é de separar o conjunto de dados que se tem em mãos em dois conjuntos: o de treino, que será usado para construir o modelo, e o de teste, que será usado como input do modelo para avaliar sua acurácia.

Após obter as previsões, deve-se usar uma ou mais métricas de erro (ver capítulo posterior) para avaliar a precisão do modelo.

```{r}

set.seed(1993)  # escolhendo uma seed para tornar os resultados previsiveis
indice <- sample(1:nrow(df_sem_out), 0.8*nrow(df_sem_out))  # row indices for training data
treino <- df_sem_out[indice, ]  # model training data
teste  <- df_sem_out[-indice, ]   # test data

# construindo os dois modelos, mantendo o teste de fora
modelo.simples <- lm(consumo ~ n_carteiras, data = treino)
modelo.mult <- lm(consumo ~ n_carteiras + imposto, data = treino)

# calcular previsao baseado no dataframe de teste
prev.simples <- predict(modelo.simples, teste)
prev.mult <- predict(modelo.mult, teste)
# uma das metricas é correlação entre previsto e real:
real = teste$consumo

cor(prev.simples, real)
cor(prev.mult, real)
# outra metrica é o MAPE
mean(abs(prev.simples - real)/real)
mean(abs(prev.mult - real)/real)

```
Os dois modelos apresentam resultados semelhantes de erro. Portanto, pelo menos para este teste, não houve um aumento significativo de acurácia no modelo ao incorporar a variável `imposto` como explanatória.

### Classificação

Suponha que estamos modelando reclamações de consumidores sobre serviços prestados por empresas no Brasil. Desejamos obter, baseado em variáveis a serem coletadas, se um dado processo aberto por um cliente será resolvido ou não resolvido.  Veja que $Y$ pode assumir apenas esses dois valores (partindo do pressuposto que nenhuma reclamação fica com o Status *Em andamento* para sempre), sendo assim definida como uma variável binária.

Em Machine Learning, esse tipo de problema é chamado de **Classificação Supervisionada**, que consiste em identificar em qual valor de uma variável categórica uma nova observação pertence. Um exemplo clássico de um problema de Classificação é tentar designar se um dado e-mail deve ser identificado como spam ou não. Outras boas referências de aplicação de técnicas de Classificação são:

* ["Why Should I Trust You?": Explaining the Predictions of Any Classifier](https://arxiv.org/abs/1602.04938): Ribeiro et al. criaram um método de explicar as variáveis explanatórias de qualquer modelo de Classificação, algo fundamental para aumentar a confiança nos resultados do modelo;  
* [A machine learning approach to automatic music genre classification](http://www.scielo.br/scielo.php?pid=S0104-65002008000300002&script=sci_arttext): Neste trabalho, diversos algoritmos de classificação são testados

Existem dezenas de modelos de classificação. Um que eu gosto bastante e que é a base para modelos mais avançados é a árvore de decisão (Decision Trees)

Algoritmos baseados em Decision Trees tentam achar maneiras de criar subsets ou subgrupos do universo dos dados, onde cada subgrupo pertence a um *node*. O objetivo do modelo é criar *nodes* onde haja uma distinção clara entre as classes previstas de forma que possa a cada *node* a probabilidade de um indivíduo pertencer a uma classe. O gráfico abaixo é um exemplo simples e didática de uma árvore de decisão:

[![](https://3qeqpr26caki16dnhd19sv6by6v-wpengine.netdna-ssl.com/wp-content/uploads/2016/02/Example-Decision-Tree.png)](https://machinelearningmastery.com/classification-and-regression-trees-for-machine-learning/)

Nesse modelo, que tenta prever o sexo de uma pessoa baseada na altura e no peso, o algoritmo de classificação funciona como uma série de regras **SE-NÃO**:  

* Se a altura for maior que 180cm, o indivíduo é um homem;  
* Se a altura é menor ou igual a 180cm e o peso é maior que 80kg, o indivíduo é homem;  
* Caso contrário, o indivíduo é mulher.  

No R, os modelos de Decision Trees são aplicados principalmente pelo pacote `rpart`.

#### Construindo uma árvore de decisão

https://rpubs.com/H_Zhu/235617

Como exemplo, criaremos um modelo para identificar se um clone do universo Star Wars, baseado em algumas de suas características, é Apto ou defeituoso.

```{r}
# url <- "https://archive.ics.uci.edu/ml/machine-learning-databases/adult/adult.data"
# 
# adult <- read.table(url, sep = ',', fill = FALSE, strip.white = TRUE)
# 
# colnames(adult) <- c('age', 'workclass', 'fnlwgt', 'education', 
#                      'educatoin_num', 'marital_status', 'occupation',
#                      'relationship', 'race', 'sex', 
#                      'capital_gain', 'capital_loss', 'hours_per_week',
#                      'native_country', 'income')
# glimpse(adult)
#library(titanic)
url <- "https://raw.githubusercontent.com/sillasgonzaga/paixaopordados-blogdown/master/data/clones.csv"
df_clones <- read.csv2(url)
glimpse(df_clones)
```

Temos então uma variável resposta ´status´ e 3 possíveis variáveis preditoras. Antes de construir qualquer modelo, é fundamental realizar uma análise exploratória:

Status ~ General Jedi

```{r}
table(df_clones$general_jedi_encarregado, df_clones$status)
```

Status ~ Estatura

```{r}
df_clones %>% 
  ggplot(aes(y = estatura_cm, x = status)) + 
  geom_boxplot()

```

Status ~ Massa

```{r}
df_clones %>% 
  ggplot(aes(y = massa_em_kilos, x = status)) + 
  geom_boxplot()
```

Qual o modelo baseline? Qual o resultado de acurácia mínimo que o modelo deve ter?

```{r}
table(df_clones$status)/nrow(df_clones)
```

Criação do modelo:

```{r}
# separacao de conjuntos de treino e teste
n = nrow(df_clones)
set.seed(123)
ind_treino <- sample(1:n, size = n * 0.7)

tr <- df_clones[ind_treino, ]
te <- df_clones[-ind_treino, ]
```


```{r}
# construcao do modelo
mod <- rpart(status ~ ., data = tr)

rpart.plot(mod)

prp(mod, type = 4, extra = 6, fallen.leaves = FALSE, varlen = 0,
    faclen = 0, box.palette = "auto")
```

Vamos avaliar a acurácia do modelo:

```{r}
ycast <- predict(mod, te, type = "class")
yreal <- te$status

caret::confusionMatrix(data = ycast, reference = yreal)

```

## Aprendizado não-supervisionado

Existem apenas as variáveis de input $X$, não há um output desejado. O objetivo é modelar a estrutura dos dados para encontrar padrões ou grupos. Isso é chamado de não-supervisionado porque não há uma resposta correta. A interpretação dos resultados é mais complexa e depende do feeling humano. Uma aplicação comum são os algoritmos de clusterização, como **k-means**, que divide um conjunto de dados em **k** grupos ou clusteres, que são homogêneos ou semelhantes entre si.

### Clusterização com k-means

K-means é um método de clusterização simples cujo objetivo é particionar um conjunto de dados em $K$ grupos distintos e excludentes. Para aplicar o k-means, deve-se primeiro especificar o número desejado de clusteres $K$. Com isso, o algoritmo irá atribuir a cada observação (ex.: linha de uma tabela) um dos $K$ clusteres.

A ideia por trás do k-means é que uma boa clusterização reduz a variância dentro dos clusteres, a chamada within-cluster variation, que é definida como uma medida da diferença entre duas diferentes observações. Uma das diversas maneiras de calcular a within-cluster variation é por meio da distância euclidiana.

O gráfico abaixo mostra o resultado da aplicação do k-means para um dataset aleatório de 100 observações e 2 variáveis para diferentes valores de $K$.


```{r}
# a 2-dimensional example
simular_cluster <- function(seed_input, seed_kmeans, k){
  # simular matriz de input
  set.seed(seed_input)
  matriz_input <- matrix(c(rnorm(100, mean = 0, sd = 1),
                           rnorm(100, mean = 1.5, sd = 0.5)),
                           ncol = 2)
  colnames(matriz_input) <- c("x", "y")
  # aplicar k-means
  set.seed(seed_kmeans)
  cl <- kmeans(matriz_input, k)
  plot(x = matriz_input,
       col = cl$cluster, main = paste("k = ", k),
       xlab = "", ylab = "")
}

par(mfrow = c(2, 2))
simular_cluster(123, 1, 2)
simular_cluster(123, 1, 3)
simular_cluster(123, 1, 4)
simular_cluster(123, 1, 5)
```

#### Estudo de caso: clusterizando países por indicadores econômicos

[Link para blog post](http://sillasgonzaga.com/post/clusterizacao-r-paises/)