最終更新:2016年1月24日
※フォントや参考文献を修正しました。

前のページで色々と理屈を並べたてましたが、理屈を知っていても実際に扱えないと意味がありません。

ここでは実際にモデル選択をしてみます。

ここで用いたRコードは、まとめてこちらから見ることができます。
コードは2015年8月30日に動作確認をしました。動かないものがあれば、ご一報いただければ幸いです。



スポンサードリンク

 

目次

1.使用データ
2.重回帰分析によるモデリング
3.検定による変数選択
4.AICによる変数選択
5.Akaike WeightによるModel Averaging

 

1.使用データ

例によって乱数を発生させます。

# シミュレーションデータの作成
set.seed(0)
N <- 100                        # サンプルサイズ
Intercept <- 5                  # 切片
B1 <- 10                        # 係数1
B2 <- 5                         # 係数2
x1 <- sort(rnorm(N, sd=2))      # 説明変数1
x2 <- rnorm(N, sd=2)            # 説明変数2
e <- rnorm(n=N, sd=3)           # 誤差

y <- Intercept + B1*x1 + B2*x2 + e

切片は5。x1の傾き が10、x2の傾きが5というのが真の値。
set.seed(0) とは、発生された乱数を定める関数です。
確認します

> # 普通に実行すると、毎回値が変わる
> rnorm(10)
 [1]  0.7635935 -0.7990092 -1.1476570 -0.2894616 -0.2992151 -0.4115108  0.2522234
 [8] -0.8919211  0.4356833 -1.2375384
> rnorm(10)
 [1] -0.22426789  0.37739565  0.13333636  0.80418951 -0.05710677  0.50360797
 [7]  1.08576936 -0.69095384 -1.28459935  0.04672617
> rnorm(10)
 [1] -0.2357066 -0.5428883 -0.4333103 -0.6494716  0.7267507  1.1519118  0.9921604
 [8] -0.4295131  1.2383041 -0.2793463
> 
> # set.seed(0)を設定すると値が固定される
> set.seed(0)
> rnorm(10)
 [1]  1.262954285 -0.326233361  1.329799263  1.272429321  0.414641434 -1.539950042
 [7] -0.928567035 -0.294720447 -0.005767173  2.404653389
> set.seed(0)
> rnorm(10)
 [1]  1.262954285 -0.326233361  1.329799263  1.272429321  0.414641434 -1.539950042
 [7] -0.928567035 -0.294720447 -0.005767173  2.404653389
> set.seed(0)
> rnorm(10)
 [1]  1.262954285 -0.326233361  1.329799263  1.272429321  0.414641434 -1.539950042
 [7] -0.928567035 -0.294720447 -0.005767173  2.404653389

何回やっても、全く同じバラバラの乱数 rnorm が発生されることになります。同じ値が出たほうが結果が比較しやすくて便利だと思ったので。

 

2.重回帰分析によるモデリング

 モデル1 x1だけを使ってyを計算

> model1 <- lm(y ~ x1) 
> model1

Call:
lm(formula = y ~ x1)

Coefficients:
(Intercept)           x1  
      4.734       10.260  

 

モデル2 x1とx2の両方を使ってモデリング

> model2 <- lm(y ~ x1 + x2) 
> model2

Call:
lm(formula = y ~ x1 + x2)

Coefficients:
(Intercept)           x1           x2  
      5.202        9.973        4.986  

 

モデル3 x1とx2の交互作用も入れてモデリング

> model3 <- lm(y ~ x1*x2) 
> model3

Call:
lm(formula = y ~ x1 * x2)

Coefficients:
(Intercept)           x1           x2        x1:x2  
     5.1839       9.9807       4.9876       0.1008  

 交互作用とは、変数単体が影響を与えているのではなく、変数同士が相乗効果を起こしていると仮定することです。
 たとえば、植物を育てる時、真っ暗闇で育てていたら、どれだけ肥料をあげても無駄かもしれません。しかし、日光という別の要因が十分に存在すると言う条 件下では、肥料はとても大きな効果を持つでしょう。こういった要因同士の相乗効果を扱うこともできます。
 が、今回発生させた乱数はそんな相乗効果を入れ込んだ物ではありません。交互作用を入れた複雑なモデルは「誤ったモデル」ということになります。

 

3.検定による変数選択

 変数選択のやり方は、とりあえずもっとも複雑なモデルを作って、それをどんどんと単純にしていくという流れを取ります。

