統計ER

統計ソフトRの使い方を中心に、統計解析方法の解説をするブログ。ありそうでなかなか見つからないサンプルサイズ計算などニッチな方法について紹介しています。

統計ソフトRのISLRパッケージAutoデータを使った重回帰モデル

ブログランキングに参加しています。
まずはぽちぽちっとお願いします。
↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓↓
にほんブログ村 科学ブログ 数学へ

燃費を推測する重回帰モデルを作成した。

交互作用項を検討したり、VIFを考慮したり。

統計ソフトRのISLRパッケージを使った課題解決例。

 

 

データの準備

最初の一回だけ、ISLRパッケージをインストール。

install.packages("ISLR")

 

ISLRパッケージを呼び出して、解析開始。

library(ISLR)

 

ISLRパッケージのAutoデータセットを用いて解析。

9変数、392例のデータ。

> str(Auto)
'data.frame':   392 obs. of  9 variables:
 $ mpg         : num  18 15 18 16 17 15 14 14 14 15 ...
 $ cylinders   : num  8 8 8 8 8 8 8 8 8 8 ...
 $ displacement: num  307 350 318 304 302 429 454 440 455 390 ...
 $ horsepower  : num  130 165 150 150 140 198 220 215 225 190 ...
 $ weight      : num  3504 3693 3436 3433 3449 ...
 $ acceleration: num  12 11.5 11 12 10.5 10 9 8.5 10 8.5 ...
 $ year        : num  70 70 70 70 70 70 70 70 70 70 ...
 $ origin      : num  1 1 1 1 1 1 1 1 1 1 ...
 $ name        : Factor w/ 304 levels "amc ambassador brougham",..: 49 36 231 14 161 141 54 223 241 2 ...
> 

 

name以外の変数を使う。

Auto1 <- Auto[1:8]

 

Miles per gallonを予測する因子の相関係数

origin以外の7つの変数で相関係数行列を作成する。

> round(cor(Auto1[1:7]),3)
                mpg cylinders displacement horsepower weight acceleration   year
mpg           1.000    -0.778       -0.805     -0.778 -0.832        0.423  0.581
cylinders    -0.778     1.000        0.951      0.843  0.898       -0.505 -0.346
displacement -0.805     0.951        1.000      0.897  0.933       -0.544 -0.370
horsepower   -0.778     0.843        0.897      1.000  0.865       -0.689 -0.416
weight       -0.832     0.898        0.933      0.865  1.000       -0.417 -0.309
acceleration  0.423    -0.505       -0.544     -0.689 -0.417        1.000  0.290
year          0.581    -0.346       -0.370     -0.416 -0.309        0.290  1.000

 

Miles per gallonを国別に比較

origin別のmpg(miles per gallon=燃費)をANOVAで分析する。

originの内容は以下の通り。

Origin of car (1. American, 2. European, 3. Japanese)

res.origin <- lm(mpg ~ factor(origin), data=Auto1) 
summary(res.origin)

 

アメ車に比べヨーロッパ車は燃費がいいし、日本車はさらにいい。

Estimateがmiles per gallon推定値で、

アメ車に+7.6がヨーロッパ車、+10.4が日本車。

> summary(res.origin)

Call:
lm(formula = mpg ~ factor(origin), data = Auto1)

Residuals:
    Min      1Q  Median      3Q     Max 
-12.451  -5.034  -1.034   3.649  18.966 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      20.0335     0.4086  49.025   <2e-16 ***
factor(origin)2   7.5695     0.8767   8.634   <2e-16 ***
factor(origin)3  10.4172     0.8276  12.588   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.396 on 389 degrees of freedom
Multiple R-squared:  0.3318,    Adjusted R-squared:  0.3284 
F-statistic:  96.6 on 2 and 389 DF,  p-value: < 2.2e-16

 

Miles per gallonを予測する重回帰モデルの作成

全部の変数でmpgを推測するモデルを作成してみる。

res1 <- lm(mpg ~ cylinders + displacement + horsepower + weight + acceleration + year + factor(origin), data=Auto1)

summary(res1)

 

Adjusted R-squaredが0.8205と高く、この重回帰モデルの説明変数は目的変数をよく説明している。

displacement (排気量), weight, year, originは統計学的有意に関連あり。

例えば、yearが1増加する(1年新しくなる)と、mpgが0.777増加(燃費向上)する。

> summary(res1)

Call:
lm(formula = mpg ~ cylinders + displacement + horsepower + weight + 
    acceleration + year + factor(origin), data = Auto1)

Residuals:
    Min      1Q  Median      3Q     Max 
-9.0095 -2.0785 -0.0982  1.9856 13.3608 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -1.795e+01  4.677e+00  -3.839 0.000145 ***
cylinders       -4.897e-01  3.212e-01  -1.524 0.128215    
displacement     2.398e-02  7.653e-03   3.133 0.001863 ** 
horsepower      -1.818e-02  1.371e-02  -1.326 0.185488    
weight          -6.710e-03  6.551e-04 -10.243  < 2e-16 ***
acceleration     7.910e-02  9.822e-02   0.805 0.421101    
year             7.770e-01  5.178e-02  15.005  < 2e-16 ***
factor(origin)2  2.630e+00  5.664e-01   4.643 4.72e-06 ***
factor(origin)3  2.853e+00  5.527e-01   5.162 3.93e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.307 on 383 degrees of freedom
Multiple R-squared:  0.8242,    Adjusted R-squared:  0.8205 
F-statistic: 224.5 on 8 and 383 DF,  p-value: < 2.2e-16

 

