マーケターのためのRobynを用いた自動Marketing Mix Modeling (MMM)

Meta(旧Facebook)社の開発したAutomated Marketing Mix Modeling (MMM)ツールの一つRobynの紹介と実践をします。

Contents

1 MMMの復習
2 Robyn実践
3 まとめ

MMMの復習

以前こちらでMMMについて扱いました。改めて復習をします。
MMMはマーケティングの費用対効果を定量化する手法の総称とされています。具体的には収益に与える各メディアの効果、それ以外の効果をそれぞれ定量化するなどです。下図はMMMのイメージです。Robynの公式HPより引用しています。

Robyn tutorial : https://facebookexperimental.github.io/Robyn/docs/analysts-guide-to-MMM

MMMの実践は複数のステップからなります¹。

MMM ステップ

MMMでは成果変数(売上やコンバージョンなど)への寄与要因を以下のように仮定しています。

大項目	項目
ペイドメディア	4マス(TV/ラジオ新聞/雑誌)
	ネット(google/Yahoo)
	SNS(Facebook/Twitter/Instagram)
非ペイドメディア	オウンドメディア
	プレスリリース
	SNSでの発信
イベント変数	キャンペーンなど
	競合キャンペーン
	スポーツや音楽イベントなど
ベースライン	トレンド
	季節性
	祝休日
	気候情報

以下は入力データの一例です。

Robyn実践

RobynはMeta(Facebook)のメンバーが作成した自動MMMツールです。本稿では簡単なデータを用いて広告効果の定量化と予算の最適化を実践していきます。色々な機能がありますが、今回は上記に必要な機能に絞って解説します。また余裕があれば別記事でオプション機能の紹介をしたいと思います。なお本実践は主に、meta公式チュートリアルを参考にしています。

注意事項) WordPressの仕様でRの代入[<-]記号が[<-]に化けてしまっています。読み替えをお願いします。

Rをインストール

この手のライブラリでは珍しく、ラッパーがRしかないとのこと。Rの最新版をここからDlしてインストール。 ²

Robynのインストールとそれに必要な準備

まず作業ディレクトリ(ここでは./MMM_01を作成しそこを作業ディレクトリとした)に対象となるデータを置いておきます。次に以下のコードをRのコンソール(またはRstudio³)上から入力し実行していきます。

# 作業ディレクト
setwd("./MMM_01")
# remotes　インストール
install.packages('remotes')
# robyn install
remotes::install_github("facebookexperimental/Robyn/R")

# 作業ディレクト

setwd("./MMM_01")

# remotes　インストール

install.packages('remotes')

# robyn install

remotes::install_github("facebookexperimental/Robyn/R")

Robynはモデリングに際してNevergra⁴という最適化のソルバーを使用します。そのためにRからPythonを使用できるようにしておきます。Pythonを使用するためにはreticulate⁵というライブラリを使います。

# Nevergradをpythonから利用するために、reticulateをインストールする。
# reticulate
install.packages("reticulate")

#dplyr インストール
install.packages("dplyr")
library(dplyr)

library(reticulate)

# Nevergradをpythonから利用するために、reticulateをインストールする。

# reticulate

install.packages("reticulate")

#dplyr インストール

install.packages("dplyr")

library(dplyr)

library(reticulate)

nevergradを実行するpython環境を作成します。

virtualenv_create("r-reticulate")
py_install("nevergrad", pip = TRUE)
use_virtualenv("r-reticulate", required = TRUE)

virtualenv_create("r-reticulate")

py_install("nevergrad", pip = TRUE)

use_virtualenv("r-reticulate", required = TRUE)

作成した環境下にあるpythonの実行ファイルのパスを以下のように指定(隠しフォルダなため注意)。

use_python("~/.virtualenvs/r-reticulate/bin/python")

1	use_python("~/.virtualenvs/r-reticulate/bin/python")

#Robyn をインストール
library(Robyn)

1 2	#Robyn をインストール library(Robyn)

入力データ

入力データは、筆者が以前書籍⁶執筆用に用意したデータを用います。人工的に作成した広告の投下費用や売上などからなる日次の時系列データです。

変数	内容	仮説
日付	日時	–
ネット広告出稿金額	その日のネット広告出稿金額	売上に比例
TV広告出稿金額	その日のTV広告出稿金額	売上に比例
平均気温	その日の平均気温	売上に比例
降水量	その日の降水量	売上に反比例
売上金額	その日の売上げ	–
週末FLG	週末であるならば1それ以外0の変数	1の場合売上増加

データ読み込み

