はじめに
今年のノーベル経済学賞を受賞したDavid Card、Joshua Angrist、Guido Imbens 3氏は、よくメディアによく取り上げられる最低賃金の他に、経済学を超えて、政治学と社会学などの他の社会科学、並びに自然科学やコンピューターサイエンスなど他の分野の学問に対しても、因果推論の手法の開発で大きく貢献しています。
そこで本文では、Guido ImbensがDonald Rubinと共著した論文Imbens, G. W. and D. B. Rubin. "Estimating Outcome Distributions for Compliers in Instrumental Variables Models." The Review of Economic Studies 64, no. 4 (1997): 555-574.の実装コードを紹介する。また、可視化の部分ではAbadie, Alberto. "Bootstrap Tests for Distributional Treatment Effects in Instrumental Variable Models." Journal of the American Statistical Association 97, no. 457 (2002): 284-292.を参考する。
手法の背景と概説
Imbens and Rubin(1997)を理解するためには、まず操作変数法という概念を明確にする必要がある。理想の実験では、処置群にアサインされた人（か物）が100%処置を受け、統制群にアサインされた人は100%処置を受けないことを条件とするが、残念ながら現実の世界でこれが実現することに大きな困難が伴う場合もある。例えばビタミン摂取の実験でビタミンを摂取しないでほしい対称群の人が、自分でサプリメントを薬局で購入したり、もしくは処置群にアサインされた人がルール通りにビタミンを摂取してくれないなどといった問題がある。
そこで、「アサイン行動」はビタミンの摂取量を完璧にコントロールできないため、理想の実験ではないが、「アサイン行動」がビタミン摂取量をある程度コントロールできれば、「アサイン行動」を操作変数として扱って、操作変数の手法で実験参加者が研究者の指示を完璧に従うわけではない状況から取ったデータから、理想の実験データを再現できる。
実験参加者の分類を考えよう。参加者の行動は以下の4つのパターンに分けられる。
| グループ | 処置群にアサインした場合 | 統制群にアサインした場合 |
|---|---|---|
| complier | 処置を受ける | 処置を受けない |
| defier | 処置を受けない | 処置を受ける |
| always taker | 処置を受ける | 処置を受ける |
| never taker | 処置を受けない | 処置を受けない |
| Complierが研究者が一番好きな参加者で、大人しく言われた通りに行動する。Defierはかなり厄介で、ビタミンを飲めと言われたら決してビタミンを飲まず、ビタミンを飲むなと言われたらビタミンを飲み始める人で、always takerとnever takerは文字通り研究者の言うことを全く無視する人である。実験参加者がどのグループに属しているかがわかれば話が早いが、「処置群にアサインした場合」と「統制群にアサインした場合」はいわゆる「反実仮想」なので、研究者は片方しか観測できないため、下の表のように、どんなにデータを集めても、必ず半分のデータが欠落する。 |
| 参加者番号 | 処置群にアサインした場合 | 統制群にアサインした場合 |
|---|---|---|
| 1 | ? | 処置を受ける |
| 2 | ? | 処置を受けない |
| 3 | 処置を受けない | ? |
| 4 | ? | 処置を受けない |
| 5 | 処置を受ける | ? |
Imbens and Rubin(1997)は、defierがいないというデータから検証できない仮定を置くと、complierグループの人が全員処置を受けた場合の結果変数の分布と全員処置を受けない場合の結果変数の分布がわかることを示した。
コードと実装
Imbens and Rubin(1997)が提案した手法は以下のコードに入っている関数で実装できる。
# Implement the method introduced in Imbens, G. W. and D. B. Rubin. "Estimating Outcome Distributions for Compliers in Instrumental Variables Models." The Review of Economic Studies 64, no. 4 (1997): 555-574. and visualization of Abadie, Alberto. "Bootstrap Tests for Distributional Treatment Effects in Instrumental Variable Models." Journal of the American Statistical Association 97, no. 457 (2002): 284-292.
# データフレームのカラム名は "outcome", "treatment"(二項変数), and "instrument"(二項変数)に設定してください。
outcome_distributions <- function(data, binwidth){
Y1_bar <- mean(subset(data, instrument == 1)$outcome)
Y0_bar <- mean(subset(data, instrument == 0)$outcome)
D1_bar <- mean(subset(data, instrument == 1)$treatment)
D0_bar <- mean(subset(data, instrument == 0)$treatment)
f00_df <- subset(data, instrument == 0 & treatment == 0)
f01_df <- subset(data, instrument == 0 & treatment == 1)
f10_df <- subset(data, instrument == 1 & treatment == 0)
f11_df <- subset(data, instrument == 1 & treatment == 1)
phi_n <- 1 - mean(subset(data, instrument == 1)$treatment)
phi_a <- mean(subset(data, instrument == 0)$treatment)
phi_c <- 1 - phi_n - phi_a
bin <- seq(min(data$outcome), max(data$outcome), by = binwidth)
F00 <- c()
F01 <- c()
F10 <- c()
F11 <- c()
for (i in 1:length(bin)){
F00[i] <- nrow(subset(f00_df, outcome < bin[i]))
F01[i] <- nrow(subset(f01_df, outcome < bin[i]))
F10[i] <- nrow(subset(f10_df, outcome < bin[i]))
F11[i] <- nrow(subset(f11_df, outcome < bin[i]))
}
g_c0 <- c()
g_c1 <- c()
for (i in 1:length(bin)){
g_c0[i] <- ((phi_n + phi_c)/(phi_c))*F00[i] - (phi_n/phi_c)*F10[i]
g_c1[i] <- ((phi_a + phi_c)/(phi_c))*F11[i] - (phi_a/phi_c)*F01[i]
}
g_c0_scaled <- c()
g_c1_scaled <- c()
for (i in 1:length(g_c0)){
g_c0_scaled[i] <- g_c0[i]/max(g_c0)
g_c1_scaled[i] <- g_c1[i]/max(g_c1)
}
return (list(Y1_bar = Y1_bar, Y0_bar = Y0_bar, D1_bar = D1_bar, D0_bar = D0_bar,
phi_n = phi_n, phi_a = phi_a, phi_c = phi_c,
F00 = F00, F10 = F10, F01 = F01, F11 = F11,
g_c0 = g_c0_scaled, g_c1 = g_c1_scaled, bin = bin))
}
詳細はまた別の記事で紹介する。
ここでは早速アメリカの職業訓練データでImbens and Rubin(1997)が提案した手法を実装する。
library(haven)
# データをhttps://cran.r-project.org/web/packages/ivdesc/ivdesc.pdfが説明するようにダウンロード
jtpa <- data.frame(read_dta("http://fmwww.bc.edu/repec/bocode/j/jtpa.dta"))
# ジェンダーと年齢でデータを抽出
jtpa_male <- subset(jtpa, sex == 1)
jtpa_female <- subset(jtpa, sex == 0)
jtpa_young <- subset(jtpa, age <= 31)
jtpa_old <- subset(jtpa, age > 31)
# outcome_distribution_function.Rに対応するための簡単な前処理
data_jtpa <- na.omit(data.frame(outcome = log(jtpa$earnings + 1), treatment = jtpa$training,
instrument = jtpa$assignmt))
data_jtpa_male <- na.omit(data.frame(outcome = log(jtpa_male$earnings + 1), treatment = jtpa_male$training,
instrument = jtpa_male$assignmt))
data_jtpa_female <- na.omit(data.frame(outcome = log(jtpa_female$earnings + 1), treatment = jtpa_female$training,
instrument = jtpa_female$assignmt))
data_jtpa_young<- na.omit(data.frame(outcome = log(jtpa_young$earnings + 1), treatment = jtpa_young$training,
instrument = jtpa_young$assignmt))
data_jtpa_old <- na.omit(data.frame(outcome = log(jtpa_old$earnings + 1), treatment = jtpa_old$training,
instrument = jtpa_old$assignmt))
# complierが全員処置を受ける場合と全員受けない場合の結果変数の分布をoutcome_distribution_function.Rのoutcome_distributions関数で推定（学習）
od_jtpa <- outcome_distributions(data_jtpa, binwidth = 0.01)
od_jtpa_male <- outcome_distributions(data_jtpa_male, binwidth = 0.01)
od_jtpa_female <- outcome_distributions(data_jtpa_female, binwidth = 0.01)
od_jtpa_young <- outcome_distributions(data_jtpa_young, binwidth = 0.01)
od_jtpa_old <- outcome_distributions(data_jtpa_old, binwidth = 0.01)
plot(x = od_jtpa$bin, y = od_jtpa$g_c1, type = "l", col = "blue",
main = "counterfactural outcome distribution\nestimated by instrumental variable",
xlab = "log(earnings)", ylab = "CDF")
lines(x = od_jtpa$bin, y = od_jtpa$g_c0, col = "red")
legend("topleft", legend = c("with training", "without training"),
col = c("blue", "red"), lty = c(1, 1))
plot(x = od_jtpa_male$bin, y = od_jtpa_male$g_c1, type = "l", col = "blue",
main = "counterfactural outcome distribution\nestimated by instrumental variable (male)",
xlab = "log(earnings)", ylab = "CDF")
lines(x = od_jtpa_male$bin, y = od_jtpa_male$g_c0, col = "red")
legend("topleft", legend = c("with training", "without training"),
col = c("blue", "red"), lty = c(1, 1))
plot(x = od_jtpa_female$bin, y = od_jtpa_female$g_c1, type = "l", col = "blue",
main = "counterfactural outcome distribution\nestimated by instrumental variable (female)",
xlab = "log(earnings)", ylab = "CDF")
lines(x = od_jtpa_female$bin, y = od_jtpa_female$g_c0, col = "red")
legend("topleft", legend = c("with training", "without training"),
col = c("blue", "red"), lty = c(1, 1))
plot(x = od_jtpa_young$bin, y = od_jtpa_young$g_c1, type = "l", col = "blue",
main = "counterfactural outcome distribution\nestimated by instrumental variable (young)",
xlab = "log(earnings)", ylab = "CDF")
lines(x = od_jtpa_young$bin, y = od_jtpa_young$g_c0, col = "red")
legend("topleft", legend = c("with training", "without training"),
col = c("blue", "red"), lty = c(1, 1))
plot(x = od_jtpa_old$bin, y = od_jtpa_old$g_c1, type = "l", col = "blue",
main = "counterfactural outcome distribution\nestimated by instrumental variable (old)",
xlab = "log(earnings)", ylab = "CDF")
lines(x = od_jtpa_old$bin, y = od_jtpa_old$g_c0, col = "red")
legend("topleft", legend = c("with training", "without training"),
col = c("blue", "red"), lty = c(1, 1))
上記のコードは何枚かの図を出力するが、ここでは全員分のデータで推定した図だけで見方を説明する。
ここで青い線はcomplierが全員処置を受けた場合の結果変数の累積分布で、赤い線はcomplierが全員処置を受けない場合の結果変数の累積分布である。例えば収入の対数値が6の時、青い線が赤い線の下にあるが、これは収入の対数値が6になった時、処置を受けない場合の結果(赤)が先に頭打ちを迎えることを意味するため、処置を受けた場合の結果の方が全体的に高い。また、中央値(y軸が0.5)を見れば、青い線(処置を受けた場合)のx軸の値が赤い線(処置を受けない場合)のx軸の値より高いこともわかる。
よって、Imbens and Rubin(1997)が提案した手法からは、アメリカの職業訓練は少なくとも実験のcomplierに対しては正の効果をもたらすと言える。