こういうメタ分析系の記事を書く時というのは大抵ネタ切れの時なんですが（汗）、最近になってこの辺のポイントでつまずいて困っているビジネスデータ分析の現場の話を聞くことがまた増えてきたので自分向けの備忘録も兼ねて記事としてまとめておきます。
そうそう、時系列分析の話って厳密にやり始めるとキリがないので、例えば単位根過程まわり（特に共和分のあたりを含めた複数時系列間の関係性の話とか）は「トレンドに注意せよ」という大きなくくりにまとめて、厳密な議論は割愛して出来る限り実務面で押さえるべきポイントに絞ろうと思います*1。悪しからずご了承あれ。
周期性のあるデータには真っ先に季節調整を
ビジネス時系列データは例えば毎日毎時の売上高とか契約数とかコンバージョン数とか、どこからどう見ても曜日変動とか24時間変動などの周期性が乗っているデータであることが多いです。にもかかわらず、その手の周期性に何の処理もせずにそのまま相関係数を計算したり回帰分析しちゃったりするケースがまま見受けられます。しかし、それをやってしまうのは非常に危険で誤った結論に陥る可能性があります。
例えば、上の図の2つの時系列。1行目に示しているのはシミュレーションデータですが、どちらも7日周期の曜日変動を模した季節成分が乗せてあります。そこで黒い時系列(x)で赤い時系列(y)を回帰すると、以下のようになります。
> summary(lm(y~x)) Call: lm(formula = y ~ x) Residuals: Min 1Q Median 3Q Max -3862 -2308 -1348 3968 6291 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 8722.004 1148.984 7.591 1.49e-11 *** x 21.105 3.183 6.631 1.57e-09 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 3210 on 103 degrees of freedom Multiple R-squared: 0.2992, Adjusted R-squared: 0.2924 F-statistic: 43.97 on 1 and 103 DF, p-value: 1.568e-09
めちゃくちゃばっちり回帰関係が見られることが分かります。しかし、実はこのデータはstl関数で7日周期の季節調整をかけると、2行目のような周期成分が検出され、残差として3行目のような時系列だけが残ります。見た目からしていかにも定常過程、すなわちある平均を中心として単にランダムに上下動しているだけの時系列っぽいですね。そこで黒い残差時系列(xn)で赤い残差時系列(yn)を回帰すると、以下のようになります。
> summary(lm(yn~xn)) Call: lm(formula = yn ~ xn) Residuals: Min 1Q Median 3Q Max -821.05 -226.83 -13.55 246.47 671.18 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 40.182 29.107 1.38 0.170 xn -0.167 1.045 -0.16 0.873 Residual standard error: 297.9 on 103 degrees of freedom Multiple R-squared: 0.0002478, Adjusted R-squared: -0.009459 F-statistic: 0.02553 on 1 and 103 DF, p-value: 0.8734
当然のように、全く何の回帰関係も見られません（笑）。というように、単に周期成分が互いに乗ってしまっただけでもこれを回帰すると途端に立派な「統計的に有意な」関係性が出てきてしまいます。もちろん世の中には「その曜日変動自体に興味があるんだ！」という人もいるかと思いますが*2、曜日変動自体に関心がないにもかかわらずこのような分析をしてしまうと、誤った結論に至る危険性があるということは知っておいて損はないと思います。
トレンドを伴うデータはまずトレンドを何かしら処理してから
このブログの「時系列分析」カテゴリ記事では、散々「トレンドのある時系列を扱う時は気を付けろ」と口を酸っぱくして言ってきました。その理由は2点あります。以下簡潔に挙げておきます。
そもそもトレンドを伴う非線形なデータに対して回帰モデルを推定するのは難しい
この過去記事で紹介したデータセットですが、記事中ではStanを用いて二階差分トレンドを分離するという手法でトレンドを処理しています。ですが、例えばそれをせずに漫然と線形回帰すると以下のようになります。
> d <- read.csv('stan_trend_issue.txt', sep='\t') > d.lm <- lm(cv~., d) > matplot(cbind(d$cv, predict(d.lm, newdata=d[,-1])), type='l', lwd=c(2,3), lty=1, ylab='') > legend('topleft', c('Data', 'Fitted'), lty=1, lwd=c(2,3), col=c(1,2), cex=1.5)
学習データに直接fitさせた結果を当てはめているにもかかわらず、全く合っていません。リンク先記事中にきちんとトレンドを分離した場合の結果がありますのでここでは割愛しますが、単純に長期のトレンドが目的変数に乗っているだけでもうまくいかないというケースが結構多いというお話です。
見せかけの回帰が起きることがある
また、トレンドのある時系列同士で漫然と回帰すると思いっきり「見せかけの回帰」(spurious regression)現象を起こします。嫌な言い方をすると、ランダムに選んできた株価の時系列をただのランダムウォーク（ホワイトノイズの累積和時系列）で回帰してもほぼ必ず有意なモデルが出来てしまう、というものです。
リンク先の過去記事でも色々対処法を書いておきましたが、最も簡単なのは「差分値を取る」です。これだけでも見え方は大きく変わってきます。
> set.seed(1001) > s <- cumsum(rnorm(1000,0,1)) > set.seed(1002) > t <- cumsum(rnorm(1000,5,2)) > par(mfrow=c(2,1)) > matplot(cbind(s,t), type='l', lty=1, ylab='Cumulative') > matplot(cbind(diff(s),diff(t)), type='l', lty=1, ylab='Differential') # 原系列のままで回帰する > summary(lm(t~s)) Call: lm(formula = t ~ s) Residuals: Min 1Q Median 3Q Max -2566.4 -1431.9 129.8 1231.8 1987.5 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2495.581 44.745 55.774 < 2e-16 *** s 34.540 6.367 5.425 7.28e-08 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1411 on 998 degrees of freedom Multiple R-squared: 0.02864, Adjusted R-squared: 0.02767 F-statistic: 29.43 on 1 and 998 DF, p-value: 7.278e-08 # 差分系列で回帰する > summary(lm(diff(t)~diff(s))) Call: lm(formula = diff(t) ~ diff(s)) Residuals: Min 1Q Median 3Q Max -6.3350 -1.2688 -0.0041 1.3234 5.9260 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 5.02499 0.06082 82.615 <2e-16 *** diff(s) 0.03900 0.06297 0.619 0.536 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.922 on 997 degrees of freedom Multiple R-squared: 0.0003846, Adjusted R-squared: -0.000618 F-statistic: 0.3836 on 1 and 997 DF, p-value: 0.5358
原系列で回帰すると非常に低いp値（偏回帰係数もモデル全体も）が出てきてしまうのに対して、差分系列で回帰すると物の見事に「無関係」としか言いようのない結果が返ってきます。このデータは最初に生成したRコードを見て一目で分かるように、実はランダムウォーク同士です。それでもこういう結果になるということで、注意が必要です。
ただし短期的な上下動が激しいnoisyなデータではむしろトレンドを見た方が良いことも
一方で、実は短期の上下動が全てただのノイズでトレンドこそが重要という面倒臭いケースもあります。これについては前に寓話的な記事を書いたことがあるので、もしかしたらご記憶の方もいらっしゃるかもしれません。
こういう場合は、いわゆる「平均への回帰」的な「昨日はたまたま上がったが今日はたまたま下がった」というようなランダムな変化の繰り返しが、あたかも自分たちのアクションに応じて上下動しているかのように見えてしまうことがあり、実際に回帰モデルを組んでみたら定常部分はアクションを表す説明変数とは何の関係もありませんでしたことになることが意外とあったりするものです。
ということで、色々な側面から評価してみても短期的な上下動が何によってもたらされているかがよく分からないようなケースでは、むしろトレンドの方を重視した方が良いのではないかと個人的には考えています。単にトレンドだけ分類するなら、色々な方法論があります。RならStanでやるまでもなく、stl関数などで簡単に季節調整をかけた上で分離できます。
ちなみにこれは極端な例ですが、わざと説明変数（3個）とは全く関係のないホワイトノイズに二階差分トレンドと適当な季節調整を加えた時系列を生成させ、これをベイジアンモデリングで推定して当てはめた結果です。当たり前ですが、トレンド以外は何も合っていません（笑）。でもこういう場面で無理くり回帰モデル組みたがる人も世の中多いように思います。
こちらから介入操作ができない時系列同士でモデリングする際は、おそらくベクトル自己回帰系のモデルの方が良い
株価同士とまでは言わないにせよ、例えばwebマーケティング界隈だと「オークション型広告の日次の投下額」と「日次のコンバージョン数」との関係を知りたいというようなケースがままあったりします。こういう場合、前者はある一定期間内の総額設定ぐらいしか介入操作ができない（＝日次の投下額そのものはコントロールできない）ので、そのまま説明変数（操作変数）としてモデリングに使うのは難しいことが多いです。
そういう時のために便利なのが、ベクトル自己回帰(VAR)系のモデルです。詳細は上記過去記事や沖本本やHamiltonやその他計量時系列のテキストに譲りますが、基本的には「相互に自己回帰する」タイプのモデルで尚且つGranger因果検定・インパルス応答関数推定・予測誤差分散分解といった時系列系因果推論に使えるというメリットがあります。これなら、例えば「オークション型広告を出稿すればk日後にコンバージョン数が増加する」と言った関係性を明らかにすることもできます。
また、前述の「トレンドがある時系列同士」の問題の一つである「見せかけの回帰」を避けるという点でも、VARモデルは有効です。そもそもVARモデルはラグ値（差分値）を説明変数に用いるという点で見せかけの回帰回避のための選択肢の一つであり*3、また大抵のVARモデルの推定メソッドにはトレンド処理のオプションがついており、さらに見せかけの回帰が起きるようなトレンド時系列同士であってもVARモデルのインパルス応答関数推定が因果推論に使える*4というメリットがあります。
ただし、トレンド時系列同士でVARモデルを組む場合は共和分に注意が必要です。詳細は上記記事や沖本本などに譲りますが、Johansenの手順などで共和分ランクを決定し、共和分関係が認められる場合はベクトル誤差修正モデル(VECM)に転換する必要があります。
外生変数を説明変数（操作変数）としてモデリングする際は、多重共線性やランク落ちやその他説明変数の偏りに要注意
これはどちらかというと単純な回帰モデル自体の注意事項でしかないのですが、時系列データのモデリング時に手に入る説明変数（操作変数）的なものの多くはその辺全く勘案されていないものだったりするので、事前に目を皿のようにしてチェックした方が良いです（笑）。
結構困るのが「ある一定期間内に限って全ての投下施策（広告など）を集中させていて、残りの期間はほとんど何もやっていないor一部の恒常的施策だけ延々とやっている」みたいなケース。データとして表現するとこんな感じになりますかね。
こういうデータで回帰モデルを推定してみると、得てして「全期間で実施されている施策」のパラメータだけが異様に大きくなり、「一定期間内に集中して実施されている施策たち」のパラメータは逆に異様に小さくなってしまいがちです。色々理由は考えられるのですが、こういうケースで真っ先に疑うべきはやはり多重共線性(multicolinearity)かと。実際、このケースで{usdm}パッケージのvif関数を使ってVIF統計量を計算してみると、1つ目の説明変数以外は60とか80とか高い数字を叩き出していて、そりゃダメだろうなぁという気になります。
もちろん、多重共線性に気を付ける以前の問題としてランク落ちを起こしているとかいう論外なケースもあったりします*5。実験計画法のような概念を全く持たずに漫然と時系列データを測定しているだけのところだとありがちな話です。
介入効果を見たかったら差分の差分法(Difference-in-Differences)で
時系列データであっても、というか時系列データだからこそ「○月○日にAという施策（介入）を行ったのでそれがKPI（売上高など）に本当にプラスになったかどうか知りたい」というような問題が持ち上がることがあったりします。
これは色々な考え方がありますが、素朴に統計的因果推論の立場を踏襲するならやはり差分の差分法(Difference-in-Differences)の枠組みに則ってやるのが一番確実でしょう。詳細については上記の記事をご覧いただくのが手っ取り早いと思います。どうやるのかをまとめて書いておくと、
- 「ある時点での介入あり」「介入なし」の2つの時系列データを用意する
- 後者の時系列データを活用して、前者で「介入がなかった場合」の反実仮想(counterfactural)データを生成する
- 反実仮想データと実際の前者の時系列データとを比較し、有意な差があれば介入効果ありと結論付ける
というものです。クロスセクションデータでまとめて比較する場合のやり方は上記の記事通りで良いのですが、時系列データでやる場合は{CausalImpact}パッケージを使うのが便利です。
> set.seed(71) > x <- rnorm(100, 1000, 50) > set.seed(72) > y <- rnorm(100, 850, 75) > set.seed(73) > inc <- rnorm(50, 100, 60) > y1 <- y + c(rep(0, 50), inc) # 11.8%の上乗せ効果を仮定した > d <- cbind(y1, x) > library(CausalImpact) > impact <- CausalImpact(d, c(1L, 50L), c(51L, 100L)) > summary(impact) Posterior inference {CausalImpact} Average Cumulative Actual 969 48456 Prediction (s.d.) 849 (13) 42470 (647) 95% CI [824, 874] [41217, 43697] Absolute effect (s.d.) 120 (13) 5986 (647) 95% CI [95, 145] [4759, 7239] Relative effect (s.d.) 14% (1.5%) 14% (1.5%) # 95%CIにちゃんと入って推定されている 95% CI [11%, 17%] [11%, 17%] Posterior tail-area probability p: 0.00103 Posterior prob. of a causal effect: 99.89701% For more details, type: summary(impact, "report") > plot(impact)
というように、2つの時系列データ同士で尚且つ差分の差分法デザインに沿ったものであれば、{CausalImpact}で適切に介入による上乗せ効果が検定・推定可能だということが分かります。
時系列データといえども外生変数メインでモデリングするなら交差検証を忘れずに
これらの問題点を全てクリアして、時系列モデルの推定がうまくいったとしましょう。特に介入操作に基づく外生変数がメインのケースでは、割と上手くいくこともそれなりに少なくないというのが個人的な経験に基づく印象です。けれども、そこで最後に待ち受ける問題があります。それは「汎化性能」(generalization)です。いかに優れたモデルであっても、未来値への当てはまりが悪ければ「過学習」(overfitting)を起こしている可能性があり、妥当なモデルとは言えません。
考え方は色々ありますが、個人的には時系列モデリングはトレンドなどのせいで適切な情報量規準が使えないこともままあるのでベタですが交差検証(cross validation: CV)をやるようにしています。時間がない時はholdout CVしかやらないことが多いですが、時間があれば上記記事のようなsliding windowによるCVを回すこともあります。
モデリングというと、時と場合とテーマによっては「パラメータ（偏回帰係数）の大小にしか興味がない」場合もあったりしますが、そういう場合であっても僕は可能な限りCVをやるようにしています。何故なら、「そのパラメータの大小に基づいて意思決定して何かをした場合、仮に過学習したモデルに基づいていたらやっぱり誤った結果になる可能性が高い」ためです。何もモデルそのもので直接未来値の予測をやらない場合であっても、モデルが提示したパラメータの大小に基づいて介入操作を行うのであれば、それは間接的にはモデルから予測値を算出しているのと同じことになるわけで、ならばやはり汎化性能は重視されるべきだと思う次第です。
最後に
ほとんどこのブログの過去記事を引用してばっかりの記事になってしまいましたが（笑）、たまにはこうやってまとめて自分でも復習した方が良いかなということで備忘録代わりに置いておきます。。。お粗末様でした。
