Regression Discontinuity

author : kei ikegami


RDは計量経済学で用いられる手法の一つです。

人道的、社会的制約のためにランダム化比較実験が行えない状況を想定します。

そのような状況で因果効果を算出するために、処置の割り当てが擬似的にランダムとなる状況を利用する一連の手法が存在します。

RDはその代表的な手法であり、ImbensやAngristの手により近年急速に発展しました。

この記事では、RDの基本的な発想を解説しつつ、Juliaを使って実際にRDによる因果効果の測定を行ってみます。


RDの理論的側面を解説します。

RDには_Sharp RD_と_Fazzy RD_の二種類が存在します。ここではSharp RDのみ解説します。(Fuzzyはそれがわかれば簡単なので、教科書を読んでください) 特に、angristなどでさっと流されるband width selectionについてそのやり方を書いたので読んどいてください。

RD、その前に

次のように具体的な状況を考えます。

経済学者の河原崎は東大に行くことが学力向上に寄与するかどうかに関心を持っています。 彼は学力の客観的な指標の作成に成功しており、そのテストをPikaと名付けました。Pikaは特殊なテストで、得点分布は正規分布に従うことがわかっています。

この夏、河原崎は東大とW大学の学生に対して大規模な学力調査を行い、各大学から大学4年生1000人分のPika点数データを入手しました。

そのデータは以下のような構造をしていました。


In [2]:
# パッケージインストール
using Gadfly
using DataFrames
using GLM

In [3]:
srand(22)
n = 1000
todai = Array{Float64}([i + 1 for i in randn(n)])
wuni = Array{Float64}([i - 0.5 for i in randn(n)])
plot(
    layer(x = todai, Geom.histogram),
    layer(x = wuni, Geom.histogram, Theme(default_color=colorant"red")),
    Guide.xlabel("Pika Score"), Guide.ylabel("Population"), Guide.title("Pika Score")
)


Out[3]:
Pika Score -20 -15 -10 -5 0 5 10 15 20 -15.0 -14.5 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0 -20 -10 0 10 20 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 -40 -30 -20 -10 0 10 20 30 40 50 60 70 -30 -29 -28 -27 -26 -25 -24 -23 -22 -21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 -30 0 30 60 -30 -28 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 52 54 56 58 60 Population Pika Score

赤がw大学のpika scoreで、青が東大のpika scoreです。

河原崎はこのデータを見て、「東大に行くと学力が向上しそうだ!」と思いました。

彼はこの仮説を検証するために、学生時代に習った回帰分析を使ってpika scoreを東大生のダミー変数と定数項に回帰してみることにしました。

回帰式は下のようになります。

$ pika_{i} = \alpha + {\beta}*todai_{i} + \epsilon_{i}$

結果は以下のようでした。


