27. 重回帰分析 10: RDD 基礎

library(broom)     # For converting models into tables
library(DT)
library(haven)
library(huxtable)  # For side-by-side regression table
library(rdrobust)  # For robust nonparametric RDD  
library(rddensity) # For nonparametric RDD
library(stargazer)
library(tidyverse) # For ggplot

• 作図用代用フォントの設定

theme_set(theme_gray(base_size = 10, base_family = "HiraginoSans-W3"))  # macOS用

theme_set(theme_gray(base_size = 10, base_family = "Meiryo"))        # Windows用

I. RDD を使ってわかること

1. RDD の定義

• 不連続回帰デザイン (RDD:regression discontinuity design) は1960年に Thistlewaite and Campbell によって提唱された研究手法
  

• 研究手法としては古くから存在しているが、実際に研究デザインとして使われ始めたのは 1999年以降
  

• Angrist and Lavy (1999)・・・小学校のクラスの大きさが成績に与える影響
  

• Lee (2008)・・・米国議会選挙における現職政治家であることの影響
  

• 有益な Rパッケージ:
  ・rdrobust・・・バンド幅の選択や処置効果の推定に関して最も進歩した Rパッケージ
  ・rddensity・・・頑強なノンパラメトリックRD推定向けのRパッケージ
  

• RDD (regression discontinuity design) とはデータの発生メカニズムに注目した準実験法　　

• 無作為割り付け (RCT) を使った分析であれば実験法と呼ばれるが、RCT を使えない時に使う分析方法なので準実験法 (pseudo experiment) と呼ばれる
  

• 局所的な範囲で RCT が成立しているとみなすことができるため「最も実験研究に近いデザイン」といわれる

• RDD は次の 4 通りで呼ばれる

・回帰不連続デザイン
・不連続回帰デザイン
・非連続回帰デザイン
・RD デザイン

• 自然法則ではない強制的なルールによって生まれる境界線を利用した分析

境界線の例：

・試験で 60 点未満なら補講に参加する
・法定飲酒年齢: 21歳（米国の場合）
・高齢者の医療費の自己負担割合: 70歳
・小学校での１クラスあたりの人数の上限: 30人
・選挙区定数: 1 人

2. 何がわかるのか？

• ここでは、RDD を使う事でどのようなことが明らかになるのかいくつか例を挙げてみよう

2.1 補講が生徒の試験点数に与える処置効果

• 前期試験で 60 点未満とった生徒が補講授業を受けた時の、補講が後期試験結果に与える処置効果を推定できる

前期試験で 60 点未満だった人：	補講を受けた・・・treatment group（処置群）
前期試験で 60 点以上だった人：	補講を受けない・・・control group（統制群）

• 「補講に参加したかどうか」という要因以外に、「後期試験結果」に影響を与える違いはないという前提が必要
  
• 前期試験で 60点の境界線の前後で「後期試験結果」に違いがあれば・・・
  →　その「後期試験結果」は「補講」によるものと見なせるはず
  
• 下の図に描かれている縦の赤い矢印が推定された処置効果の大きさであり「補講の効果あり」という結果が得られる

架空データ

2.2 飲酒年齢が若者の死亡者数に与える処置効果

• 飲酒が許される年齢の時に、飲酒が若者の死亡者数に与える処置効果を推定できる
  
• アメリカでは法定飲酒年齢は州によって異なる（18, 19, 20, or 21歳）
  
• 例えば、法定飲酒年齢が 21 歳だとする

その日に 21 歳になった人：	飲酒できる・・・treatment group（処置群）
21歳になる1日前の人：	飲酒できない ・・・control group（統制群）

• 「飲酒できるかどうか」という要因以外に、「死」に影響を与える違いはない
  →　21歳の誕生日の境界線の前後で「死亡率」に違いがあれば
  →　その「死」は「飲酒」によるものと見なせるはず
  
• 下の図に描かれている縦の赤い矢印が推定された処置効果の大きさであり「21歳という飲酒年齢が若者の死亡者数に正の処置効果を与える」という結果を得られる

Source: Angrist & Pischke (2015), Mastering Mertics, Figure 4.2.

2.3 自己負担割合の低下が健康に与える処置効果

• 70 歳以上の高齢者が通院する際の自己負担率が減ったとき (30%→10%) に高齢者の通院が増えるかどうか調べたいたいとき
  
• 70 歳以上の高齢者が通院する際の自己負担率が従来の 30%から10%に減った場合を考える

70 歳以上の高齢者：	自己負担率 10% ・・・treatment group（処置群）
70 歳未満の高齢者：	自己負担率 30% ・・・control group（統制群）

• 「自己負担率が 10%か30%か」という要因以外に、「外来患者数」に影響を与える違いはない
  →　70歳の境界線の前後で「外来患者数」に違いがあれば →　その「外来患者数」は「自己負担率」によるものと見なせるはず

• 下の図に描かれている縦の赤い矢印が推定された処置効果の大きさであり、70 歳以上の高齢者が通院する際の自己負担率が従来の 30%から10%に減ると「高齢者の通院が増える」という結果を得られる

Source: 中室・津川『「原因と結果」の経済学』(2017, p. 140)

3. RDD の基本的な考え方

• 境界線（カットオフ）を生み出すルールをランダム割り付けだと見なす
  

• 事前に決めた基準に従って処置群と統制群を割り当てる
  

• カットオフ付近に位置する人々や集団を処置群と統制群に分割する状況を作り出す

• 例えば、こんな感じに割り当てる

基準を少しだけ上回った人々	treatment group（処置群）
基準を少しだけ下回った人々	control group（統制群）

• カットオフ付近の処置群と統制群の個体はほとんど同質だと考えることができる
  
• 事前に決めた基準（カットオフ）に従って処置群と統制群を割り当てる
  → 処置群と統制群は偶然（ランダム）に決められる
  → 準実験法とみなせる
  → 自己選択バイアスを回避できる
  
• ここでの前提：処置群と統制群の特徴は似ている
  → これら 2 群を比較する
  → カットオフにおける局所的な平均処置効果 (LATE) を推定できる

ポイント ・処置群と統制群が比較できるのは、カットオフ付近だけ
・RDDでは「平均処置効果 (ATE: Average Treatment Effect)」を推定できない
・RDDで推定できるのは「局所的な平均処置効果 (LATE: Local Average Treatment Effect at the cutoff) 」

3.1 RDD で使う変数と処置効果の関係

結果変数 \(Y_i\)

• 結果を表す変数 (Outcome variable)

処置変数 \(D_i\)

・0 と 1 の値をとる
・「介入変数」とも呼ばれる

割当変数 \(R_i\)

・「強制変数」とも呼ばれる (Running variable; forcing variable)

• R は D と Y の両方に影響を与える（＝交絡変数）

• 推定したい効果：D（補講授業） が Y（後期試験の点数） に与える影響
  

• 割当変数 \(R\) の値のみで処置の割り付け \(D\) が決まる
  →　\(R\) が唯一の交絡変数
  →　\(R\) さえ回帰分析でコントロールすれば
  →　脱落変数バイアスを回避して、\(D\) の処置効果を推定できる
  

• 前述した「補講授業」の例にあてはめて図で表してみる

• \(c\) はカットオフ（＝ 60点）
  
• 境界線を事前に決める割当変数 \(R\) の値 \(R_c\) = カットオフ
  
• カットオフは cutoff、cut point、threshold、閾値（いきち）とも呼ばれる
  
• 大学の授業で前期と後期に 2 つの試験を受けるとする
  
• どちらの試験も100点満点
  
• 前期試験で60点未満の生徒は夏休みの間に補講を受けた（60点は含まない）
  
• 前期試験で60点以上の生徒は補講を受けない（60点を含む）
  
• 200人の架空生徒データ (exam.score.csv) を使って具体的に考えてみる
  
• exam.score.csv をダウンロードし、RProject フォルダー内の data フォルダに保存する
  
• データを読み取る

df0 <- read_csv("data/exam.score.csv")

• ここで使うのは次の 6 つの変数

\(R\)	前期試験点数（割当変数）	観察可能
\(D\)	「D = 1: 補講あり」, 「D = 0: 補講なし」	観察可能
\(Y\)	後期試験点数	観察可能
\(y0t\)	潜在的結果変数: \(Y_i(0)\)　生徒が補講を受けない場合の仮想現実	観察不可能（理論上の変数）
\(y1t\)	潜在的結果変数: \(Y_i(1)\)　生徒が補講を受けた場合の仮想現実	観察不可能（理論上の変数）
\(y1t-y0t = ITE\)	個体処置効果: ←これが分かれば ATEが推定可能！	観察不可能（理論上の変数）

• ここでは 2 種類の変数が表示されている
  
• ひとつは私たち「人間が観察できる変数」: \(R, D, Y\)
  
• もうひとつは私たち「人間が観察できない変数」：\(y0t, y1t, ITE\)
  
• 「私たちは過去に戻れない」→　 平均処置効果 ATE を推定することはできない
  
• 通常、私たちが実際に観察し、分析で使うことができるのは観察可能な 3 つの変数だけ: \(R\), \(D\), \(Y\)
  
• ここではドラえもんのおかげで私たちが「過去に戻れて」通常、人間が観察できない 2 つの潜在的結果変数（\(y0t\) と \(y1t\)）を同時に観察できると想定してみる
  →　2つの潜在的結果変数を同時に観察できるなら、ITE がわかるので、ATE を推定できる
  
• ここで使うデータを確認する

DT::datatable(df0)

3.2 個体処置効果 (ITE)

• まず、ITE（個体処置効果）= 個人の補講処置効果 \(y1t - y0t\) を計算する
  
• 「固体処置効果 (ITE)」を平均すれば、補講の「平均処置効果 (ATE) 」は簡単に計算できる
  
• 私たちが知りたいのは補講の処置効果・・・個体処置効果: \(y1t - y0t\)

