🐥

Metaが開発したMMM(マーケティング・ミックス・モデリング)ツール「Robyn」は数学的に何をやっているのか?

2025/02/06に公開

Robynとは?

RobynとはMetaが開発した、オープンソースのMMM(マーケティング・ミックス・モデリング)パッケージです。Robynに関しては以下の記事をご覧ください。

全体概要

Robynでの解析は以下の手順で進めます。

  1. 必要データや広告効果の数理的表現
  2. 大まかな多変量回帰モデルの構築
  3. 初期モデル構築(リッジ回帰による正則化)
  4. 進化的アルゴリズムによる多目的最適化
  5. 最終モデル選択
  6. 広告予算配分の最適化

必要データや広告効果の数理的表現

【定義】時系列データ

MMMでは、時系列データを扱うため、時系列データを以下のように定義する。

  • t \in \mathbb{N}を時間(年・月・週・日など)のインデックスとする
  • 時系列データ\{y_{t}\}_{t=1}^{T}は、時刻tにおける売上やコンバージョンなどのKPIを表す。
  • 広告出稿量データ\{X_{t,j}\}_{t=1}^{T}は、時刻tにおける広告媒体jの出稿量を表す。

以上によりMMMの回帰モデルのデータセットを次のようにする。

\mathcal{D} = \{ (y_{t}, X_{t,1}, \ldots, X_{t,J}, Z_{t,1}, \ldots, Z_{t,K}) \}_{t=1}^{T}

【定義】アドストック変換

時系列データ\{X_{t}\}_{t\in\mathbb{N}}に対し、任意のtにおけるアドストック変換後の値X_{\text{decay}, t}を以下の漸化式により定義する:

(1)一般的なアドストック変換

X_{\textrm{decay}, t} = X_{t} + \theta \cdot X_{\textrm{decay}, t-1}

ここで、

  • X_{t}は時刻tにおける広告出稿料(またはインプット変数)
  • X_{\text{decay},t}は時刻tにおけるアドストック変換後の値
  • \theta \in [0, 1]減衰率(Decay Rate) を表す定数である。

(2) 幾何アドストック変換(Geometric Adstock)

X_{\textrm{decay}, t} = \sum_{k=0}^{\infty} \theta^{k} X_{t-k}

これは自己回帰型の加重平均を表し、過去の広告出稿が指数関数的に減衰することを意味する。 またこの幾何アドストック変換の特性として次が挙げられる。

  • \thetaが大きいほど過去の広告の影響が長く続く。
  • \thetaが0に近づくと、現在の広告の影響のみが残る。

(3) ワイブル分布型アドストック変換(Weibull CDF Adstock)

\lambda, k < 0を定数として、ワイブル累積分布関数(CDF)を用いたアドストック変換は以下の形で定義される:

w(k, \lambda, t) = 1 - e^{-(t/\lambda)^{k}}
X_{\textrm{decay}, t} = \sum_{k=0}^{\infty} w(k,\lambda,t-k) X_{t-k}

このときパラメータ\lambda, kは次のような意味を持つ。

  • \lambda > 0 : 減衰速度を制御するスケールパラメータ
  • k > 0: 減衰の形状を決める形状パラメータ(Shape Parameter)
    • 0 < k < 1: 初期の影響が強く、急激に減衰
    • k > 1: 初期の影響は弱く、徐々に減衰

ワイブル分布型アドストック変換の特性として次が挙げられる。

  • パラメータ\lambda, k < 0により、幾何アドストック変換よりも柔軟に減衰形状を持つ。
  • k < 1なら急速な減衰、k>1なら初期の影響が小さく後半で遅く減衰する。

また本記事におけるアドストック(Adstock)とは、広告の遅延効果を表す概念であり、広告効果が時間とともに減衰することを示している。

【定理】アドストック変換の収束

定理1(幾何アドストック変換の収束)

|\theta| < 1かつX_{t}が適切な範囲で抑えられているとき、幾何アドストック変換は収束する。以下に証明を与える。

幾何アドストックの漸化式を展開すると次のようになる。

X_{\text{decay}, t} = X_{t} + \theta X_{t-1} + \theta^2 X_{t-2} + \cdots

