AICを使った変数選択

radvance

パッケージMASSとVR †

　以下ではstepの使い方について解説しますが、これはMASSパッケージに含まれるstepAICという関数を簡略化してstatsパッケージに再録したものです。statsパッケージはRの基本パッケージ(base)に含まれているので、MASSをインストールしなくてもstepは使うことができます。

　MASSパッケージは単体では配布されていないようです。VRという名前のパッケージに含まれているので、stepAICを使う場合はインストールする必要があります。並河が2008年以降に講義で配ったR-2.8.0には最初からインストールしてあります。

変数選択の考え方 †

参考:統計ソフトRのブログ「ステップワイズ法による変数選択」

　AICがどういうものであるかは赤池の情報量基準で簡単に触れましたが、説明というほどの説明になっていないので大学院レベルの計量経済学教科書を参照してください。学部学生の皆さんは、こうした「内部的に何をやっているのかよくわからない」関数を使わずに、何種類もlmやglmを実行して出力結果を目で読むことを薦めます。10や20のlmやglmを少しずつ変数の組み合わせを変えて実行することは、コピー&ペーストを使えば、レポートなどを書く手間に比べるとたいしたものではありませんし、そのほうが勉強になりますし、何よりもstep関数を覚えなくて済むからです。

　説明変数に使えそうな変数が10種類あるとき、採用する説明変数の組み合わせは「どれも採用しない」ことを含めて2の10乗=1024通りあります。これくらいになると、step関数を利用する値打ちがあるでしょうね。

　n種類の説明変数候補があるとき、一番「AICで測った説明力の高い」説明変数の組み合わせを選ぶためには、3つの考え方があります。以下では少し説明を簡略化していますが、stepやstepAICではもう少し細かい指定もできます。

• (1)変数増加法
  • まずn種類の説明変数候補それぞれを使って、説明変数ひとつだけの単回帰式を推定します。一番AICの高かった変数は採用です。次に残った(n-1)種類の変数と、最初に採用した変数で、説明変数ふたつの重回帰を(n-1)種類試します。
  • いちばんAICの高かった式と、最初に選んだ単回帰式を比べます。ふたつ目の説明変数を入れてAICがかえって下がってしまったなら、最初に採用した変数ひとつの単回帰式を選んで終了です。
  • ふたつ目の説明変数を入れてAICが上がったなら、ふたつ目の変数も採用です。こうして順々に採用する変数を増やしていき、どれを追加してもAICが下がるならそこでやめます。最後まで行ってしまったら、全部採用で終了です。

• (2)変数減少法
  • 変数増加法の逆です。最初にすべての説明変数候補を使った式を推定し、AICが上がるようにひとつずつ変数を抜いて行きます。どれを抜いてもAICが下がるなら、そこでストップします。仮に最後まで行くと、定数項だけの式になります。

• (3)変数増減法
  • 変数をひとつずつ増やす場合、減らす場合の両方を順に試して行って、減らしても増やしてもAICが下がるなら、そこで止めます。

変数選択スクリプト †

サンプルデータ †

まずサンプルデータファイルstepaic.csvを作ります。このサイトのほとんどのサンプルデータと同様に、RguiのあるR-2.8.0\binの下にdというフォルダをつくって、そこに置くことにします。このスクリプトとデータは2009年5月に追加したので、それ以前に並河が配布したパッケージには含まれていません。

stepaic.csvの内容はこちらのリンク先にあります。

10種類の説明変数候補と被説明変数yは、まず0〜1の一様乱数を発生させ、そのあとx9とx10にyを足して作っています。このことによって、x9、x10とyに弱い相関が生まれます。

変数同士の相関 †

サンプルスクリプトは次の通りです。。

dset <- read.csv("d/stepaic.csv", header=TRUE)
attach (dset)
cor(dset)

eq1 = lm(y~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10)
eq2 = lm(y~1)

step(eq1)
step(eq2,scope=list(upper=~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10, lower=~1))

　なるべく短い説明でstepが使えるようになることを目指して、オプションは最低限にしています。まずcor(dset)で、x9とx10にyとの相関が生まれたかどうかを見ましょう。ところで相関って何だったっけという人はこちら。

　相関行列が丸ごと出てきますがyの列だけ見ましょう。

              y
x1   0.33775521
x2   0.04382216
x3   0.06421360
x4   0.02743331
x5   0.10535508
x6  -0.25231465
x7   0.18327383
x8  -0.06560688
x9   0.34240917
x10  0.67764679
y    1.00000000

　x9とx10のほか、x1とx6もyとの相関係数が(偶然に)高いのがわかります。ところが結論を先取りすると、x1は説明変数として選ばれません。たぶんそれは、x1とx9・x10の相関が0.29〜0.36と高いのに対し、x6とx9・x10の相関は0.17〜0.25(絶対値)と少し低いからです。x9とx10が先に選ばれていると、そこにx1を加えても式全体の説明力が上がらないのです。

stepの実際 †

　まず普通の回帰分析と同様に、lm()を実行して結果に名前をつけます。eq1は説明変数候補を全部並べた回帰式、eq2は定数項だけの式(y=b0またはy=a)です。

　scopeとdirectionを指定しないで

step(eq1)

　とだけ指定すると、変数減少法を指示したことになります。だからeq1にありったけの説明変数候補を詰め込んでおく必要があります。途中経過が次のように出力されてきます。

　何もしないときが「<none>」です。変数を外すとAICの絶対値が上がるとき、いちばんAICが改善するよう変数をひとつはずします。x2、x1、x3、x7、x4、x8、x5と7つの変数がはずれたところで、残りのどれを外してもAICが下がってしまうので作業終了。x6、x9、x10でyを説明するモデルが選ばれました。

step(eq2,scope=list(upper=~x1+x2+x3+x4+x5+x6+x7+x8+x9+x10, lower=~1))

　scopeを指定してdirectionを指定しないと、変数増減法を指定したことになります。upperには説明変数候補全部を含む式、lowerには定数項だけの式を指定すると、この範囲で変数を加えたり抜いたりして試します。出発点はeq2に指定した式です。ここでは定数項だけの、説明変数がない式。途中経過は次のようになります。

　よーく見ると、変数減少法では一度外した変数のことは考えなかったのに、今度は外した変数を加えてみたらどうなるか、ステップが進むたびに試しています。これが変数増減法の中身です。でも結論は同じです。