時系列データの重回帰分析で見るべき3つの指標(R^2, VIF, DW)

f:id:ekanoh:20180811021451p:plain

時系列の重回帰分析で見ておきたい3つの指標

Rで時系列データに対する重回帰分析を行ったら、まずは以下の3つの指標を見ておきたい。

  1. 自由度調整済み決定係数
  2. VIF
  3. DW検定

具体的に、時系列解析・多変量解析の教科書などでもお馴染みのairqualityデータ(に私が手を加えたもの)に対して重回帰分析を行ったことにして結果を見ていこう。

ちなみに:airqualityデータとは

1973年5月から9月までのニューヨークの大気状態データ。N=約150。

f:id:ekanoh:20180811020344p:plain

課題設定

タスクとしては、日照強度・風速・気温などのデータを使って「オゾン層の平均濃度」を推定したい。

つまり、「日照強度200で、風速が10で、気温が華氏70度、・・・」と伝えると、「それならオゾン層の濃度はYくらいだね」と見込みを教えてくれる関数を作ることになる。

従って、以下のように変数を選択した。

  • 目的変数
  • 説明変数
    • Solar.R: 日照強度
    • Wind: 風速 [miles/h]
    • Temp: 気温 [華氏°F]
    • Rainfall: 降雨量
  • 使用しない変数
    • Month: 月
    • Day: 日

※Rainfallを見て何この変数?と思った方は正解です。後ほど説明します。

分析結果

さて、実際にRのlm()関数を利用して分析結果を出力すると以下のようになる。

> res <- lm(Ozone ~ Solar.R + Wind + Temp + Rainfall, data = df)
> summary(res)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-0.94400 -0.46757 -0.03372  0.45310  1.34176 

Coefficients:
              Estimate Std. Error  t value Pr(>|t|)    
(Intercept) -0.8538419  0.6301284   -1.355    0.178    
Solar.R      0.0009634  0.0006306    1.528    0.130    
Wind        -0.0171082  0.0192479   -0.889    0.376    
Temp        -1.9870952  0.0114558 -173.457   <2e-16 ***
Rainfall     0.1987903  0.0005077  391.518   <2e-16 ***
---
Signif. codes:  0***0.001**0.01*0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5594 on 106 degrees of freedom
  (42 observations deleted due to missingness)
Multiple R-squared:  0.9997, Adjusted R-squared:  0.9997 
F-statistic: 9.728e+04 on 4 and 106 DF,  p-value: < 2.2e-16

冒頭でも紹介した以下の3つの指標を見て、モデルを評価していこう。

  1. 自由度調整済み決定係数
  2. VIF
  3. DW検定

R2(決定係数)

f:id:ekanoh:20180809192852p:plain

自由度調整済み決定係数を説明する前に、決定係数を説明しておこう。

決定係数は、回帰式の当てはまりのよさ(goodness of fit)を示した指標で、上記の出力結果でMultiple R-squaredと書かれているのがそれだ。目的変数をどれだけ説明変数で説明できているか、の指標である。

従って、高ければ高いほど回帰式の当てはまりがよいということになる。

Adjusted R2(自由度調整済み決定係数)

f:id:ekanoh:20180809192907p:plain

しかしながら、一般に決定係数R2は、説明変数が増えれば増えるほど大きくなるという性質を持っている。変数が増えたとき、適合度が大きくなることはあっても、少なくなることは直観的にも考えにくいだろう。

とはいえ、もし決定係数がさほど変わらないなら、本当は少ない変数で目的変数を説明できている方がいいはずである。

なので、変数の数(≒自由度)を使って決定係数の値を微調整したものが自由度調整済み決定係数で、これを使えば決定係数よりも正当な評価ができる。

VIF (Variance Inflation Factor)

f:id:ekanoh:20180811002319p:plain

VIF (Variance Inflation Factor)とは、重回帰分析を行なった際に各変数に多重共線性(Multicolinearity)が無いかを調べるための指標である。

多重共線性とは、説明変数の間に相関がある状態を表し、これがあると分析が上手くいかないことが多い。

まずは、以下のOzoneとTempの散布図を見てほしい。

f:id:ekanoh:20180811015550p:plain

上記の通り、与えられたデータにおいてはオゾン層の濃度と気温の間に正の相関があるように見込まれる。

