Pythonで因果推論(4)~回帰分析を用いた効果検証~
はじめに
線形回帰による因果推論について、Pythonによる実装を交えてまとめました。回帰分析の特性の導出や理論的な背景については記述しておりません。内容について誤り等ございましたら、コメントにてご指摘いただけますと幸いです。
回帰分析の概要
回帰分析とは、説明変数
因果推論における回帰分析
因果推論の文脈では、何らかの処置を表す変数
- 帰無仮説を「
b_D = 0 D Y - パラメータ
b_D D Y
などを検討します。
処置変数
今回は、線形回帰をベースとした重回帰モデルを紹介しましたが、手元のデータによってはロジスティック回帰分析やポアソン回帰分析などの一般化線形モデル(GLM: Generalized Linear Model)などを利用することもあります。
共変量の選択
因果推論における回帰分析で、最も難しい問題の1つに「共変量の選択」があります。うまく共変量が選択できていない回帰モデルを利用してしまうと、効果
そのため、バックドア基準やCIA(後述)を考えて、共変量を選択していく必要があります。
線形回帰による因果推論
設定
とある大学の経済学部生120を対象に、夏季休暇期間に計量経済学の特別講義(特講)を開催し、その効果を考えます。特講は、あくまでも大学の単位とは無関係の任意参加であるため、介入操作が行われない限りは受講するかどうかは学生の属性に依存すると考えます。夏季休暇終了後に、特講の受講にかかわらず対象者120人に対して計量経済学の試験を実施し、その得点のデータを観測します。
手元には下記のデータが存在するとします。
- 学習意欲
X_1 - 学生の学習意欲を0以上1以下の数値で表したもの
- 一様分布(0, 1)に従う
- 特講の受講
D Y
- 理系出身ダミー
X_2 - 学生が高校時代、理系だったか文系だったかを表したダミー変数
-
X_2=1 X_2=0 - 特講の受講
D Y
- 特講受講ダミー
D - 学生が特講を受講したかどうかを表すダミー変数
-
D=1 D=0 - 学習意欲
X_1 X_2 0.6 \times X_1 + 0.2 \times X_2 + noise D=1 D=0 noise
- 試験の得点
Y - 学生の試験の得点を表し、
Y = 40 \times X_1 + 20 \times X_2 + 20 \times D + noise noise
- 学生の試験の得点を表し、
Pythonによる実装
Pythonを用いて、上記のデータを作成します。
# 必要なライブラリをインポート
import numpy as np
import pandas as pd
# シードの設定
np.random.seed(0)
# データのサイズ
size = 120
# 学習意欲を表す変数X1(交絡因子)
X1 = np.random.uniform(0, 1, size)
# 理系出身ダミーX2(交絡因子)
X2 = np.random.choice([0, 1], p=[0.5, 0.5], size=size)
# 特講受講ダミーD(処置変数)
D_noise = np.random.uniform(0, 0.2, size)
D_threshold = 0.6*X1 + 0.2*X2 + D_noise
D = np.where(D_threshold>=0.5, 1, 0)
# 試験の得点Y
Y_noise = np.random.normal(0, 10, size)
Y = np.clip(40*X1 + 20*X2 + 20*D + Y_noise, 0, 100).astype("int")
# データフレームに格納し、5つだけ出力
df = pd.DataFrame({"学習意欲": X1, "理系出身ダミー": X2, "特講受講ダミー": D, "試験の得点": Y})
df.head()
(出力結果)
理系出身者数や特講受講者数を確認してみます。
# 理系出身ダミーの数の出力
print(df["理系出身ダミー"].value_counts())
# 特講の受講者数の出力
print(df["特講受講ダミー"].value_counts())
(出力結果)
1 61
0 59
Name: 理系出身ダミー, dtype: int64
0 62
1 58
Name: 特講受講ダミー, dtype: int64
どちらもほぼほぼ半分ずつ存在していることが分かります。
次に、要約統計量と、特講受講ダミー別の試験の得点の平均を確認してみます。
# 要約統計量の算出
display(df.describe())
# 特講を受講した学生としていない学生で平均値を算出
display(df.groupby(by="特講受講ダミー").mean())
(出力結果)
出力結果の下側を見ると、特講を受講する学生
そこで、回帰分析を用いて特講受講の効果
回帰モデルの構築と評価
線形回帰
今回は、線形性を仮定し、回帰モデルを
ここで推定したいパラメータは、特講の効果、すなわち
Pythonによる実装
Pythonのstatsmodelsを用いて、回帰分析を実装します。
# 必要なライブラリのインポート
import statsmodels.api as sm
# 説明変数
X = df[["学習意欲", "理系出身ダミー", "特講受講ダミー"]]
X = sm.add_constant(X)
# 被説明変数
y = df["試験の得点"]
# モデルの設定(OLS：最小二乗法を指定)
model = sm.OLS(Y, X)
# 回帰分析の実行
results = model.fit()
# 結果の詳細を表示
print(results.summary())
(出力結果)
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.882
Model: OLS Adj. R-squared: 0.879
Method: Least Squares F-statistic: 288.7
Date: Sun, 05 Jun 2022 Prob (F-statistic): 1.27e-53
Time: 06:16:28 Log-Likelihood: -436.24
No. Observations: 120 AIC: 880.5
Df Residuals: 116 BIC: 891.6
Df Model: 3
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const -1.6042 2.125 -0.755 0.452 -5.813 2.604
学習意欲 40.4059 5.152 7.843 0.000 30.202 50.610
理系出身ダミー 20.3602 2.008 10.138 0.000 16.382 24.338
特講受講ダミー 19.9008 3.209 6.202 0.000 13.546 26.256
==============================================================================
Omnibus: 0.727 Durbin-Watson: 2.238
Prob(Omnibus): 0.695 Jarque-Bera (JB): 0.570
Skew: 0.169 Prob(JB): 0.752
Kurtosis: 3.008 Cond. No. 9.93
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
/usr/local/lib/python3.7/dist-packages/statsmodels/tsa/tsatools.py:117: FutureWarning: In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only
x = pd.concat(x[::order], 1)
warningが出ていますが、そこは触れないでおきます。
特講受講ダミーのcoefの値を確認してみると、19.9008と推定されていることが分かります。また、p値(P>|t|)も0.000と非常に小さい値なので、帰無仮説(今回の場合「特講受講の効果はない(
実際のモデル(分析者からは観測できない)においても、特講受講ダミー
脱落変数バイアス(欠落変数バイアス)
脱落変数(欠落変数)とは、本来であれば必要であるがモデルから抜け落ちている変数のことです。そして、脱落変数が存在することによって生じるバイアスのことを脱落変数バイアス(OVB: Omitted Variables Bias)と言います。
たとえば、先ほどの特講の受講
Pythonによる実装
Pythonのコードは、説明変数から理系出身ダミー
# 必要なライブラリのインポート
import statsmodels.api as sm
# 説明変数
X = df[["学習意欲", "特講受講ダミー"]]
X = sm.add_constant(X)
# 被説明変数
y = df["試験の得点"]
# モデルの設定(OLS：最小二乗法を指定)
model = sm.OLS(Y, X)
# 回帰分析の実行
results = model.fit()
# 結果の詳細を表示
print(results.summary())
(出力結果)
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.777
Model: OLS Adj. R-squared: 0.773
Method: Least Squares F-statistic: 204.1
Date: Sun, 05 Jun 2022 Prob (F-statistic): 7.04e-39
Time: 06:17:58 Log-Likelihood: -474.31
No. Observations: 120 AIC: 954.6
Df Residuals: 117 BIC: 963.0
Df Model: 2
Covariance Type: nonrobust
==============================================================================
coef std err t P>|t| [0.025 0.975]
------------------------------------------------------------------------------
const 9.8890 2.457 4.024 0.000 5.022 14.756
学習意欲 21.7001 6.577 3.299 0.001 8.674 34.726
特講受講ダミー 36.7655 3.752 9.800 0.000 29.336 44.196
==============================================================================
Omnibus: 1.347 Durbin-Watson: 2.337
Prob(Omnibus): 0.510 Jarque-Bera (JB): 1.417
Skew: 0.206 Prob(JB): 0.492
Kurtosis: 2.662 Cond. No. 8.07
==============================================================================
Warnings:
[1] Standard Errors assume that the covariance matrix of the errors is correctly specified.
/usr/local/lib/python3.7/dist-packages/statsmodels/tsa/tsatools.py:117: FutureWarning: In a future version of pandas all arguments of concat except for the argument 'objs' will be keyword-only
x = pd.concat(x[::order], 1)
またしても、warningが出ていますが、そこは触れないでおきます。
回帰分析の結果を見ると、特講受講ダミーのcoefの値は36.7655となっており、特講を受講の効果が36.7655と推定されていることが分かります。本来の効果(実際には観測できない)が20ですので、16.7655もの脱落変数バイアスが生じているようです。
このような脱落変数バイアスが生じないためには「被説明変数
参考文献
- 岩波データサイエンス刊行委員会編「岩波データサイエンスVol.3-[特集]因果推論-実世界のデータから因果を読む」岩波新書(2016)
- 安井「効果検証入門」ホクソエム社(2020)
- 西山、新谷、川口、奥井「New Liberal Arts Selection 計量経済学」有斐閣(2019)
- 山本「実証分析のための計量経済学」中央経済社(2015)
おわりに
統計学や計量経済学などを学習したことがある方にとって、回帰分析は馴染み深い分析手法だと思います。しかし、それと同時に奥が深い分析手法でもあり、因果推論の文脈になると一癖も二癖も出てきます。Pythonによるコーディング自体は簡単ですが、モデルの作成やその解釈によって、回帰分析から得られる情報は大きく変わります。一見シンプルで、ポピュラーな手法にもかかわらず、分析者の腕が試されるのもまた回帰分析の魅力かもしれません。
他にも下記のような記事を書いています。ご一読いただけますと幸いです。
- 相関関係と因果関係と疑似相関
- 反実仮想と因果効果
- 介入とランダム化比較試験
- 層別解析を用いた効果検証
- 傾向スコアを用いた効果検証
- 操作変数を用いた効果検証
- 回帰不連続デザインを用いた効果検証
- 差分の差分法を用いた効果検証
- 機械学習を用いた効果検証
また、過去にLTや勉強会で発表した資料は下記リンクにまとめてあります。ぜひ、ご一読くださいませ。
Discussion