この連載について
因果推論を「本を読んだ」で終わらせず、自分の言葉で説明でき、コードで再現できる状態まで落とす30日連載です。前回のDay 28では、隣の自治体に効果が漏れる世界を扱いました。空間的に滑らかな未観測交絡と、SUTVAを破るスピルオーバーです。今日は道具ではなく語り方を変えます。効果を1つの点と星印で報告するのをやめて、分布のまま報告したら何が言えるようになるのか。図と数値はすべて合成データで、真の効果を私たちが握った状態で答え合わせをします。
TL;DR（3行）
- ベイズにしても識別の仮定は1ミリも変わりません。 変わるのは不確実性の語り方だけです。交絡を残したままベイズにすれば、間違った値のまわりに狭い事後分布が立つだけで、むしろ質が悪くなります。
- 事後分布は「効果が◯を超える確率」に直接答えます。 今日の小さなRCTは頻度論で「効果 3.20、p = 0.040」ですが、同じデータのベイズ推定は P(τ>0) = 97.6%、P(τ>2) = 78.5% と言えます。ただし95%区間は 0.04 から 6.19 と広く、「2から4の間」と言い切れるほどのデータは実はありません。
- 小標本の自治体群には階層モデルが効きます。 12自治体の生の推定は -2.89 から +8.85 まで散らばりますが、階層推定は -0.85 から +5.30 に収まり、真値に対するRMSEは 1.87 から 0.62 へ。全部バラバラでも全部同じでもない、その中間をデータに決めさせます。
今日の問い
会議で2つの報告が並びます。「施策の効果は 3.2、p は 0.04 でした」。「効果はおそらく 2 から 4 の間です」。
来年度の予算を決めるのに使えるのは、どちらの言い方か？
前者は正しい報告ですが、聞いた側が知りたいこと（この施策に投資して割に合う確率）には答えていません。後者は答えているように聞こえますが、根拠を示さないと単なる感想です。今日のゴールは、後者を根拠つきで言えるようにすること、そしてその代償に何を支払っているのかを正確に把握することです。使うのは PyMC です。
概念① 事後分布としての因果効果：効果そのものを確率変数として扱う
頻度論では、真の効果 $\tau$ は固定された未知の定数です。確率が宿るのはデータのほうで、95%信頼区間が意味するのは「同じ手続きを無限に繰り返せば、そのうち95%の区間が真値を含む」という手続きの成績です。目の前の1本の区間が真値を含む確率は0か1のどちらかで、それを95%とは呼べません。p値も同じで、「真の効果がゼロだったとして、これほど極端な差が出る確率」であり、効果がゼロである確率ではありません。
ベイズでは、$\tau$ そのものに確率分布を置きます。手元のデータを見た後の分布、事後分布は次の形で書けます。
$$
p(\tau \mid \text{data}) ;\propto; p(\text{data} \mid \tau), p(\tau)
$$
右辺の第1項が尤度（データがどれだけその $\tau$ を支持するか）、第2項が事前分布（データを見る前に持っていた情報）です。この分布が手に入れば、「$\tau$ が2を超える確率」も「$\tau$ が更新コストを上回る確率」も、面積を測るだけで答えが出ます。意思決定の言葉にそのまま乗るのは、この形です。
Imbens & Rubin 第8章は、この見方を潜在的結果の言葉で整理しています。因果推論は結局、$Y_i(0)$ と $Y_i(1)$ の片方が必ず欠測している欠測データ問題です。欠測した側をモデルの事後予測分布から埋めれば、$\tau = \frac{1}{N}\sum_i (Y_i(1) - Y_i(0))$ という量そのものが分布を持ちます。Day 14 で、二値アウトカムのアブダクションはノイズを点まで戻せず範囲までしか絞れない、だから反事実は確率になると書きました。あの「ノイズの事後分布を丸ごと扱う」見方の、集団版がこれです。
ここで最初に釘を刺しておきます。ベイズにしても、識別の仮定は何ひとつ変わりません。 非交絡性（Day 6）、正しい調整集合（Day 11）、除外制約（Day 17）、SUTVA（Day 28）。これらはすべてモデルの外にある約束事で、事前分布に書いて解決できるものではありません。交絡が残ったままベイズ推定を回せば、バイアスした値のまわりに、いかにも精密そうな狭い事後分布が立ちます。むしろ危険です。
概念② 事前分布と階層モデル：データの外から情報を入れる2つの入り口
ベイズが持ち込む情報の入口は2つあります。1つが事前分布です。過去の実証研究、物理的な上限、担当者の経験。これらを $p(\tau)$ の形で書けば、少ないデータを補強できます。ただし補強は諸刃の剣で、標本が小さいほど事前分布が答えを支配します。だから事前分布は「置いたら終わり」ではなく、振ってみて感度を報告するまでが1セットです。手を動かす②でその振れ幅を測ります。
もう1つが階層モデルです。12の自治体でパイロットをやったとします。極端な立場が2つあります。「自治体はどれも別物だから、それぞれ独立に推定する」（プーリングなし）と、「どこも同じ施策だから、全部まとめて1つの効果を推定する」（完全プーリング）です。前者は小さな自治体で推定が暴れ、後者は本当にある自治体差を消してしまいます。
階層モデルは、この2つの間をデータに決めさせます。各自治体の効果 $\tau_j$ が共通の分布から生まれたと考えるのです。
$$
\hat{\tau}_j \mid \tau_j \sim \mathcal{N}(\tau_j,, \mathrm{se}j^2), \qquad
\tau_j \sim \mathcal{N}(\mu,, \sigma\tau^2)
$$
$\sigma_\tau$ は自治体間のばらつきで、これもデータから推定します。結果として各自治体の推定は、自分の生の推定と全体平均の加重平均におおよそ落ち着きます。
$$
\hat{\tau}_j^{\text{階層}} ;\approx; \omega_j \hat{\tau}_j + (1 - \omega_j)\mu,
\qquad
\omega_j = \frac{1/\mathrm{se}_j^2}{1/\mathrm{se}j^2 + 1/\sigma\tau^2}
$$
自分の推定が精密（$\mathrm{se}_j$ が小さい）なら $\omega_j$ は1に近づいて自前の値がほぼそのまま残り、粗いなら全体平均へ引き寄せられます。この引き寄せが部分プーリング、縮む現象が**縮約（shrinkage）**です。Day 22 では年齢5分位ごとにGATEを別々に推定しましたが、あれは群ごとの独立推定、つまりプーリングなしでした。群が小さくなるほど、その素直さは高くつきます。
手を動かす①：小さなRCTを PyMC で、P(τ>0) を直接言う
60地区を30ずつに無作為割付したパイロットを合成します。真の効果は月あたり3.0トンの節水、地区ごとのばらつきは6.0です。Day 4と5で扱ったのと同じ、小さな無作為化実験です。あのときはFisherの並べ替え検定とNeymanの繰り返しサンプリングという2つの語り方を並べました。今日は3つ目の語り方を足します。まず頻度論で検定し、次に同じデータをベイズで推定します。
import numpy as np
import pymc as pm
from scipy import stats
rng = np.random.default_rng(745)
n = 30 # 各群30地区の小さなRCT
y0 = rng.normal(50, 6, n) # 対照群：施策なし
y1 = rng.normal(53, 6, n) # 処置群：真の効果 τ = 3.0（月あたり節水量トン）
# (1) 頻度論：差を検定する
diff = y1.mean() - y0.mean()
se = np.sqrt(y1.var(ddof=1) / n + y0.var(ddof=1) / n)
p = 2 * stats.t.sf(abs(diff / se), 2 * n - 2)
half = stats.t.ppf(0.975, 2 * n - 2) * se
# (2) ベイズ：τ を確率変数として、分布ごと求める
with pm.Model():
mu0 = pm.Normal("mu0", 50, 20) # 対照群の平均（弱い事前分布）
tau = pm.Normal("tau", 0, 10) # 処置効果（弱い事前分布）
sigma = pm.HalfNormal("sigma", 10) # 地区ごとのばらつき
pm.Normal("y0_obs", mu0, sigma, observed=y0)
pm.Normal("y1_obs", mu0 + tau, sigma, observed=y1)
idata = pm.sample(draws=2000, tune=1000, chains=2, cores=1,
random_seed=29, progressbar=False)
post = idata.posterior["tau"].values.flatten() # τ の事後サンプル4000本
print("頻度論 : 効果 %.2f（p = %.3f, 95%%CI [%.2f, %.2f]）"
% (diff, p, diff - half, diff + half))
print("ベイズ : 効果 %.2f（95%%区間 [%.2f, %.2f]）"
% (post.mean(), np.percentile(post, 2.5), np.percentile(post, 97.5)))
print("P(τ > 0) = %.3f … 効果がプラスである確率" % (post > 0).mean())
print("P(τ > 2) = %.3f … 採算ライン2トンを超える確率" % (post > 2).mean())
頻度論 : 効果 3.20（p = 0.040, 95%CI [0.14, 6.25]）
ベイズ : 効果 3.15（95%区間 [0.04, 6.19]）
P(τ > 0) = 0.976 … 効果がプラスである確率
P(τ > 2) = 0.785 … 採算ライン2トンを超える確率
点推定は 3.20 と 3.15 でほぼ同じ、区間も [0.14, 6.25] と [0.04, 6.19] でほぼ重なります。弱い事前分布のもとでは当たり前の一致で、大標本なら理論的にも一致します（ベルンシュタイン・フォン・ミーゼスの定理）。違うのは言えることです。ベイズ側は「効果がプラスである確率 97.6%」「採算ライン2トンを超える確率 78.5%」と、意思決定者が聞きたい形で答えます。p = 0.040 からこの2つは出てきません。
そしてもう1つ、正直に受け止めるべき事実があります。この結果を「効果はおそらく2から4の間」とは言えません。事後分布の95%は 0.04 から 6.19 に広がっていて、ほぼゼロの世界も6トン超の世界も、まだ十分ありえます。「p < 0.05 で有意」という報告は、分布に開いてみると意外なほど何も確定させていません。事後分布のいちばん地味な効用は、この過信の解毒かもしれません。
手を動かす②：事前分布を振ってみる、効くのは小標本のときだけ
事前分布は情報を足す道具ですが、足しすぎれば答えを作ってしまいます。同じデータに3つの事前分布を当て、さらに標本サイズを30地区から1000地区に増やして、効き目がどう変わるかを見ます。
import numpy as np
import pymc as pm
def make_rct(n, seed=745):
rng = np.random.default_rng(seed)
return rng.normal(50, 6, n), rng.normal(53, 6, n) # 真の効果 τ = 3.0
def fit(y0, y1, prior_mu, prior_sd):
with pm.Model():
mu0 = pm.Normal("mu0", 50, 20)
tau = pm.Normal("tau", prior_mu, prior_sd) # ここだけ差し替える
sigma = pm.HalfNormal("sigma", 10)
pm.Normal("y0_obs", mu0, sigma, observed=y0)
pm.Normal("y1_obs", mu0 + tau, sigma, observed=y1)
idata = pm.sample(draws=2000, tune=1000, chains=2, cores=1,
random_seed=29, progressbar=False)
return idata.posterior["tau"].values.flatten()
priors = [("弱情報 N(0,10)", 0.0, 10.0), # ほとんど何も言わない
("懐疑的 N(0,1)", 0.0, 1.0), # 「効果はせいぜい±2程度だろう」
("楽観的 N(8,1)", 8.0, 1.0)] # 「8トンは減るはずだ」
for n in (30, 1000):
y0, y1 = make_rct(n)
print("n = %d地区/群（素朴な差 %.2f）" % (n, y1.mean() - y0.mean()))
for name, mu, sd in priors:
d = fit(y0, y1, mu, sd)
print(" %-14s 事後平均 %+.2f 95%%区間 [%+.2f, %+.2f]"
% (name, d.mean(), np.percentile(d, 2.5), np.percentile(d, 97.5)))
n = 30地区/群（素朴な差 3.20）
弱情報 N(0,10) 事後平均 +3.15 95%区間 [+0.04, +6.19]
懐疑的 N(0,1) 事後平均 +0.92 95%区間 [-0.70, +2.58]
楽観的 N(8,1) 事後平均 +6.66 95%区間 [+5.04, +8.35]
n = 1000地区/群（素朴な差 2.70）
弱情報 N(0,10) 事後平均 +2.70 95%区間 [+2.19, +3.22]
懐疑的 N(0,1) 事後平均 +2.53 95%区間 [+2.03, +3.02]
楽観的 N(8,1) 事後平均 +3.05 95%区間 [+2.55, +3.56]
30地区のときは悲惨です。同じデータなのに、事前分布を替えるだけで答えは 0.92 から 6.66 まで動きました。7倍以上の開きで、懐疑的な事前を置いた人は「効果はほぼ無い」と報告し、楽観的な事前を置いた人は「7トン近く減る」と報告します。しかも後者の95%区間 [+5.04, +8.35] は真値3.0を含みません。事前分布に強い主張を書き込むと、データが黙って上書きされます。
1000地区に増やすと景色が変わります。3つの事後平均は 2.53 / 2.70 / 3.05 に収まり、区間もほぼ重なりました。データが増えれば尤度が事前を圧倒し、誰が何を信じて始めたかは結論に残りません。事前分布が効くのは、標本が小さいときだけです。逆に言えば、小標本の分析でベイズを使うなら、事前分布の感度分析を出さない報告は読む価値がありません。私自身、最初はこの2枚目のグラフを作るまで「弱情報事前なら中立」と何となく思っていましたが、左のパネルを見て考えを改めました。
手を動かす③：12自治体の効果を階層モデルで束ねる
最後が今日の本命です。12の自治体でパイロットをやりました。参加地区数はバラバラで、いちばん小さいAは12地区、いちばん大きいLは300地区です。真の自治体別効果は平均2.0・SD1.5の分布から生まれたとします（この真値は私たちだけが知っています）。自治体ごとの生の推定と、階層モデルの推定を比べます。
import arviz as az
import numpy as np
import pymc as pm
NJ = np.array([12, 16, 20, 24, 30, 40, 50, 60, 80, 120, 200, 300]) # 自治体ごとの地区数
NAMES = list("ABCDEFGHIJKL")
J = len(NJ)
rng = np.random.default_rng(261)
tau_true = rng.normal(2.0, 1.5, J) # 真の自治体別効果（平均2.0・SD1.5）
hat, resid = np.zeros(J), []
for j in range(J):
w = np.repeat([0, 1], NJ[j] // 2) # 自治体の中で地区を半々に無作為割付
y = 50 + tau_true[j] * w + rng.normal(0, 5.0, NJ[j])
y1, y0 = y[w == 1], y[w == 0]
hat[j] = y1.mean() - y0.mean() # 自治体ごとの生の推定（プーリングなし）
resid.append(np.concatenate([y1 - y1.mean(), y0 - y0.mean()]))
s = np.sqrt(np.sum(np.concatenate(resid) ** 2) / (NJ.sum() - 2 * J)) # 共通の個体差SD
se = s * np.sqrt(2.0 / (NJ // 2)) # 生の推定の標準誤差
# 階層モデル：12自治体の効果は、共通の分布から生まれたと考える
with pm.Model():
mu = pm.Normal("mu", 0, 5) # 全体の平均効果
sigma_tau = pm.HalfNormal("sigma_tau", 2) # 自治体間のばらつき
z = pm.Normal("z", 0, 1, shape=J) # 非中心化パラメータ化
tau = pm.Deterministic("tau", mu + sigma_tau * z)
pm.Normal("obs", tau, se, observed=hat) # 生の推定を、誤差 se つきの観測とみなす
# step に Slice を指定するのは、実行のたびに同じ数値を出すため（つまづき③で説明）
idata = pm.sample(draws=6000, tune=2000, chains=2, cores=1,
step=pm.Slice(), random_seed=29, progressbar=False)
post = idata.posterior["tau"].values.reshape(-1, J)
shrunk = post.mean(axis=0) # 階層推定（部分プーリング）
pooled = np.sum(hat / se**2) / np.sum(1 / se**2) # 完全プーリング
print("自治体 地区数 生の推定 階層推定 真値")
for j in [0, 1, 6, 10, 11]:
print(" %s %5d %+6.2f %+6.2f %+6.2f"
% (NAMES[j], NJ[j], hat[j], shrunk[j], tau_true[j]))
print("真値に対するRMSE 生 %.2f ／ 階層 %.2f ／ 完全プーリング %.2f"
% tuple(np.sqrt(np.mean((e - tau_true) ** 2))
for e in (hat, shrunk, np.full(J, pooled))))
print("自治体間のばらつき σ_τ の事後平均 = %.2f（完全プーリングの推定値 = %.2f）"
% (float(idata.posterior["sigma_tau"].mean()), pooled))
print("収束の確認 r_hat の最大 = %.2f ／ 実効サンプル数の最小 = %.0f"
% (float(az.rhat(idata)["tau"].max()), float(az.ess(idata)["tau"].min())))
自治体 地区数 生の推定 階層推定 真値
A 12 -2.89 -0.03 +0.24
B 16 +8.85 +5.01 +4.00
G 50 +6.74 +5.30 +6.07
K 200 +4.22 +3.99 +4.04
L 300 +0.48 +0.58 +1.02
真値に対するRMSE 生 1.87 ／ 階層 0.62 ／ 完全プーリング 1.84
自治体間のばらつき σ_τ の事後平均 = 2.42（完全プーリングの推定値 = 1.83）
収束の確認 r_hat の最大 = 1.00 ／ 実効サンプル数の最小 = 4493
12地区しかないAの生の推定は -2.89、つまり「施策のせいで水が増えた」です。真値は +0.24 なので、これは丸ごと偶然のブレです。16地区のBは逆に +8.85 で、真値 +4.00 の倍以上に出ています。この2つを額面どおり受け取れば、「A自治体では逆効果だったので中止、B自治体は大成功なので横展開」という、まったく根拠のない意思決定が生まれます。
階層モデルはAを -0.03、Bを +5.01 に引き戻しました。一方で200地区のKは +4.22 から +3.99 へ、300地区のLは +0.48 から +0.58 へと、ほとんど動いていません。縮む量は自治体ごとに違い、データが薄いところほど強く縮みます。 ここが完全プーリングとの決定的な違いで、50地区あるGは生の +6.74 から +5.30 までしか縮まず、「この自治体は本当に効いている」という信号が残りました（真値 +6.07）。
真値に対するRMSEを見ると、生の推定 1.87、完全プーリング 1.84 に対して、階層推定は 0.62 です。両極端のどちらよりも良く、しかも「どれくらい混ぜるか」は $\sigma_\tau$ の推定を通じてデータが決めています。私たちが手で決めたパラメータはありません。
つまづき・誤解しやすい点
- 「ベイズにすれば交絡が消える」わけではありません。 Day 6 で、割付を決めた要因を測れていなければ、観測共変量をいくら調整しても真値ゼロの効果が +2.15 に見えることを確かめました。あのデータをベイズで解いても、+2.15 のまわりに狭い事後分布が立つだけです。事前分布は $\tau$ の値についての知識であって、割付メカニズムについての知識ではありません。バイアスそのものにパラメータを立てて事前分布を置く道はあり、それは Day 9 の感度分析をベイズで書き直したものになります。
- 「無情報事前だから客観的」ではありません。 手を動かす②で使った N(0,10) も立派な情報で、「効果が±20を超えることはまずない」と主張しています。しかも無情報らしさは変数の取り方に依存し、$\tau$ に一様分布を置くことと $\exp(\tau)$ に一様分布を置くことは別の主張です。実務でやるべきことは客観性を名乗ることではなく、置いた事前を書き、振って感度を見せることです。
-
「事後分布が出た＝収束した」ではありません。 MCMCは乱数計算なので、最低限 r_hat（1.00に近いか）と実効サンプル数は毎回見ます。今日はもう1つ実害に当たりました。手を動かす③の階層モデルは、既定のNUTSだと同じ
random_seedでも実行のたびに小数第2位が動きます。NUTSの軌道が浮動小数点の丸め誤差を指数的に拡大するためで、そのままでは記事の数値を再現できません。そこで乱数の使い方が単純な Slice を明示指定しました（結論はNUTSと一致します）。診断の詳しい作法は参考のR-hat論文に譲ります。
GISデータ実務での使い方
- 小規模自治体のパイロットを、単独で「効果なし」と結論しません。 手を動かす③のAのように、地区数が10台なら推定が符号ごとひっくり返るのは普通に起きます。判断材料にすべきは自治体ごとの生の推定ではなく、複数自治体を階層モデルで束ねた後の推定です。「この自治体のデータだけで語れることは、実はほとんどない」と最初に言えるかどうかが分かれ目になります。
- 意思決定は閾値つきの確率で報告します。 管路更新でも漏水対策でも、現場が知りたいのは「効果があるか」ではなく「投資額を上回る効果が出る確率」です。事後分布があれば P(τ > 更新コスト相当) をそのまま出せます。p値をどれだけ丁寧に説明しても、この数字は出てきません。
- 事前分布は物理と実績から書けます。 水道分野には劣化速度や漏水率の蓄積があり、「効果がマイナス20%を下回ることはない」といった制約は根拠つきで書けます。ただし手を動かす②のとおり、小標本ほど事前が答えを支配します。書いた事前は必ず3通りくらい振って、結論が変わらないことを添えます。
- 空間はそのまま階層モデルの続きになります。 「隣り合う自治体は似た効果を持つ」という前回の構造は、$\tau_j$ の事前分布に空間相関を入れる形で書けます（CARモデルやガウス過程）。部分プーリングを「全自治体で一様に」ではなく「地理的に近いほど強く」行うわけで、Day 28 の空間の話と今日の階層モデルは、同じ枠組みの中で自然につながります。
参考（本棚）
- Imbens & Rubin『Causal Inference for Statistics, Social, and Biomedical Sciences』第8章 モデルベース推論 — 因果推論を欠測データ問題として定式化し、欠測した潜在的結果をベイズで補完する。今日の下敷き
- Gelman, Carlin, Stern, Dunson, Vehtari & Rubin『Bayesian Data Analysis』第5章 階層モデル — 8 schools の例。手を動かす③のモデルはこの形をそのまま自治体に読み替えたもの
- Martin, Kumar & Lao『Bayesian Modeling and Computation in Python』 — PyMC と ArviZ の実装作法。非中心化パラメータ化や診断の読み方はここが実務的
- Congdon『Bayesian Hierarchical Models, 2nd ed.』 — 階層モデルの応用カタログ。空間・時系列への拡張は実務節の4点目に直結する
- Vehtari, Gelman, Simpson, Carpenter & Bürkner (2021) Rank-normalization, folding, and localization: An improved R̂ for assessing convergence of MCMC. Bayesian Analysis 16(2) — つまづき③で逃がした収束診断の一次情報
- 図と数値はすべて合成データです。真の効果（RCTは3.0、自治体別は平均2.0・SD1.5の分布）を私たちが握ることで、事前分布の感度と部分プーリングの利得を答え合わせできる形にしました。生成スクリプトは元ノート側に添付しています
次回予告（Day 30）
これで道具はすべて出そろいました。Rubin流の設計と仮定、Pearl流の構造と識別、計量経済のデザイン、機械学習との合流、空間、そして今日のベイズ。次回Day 30は最終回、30日で歩いた地図を1枚に描き直します。3つの流派がどの問いに答え、どの仮定を買い、どの道具に落ちるのか。手法の選び方の流れ図と、各Dayへの索引と、次に読むものを置いて締めます。コードは書きません。書くのは、30日を通して身についた判断を3つに絞った言葉です。