Pythonで因果推論(4)~回帰分析を用いた効果検証~

はじめに

線形回帰による因果推論について、Pythonによる実装を交えてまとめました。回帰分析の特性の導出や理論的な背景については記述しておりません。内容について誤り等ございましたら、コメントにてご指摘いただけますと幸いです。

回帰分析の概要

回帰分析とは、説明変数$X$を1単位増減させたときに、被説明変数$Y$がどの程度変動するかを出力する分析手法です。

回帰分析については、こちらの記事もご参照ください。

因果推論における回帰分析

因果推論の文脈では、何らかの処置を表す変数$D$(以下、処置変数)が、被説明変数$Y$にどれだけの効果を与えているかを検証する際に回帰分析を利用します。例えば

$$Y = b_0 + b_DD$$

というような回帰モデルを作成し、

• 帰無仮説を「$b_D = 0$」として検定を行い、処置変数$D$が被説明変数$Y$に影響を与えていると言えるか
• パラメータ$b_D$を推定して、処置変数$D$が被説明変数$Y$にどれくらいの影響を与えているか

などを検討します。

処置変数$D$と被説明変数$Y$の両方に影響を与える交絡因子が存在するであろう場合を考えると、上記のモデルに交絡因子と予想される共変量$X_i(1 \leq i \leq n)$を追加して、

$$Y = b_0 + b_1X_1 + ... + b_nX_n + b_DD$$

という適切な重回帰モデルを作成し、パラメータ$b_D$の検定や推定を行うことで、共変量の影響を取り除いた効果の有無やその大きさを推定することができます。

今回は、線形回帰をベースとした重回帰モデルを紹介しましたが、手元のデータによってはロジスティック回帰分析やポアソン回帰分析などの一般化線形モデル(GLM: Generalized Linear Model)などを利用することもあります。

用語の確認

処置変数

意図的にとある変数を操作(介入)して効果検証を行うという文脈では介入変数と呼ぶことが一般的ですが、必ずしも介入操作が行われていなくても構わないため、便宜上定義した造語です。因果推論のテキストでは「処置変数は原因変数、被説明変数は結果変数」と呼称されていることが多いです。本記事で処置変数・被説明変数と記述しているのは、単なる私の慣習です。

共変量

本記事では、「分析者が被説明変数へに影響するであろうと予想した、処置変数を除く変数」という文脈において、その変数のことを共変量と呼ぶことにします。
これは正式にそのような定義が存在していると言うわけではなく、あくまでも私の慣習です。

共分散分析モデル

上記の重回帰モデルにおいて、特に「$D$がカテゴリカル変数」かつ「$X$が連続変数」のとき、共分散分析(ANCOVA: ANalysis of COVAriance)モデルと呼ばれます。こちらは私の慣習ではありません。笑

共変量の選択

因果推論における回帰分析で、最も難しい問題の1つに「共変量の選択」があります。うまく共変量が選択できていない回帰モデルを利用してしまうと、効果$b_D$を正しく推定することができません。そして、さらに厄介なことに、分析者が作成した回帰モデルが正しいかどうかをデータから判断することは、基本的には不可能です。

そのため、バックドア基準や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$$
    の値が0.5以上であれば特講を受講($D=1$)、0.5未満であれば特講を未受講($D=0$)とする（ただし、$noise$は一様分布(0, 0.2)に従う）
• 試験の得点$Y$(被説明変数)
  • 学生の試験の得点を表し、
    $$Y = 40 \times X_1 + 20 \times X_2 + 20 \times D + noise$$
    で表現される(ただし、$noise$は、平均0、標準偏差10の正規分布に従う)

試験の得点$Y$のモデル$\:Y = 40 \times X_1 + 20 \times X_2 + 20 \times D + noise \:$は、分析者には分かりません。因果推論の目的は、共変量である$X_1, X_2$、処置変数$D$と、試験の得点$Y$のデータをもとに、$D$の係数を推定することにあります。

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())

(出力結果)

出力結果の下側を見ると、特講を受講する学生$D=1$の方が試験の得点$Y$の平均値は高いことは確認できます。しかし、そもそも学習意欲$X_1$が高い方や理系出身$X_2=1$の方が多く、それらの変数によって生じる試験の得点$Y$の差(セレクションバイアス)を考慮する必要がありそうです。$D=1$群の試験の得点$Y$の平均と$D=0$群の試験の得点$Y$の平均を比較して差があっても、その差が特講の受講$D$に依るものなのか、はたまた学習意欲$X_1$や理系出身ダミー$X_2$に依るものも含まれているのか、現状では判断できていないという状態なんですね。

