FuzzyRDDで少人数教室の学力向上効果を検証

# 概要
・少人数教室の学力向上の効果をFuzzyRDDを利用して推論。
・FuzzyRDDとはRDD（回帰不連続デザイン）が適用できるような状況（変数にジャンプが存在する）で、操作変数法を用いることで因果効果を推定する方法。もはや操作変数法の一種と言っても過言ではない。
・操作変数法の仮定を満たすことでバイアスの軽減しながら回帰をすることが可能。
・少人数教室は国語のテストにおいては有意な差でプラスの効果があることがわかった。←後の研究で少人数教室に優位な効果がないことが示されているので、今回の分析結果は鵜呑みしないように

コード：https://github.com/harukicode0/Estimating-the-Effectiveness_of_Small_Classrooms
taskを参照

# 前半:FuzzyRDDとは
下の図のように不連続なジャンプが発生しているようなデータに対して、操作変数法を用いて、因果効果を推定する方法のこと。下の図では60の値前後でジャンプが発生している。

![スクリーンショット 2023-04-10 18.33.05.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/199083ae-d452-4df3-9ead-5f8bd5367888.png)


前半ではRDDと操作変数法を説明し、前半の最後にFuzzyについてまとめる。

## RDD（回帰不連続デザイン）とは
例えば、ある学校で２回のテストを実施する状況を考える。１回目のテストの点数が６０点以下の場合は補講を強制的に受けなければならないとする。この時２回目のテストの点数は、１回目のテストの点数が６０点以下と６０点より大きいグループで異なると考えられる（補講の効果により）。このようなシチュエーションで、補講の効果を推定する方法がRDDである。この時、補講を受けたかどうかを介入変数と呼び、介入するか否かを決めた１回目のテストの点数をRunning Variableと呼ぶ。２回目が目的変数となる。不連続の名前の由来は今回の例で言えば、60点前後でおそらく、テストの点数にジャンプ（不連続）が発生しているところに由来する。

また、ジャンプする状況が60点で必ずしも決まっておらず、確率的な場合FuzzyRDDを利用することになる。ここでいうFuzzyとはジャンプするかしないかが確率的（曖昧）な状況のことを示す。

https://qiita.com/icchi-ioo/items/d5d456c9681817cd19d1

## RDDの３つの前提
①Continuity of Conditional Regression Functions
もし、介入がなければ、連続して同じような傾向が続いていたと仮定すること。

②non-mainpulation
ユーザーなどの変数が自らの意思でがカットオフの前後に移動し、介入グループに入ったり、入らなかったりするということがないこと。例えば、3000円以上購入する場合、割引がつくキャンペーンなどをした場合、自らの意思によりカットオフを超えるユーザーが増える。カットオフ前後でデータの密度に差がある場合、このバイアスが生じている可能性があるので注意。

③カットオフ値周辺で他に影響を与えている処置が存在しないこと。

## 操作変数法とは

以下のようなDAGを仮定したとき、交絡変数Uを観測しなくてもXからYへの因果効果を正しく推定できる手法。

![スクリーンショット 2023-04-11 18.33.54.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/1a78c23b-3a6a-9cb4-606d-bc33bd3483c0.png)

この時、

・Zを操作変数
・Xが介入変数
・Yが目的変数
・Uが未観測の交絡変数

と呼ぶ。

操作変数を用いて、２段階回帰分析を行うことで因果効果を推定できる。正し、操作変数と呼ばれるものはなんでもいい訳ではなく、以下の仮定を満たした変数を操作変数として利用できる。操作変数法の一番難しいところは、操作変数を見つけることにある。


## 操作変数法の仮定

ここでは厳密に操作変数に必要な仮定を、上記DAGを前提に説明する。

①外生性
Cov(Z,U) = 0という仮定。別の言い方をすると操作変数ZとY = b0 + b1 * X + UのUとは関係がないこと。これを数学的に証明することはできない。なぜなら、Uは観測することができない変数だから。故に、この仮定を満たしているかどうかを説明するには、論理やドメイン知識を用いて説明するしかない。

②関連性
Cov(Z,X)≠0。操作変数Zと介入変数Xは関係がある。証明するためにはXを目的変数とした重回帰分析を行う。想定されるモデルは X = b0 + b1 * Z + (その他共変量) + (誤差)。この時、操作変数の回帰係数が０でないことに加えて、t値が3.2以上、F値が１０以上であることが目安。

③除外制約
介入変数Xの値を固定すると、操作変数Zの値をいくら変更してもYの値には影響しないこと。または、操作変数ZはXを通してのみ目的変数Yに影響すること。これも論理やドメイン知識を基に、この仮定を満たしていると説明する。

以上の3つの仮定を満たすとき、操作変数法でバイアスを避けて、因果効果を推定できる。

