★新サイト完成しました!
３秒後に自動的に移動します
変わらない方は こちらからどうぞ
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の傾きが５というのが真の値。
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によるモデル選
択をした方がよいと推奨されています。
まぁ、
難しい話は置いておいて、検定によるモデル選択をしてみます。
各変数が役に立っているみなせる（「役に立っていると見えるが、それはただの誤差だ」と言える確率が５％以下）かどうかは次の計算で一発で分かります。
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 よりも複雑なので当てはまりがよいように見えるが、実はそれは誤差の範囲内なのだ」という
確率が３４％もあることが分かります。これは誤差の範囲内と見なすところです。だから、やっぱり交互作用は要らなかったと。
一
方モデル1とモデル2 の比較では、「model2 の方がmodel1 よりも複雑なので当てはまりがよいように見えるが、実はそれは誤差の範囲内なの
だ」という確率が、滅茶苦茶小さな値になっています。10の-16乗くらいの確率です。これは誤差の範囲とは言えない。
二つの変数x1とx2 は両方入れるべきということになります。
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基準での選択ということになっていますが、デフォルトではAICｃでの選択ということになっています。今回の例では正規分布を仮定しているので AICｃでもよい(ハズ)のですが、とりあえず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 関数を使った場合と一致しました。
ここ からは少し怪しげな世界に入りますが、 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
前のページ ⇒ モ
デル選択 理論編
次のページ ⇒ 重回帰分析
