統計コンサルの議事メモ

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

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

先日、Marketing Mix Modelingにおける広告の累積効果(Ad-Stock効果)の推定について以下のような記事を書いた。

ushi-goroshi.hatenablog.com

その後も推定方法について調べていたところ、以下のような記事を発見した:

www.mm-lab.jp

要するに一期前の目的変数そのもの(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を用いる方法では各変数の減衰率λの値も推定することができ、その精度も十分高い。

まとめ

上記のシミュレーションの結果をまとめると:

  1. 目的変数の一期前の数値を説明変数として用いる(パルダ・モデル)ことで、Ad-Stock効果をモデルに取り入れることができる
  2. 簡単なシミュレーションではパルダ・モデルの有効性が示され、モデルを非線形にすることなく精度の良い回帰係数を推定することができた
  3. ただしoptimによる推定の方が回帰係数の推定精度が高く、またパルダ・モデルでは各変数の減衰率λを同時に得ることができない

となった。
パルダ・モデルではモデルを非線形にする必要がないためlmにより簡単に推定でき、またその推定精度も高い。さらに過去の累積効果をまとめて説明することができるため、Marketing Mix Modelの構築において最適なLagを気にする必要がなくなることは大きなメリットである。そのため一般的なMMMの構築においてはパルダ・モデルの方が望ましいのではないか。
一方で、パルダ・モデルは各変数のλを得ることができないため詳細な考察が必要な場合には向いておらず、また分布ラグの形もコイック型に限られてしまう。例えば広告間の効率を比較するためにλを示す必要があったり、使用期限が定められたクーポンのように幾何級数的な減衰が適当でない(期限ギリギリになると効果が上昇する)ことが想定されるような場合にはoptimによる推定の方が適当であると考えられる。

したがって分析者が必要とする数値および対象領域における仮定を踏まえて推定方法を選択する必要があるだろう。