売上げを気温で回帰してみる

R 統計科学

あたまだし

掲題の件、最近出たわかりやすいデータ分析の入門書、

において

  • 売上げ v.s. 気温

として回帰分析をしていたケースがあった。これが時系列データに対するやってはいけない回帰、すなわち、みせかけの回帰に相当するケースではないか？という話を久保先生としていた。


ここでは、これを簡単なモデルのシミュレーションを通じてチェックしてみたい。

売上げ(平均値が正の正規分布)を気温(平均回帰過程)で回帰してみる

簡単のため、

としてモデリングする。コードはかつてないほどの殴り書き＆モデルパラメータは適当。
kに対するforループが上記の過程に従う１パスのモンテカルロシミュレーションに対応し、それを何回も繰り返して、回帰して、気温に対する回帰係数のp値の分布を見てみる。

size_step <- 100
sigma <- 5
dt <- 0.1
theta <- 1
mu <- 25
p_value <- numeric(10^3)
for(k in 1:10^3)
{
  temp <- rep(mu, size_step+1)
  for(i in seq_len(size_step))
  {
    temp[i+1] <- temp[i] - theta*(temp[i] - mu)*dt + sigma*sqrt(dt)*rnorm(1)
  }
  
  sale <- rnorm(size_step+1, mean=20, sd=8)
  summary(lm(sale~temp))
  p_value[k] <- coef(summary(lm(sale~temp)))[2,4]
}
hist(p_value)


ほぼ一様な分布だ。さらに、 p値が5%を下回る割合を計算すると

> sum(p_value<0.05)/10^3
[1] 0.05

となったので、p値の割合(信頼度としての意味合いとしてのp値) としてもよさそうだ。

経路の可視化

ちなみに途中の1パスを可視化すると、以下のような感じ。

#plot
library("ggplot2")
data.frame(
  date=rep(Sys.Date()+seq_len(size_step+1), 2), 
  group=rep(c("temp", "sale"), each=size_step+1), 
  value=c(temp, sale)
)　%>%
ggplot(aes(x=date, y=value, group=group, color=group)) + geom_line(size=2)

売上げ(平均回帰過程)を気温(平均回帰過程)で回帰してみる

業種によってはそうなのかもしれないが、売り上げ自体が平均回帰的な会社もあるかもしれない。そんなケースもあるはずだと思って同様の計算をしてみる。

size_step <- 100
sigma <- 5
dt <- 0.1
theta <- 1
mu_x <- 25
mu_y <- 35
p_value <- numeric(10^3)
for(k in 1:10^3)
{
  x <- rep(mu_x, size_step+1)
  y <- rep(mu_y, size_step+1)
  for(i in seq_len(size_step))
  {
    x[i+1] <- x[i] - theta*(x[i] - mu_x)*dt + sigma*sqrt(dt)*rnorm(1)
    y[i+1] <- y[i] - theta*(y[i] - mu_y)*dt + sigma*sqrt(dt)*rnorm(1)
  }
  
  p_value[k] <- coef(summary(lm(y~x)))[2,4]
}
hist(p_value, 100)

・・・明らかに有意となるケースが量産されている、これはやっちゃだめなやつだ！！！

> sum(p_value<0.05)/10^3
[1] 0.507
経路の可視化


見せかけの回帰の例

さらにちなみに、ダメなケース（ランダムウォーク同士の回帰、みせかけの回帰）だと、シードを固定しなくても、２〜３回も下記コードをまわせば即"みせかけの回帰"として有意な結果(xの回帰係数が有意)を作ることができる。

x <- cumsum(rnorm(100))
y <- cumsum(rnorm(100))

例えば以下で、xの係数が有意となっているが、これは明らかにみせかけの回帰の例だ。

> summary(lm(y~x))

Call:
lm(formula = y ~ x)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.3241 -2.2174  0.0802  2.5778  5.8766 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.53328    0.67048   6.761 9.94e-10 ***
x           -0.28741    0.05257  -5.467 3.48e-07 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.052 on 98 degrees of freedom
Multiple R-squared:  0.2337,	Adjusted R-squared:  0.2259 
F-statistic: 29.89 on 1 and 98 DF,  p-value: 3.484e-07

まとめ

解析的に計算している論文(Phillips(1986))だとランダムウォーク(正規分布を累積していくもの)同士の回帰がダメだと主張していたはずなので、それに該当しないケースはOK,あるいは個別に考えて対応となるのだろうか。

ここでは

  • 1: 正規分布 v.s. 平均回帰過程
  • 2: 平均回帰過程 v.s. 平均回帰過程

で回帰したが、2のケースは明らかにみせかけの回帰であると疑ってもよさそうだ。なので、過程(確率過程)に対するノイズが正規分布の累積(ランダムウォーク)になってる過程同士の回帰がまずそうだといえそうだ。

特に今回の売上のケースはよくよく考えると平均回帰過程でモデリングするのは筋が悪く、１のように単純に"ある想定される売り上げの水準からちょっとずれる"というような考え方をベースに持ってくるべきかなと*1。

もっと簡単にいうと、気温の過程は平均回帰過程であるとした上では、売り上げをどうモデリングするか”が肝っぽい。

  • １：売り上げ(今期) = (数値)×売り上げ(前期) + 正規分布　→　アウト
  • ２：売り上げ(今期) = ある値(平均的な売り上げ) + 正規分布　→　OK

ってかんじだ。１のケースが見せかけの回帰に相当するんで、イケイケなべんちゃー（売り上げが前期比で130%UP(上記,数値に相当)！）みたいな会社だと、１がよさそうですが、その会社さんだとみせかけの回帰出せると思ってて、「あのイケイケベンチャーの売上は・・・なんと気温に依存していた！」みたいな話になるんでねーかな、みせかけだけど。

*1:数式で考えると売上げに対するノイズ(バラツキ)が次の期の売上へと影響を与える形になるので筋が悪いなって