しかしながら、カンのいい方はすでに気付いたかもしれないが、上記のRの出力結果ではTempの係数が負になっているのである!どういうことだろうか?このような場合は、まずは多重共線性を疑おう。

RでVIFを出力してみると、以下のようになる。

> library(car)
> vif(res)
 Solar.R     Wind     Temp Rainfall 
1.161255 1.648315 4.189520 5.364369 

TempとRainfallのVIFが他と比べて高いので、少なくともこの2つの変数は多重共線性を持ってしまっている可能性が高い。(なお、一般的には10以上だと要チェックだと言われている。)

実はそれもそのはずで、Rainfallは5×Ozone+10×Temp+乱数という式で表され、OzoneやTempと相関を持つように人工的に生成した数なのである。

多重共線性が疑われる場合は、例えばどちらかの指標を説明変数から外してしまおう。Rainfallを説明変数から除外してみると、以下のようにTempの係数が正になっているのが確認できた。

> res <- lm(Ozone ~ Solar.R + Wind + Temp, data = df)
> summary(res)

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

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

多重共線性を数式的に理解してみる

たとえば、 {y = a_1x_1 + a_2x_2+b}というようなシンプルな重回帰式を考えてみる。ここで、 {x_1} {x_2}に相関がある場合を考えてみる。たとえば、 {x_2 = 2x}という関係式が成り立っているとしよう。

すると、 {y=(a_1+2a_2)x_1 + b}である。もし、本当の関係式が
 {y=3x_1}だったとすると、 {a_1=a_2=1}でも解になるのだが、 {a_1=5, a_2=-1}も解になりうる(!)

DW (Durbin-Watson Test)

f:id:ekanoh:20180811002149p:plain

DW (Durbin-Watson Test)とは、誤差項(実測値と理論値の差)間に自己相関があるかないかを判別するための指標。

回帰分析では、異なる誤差項間には相関がないことを仮定しているが、時系列データにおいては、1時点前と現時点での値が(自己)相関してしまうことが多い。

この時、回帰分析の結果がふつうよりもかなり当てはまりがよくなってしまう傾向にある。これを、見せかけの回帰などと呼ぶ。

一般に、ダービン=ワトソン比の値が2よりかなり小さいときは正の相関が、2よりかなり大きいときは負の相関が、2前後のときは相関なしと判断する。

例のモデルのDW値をrで出力してみると自己相関は低い値、ダービン=ワトソン比も2に近く、p値も高いため自己相関は無いものとして考えてよい。

> durbinWatsonTest(res)
 lag Autocorrelation D-W Statistic p-value
   1      0.03150895      1.935476   0.656
 Alternative hypothesis: rho != 0

まとめ

  1. 自由度調整済み決定係数
    • 回帰式の当てはまりのチェック
      • つまり目的変数をどれだけ説明変数で説明できているか
    • 自由度(変数の数)で調整をしている
    • 1に近いほどよい
  2. VIF
    • 多重共線性のチェック
    • 一般に10を超えると×
  3. DW検定
    • 残差の自己相関のチェック
    • 有意だと×
      • 決定係数の値がよくても見せかけの回帰の可能性などがある

おまけ:まとめて出力してしまおう

上記の3つの指標は時系列の重回帰分析においては比較的重要なため、たとえば以下のような自作関数を作ってまとめて出力してしまうのがおすすめだ。

lms <- function(formula, data){
     require(car)
     # Linear Regression
     res <- lm(formula, data=data)
     
     # Output Stats
     print(summary(res))
     print("VIF:")
     print(vif(res))
     print("Durbin-Watson Test")
     print(durbinWatsonTest(res))
     
     return(res)
}

上記の自作関数を使って重回帰分析を実行すると、自由度調整済み決定係数、VIF、DWの3つの指標がすべて出力できる。

> res <- lms(Ozone ~ Solar.R + Wind + Temp, data = df)

Call:
lm(formula = formula, data = data)

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

[1] "VIF:"
 Solar.R     Wind     Temp 
1.095253 1.329070 1.431367 
[1] "Durbin-Watson Test"
 lag Autocorrelation D-W Statistic p-value
   1      0.03150895      1.935476   0.606
 Alternative hypothesis: rho != 0

Rによるデータサイエンス データ解析の基礎から最新手法まで

Rによるデータサイエンス データ解析の基礎から最新手法まで