EconMLで実践するビジネス因果推論（第9回）その結論、信じていいのか（レビュー会の5つの質問にプラセボ検定と感度分析で答える）

手法を決め、数字が出て、資料を作りました。レビュー会でデータに詳しい役員がこう聞きます。「その422円、信じていいの？」

この問いに「頑健です」と答えてはいけません。答えるべきは「どういう条件なら覆るか」です。今回は第2回の推定（クーポン効果422円）を5つの角度から揺さぶり、数字で答えます。

結論を先に3行で書きます。

モデルを替えても、期間を分けても、施策をシャッフルしても、422円は動きませんでした。動いたのは会員ランクを外したときだけ（924円）です。
感度分析の答えは「耐荷重 0.11」。見えていない交絡が会員ランク級なら覆り、流入経路や在籍月数の程度なら耐えます。
信じていい結論とは、ある5つの質問に数字で答えられる結論のことです。

この記事を読み終えると、次のことができるようになります。

レビュー会で聞かれる5つの質問を、先回りして確かめられる
プラセボ検定と期間分割で「偽の効果」「時期のブレ」を切り分けられる
sensitivity_summary と robustness_value の数字を日本語に翻訳できる
「効果がない」と「分からない」を言い分けられる

レビュー会で詰められる5つの質問

交絡は入れ切れているか。大事な変数を1つ忘れていたら？
モデルを変えても同じか。勾配ブースティングだから出た数字では？
効くはずのないもので効果が出ないか。手順そのものが効果を作り出していないか？
期間を分けても同じか。たまたまその時期だけでは？
見えていない要因で覆らないか。データにない交絡があったら？

①〜④は「同じ結論を別の条件で出し直す」だけなので、第2回の関数を1つ書けば全部回せます。⑤だけ専用の道具（感度分析）を使います。

今回は、これまでのデータに変数month（月）が追加されたデータセットを使います。

以下からダウンロードできます。

ec_coupon.csv

クーポンを受け取った人はサイトを見に来る回数が増え、たくさん買った人ほどサイトに来ています。つまり、施策の後に、施策と成果の両方の影響を受けて動く変数です。

まず、動かしてみる：①交絡を1つずつ外す

以下、コードです。

import numpy as np
import pandas as pd
from sklearn.ensemble import (
    GradientBoostingRegressor, 
    GradientBoostingClassifier
)
from econml.dml import LinearDML

# データの読み込み
df = pd.read_csv("ec_coupon.csv")

# 変数の定義
Y = df["spend_30d"]
T = df["coupon"]
X = df[["past_freq", "days_since_last"]]
W = pd.get_dummies(df[["member_rank", "channel"]], drop_first=True)
W["tenure_months"] = df["tenure_months"]

# 第2回と同じ LinearDML を学習
def fit_ate(Y, T, X, W, model_y=None, model_t=None):
    """第2回と同じ LinearDML を学習し、推定ATEと95%区間を返す"""
    if model_y is None:
        model_y = GradientBoostingRegressor(random_state=0)
    if model_t is None:
        model_t = GradientBoostingClassifier(random_state=0)
    est = LinearDML(
        model_y=model_y, 
        model_t=model_t, 
        discrete_treatment=True,
        cv=5, random_state=0
    )
    est.fit(Y, T, X=X, W=W)
    lo, hi = est.ate_interval(X, alpha=0.05)
    return est, est.ate(X), lo, hi

# 基準（第2回）
est, ate, lo, hi = fit_ate(Y, T, X, W)
print(f"基準（第2回） \t\t{ate:.0f}円 [{lo:.0f}, {hi:.0f}]")

# ① 交絡を1つずつ外して、答えがどれだけ動くか
for key in ["member_rank", "channel", "tenure_months"]:
    cols = [c for c in W.columns if c.startswith(key)]
    _, a, lo, hi = fit_ate(Y, T, X, W.drop(columns=cols))
    print(f"① {key:<14}を外す \t{a:.0f}円 [{lo:.0f}, {hi:.0f}]")

以下、実行結果です。

基準（第2回） 		422円 [350, 493]
① member_rank   を外す 	924円 [837, 1010]
① channel       を外す 	419円 [348, 490]
① tenure_months を外す 	444円 [371, 517]

会員ランクを外すと2倍に跳ね、流入経路と在籍月数は外しても±25円です。この結論を支えているのは会員ランク1本だと分かります。「入れ切れているか」への答えは、「ランク級の変数を忘れていれば覆る。忘れていないことは第4回の仕分けで確認した」です。

②モデル差し替え、③プラセボ、④期間分割

以下、コードです。

from sklearn.ensemble import (
    RandomForestRegressor, 
    RandomForestClassifier
)

# ② 中の機械学習モデルを差し替える
rf_y = RandomForestRegressor(
    n_estimators=300, 
    min_samples_leaf=10, 
    random_state=0
)
rf_t = RandomForestClassifier(
    n_estimators=300, 
    min_samples_leaf=10, 
    random_state=0
)
_, a, lo, hi = fit_ate(Y, T, X, W, rf_y, rf_t)
print(f"② ランダムフォレスト \t{a:.0f}円 [{lo:.0f}, {hi:.0f}]")

