「見せかけの回帰」の復習 - 渋谷駅前で働くデータサイエンティストのブログ

先日のことですが、Querie.meでこんな質疑がありました。

これは非常にご尤もなご意見であり、実際この問題提起に近いシチュエーションを見かけたことは五本の指では数え切れないくらいあります。ということで、今回の記事では元々の問題意識ともいえる「見せかけの回帰」について、久しぶりにちょっと復習を兼ねて書いてみようと思います。

そもそも「見せかけの回帰」とは何か
実際に見せかけの回帰において起きること
見せかけの回帰への対処法
Rコード

そもそも「見せかけの回帰」とは何か

このブログでは11年前に沖本本の輪読記事を書いた際に「見せかけの回帰」については一通り取り上げていますので、今回はその際の説明を引用するに留めます。

なお前提知識として先に書いておくと、以下に出てくる「単位根過程」というのは平たく言えば「前の時点の値にランダムなノイズを積み増して次の時点の値が得られる」タイプのトレンドを伴う時系列データです。完全なホワイトノイズを累積すると「ランダムウォーク」になりますが、これは単位根過程の典型例です。この定義に当てはまるデータは世の中に数多く、例えばファイナンスやマーケティング分野のデータはその好例とされます。

単位根過程同士を回帰すると、本来存在しないはずの回帰関係（見せかけの回帰）が勝手に生じてしまう

既に以前の記事でも紹介済みですが、再掲しておきましょう。要は、単位根過程の被説明変数に対して同じく単位根過程の説明変数によって単なるOLS（最小二乗法）による推定をかけると、偏回帰係数 $\beta_1, \beta_2, \cdots, \beta_n$ に関するt統計量が勝手に発散してしまい、その結果としてp値も爆発的に小さくなり「有意な回帰がある」という誤った結果になってしまう現象がある、ということです。

一応、沖本本pp.126-127の記述を引用しておきます。ブラウン運動の知識があると、前回の記事の時よりは理解しやすいと思います。

見せかけの回帰の現象を定義するために、2つの独立なランダムウォークを考えよう。具体的には、 $\epsilon_1t$ と $\epsilon_2t$ を独立なiid系列として、 $x_t$ と $y_t$ を
$\begin{array} x_t=x_{t-1}+\epsilon_1t, \epsilon_1t \sim iid(0,\sigma^2_1) \\ y_t=y_{t-1}+\epsilon_2t, \epsilon_2t \sim iid(0,\sigma^2_1) \end{array}$
で定義する。このとき、 $x_t$ と $y_t$ は独立なランダムウォークであるので、
$y_t=\alpha + \beta x_t + \epsilon_t$ (6.3)
となる回帰モデルを考えたとすると、真の値は $\beta = 0$ となる。したがって、(6.3)をOLSで推定し、 $\beta = 0$ という仮説検定を行うと、 $\beta = 0$ が採択される確率が高いはずである。しかしながら、変数の非定常性が驚くべき結果をもたらすことが知られている。具体的には、(6.3)をOLSで推定したときのOLS推定量 $\hat{\alpha}$ と $\hat{\beta}$ について、
$\left( \begin{array} ~ T^{-1/2} \hat{\alpha} \\ \hat{\beta} \end{array} \right) \stackrel{L}\longrightarrow \left( \begin{array} ~ \sigma_1 h_1 \\ (\sigma_1 / \sigma_2) h_2 \end{array} \right)$ (6.4)
が成立することが知られている*1。ここで $W_1(\cdot)$ と $W_2(\cdot)$ を独立な標準ブラウン運動として、
$\left( \begin{array} ~ h_1 \\ h_2 \end{array} \right)$
$= \left( \begin{array}{cc} 1 & \int W_2(r)dr \\ \int W_2(r)dr & \int W_2^2(r)dr \end{array} \right)^{-1} \left( \begin{array} \int W_1(r)dr \\ \int W_2(r)W_1(r)dr \end{array} \right)$
である。(6.4)の結果は、 $\hat{\alpha}$ が $\sqrt{T}$ の速度で発散し、 $\hat{\beta}$ がある確率変数に収束することを示している。これは、 $\sqrt{T}$ の速度で真の値に収束していく標準的な場合や、 $T$ の速度で真の値に収束していく単位根検定統計量の場合とは非常に異なるものであり、 $\hat{\alpha}$ も $\hat{\beta}$ も一致推定量ではないことがわかる。また、このため、 $\hat{\alpha}$ と $\hat{\beta}$ に関するt統計量は発散することになる。つまり、 $T$ が大きいとき、 $\hat{\alpha}$ と $\hat{\beta}$ を用いてt検定を行うと、ほぼ確実に $\hat{\alpha}=0$ と $\hat{\beta}=0$ という帰無仮説は棄却されるのである*2。また、もう1つの驚くべき事実として、回帰の決定係数 $R^2$ が漸近的に1に収束することも知られている。

