ill-identified diary

所属組織の見解などとは一切関係なく小難しい話しかしません

トップ > 計量経済学 > [計量経済学][時事ネタ?] 非連続回帰デザイン (RDD) 実践編

[計量経済学][時事ネタ?] 非連続回帰デザイン (RDD) 実践編

計量経済学 時事ネタ 論文解説 因果推論
この記事は最終更新日から3年以上が経過しています

前回のRDD理論編を踏まえて, 前回 紹介した, Angrist and lavy(1999) "Using Maimonides' Rule to Estimate the Effect of Class Size on Scholastic Achievement" (以下 AL1999 と呼ぶ) のサンプルデータが著者の一人アングリストの所属する MIT Economics : Angrist Data Archive からダウンロードできるので, これを例にやってみよう.

ただし, RDDの例としては, AL1999 は特殊である. というのも, 前回紹介した, 2000年以降の研究では, 説明変数の1つである running variable が一定の値を超えたことに対する因果効果の大きさを求めるものだった. つまり, ゼロか1かの2値変数の係数を推定することを前提としていた. 一方 AL1999 では, クラスの人数という連続変数 (整数しかありえないので厳密には連続変数ではないが) の係数を求めるものだった. rdd パッケージでも, 2値変数の係数の推定を前提としているため, AL1999 以外にも, Lee, David S., Moretti, Enrico, and Butler, Matthew J. (2004) "Do Voters Affect or Elect Policies? Evidence from the U. S. House" (以下 LMB2004 と呼ぶ) のデータを使って実演する.

AL1999 の再現

2019/7/12 追記: 山口先生のツイートで知ったのだが, 最近の研究では AL1999の著者ら自らの追試により, AL1999 で見られたような少人数教育の効果はない, と主張している (Angrist et al., 2017).

AL1999 のダウンロードできるデータは STATA 用の .dta データ*1で, 加工用のプログラムも STATA のものなので, ここで自分が R で書きなおしたものを使う*2. データセットの読み込みには, STATA の .dta 形式のファイルを読み込むには, foreign パッケージ (現時点で ver. 0.8-61) が必要である. まず, データを加工する必要がある. 以下はデータセットと同じところでダウンロードできる STATA の .do コード*3をそのまま R用に書き直したもので, 特に複雑なこともしてないので解説は省く.

library(foreign) # foreign 0.8-61

# http://economics.mit.edu/faculty/angrist/data1/data/anglavy99
df.5 <- read.dta(file="source/final5.dta", warn.missing.labels = T) # 第5学年データを読み込み. パスは適宜変更

# ヘブライ文字が化けていたので変換 (使わないけど)
df.5$townname <- iconv(df.5$townname, from="CP862", to="UTF-8")
# 変数の説明を取得できる
(labels <- attributes(df.5)$var.labels)
colnames(df.5)

# AngristLavy_Table4.do を R用に書き換え
# 数学・言語の得点で100超のものを補正
df.5$avgmath <- ifelse(df.5$avgmath>100, df.5$avgmath-100, df.5$avgmath) 
df.5$avgverb <- ifelse(df.5$avgverb>100, df.5$avgverb-100, df.5$avgverb) 
# クラスサイズ関数による予測値計算
df.5$func1   <- df.5$c_size/(as.integer( (df.5$c_size-1)/40)+1) 
df.5$func2   <- df.5$cohsize/(as.integer(df.5$cohsize/40)+1)
# 欠損値に変換
df.5$avgverb[df.5$verbsize==0]  <- NA
df.5$passverb[df.5$varbsize==0] <- NA
df.5$avgmath[df.5$mathsize==0]  <- NA
df.5$passmath[df.5$mathsize==0] <- NA

# 一部の変な個体を落とす
df.5.2 <- df.5[1 < df.5$classize & df.5$classize < 45 & df.5$c_size >5,] # 5件脱落
df.5.2 <- df.5.2[df.5.2$c_leom ==1 & df.5.2$c_pik <3,] # 脱落なし
df.5.2 <- df.5.2[!is.na(df.5.2$avgverb),] # 5件脱落

# Discontinuity サンプルの判定用変数
df.5.2$disc <- ifelse( (df.5.2$c_size>=36 & df.5.2$c_size<=45) | (df.5.2$c_size>=76 & df.5.2$c_size<=85) | 
                         (df.5.2$c_size>=116 & df.5.2$c_size<=125),T,F)

