時変係数動的線形モデル＆Pystanを使った広告効果の推定をやってみた

#やったこと
架空のお店の売上高のトレンドの推定と、とある期間に行われた広告キャンペーンによってどれくらい売上高が変化したのかを、時変係数動的線形モデルとPystanを用いて推定してみました。もともとは[時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装](https://www.amazon.co.jp/%E6%99%82%E7%B3%BB%E5%88%97%E5%88%86%E6%9E%90%E3%81%A8%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93%E3%83%A2%E3%83%87%E3%83%AB%E3%81%AE%E5%9F%BA%E7%A4%8E-R%E3%81%A8Stan%E3%81%A7%E5%AD%A6%E3%81%B6%E7%90%86%E8%AB%96%E3%81%A8%E5%AE%9F%E8%A3%85-%E9%A6%AC%E5%A0%B4-%E7%9C%9F%E5%93%89/dp/4903814874?SubscriptionId=AKIAIM37F4M6SCT5W23Q&amp;amp;tag=lvdrfree-22&amp;amp;linkCode=xm2&amp;amp;camp=2025&amp;amp;creative=165953&amp;amp;creativeASIN=4903814874)に掲載されていた事例をPystanを使って解いてみようと思ったのがきっかけです。ちなみに、本の中ではRのKFASを使ったソリューションが載っています。
また、一部結果等が見にくいところがありますが、<https://github.com/hrkzz/Data_analysis_examples/blob/master/1_Analysing_ad_effects_by_dynamic_linear_model.ipynb>にnotebookを公開しているので、必要に応じてご覧ください。

#環境
Google Colabを使いました。環境構築の手間を考えると本当に便利ですよね。

#ライブラリのインポート
```python
#基本ライブラリ
import numpy as np
import pandas as pd

#図形描画ライブラリ
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('seaborn-darkgrid')

from matplotlib.pylab import rcParams

#統計ライブラリ
import pystan
```


#シミュレーションデータ

```python
#サンプルサイズ
N = 450

#乱数固定
np.random.seed(10)

#時間によって変化する広告効果
true_reg_coef = -np.log(np.arange(1,51))*2+10

#ランダムウォークする水準値（トレンド）
mu = np.cumsum(np.random.normal(loc=0, scale=0.5, size=N)) + 25

#状態方程式
x = mu + np.concatenate((np.repeat(0,200),true_reg_coef, np.repeat(0,200)))

#観測誤差
obs_error = np.random.normal(loc=0, scale=2, size=N)

#観測値（広告効果が入った売上データ）
sales_ad = x + obs_error

#説明変数（1なら広告あり）
dummy_ad = np.zeros(N)
dummy_ad[201:251] = 1

#データフレームで整理
data = pd.DataFrame(dict(sales_ad=sales_ad, dummy_ad=dummy_ad))
```

```python
data.head()
```

||sales_ad|	dummy_ad|
|:-:|:-:|:-:|
|0|	22.834287|	0.0|
|1|	26.212849|	0.0|
|2|	25.638334|	0.0|
|3|	26.333011|	0.0|
|4|	25.669880|	0.0|

##図示
```python
rcParams['figure.figsize'] = 15,4
plt.plot(data.sales_ad)
plt.title('Simulated data series')
plt.show()
```

![1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/2b33de9e-a8e9-cf46-f0a2-d2b70da4e44f.png)


```python
rcParams['figure.figsize'] = 15,4
plt.plot(data.dummy_ad)
plt.title('Period of ad')
plt.show()
```

![2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/cff7613f-c312-6365-f2cc-43764bc02d77.png)


```python
rcParams['figure.figsize'] = 15,4
plt.plot(true_reg_coef)
plt.title('Effects of ad')
plt.show()
```

![3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/130ae091-fc7a-5c73-5dbe-5de59e633be5.png)


#モデリング
状態空間モデルを使って推定を行っていきたいと思います。このモデルでは、回帰係数が時間よって変化する時変係数のモデリングが可能となっていて、時間とともにその効果が薄れていくことを表現することが可能です。<br>
状態空間モデルのローカルレベルモデルに、広告効果に関する時変の説明変数を付け加えた、以下のモデルを使っていきます。広告効果のダミー変数の係数の右下に$t$と付けていますが、これが時変係数を表しています。

```math
\begin{align}
  \beta_t &= \beta_{t-1} + a_t&a_t \sim Normal(0,\sigma^2_a) \\
  \mu_t &= \mu_{t-1} + w_t &w_t \sim Normal(0,\sigma^2_w)  \\
  y_t &= \mu_t + a_t*dummy\_ad_t+v_t &v_t \sim Normal(0,\sigma^2_v) 
\end{align}
```

これを別の書き方にすると、

```math
\begin{align}
  \beta_t &\sim Normal(\beta_{t-1},\sigma^2_a)  \\
  \mu_t &\sim Normal(\mu_{t-1},\sigma^2_w)  \\
  y_t &\sim  Normal(\mu_t +  a_t*dummy\_ad_t,\sigma^2_v) 
\end{align}
```

##Stanの実装
Stanの式は次のとおり。なお、Stanでは事前分布を指定しないと自動的に制限された範囲内で十分に幅の広い一様分布が使われるようです。

```python

Effect_ad_model = '''
data{
    int<lower=0> N;
    real<lower=0> sales_ad[N];
    real<lower=0, upper=1> dummy_ad[N];
}

parameters{
    real<lower=0> s_a;
    real<lower=0> s_w;
    real<lower=0> s_v;
    real beta_zero;
    real beta[N];
    real mu_zero;
    real mu[N];
}

model{
    beta[1] ~ normal(beta_zero, s_a);
    
    for(i in 2:N)
      beta[i] ~ normal(beta[i-1], s_a);
    
    mu[1] ~ normal(mu_zero, s_w);
    
    for(i in 2:N)
      mu[i] ~ normal(mu[i-1], s_w);
    
    for(i in 1:N)
      sales_ad[i] ~ normal(mu[i] + beta[i]*dummy_ad[i], s_v);
}
'''
```

上記のstanファイルを直接読み込んでもいいのですが、コンパイルと乱数発生を別にして、乱数発生を色々試したかったので、外からstanファイルを読み込むことにします。読み込み方はGoogle Colab流です。

```python
from google.colab import files
uploaded = files.upload()
```

Effect_ad_model.stan(n/a) - 475 bytes, last modified: 02/05/2019 - 100% done
 Saving Effect_ad_model.stan to Effect_ad_model.stan

まずはstanファイルのコンパイル。

```python
stanmodel = pystan.StanModel(file='Effect_ad_model.stan')
```
INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_f35367349717199a1e7e29d1b7953b86 NOW.

データをstanに優しい形式にして、実行。BUGSなどのほかのパッケージに比べてStanはMCMCの自己相関が低いらしく、thinは1でも問題ないことが多いらしいですが、今回はthinを2にしておきます。

```python
data_stan = dict(N=N, sales_ad=data.sales_ad, dummy_ad=data.dummy_ad)
fit = stanmodel.sampling(data=data_stan, iter=2000, warmup=500, chains=4, thin=2, seed=1, n_jobs=-1)
```

#推定結果の解釈
1行目はモデルのファイル名
2行目にはMCMCのchain数、iterationのステップ数、warmupのステップ数、thinningのステップ数が表示。今回はデフォルトの設定。MCMCのサンプルの長さは

```math
 \frac{chains*(iterations - warmup)}{thinning}
```
になるので、今回は3000となります。

meanはMCMCのサンプルの平均値で事後平均。se_meanはmeanの標準誤差。sdは標準偏差。2.5%～97.5%はそれぞれの分位点。n_effはStanが自己相関等あkら判断した実効的なMCMCサンプル数。著者によると100くらいが望ましいです。RhatはMCMCの収束の判断のひとつ。1.1以下だと収束とみなせます。
今回はn-effはだいたい十分で、Rhatを見ると、うまく収束しているようです。

```python
print(fit)
```

Inference for Stan model: anon_model_883d422f9f425c3ce85e7b9851134b2f.
4 chains, each with iter=2000; warmup=500; thin=2; 
post-warmup draws per chain=750, total post-warmup draws=3000.

||mean|se_mean|sd|2.5%|25%|    50%|    75%|  97.5%|n_eff|Rhat|
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|s_a|0.48|0.11|0.27|0.08|0.25|0.47|0.68|1.0|6|1.44
|s_w|0.45|7.2e-3|0.07|0.33|0.4|0.44|0.49|0.61|99|1.06
|s_v|1.87|2.0e-3|0.07|1.73|1.82|1.87|1.92|2.02|1313|1.0

残りの結果は長いので省略。

s_wの真の値が0.5、s_vの真の値が2なので、なかなかいい推定結果に。

##Traceplot
Traceplotでも収束具合を確認

```python
rcParams['figure.figsize'] = 15, 10
fit.plot()
plt.show()
```

![4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/fcbfb7ae-a6e8-e9dc-e8c7-ecf970cac7ed.png)


#推定結果の図示
```python
#結果を抽出
ms = fit.extract() 

#muの事後平均を算出
mu_mean = ms['mu'].mean(axis=0)

#広告効果の回帰係数の事後平均を算出
beta_mean = ms['beta'].mean(axis=0)

#広告効果を算出
ad_effects = beta_mean*data.dummy_ad

#95パーセンタイルを抽出
mu_5 = np.array(pd.DataFrame(ms['mu']).apply(lambda x: np.percentile(x, 5), axis=0))
mu_95 = np.array(pd.DataFrame(ms['mu']).apply(lambda x: np.percentile(x, 95), axis=0))

beta_5 = np.array(pd.DataFrame(ms['beta']).apply(lambda x: np.percentile(x, 5), axis=0))
beta_95 = np.array(pd.DataFrame(ms['beta']).apply(lambda x: np.percentile(x, 95), axis=0))
ad_effects5 = beta_5*data.dummy_ad
ad_effects95 = beta_95*data.dummy_ad

#x軸
X = data.index
```

##トレンドの推定

```python
rcParams['figure.figsize'] = 15, 5

plt.plot(X, data.sales_ad, label='observed')
plt.plot(X, mu, label='true trend', c='green')
plt.plot(X, mu_mean, label='predicted trend', c='red')
plt.fill_between(X, mu_5, mu_95, color='red', alpha=0.2)
plt.legend(loc='upper left', borderaxespad=0, fontsize=15)
```

![5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/ab413c09-fbc8-ba4e-7f5a-e5654159cedd.png)

広告期間のところだけ、観測値より低く推定が出来ています。

##広告効果の推定

```python
rcParams['figure.figsize'] = 15, 5

plt.plot(X, np.concatenate((np.repeat(0,200),true_reg_coef, np.repeat(0,200))), label='true effects', c='green')
plt.plot(X, ad_effects, label='predicted effects', c='red')
plt.fill_between(X, ad_effects5, ad_effects95, color='red', alpha=0.2)
plt.legend(loc='upper left', borderaxespad=0, fontsize=15)
```

![6.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/0b0d5d56-7ebb-4bc8-72c9-30e565017e33.png)


```python
#広告期間のみにフォーカス
rcParams['figure.figsize'] = 15, 5

plt.plot(X[201:250], np.concatenate((np.repeat(0,200),true_reg_coef, np.repeat(0,200)))[201:250], label='true effects', c='green')
plt.plot(X[201:250], ad_effects[201:250], label='predicted effects', c='red')
plt.fill_between(X[201:250], ad_effects5[201:250], ad_effects95[201:250], color='red', alpha=0.2)
plt.legend(loc='upper right', borderaxespad=0, fontsize=15)
```

![7.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/1d90e4cf-11ce-7dee-cba7-1e8f22e30137.png)

右肩下がりは示せていますが、事後平均は真の効果の半分くらいになってますね。

```python
rcParams['figure.figsize'] = 15, 5

plt.plot(X, data.sales_ad, label='observed')
plt.plot(X, mu_mean+ad_effects, label='predicted trend + predicted ad effects', c='magenta')
plt.plot(X, mu_mean, label='predicted trend', c='red')
#plt.fill_between(X, mu_5, mu_95, color='red', alpha=0.2)
plt.legend(loc='upper left', borderaxespad=0, fontsize=15)
```

![8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/246473/d471fd57-f3db-d325-ca83-ff4e41d27192.png)

トレンドと広告効果を含めるとうまく、表現できている気がします。

#おわりに
以上、走り書きのようにコードを書いて、stanを回して、という感じで作業をしたので、モデリングを中心に多々不備があると思います。こうやったらいいよという、ご意見があればコメントをお願いします！

#参考書
<div  class="amazon Default"><div  align="left" class="pictBox"><a  target="_blank" href="https://www.amazon.co.jp/%E6%99%82%E7%B3%BB%E5%88%97%E5%88%86%E6%9E%90%E3%81%A8%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93%E3%83%A2%E3%83%87%E3%83%AB%E3%81%AE%E5%9F%BA%E7%A4%8E-R%E3%81%A8Stan%E3%81%A7%E5%AD%A6%E3%81%B6%E7%90%86%E8%AB%96%E3%81%A8%E5%AE%9F%E8%A3%85-%E9%A6%AC%E5%A0%B4-%E7%9C%9F%E5%93%89/dp/4903814874?SubscriptionId=AKIAIM37F4M6SCT5W23Q&amp;tag=lvdrfree-22&amp;linkCode=xm2&amp;camp=2025&amp;creative=165953&amp;creativeASIN=4903814874"><img  class="pict" style="margin-right:10px" align="left" hspace="5" border="0" alt="時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装" src="https://images-fe.ssl-images-amazon.com/images/I/516f3jVjO2L._SL160_.jpg"></a></div><div  class="itemTitle"><a  target="_blank" href="https://www.amazon.co.jp/%E6%99%82%E7%B3%BB%E5%88%97%E5%88%86%E6%9E%90%E3%81%A8%E7%8A%B6%E6%85%8B%E7%A9%BA%E9%96%93%E3%83%A2%E3%83%87%E3%83%AB%E3%81%AE%E5%9F%BA%E7%A4%8E-R%E3%81%A8Stan%E3%81%A7%E5%AD%A6%E3%81%B6%E7%90%86%E8%AB%96%E3%81%A8%E5%AE%9F%E8%A3%85-%E9%A6%AC%E5%A0%B4-%E7%9C%9F%E5%93%89/dp/4903814874?SubscriptionId=AKIAIM37F4M6SCT5W23Q&amp;amp;tag=lvdrfree-22&amp;amp;linkCode=xm2&amp;amp;camp=2025&amp;amp;creative=165953&amp;amp;creativeASIN=4903814874">時系列分析と状態空間モデルの基礎: RとStanで学ぶ理論と実装 [単行本]</a></div><div  class="itemSubTxt">馬場 真哉</div><div  class="itemSubTxt">プレアデス出版</div><div  class="itemSubTxt">2018-02-14</div></div><br  style="clear:left" clear="left"><br /><div  class="amazon Default"><div  align="left" class="pictBox">
<a  target="_blank" href="https://www.amazon.co.jp/Stan%E3%81%A8R%E3%81%A7%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E3%83%A2%E3%83%87%E3%83%AA%E3%83%B3%E3%82%B0-Wonderful-R-%E6%9D%BE%E6%B5%A6-%E5%81%A5%E5%A4%AA%E9%83%8E/dp/4320112423?SubscriptionId=AKIAIM37F4M6SCT5W23Q&amp;tag=lvdrfree-22&amp;linkCode=xm2&amp;camp=2025&amp;creative=165953&amp;creativeASIN=4320112423"><img  class="pict" style="margin-right:10px" align="left" hspace="5" border="0" alt="StanとRでベイズ統計モデリング (Wonderful R)" src="https://images-fe.ssl-images-amazon.com/images/I/416xdqDozvL._SL160_.jpg"></a></div><div  class="itemTitle"><a  target="_blank" href="https://www.amazon.co.jp/Stan%E3%81%A8R%E3%81%A7%E3%83%99%E3%82%A4%E3%82%BA%E7%B5%B1%E8%A8%88%E3%83%A2%E3%83%87%E3%83%AA%E3%83%B3%E3%82%B0-Wonderful-R-%E6%9D%BE%E6%B5%A6-%E5%81%A5%E5%A4%AA%E9%83%8E/dp/4320112423?SubscriptionId=AKIAIM37F4M6SCT5W23Q&amp;amp;tag=lvdrfree-22&amp;amp;linkCode=xm2&amp;amp;camp=2025&amp;amp;creative=165953&amp;amp;creativeASIN=4320112423">StanとRでベイズ統計モデリング (Wonderful R) [単行本]</a></div><div  class="itemSubTxt">松浦 健太郎</div><div  class="itemSubTxt">共立出版</div><div  class="itemSubTxt">2016-10-25</div></div><br  style="clear:left" clear="left">