Chapter 16 回帰不連続デザイン

16.1 回帰不連続デザインとは

回帰不連続デザイン(Regression Discontinuity Design, RDD)も、回帰不連続デザインと同様疑似実験を用いた手法の一つである。 連続的に変化すると考えられるある変数が、政策介入や処置によって急激な変化(ジャンプ)する場合に、その急激な変化を介入や処置による効果として推定する手法である。

16.1.1 回帰不連続デザインのアイデア

お酒・アルコールを飲むのは(好きな人には)楽しいですが、当然負の側面もある。健康被害を起こしたり、また事故なども起こる可能性がある。

アルコールを接種することで、死亡率は上がるのだろうか。

アルコールを接種するかどうかは、法律によって年齢が決まっている。「お酒は二十歳になってから」と言われるように、日本では飲酒の年齢は20歳からである。

つまり、18歳や19歳ならばお酒は飲まないが、20歳以上ならば飲酒するという傾向があるのである。(当然ならば法律を守っているかどうか、という議論はあるが、すくなくても法律を守っていないのは一部で多くは20歳から飲酒するだろう)

しかし、18〜22歳ぐらいまでの若者は基本的には似た性質を持っていると考えられる。20歳になって、いきなり死亡率が増えるとは考えにくい。そのため、18歳や19歳に比べて、20歳や21歳で突然死亡率が上がっているならば、それは20歳をカットオフ(閾値)として扱い(処置)が変わる飲酒傾向によるものと考えることができるだろう。

18-19歳の人々は「対照群」、20-22歳の人を「処置群」として考えることができ、比較が可能となる。これが、基本的な回帰不連続デザインのアイデアである。

年齢によって死亡率は異なるが、少なくともその変化は滑らかであると考えられる。そこで、境界線上に大きな変化(ジャンプ)が起こっていれば、それが政策や処置の効果であると考えて推定を行うのである。

ここでは、年齢のような変数を割当変数と呼ぶ。

16.1.2 回帰不連続デザインの仮定

  • 割当変数と目的変数の間には政策がなければ滑らかに変化する関係がある
  • カットオフ前後では政策以外に大きな相違がなくほぼ同じである
  • カットオフの存在によって、割当変数を操作するようなことがない
    • 例:奨学金がもらえるカットオフが80点だとする。そこがカットオフだと知っていれば、80点にギリギリみたないような人がいつも以上に頑張って80点以上を取ろうとする。これは、政策(奨学金給付)があることで、政策がなければないようなカットオフ間の移動が起こっているため、仮定を満たさない。

16.1.3 回帰分析を使って行う方法

回帰不連続デザインを回帰分析で行う場合は、式としていくつかの方法があるが、ここでは比較グループと政策実施グループで異なる傾きを持つ分割線形回帰を紹介する。

まず回帰式は以下のようになる。

Yi=α+ρPolicyi+β1xi+β2Policyi×xi+εi ここで、Policyiはダミー変数である。Policyiはもし、該当データが、割当変数でカットオフより大きい(政策によっては小さい)個人iが、処置群に入るならば1、対照群ならば0を取る。

例えば、飲酒年齢の問題では年齢が割当変数(xi)となる。そして、20歳以上ならば、Policyiは1となるが、20歳未満はPolicyiが0となる。

ここでは、因果効果として推定されるパラメータはρ(ローと読む)である。

上の概念図に対応させると、以下のようになる。

左右の直線は、政策が適用される側(処置群)とされない側(対照群)でデータに基づいてそれぞれ推定された回帰直線である。

もし仮に、推定された結果についてβ2ρがゼロだとする。すると、回帰式は

Yi=α+β1xi+εi となり、単なる単回帰式となる。図で表すと以下のような形だ。

この場合は、当然ながら政策によって割当変数と目的変数のスムーズな関係に変化は起こっていない。

もしβ2がマイナスの値で、ρがゼロだったらどうなるか。 この場合は以下のような式になる。

Yi=α+β1xi+β2Policyi×xi+εi

この場合のxiの係数は、β1xi+β2Policyi×xi=(βi+β2Policyi)xiと解釈することができる。つまり、カットオフより下ならばPolicyi=0なので、係数はβ1になるが、カットオフより上ならばPolicyi=1なので、係数はβ1+β2となる。

もし政策効果がないならば、左右の直線はつながった形になる。もし政策がなければ、目的変数は割当変数に従って滑らかに変化する。もしかしたら比例関係ではなく、どこかで傾向が変わる関係かもしれない。そのためにβ2を推定して、傾向が変わることを許すモデルにしているのである。仮にそのような関係があったとしても、「ジャンプ」がない限りは政策効果はないと考えられる。