この奇怪な現象は、最初にGrangerとNewboldが当時始まったばかりの計算機によるモンテカルロシミュレーションによって発見し(Granger & Newbold, 1974)、その後Phillipsが解析的な証明を与えたものです(Phillips, 1986)。Phillipsの論文は非常に難解ですが、興味のある人は是非読んでみてください。

ということでまとめると、沖本本p.127では見せかけの回帰を以下のように定義しています。

定義6.1（見せかけの回帰）　単位根過程 $y_t$ を定数と $y_t$ と関係のない単位根過程 $x_t$ に回帰すると、 $x_t$ と $y_t$ の間に有意な関係があり、回帰の説明力が高いように見える現象は見せかけの回帰(spurious regression)といわれる。

実際に見せかけの回帰において起きること

11年前はわざわざ実際の株価データを取ってきて試してみましたが、面倒なので最も単純な例をお見せしようと思います。まず、以下のプロットに見られるような目的変数yと説明変数V1-8からなるデータセットがあったとします。

で、これを重回帰分析にかけた結果を以下に示します。

#R> Call:
#R> lm(formula = y ~ ., data = x1)
#R> 
#R> Residuals:
#R>     Min      1Q  Median      3Q     Max 
#R> -5.4463 -1.1911  0.0043  1.3957  4.1310 
#R> 
#R> Coefficients:
#R>               Estimate Std. Error t value Pr(>|t|)    
#R> (Intercept)  0.3850376  1.1143297   0.346  0.73049    
#R> V1           0.4233928  0.1549390   2.733  0.00755 ** 
#R> V2           0.5437312  0.1108782   4.904 4.09e-06 ***
#R> V3          -0.2474450  0.1144308  -2.162  0.03321 *  
#R> V4           0.0314953  0.1403406   0.224  0.82293    
#R> V5           0.1435052  0.1153847   1.244  0.21680    
#R> V6           0.2104091  0.0997703   2.109  0.03770 *  
#R> V7           0.0003173  0.1497070   0.002  0.99831    
#R> V8          -1.0901879  0.0764275 -14.264  < 2e-16 ***
#R> ---
#R> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#R> 
#R> Residual standard error: 2.002 on 91 degrees of freedom
#R> Multiple R-squared:  0.9526,	Adjusted R-squared:  0.9485 
#R> F-statistic: 228.7 on 8 and 91 DF,  p-value: < 2.2e-16

5個の変数が有意となっただけでなく、調整済み決定係数も0.9を超える上にモデル全体のp値も極めて低いと、言うことなしといった感じです。ところが、このデータセットを一階差分（今の時点の値から前の時点の値を差し引いた「増減分」）した上で、プロットしてみるとこうなります。

もう見るからに雲行きが怪しいですね。実際に重回帰分析にかけてみると、こうなります。

#R> Call:
#R> lm(formula = y ~ ., data = x2)
#R> 
#R> Residuals:
#R>      Min       1Q   Median       3Q      Max 
#R> -2.30262 -0.57999 -0.00243  0.58873  2.79803 
#R> 
#R> Coefficients:
#R>              Estimate Std. Error t value Pr(>|t|)  
#R> (Intercept)  0.223138   0.101767   2.193   0.0309 *
#R> V1           0.031278   0.114305   0.274   0.7850  
#R> V2          -0.012966   0.101775  -0.127   0.8989  
#R> V3          -0.030722   0.102568  -0.300   0.7652  
#R> V4           0.077622   0.097733   0.794   0.4292  
#R> V5          -0.053583   0.099928  -0.536   0.5931  
#R> V6           0.006763   0.114778   0.059   0.9531  
#R> V7           0.039528   0.101882   0.388   0.6989  
#R> V8          -0.113253   0.099598  -1.137   0.2585  
#R> ---
#R> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#R> 
#R> Residual standard error: 0.9799 on 90 degrees of freedom
#R> Multiple R-squared:  0.03118,	Adjusted R-squared:  -0.05494 
#R> F-statistic: 0.362 on 8 and 90 DF,  p-value: 0.9378