## 操作変数法がバイアスなく推定できる理由

数式を用いて説明する。

```math
Cov(Z,Y) = Cov(Z, b0 + b1 * X1 + U)  \\
= Cov(Z, b0) + Cov(Z,b1 * X1) + Cov(Z, U) \\
```

この時、b０は定数のなのでCov(Z,b0)=0となる。また、操作変数法の仮定よりCov(Z,U)=0である。したがって、

```math
=b1*Cov(Z,X1)
```

となる。
b１が推定したい値なので、式をb１＝の形に変換する。

```math
b1 = \frac{Cov(Z,Y)}{Cov(Z,X1)}
```

ここで、操作変数法の仮定より、Cov(Z,X1)≠0。以上よりb１を推定できる。

また、上記の式を若干変換する。

```math
b1 = \frac{Cov(Z,Y)}{Cov(Z,X1)} \\
= \frac{\frac{Cov(Z,Y)}{Var(Z)}}{\frac{Cov(Z,X1)}{Var(Z)}} \\
= \frac{\rho}{\Phi}
```

$\rho,\Phi$はそれぞれZをYに回帰した時の回帰係数、ZをX1に回帰した時の回帰係数である。故に２段階回帰でパラメータを推定できる。

## FuzzyRDDとは

以上のRDDと操作変数法の式を抑えた上で、FuzzyRDDを説明する。FuzzyRDDとはカットオフでのジャンプが確率的な状況において、操作変数法を用いて因果効果を推定する手法である。

他にも以下の図のような、ジャンプが確率的（曖昧）に複数発生しているようなデータを取り扱う際にも利用できる。

![スクリーンショット 2023-04-11 19.22.41.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/69f65bef-9344-0abe-92bb-8698cafaa34e.png)

手順は以下の通り。
①曖昧なジャンプまたは、ジャンプが複数回発生しているようなデータを発見する
②Running Varibaleに対して、なんらかの操作変数を発見する。ジャンプが発生している法則などに着目すると発見しやすい。
③見つけた操作変数に対して、操作変数であるための3つの仮定を満たしているか確認する。また、3つの仮定に加え、単調性（Defirが存在しない）という仮定とRDDのnon-mainpulationを満たす必要がある。
④操作変数法(操作変数を用いた２段階回帰)を用いて、因果効果を推定。

# 後半
今回はイスラエルの小学校のデータを利用して、少人数教室と国語の成績の関係をFuzzyRDDで効果検証する。

## FuzzyRDDを少人数教室問題で利用するモチベーション

イスラエルの小学校では１クラスあたりの人数を４０人以下に抑えるというルールを適用している。このルールをマイモデニウス・ルールと呼ぶ。

つまり入学者数が４０人の場合、１クラスあたりの生徒の数は４０人だが、4１人となった場合、１クラスあたりの人数は２０人と２１人となる。ここで、入学者数が４０と４１人の場合で、１クラスあたりの人数に急なジャンプが生まれているためRDDを用いるモチベーションが生まれる。

実際に、入学者数と１クラスあたりの人数（今回の目的変数）の散布図を確認する。赤い点線がマイモデニウス・ルールに従った理論的に決まる１クラスあたりの人数で、青い点が実際の値。

![スクリーンショット 2023-04-11 19.37.33.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/49048022-2399-2382-fa7b-53e47ef221fc.png)

すると、必ずしも理論値通りに１クラスあたりの人数が決まっている訳ではなく、Fuzzyであることがわかる。ここでFuzzyRDDを用いる動機が生まれる。なぜならば、マイモデニウス・ルールを数式化し、操作変数とすれば、バイアスを避けながら因果効果(ここでは少人数教室が国語の成績に与える影響)を推定することができるからだ。

故に、まずマイモデニウス・ルールを数式化する。

## 今回の操作変数：１クラスあたりの生徒数

１クラスあたりの人数（クラスサイズ）はマイモデニウス・ルールで決まる。マイモデニウス・ルールを数式下したものは下記の通り。

```math
クラスサイズの理論値(操作変数)：Z = \frac{入学者数}{int(\frac{入学者数 - 1}{40}) + 1}
```

ここで,intは小数点を切り捨てする関数。

このように、操作変数を作った訳だが、この操作変数が実際に仮定を満たしているのか確認する必要がある。

## 今回の操作変数が仮定を満たしているか確認

①外生性について
Cov（Z,U）=0を満たしているか論理を持って説明する。操作変数Zは入学者数によって自動的に決まるものなので、そのほかの変数となんらかの関連を持っているとは考えにくい。また、目的変数である国語の成績とも関係なく決まっているので、目的変数とも同様に関係がないものだと思われる。

