ステップワイズ法による変数選択

みんなのR（第二版）P374にステップワイズ法に
関する記述があるのでメモします。

ステップワイズ法とは説明変数の増減を繰り返して
AICやBICなどの情報基準量やF値を評価し、
最適な変数を探す手法です。

データセットは引き続きAmerican Community Surveyを使用します。実際は説明変数のスケールや単位が異なるので、
scale関数で説明変数を標準化した方が良いです。

#American Community SurveyのデータセットをDL
acs <- read.table("http://jaredlander.com/data/acs_ny.csv",sep=",",header=TRUE)

#dplyr　Income列を追加。15000以上はTRUEとする
acs %>% 
  mutate(Income = if_else(FamilyIncome >= 15000, true = TRUE, false = FALSE)) ->acs

#ステップワイズ法でロジスティック回帰
#family = binomialにする
fullmodel <- glm(Income~NumBedrooms +NumChildren + NumPeople + NumRooms + NumUnits + NumVehicles + NumWorkers + OwnRent +YearBuilt +
                ElectricBill+FoodStamp+HeatingFuel+Insurance+Language +1 ,data=acs,family = binomial)
nullmodel <- glm(Income~1,data=acs,family = binomial)

stepwisemodel <- step(nullmodel,scope = list(lower = nullmodel,upper = fullmodel),direction = "both")

summary(stepwisemodel)

# モデルを使用して予測
probabilities <- predict(stepwisemodel, acs, type = "response")
predicted.classes <- ifelse(probabilities > 0.5, TRUE, FALSE)

# アキュラシーを計算
observed.classes <- acs$Income
mean(predicted.classes == observed.classes)

結果は

Call:
glm(formula = Income ~ NumWorkers + FoodStamp + OwnRent + NumVehicles + 
    NumChildren + Insurance + NumRooms + NumUnits + Language + 
    NumPeople, family = binomial, data = acs)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.7877   0.0717   0.1085   0.1925   1.8713  

Coefficients:
                          Estimate Std. Error z value Pr(>|z|)    
(Intercept)             -3.022e-01  3.144e-01  -0.961 0.336547    
NumWorkers               1.292e+00  6.575e-02  19.652  < 2e-16 ***
FoodStampYes            -1.441e+00  1.036e-01 -13.912  < 2e-16 ***
OwnRentOutright          1.351e+01  1.787e+02   0.076 0.939753    
OwnRentRented           -7.896e-01  1.233e-01  -6.405 1.50e-10 ***
NumVehicles              4.064e-01  6.037e-02   6.731 1.68e-11 ***
NumChildren             -2.916e-01  5.346e-02  -5.455 4.91e-08 ***
Insurance                3.266e-04  8.860e-05   3.686 0.000227 ***
NumRooms                 6.967e-02  2.473e-02   2.818 0.004837 ** 
NumUnitsSingle attached  5.770e-01  1.859e-01   3.103 0.001915 ** 
NumUnitsSingle detached  5.592e-01  1.572e-01   3.558 0.000373 ***
LanguageEnglish          7.672e-01  2.040e-01   3.761 0.000169 ***
LanguageOther            9.153e-01  4.241e-01   2.158 0.030922 *  
LanguageOther European   5.150e-01  2.429e-01   2.120 0.034011 *  
LanguageSpanish          4.708e-01  2.295e-01   2.052 0.040217 *  
NumPeople                1.224e-01  5.176e-02   2.364 0.018065 *  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 6375.4  on 22744  degrees of freedom
Residual deviance: 4200.3  on 22729  degrees of freedom
AIC: 4232.3

Number of Fisher Scoring iterations: 15

> mean(predicted.classes == observed.classes)
[1] 0.9700154

BICを基準に変数を選択してみる。

#BICを基準にモデルを選定
stepwisemodel_BIC <- step(nullmodel,scope = list(lower = nullmodel,upper = fullmodel),direction = "both",criterion = "BIC",k=log(nrow(acs)))

summary(stepwisemodel_BIC)
coef(stepwisemodel_BIC)

#モデルを使用して予測
probabilities_BIC <- predict(stepwisemodel_BIC, acs, type = "response")
predicted.classes_BIC <- ifelse(probabilities_BIC > 0.5, TRUE, FALSE)

# アキュラシーを計算
observed.classes <- acs$Income
mean(predicted.classes_BIC == observed.classes)

結果は

Call:
glm(formula = Income ~ NumWorkers + FoodStamp + OwnRent + NumVehicles + 
    NumChildren + Insurance + NumRooms, family = binomial, data = acs)

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-3.8046   0.0731   0.1104   0.1956   1.8454  

Coefficients:
                  Estimate Std. Error z value Pr(>|z|)    
(Intercept)      9.725e-01  1.961e-01   4.958 7.12e-07 ***
NumWorkers       1.293e+00  6.359e-02  20.334  < 2e-16 ***
FoodStampYes    -1.405e+00  9.973e-02 -14.087  < 2e-16 ***
OwnRentOutright  1.361e+01  1.779e+02   0.076 0.939039    
OwnRentRented   -8.035e-01  1.223e-01  -6.569 5.07e-11 ***
NumVehicles      4.257e-01  5.805e-02   7.333 2.26e-13 ***
NumChildren     -1.972e-01  3.278e-02  -6.018 1.77e-09 ***
Insurance        3.375e-04  8.890e-05   3.797 0.000147 ***
NumRooms         9.851e-02  2.428e-02   4.058 4.96e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 6375.4  on 22744  degrees of freedom
Residual deviance: 4231.0  on 22736  degrees of freedom
AIC: 4249

Number of Fisher Scoring iterations: 15

> mean(predicted.classes_BIC == observed.classes)
[1] 0.9699714