Rで学ぶ計量経済学と機械学習 6
計量経済学２：操作変数法
安藤道人（立教大学） 三田匡能 (株式会社 GA technologies)
2023-10-28
1 本講義の目的
- 操作変数法の実施方法を学ぶ。
3 データの用意
前回と同じく、架空のデータを生成して分析を行う。
前回は個人の所得\(Y\)と学歴\(X\)・能力\(A\)との関係についての架空データを次のように生成した（ただし\(\varepsilon\)は誤差項）。
\[ Y = 200 + 10A + 500X + \varepsilon \]
- サンプルサイズは10,000人
- 能力は0から100まで均等に分布する
- 能力が1上がると所得は10万円上昇する
- 能力が80 以上の約 2,000人中から約 1,000人がランダムに選ばれて大卒となる。
- 大卒だと所得が500万円上昇する
- 切片は200万円（能力=0、非大卒の平均給与）
今回は、前回のデータに新たな変数を加える。「能力に関係なく、30％の確率でランダムに選ばれた人が学費を免除され、入学を許可される」という制度があるとし、さらに、「学費免除に選ばれた人は、50%の確率で大学を卒業できる」（ただし能力が80以上の条件で選ばれて大卒となった人はこの法則の対象外）とする。
データ生成過程は、以下のように設定できる。（難しく感じる人は、コピペでかまわない）
# 事前準備 --------------------
# 乱数の種を固定 => 毎回同じように乱数を発生させるようにする。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 <- df %>%
dplyr::anti_join(university_df,
by = c("ID","ability")) # 大卒ではない人
# no_university_dfのデータフレームに、universityという変数を作成し、すべて0とする。
no_university_df["university"] = 0
# university_dfとno_university_dfを、dplyr::bind_rowsを用いて結合してあたらしいdfとし、ID順で並べる
df_temp1 <- bind_rows(university_df, no_university_df) %>% # 両者を結合
arrange(ID) # ID順にする
# ここまで前回と同じ ------------
# ここから追加部分 --------------
# 条件３：学費免除制度とそれによる卒業
# 前回はsample_frac()を用いてランダムに大卒サンプルを抽出し、また結合するという方法を用いた。
# 今回は二項分布の乱数を発生させるrbinom()を用いて、
# 元のデータフレームに学費免除ダミー変数を追加したり、大卒ダミー変数を修正したりする。
# df_tempの個体は、能力に関係なく、30%の確率で学費免除（exemption = 1）となる
# rbinom(n, m, p): 成功確率pの試行をm回試みたときの成功回数の二項分布からn個の乱数を生成
df_temp2 <- df_temp1 %>%
dplyr::mutate(exemption = rbinom(n, 1, 0.3))
# 学費免除（exemption=1）のときには、50%の確率で卒業できる（university=1）
# ただし、能力が80以上のときには学費免除の有無に関わりなく卒業できるものとする。
# mutate()の中でcase_when()を使う
# case_when()の場合分けでの数字生成がnumericで行われるように、as.numeric()の中でrbinomを使って0,1の乱数を生成
df <- df_temp2 %>% mutate(
university = case_when(
exemption == 1 ~ as.numeric(rbinom(n, 1, 0.5)), # exemption=1なら50%の確率で1
university == 1 ~ 1, #すでにuniversity=1なら1
TRUE ~ 0 # それ以外は0
))
# 所得の生成
df <- df %>% mutate(
income =
200 +
10*ability +
500*university +
rnorm(n, mean = 0, sd = 50) # 誤差項は平均=0、SD=50
)
# 伝統的なやり方はこちら→df["income"] = 200 + 10*df["ability"] + 500*df["university"] + rnorm(n, mean = 0, sd = 50) # 誤差項は平均=0、SD=50
# 最初の6行
head(df)## # A tibble: 6 × 5
## ID ability university exemption income
## <int> <dbl> <dbl> <int> <dbl>
## 1 1 89.7 1 0 1680.
## 2 2 26.6 0 1 544.
## 3 3 37.2 0 0 595.
## 4 4 57.3 0 0 730.
## 5 5 90.8 0 0 1114.
## 6 6 20.2 1 1 905.生成したデータの記述統計をHmisc::describeで表示し、どのくらいが学費免除および大卒になったかチェックしてみる。
## df
##
## 5 Variables 10000 Observations
## --------------------------------------------------------------------------------
## ID
## n missing distinct Info Mean Gmd .05 .10
## 10000 0 10000 1 5000 3334 501 1001
## .25 .50 .75 .90 .95
## 2501 5000 7500 9000 9500
##
## lowest : 1 2 3 4 5, highest: 9996 9997 9998 9999 10000
## --------------------------------------------------------------------------------
## ability
## n missing distinct Info Mean Gmd .05 .10
## 10000 0 10000 1 50.03 33.6 4.850 9.687
## .25 .50 .75 .90 .95
## 25.321 49.656 75.644 90.397 95.246
##
## lowest : 0.0106434 0.0120656 0.0200381 0.0570522 0.0605266
## highest: 99.9195 99.9258 99.9455 99.9863 99.9931
## --------------------------------------------------------------------------------
## university
## n missing distinct Info Sum Mean Gmd
## 10000 0 2 0.517 2215 0.2215 0.3449
##
## --------------------------------------------------------------------------------
## exemption
## n missing distinct Info Sum Mean Gmd
## 10000 0 2 0.63 2997 0.2997 0.4198
##
## --------------------------------------------------------------------------------
## income
## n missing distinct Info Mean Gmd .05 .10
## 10000 0 10000 1 810.8 452.7 254.8 314.1
## .25 .50 .75 .90 .95
## 491.0 770.4 1050.5 1484.5 1605.4
##
## lowest : 56.4974 68.0837 74.5511 82.2678 90.9306
## highest: 1782.56 1792.01 1800.88 1805.26 1823.69
## --------------------------------------------------------------------------------こうして生成したデータの所得と能力・学歴・テストの点数の関係をプロットすると次のようになる。前回と異なり、abilityが80を下回っても、学費補助によって大学を卒業し、所得が上昇している人々がいることが分かる。
# 大卒か否かのラベルをデータフレームに加える
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 操作変数法の概要
効果を測りたい処置変数（学歴\(X\)）を通じてのみアウトカム変数（所得\(Y\))に影響を与える変数のことを操作変数（instrumental variable）という。そのような変数が存在することは多くないが、今回の場合は学費免除ダミー\(Z\)を操作変数として用いることができる。
前回、回帰分析などで学歴の効果を推定しようとした際には、交絡要因である能力がデータとして得られないという想定の下では 学歴の効果が過大に推定される問題が起こることを述べた。操作変数法（instrumental variable method、IV method）は操作変数から得られる情報を活用することで、そのような状況においてもよりよい推定を行う方法である。
4.1 操作変数法の求め方
ベーシックな操作変数推定値は、以下の２種類のやり方によって求めることができる（両者の推定値は一致する）。
- 二段階最小二乗法(TSLS)：学歴\(X\)の変動のうち、学費免除\(Z\)の変動による変動分\(\hat{X}\)を推計し（第一段階）、\(\hat{X}\)→所得\(Y\)の因果効果を推定する（第二段階）。
- 誘導型推定値を第一段階推定値で割る：学費免除→学歴の因果効果（第一段階）、学費免除→所得の因果効果（誘導型）はそれぞれ回帰分析で推定できるので、この二つの推定値を利用して、以下のように推定する
\[ 学歴X\to 所得Y \text{の因果効果} = \frac{学費免除Z \to 所得Y\text{の因果効果}}{学費免除Z \to 学歴X\text{の因果効果}} \]
実際には、第一段階の推定結果を開示しつも、操作変数推定値およびその標準誤差は二段階最小二乗法(TSLS)によって求めることがほとんどである。
なお二段階最小二乗法を用いる場合は、手動で二段階の推定するのではなく、下記で提示するパッケージなどを使うこと。自ら\(\hat{X}\)を推定し、それを\(Y\)に回帰して推定値と標準誤差を得ると、誤った標準誤差を算出することになる。
5 Rで実践
5.1 欠落変数を無視した回帰分析
まず、欠落変数（能力\(A\))を無視して、\(Y\)を\(X\)に回帰してみる。単純な単回帰に加えて、学費補助\(Z\)をコントロール変数に加えた重回帰分析も行ってみる。
回帰分析結果は以下のようになる。
reg_biased1 <- lm(formula = income ~ university,
data = df)
reg_biased2 <- lm(formula = income ~ university + exemption,
data = df)
stargazer(reg_biased1, reg_biased2,
type = "text")##
## ===========================================================================
## Dependent variable:
## -------------------------------------------------------
## income
## (1) (2)
## ---------------------------------------------------------------------------
## university 674.704*** 713.638***
## (6.889) (7.606)
##
## exemption -80.775***
## (6.895)
##
## Constant 661.336*** 676.921***
## (3.242) (3.484)
##
## ---------------------------------------------------------------------------
## Observations 10,000 10,000
## R2 0.490 0.497
## Adjusted R2 0.490 0.496
## Residual Std. Error 286.052 (df = 9998) 284.123 (df = 9997)
## F Statistic 9,593.290*** (df = 1; 9998) 4,930.641*** (df = 2; 9997)
## ===========================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01「大卒になることによる所得上昇の効果」は、約675（万円）および約714（万円）と過大推定される。
5.2 誘導型推定値を第一段階推定値で割る
二段階最小二乗法を実施する前に、操作変数法をより深く理解するために、第一段階推定と誘導型推定を単純にOLSで実行し、両者の推定値より操作変数推定値を求めてみよう。
それぞれ、以下の回帰式で推定できる。
- 第一段階（first stage）：\(X_i =\pi_0 + \pi_1 Z_i + v_i\)
- 誘導型（reduced form）：\(Y_i = \gamma_0 + \gamma_1 Z_i + w_i\)
求めたい因果効果は
\[ 学歴X\to 所得Y \text{の因果効果} = \frac{\gamma_1}{\pi_1} = \frac{\text{誘導型の}Z_i \text{の係数}}{\text{第一段階の}Z_i\text{の係数}} \]
となる。
第一段階の推定結果は以下のようになる。
first_stage <- lm(formula = university ~ exemption,
data = df)
stargazer(first_stage, type = "text")##
## ===============================================
## Dependent variable:
## ---------------------------
## university
## -----------------------------------------------
## exemption 0.396***
## (0.008)
##
## Constant 0.103***
## (0.004)
##
## -----------------------------------------------
## Observations 10,000
## R2 0.191
## Adjusted R2 0.191
## Residual Std. Error 0.374 (df = 9998)
## F Statistic 2,358.704*** (df = 1; 9998)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01学歴ダミーの係数推定値は、約0.4である。本モデルは線形確率モデルであるため、この結果は、学費免除になると大卒になる確率が約40%上昇することを示している。
なお、データ生成過程では学費免除者の約50%(10000×0.3×0.5で約1500人）が大卒になると設定したが、なぜここでは40%と推定されるのだろうか。
それは、学費免除者の中でも、能力が80以上のものの半分（3000×0.2×0.5で約300人）は学費免除とは関係なく大卒となるため、学費免除の「因果効果」は、1500/3000=0.5ではなく、(1500-300)/3000=0.4となるからである。第一段階推定値はこの値を正確に推定している。
一方、誘導型の推定結果は以下のようになる。
##
## ===============================================
## Dependent variable:
## ---------------------------
## income
## -----------------------------------------------
## exemption 201.839***
## (8.504)
##
## Constant 750.292***
## (4.656)
##
## -----------------------------------------------
## Observations 10,000
## R2 0.053
## Adjusted R2 0.053
## Residual Std. Error 389.599 (df = 9998)
## F Statistic 563.308*** (df = 1; 9998)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01学歴ダミーの係数推定値は、約202であり、これは、「学費免除にあたると（大卒になることを通じて）平均所得が約202万円上昇する」ことを意味している。つまり誘導形推定値は、「学費免除の所得に対する平均処置効果」である。
最後に、両者の係数を割ると、操作変数推定値を得ることができる。
## exemption
## 509.6706「学歴（大卒）が所得に与える効果」である500万を正確に推定できている。ただし、このような方法では、標準誤差を得ることはできない。
5.3 二段階最小二乗法(TSLS)を実行する
{AER}パッケージに操作変数法による回帰を行うivreg()という関数がある。これを実行すれば操作変数法で推定し、標準誤差を出すことができる。
iv_reg = AER::ivreg(
formula = income ~ university | exemption, # 被説明変数 ~ 説明変数 | 操作変数 と指定
data = df
)
summary(iv_reg)##
## Call:
## AER::ivreg(formula = income ~ university | exemption, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -641.394 -249.976 -2.886 256.761 625.671
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 697.891 4.642 150.34 <2e-16 ***
## university 509.671 16.213 31.44 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 294.1 on 9998 degrees of freedom
## Multiple R-Squared: 0.4604, Adjusted R-squared: 0.4603
## Wald test: 988.2 on 1 and 9998 DF, p-value: < 2.2e-16##
## ==========================
## Model 1
## --------------------------
## (Intercept) 697.891 ***
## (4.642)
## university 509.671 ***
## (16.213)
## --------------------------
## R^2 0.460
## Adj. R^2 0.460
## Num. obs. 10000
## ==========================
## *** p < 0.001; ** p < 0.01; * p < 0.05TSLS推定値は、約510であり、先に求めた値とも一致する。
5.4 ロバスト標準誤差
ただし、TSLSにおいても、一般的には、不均一分散に対して頑健な標準誤差（ロバスト標準誤差）を使うことが望ましい。そこで、{estimatr}パッケージのiv_robust()を使えば、ロバスト標準誤差を用いてTSLSを実施できる。実際にデータ分析する際には、この関数を活用することが望ましい。
iv_reg_robust <- estimatr::iv_robust(
income ~ university | exemption, # 被説明変数 ~ 説明変数 | 操作変数
data = df,
se_type ="HC1"
)
# summary
summary(iv_reg_robust)##
## Call:
## estimatr::iv_robust(formula = income ~ university | exemption,
## data = df, se_type = "HC1")
##
## Standard error type: HC1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|) CI Lower CI Upper DF
## (Intercept) 697.9 4.648 150.16 0.000e+00 688.8 707.0 9998
## university 509.7 16.176 31.51 7.917e-208 478.0 541.4 9998
##
## Multiple R-squared: 0.4604 , Adjusted R-squared: 0.4603
## F-statistic: 992.7 on 1 and 9998 DF, p-value: < 2.2e-16##
## =========================================
## Model 1 Model 2
## -----------------------------------------
## (Intercept) 697.891 *** 697.891 ***
## (4.642) (4.648)
## university 509.671 *** 509.671 ***
## (16.213) (16.176)
## -----------------------------------------
## R^2 0.460 0.460
## Adj. R^2 0.460 0.460
## Num. obs. 10000 10000
## RMSE 294.149
## =========================================
## *** p < 0.001; ** p < 0.01; * p < 0.05htmlで出力するには、{textreg}パッケージのhtmlreg()を使う。
6 参考文献
計量経済学応用 操作変数法 (矢内勇生氏の講義ノート)
Chapter 19 - Instrumental Variables (Nick Huntington-Klein (2021)“The Effect: An Introuction to Research Design and Causality”)