…… 実は検定によるモデル選択では、どういう順番でモデル選択をするのか(簡単なモデルから発展させるのか、複雑なモデルを簡略化するのか)というのがとても重要になってきます。順番が変わると結果として得られるモデルが変わることもあるので。
 でもAICを使うとそんな心配は不要です。そう言う意味合いでもAICは便利ですね。
 水産関連では文献[1]に詳しく載っています。こちらでも、検定によるモデル選択の結果複数のモデルが選ばれてしまったのならば、AICによるモデル選択をした方がよいと推奨されています。
 また、文献[1]にはAICだけでなくAICcなどの別の情報量基準を使った時のシミュレーション結果などが載っています。高度な内容なのでここには載せませんが、サンプルサイズが小さければAICcのほうが良いといった成果がたくさん書かれているので一度読まれることを強くお勧めします。

まぁ、 難しい話は置いておいて、検定によるモデル選択をしてみます。
各変数が役に立っているみなせる(「役に立っていると見えるが、それはただの誤差だ」と言える確率が5%以下)かどうかは次の計算で一発で分かります。

> # summary関数による検定
> summary(model3)

Call:
lm(formula = y ~ x1 * x2)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.0011 -2.0218 -0.0245  1.8009  7.2887 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   5.1839     0.3200  16.200   <2e-16 ***
x1            9.9807     0.1820  54.829   <2e-16 ***
x2            4.9876     0.1663  30.001   <2e-16 ***
x1:x2         0.1008     0.1055   0.956    0.342    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.19 on 96 degrees of freedom
Multiple R-squared:  0.9771,    Adjusted R-squared:  0.9764 
F-statistic:  1364 on 3 and 96 DF,  p-value: < 2.2e-16

 星 * が 付いていれば「役に立っている変数」で、ついて無ければ使えない変数ということになります。ここではx1:x2が交互作用を表しているのですが、これだけ役に立っていない(星が一つも付いていない)ことが分かります。なので、交互作用抜きの model2 がもっとも良さそうということに。

確認します。

> anova(model2, model3)      # モデル2とモデル3の比較
Analysis of Variance Table

Model 1: y ~ x1 + x2
Model 2: y ~ x1 * x2
  Res.Df    RSS Df Sum of Sq      F Pr(>F)
1     97 986.00                           
2     96 976.71  1    9.2907 0.9132 0.3417
> 
> 
> anova(model1, model2)      # モデル1とモデル2の比較
Analysis of Variance Table

Model 1: y ~ x1
Model 2: y ~ x1 + x2
  Res.Df   RSS Df Sum of Sq      F    Pr(>F)    
1     98 10139                                  
2     97   986  1      9153 900.45 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

  俗に言う分散分析表が出力されます。model2 model3 比較では、「model3 の方がmodel2 よりも複雑なので当てはまりがよいように見えるが、実はそれは誤差の範囲内なのだ」という 確率が34%もあることが分かります。これは誤差の範囲内と見なすところです。だから、やっぱり交互作用は要らなかったと。

一 方モデル1とモデル2 の比較では、「model2 の方がmodel1 よりも複雑なので当てはまりがよいように見えるが、実はそれは誤差の範囲内なのだ」という確率が、滅茶苦茶小さな値になっています。10の-16乗くらいの確率です。これは誤差の範囲とは言えない。
 二つの変数x1とx2 は両方入れるべきということになります。

4.AICによる変数選択

AIC を全部求めてみます。 計算方法は AIC(モデル名) です

> # AICによるモデル選択
> AIC(model1)
[1] 751.6853
> 
> AIC(model2)
[1] 520.6363
> 
> AIC(model3)
[1] 521.6896

model2 がもっともAICが小さく良いモデルということになります。

次は  step 関数を使ってやってみます。

> # step関数によるモデル選択
> model.best1 <- step(model3)
Start:  AIC=235.9
y ~ x1 * x2

        Df Sum of Sq    RSS    AIC
- x1:x2  1    9.2907 986.00 234.85
<none>               976.71 235.90

Step:  AIC=234.85
y ~ x1 + x2

       Df Sum of Sq   RSS    AIC
<none>                986 234.85
- x2    1      9153 10139 465.90
- x1    1     30599 31585 579.53
> model.best1

Call:
lm(formula = y ~ x1 + x2)

Coefficients:
(Intercept)           x1           x2  
      5.202        9.973        4.986  

検定と同じく、AIC最小モデルも「交互作用なし、変数は二つとも使うモデル」ということになりました。

MuMIn  パッケージを使って調べてみます。パッケージをインストールした後に下記を実行してください。

> # MuMInパッケージを用いた総当たり法によるモデル選択
> library(MuMIn)
> options(na.action = "na.fail")
> kekka.AIC <- dredge(model3, rank="AIC")
 Fixed term is "(Intercept)" 