仮定より、X_{t}が定数Xに収束するので、t \rightarrow \inftyのとき幾何アドストック変換は次のようになる。

\begin{align*} X_{\textrm{decay}, \infty} &= \sum_{k=0}^{\infty} \theta^{k} X_{t-k} \\ &= \sum_{k=0}^{\infty} \theta^{k} X \\ &= X \sum_{k=0}^{\infty} \theta^{k} \end{align*}

ここで幾何級数の和の公式より

\sum_{k=0}^{\infty} \theta^{k} = \frac{1}{1-\theta} \hspace{10pt} (ただし |\theta | < 1)

が成り立つので、

X_{\textrm{decay}, \infty} = X \sum_{k=0}^{\infty} \theta^{k} = \frac{X}{1-\theta}

これは有限値に収束する(証明終)

定理2(ワイブル分布型アドストック変換の収束)

ワイブル分布型アドストック変換が収束することを示す。

まずはワイブル重み関数w(k,\lambda,t)の極限を考える。

\lim_{t \rightarrow \infty} w(k,\lambda,t) = \lim_{t \rightarrow \infty} \left( 1 - e^{-(\frac{t}{\lambda})^{k}} \right)

となるが、ここでe^{-(\frac{t}{\lambda})^{k}}の項に注目すると、t \rightarrow \inftyk > 0のとき以下のようになる。

e^{-(\frac{t}{\lambda})^{k}} \rightarrow 0

したがって、ワイブル重み関数全体で見たときに次が成り立つ。

\lim_{t \rightarrow \infty} w(k,\lambda,t) = \lim_{t \rightarrow \infty} \left( 1 - e^{-(\frac{t}{\lambda})^{k}} \right) = 1

したがって、有限値に収束する(証明終)

【定義】S-Curve関数変換

時系列データ\{X_{t}\}_{t \in \mathbb{N}}に対し、広告出稿量をX_{t}としたときのSカーブ関数変換(S-Curve transformation)は以下のように定義する:

S\text{-Curve}(X_{t}) = \beta \times \frac{X_{t}^{\alpha}}{X_{t}^{\alpha} + \gamma^{\alpha}}, \hspace{10pt} X_{t} \geq 0, \hspace{10pt} \alpha, \gamma > 0

ここで、

  • S\text{-Curve}(X_{t})は変換後の広告効果
  • \beta > 0は広告の最大効果(スケールパラメータ)
  • \alpha > 0は曲線の傾きを制御するシェイプパラメータ
  • \gamma > 0は広告効果が50%になる閾値(閾値パラメータ)
  • X_{t}は時刻tにおける広告出稿量

この関数は、広告出稿量X_{t}が増加するにつれて効果が頭打ちになる性質を持つ。

【性質】S-Curve関数変換の特性

S-Curve関数変換は次のような性質を持つ。

  1. 単調増加性
    (広告投入量が増加すると、効果も増加する)

    X_{1} \leq X_{2} \hspace{10pt} \Rightarrow \hspace{10pt} S\text{-Curve}(X_{1}) \leq S\text{-Curve}(X_{2})
  2. 漸近的な飽和効果
    (広告が少ないと効果も小さく、広告が非常に大きくなると最大効果\betaに収束する)

    \lim_{X_{t} \rightarrow 0} S\text{-Curve}(X_{t}) = 0, \hspace{10pt} \lim_{X_{t} \rightarrow \infty} S\text{-Curve}(X_{t}) = \beta
  3. 効果の50%閾値
    (広告出稿量が\gammaのとき、効果は最大値の50%になる)

    S\text{-Curve}(\gamma) = \frac{\beta}{2}
  4. 凹凸変化点

    • X_{t} < \gammaのとき、S-Curveは凸関数(加速的増加)
    • X_{t} > \gammaのとき、S-Curveは凹関数(鈍化する増加)
    • X_{t} = \gammaのとき、曲率が最大。

【定理】S-Curve関数の収束性

\lim_{X_{t}\rightarrow0}S\text{-Curve}(X_{t}) = 0, \hspace{10pt} \lim_{X_{t}\rightarrow \infty}S\text{-Curve}(X_{t}) = \beta

