はじめに

こんにちは、事業会社で働く新卒データサイエンティストです。

今回は、MMM（マーケティングミックスモデリング、メディアミックスモデリング）におけるアドストック関数を考えます。

アドストック関数

一般に、当日の広告投資の効果が当日やそれにつづく数週間のうちに見られるとは限りません。広告投資の効果が現れるには一定の時間を待たねばならないこともあります。

そこでMMMでは、各メディアへの広告投資効果が減衰する様子をモデリングすることが多いです。

以下、ワイブル分布アドストックと幾何的アドストックの可視化を目指します。

準備

期間は全部で12期あり、3つのメディアが存在すると考える。変数は以下の通り。

$adstock_{m,t}\cdots$$t$期におけるメディア$m$のアドストック
$advertisement_{m,t}\cdots$ $t$期におけるメディア$m$の変数（インプレッション数や広告費など）
$maxlag \cdots$減衰の最大ラグ数。ここでは仮に6とする。

ワイブル分布アドストック

$adstock_{m,t} = \sum_{k=1}^{maxlag}adweight_{m,k} * advertisement_{m,t-k+1}
= \genfrac{(}{)}{0pt}{}{\sum_{k=1}^{maxlag}adweight_{m,k} * advertisement_{m,t-k+1} \quad for \quad (t\geq maxlag)}{\sum_{k=1}^{t}adweight_{m,k} * advertisement_{m,t-k+1} \quad for \quad (otherwise)} $

ただし、
$adweight_{m,k} = \exp\lbrace(\frac{-(t-k+1)}{\lambda})^l\rbrace$,
$lambda$は関数の形を調整するパラメータで$\lambda = \frac{1}{(-log(0.001))^{1/l}}$,
$l$は減衰パラメータで$l=1$を仮定。

これをRで実装します。

time_type <- 12
max_lag <- 6
media_type <- 3

set.seed(123)
advertisement <- matrix(runif(media_type * time_type, 0, 1), nrow = media_type, ncol = time_type)
l_values <- c(5,5,5)

# Function to compute decay weight
compute_weights <- function(l, t_values) {
  lambda <- 1 / (-log(0.001))^(1 / l)
  exp(- (t_values / lambda)^l)
}

# Prepare data frame for visualization
contributions <- tibble::tibble()
for (m in 1:media_type) {
  l <- l_values[m]  
  
  for (t in 1:time_type) {
    if (t < max_lag) {
      time_indices <- 1:t
    } else {
      time_indices <- (t - max_lag + 1):t
    }
    
    weights <- sum(compute_weights(l, rev(seq_along(time_indices))))  
    ads <- advertisement[m, time_indices]
    weighted_contributions <- weights * ads
    
    # Store data
    temp_df <- tibble::tibble(
      media_type = factor(m),
      time_type = t,
      lag = time_indices,
      contribution = weighted_contributions
    )
    
    contributions <- dplyr::bind_rows(contributions, temp_df)
  }
}

# Plot the weighted contributions over time
ggplot2::ggplot(contributions, ggplot2::aes(x = lag, y = contribution, fill = factor(time_type))) +
  ggplot2::geom_bar(stat = "identity", position = "stack") +
  ggplot2::facet_wrap(~ media_type, ncol = 1) +
  ggplot2::labs(
       x = "Lag Time",
       y = "Weighted Contribution",
       fill = "Time Step") +
  ggplot2::theme_minimal()

綺麗ではないですが、アドストックが過去の期のものを含みながら減衰していく様子が分かります。

一般に、減衰期間が長期の場合は減衰パラメータ$l$が大きくなるほど残存効果は大きい、つまり緩やかにアドストックが減衰するとされています。

幾何的アドストック

一般的には

$adstock_{m,t} =advertisement_{m,t} + \lambda_m\times adstock_{m,t-1}$
（ただし、$\lambda_m$はメディア$m$に関する減衰パラメータ）

なるアドストックが知られています。

ここでは、減衰パラメータをソフトマックス関数で表現するやり方をとります。このやり方のソースはチームメンバーの記事です。

$adstock_{m,t} = \sum_{k=1}^{maxlag}\lbrace softmax(adweight_{m,k})\times advertisement_{m, t-k+1}\rbrace = \genfrac{(}{)}{0pt}{}{\sum_{k=1}^{maxlag}\lbrace softmax(adweight_{m,k})\times advertisement_{m, t-k+1}\rbrace \quad for \quad (t\geq maxlag)}{\sum_{k=1}^{t}\lbrace softmax(adweight_{m,t})\times advertisement_{m, t-k+1}\rbrace \quad for \quad (otherwise)} $

ただし、
$adweight_{m,k} \sim N(0,1) \quad for \quad k =1$
$adweight_{m,k} \sim N(adweight_{m,k-1},\sigma_{adweight}) \quad for \quad k =\lbrace2,3,\cdots,maxlag\rbrace$
を仮定。

これをRで実装してみます。

time_type <- 12
max_lag <- 6
media_type <- 3

set.seed(123)
advertisement <- matrix(runif(media_type * time_type, 0, 1), nrow = media_type, ncol = time_type)
adweight_unnormalized <- matrix(runif(media_type * max_lag, 0, 1), nrow = media_type, ncol = max_lag)

softmax <- function(x) {
  exp(x) / sum(exp(x))
}


contributions <- tibble::tibble()
for (m in 1:media_type) {
  for (t in 1:time_type) {
    if (t < max_lag) {
      weights <- softmax(adweight_unnormalized[m, 1:t])
      ads <- advertisement[m, 1:t]
      time_indices <- 1:t
    } else {
      weights <- softmax(adweight_unnormalized[m, 1:max_lag])
      ads <- advertisement[m, (t - max_lag + 1):t]
      time_indices <- (t - max_lag + 1):t
    }
    
    weighted_contributions <- weights * ads
    
    temp_df <- data.frame(
      media_type = factor(m),
      time_type = t,
      lag = time_indices,
      contribution = weighted_contributions
    )
    
    contributions <- dplyr::bind_rows(contributions, temp_df)
  }
}

ggplot2::ggplot(contributions, ggplot2::aes(x = lag, y = contribution, fill = factor(time_type))) +
  ggplot2::geom_bar(stat = "identity", position = "stack") +
  ggplot2::facet_wrap(~ media_type, ncol = 1) +
  ggplot2::labs(
       x = "Lag Time",
       y = "Weighted Contribution",
       fill = "Time Step") +
  ggplot2::theme_minimal()

ワイブル分布アドストックのときよりも綺麗に減衰を可視化することができました。
今回は減衰パラメータをソフトマックス関数で作成しているため、事前にパラメータを設定する必要がありませんでした。

おわりに

MMMでは、自社サービスの形態に合わせて適切な形状のアドストック関数を用いることが重要だと思います。いくつかのレパートリーを持っておくことで、モデルを柔軟に構築できるのではないでしょうか。

参考文献

博報堂DYメディアパートナーズ. 2023. Marketing Mix Modeling Guidebook. https://www.hakuhodody-media.co.jp/aaas/news/mmmguidebook.html
（2024年1月27日最終取得）