> kekka.AIC
Global model call: lm(formula = y ~ x1 * x2)
---
Model selection table 
  (Int)     x1    x2  x1:x2 df   logLik   AIC  delta weight
4 5.202  9.973 4.986         4 -256.318 520.6   0.00  0.629
8 5.184  9.981 4.988 0.1008  5 -255.845 521.7   1.05  0.371
2 4.734 10.260               3 -372.843 751.7 231.05  0.000
3 5.697        5.466         3 -429.658 865.3 344.68  0.000
1 5.200                      2 -444.631 893.3 372.63  0.000
Models ranked by AIC(x) 

これも、一番上にあるAIC最小モデルは「x1 x2の両方の変数を用い、交互作用はないモデル」ということになります。AICが2番目、3番目に良いモデルも一覧として見れるので便利です。

ちなみに、rank="AIC"としているので AIC基準での選択ということになっていますが、デフォルトではAICcでの選択ということになっています。今回の例では正規分布を仮定しているので AICcでもよい(ハズ)のですが、とりあえずAICで計算してみました。

AIC 最小モデルを引っ張ってくるには下記のようにします。

> # 最適なモデルを一つ選ぶ
> best.model <- get.models(kekka.AIC, subset = 1)[1]
> best.model

Call:
lm(formula = y ~ x1 + x2 + 1)

Coefficients:
(Intercept)           x1           x2  
      5.202        9.973        4.986  

検定やstep関数を使った場合と一致しました。
 

5.Akaike Weight による Model Averaging

ここからは少し怪しげな世界に入りますが、Model Averagingをしてみます。

> # Model Averaging
> avg.model <- model.avg(get.models(kekka.AIC,  subset = delta < 4)) 
> avg.model

Call:
model.avg.default(object = get.models(kekka.AIC, subset = delta < 
    4))

Component models: 
‘12’  ‘123’

Coefficients: 
       (Intercept)       x1       x2      x1:x2
subset    5.195078 9.975798 4.986874 0.10081304
full      5.195078 9.975798 4.986874 0.03743209

本来ならば使わない交互作用も組み込まれたモデルになっていますが、モデルを平均した値を使っているので、交互作用のもつウエイトは大分と小さくなっています。summary(avg.model)  で詳しく見れます。

> summary(avg.model)

Call:
model.avg.default(object = get.models(kekka.AIC, subset = delta < 
    4))

Component model call: 
lm(formula = y ~ <2 unique rhs>)

Component models: 
    df  logLik    AIC delta weight
12   4 -256.32 520.64  0.00   0.63
123  5 -255.84 521.69  1.05   0.37

Term codes: 
   x1    x2 x1:x2 
    1     2     3 

Model-averaged coefficients: 
            Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)   5.1951     0.3197      0.3237  16.048   <2e-16 ***
x1            9.9758     0.1819      0.1842  54.154   <2e-16 ***
x2            4.9869     0.1662      0.1683  29.629   <2e-16 *** 
x1:x2         0.1008     0.1055      0.1068  0.944     0.345 

Full model-averaged coefficients (with shrinkage):
       Estimate Std. Error Adjusted SE z value Pr(>|z|)    
(Intercept)  5.19508    0.31968     0.32373   16.05   <2e-16 ***
x1           9.97580    0.18191     0.18421   54.15   <2e-16 ***
x2           4.98687    0.16620     0.16831   29.63   <2e-16 ***
x1:x2        0.03743    0.08065     0.08131    0.46    0.645    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Relative variable importance: 
                     x1   x2   x1:x2
Importance:          1.00 1.00 0.37 
N containing models:    2    2    1 

Relative variable importance:を見ると、多少変な変数を入れ込んでしまっても、その悪影響は小さくなりそうだということが分かります。モデルそのものが持つ不確実性を考慮してくれているのですね。ただし、この計算によって予測精度が上がると言う保証はありません。

 

参考文献

参考文献[1]
 庄野 宏 (2006): モデル選択手法の水産資源解析への応用―情報量規準とステップワイズ検定の取り扱い― 計量生物学, 27(1), p.55-67

参考文献[2]
   Package manual in PDF

 

モデル選択に関連するその他の本

データ解析のための統計モデリング入門――
一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)


 
大変読みやすく、内容の濃い統計モデルの入門書です。通称はみどり本。
AICについての解説や、検定との比較について載っています。
 

平均・分散から始める一般化線形モデル入門


 
手前みそではございますが、管理人の書いた統計モデルの入門書です。
モデル選択として検定やstep関数、MuMInのより詳しい使い方も載せています。
 

前のページ ⇒ モデル選択_理論編    
次のページ ⇒ 重回帰分析 



スポンサードリンク