1 Load LaLonde’s experimental (randomized employment) data & PSID potential controls

nsw_df  <- read_dta("http://users.nber.org/~rdehejia/data/nsw.dta")
psid_df <- read_dta("http://users.nber.org/~rdehejia/data/psid_controls.dta")

helper function to make all standard errors robust

robust <- function(m) coeftest(m, vcov = vcovHC(m, type = "HC1"))

1.2 Pre-treatment earnings on treatment

m12 <- lm(re75 ~ treat, data = nsw_df)
robust(m12)
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3026.683    252.350  11.994   <2e-16 ***
## treat         39.415    379.009   0.104   0.9172    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This is like summary statistics in that it tells us that treatment group’s real earnings are only $39.415 more than the control group ($3026.683) on average. This is within the standard error of $379.009.

1.3 and 1.4 1978 earnings on treatment

m13 <- lm(re78 ~ treat, data = nsw_df)
robust(m13)
## 
## t test of coefficients:
## 
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5090.05     277.43 18.3474  < 2e-16 ***
## treat         886.30     488.14  1.8157  0.06984 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The program has a positive impact on earnings ($886.30, which is larger than the standard error of 488.14), but it is only statistically significant at the 0.1 level.

For this to be causal, we need there to be exogeneity, with treatment assignment, D, independent of potential outcomes, Y. That way we can have an accurate counterfactual in our control group. This is the case here because treatment is randomized, so the difference of means estimates the ATE (and ATT).

1.5 and 1.6 Add controls for age, education, and race

m15 <- lm(re78 ~ treat + age + education + black + hispanic, data = nsw_df)
robust(m15)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  2848.033   1917.945  1.4849  0.13800  
## treat         824.395    480.482  1.7158  0.08664 .
## age            30.552     32.114  0.9514  0.34174  
## education     291.923    138.726  2.1043  0.03570 *
## black       -1827.827    772.968 -2.3647  0.01831 *
## hispanic     -148.790    982.629 -0.1514  0.87969  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Since treatment was random, all of these demographics were random too in our sample. So conditioning on controls that were already random doesn’t really change our estimate (shown as treat). These are decently bad controls, just because they are irrelevant here. A bad control can also be something that controls for the mechanism that causes the treatment effect.

1.7 Is nodegree causal? Is nodegree X treatment?

m17 <- lm(re78 ~ treat * nodegree, data = nsw_df)
robust(m17)
## 
## t test of coefficients:
## 
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     5918.61     688.66  8.5944   <2e-16 ***
## treat           1058.70    1100.32  0.9622   0.3363    
## nodegree       -1017.74     751.94 -1.3535   0.1763    
## treat:nodegree  -352.24    1225.74 -0.2874   0.7739    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

The nodegree variable can’t be interpreted as causal, because it wasn’t randomly assigned. Treatment X nodegree could be interpreted as heterogeneity for education amongst the treated people, comparing the treatment effects between the nodegree group and the people with a degree group. This can be interpreted as causal, because treatment was randomly assigned.

1.8 Difference in means VS results from 1.3

means_nsw <- nsw_df %>%
  group_by(treat) %>%
  summarize(mean_re78 = mean(re78), n = n())
means_nsw
## # A tibble: 2 × 3
##   treat mean_re78     n
##   <dbl>     <dbl> <int>
## 1     0     5090.   425
## 2     1     5976.   297
diff(means_nsw$mean_re78)
## [1] 886.3037

1.9 Treated from NSW and control from PSID

mixed_df <- bind_rows(
  nsw_df  %>% filter(treat == 1),
  psid_df %>% filter(treat == 0) %>% select(-re74)
)
table(mixed_df$data_id, mixed_df$treat)
##                 
##                     0    1
##   Lalonde Sample    0  297
##   PSID           2490    0

1.10 Difference in means from new mixed dataset VS NSW dataset

means_mixed <- mixed_df %>%
  group_by(treat) %>%
  summarize(mean_re78 = mean(re78), n = n())
means_mixed
## # A tibble: 2 × 3
##   treat mean_re78     n
##   <dbl>     <dbl> <int>
## 1     0    21554.  2490
## 2     1     5976.   297
diff(means_mixed$mean_re78)
## [1] -15577.57

1.11

mixed_df %>%
  group_by(treat) %>%
  summarize(across(c(age, education, black, hispanic, married, nodegree, re75, re78), mean))
