統計コンサルの議事メモ

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

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

背景

しつこいようですが、Marketing Mix ModelingMMM)の話題です。

先日、こんな面白い論文を見つけました。 GoogleのResearcherによるMMMの論文(彼らはMedia Mix Modelingと呼んでいます)なのですが、ヒルの式を用いて広告のShape効果(Carveture効果)を推定するということをやっています。ここでShape効果・carveture効果とは、メディアの露出量に対する目的変数の反応を示す曲線を指すようで、ヒルの式とは:

$$ H(x; K, S) = \frac{1}{1 + (\frac{x}{K})^{-S}} $$

であり、$K > 0$や$S > 0$となるパラメータによってLogやSigmoidの形状を表現することができるようです。

ヒルの式によってxがどのような形状となるか、実際に確認してみましょう。まずはヒルの式を以下のように定義します。

hill_transformation <- function(x, k, s) {
   1 / (1 + (x / k)^-s)
}

続いて、パラメータとなるksをそれぞれ以下のように設定した場合をプロットしてみましょう。なおこの数値は論文から拝借しています。

x <- seq(0, 1, by = 0.05)

plot(x, hill_transformation(x, 0.4, 4), 
     type = "l", col = 2, xlim = c(0, 1), ylim = c(0, 1),
     ylab = "Hill Transformed Value")
par(new = T)
plot(x, hill_transformation(x, 0.4, 1), 
     type = "l", col = 3, xlim = c(0, 1), ylim = c(0, 1), axes = F,
     ylab = "")
par(new = T)
plot(x, hill_transformation(x, 0.8, 4), 
     type = "l", col = 4, xlim = c(0, 1), ylim = c(0, 1), axes = F,
     ylab = "")
par(new = T)
plot(x, hill_transformation(x, 0.8, 1), 
     type = "l", col = 5, xlim = c(0, 1), ylim = c(0, 1), axes = F,
     ylab = "")
legend("topleft", legend = c("k = 0.4, s = 4", "k = 0.4, s = 1", 
                             "k = 0.8, s = 4", "k = 0.8, s = 1"), 
       col = c(2:5), lty = rep(1, 4))

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

このように、パラメータを変更することで元の値を柔軟に変換することが可能です。これまで私が書いてきた記事では、いずれも広告効果(回帰係数)そのものを推定する方法にばかり注目しており、説明変数の"効き方"については触れてきませんでした。これではちょっと検討が足りないと言われても仕方ありません。

というわけで本エントリーでは、ヒルの式により変換したXを用いてシミュレーションデータを発生させ、元のXから変数変換のためのパラメータを推定できるかを検証したいと思います。なお元論文ではAd-Stock効果としてGeometric Ad-Stock:

$$ GA(x_{t}; \alpha, L) = \frac{\sum_{l = 0}^{L}x_{t-l}\alpha^{l}}{\sum_{l = 0}^{L}\alpha^{l}} $$

を用いていますが、今回の検証ではAd-Stock効果をみておらず、説明変数の生の値を変換しています1

ライブラリの読み込み

今回の分析でも{rstan}を使用します。

library(tidyverse)
library(ggplot2)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

シミュレーションデータの生成

分析用のデータ生成ですが、まず説明変数X1runifにより作成します。次に与えたパラメータからヒルの式を用いて変換し、目的変数yを作成します。

## データ生成用の関数を定義
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
   k_01      <- pars[6]  # k for X1
   s_01      <- pars[7]  # s for X1

   X_01_raw  <- runif(n)
   X_01_conv <- hill_transformation(X_01_raw, k_01, s_01)

   error <- rnorm(n, 0, sqrt(var_e))
   
   y     <- mu + beta_01 * X_01_conv + error
   dat <- data.frame(
      "Y" = y,
      "X_01"      = X_01_raw,
      "X_01_conv" = X_01_conv
   )
   return(na.omit(dat))
}

## データ生成
set.seed(123)
pars <- c(100, 5, 0.1, 0.8, 0.7, 0.4, 4)
dat <- simulate_y(pars)

X1について、

  • 推定したい広告効果 ⇒ 0.8
  • ヒルの式のkおよびs ⇒ 0.4および4