# クラスサイズ2乗項
df.5.2$c_size2 <- df.5.2$c_size^2/100

# 線形区分関数の作成
df.5.2$trend <- ifelse(0  <= df.5.2$c_size & df.5.2$c_size <=40,  df.5.2$c_size, 0)
df.5.2$trend <- ifelse(41 <= df.5.2$c_size & df.5.2$c_size <=80,  df.5.2$c_size/2 + 20, df.5.2$trend)
df.5.2$trend <- ifelse(81 <= df.5.2$c_size & df.5.2$c_size <=120, df.5.2$c_size/3 + 100/3, df.5.2$trend)
df.5.2$trend <- ifelse(121<= df.5.2$c_size & df.5.2$c_size <=160, df.5.2$c_size/4 + 130/4, df.5.2$trend)

### データセット作成ここまで ###

多少レイアウトは違うが, FigureI, II と同等のグラフを掲載する. このグラフ作成には, ggplot2 (ver. 1.0.0), plyr (ver. 1.8.1), reshape2 (ver. 1.4) パッケージを利用している.

library(ggplot2) # ver. 1.0.0.
library(plyr) # ver. 1.8.1 諸事情により dplyr が使えない状況だったので
library(reshape2) # ver. 1.4

# Fig. I, 入学者数とクラスサイズ理論値・実際のクラスサイズ 
# クラスサイズ関数描画用のデータセット
classize_line <- data.frame(x=1:max(df.5.2$c_size))
classize_line$y <- classize_line$x/(as.integer( (classize_line$x -1)/40)+1)
ggplot(data=df.5.2) + geom_line(mapping = aes(x, y), data=classize_line, linetype=4) + geom_point(aes(c_size,classize), color="blue") + labs(x="入学者数", y="平均クラスサイズ", title="入学者数とクラスサイズ") +theme_bw() + scale_x_discrete()

f:id:ill-identified:20141218011134p:plain

元の Figure I では平均だったが, ここでは全ての点をプロットした. 一点破線はクラスサイズ関数である. クラスサイズ関数にすべて一致するわけではないが, クラスサイズ関数のまわりに点が集中していることが分かる. よって, カットオフポイントが "Fuzzy" である.

# Fig. II, 成績と入学者数のプロット
# 入学者数10人間隔で集計
avg.score <- df.5.2[,c("c_size", "avgmath","avgverb")]; avg.score$c_size <- ceiling(avg.score$c_size/10)*10
avg.score <- ddply(.data = avg.score, .(c_size), summarise, mean.verb=mean(avgverb,na.rm=T), mean.math=mean(avgmath,na.rm=T),  .progress = "text")
avg.score <- melt(data = avg.score, id="c_size",value.name = "score", measure.vars =c("mean.verb", "mean.math") )
classize_line <- data.frame(x=1:22*(10))
classize_line$y <- classize_line$x/(as.integer( (classize_line$x -1)/40)+1)

ggplot(data=avg.score) + geom_line(mapping = aes(c_size, score, group=variable, linetype=variable)) + labs(x="入学者数",y="平均得点", title="入学者数ごとの平均得点") + scale_linetype_discrete(name="", labels=c("読解力", "算数") ) + geom_line(mapping = aes(x, y+40), data=classize_line, linetype=4) + xlim(5, 170) + ylim(50,85)

f:id:ill-identified:20141218011207p:plain

次に, 入学者数と成績の関係を見る. 点数は, 入学者数10人ごとに平均を取った値である. クラスサイズ関数も併記したが, 当然ながら y軸のスケールが違うので注意. アングリストらの研究で指摘されているように, クラスサイズのカットオフポイント (クラスサイズ関数が山を描いている点) で, 点数が谷になっていることから, 点数の推移にも非連続な点が存在することが示唆される.

次に, 推定結果の再現だが, 今回は第5学年について推定した表 IV の再現をする. 推定には, 実質的には2段階最小二乗法を使えばいいが, rdd というパッケージは aer パッケージに依存し, 名前通りRDD に使える関数を用意している. rdd 固有のメソッドは AL1999 では活用できないが, 後で紹介する LMB2004 では必要である. Angrist and Lavy (1999) においては, 単なる2段階最小二乗法として推定をしているため, ivreg() を使えば良い.

