Um analista da indústria automobilística pediu sua ajuda na modelagem de dados sobre os preços dos carros novos. O interesse é em modelar preço de varejo sugerido (suggested retail price) em função de outras covariaveis para 234 carros novos. O conjunto de dados está disponível no site abaixo, no arquivo cars04.csv:
http://www.amstat.org/publications/jse/datasets/04cars.txt.

  1. Analisar (fazer analise descriptiva, ajuste, diagnostico….) o modelo: \[SuggestedRetailPrice = \beta_0 + \beta_1DealerCost + \epsilon\]

Com base na saída deste modelo, o analista concluiu o seguinte: “Uma vez que o modelo explica apenas mais de 99,8% da variabilidade no Preço de varejo sugerido ( Suggested Retail Price ) e o coefciente do custo do revendedor ( Dealer Cost ) tem um valor de t superior a 412, o modelo em questão é um modelo altamente efetivo para produzir intervalos de previsão para o preço de varejo sugerido.”

Forneça uma crítica detalhada desta conclusão.

# Carregando o dataset
carros <- read.csv("cars04.csv", sep = ",")

# Ajustando o modelo
mod1 <- lm(SuggestedRetailPrice~DealerCost ,data = carros)
summary(mod1)
## 
## Call:
## lm(formula = SuggestedRetailPrice ~ DealerCost, data = carros)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1743.52  -262.59    74.92   265.98  2912.72 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -61.904248  81.801381  -0.757     0.45    
## DealerCost    1.088841   0.002638 412.768   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 587 on 232 degrees of freedom
## Multiple R-squared:  0.9986, Adjusted R-squared:  0.9986 
## F-statistic: 1.704e+05 on 1 and 232 DF,  p-value: < 2.2e-16

A conclusão do analista foi equivocada. É fácil ver que os resíduos não são normais já que sua mediana está muito diferente de zero e também que ela é assimétrica, ferindo fortemente a normalide dos resíduos. Para confirmar observe o histograma dos resíduos e o resultado do teste de Shapiro-Wilk para normalidade:

qplot(mod1$residuals, xlab = "Resíduos", bins = 30)