以下に証明を与える。

  • X_{t} \rightarrow 0のとき:

    S\text{-Curve}(X_{t}) = \beta \cdot \frac{X_{t}^{\alpha}}{X_{t}^{\alpha} + \gamma^{\alpha}} \approx \beta \cdot \frac{X_{t}^{\alpha}}{\gamma^{\alpha}} \rightarrow 0
  • X_{t} \rightarrow \inftyのとき:

    S\text{-Curve}(X_{t}) = \beta \cdot \frac{X_{t}^{\alpha}}{X_{t}^{\alpha} + \gamma^{\alpha}} \approx \beta \cdot \frac{X_{t}^{\alpha}}{X_{t}^{\alpha}} \rightarrow \beta

したがって、S-Curve関数は[0,\infty)上で0から\betaに収束する単調増加関数である。

【定理】S-Curve関数の変曲点

S-Curve関数S\text{-Curve}(X_{t})の変曲点はX_{t} = \gammaである。

以下に証明を与える。
S-Curve関数の2階微分を計算すると:

S\text{-Curve}''(X_{t}) = \beta \cdot \alpha \cdot \gamma^{\alpha} \cdot \frac{X_{t}^{\alpha-1}(\gamma^{\alpha} - X_{t}^{\alpha})}{(X_{t}^{\alpha} + \gamma^{\alpha})^3}

これが0になるのは、

X_{t}^{\alpha} = \gamma^{\alpha} \hspace{10pt} \Rightarrow \hspace{10pt} X_{t} = \gamma

したがって、X_{t} = \gammaでS-Curveの形状が凸から凹へ変わる。

大まかな多変量回帰モデルの構築

前章で数理的表現をしたものを用いて、MMMの多変量回帰モデルを構築していきます。

【定義】MMMの多変量回帰モデル

\begin{align*} y_t &= Intercept + \beta_j \times S\textrm{-Curve}(x_{decay, t, j}) \\ &+ \beta_{hol} \cdot hol_{t} + \beta_{sea} \cdot sea_{t} + \beta_{trend} \cdot trend_{t} + \cdots + \beta_{ETC} \cdot ETC_{t} + \varepsilon \end{align*}
項目 役割
y_{t} 目的変数(売上orコンバージョン数)
\textrm{Intercept} ベースラインの売上
\beta_{j} \times S\textrm{-Curve}(x_{\textrm{decay}, t, j}) 広告の影響(Adstock & S字カーブ)
\beta_{\textrm{hol}} \cdot \textrm{hol}_{t} 休日効果(祝日による売上変動)
\beta_{\textrm{sea}} \cdot \textrm{sea}_{t} 季節性効果(年間の変動)
\beta_{\textrm{trend}} \cdot \textrm{trend}_{t} 長期トレンド(市場の成長・衰退)
\beta_{\textrm{ETC}} \cdot \textrm{ETC}_{t} 競合やプロモーションの影響
\varepsilon 誤差項(モデルに含まれない要因)
  • アドストック関数変換(Adstock transformation)
    X_{\textrm{decay}, t, j} = X_{t,j} + \theta_{j} \cdot X_{\textrm{decay}, t-1, j}
  • Sカーブ関数変換(S-Curve transformation)
    S\text{-Curve}(x,j) = \beta_{j} \times \frac{x_{\text{decay},t,j}^{\alpha}}{x_{\text{decay},t,j}^{\alpha} + \gamma^{\alpha}}

【定義】誤差項の仮定

  • 誤差項\varepsilon_{t}は、売上に影響を与えるがモデルに含まれていない要因を表す。

  • 一般的に独立同分布(i.i.d.)の正規分布に従うと仮定する

    \varepsilon_{t} \sim \mathcal{N}(0,\sigma^{2})

初期モデルの構築(リッジ回帰による正則化)

まずは初期モデルを作成するのですが、MMMは多くの場合、多重共線性の問題に直面します。これは、複数の広告を同時に運用することが一般的だからです。しかし、回帰モデルにおいて複数の独立変数が相関している場合、多重共線性の問題が発生します。これにより係数推定が不安定になり、係数の有意性の解釈が困難になりオーバーフィッティングが生じます。

そこでRobynではL2正則化(リッジ)回帰を行って、上記のような多重共線性の課題に対処しています。

なぜラッソ(L1正則化)回帰ではなくリッジ(L2正則化)回帰なのか?