# ③ プラセボ：クーポンの列をシャッフルする（効果は0のはず）
T_fake = pd.Series(
    np.random.default_rng(0).permutation(T.to_numpy()), 
    index=T.index,
    name="coupon"
)
_, a, lo, hi = fit_ate(Y, T_fake, X, W)
print(f"③ シャッフル \t\t{a:.0f}円 [{lo:.0f}, {hi:.0f}]")

# ④ 配布月で前半・後半に分けて、それぞれ推定する
for label, m in [
    ("1〜3月", df["month"] <= 3), 
    ("4〜6月", df["month"] >= 4)
]:
    _, a, lo, hi = fit_ate(Y[m], T[m], X[m], W[m])
    print(
        f"④ {label} \t\t{a:.0f}円 "
        f"[{lo:.0f}, {hi:.0f}]  n={m.sum()}"
    )

以下、実行結果です。

② ランダムフォレスト 	430円 [357, 502]
③ シャッフル 		39円 [-16, 94]
④ 1〜3月 		467円 [373, 562]  n=2451
④ 4〜6月 		365円 [259, 471]  n=2549

②は430円。線形回帰に替えても423円で、モデルには依りません。③は誰に配ったかを無作為に入れ替えた「偽のクーポン」で、効果は39円、区間は0をまたぎます。5回シャッフルしても−30〜57円でした。手順が効果を作り出してはいないということです。④は467円と365円で区間が重なり、どちらも正解413円を含みます。

⑤見えていない要因で、どこまで覆るか：感度分析

①〜④はデータの中で確かめられます。⑤はデータにない変数の話なので、確かめようがない。代わりに「どのくらい強い見落としがあれば結論が覆るか」を計算します。

以下、コードです。

# ⑤ 感度分析：見えていない交絡がどのくらい強ければ、結論が覆るか
print(est.sensitivity_summary())
print(f"耐荷重（robustness_value）: {est.robustness_value():.3f}")
for c in [0.05, 0.1, 0.2]:
    lo, hi = est.sensitivity_interval(c_y=c, c_t=c)
    print(f"c={c}: [{lo:.0f}, {hi:.0f}]")

以下、実行結果です。

Sensitivity Analysis Summary for c_y=0.05, c_t=0.05, rho=1.0
================================================
CI Lower Theta Lower Theta  Theta Upper CI Upper
------------------------------------------------
 162.056     232.303 356.61     480.918   550.98
   Robustness Values for null_hypothesis=0
==============================================
Robustness Value (Theta) Robustness Value (CI)
----------------------------------------------
                   0.137                 0.111
----------------------------------------------
耐荷重（robustness_value）: 0.111
c=0.05: [162, 551]
c=0.1: [31, 682]
c=0.2: [-258, 970]

読み方は3点です。c は「見えていない交絡が、成果と施策の残りのばらつきをどれだけ説明するか」の割合。c=0.05 を仮定すると区間は [162, 551] に広がり、c=0.11 で下限が0を割る。この0.11 が耐荷重（robustness_value）です。

Theta が422円ではなく357円なのは、感度分析が効果を「X で変わらない1本の値」として扱うためで、第4回で見た X なしの推定値と同じです。区間の広がり方を見るのが目的なので、この差は気にしなくて構いません。

0.11 が大きいのか小さいのかは、手元の交絡と比べると分かります。同じ物差しで測ると、会員ランクは0.18、流入経路は0.002、在籍月数はほぼ0でした。ランク級の変数をもう1つ見落としていれば覆る。経路・在籍級なら耐える。これが⑤への答えです。

裏で何が起きているのか：橋の耐荷重

「未観測の交絡はありませんか」と聞かれて「ありません」と答えられる人はいません。あるかないかは分からない。分かるのは、あったとして、どのくらい強ければ結論が壊れるかです。

橋の耐荷重と同じです。「絶対に落ちません」とは誰も言えませんが、「10トンまで耐えます」とは言える。そして通る車が2トンなら安心し、12トンのトラックが来るなら通行止めにする。感度分析は、見えていない交絡がどこまでの重さなら耐えるかを数字にし、手元の交絡（ランク、経路、在籍）を「通る車の重さ」の目安として並べる道具です。

耐荷重は「安心」の数字ではなく「限界」の数字です。今回の 0.11 は、業務側の人に「ランクと同じくらい配布と購入に効く要因で、データにないものはありますか」と聞く問いに変わります。あるなら、その変数を取りに行くのが次の仕事です。

レビュー会の1枚

5つの質問と確かめ方、今回の答えを1枚にしました。

この1枚を推定値と一緒に持っていきます。答えが出せない質問があれば、その分析はまだ会議に出せません。

ありがちな失敗

1つ目：区間が0をまたいだことを「効果がない」と報告する。 同じデータから500人だけ取り出して推定すると、142円 [−79, 363] でした。区間が0をまたいでいます。でも正解は413円で、効果はあります。0をまたぐのは「効果がない」ではなく「このデータでは分からない」。報告の言葉を間違えると、効く施策が止まります。

2つ目：「頑健です」で終わる。 5つ確かめて全部通ったとしても、答えは「頑健です」ではありません。「ランク級の見落としがあれば覆ります。業務側で心当たりはありますか」です。覆る条件を言える人が、いちばん信用されます。