はじめに

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

突然ですが、予測精度が一番高いモデルを選べば、一番良いモデルを選んだことになるのでしょうか？

機械学習モデルを作るとき、RMSEやMAE、AUCなどの予測精度を比較して、「一番当たるモデル」を採用することはごく自然な発想です。実際、目的が純粋な予測であれば、それで何の問題もないケースも多いと思います。

しかし、モデルを使って「この施策を実施したらどうなるのか」を推定したい場合や、さらにその結果を使って「では、どの施策を選ぶべきなのか」という意思決定まで行いたい場合、本当にそれだけでいいのでしょうか？

本記事では、この疑問についてブラックボックスで、なんでもよしなにやってくれるとよく思われるLightGBMと、データ生成過程の構造に近い線形モデルを使ったシミュレーションで検証してみます。

ちなみに今回は、高度なダブル機械学習など、機械学習による推定のバイアスを補正する手法はあえて使いません。意思決定を、Kaggleでもおなじみの機械学習モデルに、特別なチューニングもせず「LightGBMでポン！」と任せたらどうなるのか、を検証します。

先に結果を少しだけお見せすると、なかなかカオスなことになりました。

予測精度で比較すると、1000回のシミュレーションすべてでLightGBMが勝ちました。データ生成過程の構造に近い線形モデルはボロ負けです。

ところが、各個体に対して「施策を実施した場合」と「実施しなかった場合」の予測値の差を取り、真の効果をどれだけ正確に推定できているかをRMSEで比較すると、今度は1000回すべてで線形モデルが勝ちました。

「なるほど。では効果推定をしたいなら線形モデルを選べばいいんですね！」、と思いますよね。

ところが、さらに「推定された効果の符号を使って、施策を実施すべきかどうかを正しく選べたか（効果がプラスだったら処置して、マイナスだったら処置しない）」という意思決定の正確さで比較すると、話はまた変わります。

LightGBMのほうが意思決定を間違える割合が高かったのは、1000回中 44.9% にすぎませんでした。つまり、効果そのもののRMSEでは100%負けていたLightGBMが、意思決定という別の基準では線形モデルとほぼ互角、むしろ勝つケースも普通に出てきます。

予測では100%勝つ。効果推定では100%負ける。でも意思決定ではほぼ互角。

本文では他の指標も紹介しますが、もうすでにわけわからん草

一体どれを信じればいいんだよ、という話です。

そして、もう一つ本記事で注目したいのが安定性です。

線形モデルは今回のデータ生成過程にかなり近い構造を明示的に与えています。一方、LightGBMにはそのような構造をほとんど教えず、データから関係を学習してもらいます。その結果、LightGBMは線形モデルより良い意思決定をすることもある一方で、シミュレーションごとの意思決定精度のばらつきも大きくなりました。

これは、前回のガウス過程の記事でも触れた、モデル構築を柔軟なモデルにアウトソースすることにも価格があるという話につながります。

LightGBMやガウス過程のような柔軟なモデルを使えば、人間がすべての関数形や交互作用を一つずつ考えてモデルに書き込む必要はありません。これは非常に強力です。しかし、誰かにお金を払っているわけではないからといって、そのアウトソースが無料なわけではありません。

その価格は、推定結果の不安定性や、モデルが何を学習したのかを人間が完全にはコントロールできないこととして支払うことになるかもしれません。

もちろん、だから「LightGBMを使うな」「線形モデルを使え」という話ではありません。今回のシミュレーションでは、実際にLightGBMが線形モデルより良い結果を出す評価基準もあります。むしろ本記事で言いたいことは、その逆です。

万能なモデルはありません。そして、万能な評価指標もありません。

予測したいのであれば、予測精度を見るべきです。効果そのものを正確に推定したいのであれば、効果推定の誤差を見るべきです。ただし、現実のデータでは、同じ個体について「処置を受けた場合」と「受けなかった場合」の結果を同時に観測することはできません。これは因果推論の根本問題であり、今回のシミュレーションのように真の効果を直接使って精度を評価することはできないため、実務ではA/Bテストなどを通じた検証が必要になります。そして、そのモデルを使って何らかの施策を選択したいのであれば、最終的な意思決定の質を見る必要があります。

つまり、「LightGBMと線形モデルのどちらが優れているか」という問いそのものが、場合によっては間違っています。

何を解きたいのかを決めない限り、何をもって「良いモデル」と呼ぶのかも決められないからです。

これはビジネスでもアカデミアでも同じだと思います。高度な手法を使うこと自体が価値なのではなく、解くべき問題を正しく定式化し、その問題に対応した評価基準を置き、そのうえで適切なモデルを選ぶことに価値があります。場合によっては、手法選択よりもその前の何を問題として定義するかのほうが、はるかに重要かもしれません。

ということで今回は、LightGBMと線形モデルを1000回戦わせながら、

「予測精度」「効果推定精度」「意思決定精度」は、本当に同じモデルを勝者として選ぶのか？

を検証していきます。

では、シミュレーション劇場の開幕です！

シミュレーションの設定

では、今回のシミュレーションの設定を説明します。

データ生成過程は、Gelman et al. Bayesian Workflow の 7.2節 "Causal inference as generalization" で紹介されている設定を参考にしています。

こちらのページに非商用でしたらダウンロードできるPDFファイルのリンクがありますので、興味ある方は是非ご確認ください：

今回登場する変数は非常にシンプルで、

