「回帰分析から学ぶ計量経済学」をPythonで写経 ～ 第７章「反実仮想の世界」⑤傾向スコアマッチング、ランダム化比較実験

第７章「反実仮想の世界」

書籍の著者　山澤成康 先生

この記事は、書籍「回帰分析から学ぶ計量経済学」第７章「反実仮想の世界」の Python写経活動 を取り扱います。

第７章は反実仮想、因果推論に取り組みます。
今回のテーマは傾向スコアマッチングとランダム化比較実験（RCT）です。
そして今回は「回帰分析から学ぶ計量経済学」寄り道写経の最終回です！

では書籍を開いて回帰分析の旅に出発です🚀

はじめに

書籍「回帰分析から学ぶ計量経済学」のご紹介

このシリーズは書籍「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」（オーム社、「テキスト」と呼びます）の Python 写経です。

テキストは、2023年11月に発売され、副題「Excelで読み解く経済のしくみ」のとおり、主に Excel を用いて、計量経済学を平易に学べる素晴らしい書籍です。
テキストの「はじめに」に著者の先生が執筆の動機を書かれています。

私にとって計量経済学は高嶺の花ですが、このテキストでさまざまな回帰分析のアプローチを知ることができました。
また、書籍の Excel 処理を Python に置き換える「寄り道写経」の実践を通じて、回帰分析のお気持ちに少し近づけた感じがいたします。

回帰分析に慣れ親しむのに丁度良いレベル感と内容ですので、これはぜひともブログにしたい！と思って現在に至ります。
計量経済学の色を薄め、データ分析の色を濃いめに書いてまいります！

引用表記

この記事は、出典に記載の書籍に掲載された文章と配布データを引用し、適宜、掲載文章・配布データを改変して書いています。
【出典】
「回帰分析から学ぶ計量経済学: Excelで読み解く経済のしくみ」
第１版第１刷、著者 山澤成康、オーム社

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

第７章 反実仮想の世界

この記事は第７章の以下の節を取り扱います。

記事に用いるデータは、オーム社の書籍紹介サイトからダウンロードできる Excel ファイル内のデータをもとにしてCSVファイルを作成し、data フォルダに格納しています。

第７章で用いるライブラリをインポートします。

### インポート

# 数値計算
import numpy as np
import pandas as pd

# 統計処理
import scipy.stats as stats
import pingouin as pg
import statsmodels.formula.api as smf

# 効果検証
from linearmodels import IV2SLS                        # 操作変数法
import rdrobust                                        # 回帰不連続デザイン
import rddensity                                       #  同上
from causallib.estimation import PropensityMatching    # 傾向スコアマッチング
from causallib.estimation import IPW                   # 傾向スコア・IPW法
from causallib.evaluation import evaluate              # 傾向スコアの評価

# 機械学習
from sklearn.linear_model import LogisticRegression    # ロジスティック回帰

# グラフ描画
from graphviz import Digraph

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
import japanize_matplotlib

7.12 傾向スコアマッチング

■ 傾向スコアマッチングの概要 p.238～
テキストによると「傾向スコアマッチングは、似たタイプの変数をペアにして政策実施など介入にあったかどうかを見る方法」です。
「政策を実施する群」（処置群）と「実施しない群」（対照群）をランダムに割り付けすることが現実的には難しい場合に、既存のデータでランダムに抽出したと同様の効果を得ようとするのが傾向スコアマッチングとのことです。

テキストの計算手続きを要約すると・・・

• 原因変数（処置）を目的変数、共変量を説明変数にして回帰分析（ロジスティック回帰など）を実行し、目的変数＝傾向スコアを得る

• 傾向スコアが同じ（近似）する処置群と対照群をマッチングしてペアをつくる

• 処置群と対照群のペアについて、対応のある平均の差の検定を行って、平均値の差（施策の効果）が有意か調べる

■ 傾向スコアマッチングの実践 p.239

7.8節で利用した職業訓練の年収に与える影響データをここでも利用します。
データを読み込みます。

### データの読み込み
# 年収⇒万円、訓練⇒1:職業訓練を受けた、2:受けなかった

df12 = pd.read_csv('./data/07_05_jobtrain.csv').reset_index()
print('df12.shape:', df12.shape)
display(df12)

【実行結果】
年収（単位：万円）が結果変数、訓練が原因変数（処置）であり、0：訓練未実施、1：訓練実施です。
年齢と就学年数が共変量です。

