Exemplo 2: No arquivo consumo combustivel.txt (veja Gray (1989) e o livro do Prof. Gilbeto) são apresentadas as siglas dos 48 estados norte-americanos contíguos juntamente com as seguintes variáveis: taxa (taxa do combustível no estado), licenca (proporção de motoristas licenciados), renda (renda per-capita), estradas (ajuda federal para as estradas) e consumo (consumo de combustível por habitante).

O interesse nesse estudo é tentar explicar o consumo de combustível pelas variáaveis taxa, licenca, renda e estradas.

dados = read.table("https://docs.ufpr.br/~lucambio/CE225/1S2009/consumo_combustivel.txt", header = T)
head(dados)

Estudo descritivo

library(lattice)
par(mfrow=c(2,2), mar=c(3,3,1,0)+.5, mgp=c(1.6,.6,0), pch=19)
plot(Consumo ~ Taxa, data = dados)
grid()
plot(Consumo ~ Licenca, data = dados)
grid()
plot(Consumo ~ Renda, data = dados)
grid()
plot(Consumo ~ Estrada, data = dados)
grid()

Observamos que, quando a Taxa cresce o consumo diminu, quando a Licenca aumenta o consumo também mas não fica clara a influência da Renda e da Estrada.

Ajuste do modelo de regressão

ajuste = lm(Consumo ~ . - Estados, data = dados)
summary(ajuste)
## 
## Call:
## lm(formula = Consumo ~ . - Estados, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -122.03  -45.57  -10.66   31.53  234.95 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.773e+02  1.855e+02   2.033 0.048207 *  
## Taxa        -3.479e+01  1.297e+01  -2.682 0.010332 *  
## Licenca      1.336e+03  1.923e+02   6.950 1.52e-08 ***
## Renda       -6.659e-02  1.722e-02  -3.867 0.000368 ***
## Estrada     -2.426e-03  3.389e-03  -0.716 0.477999    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 66.31 on 43 degrees of freedom
## Multiple R-squared:  0.6787, Adjusted R-squared:  0.6488 
## F-statistic: 22.71 on 4 and 43 DF,  p-value: 3.907e-10

Fizemos desta maneira para não considerar os Estados. Escolhemos o modelo mais simples segundo o valor do AIC.

ajuste1 = step(ajuste)
## Start:  AIC=407.37
## Consumo ~ (Estados + Taxa + Licenca + Renda + Estrada) - Estados
## 
##           Df Sum of Sq    RSS    AIC
## - Estrada  1      2252 191302 405.94
## <none>                 189050 407.37
## - Taxa     1     31632 220682 412.80
## - Renda    1     65729 254779 419.69
## - Licenca  1    212355 401405 441.51
## 
## Step:  AIC=405.94
## Consumo ~ Taxa + Licenca + Renda
## 
##           Df Sum of Sq    RSS    AIC
## <none>                 191302 405.94
## - Taxa     1     33742 225044 411.74
## - Renda    1     69532 260834 418.82
## - Licenca  1    243586 434889 443.36
summary(ajuste1, correlation = T)
## 
## Call:
## lm(formula = Consumo ~ Taxa + Licenca + Renda, data = dados)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -110.10  -51.22  -12.89   24.49  238.77 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  307.32790  156.83067   1.960  0.05639 .  
## Taxa         -29.48381   10.58358  -2.786  0.00785 ** 
## Licenca     1374.76841  183.66954   7.485 2.24e-09 ***
## Renda         -0.06802    0.01701  -3.999  0.00024 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 65.94 on 44 degrees of freedom
## Multiple R-squared:  0.6749, Adjusted R-squared:  0.6527 
## F-statistic: 30.44 on 3 and 44 DF,  p-value: 8.235e-11
## 
## Correlation of Coefficients:
##         (Intercept) Taxa  Licenca
## Taxa    -0.69                    
## Licenca -0.74        0.29        
## Renda   -0.32       -0.06 -0.17

O modelo final, o menor AIC, não considera a variável Estrada e o ajuste está guardado no objeto R ajuste1.

Vejamos agora a importência relativa de cada uma das covariáveis para explicar a resposta.

