背景

実務における因果推論において、処置変数を二値化して傾向スコアを使用したマッチングや重み付けを行うことは一般的である
しかし、プロモーションの実施回数のように処置変数が連続値をとるケースも多く、処置の回数に応じた効果を求めたい場合がある
そこで「一般化傾向スコア（Generalized Propensity Score, GPS）」を用いると、処置が連続変数の場合でも傾向スコアの考え方を用いた因果推論が可能になる
- 例えば、用量反応関数（Dose-Response Function, DRF）を推定することで、連続処置変数の非線形な効果を確認できる
本記事では、LaLonde データセットを用いて、教育年数という「連続値の処置変数」が所得に与える非線形な効果を GPS を用いて推定する方法を紹介する

データ

LaLonde（1986）の NSW（National Supported Work）プログラムデータのうち、CPS 比較グループ（sample=2）かつ nodegree=1（高卒資格なし）のサンプルを使用する
データの詳細は過去記事で紹介済み
所得（re78, re75）は右に歪んだ分布であり、右歪みの緩和と外れ値の影響を抑えるために log1p 変換を施す。処置変数 educ は 0〜11 の整数値であり変換不要のため元のスケールのまま使用する
nodegree=1 に絞ることで教育年数の範囲は 0〜11 年となり、高卒資格を持たない層における教育効果に焦点を当てる。最終的なサンプル数は 4,731 件である。

結果

Hirano & Imbens（2004）を参考に、下記ステップで進める

GPSの推定
GPSの妥当性確認
結果変数予測モデルの構築
用量反応関数の推定

GPSの推定

Hirano & Imbens（2004）に従い、共変量 $X_i$ を所与としたときの処置変数 $T_i$（教育年数）の条件付き分布を正規分布で近似する。具体的には、以下の線形回帰モデルを OLS で推定する。

$$T_i = \alpha_0 + \alpha_1 \text{age}_i + \alpha_2 \text{black}_i + \alpha_3 \text{hisp}_i + \alpha_4 \text{married}_i + \alpha_5 \text{re75_log1p}_i + \varepsilon_i$$

残差の標準偏差 $\hat{\sigma}$ を用いて、サンプル $i$ の GPS は以下のように定義される。

$$R_i = r(T_i, X_i) = \frac{1}{\hat{\sigma}} \phi\left(\frac{T_i - \hat{T}_i}{\hat{\sigma}}\right)$$

ここで $\hat{T}_i$ は $X_i$ から予測した処置変数の予測値、$\phi(\cdot)$ は標準正規分布の密度関数である。

また、処置水準 $t$ を任意の値に固定したときの GPS は $r(t, X_i) = \phi\left(\frac{t - \hat{T}_i}{\hat{\sigma}}\right) / \hat{\sigma}$ と書け、全サンプルの共変量を固定したまま $t$ の値を変化させることで、処置水準ごとのスコアを得ることができる。

GPS 推定モデル（OLS）の結果は以下の通りである。

変数	係数	t値
切片	10.798	114.0
age	−0.057	−21.2
black	−0.638	−6.95
hisp	−1.702	−20.1
married	−0.042	−0.55
re75_log1p	0.033	3.40

決定係数は $R^2 = 0.186$ であった。age（$t = -21.2$）と hisp（$t = -20.1$）の影響が大きく、年齢が高いほど・ヒスパニック系であるほど教育年数が低い傾向がある。

GPSの可視化

下図は、2 人のサンプル（教育年数が低め：25%点付近、高め：75%点付近）について、$r(t, X_i)$ を $t$（教育年数）に対してプロットしたものである。その個人の共変量 $X_i$ を固定したまま $t$ を変化させることで、各処置水準における「この個人が観測される確率密度」を示す。

GPSの妥当性確認

Hirano & Imbens（2004）の手順に従い、バランシングを検証する。まず処置変数を $K$ 区間に分割し、各区間に属するサンプルと区間外サンプルの共変量平均差を $t$ 統計量で測る（未調整）。

GPS 調整版は以下の手順で求める。

