Econometria: exemplo de seleção de modelos, Faraway (2014)

Abstract

Exercício do cap 10 do livro do Faraway (2014), Linear models with R, 2nd edition de http://www.maths.bath.ac.uk/~jjf23/LMR/scripts2/varsel.R .

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Econometria: exemplo de seleção de modelos, Faraway (2014). Campo Grande-MS,Brasil: RStudio/Rpubs, 2019 (atualizado em 2020). Disponível em http://rpubs.com/amrofi/Faraway_backward_selection.

1 Introdução

cap 10 do livro de Faraway (2014), “Linear models with R”, 2nd edition http://www.maths.bath.ac.uk/~jjf23/LMR/scripts2/varsel.R

2 Dados

Dados de 50 estados dos 1970s coletados do U.S. Census Bureau. - ‘life expectancy’ (expectativa de vida) como variável de resposta e as demais como preditoras. Ou seja, um data.frame com 50 observations para 8 variáveis.

library(faraway)
library(knitr)
data(state)
statedata <- data.frame(state.x77, row.names = state.abb)
knitr::kable(statedata)

	Population	Income	Illiteracy	Life.Exp	Murder	HS.Grad	Frost	Area
AL	3615	3624	2.1	69.05	15.1	41.3	20	50708
AK	365	6315	1.5	69.31	11.3	66.7	152	566432
AZ	2212	4530	1.8	70.55	7.8	58.1	15	113417
AR	2110	3378	1.9	70.66	10.1	39.9	65	51945
CA	21198	5114	1.1	71.71	10.3	62.6	20	156361
CO	2541	4884	0.7	72.06	6.8	63.9	166	103766
CT	3100	5348	1.1	72.48	3.1	56.0	139	4862
DE	579	4809	0.9	70.06	6.2	54.6	103	1982
FL	8277	4815	1.3	70.66	10.7	52.6	11	54090
GA	4931	4091	2.0	68.54	13.9	40.6	60	58073
HI	868	4963	1.9	73.60	6.2	61.9	0	6425
ID	813	4119	0.6	71.87	5.3	59.5	126	82677
IL	11197	5107	0.9	70.14	10.3	52.6	127	55748
IN	5313	4458	0.7	70.88	7.1	52.9	122	36097
IA	2861	4628	0.5	72.56	2.3	59.0	140	55941
KS	2280	4669	0.6	72.58	4.5	59.9	114	81787
KY	3387	3712	1.6	70.10	10.6	38.5	95	39650
LA	3806	3545	2.8	68.76	13.2	42.2	12	44930
ME	1058	3694	0.7	70.39	2.7	54.7	161	30920
MD	4122	5299	0.9	70.22	8.5	52.3	101	9891
MA	5814	4755	1.1	71.83	3.3	58.5	103	7826
MI	9111	4751	0.9	70.63	11.1	52.8	125	56817
MN	3921	4675	0.6	72.96	2.3	57.6	160	79289
MS	2341	3098	2.4	68.09	12.5	41.0	50	47296
MO	4767	4254	0.8	70.69	9.3	48.8	108	68995
MT	746	4347	0.6	70.56	5.0	59.2	155	145587
NE	1544	4508	0.6	72.60	2.9	59.3	139	76483
NV	590	5149	0.5	69.03	11.5	65.2	188	109889
NH	812	4281	0.7	71.23	3.3	57.6	174	9027
NJ	7333	5237	1.1	70.93	5.2	52.5	115	7521
NM	1144	3601	2.2	70.32	9.7	55.2	120	121412
NY	18076	4903	1.4	70.55	10.9	52.7	82	47831
NC	5441	3875	1.8	69.21	11.1	38.5	80	48798
ND	637	5087	0.8	72.78	1.4	50.3	186	69273
OH	10735	4561	0.8	70.82	7.4	53.2	124	40975
OK	2715	3983	1.1	71.42	6.4	51.6	82	68782
OR	2284	4660	0.6	72.13	4.2	60.0	44	96184
PA	11860	4449	1.0	70.43	6.1	50.2	126	44966
RI	931	4558	1.3	71.90	2.4	46.4	127	1049
SC	2816	3635	2.3	67.96	11.6	37.8	65	30225
SD	681	4167	0.5	72.08	1.7	53.3	172	75955
TN	4173	3821	1.7	70.11	11.0	41.8	70	41328
TX	12237	4188	2.2	70.90	12.2	47.4	35	262134
UT	1203	4022	0.6	72.90	4.5	67.3	137	82096
VT	472	3907	0.6	71.64	5.5	57.1	168	9267
VA	4981	4701	1.4	70.08	9.5	47.8	85	39780
WA	3559	4864	0.6	71.72	4.3	63.5	32	66570
WV	1799	3617	1.4	69.48	6.7	41.6	100	24070
WI	4589	4468	0.7	72.48	3.0	54.5	149	54464
WY	376	4566	0.6	70.29	6.9	62.9	173	97203