• 「理論上の 2 つの変数 (\(y0t\) と \(y1t\))」が観察できるので、補講の個別処置効果 (ITE) は簡単に計算できる
  

• しかし、実際にはこれら 2 つの変数を同時には観察できない
  

• 上の図のピンクで網がけされた値は実際には観測されない
  

• 網がけされていない値だけが観察可能
  

• 例えば、上から 2 人目の生徒を例に考えてみる
  

• \(R = 62, D = 0, Y = 41, y0t = 41, y1t = 79, ITE = 38\)
  

• この生徒の「前期試験点数 (\(R\)) 」は 62 点なので、現実には補講を受けていない (\(D = 0\))
  

• この生徒の補講の個別処置効果 (ITE) を計算するためには \(y0t\) と \(y1t\) の両方の情報が必要
  

• 補講を受けていないので、現実 \(y0t = 41\) は観測できるが、仮想現実 \(y1t = 79\) は観測できない
  

• この生徒は現実では、前期試験で 62 点とれたので、補講を受けず (\(D = 0\))、実際に観察されたのは「補講を受けた場合の仮想現実 (\(y0t)\)」の値なので、観測された後期試験点数 (\(Y\)) は 41 点
  

• ここで、生徒個人の補講の処置効果 ITE は \(y1t - y0t = 79 − 41 = 38\) と求めることができる

• これらの変数の関係を可視化すると次のようになる
  

• 右側の図は「前期試験結果にかかわらず、全員が試験を受けなかった」という仮想現実
  

• 左側の図は「前期試験結果にかかわらず、全員が試験を受けた」という仮想現実

• 通常、私たちが実際に観察できるのは「前期試験 R で60点以下で、補講を受けた人」（左図赤い点●）と「前期試験 R で60点以上で、補講を受けなかった人」（右図青い点●）だけ
  
• ここでは「灰色の点●を含めたすべてのドットを観測できる」と想定しているので、個別処置効果 (ITE) が計算できる

3.3 平均処置効果 (ATE)

• 個別処置効果 (ITE) と平均処置効果 (ATE) の違い

個別処置効果 (ITE)：	ある個体が処置を「受けた場合」の結果と「受けなかった場合」の結果の差
平均処置効果 (ATE)	潜在的結果 1 (\(y1t\)) から潜在的結果 0 (\(y0t\)) の値を引いて平均したもの: = 期待値

• ここでは通常は観測できない「理論上の変数 (\(y0t\), \(y1t\)) が観測可能」という前提がある
  →　つまり「灰色の点●を含めたすべての点観を測できる」と想定している
  →　個別処置効果 (ITE) が計算できる
  →　平均処置効果 (ATE) が計算できる 　　
• 平均処置効果 (ATE) は前期試験点数にかかわらず「全員が補講を受けた時の後期試験点数 \(y1t\)の平均値」から「全員が補講を受けなかった時の後期試験点数 \(y0t\)の平均値」を引いた値

• つまり、平均処置効果 (ATE) は左図の「赤い点●と灰色の点●」の平均値（縦軸）から右図の「青い点●と灰色の点●」の平均値の差 (8.175)

• 平均処置効果 ATE とは「処置群」（左図）と「統制群」（右図）の両方を含むすべての個体に対する効果なので、次のように表現できる 　　

\[平均処置効果 (ATE) = E[Y_i(1)-Y_i(0)] \\= E[Y_i(1)] - E[Y_i(0)]\]

私たちが設定したATE 値 ・補講を受けて後期試験を受けた場合の期待値: \(E[Y_i(1)]\)

mean(df0$y1t)

[1] 73.605

・補講を受けずに後期試験を受けない場合の期待値: \(E[Y_i(0)]\)

mean(df0$y0t)

[1] 65.43

・補講の平均処置効果 (\(ATE\)) は両者の差なので

mean(df0$y1t) - mean(df0$y0t)  

[1] 8.175

・私たちが設定した平均処置効果 (ATE) は 8.175点
・ ここで平均処置効果 (ATE) が計算できたのは、ドラえもんのお陰で「過去に戻れた」と想定しているためであることに注意
・通常、観察できるのは 3 つの変数 (\(R\), \(D\), \(Y\)) だけなので「補講の平均処置効果 ATE」8.175 は計算できない
・通常、私たちが観察できるのは「前期試験 R で60点以下で、補講を受けた人」（左図赤い点●）と「前期試験 R で60点以上で、補講を受けなかった人」（右図青い点●）だけ・・・灰色の点●は観察できない

3.4 局所的な平均処置効果 (LATE)

• 通常、私たちが観察できる変数は \(R\) と \(Y\) だけ（赤い点●と青い点●だけ）
  →　補講授業の平均処置効果 (ATE) は推定できない

• ドラえもんがいないので、潜在的結果変数 \(y0t, y1t\) は観察できない
  
• どうするか？　　

解決策

• 観察できる変数 (\(R\) と \(Y\)) だけを使って、2 つの局所線形回帰モデルを使い局所的な平均処置効果 (LATE: Local Average Treatment Effect at the cutoff) を推定する
  
• 閾値の左側と右側にそれぞれ別々の線形回帰モデルを当てはめ、閾値周辺のデータだけを使って、交差項を含む共分散分析を行う

・局所的な平均処置効果 (LATE) 」の定義

共変量 \(R_i\) が閾値 \(c\) の値をとるときの \(Y_i(1)\) と \(Y_i(0)\) の差の期待値

\[E[Y_i(1)-Y_i(0)|R_i = c] = E[Y_i(1)|R_i = c] - E[Y_i(0)|R_i = c]\]

・ここでは、前期試験点数が 60 点（閾値 \(c = 60\)）の時の \(Y_i(1)\) と \(Y_i(0)\) の差のこと
・前期試験で 59点とって補講を受けた生徒と、60点とって補講を受けなかった生徒は「ほとんど同じ」だと考える
・その「ほとんど同じ人」が「補講を受けた場合」と「補講を受けない場合」の後期試験結果を比較すれば補講の処置効果がわかる、というのが基本的な考え方

平均処置効果 (ATE) と局所的な平均処置効果 (LATE) の違い

平均処置効果 (ATE)	潜在的結果 1 (\(y1t\)) から潜在的結果 0 (\(y0t\)) の値を引いて平均したもの: = 期待値
局所的な平均処置効果 (LATE)：	共変量 \(R_i\) が閾値 \(c\) の値をとるときの \(Y_i(1)\) と \(Y_i(0)\) の差の期待値

私たちが設定した LATE の真値 ・理論的には、閾値における局所的な平均処置効果 (LATE) は、共変量 \(X_i\) が閾値 \(c\) の値をとるときの \(Y_i(1)\) と \(Y_i(1)\) の差の期待値

\[E[Y_i(1)-Y_i(0)|R_i = c] = E[Y_i(1)|R_i = c] - E[Y_i(0)|R_i = c]\]

・ここでは潜在的結果変数 \(y0t, y1t\) が観察可能と想定しているので、補講を「受けた場合」と「受けなかった場合」それぞれの回帰式は次のように設定できた

\[Y_i(0) = 35.9 + 0.486R \\ Y_i(1) = 35.4 + 0.629R\]　

• 具体的には、ここでは共変量 \(X_i\)（前期試験点数 \(R_i\)）が閾値 \(c = 60\) の値をとる時の後期試験点数 \(Y_i(1)\) と \(Y_i(0)\) の差の期待値のこと
  
• 従って、次のように計算できる

\[E[Y_i(1)-Y_i(0)|R_i = c] = E[Y_i(1)|R_i = c] - E[Y_i(0)|R_i = c]　\\ = E[35.4 + 0.629R|R_i=60]-E[35.9+0.486R|R_i=60] \\ = 35.4 + 0.629*60 - 35.9 + 0.486 * 60 \\ = 73.14 - 65.06 \\= 8.08\]

・私たちが設定した局所的な平均処置効果 (LATE) は 8.08点
・ ここで局所的な平均処置効果 (LATE) が計算できたのは、ドラえもんのお陰で「過去に戻れた」と想定しているためであることに注意
・通常、観察できるのは 3 つの変数 (\(R\), \(D\), \(Y\)) だけなので「補講の局所的な平均処置効果 (LATE) 」8.08 は計算できない
・通常、私たちが観察できるのは次の 2 つだけ：
・「前期試験 R で60点以下で、補講を受けた人」（赤い点●）
・「前期試験 R で60点以上で、補講を受けなかった人」（青い点●）
・ 灰色の点●は観察できない

・これ以降、私たちが設定した局所的な平均処置効果 (LATE) 8.08点をどれだけ上手に推定できるかを試みる

3.5 やってはいけない比較

• ここでは、補講の処置効果を推定する上でやってはいけない（しかしやりがちな）「ナイーブで単純な比較」と「やってもいいのではないかと思われるが問題を含む比較」を 6 つ紹介する

3.5.1 ナイーブな前後比較

• 「後期試験の平均点」から「前期試験の平均点」を引く単純な前後比較

mean(df0$Y) # 後期試験の平均点

[1] 67.755

mean(df0$R) # 前期試験の平均点

[1] 60.69

mean(df0$Y) - mean(df0$R) # 「後期試験の平均点」ー「前期試験の平均点」

[1] 7.065

• 後期試験の平均点が前期試験の平均点より 7 点高い
  →　「この 7 点は補講の処置効果だ」と推論してしまう誤り
  
• 後期試験の点数の平均が高いのは、ただ単に後期試験の方がより簡単だったためかもしれない
  →　「後期試験の平均点」から「前期試験の平均点」を引く単純な前後比較では補講授業の処置効果は推定できない

3.5.2 ナイーブな平均値の比較

・「補講を受けた人」と「補講を受けなかった人」の「後期試験点数：Y」の平均値の比較

