重回帰分析における各説明変数の説明力をAICで把握する

今回は重回帰分析における各説明変数の説明力をAICを使って把握する方法をメモしておきます。

方法としては、尤度比検定の要領でフルモデルから1つずつ説明変数を外したモデル（片方がもう片方をネストしている）を作成し、各説明変数について、尤度比検定でなくAICの差を計算するというやり方です。
AICや尤度比検定についてはこちらを参考に。

なぜこんなことをするかというと、サンプルサイズが小さい場合はP値でもわかりやすいのですが、例えば1000以上になるような膨大なデータ数になってくると、大して影響がない説明変数まで有意差が出やすくなり、どの説明変数が影響が大きいのかわかりにくくなるからです。

また、決定係数や疑似決定係数を使うやり方もありますが、これらにしてもP値にしても手持ちのデータに対しての当てはまりぐあいを見る考え方であり、予測という意味ではAICの方があっているのではないかと考えています。

今回はボストンの住宅価格のデータを使い、このやり方を試してみます。

ボストンの住宅価格のデータ

install.packages("mlbench")
library(mlbench)
data("BostonHousing")
BH<-BostonHousing
head(BH)

> head(BH)
     crim zn indus chas   nox    rm  age    dis rad tax ptratio      b lstat medv
1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98 24.0
2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14 21.6
3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03 34.7
4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94 33.4
5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33 36.2
6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21 28.7

ボストンの住宅のデータが14項目にわたって収められており、データ数は506個です。

一番右のmedvが住宅価格なのでこれを応答変数にし、残りを説明変数にして重回帰分析してみます。

重回帰分析を実行します。

重回帰分析

fit=glm(medv~crim+zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+b+lstat,data=BH,family=Gamma(link = "log"))

さすがに多重共線性がありそうなので、（というか説明変数多くてちょっと減らしたいので）vifをチェックします。

vifのチェック

library(car)
vif(fit)

vif

    crim       zn    indus     chas      nox       rm      age      dis      rad 
1.792192 2.298758 3.991596 1.073995 4.393720 1.933744 3.100826 3.955945 7.484496 
     tax  ptratio        b    lstat 
9.008554 1.799084 1.348521 2.941491

今回は厳しめに4以上のnoxとradとtaxを外します。
本当はvifが一番大きい説明変数から順番に除いていきますがここでは省略します。

再び重回帰分析

fit=glm(medv~crim+zn+indus+chas+rm+age+dis+ptratio+b+lstat,data=BH,family=Gamma(link = "log"))
summary(fit)

重回帰分析の結果

> summary(fit)

Call:
glm(formula = medv ~ crim + zn + indus + chas + rm + age + dis + 
    ptratio + b + lstat, family = Gamma(link = "log"), data = BH)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.4190463  0.1893489  18.057  < 2e-16 ***
crim        -0.0081191  0.0013606  -5.967 4.60e-09 ***
zn           0.0010587  0.0006081   1.741  0.08232 .  
indus       -0.0052163  0.0023077  -2.260  0.02423 *  
chas1        0.1288725  0.0391119   3.295  0.00105 ** 
rm           0.1001938  0.0187985   5.330 1.50e-07 ***
age         -0.0005766  0.0005821  -0.990  0.32243    
dis         -0.0435801  0.0087120  -5.002 7.88e-07 ***
ptratio     -0.0269744  0.0053215  -5.069 5.66e-07 ***
b            0.0004086  0.0001203   3.396  0.00074 ***
lstat       -0.0292869  0.0023107 -12.675  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for Gamma family taken to be 0.04716567)

    Null deviance: 81.425  on 505  degrees of freedom
Residual deviance: 20.391  on 495  degrees of freedom
AIC: 2909.3

Number of Fisher Scoring iterations: 7

次にMuMInパッケージを使用してすべての候補モデルをAICでランク付けします。
ベストモデルを選ぶ際に一般的なやり方です。
説明変数が多いので一覧は出しません。

AICでモデルのランク付け

library(MuMIn)
options(na.action = "na.fail")
dredged_models<-dredge(fit,rank="AIC")

さあ、ここからです。
説明変数が1つだけないモデルとフルモデルだけ抜き出します。
フルモデルとのAICの差を計算します。
AICの差の大きい順に並べ替えます。