区間中央値 $t_{\text{med}}$ を固定し、全サンプルについて $r(t_{\text{med}}, X_i)$ を計算する
区間内サンプルの $r(t_{\text{med}}, X_i)$ の分位点を境界として $B$ ブロックを作り、全サンプルをそのブロックに割り当てる
各ブロック内で区間内・区間外の共変量平均差を $t$ 統計量で測り、ブロックのサンプル数で加重平均する

この操作を全区間について繰り返し、調整前後の $t$ 統計量を比較することで GPS による交絡除去の効果を評価する。

	age	black	hisp	married	re75_log1p
区間1（未調整）	21.50	2.12	9.09	14.39	5.37
区間2（未調整）	−21.50	−2.12	−9.09	−14.39	−5.37
区間1（GPS調整）	10.78	−1.34	5.97	6.73	0.92
区間2（GPS調整）	−5.30	1.41	−0.83	−3.86	−1.31

hisp の t 統計量は区間2では −9.09（未調整）から −0.83（GPS調整）、re75_log1p は −5.37 から −1.31 へと縮小しており、GPS による交絡除去の効果が確認できる。

結果変数予測モデルの構築

処置変数 $T$ と GPS $R$ の 2 次多項式で結果変数を回帰する。

$$\beta(t, r) = E[Y \mid T=t, R=r] = \gamma_0 + \gamma_1 t + \gamma_2 r + \gamma_3 t^2 + \gamma_4 r^2 + \gamma_5 t \cdot r$$

OLS の推定結果

変数	係数	t値
切片	7.084	18.02
educ	0.332	1.97
gps	−1.726	−0.27
educ²	−0.031	−2.44
gps²	−32.030	−1.60
educ × gps	1.264	2.26

用量反応関数の推定

任意の処置水準 $t$ に対する用量反応関数 $\mu(t)$ を以下の手順で推定する。

教育年数$t$ のグリッドを作成し、全サンプルの $r(t, X_i)$ を計算する
$\beta(t, r(t, X_i))$ を上記モデルで予測し、全サンプルで平均する

$$\mu(t) = \frac{1}{n} \sum_{i=1}^{n} \beta(t, r(t, X_i))$$

$\mu(t)$ の計算には、観測された処置水準での GPS（$r(T_i, X_i)$）ではなく、グリッド点 $t$ での GPS（$r(t, X_i)$）を用いて $\beta$ を評価し、$n$ サンプルで平均する点に注意する。

下記の用量反応関数が得られた。

考察

分析結果について

教育年数 3 年以下での効果が大きい
- 高卒資格を持たない層（nodegree=1）において、初等教育の有無が所得に大きな影響を与えるという直感と一致する
- ただ、教育年数 3 年以下の観測数は少ないため、分散の大きい推定結果と考えられる
教育年数 9 年以上では効果が頭打ちになる
- 教育年数 9 年を超えてもそれ以上の所得向上は望みにくいことが示唆される

手法について

GPS を用いた用量反応関数の推定によって、処置水準ごとの効果の変化を可視化でき、連続変数である処置をどの程度まで実施すればよいか、の検討をつけることができる

本分析の限界

処置変数の分布が偏っているため、教育年数 3 年以下や 9 年以上の領域では観測数が少なく、推定の不確実性が高い
結果変数予測モデルは 2 次多項式を仮定しており、真の関係が異なる関数形であった場合に推定が歪む可能性がある
結果変数モデルの精度評価を行っていない

参考文献

Hirano, K., & Imbens, G. W. (2004). The propensity score with continuous treatments. Applied Bayesian Modeling and Causal Inference from Incomplete-Data Perspectives, 226164, 73-84. https://onlinelibrary.wiley.com/doi/10.1002/0470090456.ch7
登藤直弥, 小林哲郎, 稲増一憲. (2016). ソフトニュースへの接触は政治的関心を高めるか―一般化傾向スコアを用いた因果推論―. 行動計量学, 43(2), 129-141. https://doi.org/10.2333/jbhmk.43.129

コード

環境準備

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.stats import norm
from scipy import stats

データ取得・前処理

Y_COL = 're78_log1p'
T_COL = 'educ'
X_COLS = ['age', 'black', 'hisp', 'married', 're75_log1p']
LOG1P_EDUC = False

