lalala.Rmd

---
title: "Regressão Linear, agora com lm"
author: "Nazareno Andrade"
output: 
  html_notebook:
    theme: readable
    fig_width: 7
    toc: true
    toc_float: true
---

```{r setup, message=FALSE, warning=FALSE}
library(openintro)
library(tidyverse)
library(tidymodels)
theme_set(theme_bw())
```

```{r}
data("county_complete")

counties = county_complete %>% 
  select(hs_grad = hs_grad_2017, 
         unemployment = unemployment_rate_2017,
         poverty = poverty_2017) %>% 
  filter(complete.cases(.))

glimpse(counties)
```

```{r}
dados = read_csv(here::here("data/participation-per-country.csv"), show_col_types = FALSE)

dados = dados %>% 
  drop_na()

glimpse(dados)
```

```{r}
dados %>% 
  group_by(site) %>% 
  ggplot(aes(x = EPI, fill = site, color = site, alpha = 0.3)) +
  geom_density()

dados %>% 
  group_by(site) %>% 
  ggplot(aes(x = responderam_prop, fill = site, color = site, alpha = 0.5)) +
  geom_density()

dados %>% 
  group_by(site) %>% 
  ggplot(aes(x = GNI, fill = site, color = site, alpha = 0.5)) +
  geom_density()
```


# Agora automaticamente

lm == linear model

```{r}
dados_so = dados %>% 
  filter(site == "StackOverflow")

dados_so %>% 
  ggplot(aes(y = responderam_prop, x = EPI)) + 
  geom_point(alpha = 0.4, size = 1) + 
  # geom_smooth(method = "lm", se = FALSE)
    NULL
```

```{r}
mod1 <- lm(responderam_prop ~ EPI, # poverty = b0 + b1*hs_grad
           data = dados_so)

# sintaxe base R, que não usaremos
# summary(mod)

# sintaxe broom / tidymodels: 
tidy(mod1)
glance(mod1) 
```

responderam = 0.086 + 0.0077 * EPI

Função retornada foi $poverty = 71.3 - 0.642 * hs\_grad$

## Estimativas/previsões do modelo

```{r}
mod1 %>% 
  augment(dados_so)
```

```{r}
mod1 %>% 
  augment(dados_so) %>% 
  ggplot(mapping = aes(x = EPI)) + 
  geom_point(aes(y = responderam_prop), alpha = 0.4, size = 1) + 
  geom_line(aes(y = .fitted), colour = "red")  + 
  NULL
```

## Inferência

```{r}
mod1 %>% 
    tidy(conf.int = T, conf.level = .95) %>% 
    select(-p.value)
```

## Resíduos

```{r}
mod1 %>% 
  augment(dados_so) %>% 
  ggplot(aes(EPI, .resid)) + 
  geom_point(alpha = .4, size = 1) + 
  geom_hline(yintercept = 0, colour = "blue")
```

```{r}
dados_su = dados %>% 
  filter(site == "SuperUser")

dados_su %>% 
  ggplot(aes(y = responderam_prop, x = EPI)) + 
  geom_point(alpha = 0.4, size = 1) + 
  # geom_smooth(method = "lm", se = FALSE)
    NULL
```

```{r}
mod2 <- lm(responderam_prop ~ EPI, # poverty = b0 + b1*hs_grad
           data = dados_su)

# sintaxe base R, que não usaremos
# summary(mod)

# sintaxe broom / tidymodels: 
tidy(mod2)
glance(mod2) 
```

responderam = 0.094 + 0.0027 * EPI

```{r}
mod2 %>% 
  augment(dados_su)
```
```{r}
mod2 %>% 
  augment(dados_su) %>% 
  ggplot(mapping = aes(x = EPI)) + 
  geom_point(aes(y = responderam_prop), alpha = 0.4, size = 1) + 
  geom_line(aes(y = .fitted), colour = "red")  + 
  NULL
```

## Inferência

```{r}
mod2 %>% 
    tidy(conf.int = T, conf.level = .95) %>% 
    select(-p.value)
```

```{r}
dados_so %>% 
  ggplot(aes(y = responderam_prop, x = log(GNI))) + 
  geom_point(alpha = 0.4, size = 1) + 
  # geom_smooth(method = "lm", se = FALSE)
    NULL
```

```{r}
mod3 <- lm(responderam_prop ~ log(GNI), # poverty = b0 + b1*hs_grad
           data = dados_so)

# sintaxe base R, que não usaremos
# summary(mod)

# sintaxe broom / tidymodels: 
tidy(mod3)
glance(mod3) 
```

responderam = -0.023 + 0.028 * GNI

```{r}
mod3 %>% 
  augment(dados_so)
```
```{r}
mod3 %>% 
  augment(dados_so) %>% 
  ggplot(mapping = aes(x = log(GNI))) + 
  geom_point(aes(y = responderam_prop), alpha = 0.4, size = 1) + 
  geom_line(aes(y = .fitted), colour = "red")  + 
  NULL
```