AICは小さいほどよいわけですから、説明変数を外すことで、
AICが増加＝説明力がある（増加分だけ説明力が大きい）
AICが低下＝説明力がない
となります。

フルモデルから説明変数を外した際にAICが大きくなる順に並べ替え

library(tidyverse)
#データフレーム化
dredged_models<-data.frame(dredged_models)
#列数をカウント
ncol<-ncol(dredged_models)
#説明変数が1つだけ入ってないモデルとフルモデルを抜き出す
models_NA01<-dredged_models[rowSums(is.na(dredged_models[1:ncol]))<2,]
AIC_full<-AIC(fit)
#フルモデルとのAICの差を計算
models_NA01<-mutate(models_NA01,diff_AIC=models_NA01$AIC-AIC_full)
#差の大きい順に並び替え
models_NA01<-models_NA01[order(models_NA01$diff_AIC, decreasing = T),]
#余分な項目を除く
models_NA01<-models_NA01[, !(colnames(models_NA01) %in% c("df", "AIC", "logLik","delta","weight"))]
#AICを小数点以下2桁まで表示
models_NA01$diff_AIC<-round(models_NA01$diff_AIC,2)
models_NA01

並べかえの結果

> models_NA01
     X.Intercept.           age            b chas         crim         dis         indus       lstat     ptratio
960      2.480154 -3.158457e-03 0.0006734353    + -0.011044472 -0.04855710 -0.0086802438          NA -0.02866928
768      4.246231 -4.995367e-05 0.0003369839    + -0.007925229 -0.04900019 -0.0072269465 -0.03591912 -0.03162354
1016     3.429319 -6.038610e-04 0.0005956544    +           NA -0.03598893 -0.0057830802 -0.03116513 -0.03179909
1008     3.171648  5.895371e-04 0.0003955893    + -0.006965726          NA -0.0007288015 -0.02953843 -0.03186247
896      2.878505 -5.311774e-04 0.0004155527    + -0.009312604 -0.05111813 -0.0073732075 -0.02943768          NA
1020     3.433619 -4.449071e-04 0.0004326220 <NA> -0.008280932 -0.04464148 -0.0046791019 -0.02979423 -0.02859586
1022     3.625153 -5.413808e-04           NA    + -0.009306229 -0.04297365 -0.0060754592 -0.03039514 -0.02714480
992      3.343556 -7.486436e-04 0.0004422476    + -0.008270991 -0.03606366            NA -0.02981313 -0.02933390
512      3.439869 -7.020805e-04 0.0004024952    + -0.007734034 -0.03721005 -0.0052037628 -0.02919479 -0.02968623
1024     3.419046 -5.765953e-04 0.0004086231    + -0.008119108 -0.04358005 -0.0052162769 -0.02928690 -0.02697440
1023     3.396949            NA 0.0004060369    + -0.008130869 -0.04012474 -0.0055311257 -0.03001847 -0.02693127
             rm            zn diff_AIC
960  0.22072645  0.0009588535   184.86
768          NA  0.0014930530    33.99
1016 0.09904187  0.0005396440    33.26
1008 0.10996703 -0.0002361194    28.13
896  0.11387931  0.0019381090    27.35
1020 0.10214121  0.0010597923    11.12
1022 0.09487816  0.0010039414    10.54
992  0.10658589  0.0010529749     4.13
512  0.10422294            NA     1.69
1024 0.10019380  0.0010586664     0.00
1023 0.09729282  0.0011327790    -0.77

lstat、つまり「低所得者人口の割合」が圧倒的に説明力が大きいことがわかりました。
まあ、安いところに低所得者が集まるわけですから当然といえば当然という気もします。

それ以下は大きく差があいてrm（1戸当たりの「平均部屋数」）、crim（町別の「犯罪率」）、dis（「主要施設への距離」）、ptratio（町別の「生徒と先生の比率」）、chas1（「川の隣か」）、b（「町ごとの黒人の割合」）の順に大きくなっています。

また、サンプル数も500以上とそこそこ大きいことから、差が4.13のindus以下はほとんど説明力がないと判断してよいといえます。