・「補講を受けた人 D = 1」の後期試験点数の平均値（赤い●だけの平均値）

mean(df0$Y[df0$D=="1"])  

[1] 67.15625

• 「補講を受けなかった人 D = 0」の後期試験点数の平均値（青い●だけの平均値）

mean(df0$Y[df0$D=="0"])  

[1] 68.30769

• 「補講を受けた人 D = 1」の後期試験点数の平均値から「補講を受けなかった人 D = 0」の後期試験点数の平均値を引く

mean(df0$Y[df0$D=="1"]) - mean(df0$Y[df0$D=="0"])

[1] -1.151442

• 「補講を受けたと後期試験点数が約 1 点下がる！」という結果
  
• 交絡があると、このような異常な結果を引き起こす
  → 「補講を受けた人」と「補講を受けなかった人」の「後期試験点数：Y」の平均値を単純に比較しても、補講授業の処置効果は推定できない
  
• このような単純比較法は欠測データ解析では「ペアワイズ除去法」(pairwise deletion) と呼ばれる
  
• ペアワイズ除去法は潜在的結果変数の欠測が完全に無作為でない限り、妥当な方法とはいえない

3.5.3 閾値周辺の観測値平均の比較 (1)

• 補講授業の処置効果を推定するための解決策は「観察できる変数 (\(𝑅\) と \(𝑌\)) だけを使って、局所的な平均処置効果 (LATE) を推定すること」
  
• では、やってもよさそうな比較とはどのような比較なのか？
  
• 前期試験で 59点とって補講を受けた生徒と、60点とって補講を受けた生徒の後期試験の違いを確かめてみる
  
• これらの生徒はたまたま59点とり、たまたま60点とった「ほとんど同質」の生徒と考えることができる

• 前期試験で 59点とった生徒が 10人おり、彼らの後期試験の平均は73.8点

fiftynine <- df0 |> 
  dplyr::filter(R == 59) |> 
  dplyr::select(Y) 

dim(fiftynine)

[1] 10  1

mean(fiftynine$Y)

[1] 73.8

• 前期試験で 60点とった生徒が3人おり、彼らの後期試験の平均は72.3点

sixty <- df0  |> 
  dplyr::filter(R == 60) |> 
dplyr::select(Y) 

dim(sixty)

[1] 3 1

mean(sixty$Y)

[1] 72.33333

• 両者の差 (73.8 - 72.3 = 1.47)が補講の処置効果と考えることもできる

mean(fiftynine$Y) - mean(sixty$Y)

[1] 1.466667

• この場合「補講の処置効果は1.47」ということになる

• しかし、この結果を一般化するには生徒の数が少なすぎる

3.5.4 閾値周辺の観測値平均の比較 (2)

• カットポイント 60 点の前後 5 点の生徒の後期試験点数の平均値を比較してみる

• 前期試験で 55〜59点とった生徒が44名、彼らの後期試験の平均は約73点

minus5 <- df0 |> 
  dplyr::filter(R > 54 & R < 60) |> 
  dplyr::select(Y) 

dim(minus5)

[1] 44  1

table(minus5$Y)

55 56 59 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 78 79 80 81 82 83 84 
 1  1  1  1  1  1  2  1  1  3  3  1  1  1  1  3  2  2  2  1  1  1  1  3  3  1 
85 86 90 
 1  1  2 

mean(minus5$Y)

[1] 73.15909

• 前期試験で 60〜64点とった生徒は29名、彼らの後期試験の平均は約63点

plus5 <- df0  |> 
  dplyr::filter(R > 59 & R < 65) |> 
dplyr::select(Y) 

dim(plus5)

[1] 29  1

table(plus5$Y)

41 50 51 55 56 57 59 61 62 63 64 67 68 69 70 76 77 78 83 
 1  1  1  1  2  3  3  1  3  2  2  1  1  2  1  1  1  1  1 

mean(plus5$Y)

[1] 62.58621

• 両者の差 (73.15909 - 62.58621)が補講の処置効果と考えることもできる

mean(minus5$Y) - mean(plus5$Y)

[1] 10.57288

• この場合「補講の処置効果は10.6」ということになる

どこまで前期試験点数の範囲を広げるか？ ・ここでは前期試験点数が 59 点と 60 点の場合、カットポイント 60 点の前後 5 点の場合の二種類を検討
・しかし、もし範囲をすべての範囲に広げると・・・
→「3.5.2 ナイーブな平均値の比較」にたどり着いてしまう

3.5.5 交絡を含めた重回帰分析

• 補講授業の事例では、処置変数 \(D_i\) は共変量 \(R_i\) によって確定的に割り付けられているので、\(R_i\) は交絡である
  
• 交絡をモデルに含めた共分散分析（＝重回帰分析）を実行すれば処置変数 \(D_i\) の係数を補講授業の処置効果と見なせるはず

model1 <- lm(Y ~ R + D, data = df0)

tidy(model1)

term	estimate	std.error	statistic	p.value
(Intercept)	17.2	5.91	2.92	0.00392
R	0.733	0.0837	8.75	9.53e-16
D	12.6	2.04	6.19	3.5e-09

\[Y = 17.2 + 0.733R + 12.6D\]　

• ここで得られた処置 \(D_i\) の係数は 12.6 →　補講授業の効果
  
• 補講授業の処置効果の真値 (8.175) とはかなりかけ離れている
  
• その理由・・・補講を「受けた場合」と「受けなかった場合」とで回帰式の関数系が異なるため
  
• 交絡をモデルに含めた共分散分析（＝重回帰分析）を実行しても処置変数 \(D_i\) の係数を補講授業の処置効果と見なすことは難しい

補講を「受けた場合」と「受けなかった場合」の回帰式 ・ここでは潜在的結果変数 \(y0t, y1t\) が観察可能と想定しているので、補講を「受けた場合」と「受けなかった場合」それぞれの回帰式は次のようになる

\[Y_i(0) = 35.9 + 0.486R \\ Y_i(1) = 35.4 + 0.629R\]　 ・\(Y_i(0)\) の式を求めるための回帰分析

m0 <- lm(y0t ~ R, data = df0) 
tidy(m0)

term	estimate	std.error	statistic	p.value
(Intercept)	35.9	3.42	10.5	8.14e-21
R	0.486	0.0552	8.81	6.49e-16

\[Y_i(0) = 35.9 + 0.486R\]

・\(Y_i(1)\) の式を求めるための回帰分析

m1 <- lm(y1t ~ R, data = df0)
tidy(m1)

term	estimate	std.error	statistic	p.value
(Intercept)	35.4	3.55	9.97	3.05e-19
R	0.629	0.0574	11	3.92e-22

\[Y_i(1) = 35.4 + 0.629R\]

3.5.6 変数変換・交差項を含めた重回帰分析

• 適切に変数変換したり、交差項 (\(D_i:R_i\)) を重回帰モデルに加えることで「補講」の処置効果が推定できそう
  
• しかし、この場合「外挿」が問題になる
  
• 「外挿」とはデータの存在しない部分を予測すること
  
• 下の図でいえば、赤い点●と青い点●を使って灰色の点●を予測すること

• データの存在しない部分を予測する「外挿」はできるだけ少なくする必要がある
  
• 灰色の点●を含めた潜在的結果変数 (\(y0t, y1t\))を同時に観測できれば、処置群と統制群の違いを適切にモデリングできる　　
• しかし、私たちは上図のような潜在的結果変数 (\(y0t, y1t\)) を同時に観察できない
  
• 私たちが観測できるのは、潜在的結果変数 (\(y0t, y1t\)) のいずれか片方のデータを組み合わせた下図のような散布図

• 実際に観測できる赤い点●と青い点●から 灰色の点●を含めたモデリング（＝外挿）はできない
  →　変数変換したり、交差項 (\(D_i:R_i\)) を重回帰モデルに加えることで「補講」の処置効果を適切に測定することは難しい

3.6「補講の処置効果」結果のまとめ

• 「補講の処置効果」に関して、ここまで実施した分析結果をまとめてみる

平均処置効果 (ATE)の真値	8.175
局所的な平均処置効果 (LATE)の真値	8.08

分析方法	具体的な分析内容	推定された処置効果
3.5.1 ナイーブな前後比較	「後期試験の平均点」から「前期試験の平均点」を引く単純な前後比較	7.065
3.5.2 ナイーブな平均値の比較	補講を「受けた人」と「受けなかった人」の後期試験点数の平均値の比較	-1.15
3.5.3 閾値周辺の観測値平均の比較 (1)	閾値前後 2 つの観測値の比較	1.47
3.5.4 閾値周辺の観測値平均の比較 (2)	閾値±5 の複数の観測値の比較	10.57
3.5.5 交絡を含めた重回帰分析	交絡を含めた重回帰分析	12.6
3.5.6 変数変換・交差項を含めた重回帰分析	閾値前後の観測値の比較

• 我々が設定した処置効果の真値 (8.175) と比較すると、他の 5 つの推定値がかけ離れており、正しく推定できていないことがわかる
  
• 3.5.1 ナイーブな前後比較は論外として、平均処置効果の真値 (8.175) に最も近い値は 3.5.4 閾値周辺の観測値平均の比較 (2)の 10.57
  →　重回帰分析と閾値を使った局所的な平均処置効果 (LATE) を推定する必要がある

4. RDD による局所的な平均処置効果 (LATE)

4.1 RDD の長所

(1) 回帰不連続デザインは最も実験研究に近いデザイン

• カットオフ付近に限定すれば、ランダム化比較試験 (RCT) と見なすことができる
  

• 回帰不連続デザインは、共分散分析を統計的因果推論の手法として拡張したもの
  

• 閾値（いきち・しきいち）付近の値は無作為な副標本 (subsample) と考えられる
  

• 閾値付近の値（ここでは前期試験点数 = 60点）の前後の点数をとった生徒の学力は、誤差範囲内でほぼ同じと考えることができる
  