x：0から10までの連続変数
z：0または1を取る処置変数
y：結果変数

の3つです。

まず、x は一様分布から生成します。

$$
x_i \sim Uniform(0,10)
$$

次に、処置 z は完全にランダムに割り当てるのではなく、x の値によって処置を受ける確率が変化するようにします。

$$
p_i = \frac{\exp(x_i-5)}{1+\exp(x_i-5)}
$$

$$
z_i \sim Bernoulli(p_i)
$$

つまり、x が小さい個体ほど z = 0 になりやすく、x が大きい個体ほど z = 1 になりやすい設定です。例えば x = 5 では処置確率が50%ですが、そこから離れるほど処置確率は0または1に近づいていきます。

ここが今回のシミュレーションでは重要です。

単純なA/Bテストのように全領域で z = 0 と z = 1 が同程度観測されるわけではありません。特に x が0や10に近い領域では、片方の処置しかほとんど観測されなくなります。そのためモデルは、観測された結果変数を単純に予測するだけでなく、データが相対的に少ない処置との組み合わせについても結果を推定する必要があります。

結果変数 y は、次の線形モデルから生成します。

$$
y_i = \beta_0 + \beta_1 x_i + \beta_2 z_i + \beta_3 x_i z_i + \epsilon_i
$$

$$
\epsilon_i \sim Normal(0,0.5^2)
$$

係数は、

$$
(\beta_0,\beta_1,\beta_2,\beta_3) = (0.1,0.2,-3,0.4)
$$

とします。

ここでは x と z の交互作用を入れているため、z の効果はすべての個体で一定ではありません。z = 1 にした場合と z = 0 にした場合の期待値の差を計算すると、

$$
\tau(x) = E[y \mid x,z=1] - E[y \mid x,z=0] =(-3+0.4x)
$$

となります。

したがって、真の効果は x によって変化します。

例えば x = 0 では真の効果は -3、x = 5 では -1、x = 10 では 1 です。そして、

$$
-3+0.4x=0
$$

を解くと、

$$
x=7.5
$$

となるため、x = 7.5 が意思決定の境界になります。

今回のシミュレーションでは単に「効果をどれだけ正確に当てられたか」だけでなく、推定された効果の符号を使って、

$$
\hat{\tau}(x)>0
$$

なら z = 1 を選び、

$$
\hat{\tau}(x)<0
$$

なら z = 0 を選ぶ、という意思決定も評価します。

つまり、真の効果の大きさを多少間違えていても符号が正しければ意思決定には成功します。一方、真の効果にかなり近い値を推定できていても、ちょうど0をまたいでしまえば意思決定としては失敗します。

この違いが、後ほどかなり面白くて、少し怖い結果を生みます😱

実際のRコードは以下です。

set.seed(123)

b <- c(0.1, 0.2, -3, 0.4)

train_df <- 50000 |>
  runif(0, 10) |>
  tibble::tibble(
    x = _
  ) |>
  dplyr::mutate(
    z = x |>
      purrr::map_int(
        \(x){
          rbinom(1, 1, exp(x-5)/(1 + exp(x-5)))
        }
      ),
    e = rnorm(dplyr::n(), sd = 0.5),
    y = b[1] + b[2] * x + b[3] * z + b[4] * x * z + e
  )

test_df <- 50000 |>
  runif(0, 10) |>
  tibble::tibble(
    x = _
  ) |>
  dplyr::mutate(
    z = x |>
      purrr::map_int(
        \(x){
          rbinom(1, 1, exp(x-5)/(1 + exp(x-5)))
        }
      ),
    e = rnorm(dplyr::n(), sd = 0.5),
    y = b[1] + b[2] * x + b[3] * z + b[4] * x * z + e,
    true_effect = (b[1] + b[2] * x + b[3] * 1 + b[4] * x * 1) - (b[1] + b[2] * x + b[3] * 0 + b[4] * x * 0)
  )

学習データとテストデータはそれぞれ50,000件生成します。シミュレーションなので、通常の実データでは観測できない個体ごとの真の効果 true_effect もこちらでは知ることができます。

そして、この同じデータを使って、アウトカムの真の関数形に近い構造を明示的に与えた線形モデルと、そうした関数形を明示的には与えずデータから柔軟に関係を学習する LightGBM を比較します。

ここで重要なのは、「どちらが最強のモデルなのか」を決めることではありません。

今回比較したいのは、

観測された y を正確に予測できるか
真の効果 $\tau(x)$ を正確に推定できるか
z = 0 と z = 1 のどちらを選ぶべきかを正しく判断できるか

という、似ているようで実は異なる3つのタスクです。

同じデータ、同じ2つのモデルでも、「何を正解とするか」を変えたとき、モデルの評価はどう変わるのでしょうか。

ここから、だんだん話がおかしくなっていきます笑

まずは1回、LightGBMと線形モデルを戦わせてみる

ではいきなり1000回のシミュレーションを回す前に、まずは一つのシードでLightGBMと線形モデルを推定し、何が起きるのかを確認してみましょう。

今回比較するのは、データ生成過程のアウトカム構造と同じ関数形を持つ線形ベイズモデルと、関数形を人間が明示的に指定せずにデータから柔軟に学習するLightGBMです。

線形モデル

まず線形モデルから見ていきます。

前節で説明した通り、今回の結果変数は、

$$
y_i = \beta_0 + \beta_1 x_i + \beta_2 z_i + \beta_3 x_i z_i + \epsilon_i
$$

