🐍

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

2022/06/06に公開約13,000字

はじめに

線形回帰による因果推論について、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)などを利用することもあります。

共変量の選択

因果推論における回帰分析で、最も難しい問題の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の正規分布に従う)

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という値はかなり近い値で推定できていると言えます。しかし、実際のデータ分析業務で、これほどうまく効果を推定できることは稀です。多くの場合、効果を推定する上でモデルに必要な変数が抜けていたり、余計な変数が含まれていることが多々あります。

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

脱落変数(欠落変数)とは、本来であれば必要であるがモデルから抜け落ちている変数のことです。そして、脱落変数が存在することによって生じるバイアスのことを脱落変数バイアス(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に対して相関のある変数をモデルに組み込む」必要があります。

参考文献

おわりに

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

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

Discussion

ログインするとコメントできます