• 閾値付近のデータだけに限定すれば、その局所的なデータにおいては処置が無作為に割り付けられていると考えることが可能
  →　事実上の実験研究とみなせる
  

• 回帰不連続デザインでは割当変数 \(R_i\) が唯一の交絡変数
  

• 自然実験 (natural experiment)や準実験 (quasi-experiment) と位置づけることが可能

• 但し、RDDの 3 つの仮定が満たされる場合に限る

(2) 結果を可視化して示すことができる = 結果がわかりやすい

• カットオフにおける「ジャンプ」を示す
  => 局所的な平均処置効果 (LATE)
  
• カットオフを上手に工夫すれば、様々な社会現象に応用可能

4.2 RDD の短所

(1) RDD の仮定が正しいことを証明できない

• これはどのような分析にも言えること
  
• 仮定が正しいことを証明できないので、傍証を集めて信憑性を高める

(2) 処置効果を推定できるのはカットオフ付近のみである

• 局所的な平均処置効果 (LATE) は閾値付近の個体だけに適応可能な処置効果
  →　適応範囲がかなり限定的
  →　適応範囲が限定的ではあっても、分析目的によってはかなり有益
  
• カットオフ付近のデータは全体の一部に過ぎない
  →　多くのサンプルを捨てることになる
  →　サンプル全体についての処置効果はわからない
  
• サンプルサイズが多ければ、情報量が多くなるので、不確実性の小さい推論ができる

(3) バンド幅 (bandwidth) によって推定結果が変化する

• 処置効果を推定する際には「どれだけのサンプルを捨てるか」が問題になる
  
• バンド幅を狭める　→　正確な因果推定ができるが統計的不確実性が増す
  
• バンド幅を広げる　→　全体の推定は可能だが処置効果の正確性が減少
  
• 正しい関数形がわからないとバイアスが生じる
  →　正確な因果推論ができない

4.3 RDD の仮定

• RDD で処置効果を推定するためには 3 つの仮定を満たす必要がある

4.3.1 仮定 (1)

カットオフによって処置の値が変わらない場合、結果変数の値が非連続的なジャンプをしない

• つまり

\[E[Y_i(0)|R_i] と E[Y_i(1)|R_i]が R_c で連続な関数である\]

・\(E[Y_i(0)]\): 処置が 0 の時（つまり補講を受けない時）の潜在的結果

・\(E[Y_i(1)]\): 処置が 1 の時（つまり補講を受けた時）の潜在的結果

・\(R_c\): 割当変数 \(R\) の値 = カットオフ

• 「カットオフによって処置の値が変わらない」 = 「カットオフを作るルールはない」
  
• 例えば、60点のところに「ジャンプ」があるのは、偶然なのではないかと考えてみる
  
• 「補講を受けたかどうかが決まる 60点でなければ、60点でのジャンプはない」と仮定したい
  →　そのためには例えば「50点と 70点ではジャンプがない」ことを示せばよい
  →　そうすれば説得力が増す

割当変数 \(R_i\)の非連続性を確認理由 処置群と統制群を分ける際の恣意性を排除するため
・Figure 2 はマラソンに参加した 950万人が完走するために要した時間の分布図
・時間は 30 分ごとに示している

Source: https://www.nber.org/system/files/working_papers/w20343/w20343.pdf

・3時間30分、4時間、4時間30分など区切りのいいタイム直前にゴールしている人々が多い
・区切りのいいタイム直後にゴールしている人々が少ない

このグラフから分かること

・「もう少し頑張れば 4 時間以内にフィニッシュできる！」と考えて人々は普段の実力以上に頑張る
→　本来走るのが「遅い人」が「速い人」に分類されてしまう
→　処置群と統制群を分ける際に恣意性が含まれてしまう (manipulation)
→　マラソン走者にとってはいいことだが、分析者にとってはマイナス

4.3.2 仮定 (2)

カットオフで「非連続的な変化」をするのは処置変数だけ

• カットオフ付近で、結果に影響を与える他の「処置」がないことが必要
  

• 例えば、この仮定が破られる例として・・・

• 前期試験で60点未満だった生徒たちにだけ、後期試験の予想問題が配付されたとする
  →　この場合、後期試験得点が高くなったのは「補講」のせいなのか「予想問題が配布された」せいなのか区別できない

• 後期の予想問題配布のような処置がなければ・・・
  →　カットオフで「非連続的な変化」をするのは処置変数だけという仮定は満たされる

4.3.3 仮定 (3)

割当変数の値を自分では選べない

• 生徒たちは補講に参加するかどうかということを自分で選ぶことはできない
  

• カットオフ付近で、割当変数の分布に不自然な変化がないかどうか確認する必要がある
  

• 例えば、この仮定が破られる例として・・・

• 前期試験結果において、60点の生徒の数が異常に多い場合
  →　何らかの不正が行われた可能性を否定できない
  ・カンニングペーパーの使用、試験監督への賄賂などなど
  

• カットオフ付近のヒストグラムを描いて確認する必要がある

II. 補講授業効果の推定

• ここでは架空データを使って補講授業の処置効果を推定してみる

5. 補講授業の処置効果を推定する

RDDの推定手順

Step 1: 処置群と統制群を分けるルールはあるか？	ルールの確認
Step 2: RDDデザインは sharp か fuzzy か？	散布図で確認
Step 3: カットオフ付近での結果変数 \(Y\) の非連続性	散布図とRDプロットで確認
Step 4: カットオフ付近での割当変数 \(R\) の非連続性	理論的に検討
Step 5: 割当変数 \(R\) の値を自分では選べない	ヒストグラム、rdplotdensity()関数、 rdensity()関数で検定
Step 6: 処置効果の推定（パラメトリック・ノンパラメトリック）
Step 7: 推定結果を比較する

• exam.score.csv をダウンロードし、RProject フォルダー内の data フォルダに保存する
  
• データを読み取る

df0 <- read_csv("data/exam.score.csv")

df1 <- df0  |>  
  select(R, D, Y)

• データフレーム df0 の記述統計を確認する

stargazer(as.data.frame(df1), type ="html")

Statistic	N	Mean	St. Dev.	Min	Max
R	200	60.690	12.202	27	92
D	200	0.480	0.501	0	1
Y	200	67.755	10.767	41	100

DT::datatable(df1)

5.1 Step 1: 処置群と統制群を分けるルールはあるか？

• 処置群と統制群が何らかのルールで決められているかどうかを確認する
  →　決められている:
  
• 割当変数 \(R\) の値が 60 点以下なら「補講を受けた」、60点以上なら「補講を受けず」

5.2 Step 2: RDD の種類は sharp か fuzzy か？

• RDD には sharp と fuzzy の 2 種類ある

Sharp RDD	割当変数 \(R\) の値が分かれば、処置の値が確定する場合
Fuzzy RDD	割当変数 \(R\) の値によって処置される確率が決まり、その確率がカットオフ上でジャンプする場合

• ここで使う RDD が Sharp なのか Fuzzy なのかを判断するため、\(x\) 軸に「割当変数」（前期試験点数）、\(y\) 軸に「処置変数」（補講を受けたかどうか）の散布図を描いて確かめる

theme_set(theme_bw(base_family = "HiraKakuProN-W3"))

ggplot(data = df1,
       mapping = aes(x = R, 
                     y = D, 
                     color = factor(D))
       ) +
  geom_point(alpha = 0.5, 
             alpha = 0.5,
             position = position_jitter(width = 0, height = 0.15)) +
  geom_vline(xintercept = 60, 
             color = "purple", 
             linetype = 2) +
  scale_color_manual(values = c("blue", "red"))

• カットオフ(60点) を境界線として「補講を受けた生徒数」と「補講を受けなかった生徒数」が明確に分かれている
  
• 60点以上取得した生徒は補講を受けていない
  
• 数値で確認する

df1 |> 
  group_by(D, R < 60) |> 
  summarize(count = n())

# A tibble: 2 × 3
# Groups:   D [2]
      D `R < 60` count
  <dbl> <lgl>    <int>
1     0 FALSE      104
2     1 TRUE        96

• 補講を受けた生徒数が 96 名
  
• 補講を受けなかった生徒数が 104 名
  
• オーバーラップする生徒がいない → RDD デザインは sharp

5.3 Step 3: カットオフ付近での結果変数 (Y) の非連続性

• 仮定 (1): カットオフによって処置の値が変わらない場合、結果変数の値が非連続的なジャンプをしないの確認

• 散布図と RDプロットで確認

散布図 (c=60)

• カットオフ付近での結果変数 \(Y\) の非連続性を確認
  
• ここでは非連続性（= ジャンプ）があることを期待

ggplot(df1, aes(x = R, y = Y, color = factor(D)))+
  geom_vline(xintercept = 60, color = "purple", linetype = 2) +
  geom_point(size = 1, alpha = 0.5) + 
  geom_smooth(data = filter(df1, R <= 60), method = "lm", se = FALSE) +
  geom_smooth(data = filter(df1, R > 60), method = "lm", se = FALSE) +
  scale_color_manual(values = c("blue", "red"))

• 割当変数 \(R_i\)の カットオフ付近ではジャンプが認められる
  →　補講の処置効果が期待できる
  
• 「カットオフによって処置の値が変わらない」 = 「カットオフを作るルールはない」
• 例えば、60点のところに「ジャンプ」があるのは、偶然なのではないかと考えてみる
  
• 60点でなければ、60点でのジャンプはない」と仮定したい
  →　そのためには例えば「50点と 70点ではジャンプがない」ことを示せばよい
  →　そうすれば説得力が増す

RD プロット (c=60)

• RDD では、上記のような通常の散布図に加えて RD プロットと呼ばれる特殊なグラフを使う
  