names(statedata)

[1] "Population" "Income"     "Illiteracy" "Life.Exp"   "Murder"    
[6] "HS.Grad"    "Frost"      "Area"      

# [1] 'Population' 'Income' 'Illiteracy' 'Life.Exp' 'Murder' [6] 'HS.Grad'
# 'Frost' 'Area'

3 Backward Selection

Estima-se o modelo completo - com todas as variáveis do data.frame - e depois segue retirando e atualizando o modelo anterior (lmod) por meio da função update. O pesquisador deve tomar cuidado para observar que o lmod é “atualizado” e portanto, sobreposto com os resultados do update. Caso deseje guardar os resultados de cada etapa, deverá modificar o nome do objeto lmod a cada atualização.

Fiz uma adaptação do script de Faraway (2014) a fim de obter os critérios de informação AIC e resultados completos do summary.

# MODELO COMPLETO
lmod <- lm(Life.Exp ~ ., statedata)
summary(lmod)

Call:
lm(formula = Life.Exp ~ ., data = statedata)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.48895 -0.51232 -0.02747  0.57002  1.49447 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.094e+01  1.748e+00  40.586  < 2e-16 ***
Population   5.180e-05  2.919e-05   1.775   0.0832 .  
Income      -2.180e-05  2.444e-04  -0.089   0.9293    
Illiteracy   3.382e-02  3.663e-01   0.092   0.9269    
Murder      -3.011e-01  4.662e-02  -6.459 8.68e-08 ***
HS.Grad      4.893e-02  2.332e-02   2.098   0.0420 *  
Frost       -5.735e-03  3.143e-03  -1.825   0.0752 .  
Area        -7.383e-08  1.668e-06  -0.044   0.9649    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7448 on 42 degrees of freedom
Multiple R-squared:  0.7362,    Adjusted R-squared:  0.6922 
F-statistic: 16.74 on 7 and 42 DF,  p-value: 2.534e-10

AIC(lmod)

[1] 121.7092

A partir desses resultados, inicia-se a retirada das variáveis não significativas uma a uma e observam-se os resultados do novo modelo “atualizado” de modo que se tenha menores AIC e maiores R quadrado ajustado. A sequência sugerida é retirar a variável de maior p-value a cada passo. Ou seja, a primeira a ser retirada é a variável ‘Area’, com p-value = 0.9649 na última regressão. Os valores para atentar são: Adjusted R-squared: 0.6922 e AIC = 121.7092.

3.1 Retirando Area

lmod <- update(lmod, . ~ . - Area)
summary(lmod)

Call:
lm(formula = Life.Exp ~ Population + Income + Illiteracy + Murder + 
    HS.Grad + Frost, data = statedata)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.49047 -0.52533 -0.02546  0.57160  1.50374 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.099e+01  1.387e+00  51.165  < 2e-16 ***