政策効果ρがゼロでない場合のみ、2つの直線の間に「ジャンプ」が生じる。これが政策が行われたことによる効果であると解釈される。

16.2 Rでの演習

16.2.1 スクリプトの準備

新しいスクリプトにコードを書いていく。 新しくスクリプトを作成し、rdd.Rという名前をつけて保存する。

16.2.2 パッケージの準備

以下のパッケージを使うので読み込んでおく。 インストールしてない場合は、インストールする。

# ライブラリ
library(tidyverse)
library(skimr)
library(fixest)

16.2.3 データの準備

データをダウンロードする。

# データのダウンロードとdataフォルダへの格納を自動的に行う
download.file("https://github.com/keita43a/regression_tutorial/blob/main/data/data_rdd.csv?raw=TRUE", 
              destfile = "data/data_rdd.csv")

# データを読み込む
data_rdd = read_csv("data/data_rdd.csv")

16.2.4 データの確認

このデータは生成されたデータであるが、実際の研究に使われたデータを模倣したものである。 Shigeoka (2014)では、日本の医療費の自己負担割合が70歳を境に3割から1割に下がることに着目し、この自己負担割合の変化が高齢者の死亡率を下げることに貢献しているかどうかを検証した。もし、自己負担割合が下がることで高齢者が病気や怪我のときに病院に行くようになり、それが結果として死亡率の低下に貢献しているのならばこの政策は一定の効果があると言える。

回帰不連続デザインを使う理由は70歳というカットオフが存在するからである。69歳11ヶ月の人は自己負担割合が3割だが、70歳0ヶ月になると突然1割になる。しかし、70歳になったからといって突然不健康になったり死亡率が上がるわけではないため、70歳になって急激に死亡率が下がっているのならば、これは自己負担割合の違いが死亡率の低下に貢献していると言えるというわけだ。

データの中身を見てみる。 まずは、どんなデータなのかhead()で確認する。

head(data_rdd)
## # A tibble: 6 × 3
##     age   visits  death
##   <dbl>    <dbl>  <dbl>
## 1  65.1 -0.0922  -0.281
## 2  65.2 -0.0794  -0.298
## 3  65.2 -0.115   -0.286
## 4  65.3 -0.0603  -0.282
## 5  65.3 -0.121   -0.286
## 6  65.5 -0.00707 -0.281

16.2.5 データの記述統計

