ベイズ因果推論｜因果推論

← 因果推論一覧

🎓 レベル：発展　|　重要度：B（標準）

📎 前提：回帰による調整とその限界　|　識別の仮定　|　数理：ベイズ線形回帰（ベイズ）・事前・尤度・事後・周辺尤度（ベイズ）

🔧 環境：PyMC 6.0.1 / ArviZ 1.2.0（Windows）。pm.sample は cores=1、事後の数値は idata.posterior[...].values から取り出します。

要点（BLUF）

識別はベイズでも一切変わらない。「 $\tau$ が ATE である」と言えるのは、バックドアが閉じている（条件付き交換可能性・正値性・SUTVA、識別の仮定）ときだけです。事前分布は識別を救いません。ベイズが担うのは識別の後の「推定と不確実性の表現」です。
推定：処置係数 $\tau$ に事前分布を置いてベイズ回帰を行うと、ATE の点推定でなく事後分布 $p(\tau\mid\text{data})$ が得られます。フラットに近い事前なら事後平均は OLS 点推定とほぼ一致し、加えて信用区間が不確実性をそのまま返します。
真の $\text{ATE}=2.0$ を仕込んだ交絡データで、素朴比較が外れ、ベイズ回帰の事後平均が真値・89% 信用区間が真値を覆うことを PyMC で確認します。

1. 識別が先、ベイズは後

因果効果を「点」でなく「分布」で受け取れるのがベイズの魅力ですが、最初に強調すべきは順序です。

回帰による調整とその限界で見たとおり、回帰の処置係数 $\tau$ が ATE に等しくなるのは、調整集合 $C$ がバックドア基準を満たし、

\tau \;=\; E_C\big[\,E[Y\mid X{=}1,C]-E[Y\mid X{=}0,C]\,\big]

という g公式が成り立つときだけです。これは識別の問題で、データの量や推論の流儀（頻度論かベイズか）とは無関係に決まります。未観測交絡があれば、いくら精緻な事後分布を描いても、その分布は「ATE の分布」ではなく「交絡で歪んだ係数の分布」にすぎません。ベイズの事後分布は識別の正しさを保証しないし、事前分布で交絡を消すこともできません。

❓ ここで一度立ち止まりましょう：「事後分布が真値の周りに集中したら、識別できている証拠だろうか？」違います。識別が崩れていれば、事後分布は間違った値の周りに自信たっぷりに集中します。狭い信用区間は「精度が高い」ことを意味しても「正しい」ことは意味しません。

したがって本ノートは「 $C$ で条件づければバックドアが閉じる」という前提のもとで、その先の推定をベイズで行う話に限定します。

2. ベイズ回帰モデル：処置係数に事前分布を置く

処置 $X\in\{0,1\}$ 、結果 $Y$ 、交絡 $C$ に対し、条件付き期待値が線形だとして次のモデルを置きます（尤度）。

Y_i \mid X_i, C_i \;\sim\; \mathcal{N}\!\big(\alpha + \tau X_i + \beta C_i,\ \sigma^2\big)

ここに事前分布を与えます。係数には弱情報事前、尺度には半正規分布を使います。

\alpha,\ \tau,\ \beta \;\sim\; \mathcal{N}(0,\ 5^2), \qquad \sigma \;\sim\; \text{HalfNormal}(2)

ベイズの定理より事後分布は尤度と事前の積に比例します。

p(\alpha,\tau,\beta,\sigma \mid \text{data}) \;\propto\; p(\text{data}\mid \alpha,\tau,\beta,\sigma)\;p(\alpha)p(\tau)p(\beta)p(\sigma)

私たちが知りたいのは ATE、すなわち $\tau$ の周辺事後分布 $p(\tau\mid\text{data})$ です。識別の仮定が成り立つ前提のもとで、この分布が「ATE がどのあたりにありそうか」をまるごと表します。事前・尤度・事後の関係は事前・尤度・事後・周辺尤度（ベイズ）、線形回帰のベイズ版はベイズ線形回帰（ベイズ）が土台です。

3. PyMC で事後分布を得る

真の効果を仕込んだ擬似データを作ります。交絡 $C$ が処置の受けやすさと結果の両方を押し上げ、真の処置効果は均一で $\text{ATE}=2.0$ とします。素朴比較・ベイズ回帰の順に当てます。

import numpy as np
import pandas as pd
import pymc as pm
import arviz as az

# === 真のATEを仕込んだ交絡データを作り、ベイズ回帰で処置効果の事後分布を得る ===
rng = np.random.default_rng(20260627)
n = 800

# 交絡 C: 処置の受けやすさと結果の両方を押し上げる共通原因
C = rng.normal(0.0, 1.0, size=n)

# 処置割り当て: C が高いほど処置されやすい(選択あり)
prob_treat = 1.0 / (1.0 + np.exp(-0.8 * C))
X = rng.binomial(1, prob_treat)