Population   5.188e-05  2.879e-05   1.802   0.0785 .  
Income      -2.444e-05  2.343e-04  -0.104   0.9174    
Illiteracy   2.846e-02  3.416e-01   0.083   0.9340    
Murder      -3.018e-01  4.334e-02  -6.963 1.45e-08 ***
HS.Grad      4.847e-02  2.067e-02   2.345   0.0237 *  
Frost       -5.776e-03  2.970e-03  -1.945   0.0584 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7361 on 43 degrees of freedom
Multiple R-squared:  0.7361,    Adjusted R-squared:  0.6993 
F-statistic: 19.99 on 6 and 43 DF,  p-value: 5.362e-11

AIC(lmod)

[1] 119.7116

Neste caso, observa-se a queda de AIC, passando agora para 119.7116, e melhora em Adjusted R-squared: 0.6993. Ou seja, o modelo melhorou retirando a variável.

3.2 Retirando Illiteracy

lmod <- update(lmod, . ~ . - Illiteracy)
summary(lmod)

Call:
lm(formula = Life.Exp ~ Population + Income + Murder + HS.Grad + 
    Frost, data = statedata)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4892 -0.5122 -0.0329  0.5645  1.5166 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.107e+01  1.029e+00  69.067  < 2e-16 ***
Population   5.115e-05  2.709e-05   1.888   0.0657 .  
Income      -2.477e-05  2.316e-04  -0.107   0.9153    
Murder      -3.000e-01  3.704e-02  -8.099 2.91e-10 ***
HS.Grad      4.776e-02  1.859e-02   2.569   0.0137 *  
Frost       -5.910e-03  2.468e-03  -2.395   0.0210 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7277 on 44 degrees of freedom
Multiple R-squared:  0.7361,    Adjusted R-squared:  0.7061 
F-statistic: 24.55 on 5 and 44 DF,  p-value: 1.019e-11

AIC(lmod)

[1] 117.7196

Neste caso, observa-se a queda de AIC, passando agora para 117.7196, e melhora em Adjusted R-squared: 0.7061. Ou seja, o modelo melhorou retirando a variável.

3.3 Retirando Income

lmod <- update(lmod, . ~ . - Income)
summary(lmod)

Call:
lm(formula = Life.Exp ~ Population + Murder + HS.Grad + Frost, 
    data = statedata)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.47095 -0.53464 -0.03701  0.57621  1.50683 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  7.103e+01  9.529e-01  74.542  < 2e-16 ***
Population   5.014e-05  2.512e-05   1.996  0.05201 .  
Murder      -3.001e-01  3.661e-02  -8.199 1.77e-10 ***
HS.Grad      4.658e-02  1.483e-02   3.142  0.00297 ** 
Frost       -5.943e-03  2.421e-03  -2.455  0.01802 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7197 on 45 degrees of freedom
Multiple R-squared:  0.736, Adjusted R-squared:  0.7126 
F-statistic: 31.37 on 4 and 45 DF,  p-value: 1.696e-12

AIC(lmod)

[1] 115.7326

Neste caso, observa-se a queda de AIC, passando agora para 115.7326, e melhora em Adjusted R-squared: 0.7126. Ou seja, o modelo melhorou retirando a variável.

3.4 Retirando Population

lmod <- update(lmod, . ~ . - Population)
summary(lmod)

Call:
lm(formula = Life.Exp ~ Murder + HS.Grad + Frost, data = statedata)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.5015 -0.5391  0.1014  0.5921  1.2268 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 71.036379   0.983262  72.246  < 2e-16 ***
Murder      -0.283065   0.036731  -7.706 8.04e-10 ***
HS.Grad      0.049949   0.015201   3.286  0.00195 ** 
Frost       -0.006912   0.002447  -2.824  0.00699 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7427 on 46 degrees of freedom
Multiple R-squared:  0.7127,    Adjusted R-squared:  0.6939 
F-statistic: 38.03 on 3 and 46 DF,  p-value: 1.634e-12

