何気なくR-Bloggerのタイムラインを見ていたら、"CausalImpact: A new open-source package for estimating causal effects in time series | Google Open Source Blog"という記事がシェアされていたので見に行ってみたのでした。これはもう読んで字の如く「GoogleがキャンペーンがKPIにもたらす因果的影響を時系列から推定する」ためのRパッケージの話題で、その名も{CausalImpact}という。
ということで、ちろっと触ってみたので簡単にレビューしてみようと思います。本当は色々試してみたかったんですが、ちょっと手元に良いデータがないのでヘルプの事例のみでご勘弁を。。。
インストール
追記 (Jan 29 2020)
現在はCRANからインストールできます。
install.packages("CausalImpact") library(CausalImpact)
GitHubに書いてありますが、{devtools}入れるところまで含めてたったのこれだけでおしまいです。
install.packages("devtools") library(devtools) devtools::install_github("google/CausalImpact") library(CausalImpact)
かなりあっさり入ります。が、依存パッケージがかなり多くて、特に{Rcpp}まわりでコケる人は多いかも。{Rcpp}を初めてインストールする人で何故か毎回コケてしまう人はこちらの過去記事を、gcc周り以外のエラーが出る人は{Rcpp}を一度削除してから再インストールしてみて下さい。{Rcpp}がきちんと入れば、{CausalImpact}もすんなり入るはずです。
コンセプト
GitHubの説明を見るよりもリリースブログ記事を見た方が早いと思います。元となったアイデアについては既に論文になっていて、以下のリンクから見られるようになっています。
要は「キャンペーンの影響を受けていると思しき目的変数の時系列データに対して、1) 影響を受けていない他の共変量となり得る時系列と、2) キャンペーン期間を表すベクトル、を用意して、キャンペーンの影響を受けなかったと仮定した場合の目的変数の予測値を状態空間＋MCMCで算出し、実測値と比較することでキャンペーンの因果的影響の大きさを見る」ということのようです。
多分ですが、論文中のFigure 2のグラフィカルモデルを見た方が話としては手っ取り早いかと。この辺は僕のDLMの勉強が進んでないので、そのうち勉強が済んだら後でまとめようかなと思ってます。。。
とりあえず{CausalImpact}を試してみる
そんなわけで、面倒なのでヘルプに従って例題をやってみます。ヘルプに載ってるのは以下のようなコード。
> set.seed(1) > x1 <- 100 + arima.sim(model = list(ar = 0.999), n = 100) # x1は共変量（説明変数） > y <- 1.2 * x1 + rnorm(100) > y[71:100] <- y[71:100] + 10 # yが目的変数 # 71番目以降に10だけキャンペーン効果が上乗せされている > data <- cbind(y, x1) > pre.period <- c(1, 70) # キャンペーン前 > post.period <- c(71, 100) # キャンペーン後 > impact <- CausalImpact(data, pre.period, post.period) # モデル計算を行い、オブジェクトに結果を収納する > summary(impact) Posterior inference {CausalImpact} Average Cumulative Actual 117 3511 Prediction (s.d.) 106 (0.43) 3195 (12.94) 95% CI [106, 107] [3169, 3219] Absolute effect (s.d.) 11 (0.43) 317 (12.94) 95% CI [9.7, 11] [292.4, 343] Relative effect (s.d.) 9.9% (0.41%) 9.9% (0.41%) 95% CI [9.2%, 11%] [9.2%, 11%] Posterior tail-area probability p: 0.00111 Posterior prob. of a causal effect: 99.88901% For more details, type: summary(impact, "report") > plot(impact)
まぁこんな感じで、キャンペーンの効果が真ん中の行に図示されていて、x=71以降でキャンペーンによるKPI増大が捉えられているのが分かるかと思います。
他の時系列データでやってみる
ちょっとイタズラして、こんな感じのデータを作って試してみました。イメージとしては、ソシャゲなどで毎週月曜初めとかに週1イベントがあってスパイクが入るような時系列データです。
> x1<-matrix(rnorm(300)) # ホワイトノイズ > x2<-matrix(rep(c(15,rep(0,49)),6)) # スパイク異常値が6回生じる時系列 > x<-x1+x2 # これらを足し合わせる > x[151:300]<-x[151:300]+5 # 後半だけ5の上乗せ分がキャンペーンによって生じる > impact2<-CausalImpact(cbind(x,x1),pre.period = c(1,150),post.period = c(151,300)) > plot(impact2)
ちゃんとキャンペーン効果があるという結論になっています。では、さらにイタズラしてこんな感じのデータではどうでしょうか？
> x1<-matrix(rnorm(300)) # 同じくホワイトノイズ > x2<-matrix(rep(c(15,rep(0,49)),6)) # 同じくスパイク異常値が6回生じる時系列 > x<-x1+x2 # 同じくこれらを足し合わせる > x[151:200]<-x[151:200]+5 # ごく一部の期間にだけ5の上乗せ分がキャンペーンによって生じる > impact3<-CausalImpact(cbind(x,x1),pre.period = c(1,150),post.period = c(151,300)) # 予測期間は同じにしておく > plot(impact3)
お、こうしてみるとcumulativeのところにキャンペーン効果が線形に乗らなかったことを示唆する感じの結果が見られますね（途中でサチって代わりに分散が増大している）。ということで、途中でキャンペーン効果が失われた場合でも何となく検出できそうな気がしてきました。
色々なツッコミ
これはちょっとDLMやってる人なら同じこと言いそうな感じですが、
@TJO_datasci 元論文、カルマンフィルタでいんじゃねーのこれモデルにしかみえねぇ・・・
— teramonagi (@teramonagi) 2014年9月14日
というのはまさにその通りで、そう言えば鉄板のナイル川流量データ（Rに入っている有名なサンプルデータセットで、1898年のアスワン・ダム建設で大幅に流量が変わる様子のモデリングをやるのがハイライト）でもカルマンフィルタで容易にやれますよねという。
あと、搦め手としては{MSwM}パッケージとかでマルコフ転換モデルでキャンペーン実施時期と一致するレジーム転換がEMアルゴリズムによって推定できればOK的な発想もアリかと。こんな感じです。
> set.seed(1) > x1 <- 100 + arima.sim(model = list(ar = 0.999), n = 100) > y <- 1.2 * x1 + rnorm(100) > y[71:100] <- y[71:100] + 10 > model1<-lm(y~x1) > msmModel1<-msmFit(model1,k = 2,sw = rep(T,3)) > plotProb(msmModel1,2)
この方法でも、やっぱりx=71以降でレジーム転換が起きているのが見て取れる＝キャンペーンの効果が生じているという判定ができるという。
まぁ、どう使っても良い気はしますが、例えばGranger因果なんかに比べると因果因果してるわけでもなく、DLMやカルマンフィルタよりも先進的というわけでもなく、でもキャンペーン期間が分かってる時にサクッと計算するのは簡単だなーとか、そんな感じでしょうか。ということで、お後がよろしいようで。。。