shapiro.test(mod1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod1$residuals
## W = 0.93723, p-value = 1.839e-08

Note também que o desvio padrão dos resíduos e de \(\hat\beta_0\) são muito altas, ocasinando em valores ajustados ruins. Para melhor compreenssão vamos ver o gráfico dos erros absolutos: \[e_{absoluto} = |Y - \hat Y|\]

erro.mod1 <- abs(carros$SuggestedRetailPrice-mod1$fitted.values)
erro.mod1 %>% qplot(xlab = "Erros absolutos", bins = 30)

Ainda podemos confirmar nossas suspeitas observando os gráficos de diagnóstico.

par(mfrow = c(2,2))
plot(mod1)

par(mfrow = c(1,1))

Observe que o gráfico de resíduos vs valores ajustados sugere heterocedasticidade, o qqplot confirma ainda mais nossa hipotíse anterior (não normalidade dos resíduos) e no gráfico de resíduos studentizados vs valores ajustados temos erros correlacionados.

Sendo assim, mesmo que a variável explicativa seja altamente significativa e que o R² seja alto não quer dizer que o modelo ajustado com somente estas características é o melhor. Temos que observar outros fatores.

  1. O analista ficou tão impressionado com suas respostas no item (1), que solicitou a você fazer o próximo estágio da análise de dados, a saber, uma análise dos efeitos de diferentes aspectos de um carro em seu preço de varejo sugerido. Os dados estão disponíveis para todos os 234 carros nas seguintes variáveis:

\(Y\) = SuggestedRetailPrice; \(X_1\) = EngineSize; \(X_2\) = Cylinders; \(X_3\) = Horsepower; \(X_4\) = HighwayMPG; \(X_5\) = Weight; \(X_6\) = WheelBase; \(X_7\) = Hybrid

onde

\(X_7 = 1\) se o carro é hibrido.

\(X_7 = 0\), caso contrário

Faça uma análise para esse modelo, logo responda:

  1. Decida se o modelo é válido. Dê fundamentos para apoiar sua resposta.

Estimando o modelo sugerido:

carros <- carros %>% transform(Hybrid = as.factor(Hybrid))
mod2 <- lm(SuggestedRetailPrice~EngineSize + Cylinders + Horsepower +
             HighwayMPG + Weight + WheelBase + Hybrid , data = carros)
summary(mod2)
## 
## Call:
## lm(formula = SuggestedRetailPrice ~ EngineSize + Cylinders + 
##     Horsepower + HighwayMPG + Weight + WheelBase + Hybrid, data = carros)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -17436  -4134    173   3561  46392 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -68965.793  16180.381  -4.262 2.97e-05 ***
## EngineSize   -6957.457   1600.137  -4.348 2.08e-05 ***
## Cylinders     3564.755    969.633   3.676 0.000296 ***
## Horsepower     179.702     16.411  10.950  < 2e-16 ***
## HighwayMPG     637.939    202.724   3.147 0.001873 ** 
## Weight          11.911      2.658   4.481 1.18e-05 ***
## WheelBase       47.607    178.070   0.267 0.789444    
## Hybrid1        431.759   6092.087   0.071 0.943562    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7533 on 226 degrees of freedom
## Multiple R-squared:  0.7819, Adjusted R-squared:  0.7751 
## F-statistic: 115.7 on 7 and 226 DF,  p-value: < 2.2e-16

Analisando o output acima vemos que o modelo possui algumas variáveis significativas e um R² razoável, contudo, temos o mesmo problema do item anterior, sendo que este é bem pior. Para elucidar nossas dúvidas vamos fazer uma análise de diagnóstico desse modelo.

# Normalidade dos resíduos
qplot(mod2$residuals, xlab = "Resíduos", bins = 30)

shapiro.test(mod2$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod2$residuals
## W = 0.87602, p-value = 6.997e-13
# Gráficos de diagnóstico
par(mfrow = c(2,3))
plot(mod2, which = c(1:5))
par(mfrow = c(1,1))

Investigando o gráfico de valores ajustados vs resíduos podemos observar que existe heterocedasticidade (forma de funil) e que a um padrão curvo nos rerísudos sugerindo multicolinearidade. Já o qqplot evidência oque vimos no histograma e no teste shapiro-wilk e ainda alguns possíveis pontos de alavanca nos últimos dois gráficos.

  1. O gráfico dos resíduos contra valores ajustados produz um padrão curvo. Descreva o que, se alguma coisa pode ser observada sobre o modelo neste gráfico.

Estamos diante de um problema de multicolinearidade ! Para isso vamos observar a matriz de correlação:

# carros[,c(3,5,6,7,9,10,11)] indica que estamos fazendo as operações com as colunas especificadas do data.frame carros
cor(carros[,c(3,5,6,7,9,10,11)]) %>% 
  kable("html") %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)
SuggestedRetailPrice EngineSize Cylinders Horsepower HighwayMPG Weight WheelBase
SuggestedRetailPrice 1.0000000 0.7043118 0.7623161 0.8531224 -0.5661899 0.7798208 0.6795325
EngineSize 0.7043118 1.0000000 0.9275933 0.8246932 -0.6564508 0.8447530 0.8171100
Cylinders 0.7623161 0.9275933 1.0000000 0.8405416 -0.6608553 0.8319864 0.7725715
Horsepower 0.8531224 0.8246932 0.8405416 1.0000000 -0.7189530 0.8312084 0.7283060
HighwayMPG -0.5661899 -0.6564508 -0.6608553 -0.7189530 1.0000000 -0.7605876 -0.5835760
Weight 0.7798208 0.8447530 0.8319864 0.8312084 -0.7605876 1.0000000 0.8524790
WheelBase 0.6795325 0.8171100 0.7725715 0.7283060 -0.5835760 0.8524790 1.0000000

Veja que todas as nossas variáveis explicativas são correlacionadas entre si, ocasionando este problema.

  1. Identifique quaisquer pontos de alavancagem ruins para o modelo.

Observando os gráficos de diagnóstico podemos ressaltar as seguintes observações como possíveis pontos de alavanca: 222, 223 e 67. Antes de observa-los façamos um sumário dos dados:

# Carros hibridos
carros[,c(2,3,5,6,7,9,10,11)] %>% 
  filter(Hybrid==1) %>% 
  summary
##  Hybrid SuggestedRetailPrice   EngineSize      Cylinders    
##  0:0    Min.   :19110        Min.   :1.400   Min.   :3.000  
##  1:3    1st Qu.:19625        1st Qu.:1.450   1st Qu.:3.500  
##         Median :20140        Median :1.500   Median :4.000  
##         Mean   :19920        Mean   :1.633   Mean   :3.667  
##         3rd Qu.:20325        3rd Qu.:1.750   3rd Qu.:4.000  
##         Max.   :20510        Max.   :2.000   Max.   :4.000  
##    Horsepower      HighwayMPG       Weight       WheelBase    
##  Min.   : 73.0   Min.   :51.0   Min.   :1850   Min.   : 95.0  
##  1st Qu.: 83.0   1st Qu.:51.0   1st Qu.:2291   1st Qu.: 99.0  
##  Median : 93.0   Median :51.0   Median :2732   Median :103.0  
##  Mean   : 92.0   Mean   :56.0   Mean   :2491   Mean   :101.3  
##  3rd Qu.:101.5   3rd Qu.:58.5   3rd Qu.:2811   3rd Qu.:104.5  
##  Max.   :110.0   Max.   :66.0   Max.   :2890   Max.   :106.0
# Carros nao-hibridos
carros[,c(2,3,5,6,7,9,10,11)] %>% 
  filter(Hybrid==0) %>% 
  summary
##  Hybrid  SuggestedRetailPrice   EngineSize      Cylinders     
##  0:231   Min.   : 10280       Min.   :1.500   Min.   : 4.000  
##  1:  0   1st Qu.: 19188       1st Qu.:2.200   1st Qu.: 4.000  
##          Median : 26189       Median :2.900   Median : 6.000  
##          Mean   : 29885       Mean   :2.916   Mean   : 5.541  
##          3rd Qu.: 36945       3rd Qu.:3.500   3rd Qu.: 6.000  
##          Max.   :128420       Max.   :5.500   Max.   :12.000  
##    Horsepower      HighwayMPG        Weight       WheelBase    
##  Min.   :100.0   Min.   :19.00   Min.   :2035   Min.   : 93.0  
##  1st Qu.:152.5   1st Qu.:26.00   1st Qu.:2980   1st Qu.:104.0  
##  Median :200.0   Median :28.00   Median :3416   Median :107.0  
##  Mean   :201.2   Mean   :29.05   Mean   :3324   Mean   :107.2  
##  3rd Qu.:232.0   3rd Qu.:31.00   3rd Qu.:3650   3rd Qu.:111.0  
##  Max.   :493.0   Max.   :46.00   Max.   :4474   Max.   :124.0
# Observações
carros[c(67,222,223),c(2,3,5,6,7,9,10,11)] %>% 
  kable("html") %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)