library(relaimpo)
## Warning: package 'relaimpo' was built under R version 4.0.5
## Loading required package: MASS
## Loading required package: boot
## 
## Attaching package: 'boot'
## The following object is masked from 'package:lattice':
## 
##     melanoma
## Loading required package: survey
## Warning: package 'survey' was built under R version 4.0.5
## Loading required package: grid
## Loading required package: Matrix
## Loading required package: survival
## 
## Attaching package: 'survival'
## The following object is masked from 'package:boot':
## 
##     aml
## 
## Attaching package: 'survey'
## The following object is masked from 'package:graphics':
## 
##     dotchart
## Loading required package: mitools
## Warning: package 'mitools' was built under R version 4.0.5
## This is the global version of package relaimpo.
## If you are a non-US user, a version with the interesting additional metric pmvd is available
## from Ulrike Groempings web site at prof.beuth-hochschule.de/groemping.
calc.relimp(ajuste1, rela = TRUE)
## Response variable: Consumo 
## Total response variance: 12518.44 
## Analysis based on 48 observations 
## 
## 3 Regressors: 
## Taxa Licenca Renda 
## Proportion of variance explained by model: 67.49%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##               lmg
## Taxa    0.1953570
## Licenca 0.6706823
## Renda   0.1339606
## 
## Average coefficients for different model sizes: 
## 
##                    1X           2Xs           3Xs
## Taxa     -53.10629788  -42.41256600  -29.48380857
## Licenca 1409.84211133 1388.26447384 1374.76841106
## Renda     -0.04776056   -0.05878876   -0.06802286

Significa que do R^2=0.6749, 19,5% é explicado pela variável Taxa, 67% pela Licenca e 13,5% pela Renda.

O quê aconteceria se padronizamos nossas variáveis explicativas?

Consumo.p = (dados$Consumo-mean(dados$Consumo))/sd(dados$Consumo)
Taxa.p = (dados$Taxa-mean(dados$Taxa))/sd(dados$Taxa)
Licenca.p = (dados$Licenca-mean(dados$Licenca))/sd(dados$Licenca)
Renda.p = (dados$Renda-mean(dados$Renda))/sd(dados$Renda)

par(mfrow=c(2,2), mar=c(3,3,1,0)+.5, mgp=c(1.6,.6,0), pch=19)
plot(Consumo.p ~ Taxa.p, xlab = "Taxa padronizada", ylab = "Consumo")
grid()
plot(Consumo.p ~ Licenca.p, xlab = "Licença padronizada", ylab = "Consumo")
grid()
plot(Consumo.p ~ Renda.p, xlab = "Renda padronizada", ylab = "Consumo")
grid()

Ajuste do modelo novo de regressão com as variáveis padronizadas

ajuste2 = lm(Consumo.p ~ Taxa.p + Licenca.p + Renda.p)
summary(ajuste2)
## 
## Call:
## lm(formula = Consumo.p ~ Taxa.p + Licenca.p + Renda.p)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.9840 -0.4578 -0.1152  0.2189  2.1341 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.195e-16  8.506e-02   0.000  1.00000    
## Taxa.p      -2.505e-01  8.994e-02  -2.786  0.00785 ** 
## Licenca.p    6.816e-01  9.106e-02   7.485 2.24e-09 ***
## Renda.p     -3.487e-01  8.721e-02  -3.999  0.00024 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5893 on 44 degrees of freedom
## Multiple R-squared:  0.6749, Adjusted R-squared:  0.6527 
## F-statistic: 30.44 on 3 and 44 DF,  p-value: 8.235e-11

Como interpretamos os coeficiêntes da regressão? continuam não informando a importância de cada variável explicativa na resposta.

calc.relimp(ajuste2, rela = TRUE)
## Response variable: Consumo.p 
## Total response variance: 1 
## Analysis based on 48 observations 
## 
## 3 Regressors: 
## Taxa.p Licenca.p Renda.p 
## Proportion of variance explained by model: 67.49%
## Metrics are normalized to sum to 100% (rela=TRUE). 
## 
## Relative importance metrics: 
## 
##                 lmg
## Taxa.p    0.1953570
## Licenca.p 0.6706823
## Renda.p   0.1339606
## 
## Average coefficients for different model sizes: 
## 
##                   1X        2Xs        3Xs
## Taxa.p    -0.4512803 -0.3604084 -0.2505439
## Licenca.p  0.6989654  0.6882678  0.6815767
## Renda.p   -0.2448621 -0.3014022 -0.3487442

Como percebemos, a importância relativa de cada covariável não mudou!