というデータ生成過程から作っています。

そこでベイズモデル側でも、

$$
E[Y\mid X=x,Z=z] = \beta_0+\beta_1x+\beta_2z+\beta_3xz
$$

という同じ関数形を持つモデルを作ります。事前分布はあえて特に設定せず、均一分布とします。

Stanコードは非常にシンプルです。

linear_model.stan

data {
  int N;

  vector[N] x;
  vector[N] z;
  vector[N] y;
}
parameters {
  vector[4] beta;
  real<lower=0> sigma;
}
model {
  y ~ normal(beta[1] + beta[2] * x + beta[3] * z + beta[4] * x .* z, sigma);
}

ここで一点だけ注意が必要です。

このモデルはアウトカムのデータ生成過程と同じ関数形を持っているものの、データ生成過程全体をモデル化しているわけではありません。

前節では、

$$
P(Z=1\mid X=x) = \frac{\exp(x-5)} {1+\exp(x-5)}
$$

という処置割当の仕組みも使っていましたが、今回のベイズモデルにはこの部分は入っていません。あくまで $Y\mid X,Z$ の条件付き期待値をモデル化しています。

したがって、本記事では「完全に正しく指定されたモデル」とは呼ばず、アウトカムの真の関数形を知っている、かなり有利な線形モデルくらいに考えてください。

モデルを変分推論で推定して、係数を保存します。

m_init_lm <- cmdstanr::cmdstan_model("linear_model.stan")

m_estimate_lm <- m_init_lm$variational(
  seed = 123,
  data = list(
    N = nrow(train_df),
    x = train_df$x,
    z = train_df$z,
    y = train_df$y
  )
)

m_summary_lm_beta <- m_estimate_lm$summary() |>
  dplyr::filter(stringr::str_detect(variable, "beta")) |>
  dplyr::pull(mean)

対戦相手はLightGBM

一方、LightGBMにはこのような関数形を教えません。

x と z を渡して、

「よろしく！」

とします。

m_lightgbm <- lightgbm::lightgbm(
  data = lightgbm::lgb.Dataset(
    data = train_df |>
      dplyr::select(x, z) |>
      as.matrix(),
    label = train_df$y,
    categorical_feature = c("z")
  ),
  params = list(
    objective = "regression",
    metric = "rmse"
  )
)

ハイパーパラメータの大規模な探索なども行いません。

今回検証したいのは「LightGBMを限界までチューニングしたら線形モデルに勝てるか」ではなく、アウトカムの構造を人間がかなり具体的に与えたモデルと、構造の学習を柔軟な機械学習モデルにそのまま任せた場合に、評価するタスクによって何が起きるのかだからです。

z=0とz=1の世界をそれぞれ予測する

次に、テストデータの各個体について z = 0 と z = 1 の2パターンを作り、それぞれのモデルから予測値を出します。

LightGBMについては以下のように計算します。

z_change_lightgbm_df <- c(0, 1) |>
  purrr::map(
    \(this_z){
      test_df |>
        dplyr::select(x, z) |>
        dplyr::mutate(z = this_z) |>
        as.matrix() |>
        predict(object = m_lightgbm) |>
        tibble::tibble(z = _) |>
        `colnames<-`(stringr::str_c("z_", this_z))
    }
  ) |>
  dplyr::bind_cols() |>
  dplyr::bind_cols(
    y_pred = test_df |>
      dplyr::select(x, z) |>
      as.matrix() |>
      predict(object = m_lightgbm),
    y = test_df$y
  )

線形モデルでも同じことを行います。

z_change_lm_df <- c(0, 1) |>
  purrr::map(
    \(this_z){
      test_df |>
        dplyr::select(x, z) |>
        dplyr::mutate(z = this_z) |>
        dplyr::mutate(
          y = m_summary_lm_beta[1] + m_summary_lm_beta[2] * x + m_summary_lm_beta[3] * z + m_summary_lm_beta[4] * x * z
        ) |>
        dplyr::select(y) |>
        `colnames<-`(stringr::str_c("z_", this_z))
    }
  ) |>
  dplyr::bind_cols() |>
  dplyr::bind_cols(
    y_pred = test_df |>
      dplyr::select(x, z) |>
      dplyr::mutate(
        y = m_summary_lm_beta[1] + m_summary_lm_beta[2] * x + m_summary_lm_beta[3] * z + m_summary_lm_beta[4] * x * z
      ) |>
      dplyr::pull(y),
    y = test_df$y
  )

ここで、

$$
\hat{\tau}(x) = \hat{Y}(z=1)-\hat{Y}(z=0)
$$

を計算すれば、モデルが推定した個体ごとの効果を得ることができます。

そして今回はシミュレーションなので、真の効果

$$
\tau(x)=-3+0.4x
$$

も分かっています。

これを利用して、今回は次の4つの評価指標を見ます。

1. 予測RMSE

まずは普通の予測問題です。

$$
RMSE_{pred} = \sqrt{\frac{1}{N}\sum_{i=1}^{N}(y_i-\hat{y_i)}^2}
$$

実際に観測された z のもとで、結果変数 y をどれだけ正確に予測できたかを評価します。

2. 効果推定RMSE

次に、

$$
RMSE_{effect} = \sqrt{\frac{1}{N}\sum_{i=1}^N(\tau_i-\hat\tau_i)^2}
$$

を計算します。

