Projeto Final - Inferência e Análise de Regressão

“Você consegue predizer a eficiência de combustível de um carro?”

Motivação do trabalho

O presente trabalho tenta responder à pergunta “Você consegue predizer a eficiência de combustível de um carro?”, e para isso utiliza da base de dados AutoMPG, disponibilizado pelo centro de estudos em aprendizado de máquina da UC Irvine. O site informa que esta base de dados é uma versão modificada da base de dados original disponibilizada pela StatLib.

Análise Exploratória dos dados

Primeiro, precisamos importar as bibliotecas necessárias:

library(knitr)
library(rmarkdown)
library(htmltools)

library(MASS)
library(tidyverse)
library(GGally)
library(stargazer)
library(car)

theme_set(theme_classic()) # importante para termos o mesmo layout para os gráficos

Vamos importar a base de dados e observar a organização tabular dos mesmos:

car_data <- read.table("./Auto MPG/auto-mpg.data", header=FALSE)

# Vamos adicionar os nomes às variáveis
names(car_data) <- c("mpg", "cylinders", "displacement", "horsepower", 
                     "weight", "acceleration", "model_year", "origin", 
                     "car_name")

glimpse(car_data)

## Rows: 398
## Columns: 9
## $ mpg          <dbl> 18, 15, 18, 16, 17, 15, 14, 14, 14, 15, 15, 14, 15, 14, 2…
## $ cylinders    <int> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 4, 6, 6, 6, 4, …
## $ displacement <dbl> 307, 350, 318, 304, 302, 429, 454, 440, 455, 390, 383, 34…
## $ horsepower   <chr> "130.0", "165.0", "150.0", "150.0", "140.0", "198.0", "22…
## $ weight       <dbl> 3504, 3693, 3436, 3433, 3449, 4341, 4354, 4312, 4425, 385…
## $ acceleration <dbl> 12.0, 11.5, 11.0, 12.0, 10.5, 10.0, 9.0, 8.5, 10.0, 8.5, …
## $ model_year   <int> 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 70, 7…
## $ origin       <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 3, …
## $ car_name     <chr> "chevrolet chevelle malibu", "buick skylark 320", "plymou…

paged_table(car_data)

Aqui podemos observar alguns pontos, como:

• O dataset possui 398 linhas e 9 colunas;
• A variável horsepower está representada como caractere
• As variável car_name à primeira vista não parece ser útil para a análise, já que o nome do carro não interfere no gasto de combustível

Agora também é possível observar que desejamos predizer o valor da variável dependente mpg(miles per galon) atraveś de uma regressão linear múltipla, utilizando como variáveis independentes as demais variáveis do dataset.

Na chunk abaixo nós realizamos a transformação de algumas variáveis e retiramos possíveis NA’s que apareceram no nosso dataset.

car_data <- car_data %>% 
                filter(car_data$horsepower != "?" ) %>% 
                select(c(-car_name))

car_data$horsepower <- as.numeric(car_data$horsepower)
car_data$origin <- as.factor(car_data$origin)
car_data$cylinders <- as.factor(car_data$cylinders)

Agora vamos verificar a correlação entre as variáveis duas a duas, com um gráfico pairplot:

ggpairs(car_data, 
            title = "Análise dois a dois", 
            mapping = aes(color = cylinders),
            legend = 1
        )

Observando a tabela acima, podemos notar que, aparentemente, carros com 4 cilindros tem uma eficiência melhor que carros com mais cilindros, e que as variáveis ‘weight’ e ‘horsepower’ parecem ser positivamente correlacionadas. Entretanto, a variável ‘horsepower’ aparenta estar negativamente correlacionada com a variável ‘acceleration’. Vale observar que a variável ‘mpg’ não parece ter correlação linear com nenhuma das variáveis independentes do dataset.

Ajustando o modelo de regressão

Primeiro Modelo

Para a confecção do primeiro modelo, vamos ajustar uma reta com todas as variáveis do nosso dataset:

model <- lm(mpg ~ .,
            data = car_data)
