はじめに

こんにちは、事業会社で働いているデータサイエンティストです：

今回の記事では、ディリクレ過程を含めた、新しいマーケティングミックスモデルを提案します。この新しいモデルをアメリカのマクロ経済データに応用し、各変数が失業率にどのように寄与するかを可視化します。

しかし、内生性を考慮しない結果、経済理論と常識の予想とは逆の値になっている係数があり、連邦準備制度理事会（FRB、アメリカの中央銀行に相当）がこのようなモデルで意思決定すると、とんでもないことになることを示します。

最後に、内生性の正体を少し教科書的な説明とは違う視点から説明し、データサイエンティストのあなたが内生性の問題をどこまで真剣に考える必要があるのかについて話します。

（先に言っておくと、内生性の正体は理不尽がないことだと思います）

ではまずはモデルの話に入ります！

モデルの考え方

このモデルは、独立変数の過去効果の加重平均と、従属変数自身のトレンドが、従属変数の値を決めると仮定します。

この加重平均を取る考え方はJin, Wang, Sun, Chan, Koehler（2017）のデータサイエンス界隈で有名な論文を参考にしますが、geometric decayの推定がうまくいかないので、ランダムウォークを標準化して重みを作ります。また、Hill関数の推定もうまくいかないので、Hill関数は設定していません。

これはマーケティングとマクロ経済のドメインの違いから来るものかと思われます。

トレンドに関しては、ディリクレ過程で推定します。こちらの記事と近い考えです：

要するに、トレンドが変動したりすると言っても、近い範囲内だったら線形トレンドでいいよねという発想です。その適切な範囲などをディリクレ過程に推定してもらいます。ベイズ版の局所線形回帰のようなものです。また、トレンドの数まで勝手に推定してくれます。

モデル定式化

モデルはまず数式で書くとこのようになります：

全体のディリクレ過程

$$d\_alpha\sim Gamma(0.1, 0.1)$$

$$z\sim Stick\space Breaking(d\_alpha)$$

トレンドk

$$トレンド中心点_{k}\sim Normal(0, 1)$$

$$トレンド管轄範囲_{k}\sim Gamma(0.1, 0.1)$$

$$線形トレンド係数_{k}\sim Normal(0,1)$$

説明変数mのラグ数lの累計効果ウエイト

$$ウエイト_{m,l}\sim Softmax(Normal(ウエイト_{m,l-1}, 1))$$

説明変数mの対象時間tにおける累計効果

$$累計効果_{m,t} = ウエイト_{m} * 説明変数の値_{m,(t - ラグ数 + 1):t}$$

説明変数pの係数

$$sigma_{beta_{p}}\sim Gamma(0.1, 0.1)$$

$$beta_{p} \sim Normal(0,sigma_{beta_{p}}^2)$$

他のパラメータ

$$sigma \sim Gamma(0.1, 0.1)$$

$$切片 \sim Normal(0, 1)$$

観測値iの分布

$$\eta_{i} \sim Categorical(z)$$

$$時間_{i} \sim Normal(トレンド中心点_{\eta_{i}}, トレンド管轄範囲_{\eta_{i}}^2)$$

