はじめに

マーケティング・ミックス・モデリングに興味を持ったので、簡単に触れてみます。

メモとして残しておきます。

コードを動かした結果は少しづつ追加していきます。

基本的なモデル

マーケティング・ミックス・モデリングではいくつかの効果を考えるが、繰越効果をどのようにモデリングするかが主眼になってようである。

マーケティング・ミックス・モデルの話に入る前に基本的な時系列モデルについて触れておく。

分布ラグモデ/直接ラグモデル

過去時点の説明変数をモデルのそのまま加えるモデル。

$\displaystyle{ Y_t = \alpha + \beta_0 X_t + \beta_1 X_t + \beta_2 X_{t-2} + \dots + \beta_K X_{t-K} + \epsilon_t \\ = \alpha + \sum_{k=0}^K \beta_k X_{t-k} + \epsilon_t }$

𝑘次繰越直接に表現している。

ラグ演算子を利用して書き直すと次にようになる。

$\displaystyle{ Y_t = \alpha + \beta_0 \left(1+ \frac{\beta_1}{\beta_0}L+ \frac{\beta_2}{\beta_0}L^2+\dots+ \frac{\beta_K}{\beta_0}L^K \right) X_t+ \epsilon_t }$

MA過程の形になっている。

このモデルでは、ラグの大きさをどう指定するかや多重共線性の問題が残ってしまう。

幾何級数型分布ラグモデル

$\displaystyle{ Y_t = \alpha + \beta X_t + \lambda X_t + \lambda ^2 X_{t-2} + \dots + \epsilon_t }$

ここで $\lambda$ を繰越率 retention rate と呼ぶ。

効果は1時点ごとに $\lambda$ 倍に減少することを意味する。

マーケティング活動の繰越効果の表現方法として、古くからもっともよく用いられておりコイック・ラグモデルとも呼ばれる。

利点としてはパラメータ数が 1 増えるだけ最大ラグ数を決めなくてよいことが挙げられる。

コイックモデルは自己回帰モデルに変形することができる。

$\displaystyle{ Y_t = \alpha + \beta X_t + \beta \sum_{k=1}^\infty \lambda^k X_{t-k} + \epsilon_t }$

$\displaystyle{ -\lambda Y_{t-1} = -\lambda \alpha -\lambda \beta \sum_{k=1}^\infty \lambda^k X_{t-1-k} + \epsilon_t }$

$\displaystyle{ Y_t = (1-\lambda ) \alpha + \beta X_t + \lambda X_{t-1} + (\epsilon_t -\lambda \epsilon_t ) }$

ラグ演算子を用いて書き換えるとAR(1)過程の形になる。

効果の合計は次のようになる。

$\displaystyle{ \beta \sum_{k=0}^{\infty} \lambda^k = \frac{\beta}{1-\lambda} }$

自己回帰分布ラグモデル

$\displaystyle{ Y_t = \alpha + \sum_{i=1}^p \gamma_i Y_{t-1}+ \sum_{k=0}^q\lambda^k \beta X_{t-k} + \epsilon_t }$

$\gamma_i \dots , \gamma_p$ は繰越効果のうち、マーケティング活動の種類に依存しない成分であると解釈できる。

幾何級数型分布ラグモデルでは撹乱項の効果は繰り越されないが、自己回帰分布ラグモデルでは撹乱項の効果は繰り越される。

また、ラグ演算子を利用し書き直すとARMAXモデルとなる。

$\displaystyle{ Y_t = \frac{\alpha}{1-\gamma_1 L - \dots -\gamma_p L^p} + \beta_0 \frac{\alpha}{1-\gamma_1 L - \dots -\gamma_p L^p} X_t+ \frac{1}{1-\gamma_1 L - \dots -\gamma_p L^p} \epsilon_t }$

R package

以上のモデルを実装するには次のようなpackageで実施できるそうです。

dyn パッケージ:推定できるモデル己回帰分布ラグモデル
gets パッケージ:推定できるモデル己回帰分布ラグモデル
marima パッケージ:推定できるモデル己回帰分布ラグモデル
sysid パッケージ:推定できるモデル己回帰分布ラグモデル ARMA 撹乱
dse パッケージ:推定できるモデル己回帰分布ラグモデル ARMA 撹乱
dlsem パッケージ:推定できるモデル直接ラグモデル , 負の二項分布ラグモデル

Media Mix Modeling

続いてGoogleでのMMMの取り組みを紹介している論文を読みました。

GoogleではMarketing Mix ModelingをMedia Mix Modelingと読んでいるようです。

