Rで学ぶ計量経済学と機械学習 7
計量経済学３：回帰不連続デザイン（+モンテカルロ・シミュレーション入門）
安藤道人（立教大学） 三田匡能 (株式会社 GA technologies)
2023-05-27
1 本講義の目的
- 回帰モデルを用いて不連続回帰デザインによる推定を行う方法を学ぶ。
- ループの活用やモンテカルロ・シミュレーションの基礎を学ぶ。
2 パッケージの読み込み
今回の分析で使うパッケージを読み込む。
pacman::p_load(tidyverse,
stargazer,
Hmisc,
estimatr,
texreg)3 データ
3.1 データ生成過程
前回、前々回と同様、自分で生成したデータを用いて分析していく。
すなわち、個人の所得\(Y\)と学歴\(X\)・能力\(A\)との関係についての架空データを次のように生成する。
\[ Y = 200 + 10A + 500X+ \varepsilon \]
- サンプルは1万人
- 切片は200万円
- 能力が1上がると所得は10上昇する
- 能力は0から100まで均等に分布する
- 大卒だと所得が500万円上昇する
- 能力を部分的に反映した学力テストの点数が 180点以上であれば大卒となる
- この最後の学力テスト点数が新しいデータ生成条件であり、これを回帰不連続デザインで利用する。
# 事前準備 --------------------
# 乱数の種を固定
set.seed(0)
# データの生成 ----------------
n <- 10000
# 能力は0から100まで均等に分布
ability <- runif(n, min = 0, max = 100)
# IDとabilityをデータフレームに格納する
df <- tibble(ID = 1:n, ability)
# 大卒フラグ
## 条件：能力を部分的に反映した学力テストの点数が 180 点以上であれば大卒となる
## 10%くらいが該当するように恣意的に設定
## abilityとscoreの関係は非線形なものとする
## 今回はmutate()とcase_when()を使用して作成
# 学力テストの点数 score
df <- df %>% mutate(score = 30 * log10(ability) +
rnorm(n, mean = 115, sd = 10))
# 大卒ダミー university (score 180点以上=1)
df <- df %>% mutate(university = case_when(score >= 180 ~ 1,
TRUE ~ 0))
# 所得 income
df <- df %>% mutate(income = 200 + 10*ability + 500*university +
rnorm(n, sd = 50))
## dplyrを使わないやり方だと、以下の通り。
#df["score"] = 30*log10(ability) + rnorm(n, mean = 115, sd = 10)
#df["university"] = 1*(df["score"] >= 180)
#df["income"] = 200 + 10*df["ability"] + 500*df["university"] + rnorm(n, sd = 50)
# 最初の6行
head(df)## # A tibble: 6 × 5
## ID ability score university income
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 89.7 178. 0 1076.
## 2 2 26.6 162. 0 412.
## 3 3 37.2 150. 0 615.
## 4 4 57.3 157. 0 829.
## 5 5 90.8 184. 1 1654.
## 6 6 20.2 136. 0 411.今回はデータを生成する処理を繰り返すことになるので、データの生成処理を関数にまとめる
generate_data <- function(n = 10000) {
# 能力は0から100まで均等に分布
ability <- runif(n, min = 0, max = 100)
# IDとabilityをデータフレームに格納する
df <- tibble(ID = 1:n, ability)
# 大卒フラグ
## 条件：能力を部分的に反映した学力テストの点数が 180 点以上であれば大卒となる
## 10%くらいが該当するように恣意的に設定
## abilityとscoreの関係は非線形なものとする
## 今回はmutate()とcase_when()を使用して作成
df <- df %>% mutate(
# 学力テストの点数 score
score = 30 * log10(ability) + rnorm(n, mean = 115, sd = 10),
# 大卒ダミー university (score 180点以上=1)
university = case_when(score >= 180 ~ 1, TRUE ~ 0),
# 所得 income
income = 200 + 10*ability + 500*university + rnorm(n, sd = 50),
)
return(df)
}seedを指定してからデータを生成すれば同じデータが得られる
set.seed(0)
df2 <- generate_data()
all(df == df2) # すべて同じ値であることを確認## [1] TRUE3.2 グラフと記述統計
こうして生成したデータの所得と能力・学歴・テストの点数の関係をプロットすると次のようになる。
# 塗り分けプロット
## 大卒か否かのラベルをデータフレームに加える
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")
また、データの記述統計は以下のようになる。
Hmisc::describe(df)## df
##
## 6 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.01064336 0.01206559 0.02003808 0.05705222 0.06052661
## highest: 99.91946789 99.92583303 99.94554052 99.98634593 99.99305937
## --------------------------------------------------------------------------------
## score
## n missing distinct Info Mean Gmd .05 .10
## 10000 0 10000 1 161.9 17.93 131.3 141.2
## .25 .50 .75 .90 .95
## 153.3 164.0 173.1 180.5 184.9
##
## lowest : 53.28082 55.15381 58.62861 66.52667 72.85697
## highest: 203.34014 204.18167 205.59547 205.94730 206.05713
## --------------------------------------------------------------------------------
## university
## n missing distinct Info Sum Mean Gmd
## 10000 0 2 0.289 1082 0.1082 0.193
##
## --------------------------------------------------------------------------------
## income
## n missing distinct Info Mean Gmd .05 .10
## 10000 0 10000 1 754.3 423.3 243.1 295.9
## .25 .50 .75 .90 .95
## 450.8 707.5 1007.4 1231.8 1521.0
##
## lowest : 65.29735 89.03678 90.37936 98.69741 101.73090
## highest: 1789.89534 1797.44903 1803.73766 1805.55985 1813.89118
## --------------------------------------------------------------------------------
## edu_label
## n missing distinct
## 10000 0 2
##
## Value Grad. Not grad.
## Frequency 1082 8918
## Proportion 0.108 0.892
## --------------------------------------------------------------------------------4 モンテカルロ・シミュレーション入門
4.1 (欠落変数を無視した)回帰分析
次に、欠落変数（能力\(A\))を無視して、\(Y\)を\(X\)に回帰してみる。前回同様、単純な単回帰に加えて、学力テストの点数をコントロール変数に加えた重回帰分析も行ってみる。
回帰分析結果は以下のようになる。
reg1 = lm(income ~ university,
data = df)
reg2 = lm(income ~ university + score,
data = df)
stargazer(reg1, reg2,
type = "text")##
## ===========================================================================
## Dependent variable:
## -------------------------------------------------------
## income
## (1) (2)
## ---------------------------------------------------------------------------
## university 804.364*** 481.644***
## (9.009) (8.067)
##
## score 12.202***
## (0.151)
##
## Constant 667.226*** -1,273.841***
## (2.963) (24.140)
##
## ---------------------------------------------------------------------------
## Observations 10,000 10,000
## R2 0.444 0.663
## Adjusted R2 0.444 0.663
## Residual Std. Error 279.846 (df = 9998) 217.692 (df = 9997)
## F Statistic 7,971.908*** (df = 1; 9998) 9,849.560*** (df = 2; 9997)
## ===========================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01「大卒になることによる所得上昇の効果」は、単回帰分析(reg1)だと約804（万円）と過大推定される。
scoreを含めた重回帰分析(reg2)だとおよび約482（万円）となる。
4.2 モンテカルロ・シミュレーションとは
scoreを含めた重回帰分析(reg2)における約482（万円）という係数推定値は、500（万円）に比較的近いが、これは、学力テストの点数(score)をコント―ルすることによって欠落変数バイアスが完全に除去できた結果なのだろうか？
このことを、理論的にではなく、基本的なシミュレーションにより検証してみよう。
もし上記の重回帰モデル(reg2)におけるOLS推定量が正しいのであれば、上記の482（万円）と500（万円）の差は「たまたま」生じたにすぎないということになる。
したがって、同様のデータ生成と回帰分析を何百回・何千回と実施すれば、universityの係数のOLS推定値は500(万円)を中心に正規分布するはずである。これを実際に施行してみるのが、モンテカルロ・シミュレーションである。
4.3 モンテカルロ・シミュレーションの実行
下記では、上記のデータ生成と回帰分析(reg2)をループによって1000回繰り返し、その結果を保存する。
# 再現性の確保のため、乱数の種を固定
set.seed(0)
# シミュレーション結果を入れる「箱」の作成
coef <- rep(NA, 1000)
# 同じ回帰分析(上のreg2)を1000回繰りかえす。
for(i in 1:1000){ # forによるループ。1000回繰り返す
# データの生成
df <- generate_data()
# 回帰
reg2 <- lm(income ~ university + score, data = df)
# 回帰係数推定値の保存
coef[i] = reg2$coefficients[2]
}シミュレーション結果（universityの係数推定値の平均、標準偏差、ヒストグラム）は以下のようになる。
#平均、標準誤差
mean(coef)## [1] 464.0726sd(coef)## [1] 7.627223# ヒストグラム
hist(coef)
平均は464.1、標準偏差は7.6であり、上記のreg2の回帰係数値とその標準誤差と近い。
そして、reg2を1000回施行して得たuniversityの係数推計値の平均が464.1であり、500から十分離れていることより、reg2の回帰モデルのuniversityのOLS推定量にもバイアスがあることも分かる。
5 回帰不連続デザイン
ability変数がない限り、通常の回帰分析では推定結果にバイアスが生じてしまう。そこで、「能力を部分的に反映した学力テストの点数が
180 点以上であれば大卒となる」という条件を利用した分析を考える。
まず冒頭のデータをもう一度生成させる。
# データ生成
set.seed(0)
df <- generate_data()
head(df)## # A tibble: 6 × 5
## ID ability score university income
## <int> <dbl> <dbl> <dbl> <dbl>
## 1 1 89.7 178. 0 1076.
## 2 2 26.6 162. 0 412.
## 3 3 37.2 150. 0 615.
## 4 4 57.3 157. 0 829.
## 5 5 90.8 184. 1 1654.
## 6 6 20.2 136. 0 411.X軸を学力テストの点数(score)、Y軸を所得(income)として散布図を描くと、下記のようになる。
# 塗り分けプロット
## 大卒か否かのラベルをデータフレームに加える
df <- df %>% mutate(edu_label =
case_when(university == 1 ~ "Grad.",
university == 0 ~ "Not grad."))
## 散布図を描く with 大卒ラベル
ggplot(df, aes(x = score,
y = income,
color = edu_label)) +
geom_point(alpha = 0.5) +
labs(title = "Score and Income")
テストの点数(score)の閾値（180点）を境に、所得が全体的に上振れ（ジャンプ）している。これは、学力テストで180点を大卒とおなり、所得が上昇するからである。
回帰不連続デザイン(regression discontinuity design: RD design, RDD)はこうした閾値を伴うデータ生成過程(DGP)を利用する調査設計（デザイン）で、以下のようなロジックをもとに分析を行う。
- スコアが180点になる閾値で大卒になれるかどうかという「不連続な変化」を決めるというのは人為的な制度であり、その背景にある個々人の所得の決定要因（ここでは能力\(A\)のみだが、一般的にはあらゆる決定要因）は「連続的」に分布しているはず。
- したがって、この閾値において所得に明確なジャンプ（不連続性）がある理由は、この閾値で同じく不連続に変化している「大卒」という要因以外には考えられない
5.1 RDデザインによる推定
RDデザインの推定は、より高度なRD推定方法を扱うパッケージが存在するが、lm関数を使った重回帰分析で推定できる。
基本的には閾値（ここではscore=180)での「ジャンプ」の大きさを推定すればよい。
もっとも簡単な方法は、線形回帰分析でこのジャンプを推定することである。まずは簡便に、ggplot2を使って、グラフを使ってこのジャンプを確認してみよう。
## 散布図を描く with 大卒ラベル
ggplot(df, aes(x = score,
y = income,
color = edu_label,
group = edu_label)) +
geom_point(alpha = 0.5) +
labs(title = "Score and Income") +
geom_smooth(method = lm,
se = TRUE,
color = "black") # group=edu_labelごとの回帰直線, SEあり
上図のようなRD分析を線形RD(linear RD)と呼び、以下のように推定できる。
# 閾値=0となるように割当変数(score_gap)を作り直す
df <- df %>% mutate(score_gap = score - 180)
# 線形RD: X軸にscore_gapを用いることにより、y切片の変化すなわちuniversityダミーの係数を上図の「ジャンプ」の大きさと解釈できる。
linear_RD <- lm(income ~ university + score_gap + score_gap:university,
data = df)
stargazer(linear_RD,
type = "text")##
## ================================================
## Dependent variable:
## ---------------------------
## income
## ------------------------------------------------
## university 527.334***
## (10.890)
##
## score_gap 12.304***
## (0.152)
##
## university:score_gap -8.756***
## (1.405)
##
## Constant 924.648***
## (3.919)
##
## ------------------------------------------------
## Observations 10,000
## R2 0.665
## Adjusted R2 0.665
## Residual Std. Error 217.281 (df = 9996)
## F Statistic 6,604.168*** (df = 3; 9996)
## ================================================
## Note: *p<0.1; **p<0.05; ***p<0.01universityの係数は閾値でのジャンプの大きさ、scoreの係数は閾値左側の回帰直線の傾き、score:universityの係数は閾値左側と閾値右側の回帰直線の傾きの変化である。universityの係数推定値は500に近い。
5.2 多項式RD
閾値でのジャンプの大きさを測りたいだけなのに、閾値の左右におけるscoreとincomeの関係を線形回帰モデルで捉えるのは適切ではないかもしれない。
ただし、loess曲線を用いて閾値の左右におけるscoreとincomeの関係をフレキシブルな曲線で捉えてみても、ジャンプの大きさはあまり変わらないように見える。
## 散布図を描く with 大卒ラベル
ggplot(df, aes(x = score,
y = income,
color = edu_label,
group = edu_label)) +
geom_point(alpha = 0.5)+
labs(title = "Score and Income") +
geom_smooth(method = lm,
se = TRUE,
color = "black") + # group=edu_labelごとの回帰直線, SEあり
geom_smooth(method = loess,
se = TRUE,
color = "black") # group=edu_labelごとのloess曲線, SEあり
そこで、閾値の左右におけるscoreとincomeの関係を、それぞれ異なる二次、三次、四次曲線でとらえた多項式RDを推定してみる。
# 二次多項式RD, I(score_gap^2)はscore_gap^2
secondpoly_RD <- lm(income ~ university +
score_gap +
I(score_gap^2) +
score_gap:university +
I(score_gap^2):university,
data = df)
# 三次多項式RD, I(score_gap^3)はscore_gap^3
thirdpoly_RD <- lm(income ~ university +
score_gap +
I(score_gap^2) +
I(score_gap^3) +
score_gap:university +
I(score_gap^2):university +
I(score_gap^3):university,
data = df)
# 四次多項式RD, I(score_gap^4)はscore_gap^4
fourthpoly_RD <- lm(income ~ university +
score_gap +
I(score_gap^2) +
I(score_gap^3) +
I(score_gap^4) +
score_gap:university +
I(score_gap^2):university +
I(score_gap^3):university +
I(score_gap^4):university ,
data = df)
# 一次~四次多項式RD
stargazer(linear_RD,
secondpoly_RD,
thirdpoly_RD,
fourthpoly_RD,
type = "text",
digits = 1,
df = FALSE) # 幅を狭くするため、digits=1, df=FALSEとする##
## ====================================================================
## Dependent variable:
## -------------------------------------------
## income
## (1) (2) (3) (4)
## --------------------------------------------------------------------
## university 527.3*** 461.0*** 488.3*** 523.1***
## (10.9) (14.1) (17.2) (20.6)
##
## score_gap 12.3*** 18.5*** 14.9*** 6.4***
## (0.2) (0.4) (0.7) (1.2)
##
## I(score_gap2) 0.1*** -0.02 -0.5***
## (0.01) (0.02) (0.1)
##
## I(score_gap3) -0.001*** -0.01***
## (0.000) (0.001)
##
## I(score_gap4) -0.000***
## (0.000)
##
## university:score_gap -8.8*** -14.3*** -13.1* -0.4
## (1.4) (3.7) (7.7) (14.0)
##
## university:I(score_gap2) -0.1 0.3 -0.2
## (0.2) (0.9) (2.8)
##
## university:I(score_gap3) -0.01 0.1
## (0.03) (0.2)
##
## university:I(score_gap4) -0.002
## (0.004)
##
## Constant 924.6*** 989.5*** 965.5*** 926.9***
## (3.9) (5.3) (6.7) (8.1)
##
## --------------------------------------------------------------------
## Observations 10,000 10,000 10,000 10,000
## R2 0.7 0.7 0.7 0.7
## Adjusted R2 0.7 0.7 0.7 0.7
## Residual Std. Error 217.3 213.9 213.5 212.8
## F Statistic 6,604.2*** 4,152.7*** 2,981.1*** 2,341.7***
## ====================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01どの結果も、universityの係数推定値は、だいたい500前後になる。ただし、多項式の次元を大きくしていくほど、標準誤差も大きくなる。
5.3 局所線形RD
もう一つの考え方として、閾値から遠ざかるほど、閾値でのジャンプの大きさを測るには不必要な情報となる。したがって、線形モデルでのRDを使いつつ、推定に用いるサンプルを閾値近傍の局所的（local)なものに限定すればよいのではないかと考える。
この考え方に従って、閾値前後のバンド幅(bandwidth)を設定し、そのバンド幅の内部のサンプルだけを使って線形RDを行う。
# バンド幅=30
band30_df <- df %>% filter(-30 <= score_gap & score_gap < 30)
band30_RD <- lm(income ~ university + score_gap + score_gap:university,
data = band30_df)
# バンド幅=10
band10_df <- df %>% filter(-10 <= score_gap & score_gap < 10)
band10_RD <- lm(income ~ university + score_gap + score_gap:university,
data = band10_df)
# バンド幅=5
band5_df <- df %>% filter(-5 <= score_gap & score_gap < 5)
band5_RD <- lm(income ~ university + score_gap + score_gap:university,
data = band5_df)
# バンド幅=1
band1_df <- df %>% filter(-1 <= score_gap & score_gap < 1)
band1_RD <- lm(income ~ university + score_gap + score_gap:university,
data = band1_df)
stargazer(band30_RD,
band10_RD,
band5_RD,
band1_RD,
type = "text",
digits = 1,
df = FALSE)##
## ============================================================
## Dependent variable:
## ---------------------------------------
## income
## (1) (2) (3) (4)
## ------------------------------------------------------------
## university 494.3*** 510.8*** 501.2*** 569.3***
## (11.9) (15.2) (19.7) (40.8)
##
## score_gap 14.1*** 12.3*** 15.0*** -27.4
## (0.3) (1.5) (4.4) (51.3)
##
## university:score_gap -10.6*** -8.7*** -8.9 -58.9
## (1.5) (2.9) (7.2) (70.9)
##
## Constant 957.7*** 940.7*** 946.6*** 920.9***
## (5.6) (9.5) (13.3) (29.8)
##
## ------------------------------------------------------------
## Observations 8,071 3,128 1,534 327
## R2 0.6 0.6 0.7 0.7
## Adjusted R2 0.6 0.6 0.7 0.7
## Residual Std. Error 224.7 206.1 195.5 186.3
## F Statistic 3,803.0*** 1,794.2*** 972.6*** 208.3***
## ============================================================
## Note: *p<0.1; **p<0.05; ***p<0.01バンド幅が10や5の場合はuniversity
の係数推定値は500に近いが、バンド幅を1にまで狭めると、推定値はかなり大きくなってしまっている。
また、標準誤差もバンド幅を狭めるごとに大きくなる。
一般的に、RD推定は、閾値前後のサンプルサイズが十分に大きくないと推定結果は不安定になりがちである。一方で、サンプルサイズを確保するために閾値を大きくするとバイアスが生じる可能性が高くなり、閾値を大きくすることによるバイアスを高次多項式によってコントロールしようとすると多項式により推定結果は不安定になる。
ためしに以下では、サンプルサイズを100倍すなわちn=1,000,000にして局所線形RDを推定してみる。
# データ生成
set.seed(0)
df <- generate_data(n = 1000000)
#割当変数score_gapの作成
df <- df %>% mutate(score_gap = score - 180)
# バンド幅=30
band30_df <- df %>% filter(-30 <= score_gap & score_gap < 30)
band30_RD <- lm(income ~ university + score_gap + score_gap:university,
data = band30_df)
# バンド幅=10
band10_df <- df %>% filter(-10 <= score_gap & score_gap < 10)
band10_RD <- lm(income ~ university + score_gap + score_gap:university,
data = band10_df)
# バンド幅=5
band5_df <- df %>% filter(-5 <= score_gap & score_gap < 5)
band5_RD <- lm(income ~ university + score_gap + score_gap:university,
data = band5_df)
# バンド幅=1
band1_df <- df %>% filter(-1 <= score_gap & score_gap < 1)
band1_RD <- lm(income ~ university + score_gap + score_gap:university,
data = band1_df)
stargazer(band30_RD,
band10_RD,
band5_RD,
band1_RD,
type = "text",
digits = 1,
df = FALSE)##
## ======================================================================
## Dependent variable:
## -------------------------------------------------
## income
## (1) (2) (3) (4)
## ----------------------------------------------------------------------
## university 469.0*** 495.6*** 498.1*** 502.1***
## (1.2) (1.5) (2.0) (4.5)
##
## score_gap 14.1*** 9.9*** 9.1*** -0.9
## (0.03) (0.2) (0.4) (5.4)
##
## university:score_gap -8.6*** -3.7*** -2.9*** 6.5
## (0.1) (0.3) (0.7) (7.8)
##
## Constant 959.7*** 930.0*** 927.5*** 923.7***
## (0.6) (1.0) (1.4) (3.1)
##
## ----------------------------------------------------------------------
## Observations 804,974 313,249 155,905 31,170
## R2 0.6 0.6 0.6 0.6
## Adjusted R2 0.6 0.6 0.6 0.6
## Residual Std. Error 223.1 206.7 200.2 197.1
## F Statistic 366,465.4*** 166,014.0*** 88,226.4*** 16,988.4***
## ======================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01バンド幅30だと少し下方バイアスがあるように見えるが、バンド幅を10以下にすると、ほぼ500になる。
5.4 局所多項式RD
サンプルを閾値近傍のものに限定した局所的な分析と多項式（高次元）回帰モデルを組み合わせることもできる。ただし、局所回帰と高次元回帰モデルを組み合わせると標準誤差が大きく傾向があるので、ここでは二次多項式回帰のみを用いてみる。
# データ生成
set.seed(0)
df <- generate_data(n = 10000)
#割当変数score_gapの作成
df <- df %>% mutate(score_gap = score - 180)
# バンド幅=30
band30_df <- df %>% filter(-30 <= score_gap & score_gap < 30)
band30_RD <- lm(income ~ university + score_gap + I(score_gap^2) +
score_gap:university + I(score_gap^2):university,
data = band30_df)
# バンド幅=10
band10_df <- df %>% filter(-10 <= score_gap & score_gap < 10)
band10_RD <- lm(income ~ university + score_gap + I(score_gap^2) +
score_gap:university + I(score_gap^2):university,
data = band10_df)
# バンド幅=5
band5_df <- df %>% filter(-5 <= score_gap & score_gap < 5)
band5_RD <- lm(income ~ university + score_gap + I(score_gap^2) +
score_gap:university + I(score_gap^2):university,
data = band5_df)
# バンド幅=1
band1_df <- df %>% filter(-1 <= score_gap & score_gap < 1)
band1_RD <- lm(income ~ university + score_gap + I(score_gap^2) +
score_gap:university + I(score_gap^2):university,
data = band1_df)
stargazer(band30_RD,
band10_RD,
band5_RD,
band1_RD,
type = "text",
digits = 1,
df = FALSE)##
## ================================================================
## Dependent variable:
## ---------------------------------------
## income
## (1) (2) (3) (4)
## ----------------------------------------------------------------
## university 519.3*** 510.8*** 527.5*** 556.4***
## (16.3) (22.2) (29.2) (61.5)
##
## score_gap 9.2*** 13.0** -2.0 142.9
## (1.3) (6.3) (18.0) (207.0)
##
## I(score_gap2) -0.2*** 0.1 -3.3 167.7
## (0.04) (0.6) (3.4) (197.4)
##
## university:score_gap -5.0 -10.2 -7.1 -337.4
## (4.1) (11.1) (28.0) (286.0)
##
## university:I(score_gap2) 0.1 0.04 6.6 -57.7
## (0.2) (1.1) (5.6) (276.9)
##
## Constant 931.2*** 941.9*** 931.5*** 950.7***
## (8.8) (14.7) (20.3) (46.0)
##
## ----------------------------------------------------------------
## Observations 8,071 3,128 1,534 327
## R2 0.6 0.6 0.7 0.7
## Adjusted R2 0.6 0.6 0.7 0.7
## Residual Std. Error 224.5 206.1 195.5 186.5
## F Statistic 2,288.5*** 1,075.9*** 583.7*** 124.8***
## ================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01全体としてRD推定値(universityの係数推定値）は500よりも大きくなり、とりわけもっともバンド幅が小さいときには500から大きく離れている。
ここでも、サンプルサイズを100倍(100万人)にするとどうなるか見てみよう。
# データ生成
set.seed(0)
df <- generate_data(n = 1000000)
#割当変数score_gapの作成
df <- df %>% mutate(score_gap = score - 180)
# バンド幅=50
band30_df <- df %>% filter(-30 <= score_gap & score_gap < 30)
band30_RD <- lm(income ~ university + score_gap + I(score_gap^2) +
score_gap:university + I(score_gap^2):university,
data = band30_df)
# バンド幅=10
band10_df <- df %>% filter(-10 <= score_gap & score_gap < 10)
band10_RD <- lm(income ~ university + score_gap + I(score_gap^2) +
score_gap:university + I(score_gap^2):university,
data = band10_df)
# バンド幅=5
band5_df <- df %>% filter(-5 <= score_gap & score_gap < 5)
band5_RD <- lm(income ~ university + score_gap + I(score_gap^2) +
score_gap:university + I(score_gap^2):university,
data = band5_df)
# バンド幅=1
band1_df <- df %>% filter(-1 <= score_gap & score_gap < 1)
band1_RD <- lm(income ~ university + score_gap + I(score_gap^2) +
score_gap:university + I(score_gap^2):university,
data = band1_df)
stargazer(band30_RD,
band10_RD,
band5_RD,
band1_RD,
type = "text",
digits = 1,
df = FALSE)##
## =========================================================================
## Dependent variable:
## ------------------------------------------------
## income
## (1) (2) (3) (4)
## -------------------------------------------------------------------------
## university 496.5*** 498.7*** 496.6*** 506.0***
## (1.6) (2.3) (3.1) (6.7)
##
## score_gap 8.0*** 7.9*** 11.6*** -15.4
## (0.1) (0.6) (1.8) (21.5)
##
## I(score_gap2) -0.2*** -0.2*** 0.5 -14.5
## (0.004) (0.1) (0.3) (20.7)
##
## university:score_gap -0.7 -1.2 -6.4** 12.1
## (0.4) (1.1) (2.9) (31.0)
##
## university:I(score_gap2) 0.1*** 0.1 -0.3 23.5
## (0.02) (0.1) (0.6) (30.1)
##
## Constant 927.4*** 926.3*** 929.7*** 921.3***
## (0.9) (1.5) (2.1) (4.7)
##
## -------------------------------------------------------------------------
## Observations 804,974 313,249 155,905 31,170
## R2 0.6 0.6 0.6 0.6
## Adjusted R2 0.6 0.6 0.6 0.6
## Residual Std. Error 222.8 206.7 200.2 197.1
## F Statistic 220,962.3*** 99,613.5*** 52,936.3*** 10,192.8***
## =========================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01どの推定値も500にかなり近くなる。
6 補足
上記の一連の分析からもわかるように、RD推定値は、設定としてRDデザインの適用がふさわしい場合においても、サンプルサイズ、バンド幅、多項式回帰モデルの設定などによって結果が大きくかわりうる。
RD推定の最適な方法については様々な専門的議論がなされ、Rのパッケージも提供されている（{rdrobust}パッケージなど）。
ただし上記の方法で、lm関数を使って基本的なRD推定を行うことが可能であるし、この作業によって、何をどう分析しているのかを理解することは重要である。
また標準誤差についても、ここでは修正しなかったが、ロバスト標準誤差に置き換えるほうが望ましい。
7 参考文献
安藤道人（2015) 「多重回帰分析と回帰不連続デザイン」 『日本労働研究雑誌』 No.657 pp.12-13
安藤道人「計量経済学２」あるいは「計量経済特論２」の「回帰不連続デザイン」の講義資料
計量経済学応用 回帰不連続デザイン (矢内勇生氏の講義ノート)
Chapter 20 - Regression Discontinuity (Nick Huntington-Klein (2021)“The Effect: An Introuction to Research Design and Causality”)
RD Packages (より学術的に洗練された方法で不連続回帰を実施するためのRパッケージ)