stargazer::stargazer(model, type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                 mpg            
## -----------------------------------------------
## cylinders4                   6.722***          
##                               (1.654)          
##                                                
## cylinders5                   7.078***          
##                               (2.516)          
##                                                
## cylinders6                    3.351*           
##                               (1.824)          
##                                                
## cylinders8                    5.099**          
##                               (2.109)          
##                                                
## displacement                 0.019***          
##                               (0.007)          
##                                                
## horsepower                   -0.035***         
##                               (0.013)          
##                                                
## weight                       -0.006***         
##                               (0.001)          
##                                                
## acceleration                   0.026           
##                               (0.093)          
##                                                
## model_year                   0.737***          
##                               (0.049)          
##                                                
## origin2                      1.764***          
##                               (0.551)          
##                                                
## origin3                      2.617***          
##                               (0.527)          
##                                                
## Constant                    -22.080***         
##                               (4.541)          
##                                                
## -----------------------------------------------
## Observations                    392            
## R2                             0.847           
## Adjusted R2                    0.842           
## Residual Std. Error      3.098 (df = 380)      
## F Statistic          191.118*** (df = 11; 380) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

anova(model)

## Analysis of Variance Table
## 
## Response: mpg
##               Df  Sum Sq Mean Sq  F value    Pr(>F)    
## cylinders      4 15274.5  3818.6 397.9576 < 2.2e-16 ***
## displacement   1  1098.0  1098.0 114.4233 < 2.2e-16 ***
## horsepower     1   588.0   588.0  61.2785 4.979e-14 ***
## weight         1   715.1   715.1  74.5286 < 2.2e-16 ***
## acceleration   1     7.7     7.7   0.7991    0.3719    
## model_year     1  2245.4  2245.4 234.0069 < 2.2e-16 ***
## origin         2   244.0   122.0  12.7130 4.526e-06 ***
## Residuals    380  3646.3     9.6                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Aqui podemos perceber algumas informações interessantes, como o R² de aproximadamente 0.84. Porém, ao observar o teste anova, observamos que a variável “acceleration” não é significativa, e podemos retirar do modelo.

Segundo Modelo

Sendo assim, criamos um segundo modelo retirando a variável não significativa:

model_2 <- lm(mpg ~ displacement + horsepower + weight + cylinders + origin + model_year,
              data = car_data)
stargazer::stargazer(model_2, type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                                 mpg            
## -----------------------------------------------
## displacement                  0.018**          
##                               (0.007)          
##                                                
## horsepower                   -0.037***         
##                               (0.011)          
##                                                
## weight                       -0.006***         
##                               (0.001)          
##                                                
## cylinders4                   6.784***          
##                               (1.637)          
##                                                
## cylinders5                   7.147***          
##                               (2.501)          
##                                                
## cylinders6                    3.403*           
##                               (1.813)          
##                                                
## cylinders8                    5.137**          
##                               (2.102)          
##                                                
## origin2                      1.763***          
##                               (0.551)          
##                                                
## origin3                      2.621***          
##                               (0.526)          
##                                                
## model_year                   0.736***          
##                               (0.049)          
##                                                
## Constant                    -21.623***         
##                               (4.231)          
##                                                
## -----------------------------------------------
## Observations                    392            
## R2                             0.847           
## Adjusted R2                    0.843           
## Residual Std. Error      3.094 (df = 381)      
## F Statistic          210.731*** (df = 10; 381) 
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

anova(model_2)

## Analysis of Variance Table
## 
## Response: mpg
##               Df  Sum Sq Mean Sq  F value    Pr(>F)    
## displacement   1 15440.2 15440.2 1612.999 < 2.2e-16 ***
## horsepower     1   383.5   383.5   40.062 6.921e-10 ***
## weight         1  1015.3  1015.3  106.067 < 2.2e-16 ***
## cylinders      4   836.6   209.2   21.851 2.994e-16 ***
## origin         2   309.8   154.9   16.182 1.799e-07 ***
## model_year     1  2186.5  2186.5  228.422 < 2.2e-16 ***
## Residuals    381  3647.1     9.6                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

shapiro.test(model_2$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_2$residuals
## W = 0.96713, p-value = 1.035e-07

Agora observamos um R² de 0.84, o que é bom, e a análise anova nos mostra que exitem apenas variáveis significativas no modelo. Utilizaremos a função stepAIC, muito comum para a realização do processo de ‘feature selection’, e averiguar a possibilidade de realizar alguma alteração no modelo.

backward <- stepAIC(model_2, direcion="backward", trace=FALSE)
anova(backward)

## Analysis of Variance Table
## 
## Response: mpg
##               Df  Sum Sq Mean Sq  F value    Pr(>F)    
## displacement   1 15440.2 15440.2 1612.999 < 2.2e-16 ***
## horsepower     1   383.5   383.5   40.062 6.921e-10 ***
## weight         1  1015.3  1015.3  106.067 < 2.2e-16 ***
## cylinders      4   836.6   209.2   21.851 2.994e-16 ***
## origin         2   309.8   154.9   16.182 1.799e-07 ***
## model_year     1  2186.5  2186.5  228.422 < 2.2e-16 ***
## Residuals    381  3647.1     9.6                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Agora utilizamos a função VIF(Variance Inflaction Factor), que é um preditor que nos auxilia a verificar a multicolinearidade no modelo(saiba mais). De forma geral, podemos afirmar que um VIF maior que 5 ou 10 é ruim, e o modelo apresenta problemas ao estimar os valores desejados.

vif(backward)

##                   GVIF Df GVIF^(1/(2*Df))
## displacement 22.984822  1        4.794249
## horsepower    6.947667  1        2.635843
## weight        9.028147  1        3.004687
## cylinders    15.142091  4        1.404505
## origin        2.342253  2        1.237110
## model_year    1.313734  1        1.146182

Observe que as variáveis ‘displacement’ e ‘cilynders’ tem valores bem acima de 10, então começaremos retirando a variável ‘displacement’.

Terceiro Modelo

model_3 <- lm(mpg ~ horsepower + weight + cylinders + origin + model_year,
              data = car_data)
anova(model_3)

## Analysis of Variance Table
## 
## Response: mpg
##             Df  Sum Sq Mean Sq  F value    Pr(>F)    
## horsepower   1 14433.1 14433.1 1485.829 < 2.2e-16 ***
## weight       1  2392.1  2392.1  246.254 < 2.2e-16 ***
## cylinders    4   850.5   212.6   21.888 2.795e-16 ***
## origin       2   301.3   150.6   15.508 3.344e-07 ***
## model_year   1  2131.4  2131.4  219.421 < 2.2e-16 ***
## Residuals  382  3710.7     9.7                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

vif(model_3)

##                GVIF Df GVIF^(1/(2*Df))
## horsepower 5.713676  1        2.390330
## weight     7.449175  1        2.729318
## cylinders  8.535167  4        1.307379
## origin     2.015556  2        1.191513
## model_year 1.298248  1        1.139407

Ainda encontramos um vif alto para a variávell ‘cylinders’, então também iremos retirá-la:

model_4 <- lm(mpg ~ horsepower + weight + origin + model_year,
              data = car_data)
anova(model_4)

## Analysis of Variance Table
## 
## Response: mpg
##             Df  Sum Sq Mean Sq  F value    Pr(>F)    
## horsepower   1 14433.1 14433.1 1295.250 < 2.2e-16 ***
## weight       1  2392.1  2392.1  214.669 < 2.2e-16 ***
## origin       2   307.4   153.7   13.792 1.639e-06 ***
## model_year   1  2385.3  2385.3  214.057 < 2.2e-16 ***
## Residuals  386  4301.2    11.1                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

vif(model_4)

##                GVIF Df GVIF^(1/(2*Df))
## horsepower 4.537936  1        2.130243
## weight     4.919323  1        2.217955
## origin     1.668152  2        1.136472
## model_year 1.266033  1        1.125181

Agora encontramos valores abaixo de 5 em todas as variáveis independentes, indicando que chegamos possivelmente a um modelo eficiente.

Avaliação dos pressupostos da regressão

Observe que podemos tentar linearizar mais a reta, ou seja, eliminar a não linearidade entre as variáveis; para isso, utilizaremos a transformação de BoxCox:

library(fpp)
lambda <- BoxCox.lambda(car_data$mpg, method=c("loglik"), lower=-3, upper= 3)
data_t <- BoxCox(car_data$mpg, lambda)

model_5 <- lm(data_t ~ horsepower + weight + origin + model_year,
              data = car_data)

stargazer::stargazer(model_5, type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               data_t           
## -----------------------------------------------
## horsepower                   -0.002**          
##                               (0.001)          
##                                                
## weight                      -0.0005***         
##                              (0.00003)         
##                                                
## origin2                      0.133***          
##                               (0.035)          
##                                                
## origin3                      0.134***          
##                               (0.036)          
##                                                
## model_year                   0.057***          
##                               (0.004)          
##                                                
## Constant                     1.464***          
##                               (0.287)          
##                                                
## -----------------------------------------------
## Observations                    392            
## R2                             0.871           
## Adjusted R2                    0.869           
## Residual Std. Error      0.228 (df = 386)      
## F Statistic          520.155*** (df = 5; 386)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

anova(model_5)

## Analysis of Variance Table
## 
## Response: data_t
##             Df  Sum Sq Mean Sq  F value    Pr(>F)    
## horsepower   1 104.966 104.966 2016.612 < 2.2e-16 ***
## weight       1  15.576  15.576  299.247 < 2.2e-16 ***
## origin       2   1.169   0.584   11.226 1.824e-05 ***
## model_year   1  13.661  13.661  262.461 < 2.2e-16 ***
## Residuals  386  20.091   0.052                       
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

vif(model_5)

##                GVIF Df GVIF^(1/(2*Df))
## horsepower 4.537936  1        2.130243
## weight     4.919323  1        2.217955
## origin     1.668152  2        1.136472
## model_year 1.266033  1        1.125181

E iremos avaliar novamente as pressuposições da nossa regressão:

Linearidade

Observamos que apesar de não ser o ideal, a curva se tornou mais suave nas pontas, possivelmente melhorando o teste da normalidade dos resíduos.

Normalidade dos Resíduos

shapiro.test(model_5$residuals)

## 
##  Shapiro-Wilk normality test
## 
## data:  model_5$residuals
## W = 0.9923, p-value = 0.04056

ggplot() +
    ggtitle("Gráfico QQPlot") +
    geom_qq(aes(sample = rstandard(model_5))) +
    geom_abline(color = "red") +
    coord_fixed()

Observe que obtivemos um p-valor < 0.5, indicando que há evidências estatísticas suficientes para rejeitar a hipótese nula de que os resíduos da regressão seguem uma distribuição normal. Em outras palavras, o p-valor é menor que um nível de significância comum (0.05, por exemplo), o que sugere que a distribuição dos resíduos não é normal.

Homocedasticidade

# Teste de Breusch-Pagan
bptest(model_5)

## 
##  studentized Breusch-Pagan test
## 
## data:  model_5
## BP = 13.52, df = 5, p-value = 0.01896

res_sqrt <- sqrt(abs(rstandard(model_5)))

ggplot(car_data, aes(fitted(model_5), res_sqrt)) +
    ggtitle("Verificando a Homocedasticidade") +
    geom_point() +
    geom_smooth()

## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'

Multicolinearidade

vif(model_5)

##                GVIF Df GVIF^(1/(2*Df))
## horsepower 4.537936  1        2.130243
## weight     4.919323  1        2.217955
## origin     1.668152  2        1.136472
## model_year 1.266033  1        1.125181

Conclusão

stargazer::stargazer(model_5, type = "text")

## 
## ===============================================
##                         Dependent variable:    
##                     ---------------------------
##                               data_t           
## -----------------------------------------------
## horsepower                   -0.002**          
##                               (0.001)          
##                                                
## weight                      -0.0005***         
##                              (0.00003)         
##                                                
## origin2                      0.133***          
##                               (0.035)          
##                                                
## origin3                      0.134***          
##                               (0.036)          
##                                                
## model_year                   0.057***          
##                               (0.004)          
##                                                
## Constant                     1.464***          
##                               (0.287)          
##                                                
## -----------------------------------------------
## Observations                    392            
## R2                             0.871           
## Adjusted R2                    0.869           
## Residual Std. Error      0.228 (df = 386)      
## F Statistic          520.155*** (df = 5; 386)  
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Ao final do ajuste do modelo, foi possível observar um valor R2 de 0.87, indicando que o modelo tem uma boa capacidade explicativa. As verificações das hipóteses também foram melhores após a realização da transformação de BoxCox, apesar de não serem suficientes(Não passou nos testes de normalidade e homocedasticidade).

Uma das possíveis alternativas seriam a criação de modelos não-lineares, como modelos logísticos, já que o modelo linear não foi suficiente para explicar de forma desejável o problema proposto.

Assim, o modelo final apresenta o seguinte formato:

\[\begin{equation} Y = \beta_0 + \beta_1 horsepower + \beta_2 weight + \beta_3 origin + \beta_4 modelyear + \varepsilon \end{equation}\]

Onde ‘origin’ é uma variável dummy.