ラッソ回帰を使うと特徴量が小さい変数を一気に消してしまいます。
広告効果を検証する際に、このように一気に変数を消してしまうことは実務上使いにくく、リッジ回帰だと変数を消すのではなく、大きさを減らすだけにとどめることができるからです。

Robynで行われているリッジ回帰の損失関数は以下の通りです。

【定義】リッジ回帰の損失関数

L(\beta) = \sum_{t=1}^{n}(y_{t} - f_{\beta}(x_{t}))^{2} + \lambda \sum_{j=1}^{p}\beta_{j}^{2}

ここで、

  • L(\beta) は損失関数(目的関数)。
  • y_{t}は時刻 t における実際の売上(従属変数)
  • f_{\beta}(x_{t}) は回帰モデルの予測値
  • \beta_{j}j 番目の変数に対応する回帰係数
  • \lambda > 0 は正則化強度を制御するハイパーパラメータ

【定義】リッジ回帰の目的

リッジ回帰の目的は、以下のペナルティ付き最適化問題を解くことである。

\underset{\beta}{\text{argmin}} \left[ \sum_{t=1}^{T}(y_{t} - f_{\beta}(x_{t}))^{2} + \lambda \sum_{j=1}^{p}\beta_{j}^{2} \right]

また、上記を求めることによって、以下の目的を同時に達成します。

  1. S-Curve変換後の広告効果S(X_{\text{decay}, t, j})の係数\beta_{j}を推定
  2. 外生変数や統制変数Z_{t, k}, C_{t,m}の係数\beta_{k}, \beta_{m}を推定
  3. パラメータの過大適合(過学習)を防ぐために、係数の大きさを制約
なぜベイズモデルを使用しないのか?

ベイズモデルは、現場の実情に応じたモデリングが容易なので近年は人気が高まっています。
ただRobynでは、複数の目的を同時に達成することが求められているため、ベイズモデルではなくリッジ回帰(正則化回帰)を使います。

実際に、リッジ回帰正規分布を事前分布に仮定したベイズ回帰のMAP推定と等価であり、ラッソ回帰ラプラス分布を事前分布に仮定したベイズ回帰のMAP推定と等価になります。
これに関してはまた記事にしようと思います。

進化的アルゴリズムによる多目的変数の最適化

リッジ回帰で損失関数が最小になる\betaを求め、初期モデルを作成します。その後、以下の3つの目的関数を最小にすることを目標としてハイパーパラメータを最適化します。

【定義】正則化平均二乗誤差:NRMSE

データセットを\mathcal{D} = {(X_{t}, y_{t})}_{t=1}^{n}とし、モデルの予測値を\hat{y}_{t}とすると、NRMSE(Normalized Root Mean Squared Error)は以下のように定義される。

\text{NRMSE} = \frac{\sqrt{\frac{1}{n}\sum_{t=1}^{n}(y_{t} - \hat{y}_{t})^{2}}}{\text{max}(y_{t}) - \text{min}(y_{t})}

NRMSEは、通常平均二乗誤差(RMSE)をデータの範囲で正規化することで、異なるスケール間で比較可能な誤差指標となる。

【定義】分解二乗距離の平方根: DECOMP.RSSD

各メディアjの支出シェアをS_{j}、その影響の推定値をE_{j}とすると、分解二乗距離の平方根(DECOMP.RSSD: Decomposition Root Sum of Squared Distance)は以下のように定義される。

\text{DECOMP.RSSD} = \sqrt{\sum_{j=1}^{p}(S_{j} - E_{j})^{2}}

ここで、各成分は以下の通り:

  • S_{j} = \frac{\text{広告費}_{j}}{\sum_{k=1}^{p}\text{広告費}_{k}} (広告の支出シェア)
  • E_{j} = \frac{\text{推定影響}_{j}}{\sum_{k=1}^{p}\text{推定影響}_{k}} (広告の効果シェア)

DECOMP.RSSDは、広告支出の割合と広告割合の分解割合の差を測定する指標。広告費が高いメディアが過剰に効果を持つことを防ぐことで、現実的なモデルを確保することに繋げる。

【定義】平均絶対パーセント誤差: MAPE.LIFT