②関連性について
Cov(Z,X)≠0。先ほどの散布図を確認したように、操作変数Zと実際のクラスサイズには関連があるように見える。念の為、共変量を含めた重回帰分析でt値とF値を確認する。

![スクリーンショット 2023-04-11 19.58.18.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/d3b56658-8f6a-71af-cad8-feb0fc74fa94.png)

実際のクラスサイズ：この重回帰分析の目的変数
func2：操作変数
tipuach：共変量、生徒の学校以外の教育資源（家庭環境や住んでいる場所など）がどの程度悪いのか示す指標。
c_size：共変量、入学者数。

F値1036、t値32より基準を満たしていることがわかる。

③除外制約について
介入変数である１クラスあたりの人数を固定した時、操作変数をいくら変化させても、学業成績には影響しないことは論理的にわかる。故に、除外制約を満たしている。

④単調性
被験者を介入に割り当てたときに非介入の行動をし、非介入に割り当てたときに介入の行動をする被験者のことを天邪鬼という。この天邪鬼が存在しないという仮定。少人数教室の文脈で再説明すると、入学者数が増えると、１クラスあたりの人数を大きくし、入学者数が減ると、１クラスあたりの人数を小さく学校が存在しないことを仮定するもの。学校はマイモデニウス・ルールに従うので、今回の単調性は満たしていると考えらえる。

## ２段階回帰式の構造
少人数教室の効果を予測するための２段階回帰式は以下のようになる

```math
１段階目：classize = b1 * func2
```

・classizeは１クラスあたりの人数
・func2はマイモデニウス・ルールを数式化したもの。

```math
２段階目：avgverb = b0 + b1 * c_size + b2 * tipuach + b3 * classizeの予測値
```

・avgverbが国語の点数
・c_sizeが入学者数
・tipuachが学校以外の教育資源（親や塾等）のビハインド具合をパーセンテージにしたもの。数字が大きいほど、学校以外の教育資源に乏しい。

以上の２段階回帰式を用いて、classizeの効果（少人数教室が国語の成績に与える影響）を推定する。

補足だが、１段階目の回帰式に２段階目の回帰式で利用した共変量を入れる場合もある。今回だとc_sizeとtipuachを１段階目の共変量に入れて分析することもある。

## 結果と解釈

![スクリーンショット 2023-04-11 23.45.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/79795b45-6739-f70d-41d7-10a3b0609ed9.png)

classizeがマイナスの値で有意差あり。故に、classizeを小さくするとプラスの効果があることが示された。つまり、少人数教室には国語の成績にプラスの効果があることが示された。

また、カットオフが発生しているデータ±５前後(36~45,76~85,116~125)のデータのみを使った操作変数法の結果は以下の通りとなる。

![スクリーンショット 2023-04-11 23.55.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2707093/1ff89991-0748-d8a6-9af3-0fb2a0e18db6.png)


数値は違うが全データを利用した分析と同様の結果が示された。



# 参考文献
データセット  

https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/XRSUJU  

サイト  
dtaファイルの読み込み  

https://amatosan.hatenablog.com/entry/2016/02/15/003513  

はてなブログ、少人数教室の効果推定  

https://ill-identified.hatenablog.com/entry/2014/12/18/012338  

操作変数の良し悪しを判断する方法について記述がある。  

https://www.nishiyama.kier.kyoto-u.ac.jp/2018/jugyochukei3.pdf  

pythonで２段階回帰分析  

https://py4etrics.github.io/19_IV2SLS.html

矢内先生のスライド  

http://yukiyanai.github.io/jp/classes/econometrics2/contents/slides/metrics2_topic09_slides.pdf  

日本、回帰不連続デザイン  

https://core.ac.uk/reader/229320690  

統計的因果推論と理論の実装　

https://www.amazon.co.jp/%E7%B5%B1%E8%A8%88%E7%9A%84%E5%9B%A0%E6%9E%9C%E6%8E%A8%E8%AB%96%E3%81%AE%E7%90%86%E8%AB%96%E3%81%A8%E5%AE%9F%E8%A3%85-Wonderful-R-%E9%AB%98%E6%A9%8B-%E5%B0%86%E5%AE%9C/dp/4320112458/ref=sr_1_2?__mk_ja_JP=%E3%82%AB%E3%82%BF%E3%82%AB%E3%83%8A&crid=3JN92GUK9DB1P&keywords=%E7%B5%B1%E8%A8%88%E7%9A%84%E5%9B%A0%E6%9E%9C%E6%8E%A8%E8%AB%96&qid=1681117674&sprefix=%E7%B5%B1%E8%A8%88%E7%9A%84%E5%9B%A0%E6%9E%9C%E6%8E%A8%E8%AB%96%2Caps%2C186&sr=8-2