Ad-Stock効果を推定しつつ回帰を回したい②
先日、Marketing Mix Modelingにおける広告の累積効果（Ad-Stock効果）の推定について以下のような記事を書いた。
その後も推定方法について調べていたところ、以下のような記事を発見した：
要するに一期前の目的変数そのもの（y_{t-1}）を説明変数として含めることで、過去の累積効果、すなわちAd-Stock効果を説明することができるらしい。これを発表者の名前からパルダモデルと呼ぶとのこと。
数式を用いた説明は以下で発見できた：
http://www-cc.gakushuin.ac.jp/~20130021/ecmr/chapter6.pdf
この記事において、コイック型の分布ラグを仮定したモデルはARX(1)と等価であることが解説されている。
それでは、このパルダ・モデルと先日のoptimによる推定ではどちらがより精度よく推定できるのだろうか？
シミュレーションの実行：パルダ・モデル
パルダ・モデルの推定精度を検証するためにシミュレーションを実行した。
まずは先日の記事同様にシミュレーション用のサンプルデータを発生させるための関数を作成する。真のモデルはstats::filterによって生成された説明変数にbetaを乗じたもので、このbetaを元データから推定したい。
simulate_y <- function(pars) { n <- pars[1] # num of observation mu <- pars[2] # intercept beta_01 <- pars[3] # regression coefficient of X1 to be esitmated beta_02 <- pars[4] # regression coefficient of X2 to be esitmated lambda_01 <- pars[5] # decay rate of Ad-Stock effect of X1 lambda_02 <- pars[6] # decay rate of Ad-Stock effect of X2 var_e <- pars[7] # residual variance X_01_raw <- rnorm(n, 100, 2) * ifelse(runif(n) > 0.7, 0, 1) X_01_fil <- stats::filter(X_01_raw, lambda_01, "recursive") X_02_raw <- rnorm(n, 100, 2) * ifelse(runif(n) > 0.2, 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, "Y_lag" = dplyr::lag(y, 1)) return(na.omit(dat)) }
seedを固定し、500回の反復シミュレーションを実行する。lmの中にY_lagが含められていることに注意。betaの真の値は0.5および0.3と設定した。また（今回は推定の対象ではないが）減衰率はそれぞれ0.9および0.2である。
## simulation set.seed(123) init_par <- array(c(100, 5, 0.5, 0.3, 0.9, 0.2, 5)) k <- 500 results_palda <- matrix(0, k, 4) for (i in 1:k) { dat <- simulate_y(init_par) res <- lm(Y ~ X_01 + X_02 + Y_lag, data = dat) results_palda[i, ] <- coef(res) }
> summary(results_palda) V1 V2 V3 V4 Min. :-15.4731 Min. :0.4322 Min. :0.2117 Min. :0.8300 1st Qu.: -1.8896 1st Qu.:0.4851 1st Qu.:0.2900 1st Qu.:0.8741 Median : 0.8618 Median :0.5000 Median :0.3049 Median :0.8820 Mean : 1.0207 Mean :0.5002 Mean :0.3037 Mean :0.8820 3rd Qu.: 4.0467 3rd Qu.:0.5146 3rd Qu.:0.3175 3rd Qu.:0.8907 Max. : 16.4760 Max. :0.5597 Max. :0.3624 Max. :0.9254
V2およびV3に注目すると、それぞれ0.43〜0.56、0.21〜0.36の範囲で分布しており、Medianおよび平均は両方とも一致している。どうやらパルダ・モデルは過去の累積による影響をうまく取り除き、広告の当期に対する影響を精度よく推定できるようだ。
なおV4は一期前の目的変数に対するyの回帰係数であり、先に示したpdfによるとこの値は広告の減衰率λに相当するはずである。しかし変数が2つあるためか一致していない。
シミュレーションの実行：optimによる推定
では先日と同様にoptimを使いながら減衰率を同時に推定する方法はどうか。まずは先日と同じく、optimの目的関数を作成する。
return_AIC <- function(param) { a <- param[1] b <- param[2] dat$X_01_fil <- stats::filter(dat$X_01, a, "recursive") dat$X_02_fil <- stats::filter(dat$X_02, b, "recursive") AIC(lm(Y ~ X_01_fil + X_02_fil, dat)) }
続いてシミュレーションを実行する。反復回数はパルダ・モデルと同じく500回で、betaの値も同じである。実行にやや時間がかかったので進捗を見られるように100回ごとにRoundを表示した。
set.seed(456) init_par <- array(c(100, 5, 0.5, 0.3, 0.9, 0.2, 5)) decay_par <- array(c(0.5, 0.5)) k <- 500 results_optim <- matrix(0, k, 5) for (i in 1:k) { dat <- simulate_y(init_par) res <- optim(par = optim(par = decay_par, fn = return_AIC)$par, fn = return_AIC)$par dat$X_01_fil <- stats::filter(dat$X_01, res[1], "recursive") dat$X_02_fil <- stats::filter(dat$X_02, res[2], "recursive") lm_out <- lm(Y ~ X_01_fil + X_02_fil, data = dat) results_optim[i, ] <- c(res, coef(lm_out)) if (i %% 50 == 0) cat("Round: ", i, "\n") }
> summary(results_optim) V1 V2 V3 V4 V5 Min. :0.8866 Min. :0.06162 Min. : 1.983 Min. :0.4608 Min. :0.2649 1st Qu.:0.8916 1st Qu.:0.17331 1st Qu.: 6.280 1st Qu.:0.4838 1st Qu.:0.2941 Median :0.8941 Median :0.19419 Median :35.109 Median :0.4923 Median :0.3001 Mean :0.8949 Mean :0.19289 Mean :27.462 Mean :0.4908 Mean :0.3002 3rd Qu.:0.8996 3rd Qu.:0.21262 3rd Qu.:38.388 3rd Qu.:0.4989 3rd Qu.:0.3059 Max. :0.9010 Max. :0.29186 Max. :46.506 Max. :0.5200 Max. :0.3286
V4およびV5を見ると、それぞれの範囲は0.46〜0.52および0.26〜0.33と、パルダ・モデルと比較して精度よく推定できている。加えてoptimを用いる方法では各変数の減衰率λの値も推定することができ、その精度も十分高い。
まとめ
上記のシミュレーションの結果をまとめると：
- 目的変数の一期前の数値を説明変数として用いる（パルダ・モデル）ことで、Ad-Stock効果をモデルに取り入れることができる
- 簡単なシミュレーションではパルダ・モデルの有効性が示され、モデルを非線形にすることなく精度の良い回帰係数を推定することができた
- ただしoptimによる推定の方が回帰係数の推定精度が高く、またパルダ・モデルでは各変数の減衰率λを同時に得ることができない
となった。
パルダ・モデルではモデルを非線形にする必要がないためlmにより簡単に推定でき、またその推定精度も高い。さらに過去の累積効果をまとめて説明することができるため、Marketing Mix Modelの構築において最適なLagを気にする必要がなくなることは大きなメリットである。そのため一般的なMMMの構築においてはパルダ・モデルの方が望ましいのではないか。
一方で、パルダ・モデルは各変数のλを得ることができないため詳細な考察が必要な場合には向いておらず、また分布ラグの形もコイック型に限られてしまう。例えば広告間の効率を比較するためにλを示す必要があったり、使用期限が定められたクーポンのように幾何級数的な減衰が適当でない（期限ギリギリになると効果が上昇する）ことが想定されるような場合にはoptimによる推定の方が適当であると考えられる。
したがって分析者が必要とする数値および対象領域における仮定を踏まえて推定方法を選択する必要があるだろう。
