Rで重回帰分析や多重ロジスティック回帰分析などlm()やglm()でつくった回帰式の独立変数の取捨選択を行うstep()関数をつかって、任意の一部の独立変数を強制投入し、それ以外の任意の独立変数をステップワイズでふるいにかける方法をメモしておく。
やり方としては、従属変数に対して強制投入したい独立変数だけでなるモデル式をつくって、「これが最小モデルだよ」と指定してやるとできる。
本稿の内容はstep()で任意の独立変数を強制投入してステップワイズする方法だけで、増減法がいいとか減増法がいいとか、総当たりがいいとか、そもそもステップワイズは……という点については触れない。
例として、mtcarsのデータの場合でみてみる。
燃費を従属変数としてその他の値を独立変数とし、最大モデルからステップワイズで独立変数を取捨選択する場合は次のようになる。
step(lm(mpg ~ ., data = mtcars))
## Start: AIC=70.9
## mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb
##
## Df Sum of Sq RSS AIC
## - cyl 1 0.0799 147.57 68.915
## - vs 1 0.1601 147.66 68.932
## - carb 1 0.4067 147.90 68.986
## - gear 1 1.3531 148.85 69.190
## - drat 1 1.6270 149.12 69.249
## - disp 1 3.9167 151.41 69.736
## - hp 1 6.8399 154.33 70.348
## - qsec 1 8.8641 156.36 70.765
## <none> 147.49 70.898
## - am 1 10.5467 158.04 71.108
## - wt 1 27.0144 174.51 74.280
##
## Step: AIC=68.92
## mpg ~ disp + hp + drat + wt + qsec + vs + am + gear + carb
##
## Df Sum of Sq RSS AIC
## - vs 1 0.2685 147.84 66.973
## - carb 1 0.5201 148.09 67.028
## - gear 1 1.8211 149.40 67.308
## - drat 1 1.9826 149.56 67.342
## - disp 1 3.9009 151.47 67.750
## - hp 1 7.3632 154.94 68.473
## <none> 147.57 68.915
## - qsec 1 10.0933 157.67 69.032
## - am 1 11.8359 159.41 69.384
## - wt 1 27.0280 174.60 72.297
##
## Step: AIC=66.97
## mpg ~ disp + hp + drat + wt + qsec + am + gear + carb
##
## Df Sum of Sq RSS AIC
## - carb 1 0.6855 148.53 65.121
## - gear 1 2.1437 149.99 65.434
## - drat 1 2.2139 150.06 65.449
## - disp 1 3.6467 151.49 65.753
## - hp 1 7.1060 154.95 66.475
## <none> 147.84 66.973
## - am 1 11.5694 159.41 67.384
## - qsec 1 15.6830 163.53 68.200
## - wt 1 27.3799 175.22 70.410
##
## Step: AIC=65.12
## mpg ~ disp + hp + drat + wt + qsec + am + gear
##
## Df Sum of Sq RSS AIC
## - gear 1 1.565 150.09 63.457
## - drat 1 1.932 150.46 63.535
## <none> 148.53 65.121
## - disp 1 10.110 158.64 65.229
## - am 1 12.323 160.85 65.672
## - hp 1 14.826 163.35 66.166
## - qsec 1 26.408 174.94 68.358
## - wt 1 69.127 217.66 75.350
##
## Step: AIC=63.46
## mpg ~ disp + hp + drat + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## - drat 1 3.345 153.44 62.162
## - disp 1 8.545 158.64 63.229
## <none> 150.09 63.457
## - hp 1 13.285 163.38 64.171
## - am 1 20.036 170.13 65.466
## - qsec 1 25.574 175.67 66.491
## - wt 1 67.572 217.66 73.351
##
## Step: AIC=62.16
## mpg ~ disp + hp + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## - disp 1 6.629 160.07 61.515
## <none> 153.44 62.162
## - hp 1 12.572 166.01 62.682
## - qsec 1 26.470 179.91 65.255
## - am 1 32.198 185.63 66.258
## - wt 1 69.043 222.48 72.051
##
## Step: AIC=61.52
## mpg ~ hp + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## - hp 1 9.219 169.29 61.307
## <none> 160.07 61.515
## - qsec 1 20.225 180.29 63.323
## - am 1 25.993 186.06 64.331
## - wt 1 78.494 238.56 72.284
##
## Step: AIC=61.31
## mpg ~ wt + qsec + am
##
## Df Sum of Sq RSS AIC
## <none> 169.29 61.307
## - am 1 26.178 195.46 63.908
## - qsec 1 109.034 278.32 75.217
## - wt 1 183.347 352.63 82.790
##
## Call:
## lm(formula = mpg ~ wt + qsec + am, data = mtcars)
##
## Coefficients:
## (Intercept) wt qsec am
## 9.618 -3.917 1.226 2.936
これで、最終的に車体重量と1/4マイルのタイムとトランスミッションがモデルとして選択され、車体重量が軽いほうが、1/4マイルのタイムが遅い方が、トランスミッションがマニュアルのほうが燃費がよいというモデルが抽出される。
ここで、「いや、燃費にエンジンの大きさ（気筒数）関係してるでしょ」と思って、cylを強制投入したかったら、最小モデルとしてmpg ~ cylをscopeにlowerとして指定すると、これは必ず選択される。
step(lm(mpg ~ ., data = mtcars), scope = list(lower = mpg ~ cyl))
## Start: AIC=70.9
## mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb
##
## Df Sum of Sq RSS AIC
## - vs 1 0.1601 147.66 68.932
## - carb 1 0.4067 147.90 68.986
## - gear 1 1.3531 148.85 69.190
## - drat 1 1.6270 149.12 69.249
## - disp 1 3.9167 151.41 69.736
## - hp 1 6.8399 154.33 70.348
## - qsec 1 8.8641 156.36 70.765
## <none> 147.49 70.898
## - am 1 10.5467 158.04 71.108
## - wt 1 27.0144 174.51 74.280
##
## Step: AIC=68.93
## mpg ~ cyl + disp + hp + drat + wt + qsec + am + gear + carb
##
## Df Sum of Sq RSS AIC
## - carb 1 0.4593 148.11 67.032
## - gear 1 1.3966 149.05 67.234
## - drat 1 1.6590 149.31 67.290
## - disp 1 3.7977 151.45 67.745
## - hp 1 6.7849 154.44 68.370
## <none> 147.66 68.932
## - am 1 10.4691 158.12 69.125
## - qsec 1 11.2569 158.91 69.284
## - wt 1 27.5664 175.22 72.410
##
## Step: AIC=67.03
## mpg ~ cyl + disp + hp + drat + wt + qsec + am + gear
##
## Df Sum of Sq RSS AIC
## - gear 1 0.976 149.09 65.242
## - drat 1 1.379 149.49 65.328
## <none> 148.11 67.032
## - disp 1 10.523 158.64 67.228
## - am 1 10.682 158.80 67.260
## - hp 1 12.061 160.18 67.537
## - qsec 1 14.314 162.43 67.984
## - wt 1 65.940 214.05 76.816
##
## Step: AIC=65.24
## mpg ~ cyl + disp + hp + drat + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## - drat 1 1.901 150.99 63.648
## - disp 1 9.548 158.64 65.228
## <none> 149.09 65.242
## - hp 1 11.366 160.46 65.593
## - qsec 1 13.344 162.43 65.985
## - am 1 14.387 163.48 66.190
## - wt 1 65.175 214.26 74.847
##
## Step: AIC=63.65
## mpg ~ cyl + disp + hp + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## - disp 1 8.826 159.82 63.465
## <none> 150.99 63.648
## - hp 1 10.423 161.41 63.784
## - qsec 1 12.129 163.12 64.120
## - am 1 17.694 168.69 65.194
## - wt 1 65.318 216.31 73.151
##
## Step: AIC=63.47
## mpg ~ cyl + hp + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## - hp 1 7.967 167.78 63.022
## - qsec 1 10.180 170.00 63.442
## <none> 159.82 63.465
## - am 1 15.402 175.22 64.410
## - wt 1 60.723 220.54 71.771
##
## Step: AIC=63.02
## mpg ~ cyl + wt + qsec + am
##
## Df Sum of Sq RSS AIC
## <none> 167.78 63.022
## - am 1 12.820 180.60 63.378
## - qsec 1 23.262 191.05 65.177
## - wt 1 99.706 267.49 75.947
##
## Call:
## lm(formula = mpg ~ cyl + wt + qsec + am, data = mtcars)
##
## Coefficients:
## (Intercept) cyl wt qsec am
## 14.9289 -0.3542 -3.6439 1.0126 2.4679
でcylが強制投入されたモデルがでてくる。
もちろん、初期のモデルに強制投入したい独立変数が含まれていない場合はエラーがでる。
step(lm(mpg ~ wt, data = mtcars), scope = list(lower = mpg ~ cyl))
## Start: AIC=73.22
## mpg ~ wt
## Error in factor.scope(ffac, list(add = fadd, drop = fdrop)): 下位のスコープはモデル中に含まれない項 'cyl' を持ちます
ちなみに、scopeで指定するのは直接式の形でなくても、lm()やglm()オブジェクトでもよしなにしてくれる。
model_lower <- glm(mpg ~ cyl, data = mtcars) model_full <- glm(mpg ~ ., data = mtcars) step(model_full, scope = list(lower = model_lower))
## Start: AIC=163.71
## mpg ~ cyl + disp + hp + drat + wt + qsec + vs + am + gear + carb
##
## Df Deviance AIC
## - vs 1 147.66 161.75
## - carb 1 147.90 161.80
## - gear 1 148.85 162.00
## - drat 1 149.12 162.06
## - disp 1 151.41 162.55
## - hp 1 154.33 163.16
## - qsec 1 156.36 163.58
## <none> 147.49 163.71
## - am 1 158.04 163.92
## - wt 1 174.51 167.09
##
## Step: AIC=161.74
## mpg ~ cyl + disp + hp + drat + wt + qsec + am + gear + carb
##
## Df Deviance AIC
## - carb 1 148.11 159.84
## - gear 1 149.05 160.05
## - drat 1 149.31 160.10
## - disp 1 151.45 160.56
## - hp 1 154.44 161.18
## <none> 147.66 161.75
## - am 1 158.12 161.94
## - qsec 1 158.91 162.10
## - wt 1 175.22 165.22
##
## Step: AIC=159.84
## mpg ~ cyl + disp + hp + drat + wt + qsec + am + gear
##
## Df Deviance AIC
## - gear 1 149.09 158.05
## - drat 1 149.49 158.14
## <none> 148.11 159.84
## - disp 1 158.64 160.04
## - am 1 158.80 160.07
## - hp 1 160.18 160.35
## - qsec 1 162.43 160.80
## - wt 1 214.05 169.63
##
## Step: AIC=158.05
## mpg ~ cyl + disp + hp + drat + wt + qsec + am
##
## Df Deviance AIC
## - drat 1 150.99 156.46
## - disp 1 158.64 158.04
## <none> 149.09 158.05
## - hp 1 160.46 158.41
## - qsec 1 162.43 158.80
## - am 1 163.48 159.00
## - wt 1 214.26 167.66
##
## Step: AIC=156.46
## mpg ~ cyl + disp + hp + wt + qsec + am
##
## Df Deviance AIC
## - disp 1 159.82 156.28
## <none> 150.99 156.46
## - hp 1 161.41 156.60
## - qsec 1 163.12 156.93
## - am 1 168.69 158.01
## - wt 1 216.31 165.96
##
## Step: AIC=156.28
## mpg ~ cyl + hp + wt + qsec + am
##
## Df Deviance AIC
## - hp 1 167.78 155.83
## - qsec 1 170.00 156.25
## <none> 159.82 156.28
## - am 1 175.22 157.22
## - wt 1 220.54 164.58
##
## Step: AIC=155.83
## mpg ~ cyl + wt + qsec + am
##
## Df Deviance AIC
## <none> 167.78 155.83
## - am 1 180.60 156.19
## - qsec 1 191.05 157.99
## - wt 1 267.49 168.76
##
## Call: glm(formula = mpg ~ cyl + wt + qsec + am, data = mtcars)
##
## Coefficients:
## (Intercept) cyl wt qsec am
## 14.9289 -0.3542 -3.6439 1.0126 2.4679
##
## Degrees of Freedom: 31 Total (i.e. Null); 27 Residual
## Null Deviance: 1126
## Residual Deviance: 167.8 AIC: 155.8
