ステップワイズ回帰とは説明する変数(独立変数)に何を入れれば、最も説明力が高いモデルが作れるかを自動的に考えてくれるという方法だ。日本語ではSASのJMPのページの解説がよさそうに思えた。
PCで統計パッケージを使って行えば、自動的に最もフィッティングの良いモデルの選択ができる、という触れ込みだが、実際には使うべきではない。仮説の無いデータ分析をすることがNGであることはウェブにもよく書いてあることだが、技術的にも、リッジ回帰、LASSO、Elastic Netといった洗練された方法があるので、現代でステップワイズ回帰をわざわざ行う理由が見当たらない。
ステップ関数を利用したステップワイズ回帰
ステップワイズ回帰はRの標準機能で走らせることができる。
まずはデータを読み込む。よく利用しているAERパッケージに含まれるCPS1985を今回も使って説明してみる。
データの読み込み。
library(AER) data(CPS1985)
CPS1985には下記の変数が含まれている。
- wage...賃金
- education...教育年数
- experience...就労年数
- age...年齢
- ethnicity...エスニシティ
- region...地域
- gender...性別
- 職業...occupation
- sector...職種
- union...組合加入
- married....結婚
賃金を説明するモデルをステップワイズ回帰で求めてみよう。
fit1 <- lm(formula = wage ~ education + experience + age + ethnicity
+ region + gender + occupation + sector + union + married,
data = CPS1985) # fit1に回帰モデルを格納する
res1 <- step(fit1) #格納した回帰モデルfit1をステップ関数で分析
summary(res1) #結果表示
回帰モデルをstep()に入れればそれでOKだ。非常に簡単。
Rだとフィッテングを判断するのはAICである。
Step: AIC=1564.76
wage ~ education + experience + region + gender + occupation + union
Df Sum of Sq RSS AIC
<none> 9599.4 1564.8
- region 1 50.52 9649.9 1565.6
- union 1 163.53 9762.9 1571.8
- gender 1 362.48 9961.9 1582.5
- occupation 5 684.51 10283.9 1591.5
- experience 1 587.41 10186.8 1594.5
- education 1 845.68 10445.1 1607.8
説明モデルとして採用されたのは、以下の変数である。
- education...教育年数
- experience...就労年数
- region...地域
- gender...性別
- occupation...職業
- union...組合加入
下の部分は変数を投入した場合のAICが表示されている。
<none>と書いてある行が最適だと判断されたモデルの統計量である。AICは1564.8となっている。一つ下にあるregion(地域)を入れたもののAICが1565.6になり、モデルのフィッティングが悪くなっていることがわかる。
MASSパッケージを用いてステップワイズ回帰
MASSパッケージを使ってもステップワイズ回帰ができる。ノーマルなステップワイズ回帰であれば、ステップ関数を使ったものと同じものが出力されるのであまり出番はないかもしれない。
library(MASS) res2 <- stepAIC(fit1, direction="both") res2$anova
結果表示の最後に$anovaをつけないと、変数の入れ替えでのAICの変化が表示されないので、MASSパッケージを使う際にはつける必要がある。
direction="both"はデフォルト値なので入れなくてもよい。bothは変数増減法を指定するオプションで変数増加法と変数減少法もあるが、特に変えなくても良いと思う。ステップワイズ回帰は総当たり方式ではないので、この部分の指定によって違ったモデルが最適だと出力される場合があり、どのような順番で変数を検討しているかは統計パッケージで異なるので、どこかに結果を出す場合には、ソフトウェアとどのような方法で変数選択を行ったかを書いた方が良い。
stepAIC {MASS}
https://stat.ethz.ch/R-manual/R-devel/RHOME/library/MASS/html/stepAIC.html
選択されモデルの交互作用項のリストワイズ回帰
MASSパッケージでは、交互作用の分析もできる。
さきほど、分析で選択された変数で新しい回帰モデルをつくる。
scope=list(upper=~X1 * X2 * ... * Xn, lower=~1)の中に回帰分析と同じ変数をいれる。アスタリスクは交互作用を意味している。
fit2 <- lm(formula = wage ~ education + experience +
region +gender + occupation + union,
data = CPS1985)
# 選択された変数で構成した回帰モデル
res3 <- stepAIC(fit2,scope=list(upper=~education * experience *
region * gender * occupation * union, lower=~1),
direction="both")
# 交互作用の検討
res3$anova # 結果の表示
結果は以下のようになる。
Initial Model:
wage ~ education + experience + region + gender + occupation +
union
Final Model:
wage ~ education + experience + region + gender + occupation +
union + experience:gender + education:region + education:experience +
education:union + education:occupation + experience:union +
education:experience:union
Step Df Deviance Resid. Df Resid. Dev AIC
1 523 9599.374 1564.757
2 + experience:gender 1 103.62041 522 9495.753 1560.961
3 + education:region 1 64.37861 521 9431.375 1559.328
4 + education:experience 1 57.80331 520 9373.571 1558.045
5 + education:union 1 64.71459 519 9308.857 1556.346
6 + education:occupation 5 204.93102 514 9103.926 1554.459
7 + experience:union 1 72.08144 513 9031.844 1552.214
8 + education:experience:union 1 106.19011 512 8925.654 1547.898
Final Modelが選択されたモデルである。交互作用項が7つ含まれている。
交互作用項を含めたVIFの出力
このようなモデルだと多重共線性も検討したい。VIFの出力はcarパッケージに入っているので、さきほどのステップワイズ回帰で出力されたモデルをそのまま入れてみよう。
library(car)
fit3 <- lm(formula =wage ~ education + experience + region +
gender + occupation + union + experience:gender +
education:region + education:experience +
education:union + education:occupation +
experience:union + education:experience:union,
data = CPS1985)
vif(fit3)
結果は以下のようになる。
GVIF Df GVIF^(1/(2*Df)) education 1.345057e+01 1 3.667502 experience 2.725167e+01 1 5.220313 region 2.880692e+01 1 5.367208 gender 3.444389e+00 1 1.855907 occupation 2.969113e+08 5 7.034977 union 1.419677e+02 1 11.915022 experience:gender 4.597440e+00 1 2.144164 education:region 3.378791e+01 1 5.812737 education:experience 2.331528e+01 1 4.828590 education:union 1.343270e+02 1 11.589952 education:occupation 3.895536e+08 5 7.228642 experience:union 9.909437e+01 1 9.954616 education:experience:union 8.924635e+01 1 9.447029
全体的に数値が高い。VIFは大きさだけでは判断できないとはいえ、このモデルはおそらくダメなモデルである。
ロジットモデルのステップワイズ回帰
例として、結婚を推定する変数を探索するロジットモデルのステップワイズ回帰のスクリプトを書いた。
SPSSでは(おそらく)重回帰のステップワイズ分析くらいしかできなかったと思うが→やってみたらできた。
fit4 <- glm(formula = married ~ wage+ education + experience + age + ethnicity
+ region + gender + occupation + sector + union, family=binomial,
data = CPS1985)
res4 <- step(fit4)
summary(res4)
MASSパッケージを用いる場合次のように書く。
res5 <- stepAIC(fit4, direction="both") res5$anova