## # A tibble: 2 × 9
##   treat   age education black hispanic married nodegree   re75   re78
##   <dbl> <dbl>     <dbl> <dbl>    <dbl>   <dbl>    <dbl>  <dbl>  <dbl>
## 1     0  34.9      12.1 0.251   0.0325   0.866    0.305 19063. 21554.
## 2     1  24.6      10.4 0.801   0.0943   0.168    0.731  3066.  5976.

The signs are very different. This is because our counterfactual isn’t an accurate counterfactual for this group anymore, are you can see from the summary statistics. This PSID group is wealthier. So if the control group is already making more money than the treatment group for reasons unrelated to the treatment, then the model is accidentally going to be measuring that baseline earnings gap instead of measuring the effect of the new program on earnings.

1.12 Matching (nearest neighbor)

covs <- c("age", "education", "black", "hispanic", "married", "nodegree")

X_t <- as.matrix(mixed_df[mixed_df$treat == 1, covs])
X_c <- as.matrix(mixed_df[mixed_df$treat == 0, covs])
y_t <- mixed_df$re78[mixed_df$treat == 1]
y_c <- mixed_df$re78[mixed_df$treat == 0]

nn_att <- function(D, y_t, y_c) {
  y0_hat <- apply(D, 1, function(d) mean(y_c[d == min(d)]))
  list(att = mean(y_t - y0_hat), y0_hat = y0_hat)
}

1.12.1 Euclidian matched ATE estimate

euclid <- function(A, B) {
  sqrt(pmax(outer(rowSums(A^2), rowSums(B^2), "+") - 2 * A %*% t(B), 0))
}
D_euc   <- euclid(X_t, X_c)
nn_euc  <- nn_att(D_euc, y_t, y_c)
nn_euc$att
## [1] -6781.174

1.12.2 Mahalanobis matched ATE estimate

D_mah  <- mahalanobis.dist(data.x = X_t, data.y = X_c)   # StatMatch
nn_mah <- nn_att(as.matrix(D_mah), y_t, y_c)
nn_mah$att
## [1] -3985.742

1.12.3

Matched estimates are not as large of negative numbers as with the non-matched dataset.

1.13.1

ps_logit <- glm(treat ~ age + education + black + hispanic + married + nodegree,
                family = binomial, data = mixed_df)
summary(ps_logit)
## 
## Call:
## glm(formula = treat ~ age + education + black + hispanic + married + 
##     nodegree, family = binomial, data = mixed_df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  0.64322    0.87228   0.737    0.461    
## age         -0.10149    0.01194  -8.503  < 2e-16 ***
## education   -0.04369    0.05351  -0.816    0.414    
## black        2.16986    0.22689   9.564  < 2e-16 ***
## hispanic     2.40421    0.38345   6.270 3.61e-10 ***
## married     -2.56670    0.19010 -13.502  < 2e-16 ***
## nodegree     1.30217    0.25445   5.118 3.09e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1891.12  on 2786  degrees of freedom
## Residual deviance:  890.39  on 2780  degrees of freedom
## AIC: 904.39
## 
## Number of Fisher Scoring iterations: 7

Age, black, hispanic, married, and nodegree are all predictive of treatment.

1.13.2 Add estimated propensity scores to the mixed dataset

mixed_df$p_hat <- ps_logit$fitted.values

mixed_df %>% group_by(treat) %>%
  summarize(min = min(p_hat), p25 = quantile(p_hat, .25), median = median(p_hat),
            p75 = quantile(p_hat, .75), max = max(p_hat))
## # A tibble: 2 × 6
##   treat      min     p25  median    p75   max
##   <dbl>    <dbl>   <dbl>   <dbl>  <dbl> <dbl>
## 1     0 0.000262 0.00250 0.00623 0.0384 0.873
## 2     1 0.00261  0.328   0.643   0.837  0.903
ggplot(mixed_df, aes(x = p_hat, fill = factor(treat))) +
  geom_histogram(bins = 40, alpha = 0.6, position = "identity") +
  labs(x = expression(hat(p)[i]), y = "Count", fill = "Treated") +
  theme_minimal()

There is overlap in [0.0026, 0.873], but the amount of overlap is pretty small because the IQR of the control is pretty far from the IQR of the treatment group, which has a much higher median.

1.13.3 Enforce overlap

p_min_t <- min(mixed_df$p_hat[mixed_df$treat == 1])
p_max_c <- max(mixed_df$p_hat[mixed_df$treat == 0])
c(p_min_t, p_max_c)
## [1] 0.002614915 0.872671549
overlap_df <- mixed_df %>% filter(p_hat >= p_min_t, p_hat <= p_max_c)
table(overlap_df$treat)
## 
##    0    1 
## 1854  265
ggplot(overlap_df, aes(x = p_hat, fill = factor(treat))) +
  geom_histogram(bins = 40, alpha = 0.6, position = "identity") +
  labs(x = expression(hat(p)[i]), y = "Count", fill = "Treated",
       title = "After enforcing overlap") +
  theme_minimal()