$$従属変数_{i} \sim Normal(切片 + 時間_{i} * 線形トレンド係数_{\eta_{i}} + beta '* 累計効果_{時間_{i}}, sigma^2)$$

モデル実装

Stanでの実装コードはこちらです：

mmm_dirichlet.stan

data {
  int media_type;
  int time_type;
  int P;
  int<upper=time_type> max_lag;
  int<upper=time_type> train_period;
  
  matrix[media_type, time_type] media;
  vector[time_type] normalized_time;
  array[time_type] real outcome;
}
parameters {
  real<lower=0> d_alpha;                                       // ディリクレ過程の全体のパラメータ
  vector<lower=0, upper=1>[P - 1] breaks;  // ディリクレ過程のstick-breaking representationのためのパラメータ
  vector[P] trend_middle_point;
  vector<lower=0>[P] trend_spread;
  vector[P] trend;
  vector<lower=0>[media_type] adweight_sigma;
  array[media_type] vector[max_lag] adweight_unnormalized;
  
  real<lower=0> sigma;
  real intercept;
  vector<lower=0>[media_type] beta_sigma;
  vector[media_type] beta;
}
transformed parameters {
  simplex[P] z; //
  matrix[media_type, time_type] adstock;
  for (m in 1:media_type){
    for (t in 1:time_type){
      if (t < max_lag){
        adstock[m, t] = softmax(adweight_unnormalized[m,1:t]) '* media[m, 1:t]';
      }
      else {
        adstock[m, t] = softmax(adweight_unnormalized[m, :]) '* media[m, (t - max_lag + 1):t]';
      }
    }
  }
  {
    // stick-breaking representationの変換開始
    // https://discourse.mc-stan.org/t/better-way-of-modelling-stick-breaking-process/2530/2 を参考
    z[1] = breaks[1];
    real sum = z[1];
    for (p in 2:(P - 1)) {
      z[p] = (1 - sum) * breaks[p];
      sum += z[p];
    }
    z[P] = 1 - sum;
  }
}
model {
  d_alpha ~ gamma(0.1, 0.1);
  breaks ~ beta(1, d_alpha);
  
  trend_middle_point ~ normal(0, 1);
  trend_spread ~ gamma(0.1, 0.1);
  trend ~ normal(0, 1);
  
  adweight_sigma ~ gamma(0.1, 0.1);
  
  
  for (m in 1:media_type){
    adweight_unnormalized[m, 1] ~ normal(0, 1);
    for (l in 2:max_lag)
    adweight_unnormalized[m, l] ~ normal(adweight_unnormalized[m, l - 1], adweight_sigma[m] ^ 2);
  }
  
  sigma ~ gamma(0.1, 0.1);
  intercept ~ normal(0, 10);
  beta_sigma ~ gamma(0.1, 0.1);
  
  beta ~ normal(0, beta_sigma .^ 2);
  
  for (t in 1:train_period){
    vector[P] case_vector;
    for (p in 1:P){
      case_vector[p] = log(z[p]) + 
                       normal_lpdf(normalized_time[t] | trend_middle_point[p], trend_spread[p] ^ 2) + 
                       normal_lpdf(outcome[t] | intercept + normalized_time[t] * trend[p] + beta '* adstock[:, t], sigma ^ 2);
    }
    target += log_sum_exp(case_vector);
  }
}
generated quantities {
  array[time_type] vector[P] estimated_state;
  
  for (t in 1:time_type){
    estimated_state[t] = rep_vector(0.0, P);
  }
  
  vector[time_type] predicted_trend;
  vector[time_type] predicted_outcome;
  matrix[media_type, time_type] ad_contribution;
  
  for (t in 1:time_type){
    vector[P] case_vector;
    for (p in 1:P){
      case_vector[p] = log(z[p]) + normal_lpdf(normalized_time[t] | trend_middle_point[p], trend_spread[p] ^ 2);
    }
    int state = categorical_rng(softmax(case_vector));
    estimated_state[t, state] = 1.0; 
    predicted_trend[t] = normalized_time[t] * trend[state];
    for (m in 1:media_type){
      ad_contribution[m, t] = beta[m] * adstock[m, t];
    }
    predicted_outcome[t] = normal_rng(intercept + predicted_trend[t] + beta '* adstock[:, t], sigma ^ 2);
  }
}

前処理

今回の記事では、アメリカの失業率を、物価（CPIAUCSL）、金利（FEDFUNDS）、投資（INVEST）で予測するモデルを作ります。

この記事はマクロ経済学の記事ではないため、変数定義の妥当性等は割愛します。興味ある方はこちらのサイトで変数名を検索してください：

まず、データを標準化し、差分を取ります：

data("fred_md", package = "BVAR")

fred_df <- fred_md |>
  tibble::tibble() |>
  dplyr::select(UNRATE, CPIAUCSL, FEDFUNDS, INVEST) |>
  dplyr::mutate(dplyr::across(c(UNRATE, CPIAUCSL, FEDFUNDS, INVEST), 
                              ~ (. - mean(.))/sd(.), .names = "{.col}_normalized")) |>
  
  dplyr::mutate(
    CPIAUCSL_normalized = CPIAUCSL_normalized - dplyr::lag(CPIAUCSL_normalized, 1),
    FEDFUNDS_normalized = FEDFUNDS_normalized - dplyr::lag(FEDFUNDS_normalized, 1),
    INVEST_normalized = INVEST_normalized - dplyr::lag(INVEST_normalized, 1),
  ) |>
  # 一期目がNAになるので外す
  tidyr::drop_na()

続いて、Stan用のデータを作成します：

fred_data_list <- list(
  media_type = 3,
  time_type = nrow(fred_df),
  P = 10,
  max_lag = 24,
  train_period = nrow(fred_df) - 4,
  media = fred_df |>
    dplyr::select(CPIAUCSL_normalized, FEDFUNDS_normalized, INVEST_normalized) |>
    as.matrix() |>
    t(),
  normalized_time = (1:nrow(fred_df) - mean(1:nrow(fred_df)))/sd(1:nrow(fred_df)),
  outcome = fred_df$UNRATE_normalized
)

モデル推定

では、モデル推定に入りましょう！

> m_mmm_init <- cmdstanr::cmdstan_model("mmm_dirichlet.stan")
>
> m_mmm_estimate <- m_mmm_init$variational(
     seed = 12345,
     data = fred_data_list,
     iter = 10000
 )
------------------------------------------------------------ 
EXPERIMENTAL ALGORITHM: 
  This procedure has not been thoroughly tested and may be unstable 
  or buggy. The interface is subject to change. 
------------------------------------------------------------ 
Gradient evaluation took 0.003549 seconds 
1000 transitions using 10 leapfrog steps per transition would take 35.49 seconds. 
Adjust your expectations accordingly! 
Begin eta adaptation. 
Iteration:   1 / 250 [  0%]  (Adaptation) 
Iteration:  50 / 250 [ 20%]  (Adaptation) 
Iteration: 100 / 250 [ 40%]  (Adaptation) 
Iteration: 150 / 250 [ 60%]  (Adaptation) 
Iteration: 200 / 250 [ 80%]  (Adaptation) 
Success! Found best value [eta = 1] earlier than expected. 
Begin stochastic gradient ascent. 
  iter             ELBO   delta_ELBO_mean   delta_ELBO_med   notes  
   100        -2286.309             1.000            1.000 
   200        -1970.363             0.580            1.000 
   300        -1900.413             0.399            0.160 
   400        -1805.990             0.312            0.160 
   500        -1800.496             0.250            0.052 
   600        -1806.877             0.209            0.052 
   700        -1789.804             0.181            0.037 
   800        -1798.083             0.159            0.037 
   900        -1787.790             0.142            0.010   MEDIAN ELBO CONVERGED 
Drawing a sample of size 1000 from the approximate posterior...  
COMPLETED. 
Finished in  7.6 seconds.
> 
> m_mmm_summary <- m_mmm_estimate$summary()

問題なく推定が終わりました！ガウス過程より早いですね！前のガウス過程の記事で永遠に待たされた記憶があります、、、

予測性能の可視化

m_mmm_summary |> 
  dplyr::filter(stringr::str_detect(variable, "predicted_outcome")) |>
  dplyr::mutate(
    id = purrr::map_int(
      variable,
      \(x){
        stringr::str_split(x, "\\[|\\]|,")[[1]][2] |>
          as.numeric()
      }
    )
  ) |>
  dplyr::left_join(
    tibble::tibble(
      date = as.Date(rownames(fred_md))
    ) |>
      dplyr::mutate(
        id = dplyr::row_number()
      ),
    by = "id"
  ) |>
  dplyr::bind_cols(
    answer = fred_df$UNRATE_normalized
  ) |>
  dplyr::mutate(
    status = ifelse(dplyr::row_number() <= nrow(fred_df) - 4, "学習", "検証")
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = date, y = answer, color = status)) + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = mean), color = "blue") + 
  ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::scale_x_date(date_breaks = "2 year", date_labels = "%Y") + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3") + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(title = "アメリカ失業率予測")

