Ad-Stock効果を推定しつつ回帰を回したい
最近ずっとMarketing Mix Modeling（MMM）をやっている。その中で広告効果（いわゆるROI）を推定しているのだけれど、広告の効果というものは出稿した時点だけでなく将来に渡って影響を及ぼすため、過去の広告の累積による影響（いわゆる残存効果・Ad-Stock効果）を取り上げる必要がある。ところが、このAd-Stock効果をモデルにより推定しようとすると少し難しい。
Ad-Stock効果として以下のような分布ラグ型：
を仮定するのが最も単純なのだけれど、考慮するラグの次数mによってはモデルに投入する変数の数が多くなり、また多重共線性も強くなりがち。時系列データを用いているためデータポイントが少なく、このような方法は避けたい。
そこで「m期前の広告効果はλのm乗に従って減衰する」という仮定をおき（いわゆるKoyck型）、推定すべきパラメータの数を減らしている。この場合変数の数は減らせるが、非線形な効果であるλをどのように推定するかが鍵となる。
考えられる方法として、
といったものが思いついた。しかし1は季節性や他の広告の効果などを同時に考慮できないため✖️、2は使えそうな数値がなかった。そこで今回は回帰の中でλを同時に推定できないか、以下のように調査してみた。
set.seed(123) # λに従って広告の効果が減衰するデータを生成する関数 simulate_y_01 <- function(pars) { X1 <- rnorm(100, 100, 2) * rep(c(0, rep(1, 7), 0, 1), 10) X2 <- filter(X1, pars[1], "recursive") Z1 <- rnorm(100, 5, 2) * rep(c(rep(0, 5), rep(1, 5)), 10) dat <- data.frame( "Y" = 50 + pars[2] * X2 + pars[3] * Z1 + rnorm(100, 0, pars[4]), "X1" = X1, "Z1" = Z1) return(dat) } # λを与えるとAd-Stock変数を作成して回帰を実施し、AICを返す関数 return_AIC_01 <- function(a) { dat$Tmp <- filter(dat$X1, a, "recursive") AIC(lm(Y ~ Tmp + Z1, dat)) }
Koyck型の変数の作成にはstats::filterが便利だが、使う際には{dplyr}を読み込んでいないか注意しよう。これでreturn_AIC_01を最小化した時に得られるλが、simulate_y_01で用いた数値を一致するか確認する。実際にやってみると以下の結果が得られる。
## λを小さな値とした場合 pars <- c(0.1, 0.5, 4, 5.0) n <- 500 res_01_01 <- matrix(0, n, 4) for (i in 1:n) { dat <- simulate_y_01(pars) a <- optimise(return_AIC_01, c(0, 1))$minimum X2 <- filter(dat$X1, a, "recursive") res_01_01[i, ] <- c(a, coef(lm(Y ~ X2 + Z1, data = dat))) } ## λを大きな値とした場合 pars <- c(0.9, 0.5, 4, 5.0) n <- 500 res_01_02 <- matrix(0, n, 4) for (i in 1:n) { dat <- simulate_y_01(pars) a <- optimise(return_AIC_01, c(0, 1))$minimum X2 <- filter(dat$X1, a, "recursive") res_01_02[i, ] <- c(a, coef(lm(Y ~ X2 + Z1, data = dat))) }
> summary(res_01_01) V1 V2 V3 V4 Min. :0.01545 Min. :43.81 Min. :0.4633 Min. :3.486 1st Qu.:0.08780 1st Qu.:48.73 1st Qu.:0.4918 1st Qu.:3.877 Median :0.10450 Median :49.90 Median :0.4998 Median :3.999 Mean :0.10324 Mean :49.88 Mean :0.4996 Mean :3.997 3rd Qu.:0.11992 3rd Qu.:51.02 3rd Qu.:0.5081 3rd Qu.:4.109 Max. :0.17782 Max. :55.08 Max. :0.5408 Max. :4.474 > summary(res_01_02) V1 V2 V3 V4 Min. :0.8972 Min. :40.15 Min. :0.4768 Min. :3.414 1st Qu.:0.8992 1st Qu.:47.93 1st Qu.:0.4946 1st Qu.:3.869 Median :0.9000 Median :50.04 Median :0.4995 Median :3.994 Mean :0.9000 Mean :49.94 Mean :0.5002 Mean :3.999 3rd Qu.:0.9007 3rd Qu.:52.10 3rd Qu.:0.5060 3rd Qu.:4.133 Max. :0.9031 Max. :58.68 Max. :0.5236 Max. :4.526
λが小さい時にはある程度バラついているものの、この誤差なら十分に実用的に思える。λが大きい時は全く問題がなさそうだ。
少し問題を難しくするため、2つの変数を共にAd-Stockに変更すると：
simulate_y_02 <- function(pars) { X1 <- rnorm(100, 100, 2) * rep(c(0, rep(1, 7), 0, 1), 10) X2 <- filter(X1, pars[1], "recursive") Z1 <- rnorm(100, 5, 2) * rep(c(rep(0, 5), rep(1, 5)), 10) Z2 <- filter(Z1, pars[2], "recursive") dat <- data.frame( "Y" = 50 + pars[3] * X2 + pars[4] * Z2 + rnorm(100, 0, pars[5]), "X1" = X1, "Z1" = Z1) return(dat) } return_AIC_02 <- function(param) { a <- param[1] b <- param[2] dat$TmpX <- filter(dat$X1, a, "recursive") dat$TmpZ <- filter(dat$Z1, b, "recursive") AIC(lm(Y ~ TmpX + TmpZ, dat)) } ## λを小さな値とした場合 pars <- c(0.1, 0.1, 0.5, 4, 5.0) n <- 100 res_02_01 <- matrix(0, n, 5) for (i in 1:n) { dat <- simulate_y(pars) res <- optim(par = optim(par = init_par, fn = return_AIC)$par, fn = return_AIC)$par X2 <- filter(dat$X1, res[1], "recursive") Z2 <- filter(dat$Z1, res[2], "recursive") res_02_01[i, ] <- c(res, coef(lm(Y ~ X2 + Z2, data = dat))) } ## λを大きな値とした場合 pars <- c(0.8, 0.7, 0.5, 4, 5.0) n <- 100 res_02_02 <- matrix(0, n, 5) for (i in 1:n) { dat <- simulate_y(pars) res <- optim(par = optim(par = init_par, fn = return_AIC)$par, fn = return_AIC)$par X2 <- filter(dat$X1, res[1], "recursive") Z2 <- filter(dat$Z1, res[2], "recursive") res_02_02[i, ] <- c(res, coef(lm(Y ~ X2 + Z2, data = dat))) }
計算に少し時間がかかったので、反復は100回とした。結果は以下の通り。
> summary(res_02_01) V1 V2 V3 V4 V5 Min. :0.03261 Min. :-0.02245 Min. :44.54 Min. :0.4685 Min. :3.597 1st Qu.:0.08483 1st Qu.: 0.06107 1st Qu.:48.62 1st Qu.:0.4899 1st Qu.:3.902 Median :0.09598 Median : 0.09694 Median :50.23 Median :0.4997 Median :3.997 Mean :0.09810 Mean : 0.09459 Mean :50.06 Mean :0.4997 Mean :4.017 3rd Qu.:0.11280 3rd Qu.: 0.11921 3rd Qu.:51.30 3rd Qu.:0.5081 3rd Qu.:4.122 Max. :0.14426 Max. : 0.24930 Max. :55.33 Max. :0.5419 Max. :4.606 > summary(res_02_02) V1 V2 V3 V4 V5 Min. :0.7894 Min. :0.6499 Min. :41.16 Min. :0.4630 Min. :3.701 1st Qu.:0.7961 1st Qu.:0.6833 1st Qu.:48.19 1st Qu.:0.4905 1st Qu.:3.906 Median :0.7990 Median :0.7003 Median :50.51 Median :0.5016 Median :3.984 Mean :0.7997 Mean :0.7000 Mean :49.99 Mean :0.5006 Mean :3.994 3rd Qu.:0.8030 3rd Qu.:0.7140 3rd Qu.:52.09 3rd Qu.:0.5093 3rd Qu.:4.065 Max. :0.8137 Max. :0.7496 Max. :57.14 Max. :0.5257 Max. :4.287
先ほどと変わらず、λが小さいと推定の精度が落ちるようだ。また一部のデータセットではλが負となってしまった。
ついでなので、{rBayesianOptimization}パッケージによる最適化を通常のoptimと比較してみた。
library(rBayesianOptimization) # BayesianOptimization関数では最大化を行うので、返り値となるAICに負の符号を加えてある return_AIC_BO <- function(a, b) { dat$TmpX <- filter(dat$X1, a, "recursive") dat$TmpZ <- filter(dat$Z1, b, "recursive") out <- list(Score = -AIC(lm(Y ~ TmpX + TmpZ, dat)), Pred = -AIC(lm(Y ~ TmpX + TmpZ, dat))) return(out) } pars <- c(0.1, 0.1, 0.5, 4, 5.0) dat <- simulate_y(pars) res_BO <- BayesianOptimization(return_AIC_BO, bounds = list(a = c(0, 1), b = c(0, 1)), init_points = 20, n_iter = 1, acq = "ucb", kappa = 2.576, eps = 0, verbose = TRUE) init_par <- array(c(0.5, 0.5)) res_Opt <- optim(par = optim(par = init_par, fn = return_AIC)$par, fn = return_AIC)$par
> res_BO$Best_Par a b 0.8274328 0.4757574 > res_Opt [1] 0.7959288 0.7930622
BayesianOptimizationの結果はあまり良くなく、通常のoptimの方が精度よく推定できた。以前の記事（optimってあんまり信用できないなぁ、って話 - 統計コンサルの議事メモ）ではoptimによるパラメータの推定はうまくいかなかったのだけど、lmと組み合わせているのが良いのかもしれない。
この方法でならAd-Stock効果の推定できそうだ。optimでは"L-BFGS-B"を指定すれば上限・下限を指定できたと思うので、もう少し検証してみよう。