1.13.4 Treatment effect with regression conditioned on \(\hat{p}_{i}\) (centered and uncentered)

m1134a <- lm(re78 ~ treat + p_hat, data = overlap_df)
robust(m1134a)
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)  20737.72     333.46  62.1902 < 2.2e-16 ***
## treat        -7039.79     729.45  -9.6507 < 2.2e-16 ***
## p_hat       -14388.69    1099.16 -13.0906 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
overlap_df <- overlap_df %>%
  mutate(p_c = p_hat - mean(p_hat[treat == 1]))
m1134b <- lm(re78 ~ treat * p_c, data = overlap_df)
robust(m1134b)
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)  10369.68     613.78  16.8948 < 2.2e-16 ***
## treat        -4146.51     752.71  -5.5088 4.053e-08 ***
## p_c         -20829.89    1392.63 -14.9573 < 2.2e-16 ***
## treat:p_c    18313.37    1989.96   9.2029 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
robust(lm(re78 ~ treat * p_hat, data = overlap_df))
## 
## t test of coefficients:
## 
##              Estimate Std. Error  t value  Pr(>|t|)    
## (Intercept)  21190.57     346.93  61.0799 < 2.2e-16 ***
## treat       -13660.09     955.81 -14.2916 < 2.2e-16 ***
## p_hat       -20829.89    1392.63 -14.9573 < 2.2e-16 ***
## treat:p_hat  18313.37    1989.96   9.2029 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Interacting \(\hat{p}_{i}\) tells us that earnings among NSW participants don’t vary nearly as much with \(\hat{p}_{i}\) as earnings among PSID controls do.

1.13.5

breaks <- c(p_min_t, seq(0.1, 0.9, by = 0.1), p_max_c)
breaks <- sort(unique(breaks[breaks >= p_min_t & breaks <= p_max_c]))
overlap_df$block <- cut(overlap_df$p_hat, breaks = breaks, include.lowest = TRUE)

block_tab <- overlap_df %>%
  group_by(block) %>%
  summarize(n_t = sum(treat == 1), n_c = sum(treat == 0),
            diff = mean(re78[treat == 1]) - mean(re78[treat == 0]),
            .groups = "drop")
block_tab
## # A tibble: 9 × 4
##   block           n_t   n_c    diff
##   <fct>         <int> <int>   <dbl>
## 1 [0.00261,0.1]    33  1582 -13910.
## 2 (0.1,0.2]        23    98  -7919.
## 3 (0.2,0.3]        14    44  -4356.
## 4 (0.3,0.4]        21    36  -4798.
## 5 (0.4,0.5]        29    26 -11165.
## 6 (0.5,0.6]        20    26  -6316.
## 7 (0.6,0.7]        18     9   -519.
## 8 (0.7,0.8]        44    14  -4581.
## 9 (0.8,0.873]      63    19   1362.
# keep blocks that contain both treated and controls
block_tab <- block_tab %>% filter(n_t > 0, n_c > 0)
att_block <- with(block_tab, sum(diff * n_t) / sum(n_t))              # ATT weights
ate_block <- with(block_tab, sum(diff * (n_t + n_c)) / sum(n_t + n_c)) # ATE weights
c(ATT = att_block, ATE = ate_block)
##        ATT        ATE 
##  -5200.525 -11808.186

1.13.6

The propensity score estimates weren’t very close to the original experiment estimate. This is because the mixed dataset doesn’t have the same random assignment into treatment that makes a control/treatment comparison valid.

2.1

The comparison of the impact of class size on education would be biased upwards, because wealthy families are more able to get their kids into schools that have smaller class sizes. A simple comparison would end up measuring the impact of wealth on education, not class size.

2.2

The paper randomized which students got to be in small versus larger classes within each school. The main result was that students with smaller class sizes were more likely to be enrolled in college at age 20, which had ripple effects on their lives in terms of their finances (for the better).

2.3

A binned scatterplot groups data points into bins and plots the density of the bin rather than every specific point within that bin. Basically, first you make your bins (x axis). Then you find the average x and y of each bin and plot that.

2.4.1 Binned scatter plots

source("BM_StandardErrors.R")

grade5 <- readRDS("grade5.rds")

cut <- 40.5

plot_df <- grade5 %>% filter(school_enrollment >= 20, school_enrollment <= 60)