いい感じに学習できていますね！しかも曲線はかなり滑らかで、過学習していないようです。また。信用区間が広くなっているものの、コロナ（2020）の急な失業率の上昇を当てはめることをある程度諦めて、逆に検証データにおける精度の方が良さそうという嬉しい逆説的状況も見られます。

ではデータと予測結果の散布図も確認しましょう

m_mmm_summary |> 
  dplyr::filter(stringr::str_detect(variable, "predicted_outcome")) |>
  dplyr::bind_cols(
    answer = fred_df$UNRATE_normalized
  ) |>
  dplyr::mutate(
    status = ifelse(dplyr::row_number() <= nrow(fred_df) - 4, "学習", "検証")
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_point(ggplot2::aes(x = mean, y = answer, color = status)) + 
  ggplot2::geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3")

綺麗に45度線付近で推移していますね。特に検証データの点はほぼ45度線にあるので、予測性能はかなり良いです。繰り返しになりますが、2020年の状況を学習するのを正しく諦めたのはすごいと思いました。

トレンドの可視化

次は推定されたトレンドを可視化します：

m_mmm_summary |> 
  dplyr::filter(stringr::str_detect(variable, "predicted_trend")) |>
  dplyr::mutate(
    id = purrr::map_int(
      variable,
      \(x){
        stringr::str_split(x, "\\[|\\]|,")[[1]][2] |>
          as.numeric()
      }
    )
  ) |>
  dplyr::left_join(
    tibble::tibble(
      date = as.Date(rownames(fred_md))
    ) |>
      dplyr::mutate(
        id = dplyr::row_number()
      ),
    by = "id"
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = mean), color = "blue") + 
  ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95), fill = ggplot2::alpha("blue", 0.3)) + 
  ggplot2::scale_x_date(date_breaks = "2 year", date_labels = "%Y") + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3") + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(title = "アメリカ失業率トレンド")