# 入力データのcsvを読み込む
dt_simulated = read.csv(file='sample.csv')
# change var names 日本語はplotした時に文字化けするため半角アルファベットに変換
dt_simulated = dplyr::select(dt_simulated,datetime=1, Net.Spend=2 ,Tv.Spend=3 ,temperature=4,
rain=5, revenue=6, Weekend.FLG=7,dplyr::everything())
head(dt_simulated)

#祝日用のデータも読み込み
data("dt_prophet_holidays")
head(dt_prophet_holidays)

# 入力データのcsvを読み込む

dt_simulated = read.csv(file='sample.csv')

# change var names 日本語はplotした時に文字化けするため半角アルファベットに変換

dt_simulated = dplyr::select(dt_simulated,datetime=1, Net.Spend=2 ,Tv.Spend=3 ,temperature=4,

rain=5, revenue=6, Weekend.FLG=7,dplyr::everything())

head(dt_simulated)

#祝日用のデータも読み込み

data("dt_prophet_holidays")

head(dt_prophet_holidays)

祝日用のデータはhttps://github.com/facebookexperimental/Robyn/tree/main/R/data
からダウンロードしておきます。

最後に出力先フォルダのpathを設定します。

#create object path
robyn_object &lt;- "[作業者の任意のフォルダのパス]/MyRobyn.RDS"

1 2	#create object path robyn_object <- "[作業者の任意のフォルダのパス]/MyRobyn.RDS"

input定義

#input
InputCollect &lt;- robyn_inputs(
dt_input = dt_simulated #入力する元データ
,dt_holidays = dt_prophet_holidays #祝日データ
,date_var = "datetime" # 以下のフォーマットで"2020-01-01"
,dep_var = "revenue" # 売上高やCVなどの従属変数
,dep_var_type = "revenue" # 売上高かCVフラグか
,prophet_vars = c("trend", "season") # "trend","season", "weekday" &amp; "holiday"
,prophet_country = "US"# 国名 祝日 デフォルトで日本がないため一旦USとしておく
,context_vars = c("temperature", "rain","Weekend.FLG") # イベント情報
,paid_media_spends = c("Net.Spend","Tv.Spend") # メディア投下
,paid_media_vars = c("Net.Spend","Tv.Spend") # メディア
# paid_media_vars must have same order as paid_media_spends. Use media exposure metrics like
# impressions, GRP etc. If not applicable, use spend instead.
#,organic_vars = c("newsletter") # PRなどの非広告メディア
,factor_vars = c("Weekend.FLG") # イベント因子
,window_start = "2019-04-01" # モデル構築に使用するデータの開始日
,window_end = "2020-03-31" # モデル構築に使用するデータの終了日
,adstock = "geometric" # adstockの形状
)
print(InputCollect)

#input

InputCollect <- robyn_inputs(

dt_input = dt_simulated #入力する元データ

,dt_holidays = dt_prophet_holidays #祝日データ

,date_var = "datetime" # 以下のフォーマットで"2020-01-01"

,dep_var = "revenue" # 売上高やCVなどの従属変数

,dep_var_type = "revenue" # 売上高かCVフラグか

,prophet_vars = c("trend", "season") # "trend","season", "weekday" & "holiday"

,prophet_country = "US"# 国名祝日デフォルトで日本がないため一旦USとしておく

,context_vars = c("temperature", "rain","Weekend.FLG") # イベント情報

,paid_media_spends = c("Net.Spend","Tv.Spend") # メディア投下

,paid_media_vars = c("Net.Spend","Tv.Spend") # メディア

# paid_media_vars must have same order as paid_media_spends. Use media exposure metrics like

# impressions, GRP etc. If not applicable, use spend instead.

#,organic_vars = c("newsletter") # PRなどの非広告メディア

,factor_vars = c("Weekend.FLG") # イベント因子

,window_start = "2019-04-01" # モデル構築に使用するデータの開始日

,window_end = "2020-03-31" # モデル構築に使用するデータの終了日

,adstock = "geometric" # adstockの形状

)

print(InputCollect)

今回行うデータ特有の注意事項

今回説明変数にweekendFLGがあるため、prophet_varsでは[weekday]という項目を外している。
日本の祝日がデフォルトでは使用できないため、prophet_varの設定で[holiday]を外している。日本の祝日の使い方は別途で解説するつもりである。

そのほかの細かいinputの細かい定義は公式HPをご参照してください。(どこかで整理するかも)

Prophet⁷に関連した変数があることからわかるように、backendとしてProphetを一部使用しているようです。Ptophetについての記事はこちらをご参照ください。