こちらは結果変数そのものではなく、z = 0 から z = 1 に変えたときの効果をどれだけ正確に推定できたかを見る指標です。

3. 効果推定の平均的なズレ

さらに、

$$
\frac{1}{N}\sum_i\tau_i - \frac{1}{N}\sum_i\hat\tau_i
$$

も計算します。

コード上では bias_estimation としています。これは個体ごとの推定誤差ではなく、真の平均効果と推定された平均効果のズレを見るためのものです。

4. 意思決定の間違い率

そして最後が今回かなり重要な指標です。

真の効果と推定された効果の符号が逆になった割合、

$$
\frac{1}{N}
\sum_i
I\left(
\tau_i\hat\tau_i<0
\right)
$$

を計算します。

例えば真の効果が正なのにモデルが負と推定してしまえば、本来 z = 1 を選ぶべき個体に対して z = 0 を選んでしまいます。その逆も同様です。

つまりこれは、効果の値そのものをどれだけ正確に当てたかではなく、その推定値を使って意思決定したときに、どれくらい間違えるのかを見る指標です。

実際の計算は以下です。

z_change_lightgbm_df |>
  dplyr::bind_cols(
    answer = test_df$true_effect
  ) |>
  dplyr::mutate(
    diff = z_1 - z_0
  ) |>
  dplyr::summarise(
    effect_rmse = sqrt(mean((answer - diff)^2)),
    pred_rmse = sqrt(mean((y - y_pred)^2)),
    bias_estimation = mean(answer) - mean(diff),
    wrong_choice = mean(diff * answer < 0)
  )


z_change_lm_df |>
  dplyr::bind_cols(
    answer = test_df$true_effect
  ) |>
  dplyr::mutate(
    diff = z_1 - z_0
  ) |>
  dplyr::summarise(
    effect_rmse = sqrt(mean((answer - diff)^2)),
    pred_rmse = sqrt(mean((y - y_pred)^2)),
    bias_estimation = mean(answer) - mean(diff),
    wrong_choice = mean(diff * answer < 0)
  )

結果はこちらです。

# A tibble: 1 × 4
  effect_rmse pred_rmse bias_estimation wrong_choice
        <dbl>     <dbl>           <dbl>        <dbl>
1       0.103     0.503         0.00522      0.00804
# A tibble: 1 × 4
  effect_rmse pred_rmse bias_estimation wrong_choice
        <dbl>     <dbl>           <dbl>        <dbl>
1      0.0578     0.512         -0.0561       0.0169

上がLightGBM、下が線形モデルです。

🤔

もうこの時点で、ちょっと様子がおかしいですね。

結果変数の予測RMSEではLightGBMのほうが小さいです。

ところが、効果推定RMSEでは線形モデルのほうが小さくなっています。

ここまでは「なるほど、予測に強いLightGBMと、効果推定に強い線形モデルということか」と思うかもしれません。

しかし、意思決定の間違い率 wrong_choice を見ると、また話がひっくり返ります。

このシードでは、

LightGBM：0.804%
線形モデル：1.69%

となり、効果そのものは線形モデルのほうが正確に推定できているのに、どちらの処置を選ぶべきかという意思決定ではLightGBMのほうが間違いが少なくなりました。

なぜこんなことが起きるのでしょうか。

理由の一つは、効果推定RMSEと意思決定の間違い率が、そもそも違うものを正解としているからです。

効果推定RMSEでは、$\hat\tau(x)$ が $\tau(x)$ からどれくらい離れているかが重要です。

一方、今回の意思決定では極端に言えば、$\mathrm{sign}(\hat\tau(x))$さえ正しければ、効果の大きさを多少間違えていても意思決定としては正解です。

逆に、真の効果にかなり近い値を推定できていても、意思決定境界である0をほんの少しまたいでしまえば、その個体については不正解になります。

つまり、この1回の結果だけでも、

「予測がうまい」「効果を正確に推定できる」「良い意思決定ができる」は、同じ意味ではない

ことが見えてきました。

とはいえ、たまたまシード123でこうなっただけかもしれません。

ということで、ここから1000回シードを変えて、本当にこのカオスが再現するのかを確認していきましょう。

大規模シミュレーション

先ほどはシード123の1回だけを取り出して、LightGBMと線形モデルの結果を比較しました。

しかし、1回のシミュレーションだけでは、たまたまそのシードでそうなっただけかもしれません。特に今回注目したいのは、単純な平均的な性能だけではなく、データを生成し直したときに、それぞれのモデルの性能がどの程度安定しているのかという点です。

そこでここからは、シードを1から1000まで変えながら、先ほどとまったく同じ実験を1000回繰り返します。

各シードでは、まず50,000件の学習データと50,000件のテストデータを新しく生成します。データ生成過程そのものはすべてのシードで共通で、

$$
X_i \sim Uniform(0,10)
$$

$$
Z_i \sim Bernoulli \left(\frac{\exp(X_i-5)}{1+\exp(X_i-5)}\right)
$$

$$
Y_i = 0.1+0.2X_i-3Z_i+0.4X_iZ_i+\epsilon_i
$$

$$
\epsilon_i\sim N(0,0.5^2)
$$

としています。

つまり、真の世界そのものは変えず、そこから観測されるサンプルだけを1000回作り直すことになります。

そして、それぞれの学習データに対して、前節と同じ線形モデルとLightGBMをそれぞれ推定します。

線形モデルには、

$$
E[Y\mid X,Z] = \beta_0+\beta_1X+\beta_2Z+\beta_3XZ
$$