## Inferência

```{r}
mod3 %>% 
    tidy(conf.int = T, conf.level = .95) %>% 
    select(-p.value)
```

------------------------------------------------------------------------

## Exemplo de relato

A conclusão a eque chegamos com um modelo assim é algo como:

> Regressão linear simples foi utilizada para analisar a associação entre escolaridade (medida como % de pessoas que completou o ensino médio) e pobreza (% de pessoas abaixo da linha da pobreza) nos condados. Um modelo no formato pobreza = 73 -0.642\*escolaridade explica 40% da variância da variável de resposta (R2 = 0,40). O aumento de 1% da escolaridade em um condado produz uma decréscimo de 0,642% (IC 95% [0,61; 0,67]) no valor esperado da pobreza. Portanto...


# Atividade de regressao multipla

```{r}
dados %>% 
  group_by(site) %>% 
  ggplot(aes(x = IDV, fill = site, color = site, alpha = 0.3)) +
  geom_density()
```

a. Construa um modelo 1 com a variável responderam_prop com variável de resposta e fluência em inglês da população (EPI), produto interno bruto do país (GNI) e disponibilidade de internet no país como variáveis de explicação. Comente esse modelo em termos dos coeficientes e do ajuste. Estamos interessados em fazer inferência sobre os coeficiente.

```{r}

dados_so = dados %>% 
  filter(site == "StackOverflow")

dados_so %>% 
  ggplot(aes(y = responderam_prop, x = EPI)) + 
  geom_count(alpha = 0.5, size = 2)

dados_so %>% 
  ggplot(aes(y = responderam_prop, x = GNI)) + 
  geom_count(alpha = 0.5, size = 2)

dados_so %>% 
  ggplot(aes(y = responderam_prop, x = Internet)) + 
  geom_count(alpha = 0.5, size = 2)

```

```{r}
mod1 = lm(responderam_prop ~ EPI + GNI + Internet, # poverty = b0 + b1*hs_grad
           data = dados_so)

# sintaxe base R, que não usaremos
# summary(mod)

# sintaxe broom / tidymodels: 
tidy(mod1, conf.int = TRUE)
glance(mod1) 
```

responderam_prop = 0.1605 + 0.0012 * EPI + 0.000012 * GNI - 0.0002 * Internet

```{r}
tidy(mod1, conf.int = T) %>% 
  filter(term != "(Intercept)") %>% 
  ggplot(aes(x = reorder(term, estimate) , y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_linerange() +
  geom_point() +
  xlab("term") +
  coord_flip()
```


Construa um modelo 2 que além das variáveis do modelo 1 tem também o IDV. Esse é um modelo que considera uma variável de cultura. Comparando o modelo 2 com o modelo 1, o que podemos afirmar sobre o efeito do individualismo no comportamento das pessoas de diferentes países no stackoverflow? Há um efeito relevante (lembre de considerar a inferência para a população de onde vem os dados)? O modelo é mais explicativo do que sem a variável relacionada a cultura? 

```{r}

dados_so = dados %>% 
  filter(site == "StackOverflow")

dados_so %>% 
  ggplot(aes(y = responderam_prop, x = IDV)) + 
  geom_count(alpha = 0.5, size = 2) %>% 
  intersept
  
```

```{r}
mod2 <- lm(responderam_prop ~ EPI + GNI + Internet + IDV, # poverty = b0 + b1*hs_grad
           data = dados_su)

# sintaxe base R, que não usaremos
# summary(mod)

# sintaxe broom / tidymodels: 
tidy(mod2, conf.int = TRUE)
glance(mod2) 
```

responderam_prop = 0.193 + 0.00034 * EPI + 0.000001 * GNI - 0.0004 * Internet + 0.0008 * IDV

```{r}
tidy(mod2, conf.int = T) %>% 
  filter(term != "(Intercept)") %>% 
  ggplot(aes(x = reorder(term, estimate) , y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_linerange() +
  geom_point() +
  xlab("term") +
  coord_flip()
```



c. Construa uma outra versão do modelo 2 usando agora os dados do SuperUser. Os resultados são consistentes com os do StackOverflow? Comente e mostre evidência que embase sua conclusão. 

```{r}
mod3 <- lm(responderam_prop ~ EPI + GNI + Internet + IDV, # poverty = b0 + b1*hs_grad
           data = dados_su)

# sintaxe base R, que não usaremos
# summary(mod)

# sintaxe broom / tidymodels: 
tidy(mod3, conf.int = TRUE)
glance(mod3) 
```

responderam_prop = 0.193 + 0.00034 * EPI + 0.000001 * GNI - 0.0004 * Internet + 0.0008 * IDV

```{r}
tidy(mod3, conf.int = T) %>% 
  filter(term != "(Intercept)") %>% 
  ggplot(aes(x = reorder(term, estimate) , y = estimate, ymin = conf.low, ymax = conf.high)) +
  geom_linerange() +
  geom_point() +
  xlab("term") +
  coord_flip()
```