例として, 読解力テストの成績を被説明変数にした時の RDD による推定は, 下のようになる.

library(sdd) # ver. 0.5.6
summary(res <- ivreg(formula= avgverb~classize + tipuach| func1+ tipuach , data=df.5.2, subset=df.5.2$disc) )

avgverb は読解力の平均点, classize は実際のクラスサイズ, tipuach は PD を表す. ivreg() 関数に与える forumla は, (被説明変数)~(説明変数1)+(説明変数2)+... | (操作変数1) + (操作変数2) + ... という構文になっているので, この場合は説明変数がクラスサイズとPD, | で挟んだ誘導形方程式には, 内生変数でない説明変数と操作変数を並べる. また, disc はカットオフからプラスマイナス5点の範囲内にある変数かどうかを判定するダミー変数である*4. これで元論文の表IVの列 (5) と同じ係数が出るはずである. ただし, 標準誤差については一致しない. なぜなら表IVは通常の標準誤差ではなくクラスタリング標準誤差を掲載しているからだ. これはクラス単位のデータであるため, 同一の学校から抽出されたクラスのレコードがいくつもある. よって, 説明変数からは観察できないその学校固有の特徴がある場合, 誤差項の相関が発生する. サンプルこのような個体の「クラスタ」が存在する場合, 通常の計算方法では標準誤差が過少に算出されてしまう. よって, ivpack パッケージを用いて, クラスタリング標準誤差を別途計算する. 標準誤差の詳しい話は, Cameron, Gelbach and Miller (2006) "Robust inference with multi-way clustering" や 『ほとんど無害な計量経済学』の8章などを参考に.

library(ivpack) # ver. 1.2
res.se <- cluster.robust.se(res, clusterid = df.5.2$schlcode[df.5.2$disc & !is.na(df.5.2[,"avgverb"])])[,2]
res.se

ivreg の出力と, サンプルサイズと同じ要素数の, クラスタを分類するベクトルを与えると標準誤差を計算してくれる. 今回はデータセットに学校の識別コード schlcode があるので, それを利用している.

# 表IVの推定をまとめて実行
depvars <- c("avgverb","avgmath")
result <- list()
for (dep in depvars ){
  # 1 
  result[[dep]][[1]] <- as.list(ivreg(formula= as.formula(paste(dep,"classize + tipuach| func1 + tipuach", sep="~" ) ), data=df.5.2))
  # 2
  result[[dep]][[2]] <- as.list(ivreg(formula= as.formula(paste(dep,"classize + tipuach + c_size| func1 + tipuach + c_size", sep="~" ) ), data=df.5.2))
  # 3
  result[[dep]][[3]] <- as.list(ivreg(formula= as.formula(paste(dep,"classize + tipuach + c_size + c_size2| func1 + tipuach + c_size + c_size2", sep="~" ) ), data=df.5.2))
  # 4 trend
  result[[dep]][[4]] <- as.list(ivreg(formula= as.formula(paste(dep,"classize + tipuach + trend| func1 + tipuach + trend", sep="~" ) ), data=df.5.2))
  # RDD
  result[[dep]][[5]] <- as.list(ivreg(formula= as.formula(paste(dep,"classize + tipuach| func1+ tipuach", sep="~" ) ), data=df.5.2, subset=df.5.2$disc))
  result[[dep]][[6]] <- as.list(ivreg(formula= as.formula(paste(dep,"classize + tipuach + c_size| func1 + tipuach + c_size", sep="~" ) ), data=df.5.2, subset=df.5.2$disc))
}

# 元論文では clustering-robust SE を使っているので修正
l.clu.se <- list()
for (dep in depvars){
    # (1)~(4) は通常の2SLS
  for (i in 1:4){
    l.clu.se[[dep]][[i]] <- cluster.robust.se(result[[dep]][[i]], clusterid = df.5.2$schlcode[!is.na(df.5.2[,dep])])[,2]
  }
    #RDD はサンプルを選別しているので注意
  for (i in 5:6){
    l.clu.se[[dep]][[i]] <- cluster.robust.se(result[[dep]][[i]], clusterid = df.5.2$schlcode[df.5.2$disc & !is.na(df.5.2[,dep])])[,2]
  }
}

これらを利用して, 表IVと同等の結果を得られた.

