🐍

Pythonで因果推論(3)~介入とランダム化比較試験~

2022/05/31に公開約8,100字

はじめに

介入やランダム化比較実験(RCT)について、Pythonによる実装を交えてまとめました。本記事では、グラフ的な表現や調整に関する記述はなく、介入操作の概要と(調整を必要としない)ランダム化比較実験についてのみ取り扱っています。内容について誤り等がありましたら、コメントにてご指摘いただけますと幸いです。

介入

介入とは、「因果推論をする際に、とある変数の値を変化させる操作」のことを表します。そして、多くの場合では介入操作の因果効果を推定することが、その因果推論の目的となっています。

介入操作の具体例

ここで具体例として、こちらの記事で用いた「とある大学に所属する経済学部生の、計量経済学の試験の得点Yに対する特別講義(以下、特講)受講Dの効果」を考えたいと思います。

こちらの記事では、特講の受講するかどうかDは学生個人の学習意欲Xに依存していました。すなわち、学習意欲Xが高い学生ほど特講を受講し(D=1)、学習意欲Xの低い学生は特講を受講しない(D=0)傾向にあったわけです。ここで注意してほしいのが、この場面では特講を受講する・しないというのは介入ではありません。あくまでも観察によって得られた結果です。このとき、処置群における平均処置効果ATTは

ATT = E(Y_1-Y_0|D=1)
と表されます。

それに対して、例えば「出席番号がi番未満の学生は特講を受講し(D=1)、i番以上の学生は特講を受講しない(D=0)でください」というように、(無理やり)特講の受講を割り当てる場合、特講を受講する・しないというのは介入になります。このとき、介入により研修を受講させたときの効果は、介入操作を表す

do(D=1)
という数式を用いて
E(Y_1-Y_0|do(D=1))
と表されます。

介入操作を行う目的

介入操作を行う目的は、効果を推定しやすくすることにあります。

先ほどの「特講の受講Dと試験の得点Y」を例に挙げます。

本来であれば、学習意欲Xが高い学生ほど、特講を受講し(D=1)、試験の得点Yも高いという傾向にあります。そのため、特講を受講した学生の方が試験の得点が高い傾向にあるという結果が出ても、それが特講を受講した効果なのか、学生の学習意欲に依る差なのか、正しく評価できないという問題点がありました。

それに対し、介入操作を行うと、学習意欲Xが高い学生ほど特講を受講する(D=1)という因果構造がなくなります。そのため、特講を受講した学生と特講を受講していない学生の試験の得点を比較することで、特講の効果を評価できるようになるのです。

ランダム化比較実験(RCT)

ランダム化比較実験(RCT: Randomized Controlled Trial)とは、介入をランダム(無作為)に割り振って実験し、その結果として得られたデータを比較する分析手法です。

介入の有無をランダムに割り振ることによって、介入が行われるサンプルと行われないサンプルにおけるそのほかの要因も平均的には同一になることが期待できるというのが基本的な考え方です。

先ほどの「特講の受講Dと試験の得点Y」を例に挙げると、特講の受講をランダムに割り振ることによって、特講を受講する(D=1)学生の学習意欲Xと特講を受講しない(D=0)学生の学習意欲Xの平均が同じくらいになるということですね。

Pythonによる実装

Pythonで、疑似データを作成し、介入Dをランダムに割り振り、その結果を比較してみたいと思います。

疑似データの作成

  • 学習意欲X: 一様分布(0, 100)に従う
  • 特講の受講D: ランダムに割り振る
  • 試験の得点Y: 学習意欲Xと特講の受講Dに依存する
    • Y = 0.5X + 20D + noise
    • noiseは、正規分布(0, 10)に従う
    • Yは、0 \leq Y \leq 100を満たす整数
# 必要なライブラリをインポート
import numpy as np
import pandas as pd

# シードの設定
np.random.seed(0)
# でーたのサイズ
size = 60

# 学習意欲を表す変数X
X = np.random.uniform(0, 100, size)
# 特講の受講を表す変数D
D = np.random.choice([0, 1], p=[0.5, 0.5], size=size)

# 試験の得点Y
Y_noise = np.random.normal(0, 10, size)
Y = np.clip(0.5*X + 20*D + Y_noise, 0, 100).astype("int")

# データフレームに格納し、5つだけ出力
df = pd.DataFrame({"学習意欲": X, "特講の受講": D, "試験の得点": Y})
df.head()

(出力結果)

次に介入の割り振られ方を確認してみます。

# 介入の割り振られ方
df["特講の受講"].value_counts()

(出力結果)

0    30
1    30
Name: D, dtype: int64

人数は均等に割り振られていますね。

次に、2群の平均値に注目してみます。

# 2群の平均値の算出
df.groupby(by="特講の受講").mean()

(出力結果)

学習意欲Xの平均はほぼ等しいとまではいきませんがそれほど差がないのに対し、試験の得点Yの平均は介入群(D=1)の方がかなり高いですね。

ヒストグラムも確認してみます。

# 必要なライブラリのインポート
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib

# do(D=1)
df_D1 = df[df["特講の受講"]==1]
# do(D=0)
df_D0 = df[df["特講の受講"]==0]

# ヒストグラムの描画
fig = plt.figure(figsize=(10, 3))

# 学習意欲Xのヒストグラム
ax1 = plt.subplot(121)
ax1.hist(df_D1["学習意欲"], alpha=0.4, label="特講受講")
ax1.hist(df_D0["学習意欲"], alpha=0.4, label="特講未受講")
ax1.set_title("学習意欲")

# 試験の得点Yのヒストグラム
ax2 = plt.subplot(122)
ax2.hist(df_D1["試験の得点"], alpha=0.4, label="特講受講")
ax2.hist(df_D0["試験の得点"], alpha=0.4, label="特講未受講")
ax2.set_title("試験の得点")

plt.legend()
plt.show()

(出力結果)

学習意欲Xは介入に関わらず似たような分布をしていますが、試験の得点Yは介入群の方がかなり分布が右に寄っていることが分かりますね。

次に、特講の受講Dはどれくらい試験の得点Yに効果があるかを確認します。介入の効果は

E(Y_1|do(D=1)) - E(Y_0|do(D=0))
で表されるので、これを計算します。

# 介入の効果
df_D1["特講の受講"].mean() - df_D0["特講の受講"].mean()

(出力結果)

15.466666666666665

よって、特講の受講には、試験の得点を15点上げる効果があるというように解釈できます。

今回は明らかに特講の受講に効果があると言えそうですが、一応最後に有意差検定もしておきたいと思います。この場合、試験の得点Yのデータに正規性があるとは言い切れないので、有意水準を0.05、帰無仮説を「2群間の平均値に差はない」とし、マン・ホイットニーのU検定で検定します。

# 必要なライブラリをインポート
from scipy import stats

# マン・ホイットニーのU検定
result = stats.mannwhitneyu(df_D1["試験の得点"], df_D0["試験の得点"], alternative="two-sided")
# p値を出力
print(f"p値: {result.pvalue}")

(出力結果)

p値: 0.0015851489181872057

p < 0.05なので、帰無仮説は棄却され2群間の平均には差がある、すなわち特講には試験の得点を上げる効果があるということが分かりました。

参考文献

おわりに

ランダム化比較実験(RCT)は、実際のデータ分析現場で行うということはなかなか現実的ではありませんが、介入について理解する上では非常に有用だと思います。また、介入は因果推論・効果検証を行う上で不可欠と言っても過言ではありません。介入はまだまだ押さえておく必要のある項目が山ほどあるので、少しずつまとめていけたらと思います。

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

Discussion

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