最初のステップの「ロジスティック回帰を実行して傾向スコアを算出」します。
訓練を目的変数、年齢と就学年数を説明変数にします。
statsmodels の logit を利用します。

### 傾向スコアマッチング ステップ１ 傾向スコアの算出

# ロジスティック回帰の実行
result = smf.logit(formula='訓練 ~ 年齢 + 就学年数', data=df12).fit()
# 回帰の結果の表示
display(result.summary())

【実行結果】
係数の推定値からテキストのとおり、年齢が低いほど、就学年数が少ないほど職業訓練を受ける傾向が強いことが分かります。

傾向スコアを算出します。
「訓練」の確率値を傾向スコアに用います。

## 予測値＝傾向スコアの取得
df12['傾向スコア'] = result.predict()
display(df12)

【実行結果】
傾向スコアが同一の複数行が目に入ってきます。。。

２番目のステップ「傾向スコアが同じ処置群と対照群をマッチング」します。
テキストのデータは傾向スコアが一致するペアを作れるように出来ているので、完全一致でマッチング処理をしてみます。

### 傾向スコアマッチング ステップ２ マッチング処理 ※復元なしの完全一致マッチング

## マッチング用データの作成
# 訓練あり1と訓練なし0に分割, 列名にT1,T0を付与して識別
T1 = df12[df12['訓練']==1].set_axis([f'{c}T1' for c in df12.columns], axis=1)
T0 = df12[df12['訓練']==0].set_axis([f'{c}T0' for c in df12.columns], axis=1)

## 傾向スコア・マッチング処理
# マッチしたデータを格納するデータフレームの初期化
T1_match = pd.DataFrame()
T0_match = pd.DataFrame()
# 訓練ありT1の各行ごとに傾向スコアが一致する訓練なしT0の取得を繰り返し処理
for T1idx in T1.index:
    for T0idx in T0.index:
        # 訓練ありT1の傾向スコアに一致する訓練なしT0を取得
        if T0.loc[T0idx, '傾向スコアT0'] == T1.loc[T1idx, '傾向スコアT1']:
            # マッチングした順にmatchデータフレームに追加格納
            T0_match = pd.concat([T0_match, T0[T0.index==T0idx]], axis=0)
            T1_match = pd.concat([T1_match, T1[T1.index==T1idx]], axis=0)
            break
# マッチングしたT1とT0とを１つのデータフレームに結合
matching_df = pd.concat([T1_match.reset_index(drop=True),
                         T0_match.reset_index(drop=True)], axis=1)
# 年収の差分列の追加
matching_df['年収差'] = (
            matching_df['年収T1'] - matching_df['年収T0'])

# 平均データフレームの作成
matching_df_mean = matching_df.mean(numeric_only=True).rename('平均').to_frame().T

## データフレームの表示
display(matching_df.round(2))
display(matching_df_mean.round(2))

【実行結果】
３対のペアを作れました。

訓練実施（T1）の平均年収は 533.8 万円、訓練未実施（T0）の平均年収は 453.4 万円、平均の差（訓練の効果候補）は 100.5 万円です。

最後のステップです。
平均の差が有意かどうかを対応のある平均の差の検定で調べます。
scipy.stats の ttest_rel を利用します。

### 対応する２標本の平均値の差のｔ検定 scipy.stats利用
stats.ttest_rel(matching_df['年収T1'], matching_df['年収T0'],
                alternative='greater')

【実行結果】
$${p}$$値は$${0.00005}$$であり、$${5\%}$$水準で平均の差は有意となりました。
テキスト的には、職業訓練の効果はおよそ 100 万円、ということだと思います。

ちなみに pingouin で対応のある平均の差の検定を行うとこうなります。

# pingouin利用
pg.ttest(matching_df['年収T1'], matching_df['年収T0'], paired=True,
         alternative='greater')

【実行結果】

■ causallib ライブラリで傾向スコアマッチング
因果推論用のライブラリ「causallib」を用いて傾向スコアマッチングを行い、平均処置効果（ATE）を推定します。
書籍「因果推論」を参考にしました。

データを共変量、原因変数（処置変数）、結果変数に分けます。

### causallibを用いた傾向スコア
# 参考: 書籍「因果推論」

## データの準備
X = df12[['年齢', '就学年数']]   # 共変量（説明変数）
a = df12['訓練']                 # 処置変数
y = df12['年収']                 # 結果変数