Hybrid SuggestedRetailPrice EngineSize Cylinders Horsepower HighwayMPG Weight WheelBase
67 1 19110 2.0 3 73 66 1850 95
222 0 94820 5.0 8 302 24 4085 114
223 0 128420 5.5 12 493 19 4473 114

Note que dentre os carros hibridos, a observação 67 é a que possui os valores extremos(máximos e mínimos) de todas as variáveis. No caso da observação 222 e 223, elas possuem características muito parecidas, com valores que diferem muito da média para todas as variáveis, sendo que, algumas aparecem como pontos extremos em alguns casos.

Agora, vamos tentar classificar estes pontos como bons ou ruins.

# Retirando a observao 222
mod2.222 <- lm(SuggestedRetailPrice~EngineSize + Cylinders + Horsepower +
             HighwayMPG + Weight + WheelBase + Hybrid ,
             data = carros, subset = -222)
summary(mod2.222)
## 
## Call:
## lm(formula = SuggestedRetailPrice ~ EngineSize + Cylinders + 
##     Horsepower + HighwayMPG + Weight + WheelBase + Hybrid, data = carros, 
##     subset = -222)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -16277  -4233    139   3236  35807 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -73535.431  14749.631  -4.986 1.23e-06 ***
## EngineSize   -8391.883   1471.953  -5.701 3.71e-08 ***
## Cylinders     3967.705    884.933   4.484 1.17e-05 ***
## Horsepower     178.470     14.946  11.941  < 2e-16 ***
## HighwayMPG     612.269    184.649   3.316  0.00107 ** 
## Weight          11.770      2.421   4.862 2.18e-06 ***
## WheelBase      120.117    162.501   0.739  0.46057    
## Hybrid1        420.995   5547.784   0.076  0.93958    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6860 on 225 degrees of freedom
## Multiple R-squared:  0.8059, Adjusted R-squared:  0.7998 
## F-statistic: 133.4 on 7 and 225 DF,  p-value: < 2.2e-16
mod2.223 <- lm(SuggestedRetailPrice~EngineSize + Cylinders + Horsepower +
             HighwayMPG + Weight + WheelBase + Hybrid ,
             data = carros, subset = -223)
