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 7/全30回 傾向スコアとIPW、高次元の交絡を1次元に潰す

0
Posted at

この連載について

因果推論を「本を読んだ」で終わらせず、自分の言葉で説明でき、コードで再現できる状態まで落とす30日連載です。直前のDay 6では、観察データでも「同じ$X$の層の中では割付がランダムに見える」という非交絡性の土台を置きました。今日はその続き、共変量が多くて層を切れないときの道具です。ここからは合成データではなく実データを使います。疫学の定番教材データ nhefs(禁煙と体重変化)です。


TL;DR(3行)

  • 傾向スコア $e(x) = \Pr(W=1 \mid X=x)$ は高次元の$X$を1次元に潰す要約。 同じ$e(x)$の中では共変量の分布が両群でそろう(バランシングスコア)ので、9個の共変量の代わりに$e(x)$1本で条件付けできます。nhefsの素朴な差 2.54 kg は、IPWで 3.42 kg に補正されます(今回は過小評価の向き)。
  • IPW=逆確率重み付けは、各人を $1/e(x)$(対照は $1/(1-e(x))$)で膨らませて「全員が処置を受けた世界」と「全員が受けなかった世界」の擬似母集団を作ります。 重み付け後は年齢のSMDが 0.28 → 0.01 に落ち、9共変量すべてが |SMD| < 0.03 にそろいます。
  • 重みは裾が命。 最大重み 16.7(1人で16.7人分)は安定化重みで 4.3 に縮み、上位1%をトリミングしても推定は 3.44 → 3.46 kg とほぼ動きません。トリミングで結論が大きく動くなら、それはオーバーラップ(Day 6)が薄い信号です。

今日の問い

Day 6 で「同じ$X$の層の中では割付がランダムに見える」ことを確かめた。でも現実の共変量は1個ではない。年齢、性別、人種、学歴、喫煙本数、喫煙年数…と9個あったら、層の数は掛け算で爆発し、処置群か対照群が空っぽの層だらけになる。
高次元の交絡を、層を切らずに調整する方法はあるか?

答えの中心が傾向スコアです。今日はまず「なぜ1次元の要約で足りるのか」を言葉にして、次に実データで素朴な差、傾向スコアの推定、手組みのIPW、安定化重み、重みのトリミングまでを一気にコード化します。


概念① 傾向スコア:高次元のXを1次元に潰すバランシングスコア

傾向スコアは、共変量$X$を持つ人が処置を受ける条件付き確率です。

$$
e(x) = \Pr(W_i = 1 \mid X_i = x)
$$

Rosenbaum & Rubin (1983) が示した性質が2つあります。第一に、傾向スコアで条件付けると共変量と割付が独立になります($X_i \perp!!!\perp W_i \mid e(X_i)$)。同じ$e(x)$の人たちを集めた層の中では、処置群と対照群で$X$の分布がそろうのです。この性質を持つ量をバランシングスコアと呼びます。第二に、Day 6 の非交絡性が$X$について成り立つなら、$e(x)$について条件付けても成り立ちます。

$$
\bigl(Y_i(0),, Y_i(1)\bigr) \perp!!!\perp W_i ,\bigm|, e(X_i)
$$

なぜ1次元で足りるのでしょうか。同じ$e(x)$を持つ2人は、$X$の中身が違っても「処置に回る確率が同じだった」人たちです。片方が処置を受け、片方が受けなかったのは、確率が同じだった以上、偶然の振り分けです。つまり Day 6 の「層の中はRCT」の層を、$X$の全組み合わせではなく$e(x)$の値で切ってよいわけです。9次元の条件付けが1次元に潰れます。

実務では$e(x)$は未知なので、ロジスティック回帰などで推定します。ここで大事な視点を Imbens & Rubin 第13章から持ってきます。傾向スコアの推定で目指すのは割付の予測精度ではなくバランスです。推定した$e(x)$で条件付けたときに共変量がそろうこと、それだけが役目です。この違いは後の「つまづき」で効いてきます。

概念② IPW:重みを付けると、なぜ交絡が消えるのか

**IPW(逆確率重み付け、inverse probability weighting)**は、傾向スコアの逆数を重みにする推定法です。Horvitz-Thompson 型の式はこう書けます。

$$
\hat{E}[Y(1)] = \frac{1}{N}\sum_{i=1}^{N} \frac{W_i, Y_i}{e(X_i)}, \qquad
\hat{E}[Y(0)] = \frac{1}{N}\sum_{i=1}^{N} \frac{(1-W_i), Y_i}{1-e(X_i)}
$$