何この信用区間笑ってなりますが、これはディリクレ過程をベイズ風局所線形回帰で利用した結果です。線形なのにどうして滑らかな曲線になっているの？などの質問もあると思いますが、これは管轄範囲にかぶりのある複数の線形トレンド同士の協働の結果（加重平均）として解釈してもいいと思います。

ただの感想になりますが、全体的に滑らかでトレンドっぽくて（？）いいですね。

さて、このトレンドの解釈なんですが、1973年あたりから急激な上昇が見られます。これは第1次オイルショックです。第2次オイルショックあたりでピークが来て、それ以降は2009年の世界金融危機まで線形トレンドになっています。コロナ禍前後の2019年から2022年までにも一つのピークが見られますが、2023年は下降傾向のトレンドが推定されています。

続いては、トレンドの遷移過程（変化点）を可視化します：

# ディリクレ過程に明示的にトレンドの数を教える必要がない代わりに、
# 不要なトレンドもデータとして保持してしまうため、
# モデルに利用されるトレンドしか表示しないようにする

estimated_state <- m_mmm_summary |> 
  dplyr::filter(stringr::str_detect(variable, "estimated_state")) |> 
  dplyr::pull(mean) |> 
  matrix(ncol = 10) 

used_states <- estimated_state |>
  nrow() |>
  seq_len() |>
  purrr::map_int(
    \(x){
      which.max(estimated_state[x,])
    }
  ) |>
  unique()