# 結果: 真の処置効果は均一で ATE_true = 2.0、交絡 C も正に効く
ATE_true = 2.0
Y = 1.0 + ATE_true * X + 1.5 * C + rng.normal(0.0, 1.0, size=n)

# (1) 素朴比較(調整なし) ... 交絡で過大評価
naive = Y[X == 1].mean() - Y[X == 0].mean()
print(f"真の ATE        : {ATE_true:.2f}")
print(f"(1) 素朴比較     : {naive:.2f}  (交絡で過大)")

# (2) ベイズ回帰: Y = alpha + tau*X + beta*C + 正規誤差。tau に弱情報事前を置く
with pm.Model() as model:
    alpha = pm.Normal("alpha", mu=0.0, sigma=5.0)
    tau = pm.Normal("tau", mu=0.0, sigma=5.0)      # 処置効果(知りたい量)
    beta = pm.Normal("beta", mu=0.0, sigma=5.0)    # 交絡 C の係数
    sigma = pm.HalfNormal("sigma", sigma=2.0)
    mu = alpha + tau * X + beta * C
    pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
    idata = pm.sample(
        draws=1000, tune=1000, chains=4, cores=1,
        random_seed=42, progressbar=False,
    )

# 事後は文字列で返る az.summary ではなく posterior の生値から計算する
tau_draws = idata.posterior["tau"].values.reshape(-1)
post_mean = tau_draws.mean()
lo, hi = np.quantile(tau_draws, [0.055, 0.945])   # 89%等裾区間(ETI)
n_div = int(idata.sample_stats["diverging"].values.sum())

print(f"(2) 事後平均 tau : {post_mean:.2f}")
print(f"    89%信用区間  : [{lo:.2f}, {hi:.2f}]")
print(f"    真値2.0を覆う : {lo <= ATE_true <= hi}")
print(f"    発散数        : {n_div}")

実行すると次のように印字されます。

真の ATE        : 2.00
(1) 素朴比較     : 3.01  (交絡で過大)
(2) 事後平均 tau : 2.05
    89%信用区間  : [1.92, 2.17]
    真値2.0を覆う : True
    発散数        : 0

出力の意味：素朴比較 $3.01$ は真値 $2.0$ を上回ります（ $C$ が高い個体ほど処置され、 $C$ が $Y$ を押し上げる正のバイアス）。一方、 $C$ を入れたベイズ回帰の事後平均 $\tau$ は $2.05$ と真値にほぼ一致し、89% 信用区間 $[1.92,\ 2.17]$ は真値 $2.0$ をしっかり覆います。発散（divergence）は 0 で、サンプリングは健全です。ここで得たのは1つの数ではなく、 $\tau$ の事後分布そのもの——「ATE は 2 前後で、おおむね 1.9〜2.2 の幅にありそう」という不確実性込みの答えです。

補足：事後の数値は az.summary が文字列で返すため、計算は必ず idata.posterior["tau"].values から取り出します。区間はここでは 89% の**等裾区間（ETI）**を分位点で計算しています。

4. OLS は「点」、ベイズは「分布」

ベイズの事後分布を、OLS の点推定と重ねて可視化します。同じデータに OLS（Y ~ X + C）を当て、その処置係数を縦線で、事後分布をヒストグラムで描きます。

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import statsmodels.formula.api as smf
import pandas as pd
import pymc as pm

# === OLSは「点」、ベイズは「分布」: 処置効果の不確実性を可視化する ===
rng = np.random.default_rng(20260627)
n = 800
C = rng.normal(0.0, 1.0, size=n)
prob_treat = 1.0 / (1.0 + np.exp(-0.8 * C))
X = rng.binomial(1, prob_treat)
ATE_true = 2.0
Y = 1.0 + ATE_true * X + 1.5 * C + rng.normal(0.0, 1.0, size=n)
df = pd.DataFrame({"Y": Y, "X": X, "C": C})

# OLSの点推定(処置係数)
tau_ols = smf.ols("Y ~ X + C", data=df).fit().params["X"]

# ベイズ回帰の事後サンプル
with pm.Model() as model:
    alpha = pm.Normal("alpha", mu=0.0, sigma=5.0)
    tau = pm.Normal("tau", mu=0.0, sigma=5.0)
    beta = pm.Normal("beta", mu=0.0, sigma=5.0)
    sigma = pm.HalfNormal("sigma", sigma=2.0)
    mu = alpha + tau * X + beta * C
    pm.Normal("Y_obs", mu=mu, sigma=sigma, observed=Y)
    idata = pm.sample(draws=1000, tune=1000, chains=4, cores=1,
                      random_seed=42, progressbar=False)

tau_draws = idata.posterior["tau"].values.reshape(-1)
lo, hi = np.quantile(tau_draws, [0.055, 0.945])

fig, ax = plt.subplots(figsize=(7, 4.2))
ax.hist(tau_draws, bins=40, density=True, color="tab:blue", alpha=0.55,
        label="処置効果 tau の事後分布")