ハイパーパラメータ指定

ハイパーパラメータは広告効果の逓減(saturation)を制御するパラメータとAdstockと呼ばれる広告効果の減衰を制御するパラメータの2つの種類があります。前者のパラメータは[alpha]と[gamma]です。後者は減衰のタイプが幾何型とワイブル型があり、それぞれ幾何型が[theta]を、ワイブル型が[scale][shape]というパラメータを持ちます。

使用するハイパーパラメータを表示し、その可変領域を設定していきます。この例では幾何型(geometric)の減衰を指定したので、設定するのはthetaとalpha/gammaの3つとメディアがTVとNetの2つで計6つです(2*3=6)。

#使用するハイパーパラメータ表示
hyper_names(adstock = InputCollect$adstock, all_media = InputCollect$all_media)

#上記のハイパーパラメータの可変領域を設定
hyperparameters &lt;- list(
Net.Spend_alphas = c(0.5, 3)
,Net.Spend_gammas = c(0.3, 1)
,Net.Spend_thetas = c(0, 0.3)

,Tv.Spend_alphas = c(0.5, 3)
,Tv.Spend_gammas = c(0.3, 1)
,Tv.Spend_thetas = c(0.1, 0.4)
)

InputCollect &lt;- robyn_inputs(InputCollect = InputCollect, hyperparameters = hyperparameters)
print(InputCollect)

#使用するハイパーパラメータ表示

hyper_names(adstock = InputCollect$adstock, all_media = InputCollect$all_media)

#上記のハイパーパラメータの可変領域を設定

hyperparameters <- list(

Net.Spend_alphas = c(0.5, 3)

,Net.Spend_gammas = c(0.3, 1)

,Net.Spend_thetas = c(0, 0.3)

,Tv.Spend_alphas = c(0.5, 3)

,Tv.Spend_gammas = c(0.3, 1)

,Tv.Spend_thetas = c(0.1, 0.4)

)

InputCollect <- robyn_inputs(InputCollect = InputCollect, hyperparameters = hyperparameters)

print(InputCollect)

plot_adstock(plot = TRUE)
plot_saturation(plot = TRUE)

1 2	plot_adstock(plot = TRUE) plot_saturation(plot = TRUE)

とするとハイパーパラメータごとに減衰やサチュレーションの様子を見ることができます。

上記はオプションなので、グラフを表示させる必要がない場合は「plot=FALSE」にしてください。

(参考) 各種定義の詳細⁸

モデル式

Ridge回帰モデル

Adstockやsaturationの式

幾何型Adstockは decay_t,j := x_t,j + θ * decay_t-1,j
ワイブル型は decay_t,j := 1-(exp(decay_t-1,jj/α))
saturationは c * (x^α / (x^α + γ^α )) cはメディアxに対して得られた回帰係数

モデル構築

モデルを構築していきます。マニュアルでは初期モデルの構築とあるように、この後追加データがある場合それを加えて初期モデルを修正することもできます。

https://facebookexperimental.github.io/Robyn/docs/demo-R-script にある方法とDemoのhttps://github.com/facebookexperimental/Robyn/blob/main/demo/demo.R の方法ではやり方が異なるようで、ここではdemo.Rの方を参考にしています。

## モデル構築 実行
OutputModels &lt;- robyn_run(
InputCollect = InputCollect # feed in all model specification
#, cores = NULL # default
#, add_penalty_factor = FALSE # Untested feature. Use with caution.
, iterations = 2000 # recommended for the dummy dataset
, trials = 5 # recommended for the dummy dataset
, outputs = FALSE # outputs = FALSE disables direct model output
)
print(OutputModels)

## モデル構築実行

OutputModels <- robyn_run(

InputCollect = InputCollect # feed in all model specification

#, cores = NULL # default

#, add_penalty_factor = FALSE # Untested feature. Use with caution.

, iterations = 2000 # recommended for the dummy dataset

, trials = 5 # recommended for the dummy dataset

, outputs = FALSE # outputs = FALSE disables direct model output

)

print(OutputModels)

2000回のiterationsを行う試行を5回ここでは行います。結構時間がかかります。

厳密な説明は難しいのですが、ここでのモデリングはリッジ回帰をモデリングしているようで、ハイパーパラメータと通常のパラメータを同時に学習しているようです。この時Prophetに関係する部分のlossと、通常のリッジ回帰の部分のlossに対してMulti-objective optimization⁹を行なっています。Nevergradはこの最適化に対して用いられているようです。