• RD プロットとは、散布図の横軸をいくつかのビン (bin) に分割して、そのビンの中に入る \(Y\) の平均値をドットで示した図

• rdrobustパッケージの rdplot()関数 を使って RD プロット を表示してみる

library(rdrobust)

rdplot(df0$Y, df0$R, c = 60)

[1] "Mass points detected in the running variable."

• データをスムーズ化しているので、 \(Y\) にジャンプがあるかどうか確認しやすい

• 「カットオフによって処置の値が変わらない場合」の一例として50点でのジャンプの有無を確認してみる

散布図 (c=50)

df1 |>
  mutate(threshold = as.factor(ifelse(R < 50, 1, 0))) |>
  ggplot(aes(x = R, y = Y, color = threshold)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  guides(color = FALSE) +
  geom_vline(xintercept = 50, color = "black", linetype = 2)+
  scale_color_manual(values = c("blue", "red"))

RD プロット (c=50)

rdplot(df0$Y, df0$R, c = 50)

[1] "Mass points detected in the running variable."

• 50点でのジャンプは認められない
  
• むしろマイナスにジャンプしている
  
• もう一つ「カットオフによって処置の値が変わらない場合」の例として70点でのジャンプの有無を確認してみる

散布図 (c=70)

df1 |>
  mutate(threshold = as.factor(ifelse(R < 70, 1, 0))) |>
  ggplot(aes(x = R, y = Y, color = threshold)) +
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  guides(color = FALSE) +
  geom_vline(xintercept = 70, color = "black", linetype = 2) +
  scale_color_manual(values = c("blue", "red"))

RD プロット (c=70)

rdplot(df0$Y, df0$R, c = 70)

[1] "Mass points detected in the running variable."

• やはり、70点でも大きなジャンプは認められない
  →　カットオフによって処置の値が変わらない場合（50点と70点の場合）
  →　結果変数の値が非連続的なジャンプをせず、60点の時にだけ比較的大きな「ジャンプ」が確認できた
  → 仮定 (1): カットオフによって処置の値が変わらない場合、結果変数の値が非連続的なジャンプをしないは満たされる

5.4 Step 4: カットオフ付近での割当変数 \(R\) の非連続性

• 仮定(2): カットオフで『非連続的な変化』をするのは処置変数だけの確認
  

• カットオフ付近で、結果に影響を与える他の「処置」がないことが必要
  

• 例えば、この仮定が破られる例として・・・

• 前期試験で60点未満だった生徒たちにだけ、後期試験の予想問題が配付されたとする
  →　この場合、後期試験得点が高くなったのは「補講」のせいなのか「予想問題が配布された」せいなのか区別できない

• 理論的に検討してみて、後期の予想問題配布のような処置がないなら・・・

→　仮定 (2): カットオフで「非連続的な変化」をするのは処置変数だけは満たされる

5.5 Step 5: 割当変数 \(R\) の値を自分では選べない

• 仮定 (3): 割当変数の値を自分では選べないの確認

• ヒストグラム、rdplotdensity()関数、rdensity()関数で検定

ヒストグラム

• 生徒たちは補講に参加するかどうかということを自分で選ぶことはできない
  
• カットオフ付近で、割当変数の分布に不自然な変化がないかどうか確認する必要がある
  
• 割当変数 \(R_i\)のカットオフ付近での非連続性を確認
  
• ここでは非連続性（= ジャンプ）がないことを期待
  
• 前期試験（割当変数：R）で不正が行われたとする
  
• この場合、カットオフである 60 点を少し超える点数を獲得する生徒が多くなるはず（= 非連続になるはず）
  
• 不正がなければ、カットオフである 60 点付近の生徒数はほぼ同じ（= 連続）になるはず
  
• もし、カットオフ前後で生徒数に大きな差がある（= 非連続）ようならそれは問題
  
• カットオフ付近のヒストグラムを描く

ggplot(df1, aes(x = R, fill = D)) +
  geom_histogram(binwidth = 5, boundary = 60, color = "white") +
  geom_vline(xintercept = 60, color = "purple", linetype = 2) 

• ヒストグラムを目視する限り、60点を超えたあたりの人数が多いようにはみえない

rdplotdensity()関数で検定

• 目視だけだと統計的有意性が分からないので、rdplotdensity()関数 を使って95%信頼区間を表示してみる

rdplotdensity(rdd = rddensity(df1$R, c = 60),
              X = df1$R) |> 
  summary()

        Length Class     Mode
Estl     4     lpdensity list
Estr     4     lpdensity list
Estplot 10     gg        list

• カットオフ付近での 95%信頼区間がオーバーラップしている
  → カットオフである 60 点付近の生徒数に差がない
  → 連続している
  
• カットオフである 60 点を少し超える点数を獲得する生徒は増えていない
  
• むしろ数が少ない
  →　点数の不正操作はないと思われる

rddensity()関数で検定

• 割当変数の操作の確認は rddensity パッケージの rddensity()関数を使っても確認できる

library(rddensity)

summary(rddensity(df0$R, c = 60))

Manipulation testing using local polynomial density estimation.

Number of obs =       200
Model =               unrestricted
Kernel =              triangular
BW method =           estimated
VCE method =          jackknife

c = 60                Left of c           Right of c          
Number of obs         96                  104                 
Eff. Number of obs    92                  97                  
Order est. (p)        2                   2                   
Order bias (q)        3                   3                   
BW est. (h)           26                  26                  

Method                T                   P > |T|             
Robust                -0.8813             0.3782              


P-values of binomial tests (H0: p=0.5).

Window Length / 2          <c     >=c    P>|T|
1.000                      10      11    1.0000
2.000                      18      17    1.0000
3.000                      24      24    1.0000
4.000                      34      29    0.6147
5.000                      44      39    0.6609
6.000                      49      48    1.0000
7.000                      54      54    1.0000
8.000                      56      58    0.9254
9.000                      60      64    0.7877
10.000                     65      65    1.0000

• 注目すべきところ

• ここでの帰無仮説：「カットオフ前後の平均値は連続している」
• p-value = 0.3782 →　帰無仮説は棄却されず
  → 連続している　→　割り当て変数 \(R\) に操作は行われていない

・仮定 (3): 割当変数の値を自分では選べないは満たされる

5.6 Step 6: 処置効果を推定する

• 処置効果の推定方法は次の 2 つ

パラメトリック RD	: 特定の関数形による回帰式で処置効果を推定する方法
ノンパラメトリック RD	: 閾値の周辺データだけを使って、交差項を含む共分散分析で処置効果を推定する方法

• それぞれの推定方法の特徴：

推定の方法	関数形の仮定	仮定
パラメトリック RD	仮定する	厳しい
ノンパラメトリック RD	仮定しない	緩い

一般的に、パラメトリック RD は関数形の特定が難しいので、ノンパラメトリック RD の結果を採用することが多い

パラメトリック RD

• 分析者の事前知識を仮定としてモデルに取り入れる方法
  
• 「もしこういうふうに仮定したら、どういう結果になるだろう？」と考える
  
• 関数系にも外生変数の分布にも仮定をおくアプローチ
  
• 通常は関数系には線形、外生変数の分布にはガウス分布を仮定
• 仮定をおくことでモデルが単純になったり、数理的な扱いが楽になるというのが利点
  
• 回帰式の関数形を正しく特定することが難しいというのが難点

ノンパラメトリック RD

• 関数系にも外生変数の分布にも仮定をおかないアプローチ
  
• 関数系は線形でも非線形でも何でもよい
  
• 外生変数の分布も何でもよい
  
• できるだけ仮定をおかずに、データだけでどこまで推測できるかという限界を明らかにすることを重視

ここからはそれぞれの推定方法を解説し実際に推定する

6. パラメトリック RD（補講授業）

6.1 パラメトリック RD とは

• カットオフ付近以外も含む推定
  
• ノンパラメトリック RD との違い:
  
• 仮定が多い　→　考えるべき事が多い
  
• 「\(Y\) を \(X\) に回帰する」とき、パラメトリック回帰では次のような回帰関数を想定する

\[Y_i = f(X_i)\]

• パラメトリックはさらにその先の「関数形」を考える必要あり
  
• 特定の関数形による回帰式で処置効果を推定する 　
• 回帰式の関数形を正しく特定することが必要
  →　そうでないと、処置効果を推定できない
  
• しかし、その関数形が正しいかどうかはわからない
  
• ここでは最も単純な回帰式を考えてみる

\[Y_i = α + ρD_i + βR_i + ε_i\]

• これは \(Y_i\) を \(D_i\) と \(R_i\) に回帰した単純な線形回帰

• 推定したいのは処置効果の推定値 \(ρ\)（ロー）
  

• 割当変数 \(R_i\) は交絡変数 →　コントロールする必要がある
  

• この関数形が正しければ
  →　\(ρ\)（ロー）の推定値 ＝ 処置効果の推定値

6.2 パラメトリック RD の潜在的結果

• パラメトリック RD の難しいところ：
  →　関数形だけに依存して、実際に観測していないところを予測しなければならないこと
  
• ここで想定する回帰関数:

\[E[Y_i | D_i, R_i] = α + ρD_i + βR_i\]

• \(D_i\)が 1 と 0 それぞれの場合の処置効果は

● \(D_i = 1: E[Y_i(1) | D_i = 1, R_i] = α + ρ + βR_i\)
● \(D_i = 0: E[Y_i(0) | D_i = 0, R_i] = α + βR_i\)

• 交絡変数は \(R\) だけだと分かっている
  →　条件付き交換可能性がある
  →　全ての \(R_i\) について \(ρ\)（ロー）は次のように表せる

\[ρ = E[Y_i(1)|R_i] - E[Y_i(0)|R_i]\]

• \(ρ\)（ロー）の期待値をとると

\[E[ρ] = E[E[Y_i(1)|R_i] - E[Y_i(0)|R_i]]\]
- これが平均処置効果 = \(R_i\) に条件付けられた平均処置効果

• ここで直面する問題：
  \(R_i\) に条件付けると、\(E[Y_i(1)|R_i]\) または \(E[Y_i(0)|R_i]\) のいずれか一方しか観察できない
  →　観測値のないところを推定する必要がある
  →　もし間違った関数形を使うと、推定を間違えてしまう

推定結果は関数形に依存する

• その理由は、実際には観測されていないところを関数で推定しようとするため

• パラメトリック RD でバイアスをなくすには、正しい関数形が必要
  

• しかし、関数形は無数に存在する

• 関数形の例：

\(R_i\) の次元が増える
\(Y_i = α + ρD_i + βR_i + ε_i\)
\(Y_i = α + ρD_i + β_1R_i + β_2R_i^2 + ε_i\)
\(Y_i = α + ρD_i + β_1R_i + β_2R_i^2 + β_3R_i^3 + ε_i\)
\(Y_i = α + ρD_i + β_1R_i + β_2R_i^2 + β_3R_i^3+・・・+β_kR_i^k+ ε_i\)

交差項が含まれる

\(Y_i = α + ρD_i + βR_i + γ(D_i×R_i) + ε_i\)
\(Y_i = α + ρD_i + βR_i + γ_1(D_i×R_i) + γ_2(D_i×R_i^2) + ε_i\)

• このようなモデルを仮定し、特定した関数形に現れるパラメタ (parameters) を推定するのがパラメトリック回帰

パラメトリック RD のジレンマ

• どのモデルが正しいのか事前に知ることはできない
• しかし、正しいモデルを選ばないと、正しく推定できない

関数形の違いによる推定の差に注意

曲線的な関係がありそうな時には曲線を当てはめる

• 2 次の項（曲線を想定）— \(R_i^2\)
• 交差項（傾きが異なると想定）— \(γ_1(D_i・R_i), γ_2(D_i・R_i^2)\)

\[Y_i = α + ρD_i + β_1R_i + β_2R_i^2 + γ_1(D_i・R_i) + γ_2(D_i・R_i^2) + u_i\]

→　処置効果を正しく推定

曲線的な関係がありそうなのに、直線を当てはめた場合

• 1 次の項（曲線を想定）— \(R_i\)
• 交差項を含めない（傾きが同じと想定）

\[Y_i = α + ρD_i + β_1R_i + u_i\]

→　処置効果を過小に推定してしまう

• 曲線の「フィット」を改善することが目的ではない
  
• 回帰式の多項式を増やせば、回帰曲線とデータの「フィット」は簡単に改善できる
  
• 標本サイズが 2 なら 1 次（直線）で完全にフィット可能
  
• 標本サイズが 3 なら 2 次の多項式（曲線）で完全にフィット可能
  
• 標本サイズが \(n\) なら \(n - 1\) 次の多項式（曲線）で完全にフィット可能
  
• しかし、RDD ではあまり高次の多項式を使うべきではない (Gelman & Imbens 2019)

パラメトリック RD を推定する上で必要な R コード

list() と alist()

list()

・list() を使ってリストデータを作成してみる

list_1 <- list(a = 10 + 2,
             b = 10 - 2,
             c = 10 * 2,
             d = 10 / 2)

• 長さ 4 のリストができた

is.list(list_1)

[1] TRUE

length(list_1)

[1] 4

• リストの中身を確認

list_1

$a
[1] 12

$b
[1] 8

$c
[1] 20

$d
[1] 5

• list() の場合は、リストを作るときに指定した「式の結果」がリストの中に保存されている
  
• では次に alist() を使ってリストを作ってみる

alist()

list_2 <- alist(a = 10 + 2,
             b = 10 - 2,
             c = 10 * 2,
             d = 10 / 2)

• 長さ 4 のリストができた

is.list(list_2)

[1] TRUE

length(list_2)

[1] 4

• リストの中身を確認

list_2

$a
10 + 2

$b
10 - 2

$c
10 * 2

$d
10/2

• alist() を使うと「式の結果」ではなく「式そのもの」が保存されていることがわかる
  

6.3 補講授業の処置効果（パラメトリック RD）

• ここではパラメトリック RD を使って補講授業の処置効果を推定してみる

• \(c\) はカットオフ（＝ 60点）
  

• 境界線を事前に決める割当変数 \(R\) の値 \(R_c\) = カットオフ(カットオフ)
  

• 大学の授業で前期と後期に 2 つの試験を受けるとする
  

• どちらの試験も100点満点
  

• 前期試験で60点未満の生徒は夏休みの間に補講を受けた（60点は含まない）
  

• 前期試験で60点以上の生徒は補講を受けない（60点を含む）

• 「\(Y\) を \(X\) に回帰する」とき、パラメトリック回帰では次のような回帰関数を想定する

\[Y_i = f(X_i)\]

• しかし \(f(X)\) の関数は無数にある
  
• ここでは、補講効果の回帰関数として以下の 4 つの関数形を考えて、まずパラメトリック回帰をやってみる

推定に使う関数形

\(1. Y_i = α + ρD_i + βR_i + ε_i\)
\(2. Y_i = α + ρD_i + βR_i + γ(D_i×R_i) + ε_i\)
\(3. Y_i = α + ρD_i + β_1R_i + β_2R_i^2 + ε_i\)
\(4. Y_i = α + ρD_i + β_1R_i + β_2R_i^2+ γ_1(D_i×R_i)+ γ_2(D_i×R_i^2) + ε_i\)

• 割当変数 \(R\) を中心化すれば、どの関数形においても \(ρ\) （ロー）の推定値が処置効果の推定値となる

• 200人の架空生徒データ (exam.score.csv)を使って具体的に考えてみる
  

• exam.score.csv をダウンロードし、RProject フォルダー内の data フォルダに保存する
  

• データを読み取る

df0 <- read_csv("data/exam.score.csv")

• 割当変数 (R) を中心化した変数 RC を作る
  
• RC の値を 2 乗した RC_sq を作る

df0 <- df0 |> 
  mutate(RC = R - 60) |>  # R を中心化した RC の作成
  mutate(RC_sq = RC^2)     # RC を2乗した RC_sqの作成

• データの中身を確認する

stargazer(as.data.frame(df0), type = "html")

Statistic	N	Mean	St. Dev.	Min	Max
R	200	60.690	12.202	27	92
D	200	0.480	0.501	0	1
Y	200	67.755	10.767	41	100
y0t	200	65.430	11.184	38	100
y1t	200	73.605	12.492	41	100
ITE	200	8.175	10.509	-24	53
RC	200	0.690	12.202	-33	32
RC_sq	200	148.630	223.198	0	1,089

ここでは観察可能な 5 つの変数 (\(R, D, Y, RC, RC_{sq}\)) だけを使って分析する

割当変数 \(R\)	前期試験点数
処置変数 \(D\)	0 = 「補講なし」、 1 = 「補講あり」
結果変数 \(Y\)	後期試験点数
中心化した変数 \(RC\)	\(R\)を中心化した変数
中心化した変数の 2 乗値 \(RC{sq}\)	\(R\)を中心化した変数の 2 乗

• ①から④の回帰式を lm 推定定する formula としてまとめて M_1 と名前を付ける

M_1 <- alist(
  model1 = Y ~ D + RC,
  model2 = Y ~ D * RC,
  model3 = Y ~ D + RC + RC_sq,
  model4 = Y ~ D * RC + D * RC_sq) |> 
  enframe(name = "modelname", value = "formula")

• M_1 の中身を確認

M_1

modelname	formula
model1	Y ~ D + RC
model2	Y ~ D * RC
model3	Y ~ D + RC + RC_sq
model4	Y ~ D * RC + D * RC_sq

• 4 つの回帰式を、一挙に推定
  
• 推定結果を rdd_1 として保存

rdd_1 <- M_1 |> 
  mutate(model = map(.x = formula, .f = lm, data = df0)) |> 
  mutate(pred = map(.x = model, .f = predict),   # 予測値を計算する
         result = map(.x = model, .f = tidy))

names(rdd_1)

[1] "modelname" "formula"   "model"     "pred"      "result"   

• rdd_1 に保存されている結果を1つずつ確認する

6.3.1 model1 の推定結果

• model1 による推定結果を図示してみる

p1_rdd_1 <- df0 |> 
  mutate(pred = rdd_1$pred[[1]]) |> 
  ggplot(aes(x = R, color = as.factor(D))) +
  geom_point(aes(y = Y)) +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 60, color = "gray")+
  labs(x = "前期試験結果", y = "後期試験結果") +
  theme(legend.position = "none") +
  scale_color_manual(values = c("blue", "red"))