では傾向スコアマッチングを行って、平均処置効果を推定します。
ロジスティック回帰には scikit-learn の LogisticRegression を利用します。

### causallibを用いた傾向スコアマッチング マッチングは最近隣法
# 引用: 書籍「因果推論」p82～

# 傾向スコアを算出するモデルの定義
learner = LogisticRegression(
    solver='liblinear',            # 線形モデル向けの最適化アルゴリズム
    class_weight='balanced',       # 正例と負例のバランスを取るために使用
)

# 傾向スコアの算出とマッチング
pm = PropensityMatching(learner=learner)
pm.fit(X, a, y)

# ATEの計算
outcomes = pm.estimate_population_outcome(X, a)
effect = pm.estimate_effect(outcomes[1], outcomes[0])  # ATEを推定

# 結果の表示
print(f'Average outcoms of Untreated: {outcomes[0]:7.3f}')
print(f'Average outcoms of Treated  : {outcomes[1]:7.3f}')
print(f"Average Treatment Effect    : {effect['diff']:7.3f}")

【実行結果】
平均処置効果（ATE）は 52.8 万円です。
先に求めた平均値の差 100.5 万円とは大きく離れています。。。

課題1 英語のテストと研修の効果（RCT） p.241

２つの会社の英語研修の効果を平均の差の検定で分析する課題です。
６月に条件に応じて英語研修を受け、５月と７月に英語の試験を比較して研修の効果を確認します。

■ A社の英語研修の効果
A社は５月試験後、全員が６月に英語研修を受講して７月試験を受けています。
データを読み込みます。

### データの読み込み
# 研修前: 5月試験の得点、研修後: 7月試験の得点

df13a = pd.read_csv('./data/07_ex1_a.csv')

print('df13a.shape:', df13a.shape)
display(df13a)

【実行結果】
10名の５月試験の得点（研修前）と７月試験の得点（研修後）です。
同じ人が５月にも７月にも試験を受けています。

要約統計量を確認します。

## A社の分析 要約統計量
display(df13a.describe().round(1))

【実行結果】
平均値を見ると、研修前 $${78.6}$$点、研修後$${88.0}$$点となり、研修後のほうが平均値が$${9.4}$$点高くなっています。
研修の効果があったのでしょうか！？

しばらく可視化でデータを眺めてみます。

### A社の分析 棒グラフ
plt.figure(figsize=(6, 3))
sns.barplot(data=df13a, alpha=0.7)
plt.ylabel('得点');

【実行結果】
試験の得点は研修前後で変化している感じです。

個人ごとに研修前後の得点を比較してみます。

### A社の分析 個人の比較棒グラフ
plt.figure(figsize=(6, 3))
plt.bar(df13a.index-0.25, df13a['研修前'], width=0.4, label='5月試験', alpha=0.7)
plt.bar(df13a.index+0.25, df13a['研修後'], width=0.4, label='7月試験', alpha=0.7)
plt.xlabel('個人')
plt.ylabel('得点')
plt.legend(bbox_to_anchor=(1, 1))
plt.show()

【実行結果】
研修受講後に得点が増えた人も減った人もいます。

箱ひげ図で分布を見てみます。

### A社の分析 箱ひげ図
plt.figure(figsize=(6, 3))
sns.boxplot(data=df13a, fill=False)
plt.ylabel('得点');

【実行結果】
研修後のほうが得点のバラツキが小さくなっています。
研修による「底上げ」があったのかもしれません。

では平均値の差の検定を実行します。
同じ人の得点を比べていますので、対応のある２標本の平均値の差の検定を行います。
scipy.stats の ttest_rel を利用します。

### A社の分析 対応のある２標本のt検定
# p値が大きい⇒7月試験のほうが大きい（差がある）とは言えない
stats.ttest_rel(df13a['研修後'], df13a['研修前'], alternative='greater')

【実行結果】
$${p}$$値は$${0.113}$$。$${5\%}$$水準で平均値の差は有意ではない、という結果になりました。
A社の場合、研修に効果があったとは言えない状況です。。。

■ B社の英語研修の効果
B社はランダムに２つのグループに分けて、１つ目のグループは６月に英語研修を受講し、２つ目のグループは６月に英語研修を受講しませんでした。
データを読み込みます。

### データの読み込み
# 研修前: 5月試験の得点、研修後: 7月試験の得点

