10. 回帰分析

10.1. 単回帰分析

単回帰分析は、ある1つの変数(独立変数)からもう1つの変数(従属変数)を説明する式を求めます。

図10-1 回帰分析のイメージ

具体的には $Y=β_0+β_1X_1+ε$ という回帰式を立てて $β_1$(傾き)と $β_0$(切片)にあたる数値を求めます。$ε$ は誤差と呼ばれており、従属変数の観測値と式によって求められる従属変数の予測値との差のことです。さらにその数式が独立変数と従属変数との関係をきちんと説明できるかを検討します。

扱うデータ

小売業者の取り組みに関する架空のデータ advert.csv を利用して、「広告費から売上げを予測する式」を推定します。

図10-2 広告費と売上のデータ
In [1]:
#作業ディレクトリの処理は既に行っていれば不要(復習)
getwd() #作業ディレクトリの確認
setwd("R:/WinPython/WinPython-64bit-3.6.1.0Qt5/notebooks/data") #作業ディレクトリの変更
'R:/WinPython/WinPython-64bit-3.6.1.0Qt5/notebooks'
In [2]:
data = read.csv("10-1.csv", header = TRUE) #データの挿入
attach(data)

要約統計量

In [3]:
summary(data) #データビュー
plot(data) #散布図の作成
     advert           sale      
 Min.   :14.93   Min.   : 80.5  
 1st Qu.:27.02   1st Qu.:106.1  
 Median :47.99   Median :133.0  
 Mean   :44.64   Mean   :126.8  
 3rd Qu.:58.91   3rd Qu.:145.3  
 Max.   :76.93   Max.   :152.9  

単回帰分析の実行

In [4]:
result = lm(sale ~ advert) #lm関数の実行
summary(result) #lm関数の出力確認
Call:
lm(formula = sale ~ advert)

Residuals:
     Min       1Q   Median       3Q      Max 
-19.2993  -4.3818  -0.1658   5.3570  15.7968 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  79.0058     4.8201   16.39 8.14e-14 ***
advert        1.0713     0.1001   10.70 3.46e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8.864 on 22 degrees of freedom
Multiple R-squared:  0.8389,	Adjusted R-squared:  0.8316 
F-statistic: 114.6 on 1 and 22 DF,  p-value: 3.455e-10

まず、Multiple R-squaredは $R^2$乗(決定係数)を表しています。これは、回帰式によって説明することができる観測値のばらつきの割合です。$R^2$乗の値を参照すると0.8389となっており、従属変数の分散のうち83.9%を、回帰式で説明することができることを意味します。

次にCoefficientsの表を見ます。定数=79.0058,広告費=1.0713という回帰係数が得られました。そこから、 $売上げ = 79.0058 + 1.0713 × 広告費 + ε$ という単回帰式が得られます。この式から広告費が1単位増えると、売り上げが 1.0713 増えるという解釈をすることができます。

有意確率 $Pr(>|t|)$ の項目では、「回帰係数の値は0ではない」という仮説を誤って捨てて「回帰係数の値は0である」という仮説を取る確率がどのくらいあるかを確認します。有意水準を 0.05 とすると広告費の有意確率 3.46e-10(3.46×10の-10乗)はそれより小さいため、回帰係数は統計的に有意であるという結論が導かれます。

回帰直線の作図

In [5]:
plot(data)
abline(result, col=2)

10.2 重回帰分析

単回帰分析は1つの従属変数に対して、1つの独立変数を対応させましたが、2つ以上の独立変数を対応させるには重回帰分析を行います。

具体的には $Y = β_0 + β_1 X_1 + β_2 X_2 +⋯+ β_n X_n+ε$ の式を求めることになります。

扱うデータ

ここではRのデータセットairqualityを利用します。データセットairqualityはニューヨークの大気状態を6つの変数で観測した、153件の記録です。ここで、太陽の放射量(Solar.R)、風力(Wind)、温度(Temp)からオゾンの量(Ozone)を予測する式を推定してみましょう。

