読者です 読者をやめる 読者になる 読者になる

統計コンサルの議事メモ

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

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

最近ずっとMarketing Mix Modeling(MMM)をやっている。その中で広告効果(いわゆるROI)を推定しているのだけれど、広告の効果というものは出稿した時点だけでなく将来に渡って影響を及ぼすため、過去の広告の累積による影響(いわゆる残存効果・Ad-Stock効果)を取り上げる必要がある。ところが、このAd-Stock効果をモデルにより推定しようとすると少し難しい。

Ad-Stock効果として以下のような分布ラグ型:

{ \displaystyle
y_t = \mu + \beta_{t} X_{t} + \beta_{t-1} X_{t-1} + \cdots + \beta_{t-m} X_{t-m} + \epsilon_t
}

を仮定するのが最も単純なのだけれど、考慮するラグの次数mによってはモデルに投入する変数の数が多くなり、また多重共線性も強くなりがち。時系列データを用いているためデータポイントが少なく、このような方法は避けたい。

そこで「m期前の広告効果はλm乗に従って減衰する」という仮定をおき(いわゆるKoyck型)、推定すべきパラメータの数を減らしている。この場合変数の数は減らせるが、非線形な効果であるλをどのように推定するかが鍵となる。

考えられる方法として、

  1. データ加工の段階で例えば相関係数などから推定する
  2. 過去の研究や業界での知見を導入する
  3. 非線形な効果を推定する方法を考える

といったものが思いついた。しかし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"を指定すれば上限・下限を指定できたと思うので、もう少し検証してみよう。