この連載について
因果推論を「本を読んだ」で終わらせず、自分の言葉で説明でき、コードで再現できる状態まで落とす30日連載です。前回のDay 11では、DAGの上でバックドアパスを見つけ、調整してよい変数といけない変数を規則で選び分けました。ただ、その規則に従うと「なぜ」因果効果が出るのかという証明は今日に持ち越していました。今日は介入を式で書く道具、do 演算子を導入して、その宿題を回収します。本連載でいちばんの難所のひとつなので、言葉と式と図の三方向から同じことを言います。
TL;DR(3行)
- $P(Y \mid X{=}x)$ は「たまたま $x$ だった人を眺めた」分布、$P(Y \mid do(X{=}x))$ は「全員を $x$ に書き換えた」分布。 条件付けは母集団の絞り込み、do は割付機構の書き換え。この2つは一般に別の数字になります。
- 腎結石の実データ700人では、全体の治癒率は経皮術Bの勝ち(78.0% vs 82.6%)なのに、結石の大きさで層別すると小も大も開腹手術Aの勝ち。調整化公式で計算すると P(治癒|do(A)) = 83.3% > P(治癒|do(B)) = 77.9% で、介入の正解はAです。
- 調整化公式 $P(y \mid do(x)) = \sum_z P(y \mid x, z)P(z)$ は切断された分解から1行で出ます。 $P(y \mid x)$ との違いは重みが $P(z \mid x)$ から $P(z)$ に変わるだけ。Day 11 のバックドア基準はこの式で正当化され、合成SCMの do サンプリングとも数値が一致します(+1.00 vs +1.00)。
今日の問い
「治療を受けた人の治癒率」と「治療を受けさせたときの治癒率」。言葉にするとほとんど同じに見えるのに、なぜ違う数字になるのか?
そして後者は、介入していない観察データから計算できるのか?
今日のゴールは、この2つを $P(Y \mid X)$ と $P(Y \mid do(X))$ という別の式に書き分け、「条件付け=絞り込み、do=書き換え」という対比を自分の言葉とコードにすることです。前半は700人の実データで、観察と介入の答えが実際に食い違う様子を見ます。後半は合成SCMで、do を文字どおりサンプリングとして実装し、調整化公式が本当に介入の答えを返すことを確かめます。
概念① 見ることとすること:P(Y|X) と P(Y|do(X))
$P(Y \mid X{=}1)$ は、観測されたデータの中から $X=1$ だった人を取り出して眺めた分布です。母集団を絞り込むだけで、世界には何もしていません。ここで大事なのは、誰が $X=1$ になったかを決めた機構がそのまま生きていることです。重症な患者ほど手術を受けやすいなら、「手術を受けた人」は「重症な人」の集まりでもあります。$P(Y \mid X{=}1)$ にはその偏りが焼き付いています。
$P(Y \mid do(X{=}1))$ は違います。母集団は絞り込まず、全員の $X$ を 1 に書き換えたときの $Y$ の分布です。構造因果モデル(SCM)の言葉で言うと、$X$ を決めていた構造方程式をモデルから削除して、定数 $X=1$ に置き換えます。DAG の上では $X$ に入ってくる矢印をすべて切る操作なので、graph surgery(グラフの外科手術)と呼ばれます。重症度 $Z$ の分布そのものは手つかずのまま、$Z$ が $X$ を決める経路だけが消えます。
実はこの分布、すでに知っています。RCT(Day 4)の無作為化とは、本人の意欲や重症度が割付を決める仕組みを物理的に壊して、コイン投げに置き換える行為でした。つまりRCTが測っているのは $P(Y \mid do(X))$ そのものです。do 演算子は、RCTが物理世界でやっている手術を、モデルの上で式として書けるようにした記法だといえます。
概念② 切断された分解と調整化公式:Day 11 の宿題を回収する
では、介入せずに $P(y \mid do(x))$ を計算できるのか。上の DAG($Z \to X$、$Z \to Y$、$X \to Y$)なら、同時分布は親子関係に沿って分解できます(Day 10 で見た因果マルコフ条件です)。
$$
P(z, x, y) = P(z), P(x \mid z), P(y \mid x, z)
$$
$do(X{=}x)$ は $X$ の構造方程式だけを消す操作でした。分解の言葉では、割付の因子 $P(x \mid z)$ だけを積から取り除くことに対応します。これを切断された分解(truncated factorization)と呼びます。
$$
P(z, y \mid do(x)) = P(z), P(y \mid x, z)
$$
あとは $z$ を周辺化するだけで、目的の量が出ます。
$$
P(y \mid do(x)) = \sum_z P(y \mid x, z), P(z)
$$
これが**調整化公式(adjustment formula)**です。条件付き分布と並べてみます。
$$
P(y \mid x) = \sum_z P(y \mid x, z), P(z \mid x)
$$
違いはたった1箇所、層 $z$ を混ぜ合わせる重みが $P(z \mid x)$ か $P(z)$ かだけです。条件付けは「$x$ を受けた人たちの症例構成」で平均し、do は「母集団全体の症例構成」で平均します。腎結石でいえば、開腹手術Aを受けた人の75%は大結石ですが、母集団全体では大結石は49%です。観察の数字はこの偏った重みのせいで歪み、do は全員を同じ土俵に載せ直します。「見ること」と「すること」の差の正体は、この重みの差です。
ここで Day 11 の宿題が回収できます。上の導出で使ったのは「$Z$ が $X$ と $Y$ の共通原因への道、つまりバックドアパスをすべて塞いでいる」という構造だけです。一般に、$Z$ がバックドア基準を満たす集合なら同じ式が成り立つことが証明できます(Pearl 1993)。バックドア基準とは「この集合で調整すれば $P(y \mid do(x))$ が観測量だけで書ける」ことのグラフ側の言い換えだったわけです。右辺に do が残っていないことに注目してください。すべて観察データから推定できる量です。因果の量を観測量だけの式に翻訳できることを**識別(identification)**と呼び、Day 13 以降もこの言葉が主役になります。ついでに言うと、Day 6 で層別の加重平均が真値に戻ったのも、あのとき計算していたのがまさにこの式だったからです。
手を動かす①:腎結石700人で逆転を再現する
実データで見ます。Charig et al. (1986) が報告した腎結石治療の古典的な数表で、開腹手術Aと経皮的手術Bを受けた各350人について、結石の大きさ別の治癒数が残っています。Simpson のパラドックスの定番教材です。
# 腎結石700人の数表(Charig et al. 1986): (治癒数, 症例数)
data = {
("A", "小"): (81, 87), ("A", "大"): (192, 263), # A = 開腹手術
("B", "小"): (234, 270), ("B", "大"): (55, 80), # B = 経皮的手術
}
# 条件付き確率 P(治癒 | 治療, 結石サイズ) と全体 P(治癒 | 治療)
for t in ("A", "B"):
(cs, ns), (cl, nl) = data[(t, "小")], data[(t, "大")]
print("%s: 小 %.1f%% 大 %.1f%% 全体 %.1f%%"
% (t, 100 * cs / ns, 100 * cl / nl, 100 * (cs + cl) / (ns + nl)))
# 調整化公式 P(治癒|do(治療)) = Σ_z P(治癒|治療, z) P(z)
pz = {z: (data[("A", z)][1] + data[("B", z)][1]) / 700 for z in ("小", "大")}
print("P(小)=%.2f, P(大)=%.2f" % (pz["小"], pz["大"]))
for t in ("A", "B"):
do_t = sum(data[(t, z)][0] / data[(t, z)][1] * pz[z] for z in ("小", "大"))
print("P(治癒|do(%s)) = %.3f" % (t, do_t))
A: 小 93.1% 大 73.0% 全体 78.0%
B: 小 86.7% 大 68.8% 全体 82.6%
P(小)=0.51, P(大)=0.49
P(治癒|do(A)) = 0.833
P(治癒|do(B)) = 0.779
全体ではBが 78.0% 対 82.6% で勝っています。ところが小結石だけ見るとAの勝ち(93.1% vs 86.7%)、大結石だけ見てもAの勝ち(73.0% vs 68.8%)です。部分ではすべてAが良いのに、合計するとBが良くなる。これが Simpson の逆転です。
種明かしは割付にあります。大結石343人のうち77%が開腹手術Aに回され、小結石357人のうち76%が体への負担が軽いBに回されていました。医師が難しい症例ほど実績のある開腹手術を選んだ、という自然な診療判断です。結果として、Aの全体成績78.0%は「難症例ばかり受け持ったチームの成績」であり、Bの82.6%は「楽な症例で稼いだ成績」になっています。結石の大きさ $Z$ が治療の選択と治癒の両方を動かす交絡因子で、全体の比較はこの偏りごと数字にしています。
「では治療を選ばされるとしたらどちらか」に答えるのが調整化公式です。両治療を同じ重み P(小)=0.51、P(大)=0.49 で平均し直すと、P(治癒|do(A)) = 83.3%、P(治癒|do(B)) = 77.9%。介入の答えはAの勝ちで、観察の答えと逆になりました。
手を動かす②:同じSCMから「観測」と「do」をサンプリングする
do は式の上の記号ではなく、実装できる操作です。腎結石と同じ構造($Z \to X$、$Z \to Y$、$X \to Y$)の SCM を合成データで組み、観測サンプリングと do サンプリングを両方走らせます。真の効果は $\tau = +1.0$、つまり処置は誰にとってもプラスに作っておきます。
import numpy as np
rng = np.random.default_rng(12)
N, TAU = 200_000, 1.0 # 真の効果 τ = +1.0
def f_y(x, z, rng):
"""Yの構造方程式: 重症度zがYを下げ、処置xがτだけ上げる"""
return 5.0 - 2.0 * z + TAU * x + rng.normal(0, 1, len(z))
# 観測サンプリング: 割付は重症度で決まる(重症ほど処置X=1を受けやすい)
z = rng.normal(0, 1, N) # 重症度(交絡因子)
x = rng.binomial(1, 1 / (1 + np.exp(-1.5 * z))) # 割付の構造方程式
y = f_y(x, z, rng)
# doサンプリング: 母集団(z)はそのまま、割付の式だけを定数に書き換える
y_do1 = f_y(np.ones(N), z, rng) # do(X=1): 全員を処置に
y_do0 = f_y(np.zeros(N), z, rng) # do(X=0): 全員を未処置に
print("観測: E[Y|X=1] - E[Y|X=0] = %+.2f" % (y[x == 1].mean() - y[x == 0].mean()))
print("介入: E[Y|do(X=1)] - E[Y|do(X=0)] = %+.2f (真値 +1.0)" % (y_do1.mean() - y_do0.mean()))
観測: E[Y|X=1] - E[Y|X=0] = -1.12
介入: E[Y|do(X=1)] - E[Y|do(X=0)] = +1.00 (真値 +1.0)
同じモデルから引いたのに、観測の差は -1.12 で「処置は害」に見え、介入の差は +1.00 で真値どおりです。コードを見ると、do サンプリングでやったことは1行の書き換えだけだと分かります。x を作っていた割付の式を捨てて np.ones(N) を渡した、それだけです。$Z$ の分布にも $Y$ の式にも触れていません。これが graph surgery の実装で、$Z \to X$ の矢印を切る操作がコード上では「引数の差し替え」として現れます。物々しい名前の操作がたった1つの引数の差し替えになる、この落差でようやく do 演算子が自分のものになった気がしました。観測側で差が負になるのは、処置群に重症者($Z$ が高く、もともと $Y$ が低い人)が集まるからです。腎結石で開腹手術Aが不利に見えたのと同じ機構を、今度は機構ごと自作して確かめたことになります。
手を動かす③:調整化公式の答え合わせ
仕上げに概念②の主張を数値で検証します。調整化公式の右辺は観測量だけでできているので、観測データ (z, x, y) のみから計算できるはずです。それが、モデルを書き換えないと作れないはずの do サンプリングの答えと一致するかを見ます。
import numpy as np
rng = np.random.default_rng(12)
N, TAU = 200_000, 1.0
def f_y(x, z, rng):
return 5.0 - 2.0 * z + TAU * x + rng.normal(0, 1, len(z))
# 手を動かす②と同じSCM・同じ乱数系列を再現する
z = rng.normal(0, 1, N)
x = rng.binomial(1, 1 / (1 + np.exp(-1.5 * z)))
y = f_y(x, z, rng)
y_do1, y_do0 = f_y(np.ones(N), z, rng), f_y(np.zeros(N), z, rng)
# 調整化公式 E[Y|do(x)] = Σ_z E[Y|X=x, Z=z] P(z) を層別(100分位ビン)で計算する
# 使うのは観測データ (z, x, y) だけ。doサンプルには一切触らない
edges = np.quantile(z, np.linspace(0, 1, 101))
edges[0], edges[-1] = -np.inf, np.inf
bins = np.digitize(z, edges[1:-1])
adj = {}
for xv in (0, 1):
adj[xv] = sum(y[(bins == b) & (x == xv)].mean() * (bins == b).mean()
for b in range(100))
print("調整化公式(観測データのみ): E[Y|do(1)]=%.2f, E[Y|do(0)]=%.2f, 差=%+.2f"
% (adj[1], adj[0], adj[1] - adj[0]))
print("doサンプリング(答え合わせ): E[Y|do(1)]=%.2f, E[Y|do(0)]=%.2f, 差=%+.2f"
% (y_do1.mean(), y_do0.mean(), y_do1.mean() - y_do0.mean()))
調整化公式(観測データのみ): E[Y|do(1)]=6.00, E[Y|do(0)]=5.00, 差=+1.00
doサンプリング(答え合わせ): E[Y|do(1)]=6.00, E[Y|do(0)]=5.00, 差=+1.00
一致しました。$Z$ の層ごとに $E[Y \mid X{=}x, Z{=}z]$ を取り、母集団全体の重み $P(z)$ で混ぜただけの計算が、水準も差も do サンプリングと同じ +1.00 を返しています。私たちは一度も介入していないのに、介入の答えを手に入れました。ただし無条件ではありません。この計算が正しいのは「$Z$ がバックドアを全部塞いでいる」という DAG の仮定を信じたときだけです。仮定の置き方(DAG)は Day 10、集合の選び方は Day 11、その正当化が今日、という3日がかりの分業でした。
つまづき・誤解しやすい点
-
「do はデータの操作ではなく、モデルの操作です。」 データフレームから
df[df.x == 1]と絞り込む操作は、どこまで工夫しても条件付けです。do は割付の構造方程式を書き換える操作で、観測済みデータの中には存在しません。だからこそ「do をどう観測量に翻訳するか」という識別の問題が立ち、その答えのひとつが今日の調整化公式でした。翻訳できない場合の第2の道は次回Day 13で扱います。 - 「層別すれば常に正しい、わけではありません。」 今日は層別(調整)側が正解でしたが、それは $Z$(結石の大きさ)が処置より前に決まる交絡因子だったからです。もし $Z$ が処置の後に決まる中間変数や合流点なら、同じ計算がバイアスを作ります(Day 11 の悪い調整)。Simpson の逆転のどちらの数字を信じるかはデータからは決められず、DAG が決めます。DAG は仮定の宣言だ、という Day 10 の話がここで効いてきます。
- 「予測だけなら P(Y|X) で困らないことも多いです。」 「手術を受けた患者の治癒率は78%」は、記述としては何も間違っていません。事故が起きるのは、この数字を「手術させたときの治癒率」として意思決定に使った瞬間です。predict と intervene のどちらの問いに答えたいのかを先に言葉にするのが、今日の対比のいちばん実務的な使い方です。
GISデータ実務での使い方
- 雨予測がうまいモデルは、雨を降らせません。 予測モデルは $P(Y \mid X)$ の世界の住人で、特徴量とアウトカムの相関をどれだけ深く学習しても、それは「観測したときの答え」です。気圧計の針を手で動かしても天気は変わらないように、予測に効く特徴量をこちらの都合で動かしても、アウトカムが同じように動く保証はどこにもありません。モデルの精度報告と施策提案の間には、$P(Y \mid X)$ から $P(Y \mid do(X))$ への飛躍が挟まっていると意識しておきたいところです。
- 管路更新の優先度スコアと更新効果は別物です。 漏水リスクの予測モデルが返すのは P(漏水 | 管路属性)、つまり「放置したらどの管が危ないか」です。一方で更新計画の費用対効果に必要なのは P(漏水 | do(更新))、「更新したら漏水がどれだけ減るか」です。高スコアの管を更新しても、そのスコアを押し上げた要因(地盤や水圧)は更新後も残るかもしれません。優先度づけ(予測)と効果検証(介入)は estimand が違う、と切り分けて話すだけで議論がかみ合います。
- 「導入自治体では成果が出ています」は P(Y|X) です。 施策を導入した自治体は、財政や体制の面で導入できる条件がそろっていた自治体でもあります。腎結石で開腹手術が難症例を受け持っていたのとちょうど逆向きの選択です。営業資料や事例集の数字を見たら、まず「割付を決めたのは誰か」を思い出し、規模や老朽化率で層別した数字を求めるのが今日の教訓の使いどころです。
参考(本棚)
- Pearl, J. (2009) Causal Inference in Statistics: An Overview. Statistics Surveys 3, 96-146。do演算子・切断された分解・調整化公式・バックドア基準の正当化(§3.2〜3.3)。今日の下敷き
- Pearl, Glymour & Jewell『入門 統計的因果推論』(朝倉書店)。介入とバックドア調整の教科書的展開。Simpson のパラドックスの因果的解決も第1章から一貫して扱う
- Charig et al. (1986) Comparison of treatment of renal calculi by open surgery, percutaneous nephrolithotomy, and extracorporeal shockwave lithotripsy. BMJ 292, 879-882。腎結石データの原典。本文の700人の数表はこの論文の集計値(患者単位の生データではありません)
- 後半の図は合成データです。真の効果 $\tau$ を握った SCM から観測と do の両方をサンプリングすることで、条件付けとの乖離と調整化公式の一致を再現しました。数値はすべて本文コードとシード固定で再現できます
次回予告(Day 13)
調整化公式が使えるのは、バックドアを塞ぐ変数集合が観測できているときだけです。では交絡因子がどうしても測れないとき、観察データからの因果推論は全滅なのでしょうか。実は、処置と結果の間に立つ媒介変数がうまく観測できていれば、バックドアが塞げなくても因果効果を点で識別できる状況があります。喫煙→タール→肺がんの古典例で知られる第2の道、フロントドア基準です。次回は、この式を「2回のバックドア調整の合成」として1行ずつ導出し、未観測交絡入りの合成SCMでフロントドアだけが真値を当てることをオラクル比較で確かめます。

