重回帰と変数選択
重回帰とその変数選択に関するメモ。
RにはF値による検定でのステップワイズ変数選択法はないのでしょうか？
AICによる変数選択法である、step() 関数などは見当たるが検定による変数選択は見当たらない。Web検索でもヒットしない。
AICによる変数選択と検定による変数選択、どちらが良いとかそういうのは置いておいて、検定で行った時の変数の取り込みに関する基本的性質に関してメモしておきたい。
AICは非常に興味深い指標なので、そのうち詳しく勉強できたらいいと思う。
- 作者: 樺島祥介,北川源四郎,甘利俊一,赤池弘次,下平英寿,土谷隆,室田一雄
- 出版社/メーカー: 共立出版
- 発売日: 2007/07/06
- メディア: 単行本
- 購入: 4人 クリック: 74回
- この商品を含むブログ (12件) を見る
統計の専門家でなければ、AICの原理や思想を理解している人はいないと思う。実務家であれば、分かりやすい検定によって変数を選択している場合も多いと思う。
では、検定を用いた場合の変数選択法（増加法）における変数が選ばれる順序を見ていく。
疑似データを作成。
y <- rnorm(n=50,mean=0,sd=5)
x1 <- y + rnorm(n=50,mean=0,sd=3)
x2 <- x1 + rnorm(n=50,mean=0,sd=3)
x3 <- y + rnorm(n=50,mean=0,sd=10)
標準化しておく。
y <- scale(y)
x1 <- scale(x1)
x2 <- scale(x2)
x3 <- scale(x3)
相関を見る。
X <- data.frame(y,x1,x2,x3)
cor(X)
y x1 x2 x3
y 1.0000000 0.8730817 0.7581724 0.4387364
x1 0.8730817 1.0000000 0.8529398 0.3161180
x2 0.7581724 0.8529398 1.0000000 0.2672547
x3 0.4387364 0.3161180 0.2672547 1.0000000
yとx1の相関が一番高く、次にyと相関が高いのはx2、ただしx1とx2との相関も高い。 x3はyとの相関が一番低い、そしてx1、x2との相関も低くなっている。
x1、x2、x3個別にみると、x1 → x2 → x3の順番でyに対して影響（相関）が大きい。
しかし、重回帰モデルにすべて取り込んだ場合、その影響度の順序は異なる。
結論から言うと、x1 → x3 → x2の順番となる。
full <- lm(y~x1+x2+x3)
summary(full)
Call:
lm(formula = y ~ x1 + x2 + x3)
Residuals:
Min 1Q Median 3Q Max
-1.25633 -0.24685 0.01037 0.28646 1.03247
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.456e-17 6.650e-02 0.000 1.000
x1 7.723e-01 1.307e-01 5.908 3.97e-07 ***
x2 5.107e-02 1.287e-01 0.397 0.693
x3 1.809e-01 7.081e-02 2.555 0.014 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4702 on 46 degrees of freedom
Multiple R-squared: 0.7924, Adjusted R-squared: 0.7789
F-statistic: 58.53 on 3 and 46 DF, p-value: 9.728e-16
x2の係数は小さい。（p値が大きい）
個別に結果をみると、以下のようにx1 → x2 → x3の順番でyに対して影響度が高い。
y_x1 <- lm(y~x1)
summary(y_x1)
Call:
lm(formula = y ~ x1)
Residuals:
Min 1Q Median 3Q Max
-1.30111 -0.31311 0.01479 0.33869 1.08676
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.887e-17 6.967e-02 0.00 1
x1 8.731e-01 7.038e-02 12.41 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4926 on 48 degrees of freedom
Multiple R-squared: 0.7623, Adjusted R-squared: 0.7573
F-statistic: 153.9 on 1 and 48 DF, p-value: 2.2e-16
y_x2 <- lm(y~x2)
summary(y_x2)
Call:
lm(formula = y ~ x2)
Residuals:
Min 1Q Median 3Q Max
-1.8576 -0.4300 0.0155 0.4407 1.3690
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.831e-17 9.317e-02 0.000 1
x2 7.582e-01 9.412e-02 8.056 1.81e-10 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.6588 on 48 degrees of freedom
Multiple R-squared: 0.5748, Adjusted R-squared: 0.566
F-statistic: 64.89 on 1 and 48 DF, p-value: 1.814e-10
y_x3 <- lm(y~x3)
summary(y_x3)
Call:
lm(formula = y ~ x3)
Residuals:
Min 1Q Median 3Q Max
-1.5212 -0.5209 -0.1517 0.5198 2.4679
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.116e-17 1.284e-01 0.000 1.00000
x3 4.387e-01 1.297e-01 3.383 0.00144 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.9079 on 48 degrees of freedom
Multiple R-squared: 0.1925, Adjusted R-squared: 0.1757
F-statistic: 11.44 on 1 and 48 DF, p-value: 0.001437
変数選択を行った場合、まずyと一番相関が強いx1が取り込まれる。
x1のyに対する回帰を行い、残差を求める。①
resi.y_x1 <- y_x1$residuals
ようするに、x1がyを説明した部分を取り除いたことになる。
x1のx2、x3に対する回帰を行う、残差を求める。
x1のx2に対する回帰と残差。②
x2_x1 <- lm(x2~x1)
resi.x2_x1 <- x2_x1$residuals
x1のx3に対する回帰と残差。③
x3_x1 <- lm(x3~x1)
resi.x3_x1 <- x3_x1$residuals
ようするに、x1がx2とx3を説明していた部分を取り除いたということだ。
②の①に対する回帰と③の①に対する回帰を行う。
resi.y_x1_resi.x2_x1 <- lm(resi.y_x1~resi.x2_x1)
summary(resi.y_x1_resi.x2_x1)
Call:
lm(formula = resi.y_x1 ~ resi.x2_x1)
Residuals:
Min 1Q Median 3Q Max
-1.33597 -0.31891 0.03114 0.34418 1.02701
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.679e-17 6.957e-02 0.000 1.000
resi.x2_x1 4.949e-02 1.346e-01 0.368 0.715
Residual standard error: 0.4919 on 48 degrees of freedom
Multiple R-squared: 0.002808, Adjusted R-squared: -0.01797
F-statistic: 0.1351 on 1 and 48 DF, p-value: 0.7148
resi.y_x1_resi.x3_x1 <- lm(resi.y_x1~resi.x3_x1)
summary(resi.y_x1_resi.x3_x1)
Call:
lm(formula = resi.y_x1 ~ resi.x3_x1)
Residuals:
Min 1Q Median 3Q Max
-1.22041 -0.23250 0.02731 0.30354 1.03938
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.758e-17 6.521e-02 0.000 1.0000
resi.x3_x1 1.808e-01 6.944e-02 2.604 0.0122 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.4611 on 48 degrees of freedom
Multiple R-squared: 0.1238, Adjusted R-squared: 0.1055
F-statistic: 6.78 on 1 and 48 DF, p-value: 0.01223
③の①に対する回帰の方が係数が大きくなる。変数選択によりモデルにはx3が取り込まれることになる。
ちなみに、上で推定した係数は
lm(y~x1+x2)
Call: lm(formula = y ~ x1 + x2) Coefficients: (Intercept) x1 x2 -5.836e-17 8.309e-01 4.949e-02
と
lm(y~x1+x3)
Call: lm(formula = y ~ x1 + x3) Coefficients: (Intercept) x1 x3 -5.509e-17 8.159e-01 1.808e-01
に等しくなる。
今、モデルにはx1とx3が取り込まれている。
y_x1x3 <- lm(y~x1+x3)
summary(y_x1x3)
Call:
lm(formula = y ~ x1 + x3)
Residuals:
Min 1Q Median 3Q Max
-1.22041 -0.23250 0.02731 0.30354 1.03938
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -5.509e-17 6.590e-02 0.000 1.0000
x1 8.159e-01 7.017e-02 11.628 1.98e-15 ***
x3 1.808e-01 7.017e-02 2.577 0.0132 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.466 on 47 degrees of freedom
Multiple R-squared: 0.7917, Adjusted R-squared: 0.7828
F-statistic: 89.32 on 2 and 47 DF, p-value: < 2.2e-16
x1とx3でyを回帰したときの残差を求める。④
resi.y_x1x3 <- y_x1x3$residuals
x1とx3でx2を回帰し、残差を求める。⑤
x2_x1x3 <- lm(x2~x1+x3)
resi.x2_x1x3 <- x2_x1x3$residuals
上2つの残差で回帰を行う。（目的変数：④、説明変数：⑤）
resi.y_x1x3_resi.x2_x1x3 <- lm(resi.y_x1x3~resi.x2_x1x3)
summary(resi.y_x1x3_resi.x2_x1x3)
Call:
lm(formula = resi.y_x1x3 ~ resi.x2_x1x3)
Residuals:
Min 1Q Median 3Q Max
-1.25633 -0.24685 0.01037 0.28646 1.03247
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -4.543e-18 6.510e-02 0.000 1.000
resi.x2_x1x3 5.107e-02 1.260e-01 0.405 0.687
Residual standard error: 0.4603 on 48 degrees of freedom
Multiple R-squared: 0.003412, Adjusted R-squared: -0.01735
F-statistic: 0.1643 on 1 and 48 DF, p-value: 0.687
この係数は、x1、x2、x3をモデルに取り込んだときの、x2の係数に等しい。
lm(y~x1+x2+x3)
Call: lm(formula = y ~ x1 + x2 + x3) Coefficients: (Intercept) x1 x2 x3 -5.456e-17 7.723e-01 5.107e-02 1.809e-01
ある応答に影響のある要因をざっと見つけたいといった場合があったとする。
重回帰モデルを検定による変数選択で作成し、取り込まれた変数が応答に影響がありますと言ってしまうのは少々問題がある。取り込まれた変数はそれらが同時に（重回帰モデルによって）応答をよく説明しているのであって、それら各変数が応答に強い影響（もしくは相関）を持つということではないからである。
では、多数の要因から応答へ影響のあるものを拾い上げたい、しかも相関（すなわち直線関係）といった制約なしで、といった場合、「ランダムフォレスト」はどうだろうと考えている。
使ったことはあるが、そこまで分かっていない、今後の課題としてランダムフォレストを勉強する。
