library(openintro)
library(tidyverse)
library(ggbeeswarm)
library(modelr)
library(broom)
theme_set(theme_bw())

1 Introdução

O exercício é sobre regressão com múltiplas variáveis, nele temos um conjunto de dados sobre professores e desejamos encontrar o modelo linear que melhor represente a nossa variável resposta em função de todas as outras.

2 Os dados

Os dados que usaremos estão detalhados nessa página do OpenIntro.

profs = read_csv("https://www.openintro.org/stat/data/evals.csv") %>% 
    select(score, age, gender, bty_avg, pic_outfit, pic_color)
## Parsed with column specification:
## cols(
##   .default = col_double(),
##   rank = col_character(),
##   ethnicity = col_character(),
##   gender = col_character(),
##   language = col_character(),
##   cls_level = col_character(),
##   cls_profs = col_character(),
##   cls_credits = col_character(),
##   pic_outfit = col_character(),
##   pic_color = col_character()
## )
## See spec(...) for full column specifications.
glimpse(profs)
## Rows: 463
## Columns: 6
## $ score      <dbl> 4.7, 4.1, 3.9, 4.8, 4.6, 4.3, 2.8, 4.1, 3.4, 4.5, 3.8, 4.5…
## $ age        <dbl> 36, 36, 36, 36, 59, 59, 59, 51, 51, 40, 40, 40, 40, 40, 40…
## $ gender     <chr> "female", "female", "female", "female", "male", "male", "m…
## $ bty_avg    <dbl> 5.000, 5.000, 5.000, 5.000, 3.000, 3.000, 3.000, 3.333, 3.…
## $ pic_outfit <chr> "not formal", "not formal", "not formal", "not formal", "n…
## $ pic_color  <chr> "color", "color", "color", "color", "color", "color", "col…

Acima podemos ter uma visão geral dos nossos dados observando as variáveis disponíveis para análise. As nossas variáveis são score que representa a pontuação de um professor, age que é a idade dos professores, gender que é o gênero do professor, bty_avg que é uma pontuação quanto a beleza do professor, pic_outfit que indica se a foto do professor é formal ou não e pic_color que indica se a foto é colorida ou não.

Abaixo poderemos ver algumas visualizações gráficas que relacionam o score com os outros atributos do professor e com isso poderemos ter uma idea inicial das relações entre os atributos.

profs %>% 
    ggplot(aes(x = age, y = score)) + 
    geom_count()

profs %>% 
    ggplot(aes(x = bty_avg, y = score)) + 
    geom_count()

profs %>% 
    ggplot(aes(x = pic_outfit, y = score)) + 
    geom_count()

profs %>% 
    ggplot(aes(x = pic_color, y = score)) + 
    geom_count()

Os gráficos que comparam score com a idade e com o nível de beleza respectivamente, não deixam muito claro como o score varia em função de ambas as variáveis, mas quanto ao gráfico relativo a idade, podemos destacar que a concentração de professores é maior abaixo de 60 anos. Os gráficos relativos a formalidade da foto e a cor da foto respectivamente tem resultados bem semelhantes, no gráfico que analisa a formalidade a amplitude dos scores quando a foto é informal é maior, assim como a concentração de fotos informais é maior, enquanto que no gráfico relativo a cor da foto a amplitude dos scores quando a foto é colorida é maior e a concentração de fotos coloridas também é maior.

Abaixo vamos visualizar o gráfico de densidade da variável score que na nossa regressão será nossa variável de respostas. Com a análise do gráfico abaixo poderemos ter uma ideia de por volta de quais valores se concentram os scores dos professores.

profs %>%
  ggplot(aes(x = score)) + 
  geom_density(fill = "darkred") + 
  labs(title = "Concentração dos scores", x = "Score", y = "Densidade")

Podemos ver no gráfico acima que os scores dos professores se concentram mais em torno de 4.5, que é uma nota bem alta, mais adiante vamos ver que fatores influênciam no resultado do score de um determinado professor.

Abaixo podemos ver um sumário estatístico com a média, mediana e desvio padrão da variável na qual estamos interessados: score.

profs %>%
  summarise(média_scores = mean(score), mediana_score = median(score),
            desvio_padrao = sd(score))
## # A tibble: 1 x 3
##   média_scores mediana_score desvio_padrao
##          <dbl>         <dbl>         <dbl>
## 1         4.17           4.3         0.544

Com base nos resultados acima, podemos dizer que não há muitos dados destoantes do comportamento geral da variável, a média e a mediana são muito próximas e o desvio padrão é baixo, ou seja, há pouca chance de haver pontos extremos. Além disso, podemos reforçar o que foi observado no gráfico de densidade, os scores dos professores realmente tendem a ser valores altos, se concentrando bastante acima de 4.

3 Calculando o modelo

Queremos descobrir o efeito da variação das variáveis de preditoras (age, gender, bty_avg, pic_color, pic_outfit) sobre a variável resposta que é score.

modelo = lm(score ~ age + gender + bty_avg + pic_color + pic_outfit, profs)
tidy(modelo, conf.int = TRUE) %>% select(term, estimate, conf.low, conf.high)
## # A tibble: 6 x 4
##   term                 estimate conf.low conf.high
##   <chr>                   <dbl>    <dbl>     <dbl>
## 1 (Intercept)           4.37      3.96     4.78   
## 2 age                  -0.00712  -0.0126  -0.00167
## 3 gendermale            0.223     0.121    0.326  
## 4 bty_avg               0.0473    0.0126   0.0821 
## 5 pic_colorcolor       -0.212    -0.348   -0.0765 
## 6 pic_outfitnot formal -0.0152   -0.148    0.118

Com os resultados calculado acima podemos ver graficamente os intervalos onde temos confiança que os coeficientes calculados para cada uma das variáveis preditivas está.

modelo %>% 
  tidy(conf.int = TRUE) %>%
  filter(term != "(Intercept)") %>% 
  ggplot(aes(estimate, reorder(term, estimate), 
             xmin = conf.low, 
             xmax = conf.high, 
             height = 0.1)) +
  geom_vline(xintercept = 0, color = "grey") +
  # geom_point() +
  geom_errorbarh()

Com base base no resultado podemos montar o seguinte modelo: \(score = 4,4 + 0,22.gender + 0,05.btyavg - 0,007.age - 0,21.piccolor - 0,015.picoutfit\). Com base nesse resultado podemos tirar alguma conclusões sobre as variáveis preditivas:

glance(modelo)
## # A tibble: 1 x 11
##   r.squared adj.r.squared sigma statistic p.value    df logLik   AIC   BIC
##       <dbl>         <dbl> <dbl>     <dbl>   <dbl> <int>  <dbl> <dbl> <dbl>
## 1    0.0878        0.0778 0.522      8.79 5.60e-8     6  -353.  720.  749.
## # … with 2 more variables: deviance <dbl>, df.residual <int>

Como podemos ver no resultado do \(R^2\) o nosso modelo representa muito pouco da relação, sendo \(R^2 = 0,087\), ou seja, representando cerca de 8,7% da relação que analisamos.