CTRL+ALT+N

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```



## Pacotes

```{r, message=FALSE, warning=FALSE}
library(broom)
library(readxl)
library(janitor)
library(MASS)
library(lime)
library(GGally)
library(car)
library(ISLR)
library(tidymodels)
library(tidyverse)
```



## Os dados

```{r}
# cars
glimpse(cars)
grafico_muito_pesado <- ggplot(cars) + 
  geom_point(aes(x = speed, y = dist))

# saveRDS(grafico_muito_pesado, file = "grafico_muito_pesado.rds")
grafico_muito_pesado <- readRDS("grafico_muito_pesado.rds")
grafico_muito_pesado
```

## ajuste de uma regressão linear simples

```{r}
melhor_reta <- lm(dist ~ speed , data = cars)
melhor_reta
```

## tabela com as predições junto

```{r}
broom::confint_tidy(melhor_reta)
cars_com_predicoes <- augment(melhor_reta)
glimpse(cars_com_predicoes)
cars_com_predicoes
```

## Gráfico com a reta ajustada
```{r}
ggplot(cars_com_predicoes) + 
  geom_point(aes(x = speed, y = dist)) +
  geom_line(aes(x = speed, y = .fitted))
```



## Exercício 1

Calcule o EQM da melhor reta e compare com a saída do `summary(melhor_reta)`.

```{r}
# dicas:
# yh é o y chapéu (ou valor ajustado). O nome da coluna retornado pelo augment() para o yh é o .fitted
# EQM = mean((y - yh)^2).
y <- cars_com_predicoes$dist 
yh <- cars_com_predicoes$.fitted 

EQM <- sqrt(mean((y - yh)^2))

mean((y - yh)^2) # m^2
sqrt(mean((y - yh)^2)) # m
```

```{r}
summary(melhor_reta) # compare o resultado obtido com essa saída
```



## Exercício 2 

calcule beta0 e beta1

```{r}
x <- cars$speed
y <- cars$dist

xbarra <- mean(x) #média de X
ybarra <- mean(y) #média de y

a <- 0
for(i in 1:length(x)) {
  a <- a + (x[i] - xbarra)*(y[i] - ybarra)
}

b <- 0
for(i in 1:length(x)) {
  b <- b + (x[i] - xbarra)^2
}

a <- sum((x - xbarra)*(y - ybarra))
b <- sum((x - xbarra)^2)
  
beta1 <- a/b
beta0 <- ybarra - beta1 * xbarra
```





## Exercício 3

Tire as informações do objeto `melhor_reta` para decidir se speed está associado com dist. Use a função `summary()`.

```{r}
summary(melhor_reta)
```

## Exercício 4 

Interprete o parâmetro `speed` (beta 1).

```{r}

```


## Exercício 5
Calcule o R2 para a `melhor reta` e depois compare com o valor da saída do `summary()`.

```{r}
sqr <- sum((cars_com_predicoes$dist - cars_com_predicoes$.fitted)^2)
sqt <- sum((cars_com_predicoes$dist - mean(cars_com_predicoes$dist))^2)

r2 <- 1 - sqr/sqt
  
summary(melhor_reta)
```


## Exercício 6

Calcule o R2 ajustado para a `melhor reta` e depois compare com o valor da saída do `summary()`.

```{r}
N <- nrow(cars)
r2aj <-  1 - sqr * (N - 1)/(sqt * (N - 2))
```


## Exercício 7

Estude os gráficos que saem do `plot(melhor_reta)`. Procure por outliers.

```{r}
plot(melhor_reta)
```


## Exercício 8

Dados para o exercício: Credit.

- Balance (variável resposta): Saldo no Cartão de Crédito
- Gender (variável explicativa): Masculino ou Feminino

```{r}
glimpse(Credit)
```

Interprete os parametros do modelo abaixo

```{r}
modelo_balance <- lm(Balance ~ Gender, data = Credit)
summary(modelo_balance)
```



## Exercício 9

Calcule as médias de `Balance` para cada `Gender`. Compare os resultados com o exercício anterior.

```{r}
media_f <- mean(Credit$Balance[Credit$Gender == "Female"])
media_m <- mean(Credit$Balance[Credit$Gender == " Male"])