library(stargazer) # ver. 5.1. 表に出力する
# tab 4 の texコード
res <- stargazer(list(result[[1]][[1]], 
                      result[[1]][[2]],
                      result[[1]][[3]],
                      result[[1]][[4]],
                      result[[1]][[5]],
                      result[[1]][[6]],
                      result[[2]][[1]], 
                      result[[2]][[2]],
                      result[[2]][[3]],
                      result[[2]][[4]],
                      result[[2]][[5]],
                      result[[2]][[6]]),
          se=list(l.clu.se[[1]][[1]],
                  l.clu.se[[1]][[2]],
                  l.clu.se[[1]][[3]],
                  l.clu.se[[1]][[4]],
                  l.clu.se[[1]][[5]],
                  l.clu.se[[1]][[6]],
                  l.clu.se[[2]][[1]],
                  l.clu.se[[2]][[2]],
                  l.clu.se[[2]][[3]],
                  l.clu.se[[2]][[4]],
                  l.clu.se[[2]][[5]],
                  l.clu.se[[2]][[6]]),
          align=T, no.space=T, type="latex", dep.var.labels.include=F,dep.var.labels=c("verb","math"),
          column.separate =c(4,2,4,2), column.labels=c("full sample","+/- 5 sample", "full sample", "+/- 5 sample"),
          covariate.labels=c("定数項","クラスサイズ", "PD", "入学者数", "入学者数$^2/100$", "線形区分"),
          keep.stat =c("n","rsq")
)

write(res, "tab4.tex")

表の作成には stargazer を使用している. stargazer の使用法は以前 [R] [LaTeX] R での分析結果を LaTeX 形式で出力するパッケージ比較 (後編) でも紹介したので参考にしてほしい.

f:id:ill-identified:20141218011453p:plain

実例その2

去る12月14日に衆院選挙が行われたが, 第2の例として用いる LMB2004 は, 奇しくも選挙に関する研究である. 先にこの論文の要旨を簡単に述べておく.

選挙に出る候補者は, 選挙に勝ちたい. 選挙に勝つにはより多くの有権者の支持が必要で, より多くの支持を集めるにはより多くの有権者が賛成するような政策を実施する公約を掲げれば良い. ならば, どんなで政党であっても, 候補者は得票数を最大化するためにある特定の政策を掲げることになり, 候補者の唱える政策は全員同じになるという結論が出る. つまり, 例えば2名の候補者がいて, それぞれまったく正反対のスタンスをとる政党に所属していたとしても, 得票数のために両者のちょうど中間的な方針 (中庸な政策) を選択することになる, という事を意味する. この, 候補者の政策がある一点に収束する状態を, 提唱者の名前から「ダウンズ均衡 (Downsian equilibrium)」と呼ぶ. しかし, 現実には政党によって政策は異なる. これは本来はアメリカのような2大政党制の国家を前提としたものだが, 共和党と民主党の政策は対照的で, 一致しない. この理論に対する実証上の裏付けは, まだ乏しい. LMB2004 によると, ダウンズの提唱したモデルでは, 選挙の争点が1つだけしかない. 複数の争点が存在する場合も, このモデルが成り立たない可能性がある. よって, この理論は前提条件が厳しく, 頑健に欠けているという批判があった. これはその後の研究で解消されている. 次の批判意見は, 実際には政治家・政党は独自のイデオロギーを持っているから, 一旦選挙に勝ってしまえば, もう彼らを阻むものはなく当初掲げていた中庸な政策から, 自らのイデオロギーの方向に戻っていってしまうというものである. 有権者は特定のイデオロギーを持った候補者を「選ぶ」ことしかできずに, 自ら積極的に政策に「影響を与え」られることはない.

もちろん, 現実には選挙は何度も行われるから, 露骨に公約違反をすれば信用を失い, 次回から票を集めることは難しくなるだろう. つまり, 政治家の信用とは, 方針を一貫して選挙に出馬し続けることで, 有権者から信用を獲得することができる. よって, この場合, ダウンズ均衡は発生せず, 選挙のたびに各候補者は多様な政策を唱えることになる. さらに, 候補者と有権者の信頼関係によっても政策の多様性は発生する.

LMB2004 では, ダウンズの言うように, 有権者は候補者を「選ぶ」だけしかできないのか, それとも候補者と政策にともども「影響を与え」られるのか, という検証のために RDD のフレームワークを利用している.