df13b0 = pd.read_csv('./data/07_ex1_b0.csv')
df13b1 = pd.read_csv('./data/07_ex1_b1.csv')

print('df13b0.shape:', df13b0.shape)
display(df13b0)
print('df13b1.shape:', df13b1.shape)
display(df13b1)

【実行結果】
変数名 df13b0 は研修を受けていないグループです。

変数名 df13b1 は研修を受けたグループです。

データの前処理を行います。
２つのデータを統合し、研修前後の得点差を表す「差」列、研修を受けたかどうかを表す「研修」列（1：受けた、0：受けていない）を追加します。

### B社の分析 データの前処理

# 研修前と研修後の差の列を追加
df13b0['差'] = df13b0['研修後'] - df13b0['研修前']
df13b1['差'] = df13b1['研修後'] - df13b1['研修前']

# 研修受講有無の列を追加
df13b0['研修'] = 0
df13b1['研修'] = 1

# 2つのデータフレームを統合
df13b = pd.concat([df13b0, df13b1], axis=0)
print('df13b.shape: ', df13b.shape)
display(df13b)

【実行結果】

しばらくデータを眺めます。
研修を受けたグループの要約統計量です。

### B社の分析 研修を受けたケースの要約統計量
display(df13b1.describe().round(1))

【実行結果】
研修前後の得点差の平均は $${13.7}$$点。
大幅な得点増加であり、研修の効果を期待してしまいます。
ただ、得点差の標準偏差が$${9.4}$$点であり、バラツキが大きいかもしれません。

研修を受けていないグループの要約統計量です。

### B社の分析 研修を受けなかったケースの要約統計量
display(df13b0.describe().round(1))

【実行結果】
こちらも研修前後の得点差は$${11.1}$$点であり、大きな値です。
研修を受講しなくても得点が上がっているとすると、７月の試験は易しかったのでしょうか？独学で得点を伸ばしたのでしょうか？

ちなみに研修受講グループと研修未受講グループの得点差の平均値の差は$${13.7-11.1=2.6}$$点です。

グループ別の箱ひげ図を見てみます。

### B社の分析　箱ひげ図
# 研修を受講したグループは得点差のバラツキが小さい、得点減少がほぼない

sns.boxplot(data=df13b, y='差', hue='研修', width=0.5, gap=0.3, fill=False);
plt.axhline(0, color='tab:red', ls='--')
plt.grid(lw=0.5);

【実行結果】
青色が研修未受講グループ、オレンジ色が研修受講グループです。
研修受講グループのほうが得点差のバラツキが小さく、また、得点の減少者がほぼいない、そんな感じです。

では平均値の差の検定を実行します。
同一人物の得点差を比べないので、対応のない２標本の平均値の差の検定を行います。
scipy.stats の ttest_ind を利用して、得点差を調べます。
また、等分散を仮定しないので、引数 equal_var=False とし、ウェルチのｔ検定を指定します。

### 対応のない２標本の平均値の差のｔ検定 ★書籍の自由度・p値と異なる
# p値が大きいので、研修受講したグループの得点差が大きいとは言えない
stats.ttest_ind(df13b[df13b['研修']==1]['差'], df13b[df13b['研修']==0]['差'],
                equal_var=False, alternative='greater')

【実行結果】
$${p}$$値は$${0.381}$$。$${5\%}$$水準で有意ではないということです。
B社も英語研修に効果があったとは言えない状況です。。。

なおテキストと自由度と$${p}$$値が異なります。
テキストがｔ検定に使用するExcelの分析ツールは整数値の自由度を使っていますが、scipy.stats の ttest_ind は小数を含む自由度を使っているからです。

回帰分析を用いた研修の効果推定をやってみます。
得点差を目的変数、研修の実施有無を説明変数にしました。

### 回帰の分析の実行
# 研修有無の係数のp値は大きいので、研修有無の係数が0（効果なし）を受容する

result = smf.ols(formula='差 ~ 研修', data=df13b).fit()
display(result.summary())

【実行結果】
研修（実施有無）の係数の推定値は、研修受講グループと研修未受講グループの得点差の平均値の差$${13.7-11.1=2.6}$$点と同じになっています。
係数の推定値の$${p}$$値は$${0.76}$$であり、$${5\%}$$水準で有意ではない結果になっています。

課題2 少子化対策（傾向スコアマッチング） p.242