summary(mod2.223)
## 
## Call:
## lm(formula = SuggestedRetailPrice ~ EngineSize + Cylinders + 
##     Horsepower + HighwayMPG + Weight + WheelBase + Hybrid, data = carros, 
##     subset = -223)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -16202  -3910   -115   3289  47344 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -73913.816  15322.334  -4.824 2.59e-06 ***
## EngineSize   -5568.065   1535.105  -3.627 0.000355 ***
## Cylinders     2207.195    951.775   2.319 0.021290 *  
## Horsepower     159.201     15.989   9.957  < 2e-16 ***
## HighwayMPG     561.537    192.159   2.922 0.003830 ** 
## Weight          12.689      2.517   5.041 9.50e-07 ***
## WheelBase      159.507    169.637   0.940 0.348081    
## Hybrid1        975.923   5759.176   0.169 0.865591    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7120 on 225 degrees of freedom
## Multiple R-squared:  0.7673, Adjusted R-squared:   0.76 
## F-statistic:   106 on 7 and 225 DF,  p-value: < 2.2e-16
mod2.67 <- lm(SuggestedRetailPrice~EngineSize + Cylinders + Horsepower +
             HighwayMPG + Weight + WheelBase + Hybrid ,
             data = carros, subset = -67)
summary(mod2.67)
## 
## Call:
## lm(formula = SuggestedRetailPrice ~ EngineSize + Cylinders + 
##     Horsepower + HighwayMPG + Weight + WheelBase + Hybrid, data = carros, 
##     subset = -67)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -17357  -4024    248   3336  46843 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -69971.791  16138.997  -4.336 2.19e-05 ***
## EngineSize   -7607.134   1646.639  -4.620 6.47e-06 ***
## Cylinders     3777.170    975.652   3.871 0.000142 ***
## Horsepower     177.102     16.438  10.774  < 2e-16 ***
## HighwayMPG     566.147    207.064   2.734 0.006752 ** 
## Weight          11.919      2.649   4.499 1.10e-05 ***
## WheelBase       87.773    179.277   0.490 0.624899    
## Hybrid1      -3306.840   6513.845  -0.508 0.612187    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7508 on 225 degrees of freedom
## Multiple R-squared:  0.7839, Adjusted R-squared:  0.7771 
## F-statistic: 116.6 on 7 and 225 DF,  p-value: < 2.2e-16
## Reunindo os dados para classificar os pontos de alavanca
# Reunindo os betas
coef.mod2 <- data.frame(coef.original = mod2$coefficients,
                            coef.sem222 = mod2.222$coefficients,
                            coef.sem223 = mod2.223$coefficients,
                            coef.sem67 = mod2.67$coefficients)

coef.mod2  %>% 
  kable("html") %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)
coef.original coef.sem222 coef.sem223 coef.sem67
(Intercept) -68965.79339 -73535.43079 -73913.81574 -69971.7908
EngineSize -6957.45669 -8391.88259 -5568.06463 -7607.1338
Cylinders 3564.75496 3967.70480 2207.19467 3777.1704
Horsepower 179.70211 178.46987 159.20094 177.1017
HighwayMPG 637.93940 612.26852 561.53690 566.1469
Weight 11.91076 11.76988 12.68875 11.9188
WheelBase 47.60692 120.11661 159.50672 87.7726
Hybrid1 431.75855 420.99532 975.92325 -3306.8399
# R² para cada modelo e desvio dos residuos
r2 <- data.frame(original = c(stringi::stri_paste(78.19,"%"),"7533"),
                 sem222 = c(stringi::stri_paste(80.59,"%"),"6860"),
                 sem223 = c(stringi::stri_paste(76.73,"%"),"7120"),
                 sem67 = c(stringi::stri_paste(78.39,"%"),"7508"))
rownames(r2) = c("R²","Desvio")
r2 %>% 
  kable("html") %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)
original sem222 sem223 sem67
78.19% 80.59% 76.73% 78.39%
Desvio 7533 6860 7120 7508

Investigando os modelos estimados, podemos notar que para as estimativas do intercepeto, todos os modelos sofrem alterações razoaveis. Já para a variável EngineSize sempre que é feita a retirada de alguma das observações citadas há uma expressiva diferença. Para a variável WheelBase acontece algo inusitado, ao retirarmos as observações 222 e 223 as estimativas aumentam muito, oque não era de se esperar já que estas são muito proximas do máximo. Examinando a variável Hybrid, podemos notar que ao retirar a observação 223 a um aumento gritante em sua estimativa, isso acontece pois este elemento representa o valor de máximo da variável SuggestedRetailPrice para a categoria de carros não-hibridos, já ao removermos a observação 67, o sinal da estimativa muda e ela diminui drasticamente. Uma explicação razoável pra este acontecimento é que só temos 3 obserções para esta categoria. Para os outros coeficientes nada muito gritante ocorreu.Estudando o R² e o desvio padrão dos resíduos, também não temos diferenças bruscas. Para finalizarmos o item (c) vamos explorar os gráficos de diagnóstico:

par(mfrow = c(2,2))
plot(mod2.222)

plot(mod2.223)