こんにちは、たくろー(@takuro_data )です！
線形回帰分析について、理論をまとめ、Rのサンプルデータを使って単回帰分析、重回帰分析のケーススタディを行います。
要点
- 回帰分析は、変数x（原因）が変数y（結果）に与える影響を知るための方法
- 回帰線の精度は決定係数で評価
- 回帰係数は確率変数だから検定を行う
- 残差分析で回帰分析の適切さを判断
線形回帰分析の理論
線形回帰分析のRコマンド
lm(formula, data, weights, subset, na.action)| formula | 回帰式に用いる目的変数と説明変数、定数項などモデルの形式を指定。ex) y~x |
| data | 回帰分析に用いるデータセット名。 |
| weights | 説明変数に重みを指定できる。 |
| subset | データセットの一部を用いる指定をできる。 |
| na.action | 欠損値扱いを指定する引数。 |
回帰分析に多く用いられる関数は以下の通り（関数lmの結果をcars.lmに記録したとする）。
| 関数 | 内容 | 使用例 |
|---|---|---|
| coef | 回帰係数 | coef(cars.lm), cars.lm$coef |
| fitted | 用いたデータの予測値 | fitted(cars.lm), cars.lm$fitted |
| deviance | 残差の平方和 | deviance(cars.lm) |
| anova | 回帰係数の分散分析 | anova(cars.lm) |
| predict | 新たなデータに対する予測値 | predict(cars.lm) |
| 要約より簡単な結果 | print(cars.lm) | |
| summary | 回帰分析結果の要約 | summary(cars.lm) |
| plot | 回帰診断プロット | plot(cars.lm) |
線形単回帰分析ケーススタディ
Rのサンプルデータcarsを用いて分析してみます。carsは車の速度とブレーキをかけた後、車が止まるまでの距離に関して50回の実験結果を記録したデータです。
## データ準備
#サンプルデータ読み込み
data('cars')
head(cars)
#散布図
plot(cars)
#相関係数
cor(cars$speed, cars$dist)
実行結果
> head(cars)
speed dist
1 4 2
2 4 10
3 7 4
4 7 22
5 8 16
6 9 10
> plot(cars);cor(cars$speed, cars$dist)
[1] 0.8068949
散布図からspeedとdistに相関関係があるようにみえ、相関係数は0.80です。
## 単回帰分析
#モデル当てはめ
cars.lm <- lm(dist~speed, data=cars)
summary(cars.lm)
実行結果
> summary(cars.lm)
Call:
lm(formula = dist ~ speed, data = cars)
Residuals:
Min 1Q Median 3Q Max
-29.069 -9.525 -2.272 9.215 43.201
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.5791 6.7584 -2.601 0.0123 *
speed 3.9324 0.4155 9.464 1.49e-12 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared: 0.6511, Adjusted R-squared: 0.6438
F-statistic: 89.57 on 1 and 48 DF, p-value: 1.49e-12
Residuals:が残差の四分位数、Coefficients:が回帰係数の検定結果です。回帰係数の検定結果から、切片、speedの回帰係数ともに、p値が小さく有意な回帰係数です。回帰式は以下の通りで、speedが１増えると、止まるまでの距離が3.93増えることがわかります。
dist = -17.58 + 3.93 × speed
この回帰式の説明力は、0.65（決定係数）です。元データのグラフに得られた回帰直線を描いてみます。
#グラフによる分析
plot(cars)
abline(cars.lm, lwd=2)
最後に得られた回帰式が適切か残差分析をします。
#残差分析
par(mfrow=c(2,2))
plot(cars.lm)
| 左上：残差とフィット値の散布図 | 残差の全体像を表している。無相関で問題なさそうである。 |
| 右上：残差の正規QQプロット | データが正規分布に従うと、点が直線上に並ぶ。残差は概ね直線上に並んでおり、標準正規分布に従うという仮定を満たしていそうである。 |
| 左下：残差の平方根プロット | 残差の変動状況を考察する。大きく外れている値はないようにみえる。 |
| 右下：残差と影響力プロット（梃子値とクックの距離） | 個別のデータが外れ値になり、回帰モデルに過剰な影響を与えているかどうかを判断する。49番の残差は若干大きいが、クックの距離は0.5は超えていないので問題なさそうである。 |
残差分析の結果、今回得られた回帰モデルは大きな問題はないことを確認できました。
梃子値とクックの距離
梃子値は、観測値のxからデータセット内のすべての観測値のxの平均までの距離を測定します。他のデータと比べて異常な予測値を持つ観測値を識別するために使用します。クックの距離は、i番目の観測値を使用して計算された係数と観測値を使用しないで計算された係数との間の距離に対する測度です。他のデータと比べて異常な予測値を持つ観測値や、モデルがあまり適合できない観測値を識別するために使用します。0.5以上であれば影響が大きく、1を超えると特異に大きいと言われています。
線形重回帰分析ケーススタディ
Rのサンプルデータairqualityで分析してみます。airqualityは、1973年5月から9月までのニューヨークの大気状態を、6つの変数で観測・記録した154の観測データです。
#データ確認
str(airquality)
#対散布図
pairs(airquality[,1:4])
#相関行列
mcor <- cor(airquality[,1:4],use = "complete")
mcor
library(corrplot)
corrplot(mcor)
実行結果
> str(airquality)
'data.frame': 153 obs. of 6 variables:
$ Ozone : int 41 36 12 18 NA 28 23 19 8 NA ...
$ Solar.R: int 190 118 149 313 NA NA 299 99 19 194 ...
$ Wind : num 7.4 8 12.6 11.5 14.3 14.9 8.6 13.8 20.1 8.6 ...
$ Temp : int 67 72 74 62 56 66 65 59 61 69 ...
$ Month : int 5 5 5 5 5 5 5 5 5 5 ...
$ Day : int 1 2 3 4 5 6 7 8 9 10 ...
> mcor
Ozone Solar.R Wind Temp
Ozone 1.0000000 0.3483417 -0.6124966 0.6985414
Solar.R 0.3483417 1.0000000 -0.1271835 0.2940876
Wind -0.6124966 -0.1271835 1.0000000 -0.4971897
Temp 0.6985414 0.2940876 -0.4971897 1.0000000
対散布図をみると、オゾン量は、風力(相関係数-0.61)および温度（0.70）とやや強い相関がありそうです。
#モデルあてはめ
air.lm <- lm(Ozone~Solar.R+Wind+Temp, data=airquality)
summary(air.lm) 実行結果
> 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切片、偏回帰係数とも有意水準5%で有意である。回帰式は以下の通り。
Ozone = -64.34 + 0.06 × Solar.R – 3.33 × Wind + 1.65 × Temp
この回帰式の決定係数は0.60、自由度調整済み決定係数は0.59です。
#回帰診断
plot(air.lm)
library(car)
vif(air.lm)
実行結果
> vif(air.lm)
Solar.R Wind Temp
1.095253 1.329070 1.431367117番の残差が大きくなっているものの、クックの距離0.5を超えていないので、大きな問題はなさそうです。VIF値も小さく、多重共線性の問題もないです。
交互作用モデル
上記の重回帰モデルでは、目的変数と説明変数間の相関関係のみを用いているが、説明変数間の関連性がある場合、交互作用を考慮したモデルを作成することで当てはまりを改善できることがあります。
#交互作用モデルあてはめ
air.lm2 <- lm(Ozone~(Solar.R+Wind+Temp)^2, data=airquality)
summary(air.lm2) 実行結果
> 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.69、自由度調整済み係数が0.67であり、交互作用を考慮していないモデル（決定係数0.60、自由度調整済み決定係数0.59）より当てはまりがよくなりました。ただ、偏回帰係数のp値が大きいものがあるので、よりよいモデルを構築するために、変数選択を行います。
#変数選択
air.lm3 <- step(air.lm2) 実行結果
> 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
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
38635 661.61
- Solar.R:Temp 1 2141.1 40776 665.60
- Wind:Temp 1 4339.8 42975 671.43 増減法（step関数の引数directionデフォルトは増減法）で変数選択を行いました。すべての説明変数の交互作用を考慮したモデルのAICが662.37からスタートし、Solar.R:Windを除くと、AICが661.61に減り、その他の変数を除いたモデルのAICはこれ以上小さくならないので計算が終了しています。
AIC：赤池情報量基準
統計モデルの良さを評価するための指標で、次のように計算されます。AIC=-2×(モデルの最大対数尤度)+2×(モデルのパラメータ数)
AICを用いてモデルを選択する場合は、AICの値が小さいモデルをよいモデルと評価します。
summary(air.lm3)
plot(air.lm3)
vif(air.lm3) 実行結果
> summary(air.lm3)
Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp + Solar.R:Temp + Wind:Temp,
data = airquality)
Residuals:
Min 1Q Median 3Q Max
-38.398 -10.889 -2.445 7.132 93.485
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -1.368e+02 6.414e+01 -2.133 0.035252 *
Solar.R -3.531e-01 1.750e-01 -2.018 0.046184 *
Wind 1.115e+01 4.259e+00 2.617 0.010182 *
Temp 2.451e+00 8.250e-01 2.971 0.003678 **
Solar.R:Temp 5.717e-03 2.370e-03 2.412 0.017589 *
Wind:Temp -1.863e-01 5.425e-02 -3.434 0.000852 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 19.18 on 105 degrees of freedom
(42 observations deleted due to missingness)
Multiple R-squared: 0.6828, Adjusted R-squared: 0.6677
F-statistic: 45.21 on 5 and 105 DF, p-value: 2.2e-16
> vif(air.lm3)
Solar.R Wind Temp Solar.R:Temp Wind:Temp
76.06339 68.63490 18.48100 92.31022 53.36437
VIF値が10を超えており、多重共線性が疑われるので、標準化して、もう一度モデルの当てはめを行います。
標準偏回帰係数
全ての変数を標準化して重回帰分析を行ったときの回帰係数を、標準偏回帰係数を呼びます。単位が揃うので、説明変数間で回帰係数の大きさを比較することができます。#標準化
airquality$sOzone <- scale(airquality$Ozone)
airquality$sSolar.R <- scale(airquality$Solar.R)
airquality$sWind <- scale(airquality$Wind)
airquality$sTemp <- scale(airquality$Temp)
#モデル当てはめ
air.lm4 <- lm(sOzone~(sSolar.R+sWind+sTemp)^2, data=airquality)
air.lm5 <- step(air.lm4)
summary(air.lm5)
#回帰診断
vif(air.lm5)
plot(air.lm5) 実行結果
> summary(air.lm5)
Call:
lm(formula = sOzone ~ sSolar.R + sWind + sTemp + sSolar.R:sTemp +
sWind:sTemp, data = airquality)
Residuals:
Min 1Q Median 3Q Max
-1.16400 -0.33009 -0.07413 0.21621 2.83392
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.13332 0.06169 -2.161 0.032961 *
sSolar.R 0.25174 0.06466 3.893 0.000174 ***
sWind -0.35948 0.06330 -5.679 1.22e-07 ***
sTemp 0.47603 0.06591 7.223 8.42e-11 ***
sSolar.R:sTemp 0.14774 0.06125 2.412 0.017589 *
sWind:sTemp -0.18835 0.05484 -3.434 0.000852 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.5815 on 105 degrees of freedom
(42 observations deleted due to missingness)
Multiple R-squared: 0.6828, Adjusted R-squared: 0.6677
F-statistic: 45.21 on 5 and 105 DF, p-value: 2.2e-16
> vif(air.lm5)
sSolar.R sWind sTemp sSolar.R:sTemp sWind:sTemp
1.393439 1.329312 1.432529 1.442913 1.118007 
残差の正規性が若干外れているものの、VIF値も小さく、大きな問題はないモデルを構築できました。
回帰式は以下のとおりで、決定係数は0.68、自由度調整済み決定係数は0.67です。それぞれの説明変数は正規化した値です。
Ozone* = -0.13 + 0.25 × Solar.R* – 0.35 × Wind* + 0.48 × Temp* + 0.15 × Solar.R*: Temp* – 0.18 × Wind*:Temp*
標準偏回帰係数の値の大きさを比べると、Ozoneへの正の影響が大きいのがTempで、負の影響が大きいのがWindであることがわかります。
さいごに
線形回帰分析について理論をまとめ、理論を使いながら、単回帰分析、重回帰分析をRで行いました。重回帰分析について、交互作用を考慮し、多重共線性を起こさないために変数を標準化してモデルを構築しました。当てはまりのよさはあがったものの、交互作用を考慮しないほうが解釈可能性が高いこともあります。
ビジネス目的によって、当てはまりのよさか、解釈可能性のバランスを判断して、モデル構築をしていきましょう。
以上、お読みいただきありがとうございました。