rdplot(y = plot_df$classize, x = plot_df$school_enrollment, c = cut, p = 1,
       title = "Class size at the 40-student threshold",
       x.label = "5th-grade school enrollment", y.label = "Class size")
## Warning in rdplot(y = plot_df$classize, x = plot_df$school_enrollment, c = cut,
## : Mass points detected in the running variable.

2.4.2 More binned scatter plots

rdplot(y = plot_df$avgmath, x = plot_df$school_enrollment, c = cut, p = 1,
       title = "Average math score", x.label = "5th-grade school enrollment",
       y.label = "Average math score")
## Warning in rdplot(y = plot_df$avgmath, x = plot_df$school_enrollment, c = cut,
## : Mass points detected in the running variable.

rdplot(y = plot_df$avgverb, x = plot_df$school_enrollment, c = cut, p = 1,
       title = "Average verbal score", x.label = "5th-grade school enrollment",
       y.label = "Average verbal score")
## Warning in rdplot(y = plot_df$avgverb, x = plot_df$school_enrollment, c = cut,
## : Mass points detected in the running variable.

2.4.3 More binned scatter plots

rdplot(y = plot_df$disadvantaged, x = plot_df$school_enrollment, c = cut, p = 1,
       title = "% disadvantaged", x.label = "5th-grade school enrollment",
       y.label = "% disadvantaged")
## Warning in rdplot(y = plot_df$disadvantaged, x = plot_df$school_enrollment, :
## Mass points detected in the running variable.

rdplot(y = plot_df$religious, x = plot_df$school_enrollment, c = cut, p = 1,
       title = "Religious school", x.label = "5th-grade school enrollment",
       y.label = "Share religious schools")
## Warning in rdplot(y = plot_df$religious, x = plot_df$school_enrollment, : Mass
## points detected in the running variable.

rdplot(y = plot_df$female, x = plot_df$school_enrollment, c = cut, p = 1,
       title = "Fraction female", x.label = "5th-grade school enrollment",
       y.label = "Fraction female")
## Warning in rdplot(y = plot_df$female, x = plot_df$school_enrollment, c = cut, :
## Mass points detected in the running variable.

2.4.4 Histogram of schools by enrollment, then with threshold of 40

schools <- grade5 %>% count(schlcode, school_enrollment)   # one row per school
nrow(schools)
## [1] 1002
ggplot(schools, aes(x = school_enrollment)) +
  geom_histogram(binwidth = 1, boundary = 0.5, fill = "grey60", colour = "white") +
  geom_vline(xintercept = c(40.5, 80.5, 120.5), linetype = "dashed", colour = "hotpink") +
  labs(x = "5th-grade school enrollment", y = "Number of schools") +
  theme_minimal()

ggplot(filter(schools, school_enrollment <= 80), aes(x = school_enrollment)) +
  geom_histogram(binwidth = 1, boundary = 0.5, fill = "grey60", colour = "white") +
  geom_vline(xintercept = 40.5, linetype = "dashed", colour = "hotpink") +
  labs(x = "5th-grade school enrollment", y = "Number of schools") +
  theme_minimal()

2.5 Regression Analysis

rd_df <- grade5 %>%
  filter(school_enrollment > 0, school_enrollment <= 80) %>%
  mutate(above40     = as.numeric(school_enrollment > 40),
         x_c         = school_enrollment - 40,
         x_c_above40 = x_c * above40)

rd_reg <- function(yvar) {
  f   <- as.formula(paste(yvar, "~ above40 + x_c + x_c_above40"))
  m   <- lm(f, data = rd_df)
  se  <- BMlmSE(m, clustervar = as.factor(rd_df$schlcode))$se.Stata["above40"]
  b   <- coef(m)["above40"]
  data.frame(outcome = yvar, estimate = b, se_cluster = se,
             ci_low = b - 1.96 * se, ci_high = b + 1.96 * se,
             n = nobs(m), row.names = NULL)
}

rd_results <- bind_rows(lapply(c("classize", "avgmath", "avgverb"), rd_reg))
rd_results
##    outcome   estimate se_cluster      ci_low   ci_high    n
## 1 classize -11.000749   1.182990 -13.3194101 -8.682088 1180
## 2  avgmath   3.428896   1.612143   0.2690963  6.588696 1180
## 3  avgverb   2.628131   1.497048  -0.3060818  5.562345 1180

2.6

Yes, 2.4.3 is consistent with that assumption, because there aren’t any significant discontinuities in the percent of disadvantaged students, the share of religious students, or the fraction of females. 2.4.4 has a little bit of a jump around 40, so it’s not as consistent with that assumption.

