統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

Ad-Stock効果を推定しつつ回帰を回したい③

Marketing Mix Modelingにおける広告の累積効果の推定について、以下の記事を書いた。

ushi-goroshi.hatenablog.com
ushi-goroshi.hatenablog.com

本日はこの続きで、同じような条件で作成したデータに対してstanを用いたベイズ推定を実施したので、その内容を紹介する。
なお先日の記事でoptimによりAd-Stock効果の回帰係数および減衰率を十分に精度よく推定できそうなことはわかっているが、

  • 将来的なモデルの高度化
  • 多様な推定方法の確立(手持ちの武器を増やす)

ことを目的としている。

シミュレーションデータの作成 〜 stanによるフィッティングの実行

まずは事前準備として、必要なライブラリを読み込む。当然ながらstanを使うので、事前にインストールを済ませた上で{rstan}を読み込む。

## load libraries
library(dplyr)
library(tidyr)
library(tibble)

## load rstan
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

続いてシミュレーションデータ作成用の関数を用意する。これは先日の記事で紹介したものと大きな変更はないが、拡張性を考慮して変数の順序を入れ替え、また説明変数の作成に用いていたrnormの代わりにrgammaを使用している。

set.seed(123)
simulate_y <- function(pars) {
   n         <- pars[1] # num of observation
   mu        <- pars[2] # intercept
   var_e     <- pars[3] # residual variance
   beta_01   <- pars[4] # regression coefficient of X1 to be esitmated
   lambda_01 <- pars[5] # decay rate of Ad-Stock effect of X1
   beta_02   <- pars[6] # regression coefficient of X2 to be esitmated
   lambda_02 <- pars[7] # decay rate of Ad-Stock effect of X2
   
   X_01_raw <- rgamma(n, 3) * ifelse(runif(n) > 0.7, 0, 1)
   X_01_fil <- stats::filter(X_01_raw, lambda_01, "recursive")
   
   X_02_raw <- rgamma(n, 2) * ifelse(runif(n) > 0.8, 0, 1)
   X_02_fil <- stats::filter(X_02_raw, lambda_02, "recursive")
   
   error <- rnorm(n, 0, sqrt(var_e))
   
   y     <- mu + beta_01 * X_01_fil + beta_02 * X_02_fil + error
   dat <- data.frame(
      "Y"          = y,
      "X_01"       = X_01_raw,
      "X_02"       = X_02_raw,
      "X_01_Fil"   = X_01_fil,
      "X_02_Fil"   = X_02_fil,
      "Y_lag"      = dplyr::lag(y, 1),
      "True_Error" = error)
   return(dat)
}

init_par <- array(c(100, 5, 2, 0.5, 0.6, 0.8, 0.5))
dat_Ana  <- na.omit(simulate_y(init_par))

シミュレーションデータが作成できたので、stanを実行するための準備を実行する。stanに渡すためのパラメータを定義し、事前に準備したstanのスクリプトを用いてMCMCによるフィッティングを行う。
サンプルは10,000とし、そのうち最初の5,000サンプルをwarm-upとして破棄した。この辺りの設定は結果を見ながらの試行錯誤によるもので、明確な根拠に基づくものではない。traceプロットやヒストグラムを見る限り、もう少しサンプルを減らしても問題ないように思う。

dat_Ord <- list(N    = nrow(dat_Ana),
                Y    = dat_Ana$Y,
                X_01 = dat_Ana$X_01,
                X_02 = dat_Ana$X_02)

## fitting 
fit_01 <- stan(file = './AdStock_Simulation.stan', 
               data = dat_Ord, 
               iter = 10000,
               chains = 4,
               seed = 1234,
               warmup = 5000)

stanコードの記述

stan()で参照しているstanのコードは以下のように記述した。dataブロックおよびparameterブロックは大した工夫もないが、減衰率を示すlambda_X_01およびlambda_X_02はパラメータスペースを(0, 1)に制限している。

このコードで特筆すべきはtransformed parametersブロックで、ここでlambda_X_0xの値に基づきAd-Stock変数を作成している。Rであればstats::filterで簡単に作成できるのだが、今回の分析では減衰率も推定の対象であるためどうしてもstan内で記述する必要があった。ここで作成したX_01_filおよびX_02_filをmodelブロックで用いている。

data {
   int N;
   vector[N] Y;
   vector[N] X_01;
   vector[N] X_02;
}

parameters {
   real mu;
   real<lower=0> var_error;
   real beta_X_01;
   real<lower = 0, upper = 1> lambda_X_01;
   real beta_X_02;
   real<lower = 0, upper = 1> lambda_X_02;
}

transformed parameters {
   vector[N] X_01_fil;
   vector[N] X_02_fil;

   X_01_fil[1] = X_01[1];
   X_02_fil[1] = X_02[1];
   
   ## Create Ad-stock variable
   for(j in 2:N) {
      X_01_fil[j] = X_01[j] + X_01_fil[j-1] * lambda_X_01;
      X_02_fil[j] = X_02[j] + X_02_fil[j-1] * lambda_X_02;
   }
}

model {
   ## Sampling
   for(i in 1:(N)) {
      Y[i] ~ normal(mu + beta_X_01 * X_01_fil[i] + beta_X_02 * X_02_fil[i], 
      var_error);
   }
}

結果

stanを実行した結果は以下のように取得できる。

## Plot
pars <- c("mu", "var_error", "beta_X_01", "lambda_X_01", "beta_X_02",
          "lambda_X_02")
