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法」を指定することで矩形制約を与えることができる
Bayesian t-testを書いてみたい
IBM SPSS Statisticsにベイズ推論が実装されるらしい。
これ以外にも最近は身の回りでやけにBayesian t-testに関する話題を耳にすることが多いので、t-検定をベイズ化することのどこにニーズがあるのかを考えていたのだけれど、どうやらA/B testにおいて結果を途中で確認（Peeking）したり、サンプルサイズを変更したいという状況が結構多いそうだ。
古典的な、そして計画立てられたt-検定では分析の事前に必要なサンプルサイズを求めておく。通常は実験の開始以降にデータを追加したりなどしないのだが、ビジネスではA/Bテストの結果の如何によってデータを追加したり、あるいは途中経過を見て早めに施策の実施判断を行いたいと状況がよくあるため、そういったニーズにマッチしたのだろう（実際のところ伝統的なt-検定を使う場合であってもデータの追加や結果の覗き見は行われていると思うが）。
さて、Rにおけるベイジアンモデリングのためのパッケージとしては{MCMCpack}があり、より複雑なモデルを構築したいのであればStanを呼び出す方法があげられる。しかしt-検定ぐらい自分で書いてみようと思ったので、以下にスクリプトを示す。
データの準備
# パッケージの読み込み library(tidyverse) library(ggplot2)
サンプルデータとして以下のような正規分布に従う２標本を用意する。数値に意味はないので関心があれば標本サイズなどをそれぞれ変更してみると良いが、目的は標本の平均の差の分布について確認することである。
N_01 <- 100 mu_01 <- 100 sig_01 <- 10 N_02 <- 100 mu_02 <- 103 sig_02 <- 10 set.seed(123) Pop_01 <- rnorm(N_01, mu_01, sig_01) Pop_02 <- rnorm(N_02, mu_02, sig_02)
以下のようなサンプルが得られるだろう：
> Pop_01 %>% as_data_frame() %>% mutate("Group" = "01") %>% bind_rows(as_data_frame(Pop_02) %>% mutate("Group" = "02")) %>% ggplot(aes(value, fill = Group)) + geom_histogram(alpha = 0.5, position = "identity") + theme_bw()
このデータに対して古典的なt-検定を実施すると、平均値の差が0であるという帰無仮説はp = 0.44で棄却されない。ちなみにt.testはデフォルトでは等分散性を仮定しないWelchのt-検定となる。
> t.test(Pop_01, Pop_02) Welch Two Sample t-test data: Pop_01 and Pop_02 t = -0.7674, df = 197.35, p-value = 0.4438 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.642862 1.601916 sample estimates: mean of x mean of y 100.9041 101.9245
関数の定義
このデータに対して何らかのパラメータ（ここでは平均と標準偏差）を与えた時の尤度を求める関数を以下のように定義する：
likelihood <- function(pars) { mu_01 <- pars[1] sig_01 <- pars[2] mu_02 <- pars[3] sig_02 <- pars[4] sum( dnorm(Pop_01, mu_01, sig_01, log = TRUE), dnorm(Pop_02, mu_02, sig_02, log = TRUE) ) }
データの尤度は0 ~ 1の数字の乗算となるため標本サイズが大きい場合には極めて小さな数値となり、桁あふれが発生してしまう。そのためdnormにおいてlog = TRUEとし、対数尤度の総和を求めることとする。
この関数を用いて最尤法で平均と分散を求めてみる。optimを使用するが、デフォルトは最小化なのでfnscaleを負の数として最大化を行う：
> pars <- c(100, 10, 100, 20) > optim(par = optim(par = pars, fn = likelihood, control=list(fnscale = -1))$par, + fn = likelihood, control=list(fnscale = -1)) $par [1] 100.904769 9.082697 101.923277 9.620816 $value [1] -730.8205 $counts function gradient 101 NA $convergence [1] 0 $message NULL
続いて事前分布を定義する。事前分布の設定は色々と工夫すべきポイントではあるが、簡単のため-10000 ~ 10000の一様分布とした。尤度はデータと同じく対数としてある。
prior <- function(pars) { mu_01 <- pars[1] sig_01 <- pars[2] mu_02 <- pars[3] sig_02 <- pars[4] sum( dunif(mu_01, -10000, 10000, log = TRUE), dunif(mu_02, -10000, 10000, log = TRUE), dunif(sig_01, -10000, 10000, log = TRUE), dunif(sig_02, -10000, 10000, log = TRUE) ) }
事後分布は事前分布と尤度の積であるが、対数を取ってあるためここでは和となる。
posterior <- function(pars) { likelihood(pars) + prior(pars) }
MCMCの実行
以上で準備ができたので以下のようにMCMCを実行する。なおサンプラーとしてGibbs Samplerを実装したつもりだが、Gibbs Samplingにおいてパラメータの採択確率があったかは記憶が不確かなので、この辺りは後日確認しておくとしてひとまず進めることとする。サンプルの数は30,000としたが、実際には割と早く収束しているようなのでこんなには必要ない。
N_iter <- 30000 pars <- c(50, 30, 20, 50) results <- matrix(0, nrow = N_iter, ncol = 8) results[1, ] <- c(pars, 0, 0, 0, 0) for (i in 2:N_iter) { for (j in 1:4) { candidate <- pars candidate[j] <- candidate[j] + rnorm(1, sd = 0.5) ratio <- exp(posterior(candidate) - posterior(pars)) t <- runif(1) if (t < ratio) { pars <- candidate } } results[i, ] <- c(pars, ratio, t, posterior(pars), posterior(candidate)) }
以下のような結果となる。
> par(mfrow = c(2, 2)) > plot(results[, 1]) > plot(results[, 2]) > plot(results[, 3]) > plot(results[, 4])
要約統計量も見ておこう。burn-inとして最初の15,000サンプルを切り捨て、自己回帰性を小さくするために10刻みでサンプルを取得する。
> Res_MCMC <- results[seq(15000, N_iter, 10), ] > summary(Res_MCMC[, c(1:4)]) V1 V2 V3 V4 Min. : 97.79 Min. : 7.198 Min. : 98.41 Min. : 7.997 1st Qu.:100.24 1st Qu.: 8.732 1st Qu.:101.34 1st Qu.: 9.266 Median :100.86 Median : 9.191 Median :101.98 Median : 9.753 Mean :100.86 Mean : 9.233 Mean :101.98 Mean : 9.774 3rd Qu.:101.45 3rd Qu.: 9.675 3rd Qu.:102.63 3rd Qu.:10.218 Max. :104.16 Max. :11.959 Max. :105.20 Max. :12.632
単純な推定なので精度よく推定できているようである。
平均値の差の分布
では最後に、この分析の目的であった２標本の平均の差の分布について確認してみよう。
par(mfrow = c(1, 1)) hist(Res_MCMC[, 1] - Res_MCMC[, 3])
これを見ると、平均値の差の分布が比較的中心に近い位置で0を跨いでおり、伝統的なt-検定で言う所の有意であるとは言えなさそうだ。実際に平均値の差の経験分布を確認すると、95%確信区間は-3.76 ~ 1.62である。
> quantile(Res_MCMC[, 1] - Res_MCMC[, 3], seq(0, 1, 0.025)) 0% 2.5% 5% 7.5% 10% 12.5% 15% 17.5% -5.62647895 -3.76661174 -3.35462735 -3.11314731 -2.89345678 -2.70960554 -2.56386382 -2.43985754 20% 22.5% 25% 27.5% 30% 32.5% 35% 37.5% -2.30487729 -2.18744684 -2.07270632 -1.98450427 -1.89198288 -1.77273198 -1.67356108 -1.57205147 40% 42.5% 45% 47.5% 50% 52.5% 55% 57.5% -1.47598439 -1.40161516 -1.30787797 -1.21558957 -1.12424331 -1.03598088 -0.95253546 -0.85994024 60% 62.5% 65% 67.5% 70% 72.5% 75% 77.5% -0.78004294 -0.68704363 -0.59508922 -0.49850858 -0.42349299 -0.30716666 -0.19718374 -0.10822330 80% 82.5% 85% 87.5% 90% 92.5% 95% 97.5% 0.01189818 0.14609443 0.30742495 0.45296435 0.69891710 0.92152549 1.18895282 1.62042273 100% 3.58972078
改めてt-検定の結果を見てみると、95%信頼区間として同様の結果が得られていることがわかるだろう。
> t.test(Pop_01, Pop_02) Welch Two Sample t-test data: Pop_01 and Pop_02 t = -0.7674, df = 197.35, p-value = 0.4438 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -3.642862 1.601916 sample estimates: mean of x mean of y 100.9041 101.9245
おまけ
以上でt-検定のベイズ化は終わりだが、ベイズ化した分析のメリットとして以下のような疑問に答えることが可能となる：１つ目の標本の平均が２つ目の標本の平均よりも大きい確率はどの程度だろうか？
> mean(Res_MCMC[, 3] - Res_MCMC[, 1] < 0) [1] 0.2031979
上記のような質問は、古典的な仮説検定（頻度論的アプローチ）では答えることができない。仮説検定は常に「帰無仮説が棄却されるか否か」（例えば平均値と平均値に差が0であるか否か）であり、「何%の確率で棄却される」とは表現できないのだ。しかしこれをベイズ化することでより直感的な表現が可能となる。これはビジネスへの応用において確かにメリットになると感じた。
Random Forestの結果をSQLに落としたい
Random Forestは数ある機械学習アルゴリズムの中でも、高い精度を得やすく分散処理が容易であることから頻繁に用いられるものの一つであると思う。
Rでは{randomForest}パッケージが有名だが、最近では{ranger}や{Rborist}といった新しいパッケージが出ていたり、{edarf}のようなRandom Forestによるデータの解析や解釈を主眼にしたものが出ていたりして動きとして面白い。
そんなRandom Forestであるが、せっかくRで良いモデルを作ることができたとしても実際には使えないことがある。これは、例えばSQLに落とすことが容易ではないためシステムに組み込めないからだ。
ビジネスにおいてデータをRで分析する場面は数多くあるが、Rの結果を（文字通り）そのままシステムに埋め込むことはほとんどないと言って良いと思う。Rを使って得られたモデルは、SQLなどの他の環境に合わせて再構築した上で利用することが多いだろう。そこで今回はRの{randomForest}パッケージで作成したモデルをSQLに落とすための方法を紹介したい。
なおこの記事は基本的に以下の記事を参考に作成したが、ひとまずRの{randomForest}で作成したモデルからSQLを得るための流れを示すことが目的であったため大幅に簡略化してある。またSQLのインデントも私の好みに合わせたので、元記事も参照することをお薦めする。
https://gist.github.com/shanebutler/96f0e78a02c84cdcf558
処理の流れ
Random Forestというアルゴリズムは、大雑把に言えば簡易な決定木を大量に作成し、その結果の多数決によって最終的な判断を行うというものである。よってSQLに落とす場合も手順としては以下のようになる：
ただし実際には1と2は同時に行われるだろう。
1. 決定木の結果をSQLとして書き下す
それでは実装に移ろう。まずはSQLに書き下すところから。サンプルとしてirisを用い、簡単のためrandomForestの木の数は2としたモデルを作成する。
library(randomForest) set.seed(123) res <- randomForest(Species ~ ., data = iris, ntree = 2)
この処理で得られた結果はresオブジェクトに保存されるが、特にSQLを書くために必要となる決定木の分岐はgetTreeで取り出すことができる。取り出したい木が何番目であるかを指定し、ラベルの表示をTRUEにすることで：
> getTree(res, k = 1, labelVar = TRUE) left daughter right daughter split var split point status prediction 1 2 3 Petal.Width 0.80 1 <NA> 2 0 0 <NA> 0.00 -1 setosa 3 4 5 Petal.Width 1.65 1 <NA> 4 6 7 Petal.Width 1.35 1 <NA> 5 8 9 Sepal.Length 6.05 1 <NA> 6 0 0 <NA> 0.00 -1 versicolor 7 10 11 Petal.Length 4.95 1 <NA> 8 12 13 Petal.Length 4.85 1 <NA> 9 0 0 <NA> 0.00 -1 virginica 10 0 0 <NA> 0.00 -1 versicolor 11 0 0 <NA> 0.00 -1 virginica 12 14 15 Sepal.Length 5.40 1 <NA> 13 0 0 <NA> 0.00 -1 virginica 14 0 0 <NA> 0.00 -1 virginica 15 0 0 <NA> 0.00 -1 versicolor
このようなテーブルが得られる。決定木は親となるノードから次々に分岐していく形で木が形成されるが、このテーブルはそのためのルールを示したものとなっており、「left daughter」が該当する場合に次に参照する行を、「right daughter」が該当しなかった場合に参照する行を示す。
例えばこの場合、1行目を読むことで「Petal.Widthが0.8未満であるか」を基準に分岐が発生することがわかり、また該当する場合（left daughter）には2行目を参照することでsetosaになることを示している。
このようなテーブルを各決定木について得ることができるため、SQLに落とすための流れとしては：
- テーブルをSQLにするための処理を書く
- 決定木の数だけループを回す
ことが考えられるだろう。そこでこのテーブルをSQLに落とすための関数を以下のように定義する：
Make_SQL <- function(Rule_Table, cnt, ind = 1) { Rule <- Rule_Table[cnt, ] Indent <- paste(rep(" ", ind), collapse = "") var <- as.character(Rule[, "split var"]) val <- Rule[, "split point"] if(Rule[, "status"] != -1) { cat(paste0("\n", Indent, "CASE\n", Indent, " WHEN ", var, " <= ", val, " THEN")) Make_SQL(Rule_Table, Rule[, "left daughter"], ind = (ind + 2)) cat(paste0("\n", Indent, " ELSE")) Make_SQL(Rule_Table, Rule[, "right daughter"], ind = (ind + 2)) cat(paste0("\n", Indent, "END")) } else { cat(paste0(" '", Rule[, "prediction"], "'")) } }
ここでRule_Tableは先ほどのテーブルを示している。また関数の序盤で定義しているIndentはSQLを書くときのインデントを私の好みに合わせるために修正してある。
全体的に、SQLはCASE文で形成されることがわかる。これは決定木がIf THEN Elseで表現可能なルールによって作成されることからも自明だろう。なおRule_Tableのstatusには、当該ノードが末端（ターミナルノード）であるかが示されているため、ここを参照することで予測値が決まるか（それ以上の分岐がないか）がわかる。もしターミナルノードでないなら分岐した先で同様の処理を行うため、再帰的な呼び出しが必要であることがわかる。分岐した場合、該当するものについてはleft daughterの行を参照し、そうでないものについてはright daughterを参照するようになっている。
これを実行すると、例えば以下のようなSQL文が得られる：
> Make_SQL(Rule_Table, 1, 1) CASE WHEN Petal.Width <= 0.75 THEN 'setosa' ELSE CASE WHEN Petal.Width <= 1.75 THEN CASE WHEN Sepal.Width <= 2.25 THEN CASE WHEN Petal.Length <= 4.75 THEN 'versicolor' ELSE 'virginica' END ELSE CASE WHEN Petal.Width <= 1.35 THEN 'versicolor' ELSE CASE WHEN Sepal.Width <= 2.65 THEN CASE WHEN Petal.Width <= 1.45 THEN 'virginica' ELSE 'versicolor' END ELSE 'versicolor' END END END ELSE CASE WHEN Sepal.Length <= 5.95 THEN CASE WHEN Sepal.Width <= 3 THEN 'virginica' ELSE 'versicolor' END ELSE 'virginica' END END END
これをRandom Forestのモデル作成に用いた決定木の数だけ作成すれば良いということになる。なお決定木の数を大きなものとした場合にSQL文が大きくなってしまうため、実行の際には注意が必要である。
このようにして得られたSQLはこのままだと表示して終わりであるため、何らかのテーブルに格納しておく必要がある。そこでINSERT文を前に置いておくことにする。すなわち：
cat(paste0("INSERT INTO tbl_rf \nSELECT\n id, ")) Make_SQL(Rule_Table, 1, 1)
このような書き方になる。ここでidは予測する場合の粒度を示すことになる。irisの場合は1行が1データ（個体）になるため行ごとに連番を振れば良いが、ここでは省略する。
これで一つの決定木についてSQLによる結果をテーブルに格納することができるようになったため、上記のスクリプトをループで回すと以下のようになるだろう：
for (i in 1:(res$ntree)) { Rule_Table <- getTree(res, k = i, labelVar = TRUE) cat(paste0("INSERT INTO tbl_rf \nSELECT\n id, ")) Make_SQL(Rule_Table, 1, 1) cat(paste0(" as tree", i, "\nFROM\n", " input_data", ";\n\n")) }
res$ntreeにはモデルの作成に用いた決定木の数が格納されているため、その分だけループを回すことになる。またFROM句では決定木を当てはめるためのデータ（この例ではiris）を指定することになる。
1つのテーブルに各SQLの結果を集約する
上記の処理によりtbl_rfには各決定木の結果が格納されているだろう。これを用いて最終的な判断を行うには、各決定木の結果を集約する必要がある。以下のようになるだろう：
> cat(paste0("INSERT INTO rf_predictions\n", + "SELECT\n a.id,\n a.pred\n", + "FROM\n (\n SELECT id as id, pred, COUNT(*) as cnt \n", + " FROM tbl_rf\n GROUP BY id, pred\n ) a\n", + " INNER JOIN\n (\n ", + "SELECT id, MAX(cnt) as cnt\n", + " FROM\n", + " (\n", + " SELECT id as id, pred, COUNT(*) as cnt\n", + " FROM tbl_rf\n", + " GROUP BY id, pred\n", + " )\n", + " GROUP BY id\n ) b\n", + " ON a.id = b.id AND a.cnt = b.cnt;\n\n")) INSERT INTO rf_predictions SELECT a.id, a.pred FROM ( SELECT id as id, pred, COUNT(*) as cnt FROM tbl_rf GROUP BY id, pred ) a INNER JOIN ( SELECT id, MAX(cnt) as cnt FROM ( SELECT id as id, pred, COUNT(*) as cnt FROM tbl_rf GROUP BY id, pred ) GROUP BY id ) b ON a.id = b.id AND a.cnt = b.cnt;
irisのSpeciesを用いる場合はclassificationとなるため、各決定木の判断を集約し、最も多かったクラスを最終的な予測とする。上記はそのためのクエリとなっていることがわかる。
最後に、これまでのスクリプトを一つにまとめたものを以下に示す：
library(randomForest) set.seed(123) res <- randomForest(Species ~ ., data = iris, ntree = 2) Make_SQL <- function(Rule_Table, cnt, ind = 1) { Rule <- Rule_Table[cnt, ] Indent <- paste(rep(" ", ind), collapse = "") var <- as.character(Rule[, "split var"]) val <- Rule[, "split point"] if(Rule[, "status"] != -1) { cat(paste0("\n", Indent, "CASE\n", Indent, " WHEN ", var, " <= ", val, " THEN")) Make_SQL(Rule_Table, Rule[, "left daughter"], ind = (ind + 2)) cat(paste0("\n", Indent, " ELSE")) Make_SQL(Rule_Table, Rule[, "right daughter"], ind = (ind + 2)) cat(paste0("\n", Indent, "END")) } else { cat(paste0(" '", Rule[, "prediction"], "'")) } } for (i in 1:(res$ntree)) { Rule_Table <- getTree(res, k = i, labelVar = TRUE) cat(paste0("INSERT INTO tbl_rf \nSELECT\n id, ")) Make_SQL(Rule_Table, 1, 1) cat(paste0(" as tree", i, "\nFROM\n", " input_data", ";\n\n")) } cat(paste0("INSERT INTO rf_predictions\n", "SELECT\n a.id,\n a.pred\n", "FROM\n (\n SELECT id as id, pred, COUNT(*) as cnt \n", " FROM tbl_rf\n GROUP BY id, pred\n ) a\n", " INNER JOIN\n (\n ", "SELECT id, MAX(cnt) as cnt\n", " FROM\n", " (\n", " SELECT id as id, pred, COUNT(*) as cnt\n", " FROM tbl_rf\n", " GROUP BY id, pred\n", " )\n", " GROUP BY id\n ) b\n", " ON a.id = b.id AND a.cnt = b.cnt;\n\n"))
最後に
Random ForestはDeep Learningが注目されるようになった現在においても強力なアルゴリズムの一つである。またDeep Forestのように多段構造とすることでDeep Learningを上回る性能を得る場合もある。分散処理に強く、人による解釈も容易であるためビジネスにおける実用性は相当に高いアルゴリズムであるため、実装への障害によって敬遠されるのは非常に勿体無いと思う。この記事がそれらの障害を乗り越えるための一助となれば幸いである。
現状ではDeep LearningよりもRandom Forestを使いこなせる方が実用においては有益だろう。曲線的な分類に弱いなどの欠点もあるが、まだまだ発展のポテンシャルが高いアルゴリズムだと思うので今後も注目していきたい。
ggvisによるRで可視化
面白いパッケージを発見したので紹介。可視化のための機能を提供してくれるもので、{ggvis}というもの。もしかすると{ggplot2}を超えるかもしれない。
使い方が非常にシンプルで、思想として{ggplot2}に似ているので慣れた人にはスイッチしやすいと思う。dplyrユーザーにはお馴染みの%>%演算子を使用可能というのもよい。
ひとまず試してみる。
## Library installation install.packages("ggvis") ## Load library library(ggvis)
ggvisの実行。
mtcars %>% ggvis(~mpg, ~wt, fill = ~factor(cyl)) %>% layer_points()
このように可視化に用いるデータを%>%で後続の処理に渡し、後から可視化に用いるグラフを指定する形式は{ggplot2}と同じ。
では{ggvis}の何が面白いかというと、Shinyのようなインタラクティブな操作が可能なグラフを以下のように簡単に作成できること。
以下をお試しあれ。
data.frame(x = rnorm(1000)) %>% ggvis(~x, fill := input_select(c("red", "blue", "green"))) %>% layer_histograms(width = input_slider(0.01, 3, 0.1, 0.01, label="bin_width"))
連続量のヒストグラムを作成するとき、ビンの幅を決めるのは以外に難しい。R標準のパッケージである{graphics}のhistや、{MASS}のtruehistは所定のアルゴリズム（ScottとかFDとか）に基づいて自動的に決定してくれるが、実際にデータを分析する立場からするとビンを様々に変えながら眺めたいところ。
そのような時にShinyのようにインタラクティブにビン幅を変えられると嬉しいのだが、Shinyは準備に手間がかかるのでなかなか手が出ない。それを{ggvis}では先ほどのように簡単なコマンドでインタラクティブなグラフを作成できる。
他にも
mtcars %>% ggvis(x = ~wt, y = ~mpg) %>% layer_points() %>% layer_smooths(span = input_slider(0, 1, 0.5, 0.1, label = "span"))
このような使い方ができる。
このパッケージはまだ開発の途上のようで、demoで見ようとしたいくつかの機能がエラーで見れなかったり、作成したグラフが保存できなかったりする（だからこのブログでもグラフの画像を貼っていない）。
今後の修正に大きく期待できる、面白いパッケージだと思う。
Anscombeの例
二変数間の関連性の強さを評価する上で最も良く使われる指標は、間違いなくPearsonの積率相関係数、いわゆる相関だと思う。Excelでも関数により算出でき、関連の強さを-1〜1で評価できるというお手軽さが受け入れられやすいのだと思う。
さて、この相関は直線的な関係性しか評価できないため、必ず散布図とセットで見なくてはならない。表題の「Anscombeの例」とは、統計量が同じであってもデータセットが示す内容は大きく異なることがあるということを明確に教えてくれる良い教材である。実際に見てみよう。
## Load libraries library(tidyverse) library(broom) library(forcats) library(ggplot2)
まずはRに標準で格納されているデータセットanscombeを使う。
> anscombe x1 x2 x3 x4 y1 y2 y3 y4 1 10 10 10 8 8.04 9.14 7.46 6.58 2 8 8 8 8 6.95 8.14 6.77 5.76 3 13 13 13 8 7.58 8.74 12.74 7.71 4 9 9 9 8 8.81 8.77 7.11 8.84 5 11 11 11 8 8.33 9.26 7.81 8.47 6 14 14 14 8 9.96 8.10 8.84 7.04 7 6 6 6 8 7.24 6.13 6.08 5.25 8 4 4 4 19 4.26 3.10 5.39 12.50 9 12 12 12 8 10.84 9.13 8.15 5.56 10 7 7 7 8 4.82 7.26 6.42 7.91 11 5 5 5 8 5.68 4.74 5.73 6.89
anscombeにはこのように、4つのデータセットそれぞれがxとyを持つ。
このままだと少し使いにくいので整形する。
## Anscombe dataAns <- select(anscombe, "x" = x1, "y" = y1) %>% bind_rows(select(anscombe, "x" = x2, "y" = y2)) %>% bind_rows(select(anscombe, "x" = x3, "y" = y3)) %>% bind_rows(select(anscombe, "x" = x4, "y" = y4)) %>% mutate("dataset" = c(rep("I", nrow(anscombe)), rep("II", nrow(anscombe)), rep("III", nrow(anscombe)), rep("IV", nrow(anscombe))))
出来上がったデータに対して平均と標準偏差を求める。
dataAns %>% group_by(dataset) %>% summarise_each(funs(mean, sd), vars = c(x, y)) # A tibble: 4 × 5 dataset vars1_mean vars2_mean vars1_sd vars2_sd <chr> <dbl> <dbl> <dbl> <dbl> 1 I 9 7.500909 3.316625 2.031568 2 II 9 7.500909 3.316625 2.031657 3 III 9 7.500000 3.316625 2.030424 4 IV 9 7.500909 3.316625 2.030579
平均と標準偏差が小数点以下2位まで揃っている。それぞれのデータセットについて相関を求めてみると：
dataAns %>% group_by(dataset) %>% summarise("COR" = cor(x, y)) # A tibble: 4 × 2 dataset COR <chr> <dbl> 1 I 0.8164205 2 II 0.8162365 3 III 0.8162867 4 IV 0.8165214
これもほぼ同じ値となる。平均、標準偏差、相関が等しいデータセットなのでxとyも同じような値と取っているのだろうと思うのだが、プロットしてみると：
ggplot(dataAns, aes(x = x, y = y, label = dataset)) + facet_wrap(~ dataset, nrow = 2) + geom_point() + geom_smooth(method = "lm", formula = y ~ x, se = FALSE)
このように曲線を描くものや外れ値に影響を受けているものがあることがわかる。
なお相関が同じなので当然と言えば当然だが、回帰係数も同じような値となる：
dataAns %>% split(.$dataset) %>% map(., ~ lm(y ~ x, data = .) %>% tidy()) %>% bind_rows() term estimate std.error statistic p.value 1 (Intercept) 3.0000909 1.1247468 2.667348 0.025734051 2 x 0.5000909 0.1179055 4.241455 0.002169629 3 (Intercept) 3.0009091 1.1253024 2.666758 0.025758941 4 x 0.5000000 0.1179637 4.238590 0.002178816 5 (Intercept) 3.0024545 1.1244812 2.670080 0.025619109 6 x 0.4997273 0.1178777 4.239372 0.002176305 7 (Intercept) 3.0017273 1.1239211 2.670763 0.025590425 8 x 0.4999091 0.1178189 4.243028 0.002164602
このanscombeと同様に、統計量が同程度で全く内容が異なるデータセットとしてdatasaurusというものがある。これを用いて同様に見てみよう。
## Datasaurus datasaurus <- read_table2("../Data/DatasaurusDozen.tsv", col_names = TRUE, col_types = "cnn") %>% mutate(dataset = as_factor(dataset))
> print(datasaurus) # A tibble: 1,846 × 3 dataset x y <fctr> <dbl> <dbl> 1 dino 55.3846 97.1795 2 dino 51.5385 96.0256 3 dino 46.1538 94.4872 4 dino 42.8205 91.4103 5 dino 40.7692 88.3333 6 dino 38.7179 84.8718 7 dino 35.6410 79.8718 8 dino 33.0769 77.5641 9 dino 28.9744 74.4872 10 dino 26.1538 71.4103 # ... with 1,836 more rows
このデータをプロットしてみるとどうか。
ggplot(datasaurus, aes(x = x, y = y)) + facet_wrap(~ dataset, nrow = 3) + geom_point()
なんと、全く違う！
これらの全てが本当に同様の統計量を持っているのだろうか？
datasaurus %>% group_by(dataset) %>% summarise_each(funs(mean, sd), vars = c(x, y)) # A tibble: 13 × 5 dataset vars1_mean vars2_mean vars1_sd vars2_sd <fctr> <dbl> <dbl> <dbl> <dbl> 1 dino 54.26327 47.83225 16.76514 26.93540 2 away 54.26610 47.83472 16.76982 26.93974 3 h_lines 54.26144 47.83025 16.76590 26.93988 4 v_lines 54.26993 47.83699 16.76996 26.93768 5 x_shape 54.26015 47.83972 16.76996 26.93000 6 star 54.26734 47.83955 16.76896 26.93027 7 high_lines 54.26881 47.83545 16.76670 26.94000 8 dots 54.26030 47.83983 16.76774 26.93019 9 circle 54.26732 47.83772 16.76001 26.93004 10 bullseye 54.26873 47.83082 16.76924 26.93573 11 slant_up 54.26588 47.83150 16.76885 26.93861 12 slant_down 54.26785 47.83590 16.76676 26.93610 13 wide_lines 54.26692 47.83160 16.77000 26.93790
anscombe同様に、平均と分散が小数点以下2位まで同じとなっている。
このように、データから得られる統計量が同一であっても内容が同じ傾向を示すとは限らない。そのため統計量を求める場合には必ずデータをプロットすることが重要である。平均や分散ならヒストグラムを、相関なら散布図を描くことでデータの分布について大まかに把握することができるだろう。
データの可視化というものは特に分析の後半に入ってくると忘れられがちであるが、データの特性を教えてくれる重要な手法である。分析を実施する時にはまず可視化を行うことを大切にしよう。
Ad-Stock効果を推定しつつ回帰を回したい③
Marketing Mix Modelingにおける広告の累積効果の推定について、以下の記事を書いた。
ushi-goroshi.hatenablog.com
ushi-goroshi.hatenablog.com
本日はこの続きで、同じような条件で作成したデータに対してstanを用いたベイズ推定を実施したので、その内容を紹介する。
なお先日の記事でoptimによりAd-Stock効果の回帰係数および減衰率を十分に精度よく推定できそうなことはわかっているが、
- 将来的なモデルの高度化
- 多様な推定方法の確立（手持ちの武器を増やす）
ことを目的としている。
シミュレーションデータの作成 〜 stanによるフィッティングの実行
まずは事前準備として、必要なライブラリを読み込む。当然ながらstanを使うので、事前にインストールを済ませた上で{rstan}を読み込む。
## load libraries library(dplyr) library(tidyr) library(tibble) ## load rstan library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores())
続いてシミュレーションデータ作成用の関数を用意する。これは先日の記事で紹介したものと大きな変更はないが、拡張性を考慮して変数の順序を入れ替え、また説明変数の作成に用いていたrnormの代わりにrgammaを使用している。
set.seed(123) simulate_y <- function(pars) { 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 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") 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, "X_01_Fil" = X_01_fil, "X_02_Fil" = X_02_fil, "Y_lag" = dplyr::lag(y, 1), "True_Error" = error) return(dat) } init_par <- array(c(100, 5, 2, 0.5, 0.6, 0.8, 0.5)) dat_Ana <- na.omit(simulate_y(init_par))
シミュレーションデータが作成できたので、stanを実行するための準備を実行する。stanに渡すためのパラメータを定義し、事前に準備したstanのスクリプトを用いてMCMCによるフィッティングを行う。
サンプルは10,000とし、そのうち最初の5,000サンプルをwarm-upとして破棄した。この辺りの設定は結果を見ながらの試行錯誤によるもので、明確な根拠に基づくものではない。traceプロットやヒストグラムを見る限り、もう少しサンプルを減らしても問題ないように思う。
dat_Ord <- list(N = nrow(dat_Ana), Y = dat_Ana$Y, X_01 = dat_Ana$X_01, X_02 = dat_Ana$X_02) ## fitting fit_01 <- stan(file = './AdStock_Simulation.stan', data = dat_Ord, iter = 10000, chains = 4, seed = 1234, warmup = 5000)
stanコードの記述
stan()で参照しているstanのコードは以下のように記述した。dataブロックおよびparameterブロックは大した工夫もないが、減衰率を示すlambda_X_01およびlambda_X_02はパラメータスペースを(0, 1)に制限している。
このコードで特筆すべきはtransformed parametersブロックで、ここでlambda_X_0xの値に基づきAd-Stock変数を作成している。Rであればstats::filterで簡単に作成できるのだが、今回の分析では減衰率も推定の対象であるためどうしてもstan内で記述する必要があった。ここで作成したX_01_filおよびX_02_filをmodelブロックで用いている。
data { int N; vector[N] Y; vector[N] X_01; vector[N] X_02; } parameters { real mu; real<lower=0> var_error; real beta_X_01; real<lower = 0, upper = 1> lambda_X_01; real beta_X_02; real<lower = 0, upper = 1> lambda_X_02; } transformed parameters { vector[N] X_01_fil; vector[N] X_02_fil; X_01_fil[1] = X_01[1]; X_02_fil[1] = X_02[1]; ## Create Ad-stock variable for(j in 2:N) { X_01_fil[j] = X_01[j] + X_01_fil[j-1] * lambda_X_01; X_02_fil[j] = X_02[j] + X_02_fil[j-1] * lambda_X_02; } } model { ## Sampling for(i in 1:(N)) { Y[i] ~ normal(mu + beta_X_01 * X_01_fil[i] + beta_X_02 * X_02_fil[i], var_error); } }
結果
stanを実行した結果は以下のように取得できる。
## Plot pars <- c("mu", "var_error", "beta_X_01", "lambda_X_01", "beta_X_02", "lambda_X_02") print(fit_01) stan_trace(fit_01, pars = pars) stan_hist(fit_01, pars = pars) stan_dens(fit_01, separate_chains = T, pars = pars) stan_ac(fit_01, separate_chains = T, pars = pars)
結果を見てみよう。まずは推定されたパラメータの事後分布から。
> print(fit_01) Inference for Stan model: AdStock_Simulation. 4 chains, each with iter=10000; warmup=5000; thin=1; post-warmup draws per chain=5000, total post-warmup draws=20000. mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat mu 4.70 0.01 0.63 3.41 4.29 4.72 5.13 5.91 13079 1 var_error 1.45 0.00 0.11 1.26 1.37 1.44 1.52 1.68 17554 1 beta_X_01 0.50 0.00 0.08 0.35 0.44 0.50 0.55 0.65 15821 1 lambda_X_01 0.63 0.00 0.08 0.44 0.58 0.63 0.69 0.77 12174 1 beta_X_02 0.71 0.00 0.10 0.51 0.64 0.70 0.77 0.91 13692 1 lambda_X_02 0.56 0.00 0.07 0.41 0.52 0.57 0.61 0.69 14090 1 X_01_fil[1] 4.74 0.00 0.00 4.74 4.74 4.74 4.74 4.74 2 1 X_01_fil[2] 3.51 0.00 0.40 2.63 3.27 3.55 3.79 4.18 12174 1 X_01_fil[3] 4.94 0.00 0.52 3.87 4.59 4.96 5.30 5.91 12453 1 X_01_fil[4] 9.09 0.01 0.73 7.65 8.60 9.09 9.58 10.49 12477 1 … 省略 …
X_01_filは本来推定すべきパラメータではないが、stanでは値が決まっていない変数（確率変数と見なされるもの）は全てパラメータとして扱われるため、特に分布に関心はないものの結果として含まれている。
では実際に関心のあるパラメータに着目すると、muが設定した値からやや低く推定されているものの、50%ベイズ信頼区間には真の値が含まれている。またX_01の回帰係数および減衰率も同様に設定した値に近いものを得られている。
一方で残差分散およびX_02の回帰係数・減衰率については、大きく離れてはいないものの50%信頼区間からは外れてしまっている。特に残差分散は何度か試行を繰り返したものの、そのいずれにおいても過少に推定されてしまった。この原因についてはよくわからない。
なおX_02については、回帰係数および減衰率を高めに設定することでより収束しやすい傾向があった。
プロットにより作成されるグラフは、それぞれ以下の通り。
> stan_trace(fit_01, pars = pars)
> stan_hist(fit_01, pars = pars)
> stan_dens(fit_01, separate_chains = T, pars = pars)
> stan_ac(fit_01, separate_chains = T, pars = pars)
traceプロットを見ると安定した形となっており、Rhatの数値と合わせて考えても収束していると判断して良さそうだ。分布の形をchainごとで比較しても一部で歪んでいる様子は見当たらない。
なおヒストグラムを確認すると基本的には正規分布に近い形となっているが、減衰率については左に裾を引く形の分布となっている。
まとめ
上記のシミュレーションの結果をまとめると、
- stanを用いることでAd-Stock効果をそれなりに精度よく推定できる
- 減衰率のような非線形の効果を素直にモデルに取り込むことができ、また分布の形もKoyckに限定されることなく柔軟なモデリングが可能
- しかしパルダ・モデルやoptimによる推定と比較すると精度に劣り、また処理時間およびコーディングに要する時間も比較的多くなる
と言えそうだ。分析準備に要するコストは習熟度が上がれば削減されるだろうが、一回の分析に要する処理時間が相対的に大きいため試行錯誤する場合には向いていないと思われる。
一方でベイズモデリングの大きなメリットである「モデルの柔軟性」は非常に魅力的に感じ、Marketing Mix Modelingの高度化を目指すに当たって積極的に使っていきたい技術である。
上記を踏まえてこれまでの分析結果をまとめると、
- Ad-Stockの分布としてKoyck型を仮定でき、残存の規模（減衰率）に関心はない場合 → パルダ・モデル
- 残存の規模や分布に関心があり、試行錯誤が必要 → optimによる推定
- より柔軟なモデリング（階層ベイズや状態空間モデルへの拡張）を実施したい → stanによるモデリング
といった使い分けが良さそうだ。なおstanによるモデリングは分析者のアイディアおよび対象領域に対する知見の深さにかなり依存すると思われるため、少なくともMMM構築の初期では全く適していないだろう。
逆に言えば、当該領域について十分な経験を持った分析者がいるのであれば、アイディアの限りを用いて色々なモデルを試すことが出来るため大きな価値を提供することが出来るようになるだろう。
最後に、この分析に限らずMCMCによるもの全般に言えることだが、やはりパラメータの分布という形で結果を見せることができるというのは大きな強みであると感じた。統計解析に馴染みのない相手であっても、「パラメータがバラつきを持ち、その範囲を分布として示す」というのは直感的にわかりやすいだろうと思う。頻度論的なアプローチであってもパラメータの（誤差）分布を標準誤差として示すことは出来るが、理解の容易さという点では比較にならないだろう。
唯一の懸念は、「パラメータが分布を持つ」という仮定に耐えうるか、ということだろうか。
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による推定の方が適当であると考えられる。
したがって分析者が必要とする数値および対象領域における仮定を踏まえて推定方法を選択する必要があるだろう。
