「統計学関連なんでもあり」の過去ログ--- 045

No.16843　モデルの変数選択について（全般的に）　　【波音】　2012/04/29(Sun) 17:27

重回帰分析の変数選択に関する話題です。
（以下はRにあるswissのデータを使っての例です。）

～～～変数リスト～～～
Fertility：　出生率
Examination：　軍事テストで最高点を取った被徴兵者の割合
Catholic：　カトリック教徒の割合
Infant.Mortality：　1歳未満の幼児死亡率
Agricalture：　農業従事者
Education：　被徴兵者の小学校以上の学歴


> #### フルモデルを解析
> result <- lm(Fertility ~ ., data=swiss)
Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
Examination      -0.25801    0.25388  -1.016  0.31546    
Education        -0.87094    0.18303  -4.758 2.43e-05 ***
Catholic          0.10412    0.03526   2.953  0.00519 ** 
Infant.Mortality  1.07705    0.38172   2.822  0.00734 **
Multiple R-squared: 0.7067,     Adjusted R-squared: 0.671

> #### 有意でない変数をモデルから除去
> result2 <- update(result, ~. -Examination)
Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      62.10131    9.60489   6.466 8.49e-08 ***
Agriculture      -0.15462    0.06819  -2.267  0.02857 *  
Education        -0.98026    0.14814  -6.617 5.14e-08 ***
Catholic          0.12467    0.02889   4.315 9.50e-05 ***
Infant.Mortality  1.07844    0.38187   2.824  0.00722 **
Multiple R-squared: 0.6993,     Adjusted R-squared: 0.6707

上記の例はフルモデル（result）から有意でない変数を除去（result2）するというやり方です。ちなみにstep()をデフォルトで実行してもresult2と同じモデルが最終的に得られます。

【1つ目の話題】
気になるのは，自由度調整済み決定係数（Adjusted R-squared）がresultの方が少しだけでも高いということです。私としては"最適な"モデルはresult2だと思うのですが，実際に Examinationという変数は除いても除かなくても大して影響ないのかな？と。。。

モデル作成に使ったデータにモデルを当てはめて予測値を出す際には変わらないけど，別のデータ（予測をしたいデータ）にモデルを当てはめたときにresult2の方が好ましいということでしょうか？

それにしてもExaminationのあるなしがどれだけ予測に影響するのだろうか，という疑問があります。

【2つ目の話題】
それから，この掲示板でもよく見かけることですが単変量解析の結果を寄せ集めても多変量解析の結果にはならないと思っています。だから以下のような方法は間違っていると思うのですが，どうでしょうか？

(1) Fertilityを応答変数，Agricultureを説明変数とした単回帰分析を行う
(2) 係数に対するp値を求める（そのモデルのAICを求める）
(3) これを説明変数の数だけ繰り返す。最終的に説明変数5つ分のp値とAICが求まる
(4)-1　p<0.05である変数のみを用いて重回帰分析を行う
(4)-2　AICの閾値を設定，もしくは昇順に並べて上位いくつかの変数を採用して重回帰分析を行う

Rで書くとこんな感じです↓


#### 個々の単回帰モデルについてAICを算出
my.AIC <- numeric(5)
for(i in 1:5){
	res <- lm(swiss$Fertility ~ swiss[, i+1])
	my.AIC[i] <- AIC(res)
}
data.list(matrix(my.AIC, ncol=1))
369.467485089417 
350.352503439104 
348.42229681314 
364.347891835399 
366.768286523082

#### 個々の単回帰モデルについてp値を算出
my.p <- numeric(5)
for(i in 1:5){
	res <- lm(swiss$Fertility ~ swiss[, i+1])
	my.p[i] <- summary(res)$coefficients[2, 4]
}
data.list(matrix(my.p, ncol=1))
0.0149172006147264 
9.45043734069043e-07 
3.65861696596238e-07 
0.00102852319011815 
0.00358523805395057

この方法でやると，"最適な"モデルはresultということになってしまいます。個々の単回帰モデルのp値は全て0.05より小さいので全ての変数を説明変数に入れるべきという結論になるから。

しかし具体的にこの方法では何が問題なのか，と質問されてもどのように説明すればよいのか答えに窮しています。。。（汗）

No.16844　Re: モデルの変数選択について（全般的に）　　【青木繁伸】　2012/04/29(Sun) 21:18

単純にAIC最小のモデルを選択したでではないでしょうか。
決定係数により選択するという基準からはほんのわずかの差に過ぎないのですから，それを元にしてどちらかが優れていると決定しなければならないというものでもないでしょう。
将来の予測精度についての情報もはないわけですらか，どちらが優れているという判断も出来ないでしょう。
少なくとも，得られている情報からは多変量的に考えていくつかのモデルに絞られ，そのうちでより妥当性のあるモデルはデータ以外の「実質的な状況証拠」に基づいて選択されるべきものでしょう，ということでしょう。統計学的な分析結果や実質的な状況証拠に甲乙付けがたいときには，しいて数値的な違いのみでどちらが優れているというような判断はしないほうがよいと思います。

No.16851　Re: モデルの変数選択について（全般的に）　　【波音】　2012/04/30(Mon) 11:33

回答ありがとうございます。

うーん，やはり数値的基準で判断するのが難しい場合（どちらが適切か灰色な場合）は意味的な解釈をしないといけないですね。実質的な仮説に基づいてモデルを考えなければならないというのは，なかなか難しい話ではありますが重要事項ですね。