Download the sample data by clicking here.
\[ \text{ITT} = \frac{1}{N_{1}} \sum_{i=1}^N Y_i Z_i - \frac{1}{N - N_{1}} \sum_{i=1}^N Y_i (1 - Z_i). \]
d <- haven::read_dta("practice.dta")
# ITT
ITT <- mean(d$voted[d$assigned == 1]) - mean(d$voted[d$assigned == 0])
round(ITT, digits = 2)
## [1] 0.03
# ITT SE
seDiffMeans <- function(y, tx) {
y1 <- y[tx == 1]
y0 <- y[tx == 0]
n1 <- length(y1)
n0 <- length(y0)
sqrt(((var(y1)/n1 + var(y0)/n0)))
}
round(seDiffMeans(d$voted, d$assigned), digits = 2)
## [1] 0.01
ITT_LCI <- ITT - 1.96 * seDiffMeans(d$voted, d$assigned)
ITT_UCI <- ITT + 1.96 * seDiffMeans(d$voted, d$assigned)
round(c(ITT_LCI, ITT_UCI), digits = 2)
## [1] 0.01 0.04
Use the ivreg
function in the AER package. It is used as ivreg(outcome ~ treatment + controls | instruments + controls, data = yourdata)
.
# LATE via Wald estimator
library(AER)
## Loading required package: car
## Loading required package: carData
## Loading required package: lmtest
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: survival
model <- ivreg(voted ~ treated | assigned, data = d)
LATE <- summary(model)$coefficients[2, 1]
ci <- confint(model)
LATE_LCI <- unname(ci[2, 1])
LATE_UCI <- unname(ci[2, 2])
round(LATE, digits = 2)
## [1] 0.08
round(c(LATE_LCI, LATE_UCI), digits = 2)
## [1] 0.03 0.14
ivreg
function in the AER package.model2 <- ivreg(voted ~ treated + age65 | assigned + age65 , data = d)
stargazer::stargazer(model, model2, type = "html", title = "Estimation Results")
Dependent variable: | ||
voted | ||
(1) | (2) | |
treated | 0.081*** | 0.084*** |
(0.028) | (0.027) | |
age65 | 0.220*** | |
(0.008) | ||
Constant | 0.425*** | 0.372*** |
(0.004) | (0.004) | |
Observations | 20,000 | 20,000 |
R2 | 0.004 | 0.040 |
Adjusted R2 | 0.004 | 0.040 |
Residual Std. Error | 0.494 (df = 19998) | 0.485 (df = 19997) |
Note: | p<0.1; p<0.05; p<0.01 |
Download the sample data by clicking here.
library(stargazer)
##
## Please cite as:
## Hlavac, Marek (2018). stargazer: Well-Formatted Regression and Summary Statistics Tables.
## R package version 5.2.2. https://CRAN.R-project.org/package=stargazer
library(lmtest)
library(sandwich)
d <- read.csv("child_soldiering.csv")
# Use the regression estimator
mod1 <- lm(educ ~ abd + age, data = d)
mod1HC2 <- list(coeftest(mod1, vcov = vcovHC(mod1, type = "HC2"))[, 2])
stargazer(mod1, se = mod1HC2, title = "OLS Results", type = "html", keep.stat = c("n"), notes = c("HC2 Robust SEs"))
Dependent variable: | |
educ | |
abd | -0.634*** |
(0.227) | |
age | 0.031 |
(0.022) | |
Constant | 6.781*** |
(0.416) | |
Observations | 741 |
Note: | p<0.1; p<0.05; p<0.01 |
HC2 Robust SEs |
library(ggplot2)
# What is half the magnitude of the treatment coefficient?
b.half <- mod1$coefficients[2]/2
# The range of delta is 0,1
delta <- seq(0.001, 1, by = 0.01)
# Solve for gamma
gamma <- b.half/delta
# Make plot
ggplot(data.frame(delta = delta, gamma = gamma), aes(x = delta,
y = gamma)) + geom_path() + xlab(expression(delta)) + ylab(expression(gamma)) +
ylim(-3, 0)
## Warning: Removed 11 row(s) containing missing values (geom_path).
# Add covariates as benchmarks to the plot
# Estimate delta and gamma for covariates
d_cov <- c()
g_cov <- c()
covars <- c("orphan96")
for (covar in covars) {
d_cov[covar] <- mean(d[d$abd == 1, covar]) - mean(d[d$abd ==
0, covar])
g_cov[covar] <- mean(d[d[covar] == 1, "educ"]) - mean(d[d[covar] == 0, "educ"])
}
# Create data frame for plot
obs_covar <- data.frame(covar = covars, delta = abs(d_cov), gamma = -abs(g_cov))
# take absolute value for delta, and minus absolute value for
# gamma, because our unobserved confounder has a positive
# delta and negative gamma
# Make plot
ggplot(data.frame(delta = delta, gamma = gamma), aes(x = delta,
y = gamma)) + geom_path() + xlab(expression(delta)) + ylab(expression(gamma)) +
ylim(-3, 0) + geom_point(data = obs_covar, aes(x = delta, y = gamma)) +
geom_text(data = obs_covar, aes(label = covar), hjust = 0,
vjust = c(2), size = 4)
## Warning: Removed 11 row(s) containing missing values (geom_path).
Graduate School of Economics, Waseda University, ritsu.kitagawa@fuji.waseda.jp↩