というアウトカムの真の関数形と同じ構造を与えています。一方、LightGBMには x と z だけを渡し、交互作用の形などを人間から明示的には教えずに学習させます。

そのうえで、50,000件のテストデータそれぞれについて、実際に観測された z における予測値に加えて、強制的に z = 0 とした場合と z = 1 とした場合の予測値をそれぞれ計算します。

したがって、各個体について、

$$
\hat{\tau}_i = \hat{Y}_i(1)-\hat{Y}_i(0)
$$

というモデル上の効果を計算できます。

シミュレーションなので、真の効果

$$
\tau_i=-3+0.4X_i
$$

も分かっています。

各シードについて記録するのは、前節と同じく次の4つの指標です。

pred_rmse：実際に観測された結果変数 y の予測RMSE
effect_rmse：真の効果 $\tau_i$ と推定された効果 $\hat{\tau}_i$ のRMSE
bias_estimation：真の平均効果と推定された平均効果の差
wrong_choice：真の効果と推定された効果の符号が逆になり、z = 0 と z = 1 の選択を間違えた割合

これをLightGBMと線形モデルの両方について保存するため、最終的な simulation_df は1000行 × 8指標のデータになります。

また、今回は1000回のシミュレーションを順番に実行すると時間がかかるため、future と furrr を使って20並列で計算しています。

コードはこちらです：

m_init_lm <- cmdstanr::cmdstan_model("linear_model.stan")
future::plan(future::multisession(workers = 20))
simulation_df <- 1000 |>
  seq_len() |>
  furrr::future_map(
    \(this_seed){
      set.seed(this_seed)

      b <- c(0.1, 0.2, -3, 0.4)

      train_df <- 50000 |>
        runif(0, 10) |>
        tibble::tibble(
          x = _
        ) |>
        dplyr::mutate(
          z = x |>
            purrr::map_int(
              \(x){
                rbinom(1, 1, exp(x-5)/(1 + exp(x-5)))
              }
            ),
          e = rnorm(dplyr::n(), sd = 0.5),
          y = b[1] + b[2] * x + b[3] * z + b[4] * x * z + e
        )

      test_df <- 50000 |>
        runif(0, 10) |>
        tibble::tibble(
          x = _
        ) |>
        dplyr::mutate(
          z = x |>
            purrr::map_int(
              \(x){
                rbinom(1, 1, exp(x-5)/(1 + exp(x-5)))
              }
            ),
          e = rnorm(dplyr::n(), sd = 0.5),
          y = b[1] + b[2] * x + b[3] * z + b[4] * x * z + e,
          true_effect = (b[1] + b[2] * x + b[3] * 1 + b[4] * x * 1) - (b[1] + b[2] * x + b[3] * 0 + b[4] * x * 0)
        )

      m_estimate_lm <- m_init_lm$variational(
        seed = 123,
        data = list(
          N = nrow(train_df),
          x = train_df$x,
          z = train_df$z,
          y = train_df$y
        )
      )

      m_summary_lm_beta <- m_estimate_lm$summary() |>
        dplyr::filter(stringr::str_detect(variable, "beta")) |>
        dplyr::pull(mean)

      m_lightgbm <- lightgbm::lightgbm(
        data = lightgbm::lgb.Dataset(
          data = train_df |>
            dplyr::select(x, z) |>
            as.matrix(),
          label = train_df$y,
          categorical_feature = c("z")
        ),
        params = list(
          objective = "regression",
          metric = "rmse"
        )
      )

      z_change_lightgbm_df <- c(0, 1) |>
        purrr::map(
          \(this_z){
            test_df |>
              dplyr::select(x, z) |>
              dplyr::mutate(z = this_z) |>
              as.matrix() |>
              predict(object = m_lightgbm) |>
              tibble::tibble(z = _) |>
              `colnames<-`(stringr::str_c("z_", this_z))
          }
        ) |>
        dplyr::bind_cols() |>
        dplyr::bind_cols(
          y_pred = test_df |>
            dplyr::select(x, z) |>
            as.matrix() |>
            predict(object = m_lightgbm),
          y = test_df$y
        )

      z_change_lm_df <- c(0, 1) |>
        purrr::map(
          \(this_z){
            test_df |>
              dplyr::select(x, z) |>
              dplyr::mutate(z = this_z) |>
              dplyr::mutate(
                y = m_summary_lm_beta[1] + m_summary_lm_beta[2] * x + m_summary_lm_beta[3] * z + m_summary_lm_beta[4] * x * z
              ) |>
              dplyr::select(y) |>
              `colnames<-`(stringr::str_c("z_", this_z))
          }
        ) |>
        dplyr::bind_cols() |>
        dplyr::bind_cols(
          y_pred = test_df |>
            dplyr::select(x, z) |>
            dplyr::mutate(
              y = m_summary_lm_beta[1] + m_summary_lm_beta[2] * x + m_summary_lm_beta[3] * z + m_summary_lm_beta[4] * x * z
            ) |>
            dplyr::pull(y),
          y = test_df$y
        )


      z_change_lightgbm_df |>
        dplyr::bind_cols(
          answer = test_df$true_effect
        ) |>
        dplyr::mutate(
          diff = z_1 - z_0
        ) |>
        dplyr::summarise(
          lightgbm_effect_rmse = sqrt(mean((answer - diff)^2)),
          lightgbm_pred_rmse = sqrt(mean((y - y_pred)^2)),
          lightgbm_bias_estimation = mean(answer) - mean(diff),
          lightgbm_wrong_choice = mean(diff * answer < 0)
        ) |>
        dplyr::bind_cols(
          z_change_lm_df |>
            dplyr::bind_cols(
              answer = test_df$true_effect
            ) |>
            dplyr::mutate(
              diff = z_1 - z_0
            ) |>
            dplyr::summarise(
              lm_effect_rmse = sqrt(mean((answer - diff)^2)),
              lm_pred_rmse = sqrt(mean((y - y_pred)^2)),
              lm_bias_estimation = mean(answer) - mean(diff),
              lm_wrong_choice = mean(diff * answer < 0)
            )
        )
    },
    .progress = TRUE,
    .options = furrr::furrr_options(seed = TRUE)
  ) |>
  dplyr::bind_rows()