plot(p1_rdd_1)

• 前期試験点数が 60 点の時に、後期試験点数のジャンプがある
  
• この分析は、処置変数と割当変数のみを説明変数として利用する、最も単純な パラメトリック RD
  
• そのため、予測値は「直線」であり、前期試験点数が 60 点未満と 60点以上について直線の傾きが同じ
  
• このときの処置効果の推定値と標準誤差を計算してみる

rdd_1 |> 
  unnest(cols = result) |> 
  filter(modelname == "model1", term == "D") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)

term	estimate	std.error	p.value
D	12.62	2.04	0

6.3.2 model2 の推定結果

• 処置変数と割当変数の交差項をもつモデル 2 の結果を図示してみる

p2_rdd_1 <- df0 |> 
  mutate(pred = rdd_1$pred[[2]]) |> 
  ggplot(aes(x = R, color = as.factor(D))) +
  geom_point(aes(y = Y)) +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 60, color = "gray") +
  labs(x = "前期試験結果", y = "後期試験結果") +
  theme(legend.position = "none") +
  scale_color_manual(values = c("blue", "red"))

plot(p2_rdd_1)

• 交差項を使った結果、前期試験点数が 60 点未満と 60点以上の直線の傾きが少しだけ異なる
  
• この関数形でもやはり、前期試験点数が 60 点で後期試験点数のジャンプが確認できる
  
