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法」を指定することで矩形制約を与えることができる