ax.axvspan(lo, hi, color="tab:blue", alpha=0.12, label="89%信用区間")
ax.axvline(ATE_true, color="green", lw=2.5, label="真のATE = 2.0")
ax.axvline(tau_ols, color="red", lw=2.0, ls="--", label="OLS点推定")
ax.set_xlabel("処置効果 tau")
ax.set_ylabel("事後密度")
ax.set_title("OLSは1点、ベイズは効果の分布で不確実性まで返す")
ax.legend(fontsize=9)
fig.tight_layout()
plt.show()

print(f"OLS点推定 tau_ols = {tau_ols:.2f}")
print(f"事後89%区間 = [{lo:.2f}, {hi:.2f}]")

印字は次のとおりです。

OLS点推定 tau_ols = 2.04
事後89%区間 = [1.92, 2.17]

出力の意味：OLS の点推定 $2.04$ は、事後分布のほぼ中央に落ちます。弱情報事前のもとでは、事後平均は OLS とほぼ一致するのです（最尤推定＝平坦な事前のベイズ、という対応）。違いは「返ってくるもの」です。OLS は1点（と別途の標準誤差）を返すのに対し、ベイズは $\tau$ の分布まるごとを返し、89% 信用区間や「 $\tau>2.5$ である事後確率」といった問いに直接答えられます。図では真値（緑）が事後分布の山の中に収まり、不確実性の幅が一目で分かります。

5. 事前の役割・縮小・不確実性の伝播

ベイズが OLS に対して持つ実利は3つです。

事前知識の取り込み：過去の研究で「効果はせいぜい ±5 程度」と分かっていれば $\tau\sim\mathcal{N}(0,5^2)$ のように反映できます。事前が弱ければデータが支配し、結果は OLS に近づきます。
縮小（shrinkage）による安定化：サンプルが小さい・共変量が多い・群が細かく分かれるとき、平坦な最尤推定は暴れます。事前が推定値を妥当な範囲に引き戻すことで分散が下がり、汎化が安定します。これは正則化と同じ働きで、リッジ回帰は「係数に正規事前を置いた MAP 推定」に対応します（正則化との関係（ベイズ）・縮小の数理は収縮の数理（ベイズ））。階層モデルにすれば、サブグループごとの処置効果を互いに情報共有しながら推定でき、群が小さいほど全体平均へ強く縮みます。
不確実性の伝播：効果を点でなく分布で持つので、ATE の不確実性をそのまま下流の意思決定（期待効用・損失）へ流し込めます。「効果が正である事後確率」「効果が閾値を超える確率」を、 $\tau$ の事後サンプルを数えるだけで計算できます。

ただしどれも識別とは別の話である点を繰り返します。事前は「同じデータからどう学ぶか」を変えますが、「データが何を識別できるか」は変えません。交絡が残っていれば、縮小も不確実性の伝播も、間違った中心の周りで起こります。

6. 仮定の直観的意味：なぜ事前は識別を救えないのか

事後分布は尤度 $p(\text{data}\mid\theta)$ を通してデータと結びつきます。ところが尤度はバックドアパスを区別できません。未観測交絡 $U$ があると、観測データの分布は「 $X$ の因果効果が $\tau$ 」のモデルでも「 $\tau$ は小さいが $U$ が効いている」モデルでも同じように説明できてしまいます。両者は**観測上は区別不能（observationally equivalent）**で、尤度が同じ値を返すため、どんな事前を置いてもデータは真の $\tau$ の方向に情報を与えません。これが「事前は識別を救わない」の数理的な中身です。

ベイズが本当に効くのは、識別が済んだあとです。バックドアが閉じていれば尤度は真の $\tau$ にピークを持ち、そこへ事前が穏やかな正則化を加え、事後が不確実性込みの答えを返します。識別（第1〜3章）と推定（本ノート）を必ず分ける——これがベイズ因果推論を誤用しないための鉄則です。未観測交絡が疑われるなら、事前ではなく感度分析（未観測交絡への感度分析）で「どれほどの交絡があれば結論が覆るか」を調べます。

⚠️ よくある誤解・落とし穴

「事後が狭い＝正しい」ではない。信用区間の狭さは精度であって正確さではありません。識別が崩れていれば、狭い区間が間違った値を自信たっぷりに覆います。
事前で交絡を消そうとしない。「効果はゼロに近いはず」という事前を入れても、それは識別の代用になりません。交絡は事前ではなく調整集合かデザインで対処します。
MCMC の健全性を必ず確認。発散数（sample_stats["diverging"]）が多い、 $\hat R$ が 1 から離れる、有効サンプル数が小さいときは事後分布を信用してはいけません。本ノートのモデルは線形・低次元なので発散 0 でしたが、階層モデルや非線形では要注意です（収束診断はハミルトニアンモンテカルロとNUTS（ベイズ））。
事前の感度を見る。結論が事前の選び方で大きく動くなら、データが弱いサインです。弱情報事前・無情報に近い事前で結果が安定するかを確認します。