df_org = pd.read_csv('../data/input/ec675_nsw.tab', sep='\t')
df_all = df_org.copy().query('sample==2 and nodegree==1')
df_all['re75_log1p'] = np.log1p(df_all['re75'])
df_all['re78_log1p'] = np.log1p(df_all['re78'])
if LOG1P_EDUC:
    df_all['educ'] = np.log1p(df_all['educ'])

GPS 推定

X_gps = sm.add_constant(df_all[X_COLS])
model_gps = sm.OLS(df_all[T_COL], X_gps).fit()

t_hat = model_gps.fittedvalues
sigma = np.sqrt(model_gps.mse_resid)
df_all['gps'] = norm.pdf(df_all[T_COL], loc=t_hat, scale=sigma)

バランシング検証

N_TREAT_INTERVALS = 2
N_GPS_BLOCKS = 5

_, bin_edges = pd.qcut(df_all[T_COL], q=N_TREAT_INTERVALS, retbins=True, duplicates='drop')
interval_labels = [str(i + 1) for i in range(len(bin_edges) - 1)]
df_all['treat_interval'] = pd.cut(df_all[T_COL], bins=bin_edges, labels=interval_labels, include_lowest=True)
interval_medians = df_all.groupby('treat_interval', observed=True)[T_COL].median()

def t_stat_unadjusted(df, col, interval_label):
    g1 = df.loc[df['treat_interval'] == interval_label, col]
    g0 = df.loc[df['treat_interval'] != interval_label, col]
    return stats.ttest_ind(g1, g0, equal_var=False).statistic

def t_stat_adjusted(df, col, interval_label, t_med, n_blocks):
    r_t_med = norm.pdf(t_med, loc=t_hat, scale=sigma)
    df = df.copy()
    df['r_t_med'] = r_t_med
    in_interval = df['treat_interval'] == interval_label
    block_bins = df.loc[in_interval, 'r_t_med'].quantile(np.linspace(0, 1, n_blocks + 1)).values.copy()
    block_bins[0] -= 1e-10
    df['gps_block'] = pd.cut(df['r_t_med'], bins=block_bins, labels=False)
    diffs, ses, ns = [], [], []
    for b in range(n_blocks):
        mask_block = df['gps_block'] == b
        g1 = df.loc[mask_block & in_interval, col]
        g0 = df.loc[mask_block & ~in_interval, col]
        if len(g1) < 2 or len(g0) < 2:
            continue
        diff = g1.mean() - g0.mean()
        se = np.sqrt(g1.var() / len(g1) + g0.var() / len(g0))
        n = len(g1) + len(g0)
        diffs.append(diff * n)
        ses.append(se ** 2 * n ** 2)
        ns.append(n)
    if not ns:
        return np.nan
    N = sum(ns)
    return (sum(diffs) / N) / (np.sqrt(sum(ses)) / N)

結果変数予測モデル


df_fit = df_all.dropna(subset=['gps']).copy()
t, r = df_fit[T_COL], df_fit['gps']

X_outcome = sm.add_constant(pd.DataFrame({
    'educ': t, 'gps': r, 'educ2': t ** 2, 'gps2': r ** 2, 'educ_gps': t * r,
}))
model_outcome = sm.OLS(df_fit[Y_COL], X_outcome).fit()

用量反応関数の推定・可視化

t_orig = (lambda x: np.expm1(x) if LOG1P_EDUC else x)
t_grid = np.linspace(df_all[T_COL].min(), df_all[T_COL].max(), 50)

mu = []
for t_val in t_grid:
    r_t = norm.pdf(t_val, loc=t_hat, scale=sigma)
    X_pred = pd.DataFrame({
        'const': 1, 'educ': t_val, 'gps': r_t,
        'educ2': t_val ** 2, 'gps2': r_t ** 2, 'educ_gps': t_val * r_t,
    })
    mu.append(model_outcome.predict(X_pred).mean())

df_drf = pd.DataFrame({'educ': t_grid, 'mu': mu})

fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(t_orig(df_drf['educ']), df_drf['mu'], linewidth=2)
ax.set_xlabel('教育年数 (educ)')
ax.set_ylabel('log1p(re78) の期待値 µ(t)')
ax.set_title('用量反応関数（Dose-Response Function）')
plt.tight_layout()
plt.show()