AIC(lmod)

[1] 117.9743

Neste caso, observa-se aumento de AIC, passando agora para 117.9743, e piora em Adjusted R-squared: 0.6939. Ou seja, o modelo piorou retirando a variável.

3.5 Cuidado!!!

Uma atenção deve ser dada ao fato de que as vezes a variável retirada seria significativa em outro modelo. Por exemplo, veja que “Illiteracy” agora é significativa, embora seja um modelo de AIC maior que os anteriores. Portanto, alguns autores sugerem a inclusão e eliminação de variáveis em diferentes combinações de entrada e saída do modelo. Como na próxima seção.

summary(lm(Life.Exp ~ Illiteracy + Murder + Frost, statedata))

Call:
lm(formula = Life.Exp ~ Illiteracy + Murder + Frost, data = statedata)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.59010 -0.46961  0.00394  0.57060  1.92292 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 74.556717   0.584251 127.611  < 2e-16 ***
Illiteracy  -0.601761   0.298927  -2.013  0.04998 *  
Murder      -0.280047   0.043394  -6.454 6.03e-08 ***
Frost       -0.008691   0.002959  -2.937  0.00517 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7911 on 46 degrees of freedom
Multiple R-squared:  0.6739,    Adjusted R-squared:  0.6527 
F-statistic: 31.69 on 3 and 46 DF,  p-value: 2.915e-11

AIC(lm(Life.Exp ~ Illiteracy + Murder + Frost, statedata))

[1] 124.2947

4 Alternativa

Farei aqui com armazenamento dos objetos de regressão a cada update.

# alternativa de análise
library(faraway)
data(state)
statedata <- data.frame(state.x77, row.names = state.abb)
lmod1 <- lm(Life.Exp ~ ., statedata)
lmod2 <- update(lmod1, . ~ . - Area)
lmod3 <- update(lmod2, . ~ . - Illiteracy)
lmod4 <- update(lmod3, . ~ . - Income)
lmod5 <- update(lmod4, . ~ . - Population)
alternate <- lm(Life.Exp ~ Illiteracy + Murder + Frost, statedata)
lmod1$AIC <- AIC(lmod1)
lmod2$AIC <- AIC(lmod2)
lmod3$AIC <- AIC(lmod3)
lmod4$AIC <- AIC(lmod4)
lmod5$AIC <- AIC(lmod5)
alternate$AIC <- AIC(alternate)
lmod1$BIC <- BIC(lmod1)
lmod2$BIC <- BIC(lmod2)
lmod3$BIC <- BIC(lmod3)
lmod4$BIC <- BIC(lmod4)
lmod5$BIC <- BIC(lmod5)
alternate$BIC <- BIC(alternate)
library(stargazer)
stargazer(lmod1, lmod2, lmod3, lmod4, lmod5, alternate, title = "Título: Resultados da Seleção", 
    align = TRUE, type = "text", style = "all", keep.stat = c("aic", "bic", "rsq", 
        "adj.rsq", "n"))

Título: Resultados da Seleção
=======================================================================================
                                            Dependent variable:                        
                    -------------------------------------------------------------------
                                                 Life.Exp                              
                        (1)        (2)        (3)        (4)        (5)         (6)    
---------------------------------------------------------------------------------------
Population            0.0001*    0.0001*    0.0001*    0.0001*                         
                     (0.00003)  (0.00003)  (0.00003)  (0.00003)                        
                     t = 1.775  t = 1.802  t = 1.888  t = 1.996                        
                     p = 0.084  p = 0.079  p = 0.066  p = 0.053                        
Income               -0.00002    -0.00002   -0.00002                                   
                     (0.0002)    (0.0002)   (0.0002)                                   
                    t = -0.089  t = -0.104 t = -0.107                                  
                     p = 0.930  p = 0.918  p = 0.916                                   