print(fit_01)
stan_trace(fit_01, pars = pars)
stan_hist(fit_01, pars = pars)
stan_dens(fit_01, separate_chains = T, pars = pars)
stan_ac(fit_01, separate_chains = T, pars = pars)

結果を見てみよう。まずは推定されたパラメータの事後分布から。

> print(fit_01)
Inference for Stan model: AdStock_Simulation.
4 chains, each with iter=10000; warmup=5000; thin=1; 
post-warmup draws per chain=5000, total post-warmup draws=20000.

               mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
mu             4.70    0.01 0.63   3.41   4.29   4.72   5.13   5.91 13079    1
var_error      1.45    0.00 0.11   1.26   1.37   1.44   1.52   1.68 17554    1
beta_X_01      0.50    0.00 0.08   0.35   0.44   0.50   0.55   0.65 15821    1
lambda_X_01    0.63    0.00 0.08   0.44   0.58   0.63   0.69   0.77 12174    1
beta_X_02      0.71    0.00 0.10   0.51   0.64   0.70   0.77   0.91 13692    1
lambda_X_02    0.56    0.00 0.07   0.41   0.52   0.57   0.61   0.69 14090    1
X_01_fil[1]    4.74    0.00 0.00   4.74   4.74   4.74   4.74   4.74     2    1
X_01_fil[2]    3.51    0.00 0.40   2.63   3.27   3.55   3.79   4.18 12174    1
X_01_fil[3]    4.94    0.00 0.52   3.87   4.59   4.96   5.30   5.91 12453    1
X_01_fil[4]    9.09    0.01 0.73   7.65   8.60   9.09   9.58  10.49 12477    1
…
省略
…

X_01_filは本来推定すべきパラメータではないが、stanでは値が決まっていない変数(確率変数と見なされるもの)は全てパラメータとして扱われるため、特に分布に関心はないものの結果として含まれている。

では実際に関心のあるパラメータに着目すると、muが設定した値からやや低く推定されているものの、50%ベイズ信頼区間には真の値が含まれている。またX_01の回帰係数および減衰率も同様に設定した値に近いものを得られている。
一方で残差分散およびX_02の回帰係数・減衰率については、大きく離れてはいないものの50%信頼区間からは外れてしまっている。特に残差分散は何度か試行を繰り返したものの、そのいずれにおいても過少に推定されてしまった。この原因についてはよくわからない。
なおX_02については、回帰係数および減衰率を高めに設定することでより収束しやすい傾向があった。

プロットにより作成されるグラフは、それぞれ以下の通り。

> stan_trace(fit_01, pars = pars)

f:id:ushi-goroshi:20170627001030p:plain

> stan_hist(fit_01, pars = pars)

f:id:ushi-goroshi:20170627001044p:plain

> stan_dens(fit_01, separate_chains = T, pars = pars)

f:id:ushi-goroshi:20170627001116p:plain

> stan_ac(fit_01, separate_chains = T, pars = pars)

f:id:ushi-goroshi:20170627001127p:plain

traceプロットを見ると安定した形となっており、Rhatの数値と合わせて考えても収束していると判断して良さそうだ。分布の形をchainごとで比較しても一部で歪んでいる様子は見当たらない。
なおヒストグラムを確認すると基本的には正規分布に近い形となっているが、減衰率については左に裾を引く形の分布となっている。

まとめ

上記のシミュレーションの結果をまとめると、

  1. stanを用いることでAd-Stock効果をそれなりに精度よく推定できる
  2. 減衰率のような非線形の効果を素直にモデルに取り込むことができ、また分布の形もKoyckに限定されることなく柔軟なモデリングが可能
  3. しかしパルダ・モデルやoptimによる推定と比較すると精度に劣り、また処理時間およびコーディングに要する時間も比較的多くなる

と言えそうだ。分析準備に要するコストは習熟度が上がれば削減されるだろうが、一回の分析に要する処理時間が相対的に大きいため試行錯誤する場合には向いていないと思われる。
一方でベイズモデリングの大きなメリットである「モデルの柔軟性」は非常に魅力的に感じ、Marketing Mix Modelingの高度化を目指すに当たって積極的に使っていきたい技術である。

上記を踏まえてこれまでの分析結果をまとめると、

  • Ad-Stockの分布としてKoyck型を仮定でき、残存の規模(減衰率)に関心はない場合 → パルダ・モデル
  • 残存の規模や分布に関心があり、試行錯誤が必要 → optimによる推定
  • より柔軟なモデリング(階層ベイズや状態空間モデルへの拡張)を実施したい → stanによるモデリング

といった使い分けが良さそうだ。なおstanによるモデリングは分析者のアイディアおよび対象領域に対する知見の深さにかなり依存すると思われるため、少なくともMMM構築の初期では全く適していないだろう。
逆に言えば、当該領域について十分な経験を持った分析者がいるのであれば、アイディアの限りを用いて色々なモデルを試すことが出来るため大きな価値を提供することが出来るようになるだろう。

最後に、この分析に限らずMCMCによるもの全般に言えることだが、やはりパラメータの分布という形で結果を見せることができるというのは大きな強みであると感じた。統計解析に馴染みのない相手であっても、「パラメータがバラつきを持ち、その範囲を分布として示す」というのは直感的にわかりやすいだろうと思う。頻度論的なアプローチであってもパラメータの(誤差)分布を標準誤差として示すことは出来るが、理解の容易さという点では比較にならないだろう。

唯一の懸念は、「パラメータが分布を持つ」という仮定に耐えうるか、ということだろうか。