research.google

日本語でまとめているブログも多く存在するので、読みやすいと思います。

tjo.hatenablog.com

つまり、n番煎じです。

この論文では、繰越効果にあわせ形状効果について検討しています。

繰越効果の表現

繰越効果はAdstockによってモデル化されています。

$\displaystyle{ adstock(x [t-L+1,m ], \dots, x [ t,m ] ; w_m, L) = \frac{ \sum_{l=0}^{L-1} w_m(l) x [t-l,m ]}{ \sum_{l=l}^{L-1} w_m(l)} }$

L：メディア効果の長さ

P：メディア効果のピーク/遅延、最初の露出から何週間遅れているか

D：メディアチャネルの減衰/保持率、効果

$w$ は、効果がどのように繰り越され減衰していくのかを表現する係数です。

主に2種類の減衰の方法が提案されています。

幾何減衰

$\displaystyle{ w^g_m(l; \alpha_m) = \alpha^l_m }$

遅延減衰

$\displaystyle{ w^g_m(l; \alpha_m, \theta_m) = \alpha^{(l-\theta)^2}_m }$

形状効果

投資量から得られるリターンは線形に伸びることは少なく、投資量がまだ少量の場合はリターンが得られやすく、その後のリターンは低減していき一定量におつくような場合や投資量が少ない場合はリターン効率は悪くある量を超えてくるとリターン量が伸び、その後低減するような場合がある。

このような幾つか考えられる投資量とリターンの関係性を、薬学で使われているHill関数という関数を用いて表現する。

$\displaystyle{ Hill(x_{t, m}; \cal{K}_m, \cal{S}_m) = \frac{1}{1 + (x_{t, m} / \cal{K}_m)^{- \cal{S}_m}}, ~ x_{t, m} \geq 0 }$

ここで、 $x_{t,m}$ はある期間のあるメディアへの支出、 $\cal{S} _m$ は傾きの形状に関するパラメータ、 $\cal{K} _m$ は形状に関するパラメータで飽和点の半分の量を指定する。つまり $Hill(Km)=1/2$ となる点である。

f:id:saltcooky:20210818001100p:plain

グラフの形状を見ると3つのパラメータが異なったとしても同じような形状を持つため、真値の推定が適切にできない場合が多い。そのため、 $S=1$ に固定する場合がある。

$C$ 個の共変量 $z_{t, c}$ を含めてると推定するモデルは次のようになる。

$\displaystyle{ y_t = \tau + \sum_{m=1}^{M} \beta_m Hill(x_{t, m }^{*}; \mathcal{K}_m,\mathcal{S}_m) + \sum_{c=1}^{C} \gamma_{c}z_{t, c} + \epsilon_t }$

Media Mixと呼びながらミックスしていないじゃんと思いましたが、1つの時系列が複数のメディアの出稿の結果が現れているからかと思いました。

また、交互作用なども普通に組み込むことをやっていいようです。

Rで試してみる

GoogleのMMMとそれと比較するための形状効果をlog関数とした場合のモデルの二つを生成してみました。

利用したデータはこちらのマーケティングデータで、週レベルでの販売、メディアインプレッション、およびメディア支出の4年間（209週間）の記録です。

用いたカラムは以下です。 mdsp_auddig, mdsp_audtr, mdsp_vidtr, mdsp_viddig, mdsp_so, mdsp_on,

auddig：digital audio
audtr：radio
vidtr：TV
viddig： digital video,
so：social media
on：online display

平均値が1になるようにスケールして用います。

ad_timeseries_scaled_df <-
  ad_timeseries_df %>% 
  mutate(across(starts_with("mdsp_"), ~ ./mean(.)),
               across(starts_with("mdip_"), ~ ./mean(.))) %>%
  dplyr::select(mdsp_auddig, mdsp_audtr, mdsp_vidtr, mdsp_viddig, mdsp_so, mdsp_on)

4つの場合の結果を見てみました。

形状効果をlog関数+トータルの収益を表現する場合
形状効果をHill関数+トータルの収益を表現する場合
形状効果をlog関数とした場合+各メディアのインプレッションを表現する場合
形状効果をHill関数とした場合+各メディアのインプレッションを表現する場合

形状効果をlog関数+トータルの収益を表現する場合

繰越効果にはコイックモデルの考え方を用いて次のようにします。

そして、形状効果はlog関数を用いました。

Stanコードは次のようなものになりました。

なお、ベースはこちらに記載されているコードをになります。

ichi.pro