Illiteracy             0.034      0.028                                      -0.602**  
                      (0.366)    (0.342)                                      (0.299)  
                     t = 0.092  t = 0.083                                   t = -2.013 
                     p = 0.927  p = 0.934                                    p = 0.050 
Murder               -0.301***  -0.302***  -0.300***  -0.300***  -0.283***   -0.280*** 
                      (0.047)    (0.043)    (0.037)    (0.037)    (0.037)     (0.043)  
                    t = -6.459  t = -6.963 t = -8.099 t = -8.199 t = -7.706 t = -6.454 
                    p = 0.00000 p = 0.000  p = 0.000  p = 0.000  p = 0.000  p = 0.00000
HS.Grad               0.049**    0.048**    0.048**    0.047***   0.050***             
                      (0.023)    (0.021)    (0.019)    (0.015)    (0.015)              
                     t = 2.098  t = 2.345  t = 2.569  t = 3.142  t = 3.286             
                     p = 0.042  p = 0.024  p = 0.014  p = 0.003  p = 0.002             
Frost                 -0.006*    -0.006*    -0.006**   -0.006**  -0.007***   -0.009*** 
                      (0.003)    (0.003)    (0.002)    (0.002)    (0.002)     (0.003)  
                    t = -1.825  t = -1.945 t = -2.395 t = -2.455 t = -2.824 t = -2.937 
                     p = 0.076  p = 0.059  p = 0.021  p = 0.019  p = 0.007   p = 0.006 
Area                 -0.00000                                                          
                     (0.00000)                                                         
                    t = -0.044                                                         
                     p = 0.965                                                         
Constant             70.943***  70.989***  71.066***  71.027***  71.036***   74.557*** 
                      (1.748)    (1.387)    (1.029)    (0.953)    (0.983)     (0.584)  
                    t = 40.586  t = 51.165 t = 69.067 t = 74.542 t = 72.246 t = 127.611
                     p = 0.000  p = 0.000  p = 0.000  p = 0.000  p = 0.000   p = 0.000 
---------------------------------------------------------------------------------------
Observations            50          50         50         50         50         50     
R2                     0.736      0.736      0.736      0.736      0.713       0.674   
Adjusted R2            0.692      0.699      0.706      0.713      0.694       0.653   
Akaike Inf. Crit.     121.709    119.712    117.720    115.733    117.974     124.295  
Bayesian Inf. Crit.   138.917    135.008    131.104    127.205    127.534     133.855  
=======================================================================================
Note:                                                       *p<0.1; **p<0.05; ***p<0.01

Nesta tabela, percebe-se que o melhor modelo foi o (4), ou seja, o mesmo que obtivemos na rotina anterior, contendo os regressores: Population, Murder, HS.Grad e Frost. Foram os menores Akaike e Bayesian e o melhor \(R^2\) Ajustado.

5 Escolha de Modelos baseados em critérios

A função leaps::regsubsets faz a seleção dos modelos por busca exaustiva para frente e para trás, stepwise, ou reposição sequencial. Observe que nos plots de AIC para cada modelo, deseja-se um modelo de menor AIC. O slot com as informações de quais variáveis estão em cada modelo tem nome which dentro do summary do objeto de regsubsets, neste caso chamado de “b” (o objeto resultante do summary foi chamado de rs e dentro dele estará também o resultado do BIC de cada modelo).

require(leaps)
b <- regsubsets(Life.Exp ~ ., data = statedata)
rs <- summary(b)
rs$which

  (Intercept) Population Income Illiteracy Murder HS.Grad Frost  Area
1        TRUE      FALSE  FALSE      FALSE   TRUE   FALSE FALSE FALSE
2        TRUE      FALSE  FALSE      FALSE   TRUE    TRUE FALSE FALSE
3        TRUE      FALSE  FALSE      FALSE   TRUE    TRUE  TRUE FALSE
4        TRUE       TRUE  FALSE      FALSE   TRUE    TRUE  TRUE FALSE
5        TRUE       TRUE   TRUE      FALSE   TRUE    TRUE  TRUE FALSE
6        TRUE       TRUE   TRUE       TRUE   TRUE    TRUE  TRUE FALSE
7        TRUE       TRUE   TRUE       TRUE   TRUE    TRUE  TRUE  TRUE