In [4]:
srand(22)
dummy = ones(2n)
dummy[1001:end] = 0
pika = [todai; wuni]
X = reshape([ones(2n) dummy], 2n, 2)
slope_coef = (inv(X' * X) * X') * pika
coef = ("alpha", "beta")
for i in 1:length(coef)
    c = coef[i]
    d = slope_coef[i]
    println("$c = $d")
end


alpha = -0.5207758568352845
beta = 1.5139362205338938

t値等も合わせて結果を表記すると以下のようになります。

上の回帰分析

「東大ダミーは有位に正だ!やっぱり東大に行くと学力が向上するんだ!」と河原崎は思いました。

しかし、経営学者の正田がこれにまったをかけます。

「東大生はもとからW大学の学生よりも頭がいいんだから、その結果は東大に入ったことによる学力向上が測定されているとは言えないんじゃない?」

たしかにその通りです。 しかし、河原崎はこれ以上の良い方法を知りません。

どのようにすれば彼らは東大に入ることの効果を識別することができるでしょう?

RD

2013年度東京大学文科Ⅲ類の合格最低点は347.2111点でした。つまり、347.211点を0.001点でも上回った受験者は東大に合格し、0.001点でも下回った受験者は東大には行くことはできませんでした。

RDはこの事実を使って、347.211点をわずかに超えて合格した人と、わずかに下回って不合格になった人とで従属変数を比較する手法です。

今、A君の合計得点が348点で、B君の合計得点が347点だとしましょう。

A君はめでたく東大に合格できましたが、B君は残念ながら不合格、W大学に通うことになりました。

さて、この二人について先ほど正田が述べた理屈は当てはまるでしょうか?

二人の得点差はわずかに1点。これはセンター試験の国語で1問落としたとか、世界史でたまたま直前に見てたところが出たとか、そのぐらいの偶然で十分に生まれる差です。

しかし、このわずか1点によって彼らが後の4年間を過ごす大学が変わってくるわけです。

ここで東大に入ることの学力向上への効果を見たいわけですから、東大入学を処置と考えると、 ほとんど能力に差のない二人に処置の割り当てがランダムに行われる、という状況が存在していることがわかります。

これがまさに自然実験であり、擬似的なランダム化比較実験が行われていると解釈できる状況です。

Angristでやった通り、ランダム化が行われていれば処置は潜在変数と独立となるため、処置群と対照群とで従属変数の平均値を比較すれば、それを平均的な因果効果として算出することができます。

今回の文脈で言えば、A君とB君のPika Scoreの点を比べれば、それは東大で勉強したことによる因果効果を算出している、ということです。

以下で実際のデータを模して作った擬似データを使ってこの状況をグラフにします。


In [5]:
# 前期試験の得点分布作成
srand(22)
cut = 347.211
n_1 = 500
n_2 = 300
n_3 = 150
n_4 = 50
n_array = [n_1, n_2, n_3, n_4]

pop_todai = Array(Float64)
for (i, n) in enumerate(n_array)
    pop_todai = [pop_todai ; [j*25 + cut + 25*(i-1) for j in rand(n)]]
end
pop_todai = Array{Float64}(pop_todai[2:end])

pop_w = Array(Float64)
for (i, n) in enumerate(n_array)
    pop_w = [pop_w ; [cut - j*25 - 25*(i-1) for j in rand(n)]]
end
pop_w = Array{Float64}(pop_w[2:end])

plot(
layer(x = pop_todai, Geom.histogram),
layer(x = pop_w, Geom.histogram, Theme(default_color=colorant"red")),
Guide.xlabel("Entrance Exam Score"), Guide.ylabel("Population")
)


Out[5]:
Entrance Exam Score -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650 660 670 680 690 700 -250 0 250 500 750 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700 -30 -25 -20 -15 -10 -5 0 5 10 15 20 25 30 35 40 45 50 55 -25 -24 -23 -22 -21 -20 -19 -18 -17 -16 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 -25 0 25 50 -26 -24 -22 -20 -18 -16 -14 -12 -10 -8 -6 -4 -2 0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 42 44 46 48 50 Population

In [6]:
srand(22)

# 前期試験の得点とpika scoreの散布図作成
s_1 = reshape(pop_todai,1000,1)
s_2 = reshape(pop_w,1000,1)
score = vcat(s_1, s_2)
data = hcat(X, pika)
data = hcat(data, score)
n = 1000
cutoff = 347.211

# 回帰直線
d_todai = DataFrame(X = pop_todai, Y = todai)
result_todai = lm(Y ~ X, d_todai)
pre_todai = predict(result_todai)
d_w = DataFrame(X = pop_w, Y = wuni)
result_w = lm(Y ~ X, d_w)
pre_w = predict(result_w)

# plot
todai_layer = layer(x = pop_todai, y = pre_todai, Geom.line, Theme(default_color = colorant"black"))
w_layer = layer(x = pop_w, y = pre_w, Geom.line, Theme(default_color = colorant"black"))

gap = (predict(result_todai, [1, cutoff]') - predict(result_w, [1, cutoff]'))[1]
println("the effect of attend to todai on pika score is $gap")

plot(x = score, y = pika, todai_layer, w_layer, xintercept=[347.211], 
    Scale.x_continuous(minvalue= 200, maxvalue=500),
    Geom.vline(color = colorant"orange"),
    Geom.point,
    Guide.xlabel("Entrance Exam Score"), Guide.ylabel("Pika Score"), Guide.title("Scatter Plot"),
    Theme(default_point_size = 1.5pt)
)


the effect of attend to todai on pika score is 1.6033812858523626
Out[6]:
Entrance Exam Score -200 -100 0 100 200 300 400 500 600 700 800 900 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 -500 0 500 1000 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700 720 740 760 780 800 -20 -15 -10 -5 0 5 10 15 20 -15.0 -14.5 -14.0 -13.5 -13.0 -12.5 -12.0 -11.5 -11.0 -10.5 -10.0 -9.5 -9.0 -8.5 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0 10.5 11.0 11.5 12.0 12.5 13.0 13.5 14.0 14.5 15.0 -20 -10 0 10 20 -15 -14 -13 -12 -11 -10 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 Pika Score Scatter Plot

上の散布図をみてください。

オレンジの線が合格最低点である347.211点です。

オレンジの線を境にして、左側では低めに、右側では高めにpika scoreが出ていることが見て取れるでしょう。先ほど説明した状況になっているのがわかりますか?

ここでは東大に入学するという処置を受けられるか否かの境でpika scoreが不連続となっているのですね。

後の手順は簡単です。

オレンジ線の左側のデータのみで線形回帰を行いcut off(ここでは347.211点)における予測値を算出し、右側においても同様の手順でcut off における予測値を算出します。その時の予測値の差が求めている因果効果です。

ただ、一点注意しなくてはならない点があります。

それはデータをどこまで使うのか、という問題です。

知りたい因果効果が現出するのがオレンジ線上であることを考えると、前期試験の点数がとても高い学生や、極端に低い学生のデータは今必要ではありません。

私たちが知りたいのは、左側で作った予測値のオレンジ線における極限値と右側で作った予測値のオレンジ線における極限値の差です。

それゆえにオレンジ線を含む短い区間のデータのみに注目して切片を予測するのは理にかなったことといえます。

しかし、ある程度遠くのデータまで含めないと推定値が安定しないこと、またあまりにlocalな結果となってしまうことも事実です。

RDにおいては、このジレンマの中で最適なband width、すなわちデータを用いる幅を決定することが肝となってきます。

その際に理論的な背景を与えるのがLocal Linear Regressionです。以下ではこのLocal Linear Regressionについて解説し、上のデータを少し変えて最適なbandwidthを実際に求めてみたいとおもいます。

Local Linear Regression

上のデータは面白くないので、相関をつけます。さらに人数を実際に近づけて合格者400人、不合格者400人にします。 以下のモデルでpika scoreを算出し、上と同様のプロットを行います

$ pika_{i} = 10*{\rm log}(score_{i})+ \epsilon_{i} \\ pika_{i} = 9.8 * {\rm log}(score_{i})+ \epsilon_{i}$

In [7]:
srand(22)
# データ作成
srand(22)
n = 400
cutoff = 347.211
n_1 = n/2
n_2 = 3n/10
n_3 = 3n/20
n_4 = n/20
n_array = Int64[n_1, n_2, n_3, n_4]
pop_todai = Array(Float64)
for (i, r) in enumerate(n_array)
    pop_todai = [pop_todai ; [j*25 + cutoff + 25*(i-1) for j in rand(r)]]
end
pop_todai = Array{Float64}(pop_todai[2:end])
pop_w = Array(Float64)
for (i, r) in enumerate(n_array)
    pop_w = [pop_w ; [cutoff - j*25 - 25*(i-1) for j in rand(r)]]
end
pop_w = Array{Float64}(pop_w[2:end])
s_1 = reshape(pop_todai,n,1)
s_2 = reshape(pop_w,n,1)
score = vcat(s_1, s_2)
newpika = Array(Float64, 2n)
for i in 1:2*n
    if i <= n
        newpika[i] =   10*log(score[i]) + randn()
        else newpika[i] =  9.8*log(score[i])+ randn()
    end
end

# 回帰直線
d_todai = DataFrame(X = score[1:n], Y = newpika[1:n])
result_todai = lm(Y ~ X, d_todai)
pre_todai = predict(result_todai)
d_w = DataFrame(X = score[n + 1:end], Y = newpika[n + 1 : end])
result_w = lm(Y ~ X, d_w)
pre_w = predict(result_w)

gap = (predict(result_todai, [1, cutoff]') - predict(result_w, [1, cutoff]'))[1]
println("the effect of the attend to todai on pika score is $gap")

# plot
todai_layer = layer(x = pop_todai, y = pre_todai, Geom.line, Theme(default_color = colorant"black"))
w_layer = layer(x = pop_w, y = pre_w, Geom.line, Theme(default_color = colorant"black"))

plot(x = score, y = newpika, todai_layer, w_layer, xintercept=[347.211], 
    Scale.x_continuous(minvalue= 200, maxvalue=500),
    Geom.vline(color = colorant"orange"),
    Geom.point,
    Guide.xlabel("Entrance Exam Score"), Guide.ylabel("Pika Score"), Guide.title("Scatter Plot"),
    Theme(default_point_size = 1.5pt)
)


the effect of the attend to todai on pika score is 0.9455085258287923
Out[7]:
Entrance Exam Score -200 -100 0 100 200 300 400 500 600 700 800 900 -100 -90 -80 -70 -60 -50 -40 -30 -20 -10 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 300 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 500 510 520 530 540 550 560 570 580 590 600 610 620 630 640 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 800 -500 0 500 1000 -100 -80 -60 -40 -20 0 20 40 60 80 100 120 140 160 180 200 220 240 260 280 300 320 340 360 380 400 420 440 460 480 500 520 540 560 580 600 620 640 660 680 700 720 740 760 780 800 30 35 40 45 50 55 60 65 70 75 80 85 35.0 35.5 36.0 36.5 37.0 37.5 38.0 38.5 39.0 39.5 40.0 40.5 41.0 41.5 42.0 42.5 43.0 43.5 44.0 44.5 45.0 45.5 46.0 46.5 47.0 47.5 48.0 48.5 49.0 49.5 50.0 50.5 51.0 51.5 52.0 52.5 53.0 53.5 54.0 54.5 55.0 55.5 56.0 56.5 57.0 57.5 58.0 58.5 59.0 59.5 60.0 60.5 61.0 61.5 62.0 62.5 63.0 63.5 64.0 64.5 65.0 65.5 66.0 66.5 67.0 67.5 68.0 68.5 69.0 69.5 70.0 70.5 71.0 71.5 72.0 72.5 73.0 73.5 74.0 74.5 75.0 75.5 76.0 76.5 77.0 77.5 78.0 78.5 79.0 79.5 80.0 0 20 40 60 80 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 Pika Score Scatter Plot

local linear regressionは、その名の通り一部のデータを使って線形回帰を行う手法です。RDの文脈では、cut off周囲のデータのみを使って左側と右側で二つの線形回帰を行います。推定された係数をもとにcut offにおける推定値の差を算出し、それを持って処置の効果とする手法です。

その際にどこまでを周囲とするか、すなわちband widthを決めなくてはなりません。一般にこれは恣意的に決めるものではなく、data drivenな手法で求めるものされています。

band widthの決定には主に二つの手法が存在します。

1つ目はRule of Thumbで初期のband widthを決定し、それを改善していくもの。

2つ目はcross validationを行って、平均二乗誤差を最小にするband widthを選択するものです。

今回は2つ目の手法を紹介します。

cross validationとは一般に以下のような手法を指します。

  1. 既存のデータをテスト用と学習用に分割する
  2. 学習用データを用いてパラメータを推定する
  3. 推定したパラメータをもとに推定値を算出し、テスト用のデータでその推定精度を確認する

Local Linear Regressionの文脈では、これは以下の用に用いられます。

  1. cut offを境に左側に適当なband width (h)をとる
  2. band widthに含まれるデータをN個として、1つのデータを除いたN-1個のデータで回帰分析を行う
  3. 上で得られた回帰直線をもとに、除いたデータの推定値を算出し、実データとの差を保存する
  4. 2,3の作業をband widthに含まれるすべての点に対して行い、3で得た差の二乗の平均を$CV(h)$とする
  5. 複数のband widthに対して上記の作業を行い、$CV(h)$を最小にするhをband widthとして選択する

では、以下で実際にband widthを選択してみましょう。


In [65]:
# データ作成
srand(22)
n = 400
cutoff = 347.211
n_1 = n/2
n_2 = 3n/10
n_3 = 3n/20
n_4 = n/20
n_array = Int64[n_1, n_2, n_3, n_4]
pop_todai = Array(Float64)
for (i, r) in enumerate(n_array)
    pop_todai = [pop_todai ; [j*25 + cutoff + 25*(i-1) for j in rand(r)]]
end
pop_todai = Array{Float64}(pop_todai[2:end])
pop_w = Array(Float64)
for (i, r) in enumerate(n_array)
    pop_w = [pop_w ; [cutoff - j*25 - 25*(i-1) for j in rand(r)]]
end
pop_w = Array{Float64}(pop_w[2:end])
s_1 = reshape(pop_todai,n,1)
s_2 = reshape(pop_w,n,1)
score = vcat(s_1, s_2)
newpika = Array(Float64, 2n)
for i in 1:2*n
    if i <= n
        newpika[i] =   10*log(score[i]) + randn()
        else newpika[i] =  9.8*log(score[i])+ randn()
    end
end

# cross validation
band_widthes = collect(1.0:0.5:15.0)
validation_l = Array{Float64}(length(band_widthes))
validation_r = Array{Float64}(length(band_widthes))
data = DataFrame(Y = newpika, X = squeeze(score',tuple(1)))

for (t, i) in enumerate(band_widthes)
    gaps = Float64[]
    temp_d = data[cutoff .> data[2] .> cutoff-i, :]
    temp_f = data[cutoff .> data[2] .> cutoff-i, :]
    
    for j in 1:size(temp_d)[1]
        temp_d[j, 1] = NA
        independent = temp_f[j, 2]
        dependent = temp_f[j, 1]
        temp_res = lm(Y ~ X, temp_d)
        predictval = predict(temp_res, [1, independent]')
        squared_gap = ((dependent - predictval)[1])^2
        push!(gaps, squared_gap)
        temp_d[j, 1] = dependent
    end
    
    validation_l[t] = mean(gaps)
end

for (t, i) in enumerate(band_widthes)
    gaps = Float64[]
    temp_d = data[cutoff+i .> data[2] .> cutoff, :]
    temp_f = data[cutoff+i .> data[2] .> cutoff, :]
    
    for j in 1:size(temp_d)[1]
        temp_d[j, 1] = NA
        independent = temp_f[j, 2]
        dependent = temp_f[j, 1]
        temp_res = lm(Y ~ X, temp_d)
        predictval = predict(temp_res, [1, independent]')
        squared_gap = ((dependent - predictval)[1])^2
        push!(gaps, squared_gap)
        temp_d[j, 1] = dependent
    end
    validation_r[t] = mean(gaps)
end

leftband = band_widthes[indmin(validation_l)]
rightband = band_widthes[indmin(validation_r)]

for (i, n) in enumerate(band_widthes)
    println("mean squared error for left when band = $n", " is ", validation_l[i])
    println("mean squared error for right when band = $n", " is ", validation_r[i])
end

println("selected band width for the left side is $leftband")
println("selected band width for the right side is $rightband")


mean squared error for left when band = 1.0 is 1.4150357114945884
mean squared error for right when band = 1.0 is 2.2034753660874835
mean squared error for left when band = 1.5 is 0.7831889452116898
mean squared error for right when band = 1.5 is 1.2442663673074825
mean squared error for left when band = 2.0 is 0.8725352613103292
mean squared error for right when band = 2.0 is 1.3159373512332173
mean squared error for left when band = 2.5 is 0.7831064630410849
mean squared error for right when band = 2.5 is 1.3338175166413528
mean squared error for left when band = 3.0 is 0.811762326527483
mean squared error for right when band = 3.0 is 1.110549097233108
mean squared error for left when band = 3.5 is 1.0485440330361013
mean squared error for right when band = 3.5 is 1.0498761812413278
mean squared error for left when band = 4.0 is 1.8253716456581717
mean squared error for right when band = 4.0 is 1.1436881965911814
mean squared error for left when band = 4.5 is 1.6281349443489475
mean squared error for right when band = 4.5 is 1.1116556827511102
mean squared error for left when band = 5.0 is 1.5800501351360665
mean squared error for right when band = 5.0 is 1.010266666893191
mean squared error for left when band = 5.5 is 1.391303858899397
mean squared error for right when band = 5.5 is 1.0803310236296308
mean squared error for left when band = 6.0 is 1.37980909662419
mean squared error for right when band = 6.0 is 1.0822190598652544
mean squared error for left when band = 6.5 is 1.336001780473925
mean squared error for right when band = 6.5 is 1.1041952011412464
mean squared error for left when band = 7.0 is 1.3710574628701317
mean squared error for right when band = 7.0 is 1.0931307111804276
mean squared error for left when band = 7.5 is 1.2796711109169714
mean squared error for right when band = 7.5 is 1.062109391910558
mean squared error for left when band = 8.0 is 1.2550018312693498
mean squared error for right when band = 8.0 is 1.127928312517777
mean squared error for left when band = 8.5 is 1.2550018312693498
mean squared error for right when band = 8.5 is 1.10986615238273
mean squared error for left when band = 9.0 is 1.2992994145089154
mean squared error for right when band = 9.0 is 1.139600269747505
mean squared error for left when band = 9.5 is 1.2896963112484208
mean squared error for right when band = 9.5 is 1.174524770455114
mean squared error for left when band = 10.0 is 1.26417799309661
mean squared error for right when band = 10.0 is 1.1796090747004893
mean squared error for left when band = 10.5 is 1.228425980978268
mean squared error for right when band = 10.5 is 1.1596671435285775
mean squared error for left when band = 11.0 is 1.2138857326041939
mean squared error for right when band = 11.0 is 1.109857373241189
mean squared error for left when band = 11.5 is 1.2212790759529617
mean squared error for right when band = 11.5 is 1.1627882274184382
mean squared error for left when band = 12.0 is 1.2134540479283087
mean squared error for right when band = 12.0 is 1.1286752691520257
mean squared error for left when band = 12.5 is 1.1915340027555437
mean squared error for right when band = 12.5 is 1.1122991862215226
mean squared error for left when band = 13.0 is 1.1842675997981358
mean squared error for right when band = 13.0 is 1.082143970688779
mean squared error for left when band = 13.5 is 1.2035015000354474
mean squared error for right when band = 13.5 is 1.0944982751767975
mean squared error for left when band = 14.0 is 1.1767460148056839
mean squared error for right when band = 14.0 is 1.0795560257663233
mean squared error for left when band = 14.5 is 1.1521144761696713
mean squared error for right when band = 14.5 is 1.0643953235881884
mean squared error for left when band = 15.0 is 1.2176377976945791
mean squared error for right when band = 15.0 is 1.0468007096241168
selected band width for the left side is 2.5
selected band width for the right side is 5.0

上の例では cutoffの左側ではbandwidth = 2.5が選択され、右側ではbandwidth = 5.0が選択されました。ちなみにNaNは、その範囲にデータが存在しなかったことを意味しています。

では、上記の最適なbandwidthに従ってデータを以下にプロットし、先と同様に東大進学の効果を計算してみましょう。


In [63]:
srand(22)
# 回帰直線
cutoff = 347.211
n = 400
d_todai = DataFrame(X = score[1:n], Y = newpika[1:n])
d_todai = d_todai[cutoff + rightband .>d_todai[1], :]
result_todai = lm(Y ~ X, d_todai)
pre_todai = predict(result_todai)
d_w = DataFrame(X = score[n + 1:end], Y = newpika[n + 1 : end])
d_w = d_w[d_w[1] .> cutoff-leftband , :]
result_w = lm(Y ~ X, d_w)
pre_w = predict(result_w)

data = DataFrame(Y = newpika, X = squeeze(score',tuple(1)))

gap = (predict(result_todai, [1, cutoff]') - predict(result_w, [1, cutoff]'))[1]
println("the effect of the attend to todai on pika score is $gap")

# plot
todai_layer = layer(x = d_todai[1], y = pre_todai, Geom.line, Theme(default_color = colorant"black"))
w_layer = layer(x = d_w[1], y = pre_w, Geom.line, Theme(default_color = colorant"black"))

plot(x = data[cutoff+rightband .>data[2].> cutoff-leftband, 2], y = data[cutoff+rightband .>data[2].> cutoff-leftband, 1], todai_layer, w_layer, xintercept=[347.211], 
    Scale.x_continuous(minvalue= 340, maxvalue=350),
    Geom.vline(color = colorant"orange"),
    Geom.point,
    Guide.xlabel("Entrance Exam Score"), Guide.ylabel("Pika Score"), Guide.title("Scatter Plot"),
    Theme(default_point_size = 1.5pt)
)


the effect of the attend to todai on pika score is 1.3778080765596599
Out[63]:
Entrance Exam Score 320 325 330 335 340 345 350 355 360 365 370 375 325.0 325.5 326.0 326.5 327.0 327.5 328.0 328.5 329.0 329.5 330.0 330.5 331.0 331.5 332.0 332.5 333.0 333.5 334.0 334.5