データの記述統計を確認しよう。 Chapter @ref(#rstats) で説明したように、ここではskimrパッケージのskim()関数を使った。

library(skimr)

skim(data_rdd)
Table 16.1: Data summary
Name data_rdd
Number of rows 100
Number of columns 3
_______________________
Column type frequency:
numeric 3
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
age 0 1 69.85 3.02 65.10 67.19 69.94 72.50 74.93 ▇▆▇▅▇
visits 0 1 -0.03 0.08 -0.18 -0.10 -0.03 0.03 0.13 ▃▇▆▅▃
death 0 1 -0.01 0.18 -0.30 -0.17 -0.01 0.15 0.30 ▇▆▆▆▇

16.2.6 データの加工

データには年齢は記録されているが、70歳以上かどうかを示す情報は直接は記録されていない。 年齢から、70歳以上かどうか、というダミー変数を作成しよう。 ifelse(test,true,false)という関数は1つ目の引数であるtestに条件式を入れる。その条件式を評価して、真であれば2つ目の引数であるtrueに入力された値を返す。偽であればfalseである3つ目の引数の値を返す。この場合はage >= 70が条件式なので、各行においてageが70以上であれば、1を入れ、そうでなければ0を入れる。

data_rdd <- data_rdd |>
  mutate(dummy70 = ifelse(age >= 70, 1, 0))

16.2.7 政策による行動への影響

ここでは、政策によって実際に行動が変わっているかをまず確認する。

下のコードでは、aes()によるマッピングの中に、groupという項目が入っている。これは、描画に使うデータを年齢が70歳以上か未満かdummy70で分ける機能を持つ。 こうすることで、geom_smoothを使うときに、70歳未満と以上で別々の直線を描いてくれる。 geom_vlineは縦(vertical)の直線を追加するレイヤーである。この場合は、x軸上の70歳(xintercept = 70)に縦の直線を、線の種類をダッシュ(linetype="dashed")で、色を赤で(color="red")追加する。

ggplot(data_rdd, aes(x = age, y = visits, group = dummy70)) +
  geom_point(color = 'gray') +
  geom_smooth(method = 'lm', color = 'black', se = FALSE) +
  geom_vline(xintercept = 70, linetype = "dashed", color = "red") +
  labs(x = "年齢", y = "外来患者数(対数)", title = "自己負担割合が減ると外来受診が増える") +
  theme_minimal(base_family="HiraKakuPro-W3")

これを見ると、外来患者数には70歳を境に大きなジャンプがある。 つまり、自己負担割合が低くなることで、外来患者の数が急激に増加している。 このことから、自己負担割合の変化は外来診察を受けるインセンティブを高め、高齢患者の行動を変えていると言えそうだ。

16.2.7.1 回帰分析による推定

回帰分析でこの効果を推定してみよう。

この場合は

Visitsi=α+ρDummy70i+β1agei+β2agei×Dummy70i+εi という式を推定することになる。

rdd_visit <- feols(visits ~ dummy70 + age + dummy70:age, data=data_rdd)

etable(rdd_visit)
##                              rdd_visit
## Dependent Var.:                 visits
##                                       
## Constant             -0.5376. (0.2933)
## dummy70              1.789*** (0.4178)
## 店長の年齢(歳)     0.0066 (0.0044)
## dummy70 x 店長の年齢(歳) -0.0233*** (0.0060)
## __________________ ___________________
## S.E. type                          IID
## Observations                       100
## R2                             0.72216
## Adj. R2                        0.71348
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

結果はρが1.789と推定された。Visitsは対数で表示されているため、このジャンプによって訪問患者数が179%増加したという結果になる。

繰り返しになるが、これは元の研究を模倣したデータであり、説明のために大きく省略したアプローチを取っているのでオリジナルの研究とこの数字は一致しない。

16.2.8 政策による結果への影響

今度は、政策によって重要な結果、ここでは死亡率に効果があるかを見てみよう。

ggplot(data_rdd, aes(x = age, y = death, group = dummy70)) +
  geom_point(color = 'gray') +
  geom_smooth(method = 'lm', color = 'black', se = FALSE) +
  geom_vline(xintercept = 70, linetype = "dashed", color = "red") +
  labs(x = "年齢", y = "死亡率(対数)", title = "自己負担割合が減っても死亡率は変わらない") +
  theme_minimal(base_family="HiraKakuPro-W3")

上と同様に分割して回帰直線を描いているが、あたかもつながっているかのような直線になっている。大きなジャンプは見られず、70歳になって突然死亡率が下がったという傾向はなさそうである。

16.2.8.1 回帰分析による推定

回帰分析でこの効果を推定してみよう。

この場合は

Deathi=α+ρDummy70i+β1agei+β2agei×Dummy70i+εi という式を推定することになる。

rdd_death <- feols(death ~ dummy70 + age + dummy70:age, data=data_rdd)

etable(rdd_death)
##                             rdd_death
## Dependent Var.:                 death
##                                      
## Constant           -4.283*** (0.0680)
## dummy70               0.1488 (0.0969)
## 店長の年齢(歳) 0.0612*** (0.0010)
## dummy70 x 店長の年齢(歳)   -0.0021 (0.0014)
## __________________ __________________
## S.E. type                         IID
## Observations                      100
## R2                            0.99693
## Adj. R2                       0.99683
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

結果はρが0.1488と推定された。また、統計的に有意ではないため、0であるという仮説を棄却できない。 このことから、死亡率に対する健康保険の自己負担割合の低下は効果を与えていないという結論が導かれる。

16.3 演習問題:アメリカの飲酒年齢と死亡率

  1. 以下のコードを実行してデータをダウンロードして読み込む。
# データのダウンロードとdataフォルダへの格納を自動的に行う
download.file("https://github.com/keita43a/regression_tutorial/blob/main/data/data_alcohol.csv?raw=TRUE", 
              destfile = "data/data_alcohol.csv")

# データを読み込む
data_alcohol= read_csv("data/data_alcohol.csv")

これはアメリカの誕生日によって区切られた年齢グループ別のデータである。約一ヶ月ごとに誕生日からどれぐらい経っているか(例:20歳の誕生日から1ヶ月以上2ヶ月未満)などのグループ(agecell)で区切られている。またそのグループ別に、死因別や合計(all_death)の死亡数が記録されている。

  1. データの記述統計を確認する。

  2. agecellを割当変数、policyを政策変数、all_deathを目的変数として回帰不連続デザインの図を描きなさい。アメリカの飲酒年齢は21歳である。

  1. 以下のモデルを推定して、結果をコンソールに表示する。ageagecellを21歳を中心として基準化した変数である。(例えば19歳はageでは-2となる。)

all_deathit=α+ρPolicyi+β1aget+β2Policyi×aget+εit

  1. 結果を解釈する。