重みを付けるとなぜ交絡が消えるのか。ここが今日いちばんの難所なので、1段落かけます。$e(x) = 0.1$ なのに処置を受けた人を考えてください。同じ$X$を持つ人が母集団に10人いたら、処置群に現れるのはそのうち1人だけです。残りの9人は対照群に行ってしまい、処置群からは見えません。だからその1人を $1/0.1 = 10$ 人分に膨らませると、見えなくなった9人の代役を務めてもらえます。全員にこの膨らましをすると、処置群は「処置を受けやすい人に偏ったサンプル」から「母集団の全員が処置を受けた世界の縮図」に化けます。対照群も $1/(1-e(x))$ で同じことをします。こうして作った2つの世界を擬似母集団と呼びます。擬似母集団の中では、処置を受けたかどうかと$X$の結びつきが切れているので、あとは平均を引き算するだけで因果効果になります。交絡を「消した」のではなく、交絡が最初から無い世界をデータの重み替えで作り直した、というのが正確な言い方です。

ただし、この作り直しができるのは$X$に入れた交絡だけです。$X$の外にある交絡は、どんな重みを付けても1ミリも動きません。Day 6 の「非交絡性は検証できない」はここでも一切割引なしで効いています。


手を動かす①:素朴な差と、禁煙する人の偏り

データは nhefs です。Hernán & Robins『Causal Inference: What If』の教材データで、米国の健康調査 NHANES I の追跡調査から作られています。1971年時点の喫煙者約1,600人を追跡し、1982年までに禁煙したか( qsmk )と、11年間の体重変化 kg( wt82_71 )が記録されています。問いは「禁煙は体重を何kg増やすか」です。なお1982年の追跡で体重を測れなかった63人は wt82_71 が欠測のため除外します。本当はこの欠測自体も選択バイアスとして重みで扱える話ですが、今日は交絡に絞ります。

import numpy as np
import pandas as pd
from causaldata import nhefs

d = nhefs.load_pandas().data
d = d.dropna(subset=["wt82_71"]).copy()   # 追跡不能でアウトカム欠測の63人を除外

n1, n0 = (d.qsmk == 1).sum(), (d.qsmk == 0).sum()
naive = d.loc[d.qsmk == 1, "wt82_71"].mean() - d.loc[d.qsmk == 0, "wt82_71"].mean()
print("禁煙した %d 人 / 吸い続けた %d 人" % (n1, n0))
print("素朴な差 = %.2f kg" % naive)
cols = ["age", "smokeintensity", "wt71"]
print(d.groupby("qsmk")[cols].mean().round(1))
禁煙した 403 人 / 吸い続けた 1163 人
素朴な差 = 2.54 kg
       age  smokeintensity  wt71
qsmk                            
0.0   42.8            21.2  70.3
1.0   46.2            18.6  72.4

素朴な引き算は 2.54 kg です。でも下の表を見ると、禁煙した人( qsmk = 1)は吸い続けた人より平均で3.4歳年上で、1日の喫煙本数が少なく、元の体重も重めです。禁煙は無作為に起きていません。そして年齢は体重変化と負の相関(このデータで $-0.23$)を持ちます。年を取るほど11年間で体重は増えにくいのです。禁煙群には「増えにくい年長者」が多く混ざっているので、素朴な差は真の効果を過小評価する向きに歪みます。Day 6 の例は過大評価でしたが、交絡の向きは構造次第で、どちらにも転びます。


手を動かす②:ロジスティック回帰でe(x)、そして手組みIPW

傾向スコアをロジスティック回帰で推定します。仕様は Hernán & Robins の教材と同じで、性別・人種・年齢・学歴・喫煙本数・喫煙年数・運動習慣・活動度・元の体重の9共変量、連続量には2乗項を入れます。ひとつ実務の注意です。Day 6 の予告で sklearn と書きましたが、sklearn の LogisticRegression は既定でL2正則化がかかり、$e(x)$が0.5側に縮んで重みが歪みます。今日は素の最尤推定になる statsmodels を使います(sklearn なら penalty=None を指定してください)。

import statsmodels.formula.api as smf

ps_model = smf.logit(
    "qsmk ~ sex + race + age + I(age**2) + C(education)"
    " + smokeintensity + I(smokeintensity**2) + smokeyrs + I(smokeyrs**2)"
    " + C(exercise) + C(active) + wt71 + I(wt71**2)",
    data=d).fit(disp=0)
d["ps"] = ps_model.predict(d)

w = d["qsmk"].to_numpy()
y = d["wt82_71"].to_numpy()
e = d["ps"].to_numpy()
N = len(d)