候補者の選挙前の政策の方針と, 選挙後の政策の情報を比較する. 両者に違いがある (つまり, 当選後に中庸な政策から自分のイデオロギーへ立ち返っている) なら, 候補者は当選後政策を変えたことになり, 有権者は「候補者を選んだ」だけになる. 一方, 違いが見られないなら, 有権者は政策本位の候補者選びをしている. しかし, 実際には, 候補者の持つ「本当のイデオロギー」を知る方法はないので, 工夫が必要である.

そこで, 以下の3式をそれぞれ推定することで, elect/affect の効果の大きさをそれぞれ求めることになる.

\begin{equation} \mathit{ADA}_{t+1} =\alpha_{1} + \beta_{0} P_{t} + \gamma_{1} D_{t} +\varepsilon_{t} \tag{1} \end{equation}

被説明変数は, ADA; American Democratic Action (http://www.adaction.org) という団体の発表するスコアとする (以下, ADAスコアと呼ぶ). これは20程度の社会・経済の問題に対する回答から, その候補者の政策志向の「リベラルさ」スコアづけするものである. 高いほど「リベラル」≒民主党寄りの政策を主張していることになる*5. P_{t} は選挙のあるt年の民主党候補の得票比率 (0~1) で, D_{t} はその候補が当選したかどうかのダミー変数である. つまり, 得票率 P_{t} が running variable であり, 非連続な変化の発生するカットオフポイントは 0.5 である. 選挙後の候補者のスタンスを見ることで, 落選・当選によって, その後の候補者のスタンスがどう変化したかを見ることができる. この変化の大きさが, 式 (1) の ATT である \gamma で, これは選挙による総変化効果 (overall effect) と呼べる.

次に, 被説明変数を選挙後ではなく, 選挙時点にしたものも推定する.

\begin{equation} \mathit{ADA}_{t} =\alpha_{2} + \pi_{0} P_{t} + \pi_{1} D_{t} +\varepsilon_{t} \tag{2} \end{equation}

こちらは, 当選・落選を受けての変化ではないので, \pi_{1} は elect の効果の一部を表すことができる.

さらに第3の式, 次の選挙での民主党候補の当選確率を被説明変数としたモデルを推定する.

\begin{equation} \mathit{P}_{t+1}^{D} =\alpha_{3} + \delta_{0} P_{t} + \eta_{1} D_{t} +\varepsilon_{t} \tag{3} \end{equation}

ここでの \eta はある政策を主張したときに, 選挙で勝てる確率の大きさを表す. これも elect 効果の一部を表すのに使用する.

LMB2004 はこれら3式を推定して得たパラメータから, 総変化 \gamma を elect の効果と affect の効果に分解できるとした. その分解式は, 以下の式 (4-6) で表される.

\begin{equation} \gamma= \mathrm{E}\left[ \mathit{ADA}_{t+1}| D_{t}=1 \right] -\mathrm{E}[\mathit{ADA}_{t+1} ] \\ = \pi_{0}(P_{t+1}^{\ast D} - P_{t+1}^{\ast R}) + \pi_{1}(P_{t+1}^D - P_{t+1}^R) \tag{4} \end{equation}

と,

\begin{equation} \mathrm{E}[\mathit{ADA}_{t} | D_{t}=1 ] - \mathrm{E}[ \mathit{ADA}_{t} | D_{t} = 0 ] = \pi_{1} \tag{5} \end{equation} \begin{equation} \eta = \mathrm{E}[ D_{t+1} | D_{t} =1 ] - \mathrm{E}[ D_{t+1} |D_{t}=0 ]= P_{t+1}^D - P_{t+1}^R \tag{6} \end{equation}

から導ける. 式 (4) の第1項 \pi_{0}(P_{t+1}^{\ast D} - P_{t+1}^{\ast R}) は affect 効果で, 第2項 \pi_{1}(P_{t+1}^D - P_{t+1}^R)が elect 効果. 2つの P_{t}^{\ast X} はそれぞれ, 民主党 (D), 共和党 (R) が現職候補であるときの, 民主党候補の当選確率を表す. これらは観測できないため, affect 効果は式 (4) のそれ以外の部分を推定することで, 間接的に大きさを導ける.

なお, ここではいきなり単純な線形式を提示しているが, 元の論文ではより厳密で複雑な議論を経てこの式を導いていることに注意しなければならない. 今回は RDD の応用例ということを主眼においているので, この過程は省略した.

RDestimate() 関数

ようやく rdd の使用法である. 以上から, LMB2004 のモデルは Sharp RD の手法で推定できるため, RDestimate() 関数を用いることができる. データは EITM Summer Institute 2011 から入手できる. これも .dta 形式なので, foreign パッケージが必要である.

library(foreign) # ver. 0.8-61
library(rdd) # ver. 0.5.6
df <- read.dta(file = "source/LMB Data.dta") #パスは適宜変更

gamma.result <- RDestimate(formula = score~lagdemvoteshare, data=df, cluster = df$id2, cutpoint = .5)
summary(gamma.result)
pi1.result   <- RDestimate(formula = score~demvoteshare, data = df,  cluster = df$id2, cutpoint = .5)
summary(pi1.result)
eta.result   <- RDestimate(formula = democrat~lagdemvoteshare, data=df, cluster = df$id2, cutpoint = .5)
summary(eta.result)

gamma.result, pi1.result, eta.result でそれぞれ \gamma, \pi_{1}, \eta の推定結果が得られる. 今回は他に説明変数のないモデルだったので, formula= には被説明変数と, 説明変数 (running variable) だけでよい. 説明変数を追加したり, Fuzzy RD にしたい場合は, formula=被説明変数~running variable + 操作変数 | その他の説明変数 + ... となる.

cutpoint= に与えられたカットオフポイントを元にダミー変数を判断するので, formula= にダミー変数は不要である. cutpoint= はデフォルトではゼロである.

cluster= はクラスタリングSEを計算するためのオプション引数である.

デフォルトでは Imbens and Kalyanaraman の最適バンド幅に従ってバンド幅を決定するが, bw= オプション引数に数値を与えて任意にバンド幅を決めることもできる(bw=h に対して, カットオフポイントからプラスマイナス h/2 の範囲のデータで推定する.). 設定したバンド幅に対し, バンド幅が半分の場合と, 2倍の場合の結果も自動で推定してくれる. 例えば \gamma の推定結果を summary() に与えた時の

Estimates:
           Bandwidth  Observations  Estimate  Std. Error  z value  Pr(>|z|)      
LATE       0.12883     6072         18.84     2.265        8.316   9.096e-17  ***
Half-BW    0.06442     3150         20.48     3.348        6.117   9.564e-10  ***
Double-BW  0.25767    10512         22.45     1.667       13.468   2.402e-41  ***

で, 上から順に最適バンド幅, 最適幅の半分 (Half-BW), 2倍の場合 (Double-BW) の結果が出力される. LMB2004 では, バンド幅は 0.04 に対し, \gamma=21.2 という結果を出している(元論文 Table I)*6. 多くの場合, バンド幅を狭くするほど推定値が大きくなるが, 今回は2倍バンド幅で逆に増加している. ただ, バンド幅のとり方で大きく結果が変わる, というほどでもなさそうだ.

なお, 元論文では, 候補者の収入, 学歴などの個人情報のデータも用いたり, 年代別に分けて推定するなど, より細かい分析をおこなっている.

グラフで表す

順番が前後してしまうが, REestimate() で得られる推定結果 (RD オブジェクト) を plot() に代入すると, グラフで表示してくれる. ただし, かなり簡易的な機能しかなく, 基本的にタイトルを付けたりと言ったオプションはない. ggplot のように見栄えにこだわったグラフを作成することはできない.

plot(gamma.result)
plot(pi1.result)
plot(eta.result)

\gamma のグラフ

f:id:ill-identified:20141218222659p:plain

\pi_{1} のグラフ

f:id:ill-identified:20141218222731p:plain

\eta のグラフ

f:id:ill-identified:20141218222749p:plain

正直, 融通が効かないので, 使用するデータによっては見づらくなる. 一応, 使用するデータの粒度を調整する gran= という引数があるが, 推定結果ではなく, データセットの2変数を放り込んでプロットしたほうが見やすいかもしれない.

plot(df$lagdemvoteshare, df$score, xlab="vote share at t", ylab="ADA at t+1", main=expression(gamma), pch=".")
plot(df$demvoteshare, df$score, xlab="vote share at t", ylab="ADA at t", main=expression(pi[1]), pch=".")
plot(df$lagdemvoteshare, df$democrat, xlab="vote share at t", ylab="eta at t+1", main=expression(eta), pch="." )

f:id:ill-identified:20141218224937p:plain

f:id:ill-identified:20141218224945p:plain

f:id:ill-identified:20141218224951p:plain

多項式近似の場合

もうひとつの推定方法, 多項式近似についてもやってみよう. 前回説明したとおり, 多項式近似は古典的な方法の応用のため, 直感的に分かりやすい, 非連続な変化の大きさを検定しやすい, という利点がある. つまり, 線形式の当てはめでは非連続点が存在するように見えても, 多項式近似を行えば, 実はカットオフポイントが滑らかな曲線で接続できる, ということがありうる.

f:id:ill-identified:20141218012205p:plain

多項式近似の場合は, lm() 関数を使い, AIC() で AIC を比較するだけで良い*7.

df$demvoteshare_c <- df$demvoteshare -.5 # カットオフポイントまわりに変換
df$lagdemvoteshare_c <- df$lagdemvoteshare -.5

# 各多項式について, AIC を最小化するように次数を決める

# gamma
for (i in 2:20){ 
  aic <- AIC(poly.gamma <- lm(formula = score~poly(lagdemvoteshare_c, i)  + lagdemocrat, data=df[!is.na(df$lagdemvoteshare_c),]))
  if (i == 2) aic.temp <- aic
  else if (aic <= aic.temp)  aic.temp <- aic
  else {
    if(i == 3) print("2-nd or lower order polynomial minimize the AIC.")
    else if (i < 20) print( paste(i-1,"-th polynomial has minimum the AIC"),sep="" )
    else print("20-th or higher order polynomial minimize the AIC ")
    
    break
  }
}
# 8次が選択された

#pi1
for (i in 2:20){ 
  aic <- AIC(poly.pi1 <- lm(formula = score~poly(demvoteshare_c, i)  + democrat, data=df[!is.na(df$demvoteshare_c),]))
  if (i == 2) aic.temp <- aic
  else if (aic <= aic.temp)  aic.temp <- aic
  else {
    if(i == 3) print("2-nd or lower order polynomial minimize the AIC.")
    else if (i < 20) print( paste(i-1,"-th polynomial has minimum the AIC"),sep="" )
    else print("20-th or higher order polynomial minimize the AIC ")
    
    break
  }
}
# 8次が選択された

# eta 
for (i in 2:20){ 
  aic <- AIC(poly.eta <- lm(formula = democrat~poly(lagdemvoteshare_c, i)  + lagdemocrat, data=df[!is.na(df$lagdemvoteshare_c),]))
  if (i == 2) aic.temp <- aic
  else if (aic <= aic.temp)  aic.temp <- aic
  else {
    if(i == 3) print("2-nd or lower order polynomial minimize the AIC.")
    else if (i < 20) print( paste(i-1,"-th polynomial has minimum the AIC"),sep="" )
    else print("20-th or higher order polynomial minimize the AIC ")
    
    break
  }
}
# 9次が選択された

以上, Imbens and Kalyanaraman によるバンド幅を選んだ場合と, 多項式の場合の結果を, 元論文の Table I と比較した結果が以下である. なお, 多項式でも, 3つのパラメータは全て 0.1 %水準で有意であった.

\gamma \pi_{1} \eta elect ( \pi_{1} \times \eta) affect ( \gamma -\mathit{elect})
元論文 21.2 47.6 0.48 22.84 -1.64
今回(LLR) 18.8 47.1 0.42 19.78 -0.94
今回(多項式) 18.3 46.8 0.43 20.01 -1.73

つまり, 選挙による候補者の政策の変化に占める affect 効果の大きさはかなり小さく, 大半が elect 効果によるものである, という結論となった.

余談

たまたま日本の衆院選挙と重なり時事ネタっぽくなったが, 偶然である. なお, LMB2004 の研究は, 米国の民主党・共和党が対立軸を明確にしていること, 2大政党制であることが可能にしている. 日本の場合, 自民党と多くの野党が存在する状態で, (個人的な感想として) 政策の対立軸も明確でないことが多い. また, 3つ以上の政党の対立をモデル化する場合, モデルははるかに複雑になるので (そもそも2大政党制だから説明変数を得票率に, カットオフポイントを .5 と言うふうに簡単にすることができた.), 日本のデータで応用するにはもうひとひねり必要だろう.

rdrobust パッケージ

rdrobust パッケージというのも存在する. こちらは Calonico, Cattaneo, Titiunik らが作成したもので, 作者らが Econometrica に投稿した "Robust Nonparametric Confidence Intervals for Regression-Discontinuity Designs" の最新の手法が実装されている. ただし, 説明変数に running variable しか使えないようだ. まだあまり調べられていないので, 解説は別の機会に.

参考文献

  • 黒澤昌子. (2005). 積極的労働政策の評価―レビュー. フィナンシャル・レビュー, (77).
  • 西口慶彦 (2014) 『ノンパラメトリック，セミパラメトリック計量経済分析：概観』 日本経済学会・春季大会 数量的経済・政策分析分科会「ノンパラメトリック, セミパラメトリック計量経済分析」セッション資料
  • Angrist, J. D., & Lavy, V. (1999). Using Maimonides’ Rule to Estimate the Effect of Class Size on Scholastic Achievement. Quarterly Journal of Economics, 114(2), 533–575.

  • Angrist, J. D., V. Lavy, J. Leder-Luis & A. Shany, Maimonides Rule Redux, National Bureau of Economic Research, Cambridge, MA, w23486.

  • --- & Pischke, J.-S. (2008). Mostly Harmless Econometrics: An Empricist’s Companion. Princeton University Press. 邦訳 大森義明・田中隆一・野口晴子・小原美紀 訳『ほとんど無害な計量経済学』 エヌティティ出版

  • Cameron, A. C., & Trivedi, P. K. (2005). Microeconometrics: Methods and Applications. Cambridge University Press.

  • ---, Gelbach, J. B., & Miller, D. L. (2006). Robust inference with multi-way clustering. NBER Working Paper, 0327.
  • Fan, J. (1992). Design-adaptive Nonparametric Regression. Journal of American Statistical Association, 87(420), 998–1004.
  • Hahn, J., Todd, P., & Van der Klaauw, W. (1999). Evaluating the effect of an antidiscrimination law using a regression-discontinuity design. NBER Working Paper, (7131).
  • ---, Todd, P., & Van der Klaauw, W. (2001). Identification and estimation of treatment effects with a regression‐discontinuity design. Econometrica, 69(1), 201–209.
  • Imbens, G. W., & Lemieux, T. (2008). Regression discontinuity designs: A guide to practice. Journal of Econometrics, 142(2), 615–635. doi:10.1016/j.jeconom.2007.05.001
  • ---, & Kalyanaraman, K. (2012). Optimal Bandwidth Choice for the Regression Discontinuity Estimator. Review of Economic Studies, 79(3), 933–959. doi:10.1093/restud/rdr043
  • Lee, D. S., Moretti, E., & Butler, M. J. (2004). Do Voters Affect or Elect Policies? Evidence from the U. S. House. Quarterly Journal of Economics, 119(3), 807–859. doi:10.1162/0033553041502153
  • ---, & Lemieux, T. (2010). Regression discontinuity designs in economics. Journal of Economic Literature, 48(2), 281–355.
  • McCrary, J. (2008). Manipulation of the running variable in the regression discontinuity design: A density test. Journal of Econometrics, (December). およびそのSTATAコード資料
  • Wooldridge, J. M. (2010). Econometric Analysis of Cross Section and Panel Data (2nd ed.). The MIT Press.
  • ill-identified (2014) 『[時事ネタ] 少人数教育は実際どれくらい意味があるのか』
  • --- (2014) 『[非連続回帰デザイン(Regression-Discontinuity Design)][ill2014rdd1]』

以下, R関連

*1:今回は final5.dta のみ使用

*2:R本体のバージョンは 3.1.2 である

*3:リンク先のページにある, AngristLavy_Table4.do のこと

*4:サンプルコードでは integer でなく logical 型にしている.

*5:元の論文では, これ以外の尺度でも同様の方法で推定している

*6:3本の式の推定ごとにバンド幅が異なるというのはある意味不合理だが, いまのところはこのような場合のバンド幅の決め方は研究されていない?

*7:poly 関数を使うと, 欠損値を含むデータを代入した場合エラーになるので注意. subset= で欠損値を除くのではなく, data= で代入する時点で欠損値のあるレコードを落とす必要がある. これが少し不便なのだが, もう少し楽な方法はないものか

ill-identified 読者になる

広告を非表示にする
コメントを書く