と指定しており、これらのパラメータを元のX_01から推定することが狙いです。ちなみに今回のパラメータではX_01変換後のX_01は以下のようになります:

plot(dat$X_01, dat$X_01_conv, col = 2, xlim = c(0, 1), ylim = c(0, 1),
     ylab = "")

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

optimによるフィッティング

手始めにoptimを使ってフィッティングしてみましょう。これまでの記事でも度々optimを使ってきましたが、なかなか精度良く推定することが可能です。ただしAICを計算するのが面倒なのでここではbetaを直接推定しないことにして、lmを使います。よって推定対象はk_01s_01の2つだけです。

return_AIC <- function(param, dat) {
   k_01      <- param[1]
   s_01      <- param[2]
   
   dat$X_01_conv <- hill_transformation(dat$X_01, k_01, s_01)
   AIC(lm(Y ~ X_01_conv, dat))
}
### 適当な数値を入れてみる
> return_AIC(rep(0.5, 2), dat)
[1] 60.47417

ちゃんと動くようなのでフィッティングしてみましょう。

param <- rep(1, 2)
res_optim <- optim(par = optim(par = param, fn = return_AIC, dat = dat)$par,
                   fn = return_AIC, dat = dat)

パラメータを確認してみると:

true_par    <- c(0.8, 0.4, 4)
dat$X_01_conv_opt <- hill_transformation(dat$X_01, 
                                         res_optim$par[1], res_optim$par[2])
optim_par     <- c(coef(lm(Y ~ X_01_conv_opt, dat))[2],
                   res_optim$par)
> print(cbind(true_par, optim_par), digits = 2)
              true_par optim_par
X_01_conv_opt      0.8       1.1
                   0.4       0.5
                   4.0       1.9

うーん、sがちょっと外れているようですね。推定された値でのプロットを見てみましょう。

plot(dat$X_01, hill_transformation(dat$X_01, 0.4, 4), 
     col = 2, xlim = c(0, 1), ylim = c(0, 1),
     ylab = "Hill Transformed Value")
par(new = T)
plot(dat$X_01, hill_transformation(dat$X_01, res_optim$par[1], res_optim$par[2]), 
     col = 3, xlim = c(0, 1), ylim = c(0, 1), axes = F,
     ylab = "")
legend("topleft", legend = c("k = 0.4, s = 4", "k = 0.49, s = 1.9"), 
       col = c(2:3), lty = rep(1, 2))

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

ちょっとカーブのメリハリがなく、全体として直線的になってしまっています。これが今回のサンプルで偶々発生したものなのか、それともsの推定に何らかの偏りがあるのか確かめてみましょう。同様の試行を1,000回繰り返してみます。

## 5番目のパラメータはAd-Stockのλだが今回は不使用
pars <- c(100, 5, 0.1, 0.8, 0.7, 0.4, 4)

n             <- 1000
res_optim_all <- matrix(NA, n, 3)
param         <- rep(1, 2)

try_s <- function(pars) {

   dat <- simulate_y(pars)
   res_optim <- optim(par = optim(par = param, fn = return_AIC, dat = dat)$par,
                      fn = return_AIC, dat = dat)
   dat$X_01_conv_opt <- hill_transformation(dat$X_01,
                                            res_optim$par[1], res_optim$par[2])
   return(c(coef(lm(Y ~ X_01_conv_opt, dat))[2], res_optim$par))
}

for (i in 1:n) {
   res_optim_all[i, ] <- tryCatch(try_s(pars), error = function(e) return(rep(0, 3)))
}

tryCatchを使って無理やりループを回しましたが、やはり推定が不安定なようで全体の4割ちょっとしか収束していません。その収束した解のバラツキを見てみても以下のようにかなり広がってしまっており、真値である4も中心とはなっていませんね。

> nrow(res_optim_all[res_optim_all[, 1] != 0.0, ])
[1] 434
MASS::truehist(res_optim_all[res_optim_all[, 1] != 0, 3])

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

Stanによるフィッティング

では今度はStanを使ってみましょう。以下のように指定します。

### データの再作成
set.seed(123)
pars <- c(100, 5, 0.1, 0.8, 0.7, 0.4, 4)
dat <- simulate_y(pars)