m_mmm_summary |> 
  dplyr::filter(stringr::str_detect(variable, "estimated_state")) |>
  dplyr::mutate(
    id = purrr::map(
      variable,
      \(x){
        stringr::str_split(x, "\\[|\\]|,")[[1]][2:3] |>
          as.numeric()
      }
    )
  ) |>
  tidyr::unnest_wider(id, names_sep = "_") |>
  dplyr::left_join(
    tibble::tibble(
      date = as.Date(rownames(fred_md))
    ) |>
      dplyr::mutate(
        id_1 = dplyr::row_number()
      ),
    by = "id_1"
  ) |>
  dplyr::filter(id_2 %in% used_states) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = mean, color = as.factor(id_2))) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95, fill = as.factor(id_2)),
                       alpha = 0.3) + 
  ggplot2::scale_x_date(date_breaks = "2 year", date_labels = "%Y") + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3") + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(title = "アメリカ失業率トレンド変化点検出") + 
  ggplot2::guides(color = "none", fill = "none")

この図の見方を説明すると、線と信用区間がx軸の時間がそのトレンドの所管範囲内にある確率を示しています。なので、色が変わった（所管範囲の信用区間が被った）ところを変化点として解釈できます。オイルショックから世界金融危機までは一つのトレンドが独占していることが確認できます。

全体的にいいものだからFRBがこのモデルで今後意思決定すればいいじゃん！いいえ、だめです。

破滅的な意思決定

ここではまず、MMMでよく登場する、説明変数の寄与度の図を確認しましょう：

m_mmm_summary |> 
  dplyr::filter(stringr::str_detect(variable, "contribution")) |>
  dplyr::mutate(
    id = purrr::map(
      variable,
      \(x){
        stringr::str_split(x, "\\[|\\]|,")[[1]][2:3] |>
          as.numeric()
      }
    )
  ) |>
  tidyr::unnest_wider(id, names_sep = "_") |>
  dplyr::mutate(
    id_1 = dplyr::recode(as.factor(id_1),
                         "1" = "CPIAUCSL",
                         "2" = "FEDFUNDS",
                         "3" = "INVEST")
  ) |>
  dplyr::left_join(
    tibble::tibble(
      date = as.Date(rownames(fred_md))
    ) |>
      dplyr::mutate(
        id_2 = dplyr::row_number()
      ),
    by = "id_2"
  ) |>
  ggplot2::ggplot() + 
  ggplot2::geom_line(ggplot2::aes(x = date, y = mean, color = as.factor(id_1))) + 
  ggplot2::geom_ribbon(ggplot2::aes(x = date, ymin = q5, ymax = q95, fill = as.factor(id_1)), alpha = 0.3) + 
  ggplot2::scale_x_date(date_breaks = "2 year", date_labels = "%Y") + 
  ggplot2::theme_gray(base_family = "HiraKakuPro-W3") + 
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1)) + 
  ggplot2::labs(title = "estimated contribution", color = "variable", fill = "variable")

金利の影響が一番大きいです。ではここで係数を確認しますが、（beta[2]が金利の係数です）

> m_mmm_summary |> 
     dplyr::filter(stringr::str_detect(variable, "beta"))
# A tibble: 6 × 7
  variable         mean  median    sd   mad      q5   q95
  <chr>           <dbl>   <dbl> <dbl> <dbl>   <dbl> <dbl>
1 beta_sigma[1]   1.25    1.13  0.541 0.455   0.600  2.27
2 beta_sigma[2]   4.48    4.32  1.33  1.35    2.65   6.88
3 beta_sigma[3]   1.10    0.993 0.517 0.430   0.475  2.06
4 beta[1]         0.238   0.227 0.594 0.614  -0.731  1.23
5 beta[2]       -10.4   -10.4   0.600 0.578 -11.4   -9.37
6 beta[3]         0.106   0.117 0.760 0.813  -1.14   1.31

うん？-10.4？金利を上げたら失業率が下がる？そんなわけないでしょ笑。