このシミュレーションで見たいのは、単純に

「LightGBMと線形モデル、どっちが強いの？」

というモデル同士の殴り合いではありません。

むしろ確認したいのは、

同じ2つのモデルを、予測・効果推定・意思決定という異なる基準から1000回評価したとき、それでも一貫して同じモデルを「良いモデル」と呼べるのか？

という点です。

そしてもう一つ、平均的にどちらが優れているかだけでなく、データを生成し直すたびに性能がどれくらい変動するのかにも注目します。

1回だけなら偶然かもしれません。

では、1000回やったらどうなるのでしょうか。

1000回回したら、評価指標ごとに「勝者」が変わった

では結果を見てみましょう。

simulation_df |>
  dplyr::mutate(rid = dplyr::row_number()) |>
  tidyr::pivot_longer(!rid) |>
  dplyr::mutate(
    name = name |>
      purrr::map(
        \(x){
          splitted <- stringr::str_split(x, "_")[[1]]
          return(
            c(splitted[1], stringr::str_c(splitted[2], "_", splitted[3]))
          )
        }
      )
  ) |>
  tidyr::unnest_wider(name, names_sep = "_") |>
  dplyr::rename(model = name_1, metric = name_2) |>
  ggplot2::ggplot() + 
  ggplot2::geom_density(ggplot2::aes(x = value, fill = model), alpha = 0.5) + 
  ggplot2::facet_wrap(~ metric, scales = "free")

かなり面白いことになりました。

まず、この図は1000回のシミュレーションそれぞれについて計算した4つの評価指標の分布を、LightGBMと線形モデルに分けて描いたものです。

右に行くほど値が大きくなりますが、今回の4指標はいずれも基本的にゼロに近いほど良い指標です。

そして一目見て分かる通り、どちらか一方のモデルがすべての指標で勝っているわけではありません。

せっかく1000回シミュレーションしたので、分布を見るだけでなく、各シードにおいてLightGBMと線形モデルのどちらが勝ったのかも数えてみましょう。

simulation_df |>
  dplyr::summarise(
    lightgbm_wins_effect_rmse = mean(lightgbm_effect_rmse < lm_effect_rmse), 
    lightgbm_wins_pred_rmse = mean(lightgbm_pred_rmse < lm_pred_rmse), 
    lightgbm_wins_bias_estimation = mean(abs(lightgbm_bias_estimation - 0) < abs(lm_bias_estimation - 0)), 
    lightgbm_wins_wrong_choice = mean(lightgbm_wrong_choice < lm_wrong_choice)
  ) |>
  dplyr::glimpse()

結果はこちらです。

Rows: 1
Columns: 4
$ lightgbm_wins_effect_rmse     <dbl> 0
$ lightgbm_wins_pred_rmse       <dbl> 1
$ lightgbm_wins_bias_estimation <dbl> 0.998
$ lightgbm_wins_wrong_choice    <dbl> 0.551

つまり、LightGBMの勝率は、

評価指標	LightGBMの勝率	線形モデルの勝率
予測RMSE	100.0%	0.0%
効果推定RMSE	0.0%	100.0%
平均効果の絶対バイアス	99.8%	0.2%
意思決定の間違い率	55.1%	44.9%

となりました。

100%勝ったと思ったら、評価指標を変えただけで0%になる。さらに別の指標では99.8%、意思決定では55.1%になる。

同じ1000個のデータ生成過程、同じLightGBM、同じ線形モデルを比較しているにもかかわらずです。

この数字だけでも、今回の記事で言いたいことのかなりの部分が見えてきます。

予測精度ならLightGBMが100%勝つ

まず左下の pred_rmse を見ると、これは非常に分かりやすいです。

LightGBMの分布はおおむね0.50付近、線形モデルは0.51付近にあり、LightGBMのほうが一貫して予測RMSEが小さくなっています。

実際、1000回のシミュレーションすべてにおいて、

$$
RMSE_{pred}^{LightGBM} < RMSE_{pred}^{LM}
$$

となりました。

つまり、

「未知のデータに対して $Y$ を正確に予測したい」

という問題として考えるのであれば、今回のシミュレーションではLightGBMを選ぶのが自然です。

これは少し面白い結果でもあります。

今回の線形モデルには、

$$
Y=\beta_0+\beta_1X+\beta_2Z+\beta_3XZ+\epsilon
$$

という、アウトカムの真のデータ生成過程と同じ関数形をかなり露骨に教えています。

それでも、純粋な予測精度では1000回すべてLightGBMが勝っています。

「真の構造に近いモデルを作れば、当然あらゆる意味で機械学習モデルより強い」というほど話は単純ではありません。