残差プロット。予測式の性能チェック。

Residuals vs. FittedとScale-Locationは水平な赤い線が理想。

Normal Q-Qはすべての点が斜めの線上であるのが理想。

Residuals vs. Leverageは赤い点線より外側に点がないのが理想。

Residual vs. Fittedは、厳密にいうと、heteroscedasticity(分散不均一性)に見える。

Normal Q-Qは、Quartileが大きいところでずれがある。

Scale-LocationとResiduals vs. Leverageは取り立てて問題点はない。

layout(t(matrix(c(1:4),nr=2)))
plot(res1)

f:id:toukeier:20190908213122p:plain

交互作用項を入れてみた!

例として、yearとdisplacementの交互作用項を追加して再度解析。

res2 <- lm(mpg ~ cylinders + displacement + horsepower + weight + acceleration + year + year:displacement + factor(origin), data=Auto1)

summary(res2)

 

displacementとyearの交互作用項を足してみたところ、

displacement:year: p=3.51e-12 ***

という結果を得た。この項は統計学的有意であった。

> summary(res2)

Call:
lm(formula = mpg ~ cylinders + displacement + horsepower + weight + 
    acceleration + year + year:displacement + factor(origin), 
    data = Auto1)

Residuals:
    Min      1Q  Median      3Q     Max 
-7.5860 -2.0066  0.0133  1.6529 12.1025 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -6.451e+01  7.829e+00  -8.241 2.77e-15 ***
cylinders          2.606e-02  3.103e-01   0.084   0.9331    
displacement       2.806e-01  3.642e-02   7.704 1.15e-13 ***
horsepower        -3.370e-02  1.306e-02  -2.580   0.0103 *  
weight            -6.116e-03  6.212e-04  -9.844  < 2e-16 ***
acceleration       1.048e-01  9.237e-02   1.134   0.2574    
year               1.375e+00  9.645e-02  14.261  < 2e-16 ***
factor(origin)2    2.588e+00  5.323e-01   4.861 1.71e-06 ***
factor(origin)3    2.453e+00  5.224e-01   4.695 3.73e-06 ***
displacement:year -3.560e-03  4.953e-04  -7.187 3.51e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.107 on 382 degrees of freedom
Multiple R-squared:  0.8451,    Adjusted R-squared:  0.8415 
F-statistic: 231.6 on 9 and 382 DF,  p-value: < 2.2e-16

 

多重共線性を考慮してみた!

多重共線性のチェックのためにVIF(Variance Inflation Factor)を計算。

cylindersとdisplacementに多重共線性が疑われた(VIF>10)。

> round(1/(1-cor(Auto1)^2),1)
             mpg cylinders displacement horsepower weight acceleration year origin
mpg          Inf       2.5          2.8        2.5    3.3          1.2  1.5    1.5
cylinders    2.5       Inf         10.4        3.5    5.1          1.3  1.1    1.5
displacement 2.8      10.4          Inf        5.1    7.7          1.4  1.2    1.6
horsepower   2.5       3.5          5.1        Inf    4.0          1.9  1.2    1.3
weight       3.3       5.1          7.7        4.0    Inf          1.2  1.1    1.5
acceleration 1.2       1.3          1.4        1.9    1.2          Inf  1.1    1.0
year         1.5       1.1          1.2        1.2    1.1          1.1  Inf    1.0
origin       1.5       1.5          1.6        1.3    1.5          1.0  1.0    Inf

 

mpgと相関の強いdisplacementだけ残して、再解析。

res11 <- lm(mpg ~ displacement + horsepower + weight + acceleration + year + factor(origin), data=Auto1)

summary(res11)

 

> summary(res11)

Call:
lm(formula = mpg ~ displacement + horsepower + weight + acceleration + 
    year + factor(origin), data = Auto1)

Residuals:
   Min     1Q Median     3Q    Max 
-8.729 -2.115 -0.021  1.923 13.292 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -19.247521   4.607317  -4.178 3.65e-05 ***
displacement      0.016436   0.005849   2.810  0.00521 ** 
horsepower       -0.015629   0.013629  -1.147  0.25220    
weight           -0.006837   0.000651 -10.502  < 2e-16 ***
acceleration      0.086660   0.098262   0.882  0.37837    
year              0.778280   0.051867  15.005  < 2e-16 ***
factor(origin)2   2.591975   0.566840   4.573 6.50e-06 ***
factor(origin)3   2.770110   0.550988   5.028 7.64e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.312 on 384 degrees of freedom
Multiple R-squared:  0.8231,    Adjusted R-squared:  0.8199 
F-statistic: 255.3 on 7 and 384 DF,  p-value: < 2.2e-16

 

モデル診断は情報量規準で!

AIC赤池情報量規準)、BICベイズ情報量規準)の二つでモデル診断。

当てはまりの良さをチェック。

AIC(res1, res11, res2)
BIC(res1, res11, res2)

 

値が小さいほうが当てはまりがよい。

AICBICも交互作用項を入れたモデルが最も当てはまりがよい。

> AIC(res1, res11, res2)
      df      AIC
res1  10 2060.928
res11  9 2061.299
res2  11 2013.216
> BIC(res1, res11, res2)
      df      BIC
res1  10 2100.640
res11  9 2097.041
res2  11 2056.900

 

まとめ

統計ソフトRのISLRパッケージにあるAutoデータで、重回帰分析をやってみた。

順番としては、先行研究や理論をもとにモデルを作成し、残差プロットや情報量規準で当てはまりを診断する。