線形回帰分析まとめ:理論とRを使ってケーススタディ

こんにちは、タクロウ(@takuro_109)です!

線形回帰分析について、理論をまとめ、Rのサンプルデータを使って単回帰分析、重回帰分析のケーススタディを行います。

要点
  1. 回帰分析は、変数x(原因)が変数y(結果)に与える影響を知るための方法
  2. 回帰線の精度は決定係数で評価
  3. 回帰係数は確率変数だから検定を行う
  4. 残差分析で回帰分析の適切さを判断

線形回帰分析の理論

線形回帰分析の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要約より簡単な結果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が1増えると、止まるまでの距離が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.431367

117番の残差が大きくなっているものの、クックの距離0.5を超えていないので、大きな問題はなさそうです。VIF値も10を超えるものがないので、多重共線性の問題もないです。

交互作用モデル

上記の重回帰モデルでは、目的変数と説明変数間の相関関係のみを用いているが、説明変数間の関連性がある場合、交互作用を考慮したモデルを作成することで当てはまりを改善できることがあります。

コード
#交互作用モデルあてはめ
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値も10以内に収まり大きな問題はないモデルを構築できました。

回帰式は以下のとおりで、決定係数は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で行いました。重回帰分析について、交互作用を考慮し、多重共線性を起こさないために変数を標準化してモデルを構築しました。当てはまりのよさはあがったものの、交互作用を考慮しないほうが解釈可能性が高いこともあります。

ビジネス目的によって、当てはまりのよさか、解釈可能性のバランスを判断して、モデル構築をしていきましょう。

以上、お読みいただきありがとうございました。