Rでステップワイズ回帰

ステップワイズ回帰とは説明する変数(独立変数)に何を入れれば、最も説明力が高いモデルが作れるかを自動的に考えてくれるという方法だ。日本語ではSASのJMPのページの解説がよさそうに思えた。

www.jmp.com

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