• 処置効果の推定値と標準誤差を計算してみる

rdd_1 |> 
  unnest(cols = result) |> 
  filter(modelname == "model2", term == "D") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)

term	estimate	std.error	p.value
D	12.57	2.04	0

6.3.3 model3 の推定結果

• モデル 1 に割当変数の 2 乗項を加えたモデル3の推定結果を図示してみる

p3_rdd_1 <- df0 |> 
  mutate(pred = rdd_1$pred[[3]]) |> 
  ggplot(aes(x = R, color = as.factor(D))) +
  geom_point(aes(y = Y)) +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 60, color = "gray") +
  labs(x = "前期試験結果", y = "後期試験結果") +
  theme(legend.position = "none")+
  scale_color_manual(values = c("blue", "red"))
plot(p3_rdd_1)

• 2 乗項を使った結果、予測値が直線ではなく「曲線」として推定された
  
• 交差項がないので、曲線は前期試験点数が 60 点未満と 60点以上で同じと仮定
  
• この関数形でもやはり、前期試験点数が 60 点に後期試験点数のジャンプが確認できる
  
• 処置効果の推定値と標準誤差を計算してみる

rdd_1 |> 
  unnest(cols = result) |> 
  filter(modelname == "model3", term == "D") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)

term	estimate	std.error	p.value
D	12.5	2.04	0

6.3.4 model4 の推定結果

• 最後に、交差項と 2 乗項を含むモデル4の推定結果を可視化してみる
  
• 比較のため、モデル1で推定した直線も点線として表示する　　

p4_rdd_1 <- df0 |> 
  mutate(pred0 = rdd_1$pred[[1]],
         pred = rdd_1$pred[[4]]) |> 
  ggplot(aes(x = R, color = as.factor(D))) +
  geom_point(aes(y = Y)) +
  geom_line(aes(y = pred0), linetype = "dashed") +
  geom_line(aes(y = pred)) +
  geom_vline(xintercept = 60, color = "gray") +
  labs(x = "前期試験結果", y = "後期試験結果") +
  theme(legend.position = "none") +
  scale_color_manual(values = c("blue", "red"))
plot(p4_rdd_1)

• 交差項と 2 乗項の両方を使ったので、予測値が曲線になり、前期試験点数が 60 点未満と 60点以上では異なる曲線が描かれている
  
• ここでもやはり、試験点数が 60 点に後期試験点数のジャンプがある
  
• 処置効果の推定値と標準誤差を計算してみる

rdd_1 |> 
  unnest(cols = result) |> 
  filter(modelname == "model4", term == "D") |> 
  select(term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)

term	estimate	std.error	p.value
D	11.37	2.92	0

6.4 モデルの推定結果比較

rdd_1 |> 
  unnest(cols = result) |> 
  filter(term == "D") |> 
  select(modelname, term, estimate, std.error, p.value) |> 
  knitr::kable(digits = 2)

modelname	term	estimate	std.error	p.value
model1	D	12.62	2.04	0
model2	D	12.57	2.04	0
model3	D	12.50	2.04	0
model4	D	11.37	2.92	0

• モデル 1 から 3 までの結果はほぼ同じ
  
• 最も複雑なモデル 4 の推定値だけが小さい
  
• モデル 1 から 3 か、モデル 4 のうち、少なくとも一方の推定にバイアスが生じていると思われる
  
• 私たちはどちらの関数形が正しいかを知らない
  
• バイアスがどちらにあるかわからない
  
• どちらにもバイアスがあるという可能性もある
  
• この例でも明らかなように、パラメトリック RD の推定結果は関数形に依存する
  
• 誤った関数形を使うと、誤った推定結果を得る

→　大事なこと

• 社会科学（経済学、経営学、心理学、政治学など）の先行研究成果を踏まえ、きちんとした理論に基づいた妥当な関数形（＝回帰式）を選ぶ
  
• 理論的に妥当な複数の関数形によって推定した結果がどれも実質的に同じなら
  →　推定結果に対する信頼性が高まる

7. ノンパラメトリック RD（補講授業）

7.1 ノンパラメトリック RD とは

• パラメトリック RDでは、回帰式が特定の形であると想定し、割当変数 \(R\) が最小値から最大値までの範囲全体について回帰式の予測値を計算した
  

• それに対し、ノンパラメトリック RD では回帰式に含めるべき変数の選択や関数形を考えなくてよい
  

• モデリングに関する様々な制約を回避できる
  

• ノンパラメトリック RD では推定の対象となる範囲をカットオフ値の周辺に限定し、処置群と統制群の平均値を比較
  

• 閾値の周辺データだけを使って、交差項を含む共分散分析を行う
  

• 推定の対象となる範囲を決める幅のことをバンド幅 (bandwidth) と呼ぶ
  

• バンド幅の選び方によって推定値が変わる
  

• ノンパラメトリック RD ではバンド幅の選び方が重要
  

• 処置群と統制群が比較可能なのはカットオフ値付近だけ

• バンド幅が狭いほど観察数が少なくなり、局所的な無作為割り付け (RCT) が成立する可能性が高まる
  

• RCT の確率が高まれば、推定値の偏り（バイアス）は小さくなる

• しかし、バンド幅を狭くすると、推定の対象となるサンプルサイズ（観測数）が小さくなるので、推定が不安定になる
  →　標準誤差が大きくなり、推定値の精度が落ちる
  →　推定の信頼性が落ちる

このバイアスと標準誤差のトレードオフのバランスを考えて、バンド幅を選ぶ必要がある

バンド幅・推定値の精度・推定値の偏りの関係

バンド幅	観測数	RTCの可能性	推定値の偏り	標準誤差	推定値の精度
狭い	少ない	高い	小さい	大きい	低い
広い	多い	低い	大きい	小さい	高い

• 「バンド幅が広すぎる」　→　閾値の辺の値を比較できなくなる
  →　回帰不連続デザインのい意義がなくなる

• バンド幅の決め方に関しては Imbens & Kalyanaraman (2011) を参照

• データ全体での割当変数である前期試験結果 \(R\) の分布を確認する

summary(df0$R)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  27.00   54.00   61.00   60.69   67.25   92.00 

• 前期試験結果 \(R\) には最低点（27点）から最高点（92点）までが含まれている
  
• カットオフは 60点なので、分析対象を 60 点周辺のデータに絞り込み、処置効果を推定してみる

7.2 ノンパラメトリック RD の推定

• データの中身を確認する

stargazer(as.data.frame(df0), type = "html")

Statistic	N	Mean	St. Dev.	Min	Max
R	200	60.690	12.202	27	92
D	200	0.480	0.501	0	1
Y	200	67.755	10.767	41	100
y0t	200	65.430	11.184	38	100
y1t	200	73.605	12.492	41	100
ITE	200	8.175	10.509	-24	53
RC	200	0.690	12.202	-33	32
RC_sq	200	148.630	223.198	0	1,089

• ここで使う変数は観察可能な次の 2 つの変数　　

\(R\)	前期試験点数（割当変数）	観察可能
\(Y\)	後期試験点数	観察可能

7.2.1 散布図と RD プロットを描く

• 最初に、割り付け変数 \(R\) を \(x\) 軸、結果変数 \(Y\) を \(y\) 軸とした散布図と RDプロットを描く

df0 |> 
  ggplot(aes(Y, R)) +
  geom_point() +
  stat_smooth()  # 線形ではなく LOWESS 曲線を描く  

• データをスムーズ化し、\(Y\) にジャンプがあるかどうか確認しやすい RDプロットを描いてみる

library(rdrobust)

rdplot(df0$Y, df0$R, c = 60)

[1] "Mass points detected in the running variable."

• x = 60 のカットオフでジャンプが確認できる

7.2.2 rdrobust パッケージを使った分析

• RDD を実行する前に必要な仮定の確認は終わっているので、ここでは Step 6-2 を実行する

Step 1: 処置群と統制群を分けるルールはあるか？	5.1 で確認済み
Step 2: RDDデザインは sharp か fuzzy か？	5.2 で確認済み
Step 3: カットオフ付近での結果変数 \(Y\) の非連続性	5.3 で確認済み
Step 4: カットオフ付近での割当変数 \(R\) の非連続性	5.4 で確認済み
Step 5: 割当変数 \(R\) の値を自分では選べない	5.5 で確認済み
Step 6-1: 処置効果の推定（パラメトリック）	5.2 で実行済み
Step 6-2: 処置効果の推定（ノンパラメトリック）
Step 7: 推定結果を比較する

• ここでは次の複数のバンド幅をモデリングして推定する
  ・IK バンド幅
  ・CER 最適バンド幅
  ・MSE 最適化バンド幅
  
• まず、Imbens and Kalyanarman (2012) が推奨するバンド幅 IKバンド幅 を確認してみる

IK バンド幅

IKband <- rdbwselect_2014(df0$Y, 
                          df0$R, 
                          c = 60, 
                          bwselect = "IK") # 
IKband$bws

            h        b
[1,] 5.731735 4.886543

• Imbens and Kalyanarman (2012) は 5.731735 が最適なバンド幅だと推奨している
  
• では、この推奨されたバンド幅を使って、補講授業の処置効果を推定してみる

model_IK <- rdrobust(df0$Y, 
                     df0$R, 
                     c = 60, 
                     h = IKband$bws[1,1], 
                     kernel = "triangular", # カーネルは三角形がデフォルト
                     all = TRUE)            # 他の選択肢は epanechnikov と uniform  

[1] "Mass points detected in the running variable."

summary(model_IK)

補講授業の局所的処置効果 (LATE) パラメトリック(バンド幅 IK = 5.732)

• IK によって推奨されたバンド幅 (5.732) の時の推定値は 8.4

