Seleção automática de variáveis - Laboratório 4
José Augusto Barros de Oliveira
20 de setembro de 2016
Sumário
O R apresenta diferentes soluções para as demandas estatísticas, como se sabe. Ao tentar realizar as atividades do laboratório 4 (previsto para o Minitab) no R, deparei-me com uma série de possibilidades para a seleção automática de modelos. Orientei-me inicialmente por esse material [http://www.ics.uci.edu/~jutts/st108/Model_Selection_in_R.doc]. Há uma abordagem semelhante a do Minitab (utilizando-se nível de significância e valor da estatística F parcial para inclusão/exclusão de variáveis) mas o procedimento não é automático. As opções automáticas de seleção com as funções “leaps” e “step” utilizam outros critérios de julgamento (Cp, R2 ajustado e AIC).
Uma outra abordagem foi tentada, com base nesse material [http://www-bcf.usc.edu/~gareth/ISL/] cap.6 “Subset Selection Methods”. Utilizou-se a função regsubsets(), e gráficos para a escolha do melhor modelo.
A expectativa era de que as diferentes abordagens convergissem para um mesmo modelo, coisa que não aconteceu. Então, para escolher o melhor modelo, utilizou-se uma tabela que resume as escolhas possíveis e os testes relacionados.
Seleção automática de modelos usando a função “leap”
library("leaps", lib.loc="~/R/win-library/3.2")
library("MASS", lib.loc="~/R/win-library/3.2")
cimento <- read.csv2("D:/Guto/pessoal/especializacao/regressao linear 1/EST556_CimentoPortland.csv", stringsAsFactors=FALSE)
head(cimento)## X1 X2 X3 X4 Y
## 1 7 26 6 60 78.5
## 2 1 29 15 52 74.3
## 3 11 56 8 20 104.3
## 4 11 31 8 47 87.6
## 5 7 52 6 33 95.9
## 6 11 55 9 22 109.2str(cimento)## 'data.frame': 13 obs. of 5 variables:
## $ X1: int 7 1 11 11 7 11 3 1 2 21 ...
## $ X2: int 26 29 56 31 52 55 71 31 54 47 ...
## $ X3: int 6 15 8 8 6 9 17 22 18 4 ...
## $ X4: int 60 52 20 47 33 22 6 44 22 26 ...
## $ Y : num 78.5 74.3 104.3 87.6 95.9 ...plot(cimento)
cor(cimento)## X1 X2 X3 X4 Y
## X1 1.0000000 0.2285795 -0.8241338 -0.2454451 0.7307175
## X2 0.2285795 1.0000000 -0.1392424 -0.9729550 0.8162526
## X3 -0.8241338 -0.1392424 1.0000000 0.0295370 -0.5346707
## X4 -0.2454451 -0.9729550 0.0295370 1.0000000 -0.8213050
## Y 0.7307175 0.8162526 -0.5346707 -0.8213050 1.0000000Existe forte correlação entre as preditoras e a resposta, exceto por X3. Existe também forte correlação entre X4 e X2 e entre X3 e X1.
A função “leaps” do pacote de mesmo nome procura pelo melhor conjunto de preditores usando o critério escolhido: Cp de Mallows ou R2 ajustado.
leaps(x=cimento[,1:4], y=cimento[,5], names=names(cimento[,1:4]), method="Cp")## $which
## X1 X2 X3 X4
## 1 FALSE FALSE FALSE TRUE
## 1 FALSE TRUE FALSE FALSE
## 1 TRUE FALSE FALSE FALSE
## 1 FALSE FALSE TRUE FALSE
## 2 TRUE TRUE FALSE FALSE
## 2 TRUE FALSE FALSE TRUE
## 2 FALSE FALSE TRUE TRUE
## 2 FALSE TRUE TRUE FALSE
## 2 FALSE TRUE FALSE TRUE
## 2 TRUE FALSE TRUE FALSE
## 3 TRUE TRUE FALSE TRUE
## 3 TRUE TRUE TRUE FALSE
## 3 TRUE FALSE TRUE TRUE
## 3 FALSE TRUE TRUE TRUE
## 4 TRUE TRUE TRUE TRUE
##
## $label
## [1] "(Intercept)" "X1" "X2" "X3" "X4"
##
## $size
## [1] 2 2 2 2 3 3 3 3 3 3 4 4 4 4 5
##
## $Cp
## [1] 138.730833 142.486407 202.548769 315.154284 2.678242 5.495851
## [7] 22.373112 62.437716 138.225920 198.094653 3.018233 3.041280
## [13] 3.496824 7.337474 5.000000Usando-se o critério Cp de Mallows, o 5o modelo (X1 e X2) com p = 3 e Cp = 2.678 parece ser o mais indicado.
Usando-se o R2 ajustado para a seleção do modelo, o resultado é o seguinte
leaps(x=cimento[,1:4], y=cimento[,5], names=names(cimento[,1:4]), method="adjr2")## $which
## X1 X2 X3 X4
## 1 FALSE FALSE FALSE TRUE
## 1 FALSE TRUE FALSE FALSE
## 1 TRUE FALSE FALSE FALSE
## 1 FALSE FALSE TRUE FALSE
## 2 TRUE TRUE FALSE FALSE
## 2 TRUE FALSE FALSE TRUE
## 2 FALSE FALSE TRUE TRUE
## 2 FALSE TRUE TRUE FALSE
## 2 FALSE TRUE FALSE TRUE
## 2 TRUE FALSE TRUE FALSE
## 3 TRUE TRUE FALSE TRUE
## 3 TRUE TRUE TRUE FALSE
## 3 TRUE FALSE TRUE TRUE
## 3 FALSE TRUE TRUE TRUE
## 4 TRUE TRUE TRUE TRUE
##
## $label
## [1] "(Intercept)" "X1" "X2" "X3" "X4"
##
## $size
## [1] 2 2 2 2 3 3 3 3 3 3 4 4 4 4 5
##
## $adjr2
## [1] 0.6449549 0.6359290 0.4915797 0.2209521 0.9744140 0.9669653 0.9223476
## [8] 0.8164305 0.6160725 0.4578001 0.9764473 0.9763796 0.9750415 0.9637599
## [15] 0.9735634Pelo critério do R2 ajustado, o melhor modelo seria o 11o (X1,X2 e X4), com 97,64%. Mas os modelos 5, 12, 13 e 15 também apresentam R2 ajustado acima de 97%.
Stepwise selection com critério AIC
Forward
Escolhemos X4 que é a variável que apresenta o maior coeficiente de correlação com a resposta para compor o modelo base
base <- lm(Y ~ X4, data=cimento)
full <- lm(Y ~ ., data=cimento)
step(base, scope = list( upper=full, lower=~1 ), direction = "forward", trace=TRUE)## Start: AIC=58.85
## Y ~ X4
##
## Df Sum of Sq RSS AIC
## + X1 1 809.10 74.76 28.742
## + X3 1 708.13 175.74 39.853
## <none> 883.87 58.852
## + X2 1 14.99 868.88 60.629
##
## Step: AIC=28.74
## Y ~ X4 + X1
##
## Df Sum of Sq RSS AIC
## + X2 1 26.789 47.973 24.974
## + X3 1 23.926 50.836 25.728
## <none> 74.762 28.742
##
## Step: AIC=24.97
## Y ~ X4 + X1 + X2
##
## Df Sum of Sq RSS AIC
## <none> 47.973 24.974
## + X3 1 0.10909 47.864 26.944##
## Call:
## lm(formula = Y ~ X4 + X1 + X2, data = cimento)
##
## Coefficients:
## (Intercept) X4 X1 X2
## 71.6483 -0.2365 1.4519 0.4161O resultado é coerente com a seleção por “leap” e R2 ajustado (modelo com X4, X1 e X2)
Backward
step(full, direction = "backward", trace=TRUE)## Start: AIC=26.94
## Y ~ X1 + X2 + X3 + X4
##
## Df Sum of Sq RSS AIC
## - X3 1 0.1091 47.973 24.974
## - X4 1 0.2470 48.111 25.011
## - X2 1 2.9725 50.836 25.728
## <none> 47.864 26.944
## - X1 1 25.9509 73.815 30.576
##
## Step: AIC=24.97
## Y ~ X1 + X2 + X4
##
## Df Sum of Sq RSS AIC
## <none> 47.97 24.974
## - X4 1 9.93 57.90 25.420
## - X2 1 26.79 74.76 28.742
## - X1 1 820.91 868.88 60.629##
## Call:
## lm(formula = Y ~ X1 + X2 + X4, data = cimento)
##
## Coefficients:
## (Intercept) X1 X2 X4
## 71.6483 1.4519 0.4161 -0.2365Mais uma vez, o resultado aponta para o modelo com X1, X2 e X4
Seleção manual com critério do p valor e F parcial
Forward
Vamos estabelecer alfa = 0.05 e F = 4. Usaremos como ponto de partida um modelo que contém apenas o intercepto.
modelo0 <- lm(Y ~ 1, cimento)
addterm(modelo0, scope=full, test = "F")## Single term additions
##
## Model:
## Y ~ 1
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 2715.76 71.444
## X1 1 1450.08 1265.69 63.519 12.6025 0.0045520 **
## X2 1 1809.43 906.34 59.178 21.9606 0.0006648 ***
## X3 1 776.36 1939.40 69.067 4.4034 0.0597623 .
## X4 1 1831.90 883.87 58.852 22.7985 0.0005762 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1X4 tem o maior valor de F e deve ser a primeira variável a entrar no modelo
modelo1 <- update(modelo0, .~. + X4)
addterm(modelo1, scope=full, test = "F")## Single term additions
##
## Model:
## Y ~ X4
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 883.87 58.852
## X1 1 809.10 74.76 28.742 108.224 1.105e-06 ***
## X2 1 14.99 868.88 60.629 0.172 0.6867
## X3 1 708.13 175.74 39.853 40.295 8.375e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1X1 é a próxima variável a ser incluída pois satisfaz os critérios
modelo1 <- update(modelo1, .~. + X1)
addterm(modelo1, scope=full, test = "F")## Single term additions
##
## Model:
## Y ~ X4 + X1
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 74.762 28.742
## X2 1 26.789 47.973 24.974 5.0259 0.05169 .
## X3 1 23.926 50.836 25.728 4.2358 0.06969 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Ao nível alfa = 0.05, não haveria outra variável explicativa a incluir. O modelo ficaria com X1 e X4
Backward
dropterm(full, test = "F")## Single term deletions
##
## Model:
## Y ~ X1 + X2 + X3 + X4
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 47.864 26.944
## X1 1 25.9509 73.815 30.576 4.3375 0.07082 .
## X2 1 2.9725 50.836 25.728 0.4968 0.50090
## X3 1 0.1091 47.973 24.974 0.0182 0.89592
## X4 1 0.2470 48.111 25.011 0.0413 0.84407
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1X3 deve ser retirada pois apresenta o menor valor no teste F.
modelo3 <- update(full, .~. -X3)
dropterm(modelo3, test="F")## Single term deletions
##
## Model:
## Y ~ X1 + X2 + X4
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 47.97 24.974
## X1 1 820.91 868.88 60.629 154.008 5.781e-07 ***
## X2 1 26.79 74.76 28.742 5.026 0.05169 .
## X4 1 9.93 57.90 25.420 1.863 0.20540
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1X4 deve ser retirado pois não é significante
modelo4 <- update(modelo3, .~. -X4)
dropterm(modelo4, test="F")## Single term deletions
##
## Model:
## Y ~ X1 + X2
## Df Sum of Sq RSS AIC F Value Pr(F)
## <none> 57.90 25.420
## X1 1 848.43 906.34 59.178 146.52 2.692e-07 ***
## X2 1 1207.78 1265.69 63.519 208.58 5.029e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1X1 e X2 são significantes e devem permanecer no modelo.
Seleção pelos melhores subconjuntos
A função regsubsets() retorna os melhoes subconjuntos identificando o melhor modelo que contém um determinado número de preditores, avaliando a Soma do quadrado dos resíduos.
regfull <- regsubsets(Y ~ .,cimento)
regsuma <- summary(regfull)
regsuma## Subset selection object
## Call: regsubsets.formula(Y ~ ., cimento)
## 4 Variables (and intercept)
## Forced in Forced out
## X1 FALSE FALSE
## X2 FALSE FALSE
## X3 FALSE FALSE
## X4 FALSE FALSE
## 1 subsets of each size up to 4
## Selection Algorithm: exhaustive
## X1 X2 X3 X4
## 1 ( 1 ) " " " " " " "*"
## 2 ( 1 ) "*" "*" " " " "
## 3 ( 1 ) "*" "*" " " "*"
## 4 ( 1 ) "*" "*" "*" "*"Um gráfico com SQres, R2 ajustado, Cp e BIC para todos os modelos pode ajudar na escolha do melhor modelo entre as opções dadas.
par(mfrow =c(2,2))
plot (regsuma$rss, xlab="Número de variáveis", ylab="SQres", type="l")
plot (regsuma$adjr2, xlab="Número de variáveis", ylab="R2 ajustado", type="l")
maior <- which.max(regsuma$adjr2)
points(maior,regsuma$adjr2[maior], col="red", cex=2, pch=20)
plot(regsuma$cp, xlab="Número de variáveis", ylab="Cp", type="l")
menor <- which.min(regsuma$cp)
points(menor,regsuma$cp[menor], col="red", cex=2, pch=20)
plot(regsuma$bic, xlab="Número de variáveis", ylab="BIC", type="l")
menor <- which.min(regsuma$bic)
points(menor,regsuma$bic[menor], col="red", cex=2, pch=20)
O modelo com duas explicativas (X1 e X2) parece ser o mais adequado.
Comparação entre os modelos propostos pelas diferentes abordagens
| Método | preditores | QMRes | F | p valor (F) | R2 ajust | VIFs | Cp |
|---|---|---|---|---|---|---|---|
| leaps(Cp) | X1 e X2 | 5.79 | 229 | 4.4 e-09 | 97.4% | 1.06; 1.06 | 2.7 (p=3) |
| leaps(r2) | X1, X2 e X4 | 5.33 | 166.8 | 3.3 e-08 | 97.6% | 1.06;18.7;18.9 | 3.0 (P=4) |
| step(fwd) | X1, X2 e X4 | 5.33 | 166.8 | 3.3 e-08 | 97.6% | 1.06;18.7;18.9 | 3.0 (P=4) |
| step(bwd) | X1, X2 e X4 | 5.33 | 166.8 | 3.3 e-08 | 97.6% | 1.06;18.7;18.9 | 3.0 (P=4) |
| manual(fwd) | X1 e X4 | 7.48 | 176.6 | 1.8 e-07 | 96.7% | 1.06;1.06 | 5.5 (p=3) |
| manual(bwd) | X1 e X2 | 5.79 | 229 | 4.4 e-09 | 97.4% | 1.06; 1.06 | 2.7 (p=3) |
| regsubset | X1 e X2 | 5.79 | 229 | 4.4 e-09 | 97.4% | 1.06; 1.06 | 2.7 (p=3) |
melhor <- lm(Y ~ X1 + X2, cimento)
summary(melhor)##
## Call:
## lm(formula = Y ~ X1 + X2, data = cimento)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.893 -1.574 -1.302 1.363 4.048
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 52.57735 2.28617 23.00 5.46e-10 ***
## X1 1.46831 0.12130 12.11 2.69e-07 ***
## X2 0.66225 0.04585 14.44 5.03e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.406 on 10 degrees of freedom
## Multiple R-squared: 0.9787, Adjusted R-squared: 0.9744
## F-statistic: 229.5 on 2 and 10 DF, p-value: 4.407e-09anova(melhor)## Analysis of Variance Table
##
## Response: Y
## Df Sum Sq Mean Sq F value Pr(>F)
## X1 1 1450.1 1450.08 250.43 2.088e-08 ***
## X2 1 1207.8 1207.78 208.58 5.029e-08 ***
## Residuals 10 57.9 5.79
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1O melhor modelo é, portanto, Y = 52.57 + 1.47X1 +0.66X2 A forte correlação entre as explicativas X1, X3 e X2, X4 revelada pela matriz de correlação já indicava colinearidade e consequente dificuldade de inclui-las simultaneamente no modelo.