効果検証入門 2章をPythonで書く

はじめに

効果検証入門 ～正しい比較のための因果推論／計量経済学の基礎内のソースコードをPythonで再現します。

既に素晴らしい先人の実装例がありますが、自分の勉強用のメモとして残しておきます。

この記事では、2章について記載します。
コードは、githubにも掲載しています。
なお、変数名や処理内容は、基本的に書籍内に寄せて実装します。

回帰分析

線形回帰の実装は、scikit-learnでもstatsmodelsでも可能。
しかし、各変数の統計値に関しては、statsmodelsの方が有利なので、こちらの方が便利。

scikit-learn

sklearnの回帰分析

from sklearn.linear_model import LinearRegression

# モデルの学習
X = biased_data[['treatment', 'history']]
y = biased_data['spend']
model = LinearRegression(fit_intercept=True, normalize=False).fit(X, y)

# 結果の出力
print(f'R^2: {model.score(X, y)}')
print(f'intercept: {model.intercept_}')
print(f'coefficients: {model.coef_}')

scikit-learnは、結果として得られるものが後述のstatsmodelsに比べ少ない。
学習データを基に、各統計値を算出することは出来るが、そこまでしてscikit-learnを使うメリットはないと思われる。

statsmodels

statsmodelsを使ったモデルの学習方法は複数あるようだが、Rライクなのは下記の通り。

statsmodelsの回帰分析

from statsmodels.formula.api import ols

# モデルの学習
model = ols('spend ~ treatment + history', data=biased_data).fit()

# 結果の出力
model.summary()

ちなみに、結果の出力はmodel.summary().tables[1]のようにすれば、複数の情報の中から任意の表を指定できる。

また、推定値を取得したい場合、model.paramsを使って辞書形式で一覧を参照できる。

RDataの読み込み

R形式のデータの読み込みは、rdataを使うと良い。
モジュールの使い方の詳細は、こちらを参考にした。

RDataの読み込み

import rdata

parsed = rdata.parser.parse_file('./vouchers.rda')
converted = rdata.conversion.convert(parsed)
vouchers = converted['vouchers']

まとめて回帰分析

本書では、複数の目的変数のモデルをまとめて学習しているが、Pythonで同様のことは難しそう。

まとめて回帰分析

import pandas as pd
from statsmodels.formula.api import ols

# 回帰式の定義
formula_x_base = ['VOUCH0']
formula_x_covariate = [
    'SVY', 'HSVISIT', 'AGE', 'STRATA1', 'STRATA2', 'STRATA3', 'STRATA4', 'STRATA5', 'STRATA6', 'STRATAMS',
    'D1993', 'D1995', 'D1997', 'DMONTH1', 'DMONTH2', 'DMONTH3', 'DMONTH4', 'DMONTH5', 'DMONTH6',
    'DMONTH7', 'DMONTH8', 'DMONTH9', 'DMONTH10', 'DMONTH11', 'DMONTH12', 'SEX2',
]
formula_ys = [
    "TOTSCYRS","INSCHL","PRSCH_C","USNGSCH","PRSCHA_1","FINISH6","FINISH7","FINISH8","REPT6",
    "REPT","NREPT","MARRIED","HASCHILD","HOURSUM","WORKING3",
]

# 回帰結果を受け取る関数の定義
def get_regression_result(formula, data):
  model = ols(formula, data=data).fit()
  result = pd.read_html(model.summary().tables[1].as_html(), header=0)[0]
  result.columns = ['term', 'estimate', 'std.err', 'statistic', 'p.value', '0.025', '0.975']
  result
  return result

# まとめて回帰分析を実行
results = list()
for formula_y in formula_ys:
  base_reg_formula = f'{formula_y} ~ {" + ".join(formula_x_base)}'
  base_reg_model_index = f'{formula_y}_base'
  covariate_reg_formula = f'{formula_y} ~ {" + ".join(formula_x_base+formula_x_covariate)}'
  covariate_reg_model_index = f'{formula_y}_covariate'

  base_reg_result = get_regression_result(base_reg_formula, regression_data)
  base_reg_result['model_index'] = base_reg_model_index
  results.append(base_reg_result)
  
  covariate_reg_result = get_regression_result(covariate_reg_formula, regression_data)
  covariate_reg_result['model_index'] = covariate_reg_model_index
  results.append(covariate_reg_result)

df_results = pd.concat(results).reset_index(drop=True)
df_results = df_results[['model_index', 'term', 'estimate', 'std.err', 'statistic', 'p.value', '0.025', '0.975']]

介入と目的変数の関係性のプロット

恐らくmatplotlibでエラーバー付きのプロットをするのが簡単。
信頼区間を表現するために使うエラーバーのエラーの大きさを推定値と信頼区間の差から求めるのがポイント。

(ここに掲載していない)描画に利用しているデータでは、モデルから値を取得しているが、その際の有効数字が少ないため、多少のずれが生じている可能性がある。

プロット

import matplotlib.pyplot as plt

estimate = going_private_results['estimate']
estimate_error = going_private_results['estimate'] - going_private_results['0.025']  # 信頼区間と推定値の差をエラーバーの長さとする

xmin = 0
xmax = going_private_results.shape[0] - 1

plt.errorbar(range(xmax+1), estimate, estimate_error, fmt='o')
plt.hlines(y=0, xmin=xmin, xmax=xmax, colors='k', linestyles='dashed')
plt.xlabel('model_indexe')
plt.ylabel('estimate')
plt.xticks(range(going_private_results.shape[0]), going_private_results['model_index'], rotation=45)
plt.show()

関連

• 「効果検証入門」をPythonで書いた
• RData形式のデータセットをPythonで読み込む