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.
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.
\(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:
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.
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.
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 | |
|---|---|---|---|---|
| R² | 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)
plot(mod2.67)
par(mfrow = c(1,1))
Finalizando, não existe nenhuma diferença notável na retirada desses pontos do modelo proposto. Sendo assim, mesmo se tirarmos estes pontos o modelo contínua inapropriado, contudo, a retirada destes causam algumas diferenças nas estimativas dos parâmetros. Logo, podemos categoriza-los como pontos de alavanca ruins.
Faça uma análise para esse modelo, logo responda: (a) Decida se o modelo é válido. Dê fundamentos para apoiar sua resposta.
mod3 <- lm(log(SuggestedRetailPrice)~I(EngineSize^(1/4)) + log(Cylinders) + log(Horsepower) +
I(1/HighwayMPG) + Weight + log(WheelBase) + Hybrid, data = carros)
summary(mod3)
##
## Call:
## lm(formula = log(SuggestedRetailPrice) ~ I(EngineSize^(1/4)) +
## log(Cylinders) + log(Horsepower) + I(1/HighwayMPG) + Weight +
## log(WheelBase) + Hybrid, data = carros)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42288 -0.10983 -0.00203 0.10279 0.70068
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.703e+00 2.010e+00 2.838 0.00496 **
## I(EngineSize^(1/4)) -1.575e+00 3.332e-01 -4.727 4.01e-06 ***
## log(Cylinders) 2.335e-01 1.204e-01 1.940 0.05359 .
## log(Horsepower) 8.992e-01 8.876e-02 10.130 < 2e-16 ***
## I(1/HighwayMPG) 8.029e-01 4.758e+00 0.169 0.86614
## Weight 5.043e-04 6.367e-05 7.920 1.07e-13 ***
## log(WheelBase) -6.385e-02 4.715e-01 -0.135 0.89240
## Hybrid1 6.422e-01 1.150e-01 5.582 6.78e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1789 on 226 degrees of freedom
## Multiple R-squared: 0.8621, Adjusted R-squared: 0.8578
## F-statistic: 201.8 on 7 and 226 DF, p-value: < 2.2e-16
par(mfrow = c(2,3))
plot(mod3, which = c(1:5))
par(mfrow = c(1,1))
shapiro.test(mod3$residuals)
##
## Shapiro-Wilk normality test
##
## data: mod3$residuals
## W = 0.98142, p-value = 0.003655
Observando o o outpot para o novo modelo, nota-se uma melhora gritante entre este e o modelo anterior, deixando-o bem aceitável, porém, ao análisar os gráficos de diagnóstico, ainda temos algumas discrepâncias quanto as suposições para o modelo de regressão linear. Fazendo um simples teste de normalidade para os resíduos podemos constatar a quebra de uma suposição muito forte(resíduos não são normais ao nível alfa de 5%). Sendo assim, este modelo ainda não é válido.
mod3.1 <- lm(log(SuggestedRetailPrice)~I(EngineSize^(1/4)) + log(Cylinders) + log(Horsepower) + Weight + Hybrid ,
data = carros)
summary(mod3.1)
##
## Call:
## lm(formula = log(SuggestedRetailPrice) ~ I(EngineSize^(1/4)) +
## log(Cylinders) + log(Horsepower) + Weight + Hybrid, data = carros)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42224 -0.11001 -0.00099 0.10191 0.70205
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.422e+00 3.291e-01 16.474 < 2e-16 ***
## I(EngineSize^(1/4)) -1.591e+00 3.157e-01 -5.041 9.45e-07 ***
## log(Cylinders) 2.375e-01 1.186e-01 2.003 0.0463 *
## log(Horsepower) 9.049e-01 8.305e-02 10.896 < 2e-16 ***
## Weight 5.029e-04 5.203e-05 9.666 < 2e-16 ***
## Hybrid1 6.340e-01 1.080e-01 5.870 1.53e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1781 on 228 degrees of freedom
## Multiple R-squared: 0.862, Adjusted R-squared: 0.859
## F-statistic: 284.9 on 5 and 228 DF, p-value: < 2.2e-16
# Teste F parcial
anova(mod3,mod3.1)
Fazendo o teste F parcial, vemos que não rejeitamos \(H_0\), ou seja, os parametros referentes as variáveis retiradas do modelo são iguais a zero. Sendo assim, foi uma boa estrategia retirar estas variáveis.
par(mfrow = c(2,2))
plot(mod3.1)
par(mfrow = c(1,1))
Investigando os gráficos de diagnóstico, podemos ver que o problema ainda permanece. Sendo assim, é aconselhavel que não tentemos mais tratar o problema com um modelo de regressão linear, devemos utilizar algo mais sofisticado.