ところが効果推定では線形モデルが100%勝つ

次に右上の effect_rmse を見ると、順位が完全に逆転します。

今度は線形モデルの分布が明確に左側にあり、LightGBMよりも効果推定RMSEが小さくなっています。

そして1000回のシミュレーションを個別に比較してみても、LightGBMの勝率は、0%です。

つまり1000回すべてにおいて、線形モデルの効果推定RMSEのほうが小さくなりました。

$\hat{\tau}(x)=\hat{Y}(1)-\hat{Y}(0)$ を使って、処置を変えたときの効果そのものを正確に知りたいという問題にすると、予測精度では100%勝っていたLightGBMが、今度は100%負けます。

これは今回の設定を考えれば、それほど不思議ではありません。

線形モデルには最初から $X\times Z$ という正しいアウトカム構造を与えています。一方、LightGBMはあくまで観測された $Y$ の予測を目的として学習しているだけであり、個体ごとの $\hat{Y}(1)-\hat{Y}(0)$ を正確に復元することを直接目的として最適化されているわけではありません。

予測がうまいモデルだからといって、効果の差分まで同じようにうまく推定できるとは限らない。

今回の結果はその違いをかなり極端な形で示しています。

予測精度なら1000勝0敗。

効果推定精度なら0勝1000敗。

「どちらのモデルが優れているのか？」とだけ聞かれても、もはや答えようがありません。

しかし平均効果のバイアスを見ると、また話がひっくり返る

そして左上の bias_estimation が個人的にはかなり面白いところです。

効果推定RMSEでは線形モデルが1000回すべて勝っていました。

ところが、真の平均効果と推定された平均効果の差を見ると、LightGBMはおおむね0を中心に分布しているのに対し、線形モデルは負の方向にかなり明確なズレを持っています。

そこで符号そのものではなく、

$$
|\mathrm{bias}|
$$

を比較して「どちらが0に近いか」を1000回数えると、LightGBMの勝率は 99.8% になりました。

つまり今回のシミュレーションでは、

個体ごとの効果推定RMSEでは線形モデルが100%勝っているのに、平均効果の絶対バイアスではLightGBMが99.8%勝っている

という結果になりました。

ここまで来ると、「効果推定」という一言でまとめることすら危険です。

個体ごとの効果を正確に復元したいのか。

平均的な効果を正確に知りたいのか。

この違いだけでも、モデルの評価は大きく変わります。

機械学習モデルをそのまま平均処置効果などの効果推定に利用するとバイアスが問題になる、という議論は因果推論ではよく出てきます。そのためダブル機械学習など、機械学習を利用しながらバイアスを制御するためのさまざまな方法が提案されています。

しかし、今回の結果を見て、

「だから単純な回帰モデルを使えば安全！」

と結論づけることもできません。

少なくとも今回のシミュレーションでは、単純なベイズ線形モデルの平均効果の推定値のほうが系統的にズレており、平均効果の絶対バイアスという評価基準ではLightGBMが1000回中998回勝っています。

もちろん、ここは大げさに解釈すべきではありません。

今回の設定ではその差の大きさが実質的にどれほど重要なのかは別問題ですし、このシミュレーションから「LightGBMは平均処置効果推定にも安全だ」などという一般論を導くつもりもありません。

むしろ私が怖いと思うのは逆です。

シミュレーションだからこそ、我々はこのバイアスがどちら向きに、どれくらい発生しているのかを知ることができます。

実データでは真のデータ生成過程を知りません。

真の $\tau(x)$ も知りません。

したがって、同じことが実データで起きたとき、そのバイアスが0.01なのか、0.05なのか、もっと大きいのか、正方向なのか負方向なのかを、この図のように答え合わせすることはできません。

そこがシミュレーションと実務の決定的な違いです。

そして意思決定になると、さらにカオスになる

右下の wrong_choice はさらに面白いです。

線形モデルは比較的狭い範囲に集中している一方で、LightGBMの分布はかなり広がっています。

つまりLightGBMは、非常に低い間違い率を達成するシードもある一方で、線形モデルよりかなり悪くなるシードもあります。

そして1000回について直接勝敗を数えると、LightGBMの勝率は、55.1% でした。

予測RMSEでは100%。

効果推定RMSEでは0%。

平均効果の絶対バイアスでは99.8%。

それが、結局どちらの処置を選ぶべきなのか？という意思決定になると55.1% です。

これは一見すると奇妙ですが、評価している対象が違うことを考えれば矛盾ではありません。

effect_rmse は、

$$
|\tau(x)-\hat{\tau}(x)|
$$

の大きさを評価しています。

一方、今回の wrong_choice で重要なのは、

$$
\mathrm{sign}(\tau(x)) = \mathrm{sign}(\hat{\tau}(x))
$$

かどうかです。

真の効果が $0.5$ のときに $1.0$ と推定しても、意思決定としては正解です。

反対に真の効果が $0.01$ のときに $-0.01$ と推定すれば、数値としてはほとんど外していないにもかかわらず、意思決定としては不正解になります。

だからこそ、効果推定RMSEで100%勝つことと、意思決定で100%勝つことはまったく同じではありません。

そしてLightGBMの wrong_choice の分布が広いことも無視できません。

これは前回のガウス過程の記事でも触れた、柔軟なモデルにモデリングを「外注」することの対価に近いものだと私は考えています。

