×

[PR]この広告は3ヶ月以上更新がないため表示されています。
ホームページを更新後24時間以内に表示されなくなります。

HomeRな予測モデル選択 実践編

モデル選択 実践編

★新サイト完成しました!
3秒後に自動的に移動します

変わらない方は こちらからどうぞ
http://logics-of-blue.com/%E3%83%A2%E3%83%87%E3%83%AB%E9%81%B8%E6%8A%9E_%E5%AE%9F%E8%B7%B5%E7%B7%A8/

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

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

使用データ

例に よって乱数を発生させます。
set.seed(0)
N<-100
Intercept<-5
B1<-10
B2<-5
x1<-sort(rnorm(N,sd=2))
x2<-rnorm(N,sd=2)
e<-rnorm(n=N,sd=3)
y<-Intercept+B1*x1+B2*x2+e

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

set.seed(0)
  rnorm(10)

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

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

 モデル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 

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

検定による変数選択

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

…… 実は検定によるモデル選択では、どういう順番でモデル選択をするのか(簡単なモデルから発展させるのか、複雑なモデルを簡略化するのか)というのがとても 重要になってきます。順番が変わると結果として得られるモデルが変わることもあるので。
 でもAICを使うとそんな心配は不要です。そう言う意味合いでもAICは便利ですね。
 水産関連では文献[1]に詳しく載っています。こちらでも、検定によるモデル選択の結果複数のモデルが選ばれてしまったのならば、AICによるモデル選 択をした方がよいと推奨されています。

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

summary(mdel3)

> 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

 紫 で示したところが重要。下線などは私が引きました。Rでは引いてくれません。
   星 * が 付いていれば「役に立っている変数」で、ついて無ければ使えない変数ということになります。ここでは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 は両方入れるべきということになります。

AICによる変数選択

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

> AIC(model1)
  [1] 751.6853

 > AIC(model2)
  [1] 520.6363

  > AIC(model3)
  [1] 521.6896

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

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

model.best1<-step(model3)

> model.best1

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

Coefficients:
(Intercept)           x1           x2 
      5.202        9.973        4.986 

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

MuMIn  パッケージを使って調べてみます。パッケージをインストールした後、

library(MuMIn)

 kekka.AIC<-dredge(model3,rank="AIC")


 > kekka.AIC
Global model: lm(formula = y ~ x1 * x2)
---
Model selection table
  (Int)     x1    x2  x1:x2 k   R.sq Adj.R.sq     RSS   AIC   delta weight
4 5.202  9.973 4.986        4 0.9769   0.9764   986.0 520.6   0.000  0.629
5 5.184  9.981 4.988 0.1008 5 0.9771   0.9764   976.7 521.7   1.053  0.371
2 4.734 10.260              3 0.7621   0.7596 10140.0 751.7 231.000  0.000
3 5.697        5.466        3 0.2588   0.2512 31590.0 865.3 344.700  0.000
1 5.200                     2 0.0000   0.0000 42610.0 893.3 372.600  0.000

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

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

AIC 最小モデルを引っ張ってくるには

all.model <- get.models(kekka.AIC)

best.model<-all.model[1]

> best.model
$`4`

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

Coefficients:
(Intercept)           x1           x2 
      5.202        9.973        4.986 

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

Akaike Weight による Model Averaging

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

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

> avg.model
Component models:
1+2   1+2+3

Coefficients:
(Intercept)          x1          x2       x1:x2
 5.19507788  9.97579774  4.98687373  0.03743209

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

> summary(avg.model)

Model summary:
      Deviance    AIC Delta Weight
1+2     986.00 520.64  0.00   0.63
1+2+3   976.71 521.69  1.05   0.37

Variables:
    1     2     3
   x1    x2 x1:x2

Model-averaged coefficients:
            Coefficient      SE 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
 1.00  1.00  0.37

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

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

参考文献[2]
   Package manual in PDF

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


ページトップへ
Copyright (C) 2011- 海と魚と統計解析 All Rights Reserved.