Rで重回帰分析をしてみた

重回帰分析とは、複数の説明変数をもちいて目的変数を説明する分析手法のことです。

今回はRのサンプルデータを使って重回帰分析を実行してみたいと思います。

まずは、今回必要なパッケージをインストールしておきましょう。

> devtools::install_github('jaredlander/coefplot')
> install.packages("GGally")
> install.packages("dplyr")
> install.packages("corrplot")
> install.packages("ggplot2")
> install.packages("tidyr")
> install.packages("caret")

重回帰分析の実行

今回は、MASSパッケージのBostonデータセットを使って重回帰分析を実行してみます。難しそうに思うかもしれないですが、lm関数を使うだけで重回帰分析を実行することができますので安心してください。

ちなみに、Bostonのデータセットは米国のボストン市の住宅データのことです。特徴量が14、サンプル数が506となっています。ここでは、medv(所有者の持ち家の価格:1000ドル単位)を目的変数、その他の13項目を説明変数として重回帰分析を実行してみます。

> library(MASS)
# ボストンのデータをロード
> data(Boston)
# 重回帰分析の実行
> lm.Boston <- lm(medv ~ ., data = Boston)
# 分析結果の要約
> summary(lm.Boston)

Call:
lm(formula = medv ~ ., data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.595  -2.730  -0.518   1.777  26.199 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
zn           4.642e-02  1.373e-02   3.382 0.000778 ***
indus        2.056e-02  6.150e-02   0.334 0.738288    
chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
age          6.922e-04  1.321e-02   0.052 0.958229    
dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
black        9.312e-03  2.686e-03   3.467 0.000573 ***
lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.745 on 492 degrees of freedom
Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338 
F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

「Coefficients」には、各回帰係数の推定値(Estimate)、標準偏差(Std. Error)、t-値(t value)、P-値(Pr(>|t|))が表示されています。 回帰係数の右側にはP-値の大きさに応じて(アスタリスク)が付与されています。***はP-値が0~0.001、**は0.001~0.01、は0.01から0.05、.は0.05から0.1、印がない場合は0.1から1の値の範囲にあることを示しています。

「Multiple R-squared」、「Adjusted R-squared」には、決定係数及び調整済みの決定係数が表示されています。この場合は、それぞれ0.7406、0.7338となっています。

回帰係数のプロット

回帰係数とは、前項で求めた重回帰分析の各説明変数の係数及び切片を指します。重回帰分析でも回帰分析の結果から、推定された回帰分析の情報を抽出できます。

ここでは、回帰結果のlm.Bostonオブジェクトをsummary関数で要約した後に、coef関数で回帰係数を抽出しています。coef関数は、statsパッケージで提供されていて、第1引数には回帰係数を抽出かいのうなオブジェクトを指定します。

# 回帰係数の表示
> coef(summary(lm.Boston))
                 Estimate  Std. Error      t value     Pr(>|t|)
(Intercept)  3.645949e+01 5.103458811   7.14407419 3.283438e-12
crim        -1.080114e-01 0.032864994  -3.28651687 1.086810e-03
zn           4.642046e-02 0.013727462   3.38157628 7.781097e-04
indus        2.055863e-02 0.061495689   0.33431004 7.382881e-01
chas         2.686734e+00 0.861579756   3.11838086 1.925030e-03
nox         -1.776661e+01 3.819743707  -4.65125741 4.245644e-06
rm           3.809865e+00 0.417925254   9.11614020 1.979441e-18
age          6.922246e-04 0.013209782   0.05240243 9.582293e-01
dis         -1.475567e+00 0.199454735  -7.39800360 6.013491e-13
rad          3.060495e-01 0.066346440   4.61289977 5.070529e-06
tax         -1.233459e-02 0.003760536  -3.28000914 1.111637e-03
ptratio     -9.527472e-01 0.130826756  -7.28251056 1.308835e-12
black        9.311683e-03 0.002685965   3.46679256 5.728592e-04
lstat       -5.247584e-01 0.050715278 -10.34714580 7.776912e-23

上記の結果は説明変数が多い時には見づらくなってしまうので、coefplotパッケージのcoefplot関数を使用することで結果を見やすくすることができます。

library(coefplot)
coefplot(lm.Boston)

f:id:yunosuke1107:20180422034923p:plain

上のグラフでは、縦軸に説明変数、横軸に回帰係数がプロットされています。●で表された点が推定値(Estimate)、線は推定値の区間を表しています。線の左端が推定値から標準偏差(Std. error)を引いた点、右端が推定値に標準偏差を加算した点です。

変数を選択

先ほどの例ではすべての変数を使用して重回帰分析を実行していましたが、step関数を用いると、目的変数に対する説明力が低い変数を削除して回帰を実行してくれます。

> lm.Boston.step <- step(lm.Boston)
Start:  AIC=1589.64
medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + 
    tax + ptratio + black + lstat

          Df Sum of Sq   RSS    AIC
- age      1      0.06 11079 1587.7
- indus    1      2.52 11081 1587.8
<none>                 11079 1589.6
- chas     1    218.97 11298 1597.5
- tax      1    242.26 11321 1598.6
- crim     1    243.22 11322 1598.6
- zn       1    257.49 11336 1599.3
- black    1    270.63 11349 1599.8
- rad      1    479.15 11558 1609.1
- nox      1    487.16 11566 1609.4
- ptratio  1   1194.23 12273 1639.4
- dis      1   1232.41 12311 1641.0
- rm       1   1871.32 12950 1666.6
- lstat    1   2410.84 13490 1687.3

Step:  AIC=1587.65
medv ~ crim + zn + indus + chas + nox + rm + dis + rad + tax + 
    ptratio + black + lstat

          Df Sum of Sq   RSS    AIC
- indus    1      2.52 11081 1585.8
<none>                 11079 1587.7
- chas     1    219.91 11299 1595.6
- tax      1    242.24 11321 1596.6
- crim     1    243.20 11322 1596.6
- zn       1    260.32 11339 1597.4
- black    1    272.26 11351 1597.9
- rad      1    481.09 11560 1607.2
- nox      1    520.87 11600 1608.9
- ptratio  1   1200.23 12279 1637.7
- dis      1   1352.26 12431 1643.9
- rm       1   1959.55 13038 1668.0
- lstat    1   2718.88 13798 1696.7

Step:  AIC=1585.76
medv ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + 
    black + lstat

          Df Sum of Sq   RSS    AIC
<none>                 11081 1585.8
- chas     1    227.21 11309 1594.0
- crim     1    245.37 11327 1594.8
- zn       1    257.82 11339 1595.4
- black    1    270.82 11352 1596.0
- tax      1    273.62 11355 1596.1
- rad      1    500.92 11582 1606.1
- nox      1    541.91 11623 1607.9
- ptratio  1   1206.45 12288 1636.0
- dis      1   1448.94 12530 1645.9
- rm       1   1963.66 13045 1666.3
- lstat    1   2723.48 13805 1695.0

lm.Bostonオブジェクトにstep関数を適用して変数選択を実行しています。最初に、

Start:  AIC=1589.64
medv ~ crim + zn + indus + chas + nox + rm + age + dis + rad + 
    tax + ptratio + black + lstat

と表示されています。これは、全13変数を用いて重回帰分析を実行して得られた回帰式はAICが1589.64であることを表しています。

変数選択をした後の回帰式をsummary関数で要約してみます。

> summary(lm.Boston.step)

Call:
lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad + 
    tax + ptratio + black + lstat, data = Boston)

