統計コンサルの議事メモ

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

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

これまでにMarketing Mix Modelingにおける広告の累積効果の推定について、以下のような記事を書いてきた。

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

今回は遺伝的アルゴリズムを用いてパラメータの推定を行ったので、その結果を示しておく。

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

前回同様、必要なパッケージをインストールしてから読み込む。インストールするパッケージは遺伝的アルゴリズムの実装として非常にメジャーな{GA}。

install.packages("GA")
library(GA)

続いてシミュレーションデータ作成用の関数を定義する。これは先日の記事と全く同じである。

simulate_y <- function(pars) {
   ## Simulation parameters
   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
   
   ## Create true Ad-Stock variables
   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")
   
   ## Create residuals
   error <- rnorm(n, 0, sqrt(var_e))
   
   ## Create observations   
   y     <- mu + beta_01 * X_01_fil + beta_02 * X_02_fil + error
   
   ## Return dataset
   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)
}

この関数を用いてサンプルデータを発生させる。乱数のシードを固定しているので、前回と同じデータが得られる。推定対象となるパラメータも同様に前回と同じ。

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

さらに、GAで用いる評価関数(目的関数)を最小二乗法で定義する。

OLS <- function(dat, pars) {
   X_01_Fil_E <- stats::filter(dat$X_01, pars[3], "recursive")
   X_02_Fil_E <- stats::filter(dat$X_02, pars[5], "recursive")
   Y_hat <- pars[1] + pars[2] * X_01_Fil_E + pars[4] * X_02_Fil_E
   SSE   <- t(dat$Y - Y_hat) %*% (dat$Y - Y_hat)
   return(SSE)
}

結果

上記で準備が整ったので、以下のようにGAを回してみる。

res_GA <- ga(type = 'real-valued', 
             min = c(0, 0, 0, 0, 0), 
             max = c(20, 2, 1, 2, 1),
             popSize = 500, maxiter = 500, 
             names = c('Intercept', 'b_01', 'l_01', 'b_02', 'l_02'),
             keepBest = T, 
             fitness = function(b) -OLS(dat_Ana, b))

{stan}で実行した時と同様に、{GA}を用いた場合では各パラメータに対してminとmaxを指定することでパラメータスペースに矩形制約をかけることができる。

結果を見てみよう。

> summary(res_GA)$solution
     Intercept      b_01      l_01     b_02      l_02
[1,]  4.729741 0.4938017 0.6421675 0.718373 0.5559099
> init_par[c(-1, -3)] # サンプルサイズと誤差分散を除外している
[1] 5.0 0.5 0.6 0.8 0.5

推定精度はかなり良さそうだ。
なおこれらの数値は{stan}で推定したものとほとんど同様であり、前回の記事と同じくb_02の数値の誤差が大きいようである。しかしサンプルデータを発生させた時の誤差によるかもしれないので、lmでフィットした時の結果も確認しておこう。

> lm(Y ~ X_01_Fil + X_02_Fil, dat_Ana)

Call:
lm(formula = Y ~ X_01_Fil + X_02_Fil, data = dat_Ana)

Coefficients:
(Intercept)     X_01_Fil     X_02_Fil  
     5.0846       0.5170       0.7277  

やはり、なかなかの精度で推定できていたようだ。誤差分散も見てみる。

> -summary(res_GA)$fitness / (nrow(dat_Ana) - 6)
[1] 2.064139
> init_par[3]
[1] 2

誤差分散の推定値は{GA}の返り値(SSE)をサンプルサイズ - 6で除して求めた。ここで6はパラメータ数である。誤差分散も極めて精度よく推定できているようだ。

まとめ

これまでの結果と同様に、{GA}を用いたAd-Stock効果の推定はなかなか実用的であると言えそうだ。特にMarketing Mix Modelingによって広告効果の検証を行う場合には回帰係数を0以上に制約をかけたいことが多いため、パラメータの範囲を指定できる点はパルダ・モデルやoptimと比較すると大きなメリットとなるだろう*1。{GA}のインストールも簡単で、{stan}と異なり躓くポイントも少ない。

一方、{stan}同様に分析に要する時間がパルダ・モデルやoptimと比較して長いため、分析の初期段階での使用は向いていないようだ。ある程度モデルの形が見えてきた段階でパラメータを精度よく推定したい場合に用いるのが良いだろうか。

これまで様々な手法を用いてAd-Stock効果の推定を試してきたが、いずれの方法も単純なモデルでは十分に実用的と言える精度でパラメータを推定することができていた。そのため今後はより複雑・より現実を反映したデータを用いた時に、各手法の推定精度がどの程度悪化するかを検証していく必要がある。

*1:念のために言っておくと、optimでも「L-BFGS-B法」を指定することで矩形制約を与えることができる