dat_Stan <- list(N        = nrow(dat),
                 Y        = dat$Y,
                 X_01     = dat$X_01
                 )

またStanのスクリプトは以下のようになります。今回は特別な指定はしていませんが、X_01_convを作成するためにtransformed_parametersブロックを置いています。

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

parameters {
   real mu;
   real beta_01;
   real<lower=0, upper=1> shape_k_01;
   real<lower=0, upper=5> shape_s_01;
   real<lower=0> var_error;
}

transformed parameters {
   vector[N] X_01_conv;
   for(k in 1:N) {
      X_01_conv[k] = 1 / (1 + (X_01[k] / shape_k_01)^-shape_s_01);
   }
}

model {
   // Sampling
   for(i in 1:N) {
      Y[i] ~ normal(mu + beta_01 * X_01_conv[i], var_error);
   }
}

上記のモデルを用いて、フィッティングしてみましょう。

fit_01 <- stan(file = '/YourDirectory/Hill_Transformation.stan',
               data = dat_Stan,
               iter = 10000,
               chains = 4,
               seed = 123)

結果の確認

フィッティングが終わったので、結果を見てみましょう。fit_01からサンプルを抽出して加工します。

## サンプルを抽出する
res_01 <- rstan::extract(fit_01)

## 該当するパラメータを取り出す
ests <- summary(fit_01)$summary
beta_rows  <- rownames(ests)[grep("beta", rownames(ests))]
shape_rows <- rownames(ests)[grep("shape", rownames(ests))]

まずは収束の判断です。Rhatが1.1未満であることを確認した上で、トレースプロットとヒストグラムを見てみましょう。

stan_trace(fit_01, pars = beta_rows)
stan_trace(fit_01, pars = shape_rows)
stan_hist(fit_01, pars = beta_rows)
stan_hist(fit_01, pars = shape_rows)

f:id:ushi-goroshi:20180325224627p:plain f:id:ushi-goroshi:20180325224703p:plain f:id:ushi-goroshi:20180325224757p:plain f:id:ushi-goroshi:20180325224841p:plain

うーん、kが全体的に広がっていたり、sは真値とかけ離れたところにピークがあったり、推定に失敗している様子ですね。 ひとまず回帰係数の推定結果を見てみましょう。実際の値と推定値を並べてみます。

beta_par <- 
   ests %>% 
   data.frame %>% 
   select(mean) %>% 
   mutate("Par" = rownames(ests)) %>% 
   filter(Par %in% beta_rows)
> sprintf("True = 0.8, Stan_Est = %.3f, Optim_Est = %.3f", beta_par$mean, optim_par[1])
[1] "True = 0.8, Stan_Est = 1.369, Optim_Est = 1.111"

あらら、optimの結果よりも外れてしまいました。Shape効果の方はどうでしょうか。

> ests %>% 
+     data.frame() %>% 
+     select(mean) %>% 
+     mutate(Par = rownames(ests),
+            Type = "Stan") %>% 
+     filter(Par %in% shape_rows) %>% 
+     bind_rows(data.frame(
+         mean = c(0.4, 4, optim_par[-1]),
+         Par = rep(c("shape_k_01", "shape_s_01"), 2),
+         Type = c(rep("True", 2), rep("Optim", 2))
+     )) %>% 
+     spread(Type, mean) %>% 
+     select(Par, True, Stan, Optim)
         Par True      Stan     Optim
1 shape_k_01  0.4 0.6308075 0.4965623
2 shape_s_01  4.0 1.8736434 1.9041772

こちらもoptimの方がまだ設定値に近いようですね。ただoptimstanもどちらもうまく推定できていない様子なので比較に意味はないかもしれませんが。

終わりに

というわけで、今回はあまりパッとしない結果となってしまいました。元論文はちゃんと読んでいないのですが、彼らは階層ベイズで解いているようなので、もう少しパラメータに対して制約をかけてあげないとダメなのかもしれません。

本当なら良い推定値が得られるようになるまで頑張りたいのですが、ひとまず思いつくことは試した上でこの結果なので、いったん諦めて公開することとしました。何かアイディアが浮かべば再挑戦したいと思います。


  1. 要するにAd-Stock効果のパラメータの同時推定がうまく行かなかったんです。。。