リフトテスト(キャンペーンの因果効果)を考慮したMAPE(Mean Absolute Percentage Error)は以下のように定義される。

\text{MAPE.LIFT} = \frac{1}{n}\sum_{t=1}^{n} \left| \frac{y_{t}^{\text{lift}} - \hat{y}_{t}}{y_{t}^{\text{lift}}} \right|

ここで、

  • y_{t}^{\text{lift}}はリフトテストによる実測された因果効果。
  • \hat{y}_{t}はモデルによる予測値。

MAPE.LIFTは、予測値がリフトテストの結果とどれだけ一致するかを測る指標で、モデルの予測値が実験的な因果推定と一致するように調整します。

多目的最適化のハイパーパラメータ

ここまでの内容を踏まえ、多目的最適化には以下の4つのハイパーパラメータを設定します。

  • Adstocking: Geometric広告効果を選択する場合の\theta、またはWeibull広告効果を選択する場合の\lambdak
  • Saturation: S-Curve関数における \alpha および \gamma
  • Regularization: リッジ回帰のペナルティ項に対する \lambda
  • Varidation: トレーニングデータのサイズ

【定義】進化的アルゴリズム

進化的アルゴリズムは、機械学習における最適化手法の一つで、自然選択の概念に基づく。パラメータの組み合わせを参考とし、複数の世代の反復を通じて良好な解を生成し、不適合な解を排除する。進化的アルゴリズムの手法なステップは次の通りである。

  1. 初期化(Initialization): \lambda, \alpha, \gammaなどのパラメータをランダムに初期化。
  2. 適合度評価(Fitness Evalutian): NRMSEやDECOMP.RSSDを評価。
  3. 選択と交叉(Selection & Crossover): 適合度の高い解を組み合わせ、新しいハイパーパラメタセットを生成。
  4. 突然変異(Mutation): 随機に小さな変更を加え、多様性を保つ。
  5. 反復(Iteration): 上記のプロセスを複数回繰り返し、最適化された解に近づける。

数学的に表現すると以下のようになる。

\theta_{t+1} \sim \text{Mutate}(\text{Crossover}(\text{Select}(\theta_{t})))

ここで、\theta_{t}t期目のパラメータの組み合わせを意味する。

【定義】反復回数

Robynでは、最適化された解を得るために2000回以上の反復を定めている。
反復ごとに次の世代のパラメータ\thetaに対してフィードバックを行い、\alpha, \gamma, \thetaの最適値を近似する。

【定義】反復の収束条件

Robynは、次の条件を満たす場合、モデルが収束したとみなす。

  1. 最後の四分位範囲の標準偏差が、最初の3つの四分位範囲の標準偏差の平均より小さいこと。

    \sigma_{Q4}^{\text{last}} < \frac{1}{3}\sum_{i=1}^{3}\sigma_{Q4}^{(i)}
  2. 最後の四分位範囲の絶対中央値が、最初の四分位範囲の絶対中央値-最初の3つの四分位範囲の標準偏差の平均*2よりも小さいこと。

    | M_{Q4}^{\text{last}} | < | M_{Q4}^{(1)} | - 2 \cdot\frac{1}{3}\sum_{i=1}^{3}| \sigma_{Q4}^{(i)} |

【定義】パレート最適解

目的関数の集合f_{1}, f_{2}, \ldots, f_{k}に対して、決定変数の集合\Theta上でパレート最適解は以下を満たす。