In [6]:
(data = airquality) #データの代入
summary(data) #要約量軽量の算出
attach(data)
result = lm(Ozone~Solar.R+Wind+Temp) #lm関数による重回帰分析の実行
summary(result) #結果の出力
OzoneSolar.RWindTempMonthDay
41 190 7.467 5 1
36 118 8.072 5 2
12 149 12.674 5 3
18 313 11.562 5 4
NA NA 14.356 5 5
28 NA 14.966 5 6
23 299 8.665 5 7
19 99 13.859 5 8
8 19 20.161 5 9
NA 194 8.669 5 10
7 NA 6.974 5 11
16 256 9.769 5 12
11 290 9.266 5 13
14 274 10.968 5 14
18 65 13.258 5 15
14 334 11.564 5 16
34 307 12.066 5 17
6 78 18.457 5 18
30 322 11.568 5 19
11 44 9.762 5 20
1 8 9.759 5 21
11 320 16.673 5 22
4 25 9.761 5 23
32 92 12.061 5 24
NA 66 16.657 5 25
NA 266 14.958 5 26
NA NA 8.057 5 27
23 13 12.067 5 28
45 252 14.981 5 29
115 223 5.779 5 30
..................
96 167 6.991 9 1
78 197 5.192 9 2
73 183 2.893 9 3
91 189 4.693 9 4
47 95 7.487 9 5
32 92 15.584 9 6
20 252 10.980 9 7
23 220 10.378 9 8
21 230 10.975 9 9
24 259 9.773 9 10
44 236 14.981 9 11
21 259 15.576 9 12
28 238 6.377 9 13
9 24 10.971 9 14
13 112 11.571 9 15
46 237 6.978 9 16
18 224 13.867 9 17
13 27 10.376 9 18
24 238 10.368 9 19
16 201 8.082 9 20
13 238 12.664 9 21
23 14 9.271 9 22
36 139 10.381 9 23
7 49 10.369 9 24
14 20 16.663 9 25
30 193 6.970 9 26
NA 145 13.277 9 27
14 191 14.375 9 28
18 131 8.076 9 29
20 223 11.568 9 30
     Ozone           Solar.R           Wind             Temp      
 Min.   :  1.00   Min.   :  7.0   Min.   : 1.700   Min.   :56.00  
 1st Qu.: 18.00   1st Qu.:115.8   1st Qu.: 7.400   1st Qu.:72.00  
 Median : 31.50   Median :205.0   Median : 9.700   Median :79.00  
 Mean   : 42.13   Mean   :185.9   Mean   : 9.958   Mean   :77.88  
 3rd Qu.: 63.25   3rd Qu.:258.8   3rd Qu.:11.500   3rd Qu.:85.00  
 Max.   :168.00   Max.   :334.0   Max.   :20.700   Max.   :97.00  
 NA's   :37       NA's   :7                                       
     Month            Day      
 Min.   :5.000   Min.   : 1.0  
 1st Qu.:6.000   1st Qu.: 8.0  
 Median :7.000   Median :16.0  
 Mean   :6.993   Mean   :15.8  
 3rd Qu.:8.000   3rd Qu.:23.0  
 Max.   :9.000   Max.   :31.0  
                               
Call:
lm(formula = Ozone ~ Solar.R + Wind + Temp)

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

まず単回帰分析の場合とは違い「調整済み$R^2$乗(Adjusted R-squared)」の数値を確認します。調整済みR2乗は単回帰分析で確認した$R^2$乗の数値と同じ意味を持っています。今回の例では、従属変数の分散のうち59.48%を、回帰式で説明することができることを意味します。

次にCoefficientsの表から、定数 = -64.34208,太陽の放射量(Solar.R)= 0.05982、風力(Wind)= -3.33359、温度(Temp)= 1.65209という係数が得られました。これらの係数から、 $オゾンの量(Ozone)=0.05982×太陽の放射量(Solar.R)-3.33359×風力(Wind)+ 1.65209×温度(Temp)$ という重回帰式が得られます。

有意確率は単回帰分析の時と同じように解釈します。今回の例では、有意水準を 0.05 とすると太陽の放射量(Solar.R)、風力(Wind)、温度(Temp)の有意確率は 0.05 より小さいため、太陽の放射量(Solar.R)、風力(Wind)、温度(Temp)の回帰係数はそれぞれ統計的に有意であるという結論が導かれます。

多重共線性

最後に多重共線性を確認します。多重共線性とは、独立変数の間に強い相関関係がある状態です。独立変数の間に相関関係があると、回帰係数の有意性が下がり、正確な分析を行うことができなくなります。

多重共線性を確認するには、VIF(分散拡大係数)の数値を確認します。一般的に、VIFの数値が10以上の独立変数は多重共線性によって対応する回帰係数の値が信頼できないと判断します。

RでVIFを求めるには、carパッケージのvif関数を使用します。

In [ ]:
!install.packages("car", repos = c(CRAN ='https://cran.ism.ac.jp/'))
In [ ]:
library(car) #carパッケージを使用
vif(result) # VIF>10だと多重共線性に問題あり。

今回の例では、太陽の放射量(Solar.R)、風力(Wind)、温度(Temp)のVIFは10以下なので、多重共線性の問題はないと考えられます。

10.3. 関数オプション

関数 目的
lm(y ~ x)  単回帰分析
plot(data) 散布図の作図
lm(y~x1+x2+...+xn) 重回帰分析