Rで学ぶ計量経済学と機械学習 5
計量経済学１：差の検定/回帰分析/マッチング
安藤道人（立教大学） 三田匡能 (株式会社 GA technologies)
2023-05-21
1 本講義の目的
- 差の検定と回帰分析とマッチング法の基礎を学ぶ。
- データを自分で生成する方法を学ぶ。
2 平均処置効果
今回から、データの背後に潜む因果関係を推定する手法について学んでいく。
計量経済学においてよく用いられる考え方は潜在的結果の枠組み（potential outcome framework）である。 今回からこの枠組みに関する用語が登場するため、ここで軽く触れておく。
潜在的結果の枠組みのもとでは、因果関係の強さは「同じ個体が、処置を受けた時と処置を受けなかったときの結果の差」として考える。
2.1 潜在的結果
個体\(i\)が「処置を受けた場合（\(X_i=1\)）」と「処置を受けなかった場合（\(X_i=0\)）」の結果\(Y_i\)を、実際に処置があったかどうかに関わらず議論するために定めたものを潜在的結果（potential outcome）といい、それぞれ\(Y_i(1), Y_i(0)\)と書く。
実際の観察結果は\(Y_i\)とする
\[ Y_i = \left\{ \begin{array}{l} Y_i(0) \ \text{if } X_i = 0 \\ Y_i(1) \ \text{if } X_i=1 \end{array} \right. \]
2.2 処置効果
個体\(i\)が処置を受けた時と受けなかったときの結果の差を個体処置効果 （individual treatment effect: ITE）という。
\[ \tau_i = Y_i(1)-Y_i(0) \]
処置を受けたときの結果と処置を受けなかったときの結果の両方を同時に観測することはできないため、ITE \(\tau_i\)は観測することはできない。
そこで、処置を受けた群とそうでない群における結果の平均の差を用いて処置効果を推定することを考える。
\[ \tau = E[Y(1) - Y(0)] = E[Y(1)] - E[Y(0)] \]
これを平均処置効果（average treatment effect: ATE）という。
以下の節では、データを分析しながらATEの推定方法について学んでいく。
3 データの用意
データから因果関係を推定するということは、データ生成過程（Data generating process: DGP）を理解しようという試みでもある。そこで今回は、データ生成過程があらかじめ分かっている架空のデータを自分で生成し、それを分析していく。
個人の所得\(Y\)と学歴\(X\)（大卒か否か）と能力\(A\)との関係について、以下のようなDGPを持つ架空データを次のように生成する（ただし\(\varepsilon\)は撹乱項）。このDGPは現実的なものではないが、理解しやすいようにできるだけ単純に設定したものを用いる。
\[ Y = 200 + 10A + 500X+ \varepsilon \]
- サンプルサイズは10,000人
- 能力は0から100まで均等に分布する
- 能力が1上がると所得は10万円上昇する
- 能力が80 以上の約 2,000人中から約 1,000人がランダムに選ばれて大卒となる。
- 大卒だと所得が500万円上昇する
- 切片は200万円（能力=0、非大卒の平均給与）
ここで、学歴\(X\)が所得\(Y\)に与える効果に興味があるとき、\(X\)を処置変数、\(Y\)をアウトカム変数と呼ぶ。また、\(X\)と\(Y\)に影響を及ぼすことが想定される能力\(A\)のような変数を共変量/交絡要因/コントロール変数と呼ぶ（分野によって呼び方やニュアンスが若干異なるが、ほぼ同じ意味で用いて問題ない）。
まず、上記のDGPに従うデータを生成する。なお、データ生成は、これまでの学習した内容と比較してやや高度なので、とりあえずは、ただコピペしてデータ生成するだけでも構わない。
# 事前準備 --------------------
# パッケージの読み込み
pacman::p_load(tidyverse, stargazer, Matching)
# 乱数の種を固定（同じ種からは常に同じ乱数が生成される。0を他の数値に変えると異なる乱数となる。）
set.seed(0)
# データの生成 ----------------
# n:サンプルサイズ
n <- 10000
# 能力は0から100まで均等に分布。#runifは一様分布を発生させる関数。標本規模n、最小値0、最大値100
ability <- runif(n, min = 0, max = 100)
# IDとabilityをデータフレームに格納する。
# 以下の"tibble()"はtidyverseにおけるデータフレームを作成する関数。
# 代わりに"data.frame()"を用いても構わない。
df <- tibble(ID = 1:n, ability)
# 大卒ダミーの作成
# 能力が 80 以上の約 2000人の中から約 1000人をランダムに選ばれて大卒にする。
# dfからdplyr::filter()で抽出し、sample_fram()でさらに半分をランダムに抽出する。
university_df <- df %>% filter(ability >= 80) %>%
sample_frac(0.5) # 大卒の人
# university_df のデータフレームに、universityという変数を作成し、すべて1とする。
university_df <- university_df %>%
dplyr::mutate(university = 1)
## なお、university_df["university"] = 1でも作成可能
# dfからdplyr::anti_join()を用いて"university_df"とマッチしなかった人を抽出する。
no_university_df <- anti_join(df, university_df,
by = c("ID","ability")) # 大卒ではない人
# no_university_dfのデータフレームに、universityという変数を作成し、すべて0とする。
no_university_df <- no_university_df %>%
dplyr::mutate(university = 0)
## なお、no_university_df["university"] = 0でも作成可能
# university_dfとno_university_dfを、dplyr::bind_rowsを用いて結合してあたらしいdfとし、ID順で並べる
df <- bind_rows(university_df, no_university_df) %>% # 両者を結合
arrange(ID) # ID順に並べる
# 所得変数（income)の作成：上記のDGP(母集団回帰モデル）に基づいて生成
df <- df %>% mutate(income = 200 + 10*ability + 500*university + rnorm(n, mean = 0, sd = 50)) # 誤差項は平均=0、標準偏差=50
# 最初の6行
head(df)## # A tibble: 6 × 4
## ID ability university income
## <int> <dbl> <dbl> <dbl>
## 1 1 89.7 1 1558.
## 2 2 26.6 0 423.
## 3 3 37.2 0 502.
## 4 4 57.3 0 710.
## 5 5 90.8 0 1015.
## 6 6 20.2 0 393.生成したデータを用いて、能力(ability)と所得(income)との関係をプロットすると次のようになる。
# 大卒か否かのラベルをデータフレームに加える
df <- df %>% mutate(edu_label =
case_when(university == 1 ~ "Grad.",
university == 0 ~ "Not grad."))
# 散布図を描く with 大卒ラベル
ggplot(df, aes(x = ability, y = income, color = edu_label)) +
geom_point(alpha = 0.5)+
labs(title = "Ability and Income")
4 差の検定
以下では、上記のように生成したデータを、（データ生成過程を知らない）観察データとみなして分析する。分析目標は、「大学を卒業することが、所得をどれだけ上昇させるのか」という「因果効果」の推定である。
まず、学歴と所得分布の関係を箱ひげ図で表示すると、以下のようになる。
# plot（箱ひげ図）
ggplot(df, aes(x = edu_label, y = income)) +
geom_boxplot(alpha = 0.7) +
labs(title = "Income distribution by university education")
大卒と非大卒の所得の平均には統計学的に有意な差があるかどうかを\(t\)検定によって調べる。
# 大卒
grad_df <- df %>% filter(university == 1)
# 非大卒
not_grad_df <- df %>% filter(university == 0)
# それぞれのincomeを比較
t.test(x = grad_df$income, y = not_grad_df$income)##
## Welch Two Sample t-test
##
## data: grad_df$income and not_grad_df$income
## t = 253.95, df = 5215.2, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## 938.3733 952.9739
## sample estimates:
## mean of x mean of y
## 1600.458 654.784大卒の平均所得が約1600万円、非大卒の平均所得が約655万円であり、1%有意水準でも有意な差があることがわかった。
ただし、ここで「大卒になることで、所得が1600 - 655 = 945万円上昇する」と結論付けることはできない。例えば、学歴と所得の両方に影響を与える交絡要因がある場合、大卒と非大卒の所得差には、このような交絡要因の影響が反映されてしまう。
ここで、自分がこのデータの「真のデータ生成過程」を知っていることを思い起こそう。本データにおいては、所得\(Y\)が上がる要因には学歴\(X\)以外にも能力\(A\)という要因があった。そして、能力\(A\)は所得\(Y\)だけでなく、学歴\(X\)にも影響を与えていた。したがって、945万円という所得差は、学歴（大卒）による所得上昇効果だけでなく、能力による所得上昇効果も混在している。
5 （参考）セレクション・バイアスと処置割り当ての独立性
5.1 セレクション・バイアス
前節で行った方法が潜在的結果の枠組みにおいてどう表されるか示す。 前節の方法は、観測できた各群の結果の平均の差
\[ E[Y(1)|X = 1] - E[Y(0)|X = 0] \]
によってATEを推定する方法である。
これは\(E[Y(0)|X = 1]\)を足して引いて整理すると、次のようになる
\[ \begin{align} E[Y(1)|X = 1] - E[Y(0)|X = 0] &= E[Y(1)|X = 1] - E[Y(0)|X = 1] + E[Y(0)|X = 1] - E[Y(0)|X = 0]\\ &= \underbrace{E[Y(1) - Y(0)|X = 1]}_{\text{ATT}} + \underbrace{E[Y(0)|X = 1] - E[Y(0)|X = 0]}_{\text{Selection Bias}}\\ \end{align} \]
\(\text{ATT}\)は処置群における平均処置効果（average treatment effect on the treated）である。話を単純化するために処置効果が処置群と対照群で平均的に等しい（\(\text{ATT}=\text{ATE}\)）と仮定すると、上記の単純な差による推定量は\(\text{ATE} + \text{Selection Bias}\)であり、セレクション・バイアスの分だけATEから乖離する推定量であることがわかる。
5.2 ランダム化比較試験による推定
もし処置の割当\(X\)が無作為（ランダム）であれば潜在的結果\(Y(0), Y(1)\)とは独立
\[ \newcommand{\indep}{\mathop{\hspace{0.1em} \perp\!\!\!\perp \hspace{0.1em}}} (Y(0), Y(1)) \indep X \]
となる。その場合、
\[ E[Y(0)| X = 0] = E[Y(0)]\\ E[Y(1)| X = 1] = E[Y(1)] \]
になるため、各群の平均の差\(E[Y(1)|X = 1] - E[Y(0)|X = 0]\)でATEを推定することができる。
そのため、もしも分析対象のデータが処置をランダムに割り振った実験によって得られたデータである場合、両群の平均値の差をATEの推定値に使うことは妥当であると考えられる。 しかし、今回のデータは処置の割り当てがランダムではなく条件を満たさないため、この方法を使用することはできない。
5.3 交絡要因の統制による推定
交絡要因（confounder）\(A\)を条件付けた下で、処置の割当\(X\)と結果\(Y(0), Y(1)\)が独立
\[ \newcommand{\indep}{\mathop{\hspace{0.1em} \perp\!\!\!\perp \hspace{0.1em}}} (Y(0), Y(1)) \indep X | A \]
であるという条件を、無交絡性（unconfoundedness）という。無交絡性が成立するとき、平均での独立性
\[ E[Y|X=0, A] = E[Y(0)|X, A] = E[Y(0)|A]\\ E[Y|X=1, A] = E[Y(1)|X, A] = E[Y(1)|A] \]
が成立する。
そこから共変量について期待値をとればATEが推定できる。
\[ E_{A}( E[Y(1) - Y(0)|A] ) = E[Y(1) - Y(0)] \]
条件付き期待値を推定する方法は次節以降で学んでいく。
6 回帰分析
6.1 回帰分析と条件付き平均値
平均値の差の検定では、条件付き期待値\(E(変数|条件)\)の記法を使って表現すると、大卒者の所得の期待値\(E(Y|X=1)\)と非大卒者の所得の期待値\(E(Y|X=0)\)の差
\[ E(Y|X=1) - E(Y|X=0) \]
を推定していた。
しかしこの「差」においては、能力がどんな人であれ、大卒であるかどうかという点のみを条件に比較することになり、それが交絡要因である能力\(A\)による推定バイアスを生じさせていた。
そこで、イメージとしては、
- 能力\(A=80\)の人の中で、大卒(\(X=1\))の所得の条件付き期待値\(E(Y|X = 1,A = 80)\)と非大卒(\(X=0\))の所得の条件付き期待値\(E(Y|X = 0,A = 80)\)を比べる
- 能力\(A=81\)の人の中で、大卒(\(X=1\))の所得の条件付き期待値\(E(Y|X = 1,A = 81)\)と非大卒(\(X=0\))の所得の条件付き期待値\(E(Y|X = 0,A = 81)\)を比べる
- 能力\(A=82\)の人の中で…
- 能力\(A=83\)の人の中で…
といったように、新たに条件を追加して能力\(A\)の値を一定にした下で大卒と非大卒の所得を比較していけばよいと考えられる。
回帰分析の予測値\(\hat{Y}\)は、説明変数を\(X,A\)とすると、これらの説明変数の条件付き期待値\(E(Y|X,A)\)の推定値である。そのため、\(E(Y|X = x,A = a)\)のように、複数の説明変数で条件付けた下での期待値を推定することができる。これを利用して学歴\(X\)が所得\(Y\)に与える効果を推定することができる。
6.2 すべての変数が入手可能な場合
もし、手元に所得\(Y\)、能力\(A\)、学歴\(X\)のデータがあるならば、
\[ Y = \beta_0 + \beta_1 A + \beta_2 X + u \]
という回帰モデルによって係数\(\beta_0,\beta_1,\beta_2\)を推定できる。
OLS_true_model <- lm(formula = income ~ ability + university,
data = df)
stargazer(OLS_true_model, type = "text")##
## =================================================
## Dependent variable:
## -----------------------------
## income
## -------------------------------------------------
## ability 9.951***
## (0.019)
##
## university 501.949***
## (1.857)
##
## Constant 202.780***
## (1.028)
##
## -------------------------------------------------
## Observations 10,000
## R2 0.983
## Adjusted R2 0.983
## Residual Std. Error 49.974 (df = 9997)
## F Statistic 297,084.600*** (df = 2; 9997)
## =================================================
## Note: *p<0.1; **p<0.05; ***p<0.01能力(ability)、学歴（university)、ｙ切片（constant)の係数推定値は、冒頭のDGPで設定したパラメータとほとんど一致する。すなわち、バイアスのない推定値となっている。
ここで、上述したような、「能力\(A=a\)の人の中で、大卒(\(X=1\))の所得の条件付き期待値\(E(Y|X = 1,A = a)\)と非大卒(\(X=0\))の所得の条件付き期待値\(E(Y|X = 0,A = a)\)を比べる」というのは、下図の二つの回帰直線の「縦の差を測る」ことに相当する。一定の条件の下で、ここでの「縦の差」は、学歴が所得に与えた「平均処置効果」の推定値とみなすことができる。
# 散布図と回帰直線を描く with 大卒ラベル
ggplot(df, aes(x = ability, y = income,
color = edu_label,
group = edu_label)) + # groupごと
geom_point(alpha = 0.5)+
geom_smooth(method = "lm", color = "black") + # 回帰直線
labs(title = "Ability and Income")
6.3 入手不可能な変数がある場合
実際には、学歴\(X\)と所得\(Y\)に影響を及ぼす能力\(A\)を正確に観察・測定することは難しい。そこで、能力\(A\)の変数が入手できない場合を想定し、上記の回帰式から\(A\)の項を除いた単回帰モデル
\[ Y = \beta_0 + \beta_1 X + u \]
を推定するとどうなるだろうか。
結果は以下のようになる。
OLS_missing_A <- lm(formula = income ~ university,
data = df)
stargazer(OLS_missing_A, type = "text")##
## ================================================
## Dependent variable:
## ----------------------------
## income
## ------------------------------------------------
## university 945.674***
## (8.580)
##
## Constant 654.784***
## (2.756)
##
## ------------------------------------------------
## Observations 10,000
## R2 0.549
## Adjusted R2 0.548
## Residual Std. Error 261.022 (df = 9998)
## F Statistic 12,147.930*** (df = 1; 9998)
## ================================================
## Note: *p<0.1; **p<0.05; ***p<0.01この推定結果を因果的に解釈すると、「大卒になることによる所得上昇の効果」は、約946（万円）と推定されたことになる（なおこの値は、上記の「差の検定」における大卒と非大卒の平均所得差と等しい）。すなわち、データ生成時に設定した真の値500よりも過大に推定されている。
このように、正しく推定するために必要な変数（能力\(A\)）が考慮されないことによる推定量・推定値の偏りを欠落変数バイアス（omitted variable bias）という。
7 マッチング法
前節において、交絡要因である能力\(A\)の影響を除去したいときには、
- 能力\(A=a\)の人の中で、大卒(\(X=1\))の所得の条件付き期待値\(E(Y|X = 1,A = a)\)と非大卒(\(X=0\))の所得の条件付き期待値\(E(Y|X = 0,A = a)\)を比べる
という作業を行えばよいと述べた。そして、能力\(A\)を説明変数に加えた回帰分析は、このような作業を行っていると解釈できることも説明した。
マッチング法は、このような「条件(\(A\))を揃えた上での処置群(\(X=1\))と対照群(\(X=0\))の間でのアウトカム(\(Y\))の比較」を、いわば、そのまま行う分析手法である。言い換えると、マッチング法とは、共変量ベクトルが同一あるいは近しい個体を、処置群と対照群から「マッチング」させた上で両者のアウトカムの差を計算し、その全体での平均値を平均処置効果の推定値とみなす方法である。
なお、共変量ベクトルと処置変数のデータを用いて「処置を受ける確率」（傾向スコア）を推定した上で、この推定された傾向スコアを用いてマッチングする方法を「傾向スコアマッチング」という。この手法はマッチング法以上に広く使用されてきたが、ここでは扱わない。
マッチング法は、Rでは{Matching}パッケージで実行できる。
マッチングも回帰分析と同じで、必要な共変量が入手できるかどうかによって、推定バイアスの有無や程度に違いがでる。
今回の場合、能力\(A\)が入手できる場合は、これを唯一の共変量として、次のようにマッチング推定を行うことができる。
# 最近傍マッチング（マッチングの方法の１つ）
match_nn <- Matching::Match(Y = df$income, # 被説明変数・結果変数
Tr = df$university, # 処置変数
X = df$ability, # 共変量
version = "fast", # 一部の処理を省略して計算量を減らす設定
Weight = 1) # ユークリッド距離に基づくマッチング
# 結果の表示
summary(match_nn)##
## Estimate... 502.52
## SE......... 2.1998
## T-stat..... 228.44
## p.val...... < 2.22e-16
##
## Original number of observations.............. 10000
## Original number of treated obs............... 1032
## Matched number of observations............... 1032
## Matched number of observations (unweighted). 9874推定値は、設定した500（万円）に非常に近い値となっている。
8 参考文献
計量経済学応用 回帰分析 (矢内勇生氏の講義ノート)
計量経済学応用 傾向スコア (矢内勇生氏の講義ノート)
Chapter 13 - Regression (Nick Huntington-Klein (2021)“The Effect: An Introuction to Research Design and Causality”)
Chapter 14 - Matching (Nick Huntington-Klein (2021)“The Effect: An Introuction to Research Design and Causality”)
なお、今回用いたデータ生成過程（DGP)は、下記のエッセイのものと同一である。
- 安藤道人（2015) 「多重回帰分析と回帰不連続デザイン」 『日本労働研究雑誌』 No.657 pp.12-13