\forall \theta' \in \Theta, \hspace{10pt} f_{i}(\theta) \leq f_{i}(\theta') \hspace{10pt} \text{かつ} \hspace{10pt} \exists j \hspace{5pt} \text{s.t.} \hspace{5pt} f_{j}(\theta) < f_{j}(\theta')

これは、いずれの目的関数も改善できず、一方の関数を改善しようとすると他の関数が悪化する点である。

【定義】パレートフロント

全てのパレート最適解の目的関数値の集合をパレートフロントと呼ぶ。

\mathcal{P}^{\ast} = \{ f(\theta^{\ast}) \mid \theta^{\ast} \in \Theta^{\ast} \}

ここで、\mathcal{P}^{\ast}はパレートフロント(今回だとNRMSE・DECOMP.RSSD・MAPE.LIFTなどの目的関数の値の集合値)、\Theta^{\ast}はパレート最適解の集合(ハイパーパラメータの組み合わせ)を指す。

【定理】パレートフロントの収束

進化的アルゴリズムの進化反復により、複数回数kを増やすと、最適解の集合はパレートフロントに収束する。

\lim_{k \rightarrow \infty} \mathcal{P}_{k} = \mathcal{P}^{\ast}

最終モデルの選択

【定義】最終モデル選択

収束後に得られる複数のパレート最適解 \theta^{\ast} \in \Theta^{\ast} から、実務上の直感や現実的な制約に基づいて最終モデル \theta^{\text{final}} を選択する。

\theta^{\text{final}} = \text{arg}\underset{\theta \in \Theta^{\ast}}{\text{min}} \hspace{3pt} \text{BusinessPenalty}(\theta)

ここで、\text{BusinessPenalty}(\theta) は選択されたモデルがビジネスの現実にどれだけ適合するかを定量化する関数である。

予算配分の最適化

Robynの予算割り当てには2つのシナリオが存在します。具体的には以下の2つです。

  1. 設定したゴール(ROASもしくはCPA)を最大化する広告予算配分
  2. 設定したゴール(ROASもしくはCPA)を到達するために必要な最小予算配分

これらをまずは大雑把に最適化を行います(グローバル最適化)
そしてその後、細かい最適化を行います(ローカル最適化)

グローバル最適化には拡張ラグラジアン(AUGLAG)、ローカル最適化には逐次二次計画法(SLSQP) という手法によって求められます。以下に最適化したいシナリオを定義として書いておきます。

【定義】最大レスポンス最適化

x_{j} を任意のメディア j の投入広告費用とする。総予算 B とメディアレベルの制約 x_{j} \leq u_{j} が与えられた場合、ROASまたはCPAのレスポンスを最大化する最適なクロスメディア予算割り当ては次のように定義される。

\underset{x}{\max}\hspace{3pt}R(x) \hspace{5pt} \text{subject to} \hspace{5pt} \sum_{j=1}^{p}x_{j} \leq B, \hspace{10pt} x_{j} \leq u_{j} \hspace{3pt} \forall j

ここで、

  • R(x) 総レスポンス(ROASまたはCPA)
  • u_{j} は各メディアチャネルの上限予算

【定義】ターゲット効率最適化

x_{j} を任意のメディア j の投入広告費用とする。ターゲットROASまたはCPAが与えられた場合、レスポンス( R )を最大化する最適なクロスメディア予算割り当ては次のように定義される。

\underset{x}{\min} \sum_{j=1}^{p}x_{j} \hspace{5pt} \text{subject to} \hspace{5pt} R(x) \geq T, \hspace{5pt} x_{j} \leq u_{j} \hspace{3pt} \forall j

ここで、

  • T は達成すべきターゲットROASまたはCPA
  • u_{j} は各メディアチャネルの上限予算

さいごに

思いのほか記事が長くなってしまいましたが、Robynでやっていることを数学的に表現できたかなと思います。本記事の冒頭にも書きましたが、もし間違えている点などございましたら教えていただけると嬉しいです。ここからはいくつかの観点で個人的な感想を交えて、本記事を終わりたいと思います。

Robynで面白いと感じた点

❶パレートフロントを用いたモデルの選択プロセス

数理モデル的に最適なものを選ぶことはもちろん可能ではあるが、実務を行う人からすると直感に反する部分があったりどうしても制約条件を反映させたりすることが困難であることが課題でした。
しかしRobynではビジネスの知見を織り交ぜることができるよう、一度ヒトがモデルを見て実態に合った最良のモデルを選択するというプロセスが含まれているという点が面白いなと感じました。

❷独自の評価指標「DECOMP.RSSD」

モデルから導かれた投下コストシェアと、実データでのコストシェアの二乗和の平方根を表す指標ですが、これはモデルと実際の値が乖離しすぎる状況を防ぐ役割を担っています。
よくモデルの精度(NRMSE)だけを用いて評価することがあるのですが、DECOMP.RSSDという考え方を用いることで、統計的な妥当性と実務上の妥当性のバランスを取れる 非常に良くできた指標だと思いました。

Discussion