[R]回帰分析 - 線形重回帰分析

今回は、線形重回帰分析。内容は、Rによるデータサイエンス - データ解析の基礎から最新手法までの7.3節に沿ってやります。

重回帰分析とは

説明変数が複数ある回帰分析を、重回帰分析と呼びます。

説明変数のデータをX、目的変数をY、誤差をE、係数をAと表すと、回帰モデルは、

となり、係数は、最小２乗法などで求めることが出来ます。

用いるデータ

Rに用意されている、airqualityというデータセットを使います。1973年5月から9月までのニューヨークの大気状態を、6つの変数で観測し、記録した154の観測値です。

pairs(airquality[,1:4])

解析その1

とりあえず、Ozoneを目的変数として、Solar.R、Wind、Tempを説明変数とした重回帰分析を行ってみます。

air.lm <- lm(Ozone~Solar.R+Wind+Temp, data=airquality)
summary(air.lm)

要約

Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp, data = airquality)

Residuals:
    Min      1Q  Median      3Q     Max 
-40.485 -14.219  -3.551  10.097  95.619 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -64.34208   23.05472  -2.791  0.00623 ** 
Solar.R       0.05982    0.02319   2.580  0.01124 *  
Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
Temp          1.65209    0.25353   6.516 2.42e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 21.18 on 107 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-squared: 0.6059,	Adjusted R-squared: 0.5948 
F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16 

AIC

モデルの評価として、AICというものがあり、これは小さいほど良いとされています。
Rには、AICを求める関数としてextractAIC関数が用意されており、モデルのパラメータの数とAIC値を返します。

extractAIC(air.lm)

[1]   4.0000 681.7127

パラメータの数が4つで、AICが681.7127と出ました。

解析その2

散布図を見ると、説明変数間に、この図で言えば、WindとTempに相関がある事が分かります。
このような場合、説明変数一個に対する影響が、目的変数だけでなく、もう一方の説明変数、Windの変化の影響がTempにも出てしまい、適切なモデルとは言えません。
そこで、この説明変数間の相関関係（相互作用と言います）を考慮したモデルを構築してみます。
Rでは、「^2」をつけるだけで、これが出来ます。

air.lm2 <- lm(Ozone~(Solar.R+Wind+Temp)^2, data=airquality)
summary(air.lm2)

要約

Call:
lm(formula = Ozone ~ (Solar.R + Wind + Temp)^2, data = airquality)

Residuals:
    Min      1Q  Median      3Q     Max 
-38.685 -11.727  -2.169   7.360  91.244 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -1.408e+02  6.419e+01  -2.193  0.03056 * 
Solar.R      -2.260e-01  2.107e-01  -1.073  0.28591   
Wind          1.055e+01  4.290e+00   2.460  0.01555 * 
Temp          2.322e+00  8.330e-01   2.788  0.00631 **
Solar.R:Wind -7.231e-03  6.688e-03  -1.081  0.28212   
Solar.R:Temp  5.061e-03  2.445e-03   2.070  0.04089 * 
Wind:Temp    -1.613e-01  5.896e-02  -2.735  0.00733 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Residual standard error: 19.17 on 104 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-squared: 0.6863,	Adjusted R-squared: 0.6682 
F-statistic: 37.93 on 6 and 104 DF,  p-value: < 2.2e-16 

調整済み決定係数が、0.5948から0.6682にあがりました！

AIC

AICも確認してみまると、こちらも良くなったことが分かります。

extractAIC(air.lm2)

[1]   7.0000 662.3732

解析その3

解析その2の結果を見ると、モデルに用いた説明変数のp値が、やや大きいものがあります。
より良いモデルを作成するためには、変数を何らかの基準で選択して用いる必要があります。
Rでは、AICを用いてモデル•å¤‰æ•°ã‚’選択する関数stepが用意されています。

air.lm3 <- step(air.lm2)

Start:  AIC=662.37
Ozone ~ (Solar.R + Wind + Temp)^2

               Df Sum of Sq   RSS    AIC
- Solar.R:Wind  1    429.42 38635 661.61
<none>                      38205 662.37
- Solar.R:Temp  1   1574.75 39780 664.86
- Wind:Temp     1   2748.20 40954 668.08

Step:  AIC=661.61
Ozone ~ Solar.R + Wind + Temp + Solar.R:Temp + Wind:Temp

               Df Sum of Sq   RSS    AIC
<none>                      38635 661.61
- Solar.R:Temp  1    2141.1 40776 665.60
- Wind:Temp     1    4339.8 42975 671.43

結果、AICは661.61になりました。

また、モデルは、Solar.R:Windの相互作用を削除した次のモデルになります。

• 係数を出力

round(coefficients(air.lm3),2)

 (Intercept)      Solar.R         Wind         Temp Solar.R:Temp    Wind:Temp 
     -136.81        -0.35        11.15         2.45         0.01        -0.19 

• 回帰モデル

Ozone = -136.81 -0.35*Solar.R + 11.15*Wind + 2.45*Temp + 0.01*Solar.R*Temp -0.19*Wind*Temp

回帰診断図

最後に回帰診断図。

正規Q-Qプロットを見ると、両端が直線から外れており、残差が正規分布に従うという仮定が成り立ってるとは言い難いことが分かります。
さらに当てはまりを良くするためには、残差の分布を正規分布以外で扱える他のモデルが必要になってきます。ここら辺については、また別エントリーで書きます。

まとめ

• Rで線形重回帰分析を行いました。
• 説明変数間の相関関係を相互作用といい、lm関数のformulaで「^2」を付ける事でこれを考慮したモデルを作れます。
• また、step関数で、AICをもとにモデル•å¤‰æ•°é¸æŠžãŒã§ãã€ä¾¿åˆ©ã ãªã¨æ€ã„ました。