media_m - media_f
```


## Exercício 10

Explore como usar a função model.matrix(). Use a função `model.matrix()` no lugar de `lm()` do exercício 8 e veja sua saída.

```{r}
model.matrix(Balance ~ Gender, data = Credit) %>% head()
Credit %>% select(Balance, Gender)

y = b0*1 + b1*0 = b0



y = b0*1 + b1*1 = 509 + 19
```


## Exercício 11

Repita os exercícios 8, 9 e 10, mas agora usando Ethnicity no lugar de Gender.

```{r}
modelo_balance <- lm(Balance ~ Ethnicity, data = Credit)
summary(modelo_balance)
```

```{r}
Credit %>%
  group_by(Ethnicity) %>%
  summarise(
    media = mean(Balance)
  )
531.0000	+ (0)
531.0000	+ (0)
```
```{r}

model.matrix(Balance ~ Ethnicity, data = Credit) %>% head()
```


## Exercício 12

Crie um boxplot de Balance para cada Ethnicity e avalie se a análise visual é compatível com o que os valores-p indicam.

```{r}
Credit %>%
  ggplot() +
  geom_boxplot(aes(x = Balance, y = Ethnicity, fill = Ethnicity))
```


## Exercício 13: Use os dados simulados `y_x` abaixo:

**dados para o exercício**

```{r}
set.seed(1)
y_x <- tibble(
  x = runif(60),
  y = 10 + 0.5*log(x) + rnorm(30, sd = 0.1)
) %>%
  mutate(
    x2 = log(x)
  )

ggplot(y_x) + geom_point(aes(x = x2, y = y))
```

Use `lm()` para ajustar os seguintes modelos:

1) y ~ x, data = y_x
2) y ~ log10(x), data = y_x

Avalie qual modelo é melhor quanto ao EQM e quanto ao R^2.

```{r}
y_x$log_x <- log(y_x$x)

modelo_lin <- lm(y ~ x, data = y_x)
modelo_log <- lm(y ~ BoxCox(x), data = y_x)

summary(modelo_log)
summary(modelo_lin)

plot(modelo_lin)
```


## Exercício 14 

**dados para o exercício:**

```{r}
set.seed(1)
y_x_poly <- tibble(
  x = runif(30, 0, 20),
  y = 500 + 0.4 * (x-10)^3 + rnorm(30, sd = 50)
)
```

Explore a saída de model.matrix() e lm() utlizando as fórmulas:

a) `y ~ x`
b) `y ~ x + I(x^2) + I(x^3)`
c) `y ~ poly(x, 2, raw = TRUE)`
d) `y ~ poly(x, 3, raw = TRUE)`

```{r}
# Resolução do item a)
a_mm <- model.matrix(y ~ x, data = y_x_poly)
a_lm <- lm(y ~ x, data = y_x_poly)
summary(a_lm)
```

```{r}
# Item b)
b_mm <-  model.matrix(y ~ poly(x, 2, raw = TRUE), data = y_x_poly)
b_lm <-  lm(y ~ poly(x, 2, raw = TRUE), data = y_x_poly)
summary(b_lm)
head(b_mm)
```

```{r}
# Item c)
c_mm <- model.matrix(y ~ x + I(x^2) + I(x^3) , data = y_x_poly)
c_lm <- lm(y ~ x + I(x^2) + I(x^3), data = y_x_poly)
summary(c_lm)
```

```{r}
# Item d)
d_mm <- model.matrix(y ~ poly(x, 3, raw = TRUE), data = y_x_poly)
d_lm <- lm(y ~ poly(x, 3, raw = TRUE), data = y_x_poly)
summary(d_lm)

head(d_mm)
```

Compare o número de parâmetros com o número de colunas. Qual conclusão?









