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による推定の方が適当であると考えられる。
したがって分析者が必要とする数値および対象領域における仮定を踏まえて推定方法を選択する必要があるだろう。
Keras for R
RstudioがR上でKerasによるディープラーニングのモデルを構築するためのライブラリ{keras}を公開した。
以前から{tensorflow}を使えばtensorflow::import(module = "keras")でKerasを導入することができたようだが、{keras}を先にインストールすることでpythonさえ入っていればtensorflowのインストールもできるらしい。随分と楽にtensorflowを動かせることができそうなので早速試してみた。なお類似のプロジェクトがRstudio以外でも行われているようで、以下でもKerasをRから動かすためのライブラリを提供している。
実行環境はAWS上のWindowsマシン。Rのバージョンは3.3.2だけれど、Microsoft R Clientがインストールされている。
実行したスクリプトは上で紹介したページの内容なのだけれど、サンプルデータであるx_trainやx_testなどは自前で用意する必要があったのでirisを使って適当に作成した。なおこちらのページはKerasの公式ドキュメントをR用に書き換えただけの様子。
{Keras}のインストール
早速{keras}のインストールとtensorflowのインストールを実行してみる。なお私の環境では{keras}のインストールの際に{reticulate}が一緒に入らず、{keras}のインストールに失敗したため{reticulate}だけ先にインストールしている。
## library install & load devtools::install_github("rstudio/reticulate") devtools::install_github("rstudio/keras") library(keras) library(dplyr) # for data handling ## install tensorflow install_tensorflow() Error: Installing TensorFlow requires a 64-bit version of Python 3.5 Please install 64-bit Python 3.5 to continue, supported versions include: - Anaconda Python (Recommended): https://www.continuum.io/downloads#windows - Python Software Foundation : https://www.python.org/downloads/release/python-353/ Note that if you install from Python Software Foundation you must install exactly Python 3.5 (as opposed to 3.6 or higher).
{keras}のインストールがうまくいけばinstall_tensorflow()という関数一つでtensorflowをインストールすることができる。とはいえpythonのインストールまではやってくれないようで私の環境ではエラーが出てしまった。指示に従い、Anaconda経由でpythonをインストールする。なお現時点でのpythonのバージョンは3.6だった。
Anacondaのインストールが済んだので再度tensorflowのインストールを試してみる。
## install tensorflow install_tensorflow()
今度はうまくいった。
サンプルデータの作成
サンプルスクリプトでは分析用のデータが提供されていなかったため、irisを使ってサンプルデータを作成する。特別なことはしておらず、sampleによって乱数を振った上で学習用とテスト用にデータを分割した。なおデータはmatrixにする必要がある。
set.seed(123) learn_id <- sample(1:nrow(iris), 100, replace = FALSE) train_data <- iris %>% filter(row.names(.) %in% learn_id) test_data <- iris %>% filter(!row.names(.) %in% learn_id) ## training data x_train <- train_data %>% select(-Species) %>% as.matrix() y_train <- as.matrix(model.matrix(~ train_data$Species - 1)) ## test data x_test <- test_data %>% select(-Species) %>% as.matrix() y_test <- as.matrix(model.matrix(~ test_data$Species - 1))
Kerasによるmodelオブジェクトの生成
サンプルデータが出来上がったので、{keras}によるフィッティングを行う。まずはmodelオブジェクトの生成から。
## generate model object model <- keras_model_sequential() > Model Model _______________________________________________________________________________ Layer (type) Output Shape Param # =============================================================================== Total params: 0 Trainable params: 0 Non-trainable params: 0 _______________________________________________________________________________
入力層 → 隠れ層 → 出力層という簡単なニューラルネットワークを構築する。layer_denseで簡単にレイヤーを積み上げることができるようで、入力層のみinput_shapeを指定する必要がある。これはpython上でKerasを実行する際のmodel.addに相当するものと思われる。
model %>% layer_dense(units = 10, input_shape = 4) %>% # Output & Input dimension layer_activation(activation = 'relu') %>% layer_dense(units = 3) %>% layer_activation(activation = 'softmax') > model Model _______________________________________________________________________________ Layer (type) Output Shape Param # =============================================================================== dense_3 (Dense) (None, 10) 50 _______________________________________________________________________________ activation_3 (Activation) (None, 10) 0 _______________________________________________________________________________ dense_4 (Dense) (None, 3) 33 _______________________________________________________________________________ activation_4 (Activation) (None, 3) 0 =============================================================================== Total params: 83.0 Trainable params: 83.0 Non-trainable params: 0.0 _______________________________________________________________________________
パラメータの数は層ごとに(入力データの次元数+1) * 出力データの次元で決まる。+1は切片。
ここで出来上がったmodelオブジェクトの属性や構造などを見てみようとstrを実行してみたのだが、print(model)と結果は変わらなかった。やはりR的なオブジェクトではない様子。
ニューラルネットワークの構造が決まれば、最適化の際の条件を指定する。
## model parameter model %>% compile( loss = 'categorical_crossentropy', optimizer = optimizer_sgd(lr = 0.02), metrics = c('accuracy'))
フィッティング
以上で準備が整ったため、フィッティングを行う。
## fitting model %>% fit(x_train, y_train, epochs = 1000, batch_size = 32) ## 省略 # Epoch 990/1000 # 100/100 [==============================] - 0s - loss: 0.0432 - acc: 0.9800 # Epoch 991/1000 # 100/100 [==============================] - 0s - loss: 0.0442 - acc: 0.9800 # Epoch 992/1000 # 100/100 [==============================] - 0s - loss: 0.0455 - acc: 0.9800 # Epoch 993/1000 # 100/100 [==============================] - 0s - loss: 0.0508 - acc: 0.9800 # Epoch 994/1000 # 100/100 [==============================] - 0s - loss: 0.0602 - acc: 0.9800 # Epoch 995/1000 # 100/100 [==============================] - 0s - loss: 0.0437 - acc: 0.9800 # Epoch 996/1000 # 100/100 [==============================] - 0s - loss: 0.0433 - acc: 0.9800 # Epoch 997/1000 # 100/100 [==============================] - 0s - loss: 0.0432 - acc: 0.9800 # Epoch 998/1000 # 100/100 [==============================] - 0s - loss: 0.1129 - acc: 0.9600 # Epoch 999/1000 # 100/100 [==============================] - 0s - loss: 0.0529 - acc: 0.9700 # Epoch 1000/1000 # 100/100 [==============================] - 0s - loss: 0.0647 - acc: 0.9600
サンプルデータがものの100行程度しかないため一瞬で終わる。Accuracyが0.96とまずまず。続いてこちらのモデルを用いて、テストデータに対する精度評価を実行する。
## test loss_and_metrics <- model %>% evaluate(x_test, y_test, batch_size = 128) > loss_and_metrics [[1]] [1] 0.1853937 [[2]] [1] 0.96
テストデータに対しても精度が0.96と高く、過学習を起こしている様子はない。
テストデータに対するラベルの予測は以下のように実行する。
## predict classes <- model %>% predict(x_test, batch_size = 128) > classes [,1] [,2] [,3] [1,] 9.983064e-01 1.693639e-03 6.461022e-20 [2,] 9.971897e-01 2.810343e-03 7.951551e-19 [3,] 9.983315e-01 1.668578e-03 3.816257e-20 [4,] 9.984826e-01 1.517428e-03 4.521292e-20 [5,] 9.994913e-01 5.087039e-04 3.988174e-23 ## 省略
サンプルスクリプトは以上である。
注意点として、パイプ演算子(%>%)が使えるからといってdplyrが読み込まれている訳ではないので、データ加工を行うのであれば別途読み込む必要がある。
Rでディープラーニングを簡単に実行できる環境が整ってきているので、これを機会にRユーザーが増えると良いな。
帰無仮説は採択できない
統計的仮説検定は通常、以下のような手順に従って行われる：
- 帰無仮説を設定する
- 対立仮説を設定する
- 検定統計量を設定する
- 有意水準を設定する
- 実験やデータ解析によって検定統計量を求める
- 帰無仮説を真と仮定した時の検定統計量の得られる確率を求める
- 設定した有意水準に従い、判断を下す
このうち７番で行われる判断において、帰無仮説が棄却された場合は「対立仮説を採択する」という表現がなされる。では対立仮説が棄却された場合は「帰無仮説を採択する」と言うのかといえばそうではない。帰無仮説は採択できないからだ。
そもそも統計的仮説検定は「帰無仮説が真であるとき、この統計量が偶然に得られる確率はどの程度であろうか？」を定量的に評価するための手続きである。大事なポイントとしては「帰無仮説が真である」という点で、この前提の上でどの程度の検定統計量が得られれば帰無仮説を棄却できるか、ということを見ている。帰無仮説と対立仮説は排反であるため、帰無仮説が棄却された場合は自動的に対立仮説が採択される。
では帰無仮説が棄却されなかった場合はと言うと、これがややこしいところで、「帰無仮説を採択する」ことはできないのである。「棄却できない」と「採択する」は異なるからだ。
先にも述べたように、統計的仮説検定はあくまで「帰無仮説を棄却できるか」を検証するための手続きである。この手続きから得られる結果は「帰無仮説を棄却できる」「できない」の２つだけであり、採択に関する考慮は一切ない。そのため、例えば回帰分析を行った際にβ = 1.5, P < 0.01という回帰係数が得られたとしても、ここで言えることは「回帰係数が0であること（=帰無仮説）が棄却されたので、回帰係数は0ではない」ということだけであり、1.5という数値そのものに対しては何も言えないのである。
ちなみに検定統計量（t値やχ2値など）の大小をもって「極めて有意」「わずかに有意」などといった定量的は判断が下されることがあるが、厳密にはこれはよろしいものではない。伝統的な手続きに則った場合、棄却できるか否か（1か0か）の判断だけが可能なのである。
（念のため付け加えると、もちろん「t値が大きく、めったに起こらないことが発生した」といった表現は問題なく、「棄却できるか否か」について定量的な表現がよくない、ということである。最初に設定した有意水準を超えるか否かだけが焦点であり、5%と設定したのであれば、P値が0.04であろうと0.000001であろうと等しく「棄却できる」のである。）
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"を指定すれば上限・下限を指定できたと思うので、もう少し検証してみよう。