最後に作成したモデル群の中からパレート最適な物を選ぶ処理?を行うようです。正直ここは少し理解できていません。

OutputCollect &lt;- robyn_outputs(
InputCollect, OutputModels
, pareto_fronts = 1
# , calibration_constraint = 0.1 # range c(0.01, 0.1) &amp; default at 0.1
, csv_out = "pareto" # "pareto" or "all"
, clusters = TRUE # Set to TRUE to cluster similar models by ROAS. See ?robyn_clusters
, plot_pareto = TRUE # Set to FALSE to deactivate plotting and saving model one-pagers
, plot_folder = robyn_object # path for plots export
)
print(OutputCollect)

OutputCollect <- robyn_outputs(

InputCollect, OutputModels

, pareto_fronts = 1

# , calibration_constraint = 0.1 # range c(0.01, 0.1) & default at 0.1

, csv_out = "pareto" # "pareto" or "all"

, clusters = TRUE # Set to TRUE to cluster similar models by ROAS. See ?robyn_clusters

, plot_pareto = TRUE # Set to FALSE to deactivate plotting and saving model one-pagers

, plot_folder = robyn_object # path for plots export

)

print(OutputCollect)

pareto_fronts = 1とするとパレート最適なフロンティアのベストソリューション(つまり2つのlossのトレードオフが最小になる)を選ぶようです。またclusters=TRUEとすると、ソリューション群の中でROIが似た者同士を同一のクラスタとしてまとめます。kは自動で決めてくれるようですが、３になる場合が多いです。

[print(OutputCollect)]の後に以下のようにソリューション群がリストアップされます。

はじめに定義したMyRobyn.RDSのフォルダパスと同じ場所に以下のファイルが出力されます。X_YYY_Z.pngの文字列からなるファイルはそのモデルのone-pager repotです。パレート最適なクラスタの代表モデルごとにこのファイルは存在するため3つあります¹⁰(つまりクラスタ数が5つならばこのファイルは5つあります)。

ne-pager report(ここでは上図の例として2_399_4の)を見てみます。

左上の図:ウォーターフォール図
成果変数への全体寄与を1とした場合の各要素の寄与割合を示した図です。一番寄与が高いのが「平均気温」、次が「trend」、その次が「週末かどうか」です。ここではTVの寄与はNetより大きいと示されています。メディアの寄与の比較はあくまで投下量ベースであるため、効率は考慮されていません。
右上の図:時系列折れ線グラフ
実際の値と予測値のプロット図です。概ね実績値と予測値が近い動きをしていることがわかります。
真ん中左の図:横棒グラフとROIの折れ線グラフ
メディアの実際の投下シェアと、成果変数のシェアのグラフです。赤い折れ線はROIです。全体比の寄与ではTVの方が大きいのですが、ROIではNetの方が大きいことが分かります。成果指標をconversionにすると、ROIの代わりにCPAが出力されるようです。
真ん中右の図:saturationのプロットグラフ
saturation効果をそれぞれのメディアで記述しています。Netの方が飽和しにくい、つまり効果の逓減が小さいことが示されています。
左下の図:広告効果の減衰グラフ
広告効果の減衰の程度をそれぞれのメディアごとに示しています。幾何的減衰(geometric)の場合、θの値そのもののようです。
右下の図:残差プロット
実績と予測値の残差プロットです。一般論として残差に偏りがあれば説明変数だけで説明できていない要因があると考えられます。

本ツールで算出しているROIの定義に関しては調査中ですが一般的な分子から投下量を控除する定義ではなく、成果/投下量で定義しているようです¹¹。

予算の最適化

最後に構築したモデル(2_399_4)で予算の最適化を行います。

# 2_339_4を選びます。
select_model &lt;- "2_399_4" # select one from above
robyn_save(robyn_object = robyn_object # model object location and name
, select_model = select_model # selected model ID
, InputCollect = InputCollect # all model input
, OutputCollect = OutputCollect # all model output
)

# 2_339_4を選びます。

select_model <- "2_399_4" # select one from above

robyn_save(robyn_object = robyn_object # model object location and name

, select_model = select_model # selected model ID

, InputCollect = InputCollect # all model input

, OutputCollect = OutputCollect # all model output

)

#選択したモデルのメディアサマリー
OutputCollect$xDecompAgg[solID == select_model &amp; !is.na(mean_spend)
, .(rn, coef,mean_spend, mean_response, roi_mean
, total_spend,
total_response=xDecompAgg, roi_total, solID)]

#選択したモデルのメディアサマリー

