今回はR言語でステップワイズ法を実行する方法を紹介します。
AIC（赤池情報量基準）に基づくステップワイズ法の実行方法や実際の解析例をまとめました。
この記事では重回帰分析や一般化線形モデルに対する変数選択の実行方法を紹介します。
重回帰分析や一般化線形モデルの実行方法については、それぞれ以下の記事で紹介しています。
-
【R言語】関数lmによる線形回帰 単回帰分析・重回帰分析 回帰分析で予測する
R言語の回帰分析の実行方法を紹介していきます。回帰分析の実行方法だけでなく、回帰分析により作ったモデルでテストデータの予測も行っていきます。 実行方法では、回帰分析の結果の見方やその抽出方法を詳しく解 ...
続きを見る
-
【R言語】一般化線形モデル 関数glmの使いかた
一般化線形モデルを行う方法を解説していきます。 R言語で回帰分析を行う方法については、以下の記事で紹介しています。 一般化線形モデルを用いることで、母分布が正規分布でなくても線形モデルを構築することが ...
続きを見る
この記事で紹介しているプログラミングコードは以下のRスクリプトに保存しています。
ポチップ
ポチップ関数step
関数stepの使用法やその引数についてまとめました。実行例については、この後で紹介しています。
direction = c("both", "backward", "forward"),
trace = 1, keep = NULL, steps = 1000, k = 2, …)
| object | 主にlmやglmによって作られるモデルを表すオブジェクト。ステップワイズ法を適用するの一番初めのモデルとして使われる。 |
| scope | ステップワイズ法において探索されるモデルの範囲を指定する。 |
| scale | モデル選択においてAIC統計量の定義で用いられるスケール。 |
| direction | "both"で変数増減法、"backward"で変数減少法、"forward" で変数増加方を指定できる。デフォルトで"both"。 |
| trace | 正（>0）の場合、stepを実行中にステップワイズの途中経過がコンソールにプリントされる。さらに大きい数の場合、より詳細な結果がプリントされる。 |
| keep | 変数選択においてモデルに保持しておきたい変数。 |
| steps | 最大ステップ数。この数を超えると実行が止まる。 |
| k | ペナルティー項の自由度の指定に用いられる。 |
関数stepの使用例
重回帰分析
まず、重回帰分析における関数stepを実行例をみていきます。データセットとしてBostonの住宅価格のデータを用いたいので、次のパッケージMASSをインストールし、ライブラリに追加します。
1 | library(MASS) |
上記を実行したら、以下のようにBostonの住宅価格のデータフレームが使用可能になります。
1 2 3 4 5 6 7 8 9 | > data <- Boston > head(data) crim zn indus chas nox rm age dis rad tax ptratio black lstat medv 1 0.00632 18 2.31 0 0.538 6.575 65.2 4.0900 1 296 15.3 396.90 4.98 24.0 2 0.02731 0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14 21.6 3 0.02729 0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03 34.7 4 0.03237 0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94 33.4 5 0.06905 0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33 36.2 6 0.02985 0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21 28.7 |
早速、住宅価格medvに関する重回帰分析を実行します。次のように、関数lmを用いてmedvをそれ以外の変数で予測します。
1 2 3 4 5 6 7 8 9 10 11 | > regression <- lm(medv ~., data = data) > regression Call: lm(formula = medv ~ ., data = data) Coefficients: (Intercept) crim zn indus chas nox rm age 3.646e+01 -1.080e-01 4.642e-02 2.056e-02 2.687e+00 -1.777e+01 3.810e+00 6.922e-04 dis rad tax ptratio black lstat -1.476e+00 3.060e-01 -1.233e-02 -9.527e-01 9.312e-03 -5.248e-01 |
上記のようにregressionには重回帰分析の結果が入っており、切片パラメータを合わせると回帰係数が14個もあることが分かります。
重回帰分析の問題点
重回帰分析だけの場合、これらの回帰係数に関する検定\(H_0:\ \beta_i = 0,\ i = 0,\ldots ,13\)を行うと、合計14回も仮説検定を行うことになり、検定の多重性が問題になってきます。
これを回避するために、次のように関数stepを用いて変数選択を実行します。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 | > step.reg <- step(regression) Start: AIC=1589.64 medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + tax + ptratio + black + lstat Df Sum of Sq RSS AIC - age 1 0.06 11079 1587.7 - indus 1 2.52 11081 1587.8 <none> 11079 1589.6 - chas 1 218.97 11298 1597.5 - tax 1 242.26 11321 1598.6 - crim 1 243.22 11322 1598.6 - zn 1 257.49 11336 1599.3 - black 1 270.63 11349 1599.8 - rad 1 479.15 11558 1609.1 - nox 1 487.16 11566 1609.4 - ptratio 1 1194.23 12273 1639.4 - dis 1 1232.41 12311 1641.0 - rm 1 1871.32 12950 1666.6 - lstat 1 2410.84 13490 1687.3 Step: AIC=1587.65 medv ~ crim + zn + indus + chas + nox + rm + dis + rad + tax + ptratio + black + lstat Df Sum of Sq RSS AIC - indus 1 2.52 11081 1585.8 <none> 11079 1587.7 - chas 1 219.91 11299 1595.6 - tax 1 242.24 11321 1596.6 - crim 1 243.20 11322 1596.6 - zn 1 260.32 11339 1597.4 - black 1 272.26 11351 1597.9 - rad 1 481.09 11560 1607.2 - nox 1 520.87 11600 1608.9 - ptratio 1 1200.23 12279 1637.7 - dis 1 1352.26 12431 1643.9 - rm 1 1959.55 13038 1668.0 - lstat 1 2718.88 13798 1696.7 Step: AIC=1585.76 medv ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat Df Sum of Sq RSS AIC <none> 11081 1585.8 - chas 1 227.21 11309 1594.0 - crim 1 245.37 11327 1594.8 - zn 1 257.82 11339 1595.4 - black 1 270.82 11352 1596.0 - tax 1 273.62 11355 1596.1 - rad 1 500.92 11582 1606.1 - nox 1 541.91 11623 1607.9 - ptratio 1 1206.45 12288 1636.0 - dis 1 1448.94 12530 1645.9 - rm 1 1963.66 13045 1666.3 - lstat 1 2723.48 13805 1695.0 |
実行するとコンソール上に上の結果が出力されることが分かります。長くなって分かりにくいですが、Step: AIC = . . .の部分を見るとわかるように2ステップ目のモデルのAICが最小であることが分かります。
次のようにstep.regを関数summaryによって変換し、summaryStep.regを参照すると上でAICが最小だったモデルが選択されていることが確認できます。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 | > summaryStep.reg <- summary(step.reg) > summaryStep.reg Call: lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat, data = data) Residuals: Min 1Q Median 3Q Max -15.5984 -2.7386 -0.5046 1.7273 26.2373 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 36.341145 5.067492 7.171 2.73e-12 *** crim -0.108413 0.032779 -3.307 0.001010 ** zn 0.045845 0.013523 3.390 0.000754 *** chas 2.718716 0.854240 3.183 0.001551 ** nox -17.376023 3.535243 -4.915 1.21e-06 *** rm 3.801579 0.406316 9.356 < 2e-16 *** dis -1.492711 0.185731 -8.037 6.84e-15 *** rad 0.299608 0.063402 4.726 3.00e-06 *** tax -0.011778 0.003372 -3.493 0.000521 *** ptratio -0.946525 0.129066 -7.334 9.24e-13 *** black 0.009291 0.002674 3.475 0.000557 *** lstat -0.522553 0.047424 -11.019 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 4.736 on 494 degrees of freedom Multiple R-squared: 0.7406, Adjusted R-squared: 0.7348 F-statistic: 128.2 on 11 and 494 DF, p-value: < 2.2e-16 |
また、重回帰分析で関数lmから作られたオブジェクトと同様に、次のように回帰係数や決定係数など結果はリストの要素として格納されています。
1 2 3 | coef <- summaryStep.reg$coefficients #回帰係数 r.squared <- summaryStep.reg$r.squared #決定係数 adj.r.squared <- summaryStep.reg$adj.r.squared #修正決定係数 |
以下のようにして、関数cbindを用いて回帰係数、決定係数、修正決定係数をデータフレームに格納しresultTableとし、関数write.csvを用いることで、ステップワイズ法で選択されたモデルの結果をcsvファイルへと出力することができます。
1 2 | resultTable <- cbind(coef, r.squared = r.squared, adj.r.squared = adj.r.squared) write.csv(resultTable, "ステップワイズ法（重回帰分析）.csv") |
作業ディレクトリを確認すると、次の画像の結果が保存されていることが分かります。
一般化線形モデル
一般化線形モデルに対しても重回帰分析と同様に、ステップワイズ法により変数選択ができます。データセットとしDefaultを用いたいので次のパッケージLSLRをインストールしライブラリに追加します。
1 | library(ISLR) |
Defaultは以下の4つのデータdefault、student、balance、incomeから成ります。
1 2 3 4 5 6 7 8 9 | > data <- Default > head(data) default student balance income 1 No No 729.5265 44361.625 2 No Yes 817.1804 12106.135 3 No No 1073.5492 31767.139 4 No No 529.2506 35704.494 5 No No 785.6559 38463.496 6 No Yes 919.5885 7491.559 |
一般化線形モデルで紹介したように、1列目のdefaultについてロジスティック回帰分析を行っていきます。まず準備として、defaultの水準"Yes"、"No"をそれぞれ1、0に変換しnumeric型に変換してあげます。
1 2 3 4 | data$default <- as.character(data$default) data$default <- replace(data$default, data$default == "Yes", 1) data$default <- replace(data$default, data$default == "No", 0) data$default <- as.numeric(data$default) |
上を実行し終わったら、次のように関数glmを用いてdefaultについてのロジスティック回帰分析を行います。2行目を見ればわかるように、lmのときと同じようにして変数選択できます。
1 2 3 | regression <- glm(default ~., data = data, family = binomial(logit)) step.reg <- step(regression) summaryStep.reg <- summary(step.reg) |
関数summaryをglmのオブジェクトに適用すると、ステップワイズ法で選ばれたモデルの回帰係数や決定係数などの結果を得ることができます。
1 2 3 | coef <- summaryStep.reg$coefficients r.squared <- summaryStep.reg$r.squared adj.r.squared <- summaryStep.reg$adj.r.squared |
これらの結果をcsvファイルへ出力したい場合は、以下のようにデータフレームに結果を保存しwrite.csvでcsvファイルへと出力します。
1 2 | resultTable <- cbind(coef, r.squared = r.squared, adj.r.squared = adj.r.squared) write.csv(resultTable, "ステップワイズ法（ロジスティック回帰分析）.csv") |
ステップワイズ法の途中結果を得たい場合
最後に、ステップワイズ中の途中経過を取得する方法について解説します。
関数stepの引数traceについて
関数stepには引数としてtraceが用意されていますが、どれだけtraceの値を大きくしてもコンソールに表示される途中結果の詳細が増えるだけで、stepによって作られるオブジェクトが途中結果に関する回帰係数などの結果を持つようになるわけではありません。
ステップワイズ中の途中経過を得る方法を探しても具体的な例が無かったため、ここにまとめておきます。
データセットとして、再びつぎのBostonの住宅価格のデータを用います。重回帰分析のときと同様に以下のようにして重回帰分析を実行します
1 2 3 | data <- Boston regression <- lm(medv ~., data = data) |
ステップワイズの途中結果を得るために次のように、関数capture.outputを用いてコンソール上にプリントされた結果を文字列のベクトルとして保存します。
1 | stepResults <- capture.output(step(regression)) |
次にstepResultsには各ステップのモデルの式（\(Y \sim X_1 + \cdots + X_k\)）が保存されているため、各式に対応するインデックスをstart、endとして抽出します。
1 2 3 4 5 | indicesEachResultStarts <- grep("AIC=", stepResults) indicesEachResultEnds <- stepResults == "" start <- indicesEachResultStarts + 1 indicesEnds <- seq_len(length(stepResults))[indicesEachResultEnds] end <- indicesEnds[seq(1, length(indicesEnds), 2)] - 1 |
最後に、for文を用いて各ステップのモデル式に関する重回帰分析を繰り返し実行することですべてのステップの結果を保存することが可能です。
1 2 3 4 5 6 7 8 9 10 11 | for (i in seq_len(length(start))) { eachFormula <- as.formula(paste(stepResults[start[i] : end[i]], collapse = "")) eachReg <- lm(eachFormula, data = data) eachSummaryReg <- summary(eachReg) coef <- eachSummaryReg$coefficients r.squared <- eachSummaryReg$r.squared adj.r.squared <- eachSummaryReg$adj.r.squared eachResultTable <- cbind(coef, r.squared = r.squared, adj.r.squared= adj.r.squared) write.table(eachResultTable, "ステップワイズ途中経過.csv", append = TRUE, quote = FALSE, sep = ",", row.names = TRUE, col.names = NA) write.table("", "ステップワイズ途中経過.csv", append = TRUE, quote = FALSE, sep = ",", row.names = FALSE, col.names = FALSE) } |
作業ディレクトリを確認すると次の画像のような各ステップの重回帰分析の結果が保存されていることが分かります。
まとめ
R言語でステップワイズ法を実行する方法を紹介しました。
重回帰分析や一般化線形モデルで紹介した関数lm、glmで作られたオブジェクトに対し、関数stepを適用することで簡単にステップワイズ法による変数選択が可能です。