そこで、回帰分析を用いて特講受講の効果$b_D$の推定を行います。

回帰モデルの構築と評価

線形回帰

今回は、線形性を仮定し、回帰モデルを

$$Y = b_0 + b_1X_1 + b_2X_2 + b_DD$$

とします。

ここで推定したいパラメータは、特講の効果、すなわち$D$の係数にあたる$b_D$です。これを、最小二乗法(OLS)を用いて推定します。

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と非常に小さい値なので、帰無仮説(今回の場合「特講受講の効果はない($b_D=0$)」)を棄却できます。したがって、回帰分析の結果、特講を受講することで試験の得点は平均19.9008ほど増加すると解釈できます。

実際のモデル(分析者からは観測できない)においても、特講受講ダミー$D$の係数は20であることから、19.9008という値はかなり近い値で推定できていると言えます。しかし、実際のデータ分析業務で、これほどうまく効果を推定できることは稀です。多くの場合、効果を推定する上でモデルに必要な変数が抜けていたり、余計な変数が含まれていることが多々あります。

回帰分析の結果を見ると、つい

$$Y = -1.6042 + 40.4059X_1 + 20.3602X_2 + 19.9008D$$

として、回帰モデルを解釈してしまいそうですが、因果推論における回帰分析では$\hat{b_D}$以外の推定結果には議論しません。その理由は、基本的に知りたいことは$D$の効果$\hat{b_D}$であり、その他のパラメータの値が$D$の効果の解釈を助けるということはないからです。

脱落変数バイアス(欠落変数バイアス)

脱落変数(欠落変数)とは、本来であれば必要であるがモデルから抜け落ちている変数のことです。そして、脱落変数が存在することによって生じるバイアスのことを脱落変数バイアス(OVB: Omitted Variables Bias)と言います。

たとえば、先ほどの特講の受講$D$と試験の得点$Y$において、本来であれば

$$Y = b_0 + b_1X_1 + b_2X_2 + b_DD$$

であるところを、理系出身ダミー$X_2$を抜いて

$$Y = b_0 + b_1X_1 + b_DD$$

として特講の効果$b_D$を推定し、脱落変数バイアスが生じることによって推定結果にどのような影響を与えるのかを確認します。

Pythonによる実装

Pythonのコードは、説明変数から理系出身ダミー$X_2$を抜くだけなので、先ほどとほぼ同じです。

# 必要なライブラリのインポート
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もの脱落変数バイアスが生じているようです。

このような脱落変数バイアスが生じないためには「被説明変数$Y$と処置変数$D$に対して相関のある変数をモデルに組み込む」必要があります。

因果推論における回帰分析では、脱落変数バイアスが0になるように共変量を選択することが理想です。この状態は、CIA(Conditional Independence Assumption)と呼ばれ、数式では

$${Y^i_1,Y^i_0} \perp D^i | X^i$$

と表現されます。

参考文献

• 岩波データサイエンス刊行委員会編「岩波データサイエンスVol.3-[特集]因果推論-実世界のデータから因果を読む」岩波新書(2016)
• 安井「効果検証入門」ホクソエム社(2020)
• 西山、新谷、川口、奥井「New Liberal Arts Selection 計量経済学」有斐閣(2019)
• 山本「実証分析のための計量経済学」中央経済社(2015)

おわりに

統計学や計量経済学などを学習したことがある方にとって、回帰分析は馴染み深い分析手法だと思います。しかし、それと同時に奥が深い分析手法でもあり、因果推論の文脈になると一癖も二癖も出てきます。Pythonによるコーディング自体は簡単ですが、モデルの作成やその解釈によって、回帰分析から得られる情報は大きく変わります。一見シンプルで、ポピュラーな手法にもかかわらず、分析者の腕が試されるのもまた回帰分析の魅力かもしれません。

他にも下記のような記事を書いています。ご一読いただけますと幸いです。

• 相関関係と因果関係と疑似相関
• 反実仮想と因果効果
• 介入とランダム化比較試験
• 層別解析を用いた効果検証
• 傾向スコアを用いた効果検証
• 操作変数を用いた効果検証
• 回帰不連続デザインを用いた効果検証
• 差分の差分法を用いた効果検証
• 機械学習を用いた効果検証

また、過去にLTや勉強会で発表した資料は下記リンクにまとめてあります。ぜひ、ご一読くださいませ。
https://speakerdeck.com/s1ok69oo