AIC <- 50 * log(rs$rss/50) + (2:8) * 2
plot(AIC ~ I(1:7), ylab = "AIC", xlab = "Number of Predictors")

rs$bic

[1] -39.22051 -42.62472 -46.70678 -47.03640 -43.13738 -39.23342 -35.32373

# plot do BIC
plot(rs$bic ~ I(1:7), ylab = "BIC", xlab = "Número de Preditores")

plot(2:8, rs$adjr2, xlab = "No. of Parameters", ylab = "Adjusted R-square")

which.max(rs$adjr2)

[1] 4

plot(2:8, rs$cp, xlab = "No. of Parameters", ylab = "Cp Statistic")
abline(0, 1)

Em todos os casos, o melhor modelo foi o (4), com os preditores Population, Murder, HS.Grad e Frost.

6 Stepwise

O Stepwise Regression é uma combinação de eliminação backward (para trás) e seleção forward (para frente - inclusão). É a situação em que variáveis são adicionadas e removidas no processo e a cada estágio existem variações diversas de como proceder. O usual é minimizar AIC ou BIC. A função será step de um modelo de regressão, dentro do pacote stats que já vem no R básico. A função step retornará os vários modelos estimados e respectivos AIC até otimizar.

lmod <- lm(Life.Exp ~ ., data = statedata)
step(lmod)

Start:  AIC=-22.18
Life.Exp ~ Population + Income + Illiteracy + Murder + HS.Grad + 
    Frost + Area

             Df Sum of Sq    RSS     AIC
- Area        1    0.0011 23.298 -24.182
- Income      1    0.0044 23.302 -24.175
- Illiteracy  1    0.0047 23.302 -24.174
<none>                    23.297 -22.185
- Population  1    1.7472 25.044 -20.569
- Frost       1    1.8466 25.144 -20.371
- HS.Grad     1    2.4413 25.738 -19.202
- Murder      1   23.1411 46.438  10.305

Step:  AIC=-24.18
Life.Exp ~ Population + Income + Illiteracy + Murder + HS.Grad + 
    Frost

             Df Sum of Sq    RSS     AIC
- Illiteracy  1    0.0038 23.302 -26.174
- Income      1    0.0059 23.304 -26.170
<none>                    23.298 -24.182
- Population  1    1.7599 25.058 -22.541
- Frost       1    2.0488 25.347 -21.968
- HS.Grad     1    2.9804 26.279 -20.163
- Murder      1   26.2721 49.570  11.569

Step:  AIC=-26.17
Life.Exp ~ Population + Income + Murder + HS.Grad + Frost

             Df Sum of Sq    RSS     AIC
- Income      1     0.006 23.308 -28.161
<none>                    23.302 -26.174
- Population  1     1.887 25.189 -24.280
- Frost       1     3.037 26.339 -22.048
- HS.Grad     1     3.495 26.797 -21.187
- Murder      1    34.739 58.041  17.456

Step:  AIC=-28.16
Life.Exp ~ Population + Murder + HS.Grad + Frost

             Df Sum of Sq    RSS     AIC
<none>                    23.308 -28.161
- Population  1     2.064 25.372 -25.920
- Frost       1     3.122 26.430 -23.877
- HS.Grad     1     5.112 28.420 -20.246
- Murder      1    34.816 58.124  15.528

Call:
lm(formula = Life.Exp ~ Population + Murder + HS.Grad + Frost, 
    data = statedata)

Coefficients:
(Intercept)   Population       Murder      HS.Grad        Frost  
  7.103e+01    5.014e-05   -3.001e-01    4.658e-02   -5.943e-03  

