【重回帰モデル】回帰の診断（2）――変数はすべて適切か？ - Men talking over coffee with smoking Ark Royal.

※重回帰モデルについての、小暮厚之先生の講義ノート（2012年）を基にした覚書です。

この記事は、「変数をとりあえず全部入れて回帰分析をして、これとそれが有意だったので、意味のある影響を持つ変数はこれとそれで〜」のような分析を駆逐したいという強い思いから書きました。
赤池先生の情報量基準を使って変数選択をして、真のモデルにできるだけ近づけたモデルから議論をしましょう！

回帰の診断：回帰モデルを推定した後に・・・

異常な観測値はないか？（外れ値の検出）←前回の内容
変数はすべて適切か？（説明変数の選択[AIC, BIC]、多重共線性）←今回はここをやります。

使用したデータなどは、前回（【重回帰モデル】回帰の診断（1）――異常な観測値はないか？）を参照してください。

AICとその修正版――回帰分析における変数選択の基準
※R2乗、自由度調整済みR2乗は基準足りえる根拠がない。

AIC(BIC)が小さいモデルほどよいモデルとみなす。

AIC：「将来の予測」に適する

$AIC=n\log(\frac{RSS}{n})+2K$
RSS: 残差二乗和, K=説明変数の個数

BIC：「現在及び過去の説明」に適する

$BIC=n\log(\frac{RSS}{n})+\log(n)K$

■AICによる説明変数の選択

> step(Term1.kaiki3, criterion="AIC", direction="forward")
Start:  AIC=228
LNFACE ~ LNINCOME + EDUCATION + NUMHH + d197


Call:
lm(formula = LNFACE ~ LNINCOME + EDUCATION + NUMHH + d197, data = Term1)

Coefficients:
(Intercept)     LNINCOME    EDUCATION        NUMHH         d197  
     2.0494       0.5749       0.1843       0.2824       5.0623

AICは228で、先ほどのモデルと同じ変数が選択されました。

では、次に、TermLife.csvに入っていたすべての変数から変数選択を行います。

> names(Term2)
 [1] "GENDER"             "AGE"                "MARSTAT"           
 [4] "EDUCATION"          "ETHNICITY"          "SMARSTAT"          
 [7] "SGENDER"            "SAGE"               "SEDUCATION"        
[10] "NUMHH"              "INCOME"             "TOTINCOME"         
[13] "CHARITY"            "FACE"               "FACECVLIFEPOLICIES"
[16] "CASHCVLIFEPOLICIES" "BORROWCVLIFEPOL"    "NETVALUE"          
[19] "LNFACE"             "LNINCOME"          

> x = names(Term2[-c(11, 14, 19)])
> # xはTerm2からINCOME, FACE, LNFACEを除いたもの
> fmla = as.formula(paste("LNFACE~", paste(x, collapse="+")))
> fmla
LNFACE ~ GENDER + AGE + MARSTAT + EDUCATION + ETHNICITY + SMARSTAT + 
    SGENDER + SAGE + SEDUCATION + NUMHH + TOTINCOME + CHARITY + 
    FACECVLIFEPOLICIES + CASHCVLIFEPOLICIES + BORROWCVLIFEPOL + 
    NETVALUE + LNINCOME
> fmla #回帰分析のformulaを作成
LNFACE ~ GENDER + AGE + MARSTAT + EDUCATION + ETHNICITY + SMARSTAT + 
    SGENDER + SAGE + SEDUCATION + NUMHH + TOTINCOME + CHARITY + 
    FACECVLIFEPOLICIES + CASHCVLIFEPOLICIES + BORROWCVLIFEPOL + 
    NETVALUE + LNINCOME

> Term2.kaiki = lm(fmla, data=Term2)
> summary(Term2.kaiki)

Call:
lm(formula = fmla, data = Term2)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.4824 -0.8606  0.1178  0.9367  4.0784 

Coefficients:
                     Estimate Std. Error t value Pr(>|t|)    
(Intercept)         5.111e+00  1.206e+00   4.239 3.14e-05 ***
GENDER              7.466e-01  3.916e-01   1.906 0.057716 .  
AGE                -1.635e-02  1.193e-02  -1.371 0.171575    
MARSTAT            -6.695e-01  4.624e-01  -1.448 0.148894    
EDUCATION           1.709e-01  4.674e-02   3.657 0.000310 ***
ETHNICITY          -3.945e-02  6.655e-02  -0.593 0.553799    
SMARSTAT            2.040e-01  1.998e-01   1.021 0.308321    
SGENDER            -3.914e-01  5.771e-01  -0.678 0.498272    
SAGE                1.103e-02  1.424e-02   0.775 0.439325    
SEDUCATION          4.239e-02  5.153e-02   0.823 0.411484    
NUMHH               2.448e-01  7.773e-02   3.149 0.001829 ** 
TOTINCOME           3.918e-09  2.844e-08   0.138 0.890536    
CHARITY             6.300e-06  2.452e-06   2.569 0.010769 *  
FACECVLIFEPOLICIES  4.933e-08  3.830e-08   1.288 0.198908    
CASHCVLIFEPOLICIES  1.175e-06  1.081e-06   1.087 0.278078    
BORROWCVLIFEPOL     5.244e-02  4.726e-02   1.110 0.268226    
NETVALUE            2.155e-02  2.412e-01   0.089 0.928890    
LNINCOME            3.342e-01  8.553e-02   3.907 0.000119 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.483 on 257 degrees of freedom
Multiple R-squared: 0.4104,	Adjusted R-squared: 0.3714 
F-statistic: 10.52 on 17 and 257 DF,  p-value: < 2.2e-16