CER 最適バンド幅

• 次に CER が推奨するバンド幅を使って、補講授業の処置効果を推定してみる

model_cer <- rdrobust(df0$Y, 
                      df0$R, 
                      c = 60, 
                      bwselect = "cerrd", 
                      all = TRUE)

[1] "Mass points detected in the running variable."

summary(model_cer)

Call: rdrobust

Number of Obs.                  200
BW type                       cerrd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                   96          104
Eff. Number of Obs.              54           54
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   7.087        7.087
BW bias (b)                  14.240       14.240
rho (h/b)                     0.498        0.498
Unique Obs.                      26           28

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional   -11.147     3.933    -2.834     0.005   [-18.856 , -3.438]    
Bias-Corrected   -11.044     3.933    -2.808     0.005   [-18.753 , -3.334]    
        Robust   -11.044     4.493    -2.458     0.014   [-19.850 , -2.237]    
=============================================================================

MSE 最適化バンド幅

• 最後にMSE が推奨するバンド幅を使って、補講授業の処置効果を推定してみる

model_mse <- rdrobust(df0$Y, df0$R, c = 60, bwselect = "mserd", all = TRUE)

[1] "Mass points detected in the running variable."

summary(model_mse)

Call: rdrobust

Number of Obs.                  200
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                   96          104
Eff. Number of Obs.              60           64
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   9.236        9.236
BW bias (b)                  14.240       14.240
rho (h/b)                     0.649        0.649
Unique Obs.                      26           28

=============================================================================
        Method     Coef. Std. Err.         z     P>|z|      [ 95% C.I. ]       
=============================================================================
  Conventional   -12.306     3.502    -3.514     0.000   [-19.170 , -5.443]    
Bias-Corrected   -12.114     3.502    -3.459     0.001   [-18.978 , -5.250]    
        Robust   -12.114     4.256    -2.847     0.004   [-20.455 , -3.773]    
=============================================================================

7.2.3 カーネル密度推定

• ノンパラメトリックなカーネル密度推定 (kernel density estimate) とは、ヒストグラムをスムーズ化したもの
  
• 観察されたデータから確率密度関数を推定する方法
  
• カーネル関数には次の 4 種類がある

layout(matrix(1:4, 2, 2, byrow=TRUE))
plot(density(0, bw=1, kernel = "gaussian"))
plot(density(0, bw=1, kernel = "rectangular"))
plot(density(0, bw=1, kernel = "triangular"))
plot(density(0, bw=1, kernel = "epanechnikov"))

4 種類のカーネル関数

• 最も一般的なのがガウス関数（＝正規分布）
  
• 閾値 \(c\)、バンド幅を決める定数 \(h\) とすると、\(c±h\) がバンド幅
  
• このバンド幅の中に入るデータ \(x_i\) に重み付けをする方法
  
• それぞれの関数の特徴は次のとおり
  
• いずれのカーネル関数においても、閾値 \(c\) の位置で重みが最大　　

ガウス (Gaussian)	最も一般的なガウス関数（＝正規分布）
矩形・くけい (rectangular, uniform)	バンド幅内の全ての値の重みが等しい
三角形 (traiangular)・・・rddensity()関数ではデフォルト	閾値から遠ざかるにつれて、左右対称で線形的に重みが小さくなる
エパネチニコフ (Epanechnikov)	バンド幅内の観測値に対して曲線を描くように重みが小さくなる

• それぞれの関数の重み付けの方法を可視化するとこのようになる

3 種類のカーネル関数の重み付け方法

• 三角形 (triangular) の場合、真ん中の尖っているところが黄色で最も重く、左右に離れるにつれて軽くなる
  
• 矩形 (uniform) の場合、バンド幅内のデータの重みは等しい（＝全て同じ紫色）ことがわかる

どのカーネルを選択するかはあまり大きな問題ではない

• カーネルは次のいずれかの指定が可能

・kernel = "triangular"・・・rdrobust() 関数ではデフォルト設定
・kernel = "epanechnikov"
・kernel = "uniform"

7.2.4 分析結果のまとめ3

• 「補講の処置効果」に関して、ここまで実施した分析結果をまとめてみる
  
• 「平均処置効果 (ATE)」 と「局所的な平均処置効果 (LATE) 」は異なる
  ・私たちが設定した ATE の真値は 8.175
  ・ATE の求め方は 3.3 平均処置効果 (ATE)を参照
  ・ここで私たちが推定しようとしているのは LATE
  ・平均処置効果 (ATE) を推定するためには、\(X\) の定義域全体に対して「外挿」しなければならないため、ATE の推定は諦めた
  ・大域的な推定方法の代わりに、閾値における局所的な平均処置効果 (LATE) を推定することにした

・私たちがもともと設定した LATE の真値は 8.08

私たちが設定した LATE の真値 ・理論的には、閾値における局所的な平均処置効果 (LATE) は、共変量 \(X_i\) が閾値 \(c\) の値をとるときの \(Y_i(1)\) と \(Y_i(1)\) の差の期待値

\[E[Y_i(1)-Y_i(0)|R_i = c] = E[Y_i(1)|R_i = c] - E[Y_i(0)|R_i = c]\]

・ここでは潜在的結果変数 \(y0t, y1t\) が観察可能と想定しているので、補講を「受けた場合」と「受けなかった場合」それぞれの回帰式は次のように設定できた

\[Y_i(0) = 35.9 + 0.486R \\ Y_i(1) = 35.4 + 0.629R\]　

• 具体的には、ここでは共変量 \(X_i\)（前期試験点数 \(R_i\)）が閾値 \(c = 60\) の値をとる時の後期試験点数 \(Y_i(1)\) と \(Y_i(0)\) の差の期待値のこと
  
• 従って、次のように計算できる

\[E[Y_i(1)-Y_i(0)|R_i = c] = E[Y_i(1)|R_i = c] - E[Y_i(0)|R_i = c]　\\ = E[35.4 + 0.629R|R_i=60]-E[35.9+0.486R|R_i=60] \\ = 35.4 + 0.629*60 - 35.9 + 0.486 * 60 \\ = 73.14 - 65.06 \\= 8.08\]

推定結果（補講授業の処置効果）

手法	推定値	標準誤差	p-value	95%信頼区間
単純な前後比較	7.065
後期試験点数の平均値比較	-1.15
交絡を含めた重回帰分析	12.6
閾値前後の観測値の比較	1.47
閾値±5 の観測値の比較	10.57
パラメトリックmodel1	12.62	2.04	0
パラメトリックmodel2	12.57	2.04	0
パラメトリックmodel3	12.50	2.04	0
パラメトリックmodel4	11.37	2.92	0
ノンパラIK バンド幅	8.384	4.336	0.053	-16.884, 0.115
ノンパラCER	11.147	3.933	0.005	-18.856 , -3.438
ノンパラCERバイアス修正	11.044	3.933	0.005	-18.753 , -3.334
ノンパラCERロバスト	11.044	4.493	0.014	-19.850 , -2.237

パラメトリック model1 の推定結果： LATE = 12.62

• 私たちがもともと設定した LATE の真値 (8.08) に最も近いのはノンパラメトリック IK バンド幅 (8.384)
  
• ノンパラメトリック RD の局所的処置効果 (LATE) の推定値は、パラメトリック RD の推定結果よりも少しだけ LATE の真値に近い
  
• パラメトリック RD のモデルを見ると、過大推定バイアスがあったのではなないかと推測できる
  
• 1 つのデータを多角的に検討することが重要
  
• 関数形がわからないときは「境界線付近」以外のサンプルを使うことは諦めて（パラメトリック RD は諦めて）ノンパラメトリック RD でバイアスの小さな推定値を得ることを目指すほうが良い
  
• 結論としては、例えば、バイアスが小さなノンパラメトリック RD で得られた 4 つの推定値の平均を示すことも可能（すべて統計的に有意）

(8.384 + 11.147 + 11.044 + 11.044)/4

[1] 10.40475

補講授業の処置効果の推定結果 補講授業を受けると後期試験結果を 10点程度高める

7.2.5. 共変量がある場合:

• RDDデザインによって局所的な無作為割り付けができている
  →　局所的には実験研究
  →　共変量をモデルに含める必要はないのでは？
  
• 変量をモデルに含めることで、推定の精度を向上させることができる
  
• RDD における共変量は、交絡変数とは異なり、関数系を気にせず共変量をモデルに追加するだけでよい
  
• RDD では、共変量を追加してもしなくても、パラメータ推定値の一致性に影響はない

Reference

Angrist, Joshua D., and Jörn-Steffen Pischke. 2009. Mostly Harmless Econometrics: An Empiricist’s Companion. Princeton: Princeton University Press.

Guido Imbens and Karthik Kalyanaraman (2012), Optimal Bandwidth Choice for the Regression Discontinuity Estimator Get access Arrow, The Review of Economic Studies, Volume 79.

Andrew Gelman & Guido Imbens (2019), Why High-Order Polynomials Should Not Be Used in Regression Discontinuity Designs, Journal of Business & Economic Statistics, Volume 37

PMAP 8521 • (10) Regression discontinuity I: (2) Drawing lines and measuring gaps

Philipp Leppert, R Tutorial: Regression Discontinuity Design

矢内勇生「KUT計量経済学応用」

高橋将宣『統計的因果推論の理論と実装』、共立出版、2022年

清水昌平『統計的因果探索』、講談社、2017年

宋財泫 (Jaehyun Song)・矢内勇生 (Yuki Yanai)「私たちのR: ベストプラクティスの探究」

浅野正彦・ 矢内勇生『Rによる計量政治学』オーム社、2018年

松林哲也 『政治学と因果推論』、2021年

Winston Chang, R Graphics Cookbook, O’Reilly Media, 2012.

Kosuke Imai, Quantitative Social Science: An Introduction, Princeton University Press, 2017