読者です 読者をやめる 読者になる 読者になる

東京に棲む日々

データ分析、統計、ITを勉強中。未だ世に出ず。

重回帰と変数選択

重回帰とその変数選択に関するメモ。

RにはF値による検定でのステップワイズ変数選択法はないのでしょうか?
AICによる変数選択法である、step() 関数などは見当たるが検定による変数選択は見当たらない。Web検索でもヒットしない。

AICによる変数選択と検定による変数選択、どちらが良いとかそういうのは置いておいて、検定で行った時の変数の取り込みに関する基本的性質に関してメモしておきたい。

AICは非常に興味深い指標なので、そのうち詳しく勉強できたらいいと思う。

赤池情報量規準AIC―モデリング・予測・知識発見

赤池情報量規準AIC―モデリング・予測・知識発見


統計の専門家でなければ、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  

  

ある応答に影響のある要因をざっと見つけたいといった場合があったとする。

重回帰モデルを検定による変数選択で作成し、取り込まれた変数が応答に影響がありますと言ってしまうのは少々問題がある。取り込まれた変数はそれらが同時に(重回帰モデルによって)応答をよく説明しているのであって、それら各変数が応答に強い影響(もしくは相関)を持つということではないからである。

では、多数の要因から応答へ影響のあるものを拾い上げたい、しかも相関(すなわち直線関係)といった制約なしで、といった場合、「ランダムフォレスト」はどうだろうと考えている。
使ったことはあるが、そこまで分かっていない、今後の課題としてランダムフォレストを勉強する。