# Horvitz-Thompson 式をそのまま書く
ey1 = np.sum(w * y / e) / N
ey0 = np.sum((1 - w) * y / (1 - e)) / N
print("E[Y(1)] = %.2f kg, E[Y(0)] = %.2f kg" % (ey1, ey0))
print("IPWのATE = %.2f kg(素朴な差は 2.54)" % (ey1 - ey0))
print("傾向スコアの範囲: %.3f から %.3f" % (e.min(), e.max()))
E[Y(1)] = 5.20 kg, E[Y(0)] = 1.78 kg
IPWのATE = 3.42 kg(素朴な差は 2.54)
傾向スコアの範囲: 0.051 から 0.777

概念②の式をそのまま2行書くだけで、推定値は 2.54 から 3.42 kg に動きました。「全員が禁煙した世界」の平均体重変化が +5.20 kg、「全員が吸い続けた世界」が +1.78 kg、差が禁煙の効果です。教科書(Hernán & Robins 第12章)の答え約3.4 kgと一致します。下の図の左が傾向スコアの分布です。禁煙群がやや高スコア側に偏っていて、これが9共変量の交絡の1次元要約です。右が重み付け後で、両群の分布がそろい、どちらの擬似母集団もほぼ元の1,566人分になっています。

Fig 1: 傾向スコアの分布(左:生データでは禁煙群が高スコア側に偏る、右:逆確率重み付け後は両群がそろい母集団全体の縮図になる)

次に安定化重みです。素の重み $1/e(x)$ は、処置の周辺確率 $\Pr(W=1)$ を分子に乗せた $\Pr(W=1)/e(x)$ に置き換えられます(対照は $\Pr(W=0)/(1-e(x))$)。分母が小さい人は分子も小さくなるので、重みの暴れが抑えられます。

p = w.mean()                                       # 処置の周辺確率 Pr(W=1)
hw = np.where(w == 1, 1 / e, 1 / (1 - e))          # 素の逆確率重み
sw = np.where(w == 1, p / e, (1 - p) / (1 - e))    # 安定化重み

sy1 = np.sum(w * sw * y) / np.sum(w * sw)
sy0 = np.sum((1 - w) * sw * y) / np.sum((1 - w) * sw)
print("素の重み  : 平均 %.2f, 最大 %.1f, 合計 %.0f (N=%d)"
      % (hw.mean(), hw.max(), hw.sum(), N))
print("安定化重み: 平均 %.2f, 最大 %.1f, 合計 %.0f" % (sw.mean(), sw.max(), sw.sum()))
print("安定化IPWのATE = %.2f kg" % (sy1 - sy0))
素の重み  : 平均 2.00, 最大 16.7, 合計 3126 (N=1566)
安定化重み: 平均 1.00, 最大 4.3, 合計 1564
安定化IPWのATE = 3.44 kg

素の重みの合計 3126 は $2N$ とほぼ一致します。処置側と対照側にそれぞれ$N$人分の擬似母集団を作っているからです。安定化重みは合計が$N$に戻り、最大重みが 16.7 から 4.3 に縮みます。推定値は 3.44 kg で実質同じです(重みの合計で割る自己正規化の分だけ、素のHT式の3.42と小数点以下がずれます)。安定化は答えを変えずに計算を安定させる保険です。そして重み付けが仕事をした証拠がバランスです。重み付け前は年齢のSMD(標準化平均差)が 0.28 ありましたが、重み付け後は 0.01。9共変量すべてが |SMD| < 0.03 に収まります。

Fig 2: 重み付け前後の共変量バランス(年齢のSMD 0.28が重み付け後0.01へ、9共変量すべて0.03未満にそろう)


手を動かす③:重みの裾を見る、オーバーラップの実務版

Day 6 で「オーバーラップが崩れた領域は推定できない」ことを見ました。IPWではこの問題が極端な重みという形で顔を出します。$e(x)$が0や1に近い人の重みは爆発し、その数人が推定値を支配します。だからIPWを使ったら重みの分布、特に裾を必ず見ます。定番の処方がトリミングで、極端な重みを閾値で頭打ちにします(傾向スコア側を0.01から0.99などで切る流儀もあります)。

def wmean_diff(wt):
    m1 = np.sum(w * wt * y) / np.sum(w * wt)
    m0 = np.sum((1 - w) * wt * y) / np.sum((1 - w) * wt)
    return m1 - m0

cap = np.quantile(hw, 0.99)              # 素の重みの99%点で頭打ちにする
print("トリミング閾値(99%%点)= %.1f, 頭打ちになる人数 = %d" % (cap, (hw > cap).sum()))
print("トリミング前 = %.2f kg / 後 = %.2f kg"
      % (wmean_diff(hw), wmean_diff(np.minimum(hw, cap))))
print("e(x) が 0.01 未満または 0.99 超の人数 = %d" % (((e < 0.01) | (e > 0.99)).sum()))
トリミング閾値(99%点)= 7.4, 頭打ちになる人数 = 16
トリミング前 = 3.44 kg / 後 = 3.46 kg
e(x) が 0.01 未満または 0.99 超の人数 = 0

