[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
解析その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にあがりました！
解析その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をもとにモデル•変数選択ができ、便利だなと思いました。