人間が関数形を考える仕事を減らし、データから複雑な関係を学習してくれるモデルに任せる。

もちろん、それ自体には大きな価値があります。

しかし、その柔軟性は無料ではありません。

お金を誰かに払っているわけではなくても、推定結果の不安定性や、想定していなかった挙動という形で代金を払う可能性があります。

今回のLightGBMの wrong_choice の広い分布は、その一例として見ることができます。

「どのモデルが一番強いですか？」という質問自体が間違っているかもしれない

ここまでの結果をまとめると、かなりカオスです。

予測RMSEを見る → LightGBMが100%勝つ
効果推定RMSEを見る → LightGBMが0%勝つ
平均効果の絶対バイアスを見る → LightGBMが99.8%勝つ
意思決定の間違い率を見る → LightGBMが55.1%勝つ

同じデータ生成過程です。

同じ2つのモデルです。

変えたのは、「何をもって良いモデルとするのか」だけです。

だからこそ、この記事で言いたいのは「LightGBMはダメ」とか「線形モデル最強！」という話ではありません。

むしろ逆です。

手法を選ぶ前に、何を良くしたいのかを決めなければならない。

予測したいのであれば予測精度を見るべきです。

効果そのものを正確に推定したいのであれば、効果推定の誤差を見るべきです。

平均的な効果を知りたいのであれば、その平均量に対する誤差を見る必要があります。

そしてモデルを使って最終的に何らかの施策を選ぶのであれば、我々が本当に気にするべきなのは、その意思決定によって得られる結果です。

100%勝つモデルが、評価指標を変えた瞬間に0%勝つモデルになる。

この結果は、モデル選択以前に問題設定が重要であることをかなり強烈に示していると思います。

どれだけ高度なモデルを知っているかよりも、

「そもそも何を正解としてモデルを評価するべきなのか」を正しく定義できること。

実務でも研究でも、こちらのほうがはるかに重要な能力なのではないかと私は思います。

おわりに

いかがでしたか？

今回は、LightGBMと、アウトカムの真の関数形に近い構造を持つ線形ベイズモデルを1000回戦わせてみました。

結果をもう一度まとめると、

予測RMSE：LightGBMが1000回中1000回勝利
効果推定RMSE：LightGBMが1000回中0回勝利
平均効果の絶対バイアス：LightGBMが99.8%勝利
意思決定の間違い率：LightGBMが55.1%勝利

となりました。

もうめちゃくちゃです、意味わかんない。どのモデルがどの指標で勝ったかを期末テストで聞かれたら泣く笑

しかし、この記事で一番伝えたかったのは、まさにこの「めちゃくちゃさ」です。

「どちらのモデルが優れているのか？」という質問には、何をしたいのかを決めない限り答えられません。

$Y$を当てたいのであれば、予測RMSEを見るべきです。

個体ごとの効果 $\tau(x)$ を正確に推定したいのであれば、効果推定RMSEを見るべきです。

平均的な効果を知りたいのであれば、その平均量に対するバイアスを見る必要があります。

そして、そのモデルを使って最終的に z = 0 と z = 1 のどちらを選ぶのかを決めたいのであれば、その意思決定がどれだけ正しいのかを評価しなければなりません。

同じモデルでも、

評価するタスクが変われば、100%勝っていたモデルが0%勝つモデルになります。

これは「LightGBMはダメ」「線形モデルが正義」という話ではありません。

むしろ、そんな単純な結論を出せないことこそが今回の結論です。

また、今回LightGBMの意思決定精度は線形モデルよりばらつきが大きくなりました。

私は前回のガウス過程の記事でも、柔軟なモデルを使うことを、モデリング作業の一部をモデルにアウトソースすることだと表現しました。

人間が関数形や交互作用を一つずつ考えなくても、LightGBMのようなモデルはデータから非常に複雑な関係を学習してくれます。

これはものすごく便利です。

しかし、無料ではありません。

誰かに請求書を渡されるわけではなくても、推定結果のばらつき、想定していない相互作用、説明しづらい挙動、そして意思決定の不安定性という形で代金を払う可能性があります。

一方で、構造を人間が与えるモデルも無料ではありません。

今回はこちらが真のアウトカム関数を知っているという、線形モデルにかなり有利な世界を作っています。現実ではそんな答え合わせはできません。

人間が考えた構造そのものが間違っていれば、今度は間違った構造を自信満々にモデルへ教えてしまうリスクがあります。

だから私は、柔軟なモデルか構造モデルか、ベイズか機械学習か、といった「宗派」を決めること自体にはあまり興味がありません。

それよりも、

何を知りたいのか。

何を改善したいのか。

最終的にどの意思決定を良くしたいのか。

を先に考えることのほうが重要だと思います。

高度な手法を知っていることにはもちろん価値があります。

しかし、ビジネスでも研究でも、価値を生み出すのは手法そのものではなく、正しい問題を定義し、その問題に対応した評価基準を置き、その目的に適した手法を選ぶことです。

今回のシミュレーションを一言でまとめるなら、

万能な手法はない。万能な評価指標もない。

そして最後に、私が最近ますます大事だと思っていることを書いて終わります。

手法にコミットするのではなく、より良い意思決定にコミットする。

モデルが1000勝0敗でも、評価する問題を変えれば0勝1000敗になる世界では、それくらいの覚悟でちょうどいいのではないでしょうか。

最後に、私たちと一緒に働きたい方はぜひ下記のリンクもご確認ください：