O melhor modelo foi também o de Population, Frost, HS.Grad e Murder.

A partir da função lm.influence pode-se encontrar as observações mais influentes no modelo.

h <- lm.influence(lmod)$hat
names(h) <- state.abb
rev(sort(h))

        AK         CA         HI         NV         NM         TX         NY 
0.80952226 0.40885689 0.37876168 0.36524576 0.32472157 0.28416353 0.25694978 
        WA         OR         ND         LA         CT         UT         RI 
0.22268196 0.22183076 0.21968606 0.19494860 0.19362508 0.19097575 0.17082010 
        MD         AL         AZ         MS         FL         IL         PA 
0.16407127 0.16110483 0.15890792 0.15571602 0.14856937 0.13743337 0.13209846 
        NJ         SD         ME         SC         MI         WY         VT 
0.12993167 0.12528365 0.12181903 0.11757971 0.11725238 0.11638093 0.11504122 
        MA         GA         MO         KY         AR         CO         WV 
0.11274391 0.11028889 0.11026698 0.10909216 0.10447656 0.10250875 0.09964990 
        NC         DE         NH         ID         OH         MT         MN 
0.09361765 0.09322149 0.08980829 0.08755753 0.08751277 0.08638978 0.07617370 
        TN         WI         VA         OK         IA         NE         KS 
0.07006383 0.06835416 0.06393987 0.06349971 0.06200255 0.05749338 0.05537644 
        IN 
0.05198210 

Neste caso, fazendo uso do resultado da observação influente, percebe-se que a combinação ótima de variáveis é diferente daquela obtida com backward selection padrão. Aqui, o melhor modelo foi olhando a observação AK = Alaska. O resultado indicou a presença dos regressores Population, Murder, HS.Grad, Frost e Area. Observar que o Alaska é um grande (extenso) estado - o dobro do tamanho do Texas (segundo em area).

b <- regsubsets(Life.Exp ~ ., data = statedata, subset = (state.abb != "AK"))
rs <- summary(b)
rs$which[which.max(rs$adjr), ]

(Intercept)  Population      Income  Illiteracy      Murder     HS.Grad 
       TRUE        TRUE       FALSE       FALSE        TRUE        TRUE 
      Frost        Area 
       TRUE        TRUE 

O stripchart é uma forma de avaliar as estimativas, embora a nuvem de pontos possa gerar alguma confusão na interpretação.

stripchart(data.frame(scale(statedata)), method = "jitter", las = 2, vertical = TRUE)

Agora usaremos uma selação para trás com algumas variáveis em logs.

b <- regsubsets(Life.Exp ~ log(Population) + Income + Illiteracy + Murder + HS.Grad + 
    Frost + log(Area), statedata)
rs <- summary(b)
rs$which[which.max(rs$adjr), ]

    (Intercept) log(Population)          Income      Illiteracy          Murder 
           TRUE            TRUE           FALSE           FALSE            TRUE 
        HS.Grad           Frost       log(Area) 
           TRUE            TRUE           FALSE 

rs$bic

[1] -39.22051 -42.62472 -47.17452 -47.87315 -44.43397 -40.64452 -36.73638

rs$adjr2

[1] 0.6015893 0.6484991 0.6967729 0.7173392 0.7136360 0.7076938 0.7007574

Neste caso, a melhor estimação foi com Population, Murder, HS.Grad e Frost. O BIC para este modelo foi de -47.87315, contra -47.03640 no modelo sem logs.

7 Considerações finais

A análise considera a estimação e os valores de critérios. Ou seja, não foram checados os pressupostos clássicos de regressão, o que implica em observação de outras características residuais e eventuais outliers.

Referências

FARAWAY, Julian J. Linear Models with R. London: CRC Press/Taylor and Francis Group. 2014. Second Edition. Code available at: https://people.bath.ac.uk/jjf23/LMR/