nhefs はオーバーラップが良好です。傾向スコアは 0.051 から 0.777 に収まり、0.01未満や0.99超の人はゼロ。最大重み16.7の人を含む上位1%(16人)を頭打ちにしても、推定値は 3.44 から 3.46 kg へ 0.02 kg しか動きません。この「動かなさ」自体が診断です。逆にトリミングで結論が大きく動いたら、推定値が裾の数人に支配されている、つまり比較相手がほとんどいない領域に頼っている赤信号です。なおトリミングには代償もあります。重みを切った推定はもう厳密には母集団全体のATEではなく、「比較できる人たち」の効果に対象がすり替わります。切ったこと、切った後の対象が誰なのかは、明示して報告します。

Fig 3: 重みの分布と裾(左:素の逆確率重みは最大16.7まで伸び、赤の99%点7.4でトリミング、右:安定化重みは最大4.3に縮む)


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

  • 傾向スコアは「当てるほど良い」ではありません。 処置だけをよく説明してアウトカムに効かない変数を足すと、交絡は1つも減らないのに重みの分散だけ増えます。極端な話、割付を完璧に予測できたら$e(x)$は0か1に貼りつき、比較相手が消えます(Day 6 のオーバーラップ崩壊)。目指すのは予測精度ではなく、Fig 2 で確認したバランスです。
  • IPWが処理できるのは、$X$に入れた交絡だけです。 SMDが全部ゼロになっても、それは観測した共変量がそろった証拠にすぎません。測っていない交絡(禁煙の例なら「健康への意識の変化」など)は重みでは動かず、バランス表からも見えません。どれだけ強い隠れ交絡なら結論が覆るかは、Day 9 の感度分析で定量化します。
  • 推定した傾向スコアの誤差は、最終の標準誤差に響きます。 今日の手組みは点推定に絞りましたが、$e(x)$を推定してから重みに使う2段構えなので、素朴に重み付き平均の分散を計算すると不確かさを見誤ります。実務ではブートストラップか専用の頑健標準誤差を使います。機械学習で$e(x)$を推定する場合の正しい扱いは Day 24(Double ML)で戻ってきます。

実務での使いどころ(空間データ/天地人の文脈)

  • 自治体の施策評価は Day 6 の続きをそのまま実装できます。 衛星データ解析を導入した自治体は、人口・財政力・管路の老朽化率・データ整備度の点で「導入しやすい自治体」に偏っています。人口や財政などから導入確率$e(x)$を推定して重みでならせば、共変量が9個でも層を切らずに比較できます。
  • 最初に描くのは Fig 1 です。 傾向スコアの分布図は「どの自治体に比較相手がいるか」の地図になります。大都市が全部高スコア側、小規模町村が全部低スコア側に貼りついていたら、その端の効果はデータからは語れません。
  • 重みの爆発は分析の失敗ではなく情報です。 1自治体で数十自治体分の重みが出たら、そこは比較不能な領域を無理に埋めようとしている場所です。トリミングして「効果を語れるのはこの範囲の自治体」と正直に線を引く方が、結論は弱くなっても信頼されます。
  • 隣の自治体に効果が漏れる状況は、IPWの守備範囲外です。 重みは各自治体を独立に膨らませるだけなので、広域連携のような波及(SUTVA違反)は扱えません。空間の干渉は Day 28 でまとめて掘ります。

参考(本棚)

  • Imbens & Rubin『Causal Inference for Statistics, Social, and Biomedical Sciences』第13章 傾向スコアの推定。「目標は予測精度ではなくバランス」という仕様選びの規律で、今日の下敷き
  • Matheus Facure『Causal Inference in Python』第5章。IPW・安定化重み・擬似母集団・ポジティビティの実装視点で、今日のもう1枚の下敷き
  • Hernán & Robins『Causal Inference: What If』第12章。nhefsデータの本家で、IPWで約3.4 kgという答え合わせもこの章と一致
  • データ: nhefs( causaldata パッケージ経由、MITライセンス、取得日2026-08-15)。米国 NHANES I の追跡調査 NHEFS から Hernán & Robins が教材用に公開した抽出データ

次回予告(Day 8)

今日は傾向スコアで「重みを付けてならす」方向に進みました。同じ$e(x)$のもう1つの使い道が、処置を受けた1人ひとりに「そっくりな未処置の人」を探してくるマッチングです。次回は因果推論の歴史で有名なLaLondeの職業訓練データを使います。このデータには無作為化実験の答えという金標準があり、観察データの対照群とマッチングでどこまで実験の答えに迫れるかを検証できます。今日のSMDとラブプロットが、そのままバランス診断の道具として再登場します。

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?