Residuals:
     Min       1Q   Median       3Q      Max 
-15.5984  -2.7386  -0.5046   1.7273  26.2373 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
crim         -0.108413   0.032779  -3.307 0.001010 ** 
zn            0.045845   0.013523   3.390 0.000754 ***
chas          2.718716   0.854240   3.183 0.001551 ** 
nox         -17.376023   3.535243  -4.915 1.21e-06 ***
rm            3.801579   0.406316   9.356  < 2e-16 ***
dis          -1.492711   0.185731  -8.037 6.84e-15 ***
rad           0.299608   0.063402   4.726 3.00e-06 ***
tax          -0.011778   0.003372  -3.493 0.000521 ***
ptratio      -0.946525   0.129066  -7.334 9.24e-13 ***
black         0.009291   0.002674   3.475 0.000557 ***
lstat        -0.522553   0.047424 -11.019  < 2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 4.736 on 494 degrees of freedom
Multiple R-squared:  0.7406, Adjusted R-squared:  0.7348 
F-statistic: 128.2 on 11 and 494 DF,  p-value: < 2.2e-16

以上の結果を見ると、最終的に選択された変数はcrim、zn、chas、nox、rm、dis、rad、tax、ptratio、black、lstatの11個であることがわかります。 この結果はもともとの回帰式でこれら2つの説明変数のP-値が相対的に高く、回帰係数が0であるとする帰無仮説を棄却できなかったことからも納得のいく結果です。 変数選択を実行した後の回帰式の決定係数、調整済み決定係数は、それぞれ0.7406、0.7348となっています。これらの値は変数選択前に比べてそれほど変わっていません。。。 このことからも、切り捨てられた変数ageやindusのキヨがもともとそれほどなかったことがわかりますね。

各変数の相関関係