OutputCollect$xDecompAgg[solID == select_model & !is.na(mean_spend)

, .(rn, coef,mean_spend, mean_response, roi_mean

, total_spend,

total_response=xDecompAgg, roi_total, solID)]

予算最適の方法は2種類あり、”max_historical_response”と”max_response_expected_spend”です。前者は構築データと同じトータルメディア投下量で最大の売上を得る、メディアミックスを求めます。後者はメディア投下量を与えて、そのもとで最適な予算の編成を求めます。

“max_historical_response”の場合を実行してみます。”max_response_expected_spend”の公式チュートリアルを参照してください。

AllocatorCollect &lt;- robyn_allocator(
InputCollect = InputCollect
, OutputCollect = OutputCollect
, select_model = select_model
, scenario = "max_historical_response"
, channel_constr_low = c(0.01, 0.01) #メディア数と同じ長さが必要。元の投下量に対して何%まで変化を許容するかの下限
, channel_constr_up = c(10, 10) #メディア数と同じ長さが必要。元の投下量に対して何%まで変化を許容するかの上限

)
print(AllocatorCollect)
AllocatorCollect$dt_optimOut

AllocatorCollect <- robyn_allocator(

InputCollect = InputCollect

, OutputCollect = OutputCollect

, select_model = select_model

, scenario = "max_historical_response"

, channel_constr_low = c(0.01, 0.01) #メディア数と同じ長さが必要。元の投下量に対して何%まで変化を許容するかの下限

, channel_constr_up = c(10, 10) #メディア数と同じ長さが必要。元の投下量に対して何%まで変化を許容するかの上限

)

print(AllocatorCollect)

AllocatorCollect$dt_optimOut

先のone-pager reportと同じフォルダに、「2_399_4_reallocated_hist.png」というファイルができたはずです。これに今の最適化の結果が出力されています。

左上の図は最適化前と最適化後の広告費の構成比率を示しています。最適化後はほとんどネットに投下した方が良いという結果になっています。余談ですがこれは弊社標準モデルで計算した場合とほぼ同様の結果になっています¹²。
右上の図は最適化前と最適化後の売上への寄与のをそれぞれのメディアで比較しています。
下の図はそれぞれのメディアのsaturation曲線と、最適化前後の投下量のポイントを示しています。

まとめ

meta社のRobynを用いたAutomated MMMの実践を行いました。入力データとハイパーパラメータを与えることで売上とメディアの関係をモデリングすることができます。そして得られたモデルを用いて最適化予算のアロケーションを行うことが可能です。

更にここでご紹介した方法以外にも、追加データでモデルを構築するrefreshやおよモデルのcaribrationなどの機能があります。ここでご紹介できなかったものに関しては、今後どこかで解説記事を書きたいと思います。

一方でAutoと言いつつもかなり、細かいパラメータの設定を作業者が行う必要性があります。
Robynを始めとしたMMMに関してのご相談は是非弊社までお願いいたします。

当ツールは今でも改良が行われているようなので、最新版に関しての情報は公式HPをご確認ください。

投稿者プロフィール

株式会社Crosstab 代表漆畑充

当ブログは【アナリティクスやデータの総合的な活用を通じて「未知の知」の獲得に貢献する】株式会社Crosstabが作成しています。

筆者について
2007年より金融機関向けデータ分析業務に従事。与信及びカードローンのマーケテイングに関する数理モデルを作成。その後大手ネット広告会社にてアドテクノロジーに関するデータ解析を行う。またクライアントに対してデータ分析支援及び提言/コンサルティング業務を行う。
統計モデルの作成及び特にビジネスアウトプットを重視した分析が得意領域である。
お問い合わせは株式会社Crosstabまでお願いいたします。

マーケターのためのRobynを用いた自動Marketing Mix Modeling (MMM)

MMMの復習

Robyn実践

Rをインストール

Robynのインストールとそれに必要な準備

入力データ

データ読み込み

input定義

ハイパーパラメータ指定

(参考) 各種定義の詳細⁸

モデル式

Adstockやsaturationの式

モデル構築

予算の最適化

まとめ

投稿者プロフィール

最新の投稿

いいね:

関連

MMMの復習

Robyn実践

Rをインストール

Robynのインストールとそれに必要な準備

入力データ

データ読み込み

input定義

ハイパーパラメータ指定

(参考) 各種定義の詳細8https://facebookexperimental.github.io/Robyn/docs/variable-transformations

モデル式

Adstockやsaturationの式

モデル構築

予算の最適化

まとめ

投稿者プロフィール

最新の投稿

いいね:

関連

(参考) 各種定義の詳細⁸