一般的なマクロ経済学の予想では、金利が下がったら、失業率も下がるはずです。金利が下がったら、企業がより積極的に投資でき、その分より多くの人を雇うはずです。企業に金銭的な余裕ができたのにリストラに踏み切るのは経済学の直感に反します。

モデルが推定した係数が経済学の直感に反する現象の理由は、内生性にあると思われます。要するに、利上げと利下げはランダムに決められるものではなく、FRB・中央銀行がいろいろ予想して決めるからです。

「来月の失業率が上がりそう！じゃあとりあえず今月はもう利下げしよ！」という意思決定を実施できるから、我々観測者が「失業率が上がるときは金利が下がりそう、マイナス相関がありそう」という結論に至ってしまいます。因果推論的にいうと、利下げしなかった世界は観測できません。FRBはそんな意思決定をしないからです。

ただ、反実仮想的に考えて、金利を下げなかったら失業率がもっとひどくなっていたはずなので、金利の失業率に対する真の因果関係は正です。

だから、FRBがこのモデルを信じて、今後大不況で失業者が溢れる状態で「よっしゃ！マイナス相関だから金利を上げて失業率を下げるぞ！」のような意思決定をしたら、とんでもないことになります。

このモデルはFRBの意思決定に使えません。

内生性の正体

この内生性はつまりどういうこと？因果推論以外は信用するなってこと？そういうわけではないです。

私は、理不尽がないことが内生性の正体だと理解しています。

日本の中央銀行になりますが、日本銀行法第三条を確認しましょう：

日本銀行の通貨及び金融の調節における自主性は、尊重されなければならない

自主性を尊重しないといけませんということです。

これは見方を変えると、要するに理不尽なことが比較的少ないということです。日銀が金利を上げようと思えば上げられるし、下げようと思えば下げられる。そんな自由な世界です。もちろん、円滑な政策実施のため、実際は日本政府ないし外国政府・中銀とコミュニケーションする必要がありますが、それでも理不尽のない世界です。

ABテストと比較しましょう。

くじを引きました。統制群です。薬が飲めません、新しいフロント画面が見られません。理不尽しかないです。統制群くじを引いただけなのにどうしてこんな目に遭わないといけないの、、、？

あなたの手元の因果推論・効果検証の本を確認してみてください。因果推論の手法（傾向スコアマッチング、DID、回帰不連続デザインなど）はどれもこのような説明変数の変な理不尽の存在を前提としています。

なぜ説明変数のの理不尽がないと因果推論がないというと、理不尽がランダム性のそのものだからです。

MMMは使えないの？

話をMMMに戻しますが、会社のマーケティング部門は、中央銀行ほどの自主性を持っていません。

元々TVCMのリリースが来週だったのに、商品企画部門から急にTVCMの内容を変更してくださいと要望され、リリースをやむをえず3ヶ月後に変更した。

来月からOOHを出稿する予定だったのに、枠がないと代理店に言われた

どれも理不尽です。でも仕方ないです。

このように、実際の企業内のデータは、そんなに内生性まみれではなく、説明変数にある程度の外生的な変動が担保されると思います。なので、MMMが全く使えないわけではなく、まずは関連部門に説明変数の値がどのように生まれたのかを確認して、そこから判断しましょう。

結論

いかがでしたか？

記事内の内生性に関する考察は、どれも大学院在学中ではなく、データサイエンティストになった後いろいろ見たり考えたりして得た結論です。

大学院生の皆さんもぜひアカデミアの他に、民間企業への就職を考えてみてください。

アカデミアでは手に入れない知識が得られます！

ご応募お待ちしております！

最後に、この記事のレポジトリのリンクを貼ります：

（多分M1 Mac周りの問題でbuildがうまくいけてないのでまだdevelopブランチにマージしたくないです）

参考文献

Jin, Yuxue, Yueqing Wang, Yunting Sun, David Chan and Jim Koehler. "Bayesian methods for media mix modeling with carryover and shape effects." Google Research (2017).