0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

因果推論 Day 25/全30回 Causal Forest、木で異質効果を信頼区間つきで切る

0
Posted at

この連載について

因果推論を「本を読んだ」で終わらせず、自分の言葉で説明でき、コードで再現できる状態まで落とす30日連載です。前回のDay 24では、アウトカムと処置をそれぞれ機械学習で予測して残差を取り、残差同士の回帰と交差適合で「機械学習の当てはめ癖」を打ち消すDouble MLを組みました。ただしDMLが返すのは平均効果です。今日はその残差化を木と組み合わせ、Day 23で点推定止まりだった条件付き処置効果 $\tau(x)$ に信頼区間を付けるCausal Forestを組みます。図はすべて合成データで、真の $\tau(x)$ を私たちが握った状態で答え合わせをします。


TL;DR(3行)

  • 因果木は分割規準が違う。 回帰木は予測誤差、因果木は子ノード間の効果差で割ります。ベースラインの段差を $x=0$、効果の段差を $x=1$ に仕込むと、回帰木は x = 0.00、因果木は x = 1.02 を選びました。予測に良い分割≠効果を分ける分割です。
  • CausalForestDML は $\tau(x)$ を信頼区間つきで返す。 非線形な真の $\tau(x)$ に対し RMSE 0.123(X-learnerは0.242)、95%区間が真値を含む割合はオラクル検証で 96.2% でした。
  • 区間の生命線は honesty=分割用と推定用の標本分割。 分割を選んだのと同じデータで区間まで作ると、95%のはずのカバレッジが 81.3% に落ちます。標本を分ければ 94.8% に戻ります。

今日の問い

Day 23のメタラーナーは $\tau(x)$ の点推定を返すだけだった。「この層の効果は1.5」と言われても、それが $1.5 \pm 0.2$ なのか $1.5 \pm 2$ なのか分からなければ、意思決定には怖くて使えない。
$\tau(x)$ に信頼区間を付けるには、森の何を変える必要があるのか?

答えは2つの発明です。1つは分割規準を「予測誤差」から「効果の差」に取り替えること。もう1つは、分割を選ぶデータと効果を推定するデータを分けること(honesty)です。今日のゴールは、この2つがそれぞれ何を解決しているのかを言語化し、econmlのCausalForestDMLで $\tau(x)$ を区間ごと推定して、その区間が本当に約束どおり当たるのかをオラクル検証することです。


概念① 因果木:分割規準を「効果の差」に取り替える

ランダムフォレストの部品である回帰木は、葉の中の予測誤差(二乗誤差)が最も減る分割を貪欲に選びます。目的が $y$ の予測ならこれで正解です。でも私たちが欲しいのは、葉の中で効果 $\tau(x)$ が均質になる分割です。$y$ を大きく動かす変数と、$\tau$ を動かす変数は、一般に別物です。ベースラインの段差がどれだけ大きくても、処置群と対照群が同じだけ持ち上がるなら効果は一定で、そこで割る意味はありません。

では「葉の中で $\tau$ が均質になるように割る」を、どう実装するか。困ったことに、個人の効果 $\tau_i = Y_i(1) - Y_i(0)$ は誰についても観測できません(Day 2の根本問題、Day 23で「教師ラベルが無い」として再来した壁です)。Athey & Imbens (2016) の因果木はここを裏返して、子ノード間で推定効果の差が最大になる分割を選びます。

$$
\max_{s}\ \frac{n_L , n_R}{n^2} \bigl(\hat\tau_L - \hat\tau_R\bigr)^2
$$

$\hat\tau_L, \hat\tau_R$ は分割後の左右の子ノードで計算した処置群と対照群の平均差です。ラベル $\tau_i$ が無くても、ノード単位の効果推定なら計算できます。その差が大きい分割ほど異質性を捉えている、という理屈です。この規準の取り替えが第1の発明です。

Causal Forestはこの因果木をサブサンプリングで多数束ねたもので、新しい点 $x$ の効果は「$x$ と同じ葉に落ちた標本ほど重く数える」適応的な近傍平均として推定されます。Wager & Athey (2018) は、この推定量が漸近正規($\hat\tau(x)$ が正規分布に近づく)であることを示し、分散推定と合わせて点ごとの信頼区間を出せるようにしました。今日使う econml の CausalForestDML は、その前段にDay 24の残差化を挟みます。$y$ と $w$ をそれぞれニュイサンスモデルで予測して残差 $\tilde y, \tilde w$ を作り、交絡の掃除はDMLに任せ、森は残差の世界で異質性だけを担当する分業です。

概念② honesty:分割を選んだデータに、区間まで作らせない