豈に図らんや、今度は切片以外の全ての変数のp値が高くなってしまい、有意なものは一つもなし。モデル全体を見ても調整済み決定係数はほぼ0であり、p値も0.9以上と「何も説明できていない」に等しい結果となりました。

そう、実はこれらの変数は元々は2番目のセットアップの方が真であり、本来なら全ての変数はただのホワイトノイズ。1番目の方は全ての変数を累積値にただ変換して、全てランダムウォークに変えただけのものなのです。にもかかわらず、重回帰分析の結果はまるっきり異なります。これこそが、「見せかけの回帰」現象なのです。

見せかけの回帰への対処法

ということで、単位根過程のようなトレンド的なうねりのある時系列データ同士で回帰分析を行うと、とんでもないことになることが分かりました。では、そのようなデータセットに対して回帰分析を行いたい場合は、どう対処すれば良いのでしょうか？　そのアプローチについては11年前の記事でも取り上げましたが、近年普及してきた手法も幾つかありますのでここでまとめて挙げておきます。

なお、本来ならば全ての系列（変数）に対して単位根検定を行って白黒付けるべきではありますが、特にマーケティング分野では大抵の時系列データは単位根過程なので、いきなり下記の対処法をとっても大きな問題になることはあまりないと思われます。

差分系列に変換する

これは先述したシミュレーションのやり方そのままで、単に「今の時点の値から前の時点の値を差し引いて『増減分』（差分）を得る」というものです。基本的には、これさえやれば見せかけの回帰は確実に防ぐことができます。ただし、このやり方だと個々の変数の情報量が落ちてしまうことと*3、長期トレンドへの対処という意味では不完全な気がしています。

VARモデルを使う

これも前掲の11年前の記事でも取り上げていたやり方で、要は「漫然と回帰分析するのではなく時系列回帰手法であるVARを使え」というものです。勿論これも有効なアプローチではありますが、外生変数が多いケースでは使いづらい*4のでマーケティング分野ではあまり一般的ではないように見受けられます。

動的線形（状態空間）モデルやベイズ構造時系列モデルを使う

ということで、近年マーケティング分野で広く用いられているのがこちらです。即ち「目的変数からトレンド成分含めて系列相関要素を分離する」「その上で回帰で説明できる部分に絞って目的変数にアプローチする」やり方です。このやり方であれば、説明変数の情報量を減らし過ぎずに、尚且つ目的変数に対する影響度を定量化することができます。冒頭のQuerieでの質疑にも挙がっていたMMM (Media/Marketing Mix Models)も、どの実装であれ基本的にはこのやり方をとっています。

ただ、このやり方であってもトレンドの扱いには苦慮させられることが多く*5、決して万能とは言い難い点には注意が必要です。そもそも本来の回帰分析においても不均一分散の扱いは厄介ですし、トレンドといえばその不均一分散の代表みたいなものです。もう少し本質的な解決策があったら良いなと個人的にはいつも思うのですが、そんなものがあったら立派な計量経済学の論文が書けてしまうと思うと、まぁ難しいのかなと……ということで、お後がよろしいようで。

Rコード

老害も良いところの僕にはもはやモダンなRコードを書くのは無理ということで、カビが生えた感の酷いベタベタしたクソコードである点何卒ご容赦ください（泣）。

x1 <- c()

for (i in 1:9){
  set.seed(500 + i)
  tmp <- cumsum(rnorm(100, 0, 1))
  x1 <- cbind(x1, tmp)
}

x1 <- as.data.frame(x1)
nm <- c()

for (i in 1:9){
  nm <- c(nm, paste0('V', i - 1))
}

nm[1] <- 'y'
names(x1) <- nm
summary(lm(y ~ ., x1))

x2 <- diff(as.matrix(x1))
x2 <- as.data.frame(x2)
summary(lm(y ~ ., x2))

plot.ts(x1)
plot.ts(x2)