少子化対策を「実施した10自治体」と「実施しなかった10自治体」のデータを用いて、傾向スコアマッチングによる少子化対策の効果を分析します。

データを読み込みます。

### データの読み込み

df14T1 = pd.read_csv('./data/07_ex2_T1.csv')
df14T0 = pd.read_csv('./data/07_ex2_T0.csv')

print('df14T1.shape:', df14T1.shape)
display(df14T1)
print('df14T0.shape:', df14T0.shape)
display(df14T0)

【実行結果】
変数名 df14T1 は少子化対策を実施した自治体です。
既に傾向スコアは算出されています。

変数名 df14T0 は少子化対策を実施しなかった自治体です。
既に傾向スコアは算出されています。

少子化対策実施自治体と少子化対策未実施自治体を傾向スコアを用いてマッチングします。

### 傾向スコアマッチング マッチング処理 ※復元なしの完全一致マッチング

## マッチング用データの作成
# 訓練あり1と訓練なし0に分割, 列名にT1,T0を付与して識別
T1 = df14T1.set_axis([f'{c}T1' for c in df14T1.columns], axis=1)
T0 = df14T0.set_axis([f'{c}T0' for c in df14T0.columns], axis=1)

## 傾向スコア・マッチング処理
T0_match = pd.DataFrame()
T1_match = pd.DataFrame()
# 訓練ありT1の各行ごとに傾向スコアが一致する訓練なしT0の取得を繰り返し処理
for T1idx in T1.index:
    for T0idx in T0.index:
        # 訓練ありT1の傾向スコアに一致する訓練なしT0を取得
        if T0.loc[T0idx, '傾向スコアT0'] == T1.loc[T1idx, '傾向スコアT1']:
            # マッチングした順にmatchデータフレームに追加格納
            T0_match = pd.concat([T0_match, T0[T0.index==T0idx]], axis=0)
            T1_match = pd.concat([T1_match, T1[T1.index==T1idx]], axis=0)
            break
# マッチングしたT1とT0とを１つのデータフレームに結合
matching_df = pd.concat([T1_match.reset_index(drop=True),
                         T0_match.reset_index(drop=True)], axis=1)
# 政策効果の大きさの差分列の追加
matching_df['政策効果の大きさ差'] = (
            matching_df['政策効果の大きさT1'] - matching_df['政策効果の大きさT0'])
# 平均行の追加
# matching_df.loc['平均'] = matching_df.mean(numeric_only=True)

## データフレームの表示
display(matching_df.round(2))

【実行結果】
７自治体がマッチングできました。
変数名に T1 が付いているのは少子化対策実施自治体、T0 が付いているのは未実施自治体です。
「政策効果の大きさ差」はマッチングした自治体間の政策効果の大きさの差です。

では平均値の差の検定を実行します。
「政策効果の大きさ」の平均値について、少子化対策自治体と未実施自治体との間で差があるかを調べます。
マッチングしていますので、対応のある２標本の平均値の差の検定を行います。
scipy.stats の ttest_rel を利用します。

### 対応のある２標本のｔ検定
# p値が極小であり、帰無仮説「少子化対策効果の差は０である」は棄却され、
# 対策に正の効果があると言えます

stats.ttest_rel(matching_df['政策効果の大きさT1'],
                matching_df['政策効果の大きさT0'],
                alternative='greater')

【実行結果】
$${p}$$値は$${0.0000}$$であり、$${5\%}$$水準で有意です。
少子化対策の効果があったと言えます。

効果の大きさを推定してみます。
「政策効果の大きさ差」の平均値を効果として見てみます。

### 効果の概算・・・政策効果の大きさ差の平均をとる
display(matching_df['政策効果の大きさ差'].describe().round(3).to_frame())

【実行結果】
マッチングした少子化対策の実施自治体・未実施自治体の「政策効果の大きさ差の平均値は$${2.286}$$」です。

今回の写経は以上です。

第７章の課題は書籍の最後の課題ですので、この記事が書籍「回帰分析から学ぶ計量経済学」寄り道写経シリーズの最終回となります。
長期間お付き合いくださり、ありがとうございました！

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！

１．のんびり統計

統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。

２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で

書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！

３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で

書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！

４．楽しい写経 ベイズ・Python等

ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀

５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で

書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。

６．データサイエンスっぽいことを綴る

統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。

７．Python機械学習プログラミング実践記

書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。