第2の発明は地味に見えて、信頼区間の生命線です。因果木の規準は「効果の差が最大に見える」分割を探します。ということは、たとえ真の効果が完全に一様でも、ノイズがたまたま作った差を必ずどこかで見つけてしまいます。その同じデータで葉の効果を推定すると、選ばれた分割の「見かけの差」をそのまま真に受けることになり、葉の推定値は外側へ偏ります。勝者の呪い(winner's curse)です。点推定が偏れば、その周りに張った95%信頼区間は真値を外し、カバレッジは名目を割ります。

honestyは、これを標本分割で断ち切ります。標本の半分で木の構造(どこで割るか)だけを決め、残りの半分で葉の中の効果を推定します。構造の選択と推定が独立になるので、「選んだから偏る」が起きません。Wager & Athey の漸近正規性の証明も、このhonestyを前提に組み立てられています。代償は、それぞれの仕事にデータを半分しか使えないことによる分散の増加です。この損得は「手を動かす③」で数値にします。


手を動かす①:回帰木と因果木、同じデータで最良分割を比べる

概念①を最小の実験で確かめます。共変量は1つ、割付は完全ランダム(RCT)。ベースラインは $x=0$ で大きく段差を作り(両群とも+4、効果は不変)、処置効果は $x=1$ で初めて0から2に変わる、という意地悪なデータです。深さ1の「切り株」を、2つの規準で育てます。

import numpy as np

rng = np.random.default_rng(25)
n = 2000
x = rng.uniform(-2, 2, n)              # 共変量は1つだけ
w = rng.binomial(1, 0.5, n)            # RCT(割付は完全ランダム)
mu = 4.0 * (x > 0)                     # ベースライン:x=0 に大きな段差
tau = 2.0 * (x > 1)                    # 処置効果:x=1 で初めて変わる(真のCATE)
y = mu + tau * w + rng.normal(0, 1, n)

def sse(v):
    return ((v - v.mean()) ** 2).sum()

splits = np.linspace(-1.8, 1.8, 181)
for name, causal in [("回帰木(予測誤差で分割)", False), ("因果木(効果の差で分割)", True)]:
    scores = []
    for s in splits:
        L, R = x <= s, x > s
        if min(L.sum(), R.sum()) < 50:
            scores.append(-np.inf)
            continue
        if not causal:
            scores.append(sse(y) - sse(y[L]) - sse(y[R]))   # 予測誤差の減少量
        else:
            tL = y[L & (w == 1)].mean() - y[L & (w == 0)].mean()
            tR = y[R & (w == 1)].mean() - y[R & (w == 0)].mean()
            scores.append(L.sum() * R.sum() / n**2 * (tL - tR) ** 2)  # 子ノード間の効果差
    print("%sの最良分割点: x = %.2f" % (name, splits[int(np.argmax(scores))]))
回帰木(予測誤差で分割)の最良分割点: x = 0.00
因果木(効果の差で分割)の最良分割点: x = 1.02

回帰木は x = 0.00 を選びました。$y$ の段差が最大の場所です。でもそこはベースラインが動くだけで、効果は両側とも同じです。効果の異質性を知りたい私たちにとって、この分割は1歩も前進していません。因果木は x = 1.02 を選びました。真の効果が0から2に変わる場所、まさに知りたかった段差です。同じデータ、同じ木のアルゴリズムでも、規準を取り替えるだけで見るものが変わります。Day 23でS-learnerの正則化が $\tau$ を潰した話と根は同じで、予測の目的関数は $\tau$ に興味がないのです。

Fig 1: 分割規準の対比(左:yの段差はx=0、効果の段差はx=1にあるデータ。右:回帰木の規準はx=0.00、因果木の規準はx=1.02にピーク)


手を動かす②:CausalForestDMLで τ(x) を推定し、区間ごとオラクル検証する

本番です。共変量5つの合成データで、真の効果は $x_1$ のS字カーブ $\tau(x) = 1 + 2/(1+e^{-3x_1})$(1から3へ滑らかに変化)。しかも $x_1$ は傾向スコアにも効く交絡因子にします。RCTではない観察データで、前段のDML残差化が仕事をする設定です。真の $\tau$ を知っている評価用データ2,000点で、点推定の精度(RMSE)と95%区間のカバレッジを答え合わせします。X-learner(Day 23)も同じデータで比較します。実行は手元で20秒ほどです。

import numpy as np
from econml.dml import CausalForestDML
from econml.metalearners import XLearner
from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor
from sklearn.linear_model import LogisticRegression

def make_data(n, rng):
    X = rng.uniform(-2, 2, (n, 5))              # 共変量5つ(効くのは一部だけ)
    tau = 1 + 2 / (1 + np.exp(-3 * X[:, 0]))    # 真の効果:x1のS字(1〜3で変化)
    e = 1 / (1 + np.exp(-X[:, 0]))              # 傾向スコア:x1が割付にも効く交絡
    w = rng.binomial(1, e)
    y = 2 * X[:, 0] + X[:, 1] + tau * w + rng.normal(0, 1, n)
    return X, w, y, tau

rng = np.random.default_rng(25)
X, w, y, _ = make_data(4000, rng)             # 学習用
Xte, _, _, tau_te = make_data(2000, rng)      # 評価用(真のτで答え合わせ)

rf = dict(n_estimators=200, min_samples_leaf=20, random_state=0)
cf = CausalForestDML(
    model_y=RandomForestRegressor(**rf), model_t=RandomForestClassifier(**rf),
    discrete_treatment=True, n_estimators=500, min_samples_leaf=10,
    cv=2, random_state=25)
cf.fit(y, w, X=X)                             # 前段:残差化(Day 24)→ 後段:因果森
tau_cf = cf.effect(Xte)                       # τ(x) の点推定
lo, hi = cf.effect_interval(Xte, alpha=0.05)  # 95%信頼区間

xl = XLearner(models=RandomForestRegressor(**rf), propensity_model=LogisticRegression())
xl.fit(y, w, X=X)

def rmse(t):
    return np.sqrt(((t - tau_te) ** 2).mean())

cover = ((lo <= tau_te) & (tau_te <= hi)).mean()
print("Causal Forest: RMSE = %.3f、95%%区間の平均幅 = %.2f" % (rmse(tau_cf), (hi - lo).mean()))
print("X-learner    : RMSE = %.3f(区間なし)" % rmse(xl.effect(Xte)))
print("オラクル検証 : 95%%区間が真のτ(x)を含む割合 = %.1f%%(評価2000点)" % (100 * cover))
Causal Forest: RMSE = 0.123、95%区間の平均幅 = 0.59
X-learner    : RMSE = 0.242(区間なし)
オラクル検証 : 95%区間が真のτ(x)を含む割合 = 96.2%(評価2000点)

点推定の精度はCausal Forestが RMSE 0.123、X-learnerが 0.242 でした。下の図で見ると、どちらも真のS字を追えていますが、X-learnerの曲線は細かく暴れ、Causal Forestは滑らかに真値をなぞります。交絡のある観察データなので、前段の残差化が効いている形です。そして本題の区間。95%を約束した信頼区間は、評価2,000点のうち 96.2% で真の $\tau(x)$ を捕まえました。名目どおりです。「この点の効果は $2.1 \pm 0.3$」という言い方が、合成データの上ではきちんと担保されて出てくるわけです。これがメタラーナーとの決定的な違いです。

Fig 2: 真のτ(x)・Causal Forest(95%区間帯つき)・X-learnerの重ね描き(x1の断面。CFはRMSE 0.123で区間カバレッジ96.2%、X-learnerはRMSE 0.242で区間なし)


手を動かす③:honestyを外すと、95%の区間は81%しか当たらない

区間が当たったのはhonestyのおかげだ、と概念②で述べました。それを裏返しで確かめます。econmlの内部を開く代わりに、深さ1の因果木を自分で組んで、「分割も推定も同じデータ」(honestyなし)と「半分で分割を選び、残り半分で推定」(honestyあり)を2,000回のモンテカルロで比べます。意地悪なのはデータの作り方で、真の効果は全員1.0、異質性はゼロです。それでも因果木の規準は、ノイズがたまたま作った「効果の差が最大に見える」分割を毎回どこかに見つけます。実行は30秒ほどです。

import numpy as np

def best_split(x, w, y, splits, min_leaf=50):
    best_s, best_v = None, -np.inf
    for s in splits:
        L, R = x <= s, x > s
        if min(L.sum(), R.sum()) < min_leaf:
            continue
        tL = y[L & (w == 1)].mean() - y[L & (w == 0)].mean()
        tR = y[R & (w == 1)].mean() - y[R & (w == 0)].mean()
        v = L.sum() * R.sum() / len(x)**2 * (tL - tR) ** 2
        if v > best_v:
            best_v, best_s = v, s
    return best_s

def leaf_ci(x, w, y, s):
    out = []
    for m in (x <= s, x > s):
        y1, y0 = y[m & (w == 1)], y[m & (w == 0)]
        t = y1.mean() - y0.mean()
        se = np.sqrt(y1.var(ddof=1) / len(y1) + y0.var(ddof=1) / len(y0))
        out.append((t - 1.96 * se, t + 1.96 * se))
    return out

rng = np.random.default_rng(2525)
R, n = 2000, 400
splits = np.linspace(-1.6, 1.6, 33)
hit = {"adaptive": [], "honest": []}
for _ in range(R):
    x = rng.uniform(-2, 2, n)
    w = rng.binomial(1, 0.5, n)
    y = 1.0 * w + rng.normal(0, 1, n)      # 真の効果は全員 1.0(異質性なし)

    s = best_split(x, w, y, splits)        # 分割も推定も同じデータ
    hit["adaptive"] += [lo <= 1.0 <= hi for lo, hi in leaf_ci(x, w, y, s)]

    half = rng.permutation(n) < n // 2     # 半分で分割を選び、残り半分で推定
    s2 = best_split(x[half], w[half], y[half], splits)
    hit["honest"] += [lo <= 1.0 <= hi
                      for lo, hi in leaf_ci(x[~half], w[~half], y[~half], s2)]

print("honestyなし(同じデータで分割と推定): 95%%区間のカバレッジ = %.1f%%"
      % (100 * np.mean(hit["adaptive"])))
print("honestyあり(分割用と推定用を分ける): 95%%区間のカバレッジ = %.1f%%"
      % (100 * np.mean(hit["honest"])))
honestyなし(同じデータで分割と推定): 95%区間のカバレッジ = 81.3%
honestyあり(分割用と推定用を分ける): 95%区間のカバレッジ = 94.8%

差は歴然です。honestyなしでは、95%を名乗る区間が 81.3% しか当たりません。分割の規準が「差が最大に見える場所」を選んだ時点で、左右の葉の推定値はノイズの分だけ外側へ引き離されていて、その偏った点推定を中心に張った区間は真値1.0を外しやすくなります。しかも区間の計算式自体は教科書どおり正しいので、壊れていることが式からは見えません。標本を分けるだけでカバレッジは 94.8% に戻ります。推定に使える標本は半分になっているのに、です。区間を読むつもりがあるなら、honestyは省略できる工程ではありません。econmlのCausalForestDMLは honest=True が既定値で、手を動かす②の96.2%はこの仕組みの上に立っています。


つまづき・誤解しやすい点

  • 「予測精度の高い森=良い因果森」ではありません。 $\tau_i$ は観測できないので、テストデータの精度で因果モデルは選べません(Day 23のPEHEは真値を知る合成データ専用の物差しです)。現実データでの選択は、反証テスト(Day 15)や、次回Day 26のアップリフト評価のような間接的な物差しに頼ります。今日のRMSEやカバレッジの答え合わせは、真値を握ったオラクルだからできた贅沢だと覚えておいてください。
  • honestyはタダではありません。 分割にも推定にも標本の半分しか使えないぶん、点推定は粗くなります。区間が不要で点推定の精度だけが欲しい場面なら、honestyなしの適応的な木にも一理あります。ただしその森の effect_interval を読んではいけません。手を動かす③のとおり、区間の顔をした別の何かが出てきます。
  • 区間は点ごと(pointwise)で、しかも標本誤差しか語りません。 96.2%は「各点の区間がそれぞれの真値を含む割合」で、曲線全体を同時に覆う帯ではありません。そして区間の幅に入っているのはサンプリングの不確かさだけです。非交絡性(Day 6)が破れていれば、区間ごと系統的にズレます。隠れ交絡への耐性は信頼区間ではなく感度分析(Day 9)の仕事です。

GISデータ実務での使い方

  • 「どの区画に効くか」を区間つきで語れます。 漏水対策や管路更新の効果は、管路材質・布設年代・土壌・水圧で大きく異質なはずです。平均効果(Day 24)で施策の是非を決めたあと、Causal Forestで $\tau(x)$ を推定すれば、「効果の下限がプラスと言い切れる区画」から優先着手する、という守りの効いた意思決定に進めます。これは次回Day 26のターゲティングの入口です。
  • 区間の広さは「次にデータを取る場所」を教えてくれます。 区間が広い区画は「効果が無い」のではなく「まだ分からない」区画です。効果の点推定だけを地図に塗ると両者の区別が消えます。区間幅の地図を並べて、次の実証実験の対象選びに使うのが誠実な運用です。
  • 空間データではhonestyでも足りないことがあります。 honestyが守ってくれるのは「同じデータで選んで推定する」偏りです。隣接する区画のデータが空間相関で似ている場合、ランダムに半分に割っても2つの標本は独立になりません。区間は楽観的に狭くなりえます。分割をランダムでなく空間ブロック単位にする発想(空間CVと同じ処方)が必要で、ここはDay 28でまとめて掘ります。

参考(本棚)

次回予告(Day 26)

$\tau(x)$ を区間つきで推定できるようになりました。でも実務の問いはもう一歩先にあります。クーポンを配る予算が限られているとき、誰に配れば増分が最大になるか。推定した効果で顧客を並べ替え、その並べ替えの良し悪しを測る道具がQini曲線です。次回Day 26はアップリフトモデリング。「反応率が高い人」と「効果が高い人」の違いから始めて、4象限(説得可能・鉄板・無関心・天邪鬼)とターゲティングの価値を数値にします。

0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?