functions {
  real Adstock(vector t, row_vector weights) {
    return dot_product(t, weights) / sum(weights);
  }
}

data {
  int<lower=1> N;
  int<lower=1> T;
  vector[T] y;
  matrix[T,N] X;
  int<lower = 0> max_lag;
}

parameters {
  real<lower=0> V;
  vector<lower=0>[N] beta;
  real<lower=0> beta0;
  real<lower=0, upper=1> alpha[N];
  vector<lower=0,upper=ceil(max_lag/2)>[N] peak;
}

transformed parameters {
  row_vector[max_lag] lag_weights;
  matrix[T-max_lag+1, N] X_adstocked;
  real<lower = 0> cum_effect;
  
  for (t in 1:(T-max_lag+1)) {
     for (n in 1:N) {
       for (l in 1:max_lag) {
          lag_weights[max_lag-l+1] = alpha[n]^l;
       }
       cum_effect = Adstock(sub_col(X, t, n, max_lag), lag_weights);
       X_adstocked[t, n] = log(cum_effect);
     }
  } 
}

model {
  beta0 ~ inv_gamma(1, 0.1); 
  beta ~ inv_gamma(1, 0.1); 
  V ~ inv_gamma(1, 0.1); 
  alpha ~ beta(2,2);
  peak ~ inv_gamma(1, 1); 
  for(t in max_lag:T){
    y[t] ~ normal(dot_product(X_adstocked[t-max_lag+1], beta), sqrt(V));
  }
}

このコードをキックします。

data_list <- list(T = nrow(ad_timeseries_scaled_df), 
                  N = ncol(ad_timeseries_scaled_df),
                  y = log10(ad_timeseries_scaled_df$sales), 
                  X = ad_timeseries_scaled_df,
                  slope = 1,
                  max_lag = 14)

mmm_log_res <- stan(file = "model/media_mix_model_log.stan", 
                    data = data_list,
                    control = list(adapt_delta = 0.9, max_treedepth = 15),
                    iter = 4000,
                    warmup = 1000,
                    chains = 3)

結果を確認してみます。

mcmc_areas_ridges(mmm_log_res, regex_pars = "^beta\\[", prob = 0.95, prob_outer = 1, point_est = "mean")

digital audio

audtr：radio
vidtr：TV
viddig： digital video, 上からdigital audio, radio, TV, digital video, social media, online displayの係数です。

f:id:saltcooky:20210814155039p:plain

形状効果をHill関数+トータルの収益を表現する場合

functions {
  real Adstock(vector t, row_vector weights) {
    return dot_product(t, weights) / sum(weights);
  }  
  real Hill(real t, real ec, real slope) {
    return 1 / (1 + (t / ec)^(-slope));
  }
}

data {
  int<lower=1> N;
  int<lower=1> T;
  vector[T] y;
  matrix[T,N] X;
  int<lower = 0> max_lag;
}

parameters {
  real<lower=0, upper=100> V;
  vector<lower=0>[N] beta_hill;
  real<lower=0, upper=100> ec[N];
  real<lower=0, upper=100> mu0;
  real<lower=0> slope[N];
  real<lower=0, upper=max_lag-1> peak[N];
  real<lower=0, upper=1> alpha[N];

}

transformed parameters {
  row_vector<lower=0>[max_lag] lag_weights;
  matrix[T-max_lag, N] X_adstocked;
  real<lower = 0> cum_effect;
  
  for (t in 1:(T-max_lag)) {
     for (n in 1:N) {
       for (l in 1:max_lag) {
          lag_weights[max_lag-l+1] = pow(alpha[n], (l-1-peak[n])^2);
       }
       cum_effect = Adstock(sub_col(X, t, n, max_lag), lag_weights);
       X_adstocked[t, n] = Hill(cum_effect, ec[n], slope[n]);
     }
  }   
}

model {
  slope ~ gamma(3, 1);
  ec ~ inv_gamma(1, 0.01); 
  peak ~ uniform(0, max_lag - 1);
  beta_hill ~ inv_gamma(1, 0.01); 
  alpha ~ beta(10, 10);
  V ~ inv_gamma(1, 0.01); 
  mu0 ~ inv_gamma(1, 0.01); 
  for(t in 1:(T-max_lag)){
    y[t] ~ normal(mu0+dot_product(X_adstocked[t], beta_hill), sqrt(V));
  }
}

推定されたパラメータ

上からdigital audio, radio, TV, digital video, social media, online displayの係数です。

f:id:saltcooky:20210823230131p:plain