全部入れモデルでは、EDUCATION, NUMHH, CHARITY, LNINCOMEが有意になっています。これらの変数は選択されるのでしょうか。また、これら以外の変数は選択されないのでしょうか。
ここではもうひとつの情報量基準、BICを用いて変数選択をしてみます。

> step(Term2.kaiki, criterion="BIC", direction="backward")
Start:  AIC=234.2
LNFACE ~ GENDER + AGE + MARSTAT + EDUCATION + ETHNICITY + SMARSTAT + 
    SGENDER + SAGE + SEDUCATION + NUMHH + TOTINCOME + CHARITY + 
    FACECVLIFEPOLICIES + CASHCVLIFEPOLICIES + BORROWCVLIFEPOL + 
    NETVALUE + LNINCOME

                     Df Sum of Sq    RSS    AIC
- NETVALUE            1     0.018 565.40 232.21
- TOTINCOME           1     0.042 565.42 232.22
- ETHNICITY           1     0.773 566.15 232.57
- SGENDER             1     1.012 566.39 232.69
- SAGE                1     1.320 566.70 232.84
- SEDUCATION          1     1.489 566.87 232.92
- SMARSTAT            1     2.292 567.67 233.31
- CASHCVLIFEPOLICIES  1     2.599 567.98 233.46
- BORROWCVLIFEPOL     1     2.708 568.09 233.51
- FACECVLIFEPOLICIES  1     3.649 569.03 233.97
<none>                            565.38 234.20
- AGE                 1     4.135 569.51 234.20
- MARSTAT             1     4.611 569.99 234.43
- GENDER              1     7.995 573.37 236.06
- CHARITY             1    14.517 579.90 239.17
- NUMHH               1    21.821 587.20 242.61
- EDUCATION           1    29.415 594.79 246.15
- LNINCOME            1    33.584 598.96 248.07

# …中略…

Step:  AIC=221.35
LNFACE ~ GENDER + EDUCATION + NUMHH + CHARITY + FACECVLIFEPOLICIES + 
    BORROWCVLIFEPOL + LNINCOME

                     Df Sum of Sq    RSS    AIC
<none>                            580.28 221.35
- BORROWCVLIFEPOL     1     4.303 584.58 221.38
- FACECVLIFEPOLICIES  1     6.770 587.05 222.54
- GENDER              1    15.116 595.39 226.42
- CHARITY             1    18.088 598.36 227.79
- NUMHH               1    32.243 612.52 234.22
- LNINCOME            1    37.571 617.85 236.60
- EDUCATION           1    65.242 645.52 248.65

Call:
lm(formula = LNFACE ~ GENDER + EDUCATION + NUMHH + CHARITY + 
    FACECVLIFEPOLICIES + BORROWCVLIFEPOL + LNINCOME, data = Term2)

Coefficients:
       (Intercept)              GENDER           EDUCATION  
         3.611e+00           7.399e-01           2.082e-01  
             NUMHH             CHARITY  FACECVLIFEPOLICIES  
         2.446e-01           6.663e-06           5.732e-08  
   BORROWCVLIFEPOL            LNINCOME  
         6.133e-02           3.427e-01

GENDER, EDUCATION, NUMHH, CHARITY, FACECVLIFEPOLICIES, BORROWCVLIFEPOL, LNINCOMEを用いるモデルが選択されました。最初の全部入れモデルよりもAICが234.2から221.35へと減少しました。

では、このモデルではどの変数が有意となるのでしょうか。

> Term2.BIC = lm(LNFACE ~ GENDER + EDUCATION + NUMHH + CHARITY + FACECVLIFEPOLICIES + BORROWCVLIFEPOL + LNINCOME, data=Term2)
> summary(Term2.BIC)

Call:
lm(formula = LNFACE ~ GENDER + EDUCATION + NUMHH + CHARITY + 
    FACECVLIFEPOLICIES + BORROWCVLIFEPOL + LNINCOME, data = Term2)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.5707 -0.9019  0.1245  0.9090  4.1471 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        3.611e+00  8.877e-01   4.069 6.23e-05 ***
GENDER             7.399e-01  2.806e-01   2.637 0.008848 ** 
EDUCATION          2.082e-01  3.799e-02   5.479 9.89e-08 ***
NUMHH              2.446e-01  6.350e-02   3.852 0.000147 ***
CHARITY            6.663e-06  2.310e-06   2.885 0.004234 ** 
FACECVLIFEPOLICIES 5.732e-08  3.247e-08   1.765 0.078710 .  
BORROWCVLIFEPOL    6.133e-02  4.359e-02   1.407 0.160548    
LNINCOME           3.427e-01  8.241e-02   4.158 4.33e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 1.474 on 267 degrees of freedom
Multiple R-squared: 0.3949,	Adjusted R-squared: 0.379 
F-statistic: 24.89 on 7 and 267 DF,  p-value: < 2.2e-16

全部入れモデルでは、EDUCATION, NUMHH, CHARITY, LNINCOMEが有意でしたが、変数選択後のモデルでは、GENDER, EDUCATION, NUMHH, CHARITY, LNINCOMEが有意で、新しくGENDERが有意となりました。
このGENDERは全部入れモデルでは他の変数の影響を受けて、ゴチャゴチャになり埋もれていたと考えられます。変数選択は絶対ではありませんが、変数選択後のモデルを用いたほうが「真のモデル」により近い、誠実な結果が得られるはずです。
これを機に全部入れモデルだけから、「これとそれが有意だった」と結論付けるやり方はぜひ辞めてください。お願いします。

※多重共線性については次回に見送ります。