2.7

I would expect the class size to drop by 20 students, but it only drops by 11. This is probably because the “under 40” rule isn’t the most enforceable or practical in the real world, so there are a lot of exceptions. The cutoff is “fuzzy”.

2.8

b_class <- rd_results$estimate[rd_results$outcome == "classize"]
per_student <- rd_results %>%
  filter(outcome != "classize") %>%
  mutate(effect_per_student = estimate / b_class,
         effect_minus5      = -5 * effect_per_student)
per_student
##   outcome estimate se_cluster     ci_low  ci_high    n effect_per_student
## 1 avgmath 3.428896   1.612143  0.2690963 6.588696 1180         -0.3116966
## 2 avgverb 2.628131   1.497048 -0.3060818 5.562345 1180         -0.2389048
##   effect_minus5
## 1      1.558483
## 2      1.194524

Reducing class size to 35 would improve scores!

2.9 McCrary density test

b_class <- rd_results$estimate[rd_results$outcome == "classize"]
per_student <- rd_results %>%
  filter(outcome != "classize") %>%
  mutate(effect_per_student = estimate / b_class,      # Wald / fuzzy-RD ratio
         effect_minus5      = -5 * effect_per_student)  # reducing size by 5
per_student
##   outcome estimate se_cluster     ci_low  ci_high    n effect_per_student
## 1 avgmath 3.428896   1.612143  0.2690963 6.588696 1180         -0.3116966
## 2 avgverb 2.628131   1.497048 -0.3060818 5.562345 1180         -0.2389048
##   effect_minus5
## 1      1.558483
## 2      1.194524

The right side of enrollment 40 has a higher density than the left side of 40. It’s like two separate lines/curves, with a high density right above 40 students. This is not a perfect RD with perfect continuity.

2.10

No, because there is no cutoff/RD/threshold for us to observe, unlike with the 40 student cutoff.

2.11

Chetty et al (2011)’s effect is larger than mine in 2.8, and that might be because their design was a randomized experiment, and not a fuzzy RD. I was pretty close to the Fridriksson et al. (2013) paper, likely because they also used a fuzzy RD.

2.12

No, because in the Fredriksson paper, it shows long term effects on education and mid-life earnings too.

2.13

My hometown is in rural Texas, and the classes are already small. I don’t think this data is applicable to that setting. I’m not sure what would be plausible and appropriate for that setting. Maybe some kind of similar design with a RD/IV about class size if that is relevant for small towns as well… There would also need to be some kind of survey or tax record data for earnings information and educational attainment.

3.1

The IV can be which therapist the person happens to be assigned based on availability, and whether that therapist tends to say that their clients are mentally healthy or not.

3.2, 3.3, 3.4

SUTVA: The person’s treatment status into mental health court depends only on the insrument of the therapist assignment (not anyone else’s). This would be violated if the therapist’s mood from client impacts how they treat their next client. Independence: The therapist assignment is as good as random (this allows use to find the LATE with a good counterfactual). This could be violated by seletion into therapist shifts based on when a person was arrested (Saturday night arrests might look different from Tuesday morning arrests). Exclusion: any affect of therapist assignment on recidivism must be via the effect of therapist assignment into mental health court. This would be violated if the therapists notes could be seen by other people in the legal process like their parole officer, because then there is another channel to impact the outcome (behavioral bias based on what a parole officer saw in the note). Nonzero first stage: The therapist you happened to draw does impact your probability of being sent to mental health court. This is testable. This would be violated if a therapist’s identity doesn’t really impact meantal health court attendence. Monotonicity: All those affected by the instrument are affected in the same direction (you can’t be diagnosed with a mental illness and then not go to mental health court). Monotonicity would be violated if therapists disagree on which cases are bad enough to warrant meantal health court, which might be different than how lenient the therapist is.

3.5

First stage \[D_i = \pi_0 + \pi_1 Z_i + X_i'\pi_2 + \nu_i\]

Second stage \[Y_i = \beta_0 + \beta_1 \hat{D}_i + X_i'\beta_2 + \varepsilon_i\]

3.6

\[\beta_1 = E\left[Y_i(1) - Y_i(0) \mid \text{complier}\right]\]

3.7

We can test weak instruments using a first-stage F-stat. Reporting Anderson-Rubin confidence intervals is reccommended to support the F-stat. There is the Montiel Olea and Pflueger effective F statistic that accounts for heteroskedasticity. There is also the Kleibergen and Paap Wald test.