本来なら重回帰分析する前に各変数の相関関係などを把握しておくべきですが、、、。 一応今回使用したBostonデータセットの中身を確認してみましょう。 まずは相関係数を見てみます。

> library(dplyr)
> cor.Boston <- Boston %>% cor
> cor.Boston
               crim          zn       indus         chas         nox          rm         age
crim     1.00000000 -0.20046922  0.40658341 -0.055891582  0.42097171 -0.21924670  0.35273425
zn      -0.20046922  1.00000000 -0.53382819 -0.042696719 -0.51660371  0.31199059 -0.56953734
indus    0.40658341 -0.53382819  1.00000000  0.062938027  0.76365145 -0.39167585  0.64477851
chas    -0.05589158 -0.04269672  0.06293803  1.000000000  0.09120281  0.09125123  0.08651777
nox      0.42097171 -0.51660371  0.76365145  0.091202807  1.00000000 -0.30218819  0.73147010
rm      -0.21924670  0.31199059 -0.39167585  0.091251225 -0.30218819  1.00000000 -0.24026493
age      0.35273425 -0.56953734  0.64477851  0.086517774  0.73147010 -0.24026493  1.00000000
dis     -0.37967009  0.66440822 -0.70802699 -0.099175780 -0.76923011  0.20524621 -0.74788054
rad      0.62550515 -0.31194783  0.59512927 -0.007368241  0.61144056 -0.20984667  0.45602245
tax      0.58276431 -0.31456332  0.72076018 -0.035586518  0.66802320 -0.29204783  0.50645559
ptratio  0.28994558 -0.39167855  0.38324756 -0.121515174  0.18893268 -0.35550149  0.26151501
black   -0.38506394  0.17552032 -0.35697654  0.048788485 -0.38005064  0.12806864 -0.27353398
lstat    0.45562148 -0.41299457  0.60379972 -0.053929298  0.59087892 -0.61380827  0.60233853
medv    -0.38830461  0.36044534 -0.48372516  0.175260177 -0.42732077  0.69535995 -0.37695457
                dis          rad         tax    ptratio       black      lstat       medv
crim    -0.37967009  0.625505145  0.58276431  0.2899456 -0.38506394  0.4556215 -0.3883046
zn       0.66440822 -0.311947826 -0.31456332 -0.3916785  0.17552032 -0.4129946  0.3604453
indus   -0.70802699  0.595129275  0.72076018  0.3832476 -0.35697654  0.6037997 -0.4837252
chas    -0.09917578 -0.007368241 -0.03558652 -0.1215152  0.04878848 -0.0539293  0.1752602
nox     -0.76923011  0.611440563  0.66802320  0.1889327 -0.38005064  0.5908789 -0.4273208
rm       0.20524621 -0.209846668 -0.29204783 -0.3555015  0.12806864 -0.6138083  0.6953599
age     -0.74788054  0.456022452  0.50645559  0.2615150 -0.27353398  0.6023385 -0.3769546
dis      1.00000000 -0.494587930 -0.53443158 -0.2324705  0.29151167 -0.4969958  0.2499287
rad     -0.49458793  1.000000000  0.91022819  0.4647412 -0.44441282  0.4886763 -0.3816262
tax     -0.53443158  0.910228189  1.00000000  0.4608530 -0.44180801  0.5439934 -0.4685359
ptratio -0.23247054  0.464741179  0.46085304  1.0000000 -0.17738330  0.3740443 -0.5077867
black    0.29151167 -0.444412816 -0.44180801 -0.1773833  1.00000000 -0.3660869  0.3334608
lstat   -0.49699583  0.488676335  0.54399341  0.3740443 -0.36608690  1.0000000 -0.7376627
medv     0.24992873 -0.381626231 -0.46853593 -0.5077867  0.33346082 -0.7376627  1.0000000

算出した相関係数を可視化してみます。今回はcorrplotパッケージを使います。

> library(corrplot)
> corrplot(Boston %>% cor, addConf.col = TRUE)[f:id:yunosuke1107:20180422042437p:plain]

f:id:yunosuke1107:20180422042437p:plain

▼分散拡大要因のヒートマップ

library(tidyr)
library(ggplot2)
# 分散拡大要因の算出
vif <- 1/(1 - (Boston %>% dplyr:::select(-medv) %>% cor)^2)
# 縦長形式への変換
vif.l <- vif %>%
  as.data.frame %>%
  mutate(item1 = rownames(.)) %>%
  gather(item2, vif, -item1)
p <- ggplot(data = vif.l, aes(x = item1, y = item2, fill = vif)) +
  geom_tile() + geom_text(label = sprintf("%.2f", vif)) +
  scale_fill_gradient(low = "white", high = "red")
print(p)

f:id:yunosuke1107:20180422043855p:plain

広告を非表示にする