S-PLUS mini course 第31回

S-PLUS mini course の31回目をお届けします。

さて、前回は関数 step を使った変数選択についてご説明しました。
今回はその step の細かい点について、触れたいと思います。

・step は総当りではない
前回もちょっと触れましたが、総当りをして、AIC をすべて比較するということは行いません。
例えば、前回のメールでのコマンド

options(contrasts=c("contr.treatment", "contr.poly"))
# ダミー行列を分かりやすいように、通常の1,0 にかえる
fitall <- lm(Fuel ~(Weight + Disp. + Type)^2, data=fuel.frame)
# Weight, Disp. Type(因子）とその2次の交互作用項までを含む回帰を行う

であれば、可能なモデルは
主効果のみ：
変数１個　3モデル
変数２個 3*2/2=3モデル
変数３個 1モデル、合計７モデルが考えられる

主効果と交互作用：
上の全モデル（７）に、２次の交互作用の組み合わせ（多分、７通り）をかけた、最大４９モデルについて、回帰を行い、AIC を求めるということをしなければなりません。変数が３でこの程度ですから、変数が増えると大変面倒なことになってしまいます。

そこで、S-PLUS は、次にどの変数を減らしたらいいか、増やしたらいいか、を統計量 Cp を用いて調べ、その変数を増やした（あるいは減らした）ものと、元のモデルのAICを比べて、その変数の増減によって新しいモデルのAICが大きくなるようなら、その１つ前のモデルでストップします。

では、

fits <- step(fitall)

では、何をしているでしょうか？これは、まず、２次のフルモデルなので、変数増減が許されていても、減らす方しか、考えられません。そのとき、落とすことができるのは、この場合ならば、２次の交互作用項のみです（主効果から落ちて、交互作用だけが残る変数は想定していないためです）。
そのため、まず、Weight:Type（交互作用）を落とす、ということが試みされ元のモデルより AIC が小さくなりますから、さらにモデル増減を進めます。

次に落とすことができるのは、残りの２つの交互作用項で、Weight:Disp.を落とす、ということがやはり試されます。
すると、また AIC が小さくなりますから、さらに進めます。
今度は、可能性として、残った交互作用 Type:Disp. と、先ほどの変数選択で交互作用が全部落ちたWeight の主効果も候補になります。
そして、主効果 Weight を落としたらよいという結論に落ち着き、続いてnone と Disp.:Type を比較すると、none のほうがよいと推測されるので、もう変数を落とすことをしません。つまり、ここで収束します。

デフォルトで変数増減法（direction="both") なのに、どうして、変数を増やすほうがされていないか、この場合ならば、「主効果のない、交互作用項が考えられていない」ためです。お手元のデータで変数を増やして試していただけると、多分、増減が発生していると思います。

次回はこの変数選択の続きと、他のモデル（樹形モデル）を使った変数選択についてご紹介したいと思います。
お楽しみに！

お問い合わせ／資料請求はお気軽に。