Chapter 8

Code Call-out 8.1: Quantile Effects with an Exogenous Treatment

In this code call-out we will work with results and data discussed in Firpo (2007), which require us to return once again to the LaLonde (1986) data we have discussed at length in Chapter 3. However, here rather than focussing on average treatment effects on the treated, we will focus on a range of quantile treatment effects (QTEs) and quantile treatment effects on the treated (QTTs).

Estimating QTEs with Experimental Interventions

We will begin by exploring QTEs when treatment assignment is random, and does not require the incorporation of any controls. In these cases, at least for the estimates themselves, it is sufficient to simply directly calculate quantiles among the treated and control units. Let’s first load these data, and ensure that we know which our outcome and treatment variables are:

library(haven)
library(dplyr)
library(ggplot2)
library(quantreg)

df <- read_dta("data/Dehejia_Wahba_2002.dta")

# Describe data structure
glimpse(df)

Rows: 16,437
Columns: 11
$ data_id   <chr> "Dehejia-Wahba Sample", "Dehejia-Wahba Sample", "Dehejia-Wah…
$ treat     <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ age       <dbl> 37, 22, 30, 27, 33, 22, 23, 32, 22, 33, 19, 21, 18, 27, 17, …
$ education <dbl> 11, 9, 12, 11, 8, 9, 12, 11, 16, 12, 9, 13, 8, 10, 7, 10, 13…
$ black     <dbl> 1, 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ hispanic  <dbl> 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ married   <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, …
$ nodegree  <dbl> 1, 1, 0, 1, 1, 1, 0, 1, 0, 0, 1, 0, 1, 1, 1, 1, 0, 1, 0, 0, …
$ re74      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ re75      <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, …
$ re78      <dbl> 9930.0459, 3595.8940, 24909.4492, 7506.1460, 289.7899, 4056.…

summary(df$re78)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0    5089   15962   14588   25565   60308 

table(df$treat)

    0     1 
16252   185 

# Group means
df %>% group_by(treat) %>%
  summarise(mean_re78 = mean(re78, na.rm = TRUE),
            n         = n())

# Keep Dehejia-Wahba experimental sample
df <- df %>% filter(data_id == "Dehejia-Wahba Sample")

Above we have kep the data corresponding to the “Dehejia-Wahba Sample” which is the original 445 treatment and control observations in the experimental implementation discussed in Dehejia and Wahba (1999). We can also exmaine the density of outcomes for both treatment and control units below:

# Kernel density plot
ggplot(df, aes(x = re78, colour = factor(treat))) +
  geom_density(linewidth = 1.2) +
  scale_colour_manual(values = c("red","blue"),
                      labels = c("Control","Treated")) +
  labs(x = "re78", colour = "Group") +
  theme_minimal()

Examining densities of treatment and control

With these experimental data, let’s calculate some QTE. In particular, let’s calculate the quantity \(\tau_{QTE(0.8)}\), which is the effect at quantile 80. We can do this simply by calculating the quantiles in each group:

q80_treated <- quantile(df$re78[df$treat == 1], 0.80)
q80_control <- quantile(df$re78[df$treat == 0], 0.80)
cat("QTE(0.8) =", q80_treated - q80_control, "\n")

QTE(0.8) = 2273.017 

Of course, we are not limited to doing this at the 80^th quantile. We can examine QTEs at multiple points of the distribution, below examining quantiles 25, 50 and 75 (and comparing these to ATEs themselves):

quantiles <- c(0.25, 0.5, 0.75)

for (q in quantiles) {
  qt <- quantile(df$re78[df$treat == 1], q, na.rm = TRUE)
  qc <- quantile(df$re78[df$treat == 0], q, na.rm = TRUE)
  cat(sprintf("QTE at quantile %.2f = %.2f\n", q, qt - qc))
}

QTE at quantile 0.25 = 485.23
QTE at quantile 0.50 = 1093.51
QTE at quantile 0.75 = 2354.58

mean_treated <- mean(df$re78[df$treat == 1], na.rm = TRUE)
mean_control <- mean(df$re78[df$treat == 0], na.rm = TRUE)
cat("Mean Treated =", mean_treated, "\n")

Mean Treated = 6349.144 

cat("Mean Control =", mean_control, "\n")

Mean Control = 4554.801 

cat("ATE =", mean_treated - mean_control, "\n")

ATE = 1794.342 

If we inspect these QTEs and compare them to the values reported in the Supplementary Table II of Firpo (2007), we can see that they are very similar, with some marginal differences at certain quantiles likely owing to differences in manners of estimating quantiles (there are a number of ways to estimate quantiles, including different way with dealing with ties, interpolations, and so forth). We can go also further, and inspect effects across the entire distribution, which we do below at each centile:

centiles  <- 1:99
FY1       <- quantile(df$re78[df$treat == 1], centiles / 100)
FY0       <- quantile(df$re78[df$treat == 0], centiles / 100)
diff_cdf  <- FY1 - FY0

cdf_df <- data.frame(quantile = centiles, FY1 = FY1, FY0 = FY0, diff = diff_cdf)

ggplot(cdf_df) +
  geom_line(aes(x = FY1, y = quantile, colour = "Treated"), linewidth = 1) +
  geom_line(aes(x = FY0, y = quantile, colour = "Control"),
            linewidth = 1, linetype = "dashed") +
  scale_colour_manual(values = c("Treated" = "blue", "Control" = "red")) +
  labs(x = "Earnings in 1978", y = "Quantile", colour = NULL) +
  theme_minimal()

Empirical CDFs of 1978 Earnings by Treatment Status

ggplot(cdf_df, aes(x = quantile, y = diff)) +
  geom_line(linewidth = 1) +
  labs(x = "Quantile", y = "QTE") +
  theme_minimal()

QTEs Across the Distribution

Typically, inference on quantile treatment effects will proceed using a bootstrap. However, we can also estimate these quantile treatment effects (along with standard errors) using quantile regression in this case with unconditional unconfoundedness. It is important to note, however, that we may see minor differences between quantile regression and quantile treatment effects as calculated “by hand” above depending on the way that percentiles are calculated. We can see this below where we estimate a simple quantile regression at the 80^th percentile:

df <- read_dta("data/Dehejia_Wahba_2002.dta")
df <- df %>% filter(data_id == "Dehejia-Wahba Sample")

qreg_80 <- rq(re78 ~ treat, data = df, tau = 0.8)

summary(qreg_80)

Call: rq(formula = re78 ~ treat, tau = 0.8, data = df)

tau: [1] 0.8

Coefficients:
            coefficients lower bd   upper bd  
(Intercept)  8551.5332    7352.9248 10104.0272
treat        2195.8164     333.5104  4240.7364

Here, while similar, this value of 2195 is different to the unconditional QTE of 2273 reported above. However, if we dig and explore this particular case, we can see that this simply owes to the fact that in our control group, there is an even number of units, and the 80^th percentile happens to fall between two units, and in these cases decisions on interpolation differ between commands. We can confirm to ourselves that there is indeed equivalence between these two settings provided that quantiles are smooth if we simply remove one observation so that percentiles in each group now fall precisely on a specific unit for the values below:

# Sort and drop one observation (with re78 = 0)
df_drop <- df %>% arrange(treat, re78) %>% slice(-1)

for (p in c(0.25, 0.50, 0.70)) {
  qreg_p <- rq(re78 ~ treat, data = df_drop, tau = p)
  cat(sprintf("\n--- Quantile %.2f ---\n", p))
  print(summary(qreg_p, se = "boot", R = 200)$coefficients)

  q1 <- quantile(df_drop$re78[df_drop$treat == 1], p)
  q0 <- quantile(df_drop$re78[df_drop$treat == 0], p)
  cat(sprintf("Direct: treated = %.4f, control = %.4f, effect = %.4f\n",
              q1, q0, q1 - q0))
}

--- Quantile 0.25 ---
               Value Std. Error  t value  Pr(>|t|)
(Intercept)   0.0000     0.0000      NaN       NaN
treat       485.2298   305.7044 1.587252 0.1131708
Direct: treated = 485.2298, control = 0.0000, effect = 485.2298

--- Quantile 0.50 ---
               Value Std. Error  t value     Pr(>|t|)
(Intercept) 3194.010   669.9197 4.767751 2.533056e-06
treat       1038.299   913.4606 1.136665 2.562941e-01
Direct: treated = 4232.3091, control = 3194.0100, effect = 1038.2991

--- Quantile 0.70 ---
               Value Std. Error   t value   Pr(>|t|)
(Intercept) 6378.720   622.6505 10.244463 0.00000000
treat       1795.188   933.3814  1.923317 0.05508234
Direct: treated = 8164.0695, control = 6368.9097, effect = 1795.1599

In this case, we see that both the QTEs and the means among treated and control units are equivalent to the quantile-regression based calculation.

Estimating QTEs and QTTs with non-Experimental Interventions

If, rather than believing that unconditional unconfoundedness hold, we instead believe that conditional unconfoundedness is the appropriate assumption, there are two potential ways forward. To look at this, we will use the same LaLonde (1986) data, but now with the experimental treatment sample, and the non-experimental subsample as potential control units. We will control for the same factors we have discussed in code call-outs in Chapter 3. Let’s begin by loading the requisite data, which is the observational (CPS) subsample, along with treated units, and generating a number of required covariates:

df <- read_dta("data/Dehejia_Wahba_2002.dta")
df <- df %>% filter(data_id == "CPS1" | treat == 1)

df <- df %>%
  mutate(
    age2    = age * age,
    age3    = age * age * age,
    educ2   = education * education,
    u74     = as.integer(re74 == 0),
    u75     = as.integer(re75 == 0),
    edure74 = education * re74
  )

xvars <- c("age","age2","age3","education","educ2","black","hispanic",
           "married","re74","re75","u74","u75","edure74")

Conditional QTEs

Let’s start by examining how we estimate QTEs in this observational setting. A first, simpler, case exists if we are happy to report conditional QTEs. We discuss the potential drawbacks of such QTEs in the book, but for now we can just note that if we are happy to report such conditional QTEs and invoke the required assumptions discussed in Section 8.2.2.1, we can do this very simply, by just estimating a quantile regression with covariates. We do this below (with an identical group of covariates to those used in code call out 3.1), reporting conditional quantile effects at the 25^th, 50^th, 75^th, and 80^th quantiles:

fml <- as.formula(paste("re78 ~ treat +", paste(xvars, collapse = " + ")))

# Simultaneous quantile regression at four quantiles (equivalent to Stata's sqreg)
taus <- c(0.25, 0.50, 0.75, 0.80)
sqreg_fit <- rq(fml, data = df, tau = taus)
summary(sqreg_fit, se = "boot", R = 100)

Call: rq(formula = fml, tau = taus, data = df)

tau: [1] 0.25

Coefficients:
            Value       Std. Error  t value     Pr(>|t|)   
(Intercept) 11316.04186  1317.44539     8.58938     0.00000
treat         336.96883   256.42240     1.31412     0.18883
age         -1095.42275   112.69150    -9.72054     0.00000
age2           29.45193     3.20214     9.19757     0.00000
age3           -0.26046     0.03056    -8.52305     0.00000
education      -6.59781    54.59238    -0.12086     0.90381
educ2           0.38214     1.99997     0.19108     0.84847
black        -288.18823   157.79592    -1.82634     0.06782
hispanic     -278.95136   214.99213    -1.29750     0.19448
married       150.56290   106.26067     1.41692     0.15653
re74            0.20193     0.03600     5.60964     0.00000
re75            0.72535     0.03079    23.55571     0.00000
u74          1143.29539   141.97372     8.05287     0.00000
u75           802.64365   145.02619     5.53447     0.00000
edure74         0.00845     0.00114     7.44059     0.00000

Call: rq(formula = fml, tau = taus, data = df)

tau: [1] 0.5

Coefficients:
            Value       Std. Error  t value     Pr(>|t|)   
(Intercept) 18840.38606   985.25418    19.12236     0.00000
treat         691.35408   374.99667     1.84363     0.06526
age         -1312.07757    74.71936   -17.56007     0.00000
age2           31.10551     1.87826    16.56080     0.00000
age3           -0.24300     0.01567   -15.50663     0.00000
education     125.67497    38.22216     3.28801     0.00101
educ2          -0.40390     0.80465    -0.50196     0.61570
black        -174.04308    86.24114    -2.01810     0.04360
hispanic      -21.65102    70.96513    -0.30509     0.76030
married        51.21724    43.79337     1.16952     0.24221
re74            0.33911     0.03040    11.15493     0.00000
re75            0.63622     0.02699    23.57590     0.00000
u74            23.90396   200.00398     0.11952     0.90487
u75         -2110.50757   240.95091    -8.75908     0.00000
edure74        -0.00445     0.00139    -3.19831     0.00139

Call: rq(formula = fml, tau = taus, data = df)

tau: [1] 0.75

Coefficients:
            Value       Std. Error  t value     Pr(>|t|)   
(Intercept) 15936.34282  1716.90809     9.28200     0.00000
treat        3210.62942  1296.04367     2.47725     0.01325
age          -735.41179   119.12064    -6.17367     0.00000
age2           14.49999     3.01117     4.81541     0.00000
age3           -0.09685     0.02489    -3.89062     0.00010
education     262.45127    92.22226     2.84586     0.00443
educ2          19.45432     4.02027     4.83905     0.00000
black        -352.27268   137.45953    -2.56274     0.01039
hispanic      132.27756   200.04516     0.66124     0.50847
married       350.25084    84.77779     4.13140     0.00004
re74            0.65315     0.04233    15.43113     0.00000
re75            0.35694     0.01910    18.68710     0.00000
u74         -1227.02233   364.49071    -3.36640     0.00076
u75         -4516.98513   451.64713   -10.00114     0.00000
edure74        -0.03182     0.00265   -11.99063     0.00000

Call: rq(formula = fml, tau = taus, data = df)

tau: [1] 0.8

Coefficients:
            Value       Std. Error  t value     Pr(>|t|)   
(Intercept) 12357.80542  2030.45609     6.08622     0.00000
treat        3558.44701  1114.57868     3.19264     0.00141
age          -440.53281   165.83821    -2.65640     0.00791
age2            7.92283     4.24992     1.86423     0.06231
age3           -0.04923     0.03552    -1.38594     0.16579
education     427.87446    72.90011     5.86933     0.00000
educ2          18.80605     3.70859     5.07094     0.00000
black        -468.69424   150.45074    -3.11527     0.00184
hispanic      221.27140   167.01351     1.32487     0.18523
married       391.43314   112.21584     3.48822     0.00049
re74            0.71395     0.03819    18.69222     0.00000
re75            0.28623     0.01377    20.79034     0.00000
u74         -1506.82977   377.85806    -3.98782     0.00007
u75         -4358.40902   445.23114    -9.78909     0.00000
edure74        -0.03907     0.00263   -14.86184     0.00000

The above covariate structure is slightly different to that described in the Supplementary Materials of Firpo (2007), but is useful in providing an identical comparison to previous covariate specifications. We can also examine this over the entire distribution of (conditional) quantiles. Below we loop through quantiles 20-95 (ie quantiles where we observed non-zero salaries in the experimental case). We do this as a simple loop, in each iteration saving the estimated quantile treatment effect.

qte_cond <- data.frame(
  quantile = numeric(),
  qte      = numeric(),
  qte_se   = numeric()
)

for (q in 20:95) {
  tau <- q / 100
  fit <- rq(fml, data = df, tau = tau)
  s   <- summary(fit, se = "boot", R = 100)$coefficients
  qte_cond <- rbind(qte_cond, data.frame(
    quantile = tau,
    qte      = s["treat", "Value"],
    qte_se   = s["treat", "Std. Error"]
  ))
}

Finally, we can plot the resulting distribution of quantile treatment effects. We do this below, and note that we see a reasonable correspondence with the experimentally estimated QTEs laid out previously, however with some important differences particularly at the upper end of the distribution.

qte_cond <- qte_cond %>%
  mutate(
    qte_upper = qte + qnorm(0.975) * qte_se,
    qte_lower = qte + qnorm(0.025) * qte_se
  )

ggplot(qte_cond, aes(x = quantile, y = qte)) +
  geom_line(linewidth = 0.8) +
  geom_line(aes(y = qte_upper), linetype = "dashed") +
  geom_line(aes(y = qte_lower), linetype = "dashed") +
  labs(x = "Quantile", y = "Earnings in 1978") +
  scale_x_continuous(breaks = seq(0.2, 1, 0.1),
                     labels = sprintf("%.1f", seq(0.2, 1, 0.1))) +
  theme_minimal()

Conditional Quantile Treatment Effects

Unconditional QTEs

While the previously implemented methods allow for the calculation of conditional quantile effects, we can use these tools to calculate unconditional QTEs. Below, we implement Firpo (2007)’s reweighted estimator “by hand”. You may note that while there is a canned routine which can be used to implement these methods available in Stata (the ivqte library), no such implementation exists for R. Nevertheless, the estimator is relatively straightforward to implement directly from the description in Firpo (2007).

The procedure works in two steps. First, we estimate a propensity score \(\hat{p}(X)\). Second, we use this propensity score to construct IPW weights which rebalance the observed outcome distributions to recover the unconditional counterfactual distributions \(F_{Y(1)}\) and \(F_{Y(0)}\). Specifically, following Firpo (2007), the weight for each observation is: \[ \hat{w}_i = \frac{D_i - \hat{p}(X_i)}{\hat{p}(X_i)(1-\hat{p}(X_i))} - \bar{w} \] where \(\bar{w}\) is the sample mean of the uncentred weights and \(P_c = E[D_i \hat{w}_i]\) is a normalising constant. These weights are used to build reweighted empirical CDFs for the treated and control potential outcome distributions by accumulating weighted running sums over the sorted outcome values. The unconditional quantile treatment effect at quantile \(\tau\) is then obtained by inverting each reweighted CDF at \(\tau\) and taking the difference, i.e. finding the \(\tau\)^th quantile of \(F_{Y(1)}\) and subtracting the \(\tau\)^th quantile of \(F_{Y(0)}\). Propensity scores that are very close to 0 or 1 receive extreme weights and are trimmed prior to estimation, with the default trimming threshold of 0.001 matching the default in Stata’s ivqte. If you refer to the Stata section of this code call-out, you can see the our implementation by hand is virtually identical, while the Stata implementation also provides calculations for the standard errors and confidence intervals.

# Estimate propensity score via logit
ps_fml  <- as.formula(paste("treat ~", paste(xvars, collapse = " + ")))
ps_fit  <- glm(ps_fml, data = df, family = binomial("logit"))
df$pscore <- predict(ps_fit, type = "response")

# Define Weighting function
firpo_qte <- function(data, tau, trim = 0.001) {
  d   <- data %>% filter(pscore >= trim & pscore <= 1 - trim)
  y   <- d$re78
  ps  <- d$pscore
  tr  <- d$treat
  n   <- length(y)

  w   <- (tr - ps) / (ps * (1 - ps))
  w   <- w - mean(w)
  Pc  <- mean(tr * w)

  ord   <- order(y)
  y_s   <- y[ord]
  tr_s  <- tr[ord]
  w_s   <- w[ord]

  temp1 <- cumsum(tr_s * w_s)       / Pc / n
  temp0 <- cumsum((tr_s - 1) * w_s) / Pc / n

  ys    <- sort(unique(y))
  dist1 <- sapply(ys, function(yv) temp1[max(which(y_s <= yv))])
  dist0 <- sapply(ys, function(yv) temp0[max(which(y_s <= yv))])

  q1 <- ys[max(1, sum(dist1 <= tau))]
  q0 <- ys[max(1, sum(dist0 <= tau))]
  q1 - q0
}

cat("Unconditional QTEs (Firpo 2007, default trimming):\n")

Unconditional QTEs (Firpo 2007, default trimming):

for (tau in c(0.25, 0.50, 0.75, 0.80)) {
  cat(sprintf("tau = %.2f: coef = %.4f\n", tau,
              firpo_qte(df, tau)))
}

tau = 0.25: coef = 3994.2946
tau = 0.50: coef = -1352.1201
tau = 0.75: coef = -3639.6406
tau = 0.80: coef = -2850.5400

cat("\nUnconditional QTEs (no trimming):\n")

Unconditional QTEs (no trimming):

for (tau in c(0.25, 0.50, 0.75, 0.80)) {
  cat(sprintf("tau = %.2f: coef = %.4f\n", tau,
              firpo_qte(df, tau, trim = 0)))
}

tau = 0.25: coef = -762.5073
tau = 0.50: coef = -10771.1519
tau = 0.75: coef = -15471.7910
tau = 0.80: coef = -13003.9912

# With more trimming
cat("\nUnconditional QTEs (trim at 0.05, 0.95):\n")

Unconditional QTEs (trim at 0.05, 0.95):

for (tau in c(0.25, 0.50, 0.75, 0.80)) {
  cat(sprintf("tau = %.2f: coef = %.4f\n", tau,
              firpo_qte(df, tau, trim = 0.05)))
}

tau = 0.25: coef = 626.5561
tau = 0.50: coef = 1774.8838
tau = 0.75: coef = 712.2637
tau = 0.80: coef = -53.7383

In this case, we see that the estimator performs quite poorly without trimming, suggesting that the conditional unconfoundedness assumption is likely unreasonable when considering all observations. However, when using more judicious trimming, we observe estimates which are at least broadly positive. Firpo (2007) documents results which are more broadly similar to those in the experimental sub-sample when using a much richer specification for the propensity score, pointing to the importance of appropriately modelling the propensity score.

Code Call-out 8.2: Extrapolating Regression Discontinuity Parameters Away from the Cut-off

Londoño-Vélez, Rodríguez, and Fabio Sánchez (2020) study the impact of receiving financial aid for post-secondary education on rates of enrollment among low-income individuals in Colombia. Specifically, they take advantage of program eligibility rules based on cut-off scores in a wealth index to isolate effects of financial aid eligibility. These cut-off scores in the wealth index imply that for individuals whose family wealth index is below a specific explicitly designed score, they are eligible to receive financial aid provided that they meet test score requirements. However, comparable individuals with scores in the wealth index even marginally above this cut-off, are not eligible to receive financial aid.

While this suggests a standard regression discontinuity design, one novelty of these wealth cut-off scores is that they are not fixed nation-wide, but rather vary by location. For individuals living in rural areas, individuals with a score below 40.75 are eligible for financial aid, while in large metropolitan areas, this score is 57.21 (refer to Londoño-Vélez, Rodríguez, and Fabio Sánchez (2020), page 201 for full details). This suggests that we can extrapolate findings away from specific cut-offs to consider the generalisability of any treatment effect local to specific cut-off scores, following Cattaneo et al. (2021).

Confirming Discontinuities in Treatment Eligibility

To see the broad context of the study, we will begin by confirming that there is a discontinuity in financial aid eligibility rates around the test score cut-off. We begin by opening the data from Londoño-Vélez, Rodríguez, and Fabio Sánchez (2020) and working with the sample they use in the paper. Specifically here, we will also impose the restriction eligible_saber11==1 which implies that all individuals in the sample meet educational criteria for financial aid, and so the only discontinuity which exists is that owing to the wealth eligibility criteria:

library(haven)
library(dplyr)

# Load data
data <- read_dta("data/LondonoVelez_et_al_2020.dta")

data <- data %>%
  filter(icfes_per == 20142,
         eligible_saber11 == 1)

Now let’s visualise the discontinuity in eligibility in the entire sample. Here we will work with a re-centred running variable which for each individual defines their distance to the area-specific eligigility threshold which applied to them. In order to set up a simple visualisation we will use an arbitrary definition setting bins from 50 points below the threshold up to 50 points above the threshold, in increments of 2 points. For a discussion of optimal bins in this setting refer to the discussion in Chapter 6 of the book. Within each bin, we will generate average scores as bin, and will plot a single point for each average score. We will then overlay a linear fit on either side of the cut-off using the original data.

library(ggplot2)

data <- data %>%
  filter(abs(running_sisben) < 50) %>%
  mutate(cut_int = floor(running_sisben / 2) * 2)

#Create bins
bin_means <- data %>%
  group_by(cut_int) %>%
  summarise(bin     = mean(beneficiary_spp, na.rm = TRUE),
            .groups = "drop") %>%
  mutate(midpoint = cut_int + 1)

ggplot(data, aes(x = running_sisben, y = beneficiary_spp)) +
  geom_smooth(data = filter(data, running_sisben <  0),
              method = "lm", se = TRUE,
              colour = "red", fill = "grey70", alpha = 0.5) +
  geom_smooth(data = filter(data, running_sisben >= 0),
              method = "lm", se = TRUE,
              colour = "red", fill = "grey70", alpha = 0.5) +
  geom_point(data = bin_means, aes(x = midpoint, y = bin),
             colour = "black", size = 3, shape = 1) +
  labs(x = "Distance to eligibility threshold",
       y = "Proportion beneficiaries") +
  scale_y_continuous(labels = scales::label_number(accuracy = 0.1)) +
  theme_minimal()

Above we see clear evidence of a sharp increase in eligibility when individuals fall just below the cut-off point. While no one with an above cut-off score is found to be eligible, this immediately jumps to around 60% eligibility among those with a below cut-off score. In what remains of the code call out below we will focus on the impact of falling below this cut-off, rather than the eligibility criteria itself. In effect, we wil consider a sharp design rather than a fuzzy design, though note that discussion in Cattaneo et al. (2021) points to how we could generalise this for a fuzzy design.

Visualising Multiple Treatment Cut-offs

Before we consider the process of extrapolating across cut-offs, let’s begin by confirming that we do indeed see multiple cut-offs owing to the differential treatment thresholds. We do this below, using the sisben_area variable, which takes 1 for large metropolitan areas which have a cut-off of 57.21 points, and 3 for rural areas with a cut-off of 40.75 points. There is actually also a third group (other urban areas) which takes a value of 2 and has a cut-off of 56.32 points, but because this is very close to group 1, we will only focus on rural and metropolitan areas.

We will generate two graphs, simply seeking to confirm that we see a sharp cut-off for each group at the point where, theoretically, such a cut-off should appear. We will do this with the precise score on the wealth index (sisben_score).

library(patchwork)

plot_cutoff <- function(df, area, cutoff, title_label) {
  d <- df %>%
    filter(sisben_area == area) %>%
    mutate(cut_int = floor(sisben_score / 2) * 2)

  bins <- d %>%
    group_by(cut_int) %>%
    summarise(bin = mean(beneficiary_spp, na.rm = TRUE), .groups = "drop")
  
  bins$midpoint <- bins$cut_int + 1

  ggplot(d, aes(x = sisben_score, y = beneficiary_spp)) +
    geom_smooth(data = filter(d, sisben_score <  cutoff),
                method = "lm", se = TRUE,
                colour = "red", fill = "grey70", alpha = 0.5) +
    geom_smooth(data = filter(d, sisben_score >= cutoff),
                method = "lm", se = TRUE,
                colour = "red", fill = "grey70", alpha = 0.5) +
    geom_point(data = bins, aes(x = midpoint, y = bin),
               colour = "black", size = 1) +
    geom_vline(xintercept = cutoff, colour = "red") +
    labs(title = title_label, y = "SPP Beneficiary", x = "sisben_score") +
    scale_y_continuous(labels = scales::label_number(accuracy = 0.01)) +
    theme_minimal() +
    theme(legend.position = "none")
}

plot_cutoff(data, area = 3, cutoff = 40.75, title_label = "Rural") +
plot_cutoff(data, area = 1, cutoff = 57.21, title_label = "Metropolitan")

These graphs above are really only used to show descriptive patters, as we have built bins in average of 2 points, and so the final bin on the right hand side will be slightly contaminated with above-threshold points, but we can correct this by generating different cut-points easily enough, and the graphs above make clear that there are clear discontinuities at the points which correspond to the specific wealth cut-off which binds for each group.

Let’s now consider some outcome of interest, and the extrapolation of treatment effects which we seek to achieve. In particular, let’s work with the variable spadies_any which indicates whether an individual studies any type of post-secondary education. And let’s visualise the mean outcome and a local polynomial fit of the outcomes for each group. Below we do this in a single graph. To do so we generate two variables which we can use to generate means (round1 and round3, which are centered on the cut-off), and which we arbitrarily set in terms of 2 point bins. We then use group_by(...) to generate mean scores in each of these groups simultaneously, before plotting these means as well as local polynomial fits on each side. Note that we generate group means without collapsing so that we can generate local polynomial fits based on the original microdata, rather than collapsed means.

library(KernSmooth)

data <- data %>%
  mutate(
    round1 = ifelse(sisben_area == 1, floor(running_sisben / 2), NA_integer_),
    round3 = ifelse(sisben_area == 3, floor(running_sisben / 2), NA_integer_)
  ) %>%
  group_by(round1, round3, sisben_area) %>%
  mutate(n = row_number(), bin_group = mean(spadies_any, na.rm = TRUE)) %>%
  ungroup()

# Helper function to generate local polynomial fit (local linear, Epanechnikov, bwidth=10)
lp_df <- function(df, bw = 10) {
  keep <- complete.cases(df$sisben_score, df$spadies_any)
  if (!any(keep)) return(tibble(x = numeric(0), y = numeric(0)))
  lp <- locpoly(df$sisben_score[keep], df$spadies_any[keep],
                bandwidth = bw, degree = 1, kernel = "epanechnikov",
                gridsize = 200, range.x = range(df$sisben_score[keep]))
  tibble(x = lp$x, y = lp$y)
}

fit_low_left   <- lp_df(filter(data, sisben_area == 3, running_sisben < 0))
fit_low_right  <- lp_df(filter(data, sisben_area == 3, running_sisben >= 0))
fit_high_left  <- lp_df(filter(data, sisben_area == 1, running_sisben < 0))
fit_high_right <- lp_df(filter(data, sisben_area == 1, running_sisben >= 0))

ggplot() +
  geom_point(data = filter(data, n == 1, sisben_area == 1),
             aes(x = sisben_score, y = bin_group, shape = "High cut-off", color = "High cut-off"),
             size = 3) +
  geom_point(data = filter(data, n == 1, sisben_area == 3),
             aes(x = sisben_score, y = bin_group, shape = "Low cut-off", color = "Low cut-off"),
             size = 3) +
  geom_line(data = fit_low_left,  aes(x, y, linetype = "Fit (low)",  color = "Low cut-off"), linewidth = 0.8) +
  geom_line(data = fit_low_right, aes(x, y, linetype = "Fit (low)",  color = "Low cut-off"), linewidth = 0.8) +
  geom_line(data = fit_high_left, aes(x, y, linetype = "Fit (high)", color = "High cut-off"), linewidth = 0.8) +
  geom_line(data = fit_high_right,aes(x, y, linetype = "Fit (high)", color = "High cut-off"), linewidth = 1.2) +
  scale_shape_manual(name = NULL, values = c("High cut-off" = 1, "Low cut-off" = 0)) + # 1 círculo hueco, 0 cuadrado hueco
  scale_color_manual(name = NULL, values = c("High cut-off" = "navy", "Low cut-off" = "red")) +
  scale_linetype_manual(name = NULL, values = c("Fit (low)" = "dashed", "Fit (high)" = "longdash")) +
  labs(x = "SISBEN wealth index", y = "Studying any tertiary education") +
  theme_minimal() +
  theme(legend.position = "bottom")

In the above plot we can quite easily see the idea of what we wish to do when extrapolating treatment effects away from the cut-off. Specifically, we wish to consider the first cut-off, here at 40.75 points. We wish to calculate the treatment effect at this cut-off by comparing outcomes among exposed units just at the left to those just at the right in the spirit of an RDD. And then we wish to consider what the treatment effect would look like at specific points above 40.75 if we use the trend among units exposed at a higher cut-off point (those with blue circles) to extrapolate means in the below-cut-off group with red squares, before finally comparing extrapolated means with actual observed rates among those with red squares to the right of the treatment cut-off.

Extrapolating Treatment Effects Away from the Cut-off

In order to do such an extrapolation, we require a way to estimate local polynomial fits at various points (including at end-points just before treatment cut-offs), as well as the variance of these local polynomial estimates. Fortunately, in relation to the work of Cattaneo et al. (2021), the authors developed software for such local polynomial fits, incorporating elements such as robust bias correction. This is avaialable (in Stata and R) as nprobust and we will work with this package below to estimate the required quantities.

Let’s begin by imagining that we wish to extrapolate treatment effects from the true cut-off of 40.75 up a higher point on the SISBEN wealth index (50 points). To do so, we need four quantities. Firstly, we need to calculate the end point of the low-cut-off group precisely at 40.75 points (ie the point just before the discontinuity kicks in). Cattaneo et al. (2021) call this first quantity \(\mu_{0,\ell}(\ell)\) Secondly, we need to calculate the mean values in the high-cut-off group at both 40.75 and 50 points, which allows us to calcualte any trend over this range. Cattaneo et al. (2021) refer to these as \(\mu_{0,h}(\ell)\) and \(\mu_{0,h}(\bar{x})\) respectively. And finally, we wish to calculate the mean among low-cut-off group outcomes at 50 points, which Cattaneo et al. (2021) refer to as \(\mu_{1,\ell}(\bar{x})\). Once we have these points in hand, as discussed in Section 8.4.2.2 of the book, we can simply extrapolate our treatment effect as: \[ \widehat\tau_\ell(\bar{x})=\widehat\mu_{1,\ell}(\bar{x}) -[\widehat\mu_{0,h}(\bar{x})+\widehat\mu_{0,\ell}(\ell)-\widehat\mu_{0,h}(\ell)]. \]

Let’s do this below, using `lprobust’ to generate key quantities. We will also take the variance at each point, allowing us to calculate the standard error of the extrapolated treatment effect as the square root of the total variance.

library(nprobust)

# Low cut-off (Rural) - left of 40.75
c0 <- 40.75
res0 <- lprobust(y = data$spadies_any[data$sisben_area == 3 & data$sisben_score < 40.75],
                 x = data$sisben_score[data$sisben_area == 3 & data$sisben_score < 40.75],
                 eval = c0)
mu_0_l_l <- res0$Estimate[1, 5]
v_0_l_l  <- (res0$Estimate[1, 8])^2

# Low cut-off (Rural) - right of 40.75
c1 <- 50
res1 <- lprobust(y = data$spadies_any[data$sisben_area == 3 & data$sisben_score >= 40.75],
                 x = data$sisben_score[data$sisben_area == 3 & data$sisben_score >= 40.75],
                 eval = c1)
mu_1_l_x <- res1$Estimate[1, 5]
v_1_l_x  <- (res1$Estimate[1, 8])^2

# High cut-off (Metropolitan) - left of 57.21
c2 <- c(40.75, 50)
res2 <- lprobust(y = data$spadies_any[data$sisben_area == 1 & data$sisben_score < 57.21],
                 x = data$sisben_score[data$sisben_area == 1 & data$sisben_score < 57.21],
                 eval = c2,
                 covgrid = TRUE, bwselect = "mse-dpi")
mu_0_h_l <- res2$Estimate[1, 5]
mu_0_h_x <- res2$Estimate[2, 5]
v_0_h_l  <- (res2$Estimate[1, 8])^2
v_0_h_x  <- (res2$Estimate[2, 8])^2
cov      <- res2$cov.rb[2, 1]

# Effect and variance
effect   <- mu_1_l_x - (mu_0_h_x + mu_0_l_l - mu_0_h_l)
variance <- v_0_l_l + v_1_l_x + v_0_h_l + v_0_h_x - 2 * cov

cat("Effect at 50 is:", effect, "\n")

Effect at 50 is: -0.261186 

cat("Variance at 50 is:", variance, "\n")

Variance at 50 is: 0.01230848 

While this is a single extrapolation, we can of course extrapolate more widely up to any point below the second cut-off, at which point no untreated units remain. Below we conduct a similar process now extrapolating across a range of values. This simply replicates the code above, but incorporates a loop for various values \(\bar{x}\) used to extrapolate.

# Grid of evaluation points
grid_vals <- seq(41, 57.2, by = 0.8)

results <- data.frame(
  runvar = numeric(),
  estimate = numeric(),
  std_error = numeric()
)

for (num in grid_vals) {
  
  # Left of 40.75 (rural)
  res1 <- lprobust(y = data$spadies_any[data$sisben_area == 3 & data$sisben_score < 40.75],
                   x = data$sisben_score[data$sisben_area == 3 & data$sisben_score < 40.75],
                   eval = 40.75)
  b1 <- res1$Estimate[1, 5]
  v1 <- (res1$Estimate[1, 8])^2
  
  # Right of 40.75 (rural)
  res2 <- lprobust(y = data$spadies_any[data$sisben_area == 3 & data$sisben_score >= 40.75],
                   x = data$sisben_score[data$sisben_area == 3 & data$sisben_score >= 40.75],
                   eval = num)
  b2 <- res2$Estimate[1, 5]
  v2 <- (res2$Estimate[1, 8])^2
  
  # High cut-off (metro) left and right
  res3 <- lprobust(y = data$spadies_any[data$sisben_area == 1 & data$sisben_score < 57.21],
                   x = data$sisben_score[data$sisben_area == 1 & data$sisben_score < 57.21],
                   eval = c(num, num),
                   covgrid = TRUE, bwselect = "mse-dpi")
  b3 <- res3$Estimate[1, 5]
  b4 <- res3$Estimate[2, 5]
  v3 <- (res3$Estimate[1, 8])^2
  v4 <- (res3$Estimate[2, 8])^2
  cov <- res3$cov.rb[2, 1]
  
  effect   <- (b2 - b4) - (b1 - b3)
  variance <- v1 + v2 + v3 + v4 - 2 * cov
  
  results <- rbind(results, data.frame(
    runvar = num,
    estimate = effect,
    std_error = sqrt(variance)
  ))
}

# Confidence intervals
results <- results %>%
  mutate(
    LB = estimate + qnorm(0.025) * std_error,
    UB = estimate + qnorm(0.975) * std_error
  )
print(results)

         runvar   estimate std_error         LB           UB
tau.us     41.0 -0.2162904 0.1335275 -0.4779996  0.045418712
tau.us1    41.8 -0.2343613 0.1179481 -0.4655353 -0.003187306
tau.us2    42.6 -0.2303900 0.1141077 -0.4540370 -0.006743101
tau.us3    43.4 -0.2320759 0.1119424 -0.4514791 -0.012672824
tau.us4    44.2 -0.2346076 0.1123994 -0.4549064 -0.014308854
tau.us5    45.0 -0.2367794 0.1123826 -0.4570453 -0.016513593
tau.us6    45.8 -0.2388162 0.1109337 -0.4562423 -0.021390194
tau.us7    46.6 -0.2409769 0.1107751 -0.4580921 -0.023861622
tau.us8    47.4 -0.2420735 0.1105644 -0.4587758 -0.025371260
tau.us9    48.2 -0.2429924 0.1100148 -0.4586173 -0.027367381
tau.us10   49.0 -0.2434295 0.1101239 -0.4592685 -0.027590576
tau.us11   49.8 -0.2429314 0.1103221 -0.4591588 -0.026704060
tau.us12   50.6 -0.2419040 0.1105902 -0.4586569 -0.025151092
tau.us13   51.4 -0.2414686 0.1110274 -0.4590784 -0.023858830
tau.us14   52.2 -0.2407317 0.1115256 -0.4593180 -0.022145473
tau.us15   53.0 -0.2390456 0.1122224 -0.4589976 -0.019093663
tau.us16   53.8 -0.2364377 0.1130101 -0.4579335 -0.014941881
tau.us17   54.6 -0.2345768 0.1134266 -0.4568889 -0.012264772
tau.us18   55.4 -0.2334605 0.1132280 -0.4553834 -0.011537695
tau.us19   56.2 -0.2303931 0.1132518 -0.4523625 -0.008423752
tau.us20   57.0 -0.2266534 0.1133982 -0.4489099 -0.004396917

# Plot
ggplot(results, aes(x = runvar, y = estimate)) +
  geom_ribbon(aes(ymin = LB, ymax = UB), fill = "grey50", alpha = 0.3) +
  geom_point(color = "red") +
  labs(y = "Extrapolated ATT", x = "Running variable") +
  scale_y_continuous(labels = scales::label_number(format = "%03.1f")) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "black") +
  theme_minimal() +
  theme(legend.position = "bottom")

Figure 1: Non-parametric fits with multiple cut-offs

We can see what the extrapolated effects look like across a range of values, presenting these estimates along with 95% confidence intervals. In this specific case we note that extrapolated effects are quite flat, which makes sense given that both relevant non-parametric fits (shown in Figure 1) in the relevant range between 40 and 57 are quite flat.

Code Call-out 8.3: Marginal Treatment Effects

In this code call-out we consider a number of elements related to the estimation of marginal treatment effects (MTEs). To do so, we work with data from Carneiro, Lokshin, and Umapathi (2017), who estimate the returns to upper secondary schooling in Indonesia and examine how those returns vary across individuals with different likelihoods of enrolling. This setting is well-suited to consider the MTE framework. If individuals sort into schooling partly based on their own anticipated returns, then the returns among those who enrol only when schooling is made more accessible (the compliers of any given instrument) may differ substantially from returns among those who would always enrol regardless. The MTE surface we consider here traces out exactly this variation.

Following Carneiro, Lokshin, and Umapathi (2017), data from the 2000 wave of the Indonesia Family Life Survey (IFLS) is used, restricted to men aged 25-60 who are employed and report non-missing wages and schooling, along with non-missing information on an instrument, key covariates used for calculating a propensity score. This yields a sample of 2,608 working-age males. We load these data below and subset to observations with non missing observations, examining summary statistics for log earnings, the key outcome in the estimation sample.

library(tidyverse)
library(haven)
library(caret)
library(dplyr)
library(ggplot2)
library(patchwork)


data <- read_dta("data/Carneiro_et_al_2017.dta")

X <- c("age", "age2",
  "r_protest", "r_cathol", "r_other",
  "elem_f", "jsec_f", "edumiss_f",
  "elem_m", "jsec_m", "edumiss_m",
  "rural", "kmsd",

  # provinces
  "prov_NSUM", "prov_WSUM", "prov_SSUM", "prov_LAMP",
  "prov_JAKA", "prov_CJAV", "prov_YOGI", "prov_EJAV",
  "prov_BALI", "prov_WNUSA", "prov_SKALI", "prov_SSUL"
)

INT_vars <- grep("^INT_", names(data), value = TRUE)

vars_all <- c("learnhr00", "dschool", "kmsmp", X, INT_vars)

data_sample <- data %>%
  mutate(
    touse = if_else(
      if_all(all_of(vars_all), ~ !is.na(.)),
      1, 0
    )
  ) %>%
  filter(touse == 1)

In this context, our “treatment” variable of interest is the binary measure of whether an individual attended upper secondary school (dschool). The authors also propose an instrument (kmsmp), which is a measure of the distance from each individual’s community head’s office to the nearest secondary school. In this call out we do not discuss the assumptions related to this IV, but rather follow the authors in using this, given the importance of instruments in this setting to trace out the MTE over a meaningful range of the propensity score.

Propensity Score Estimation and Common Support

Given the central nature of the propensity score in the MTE framework, we begin with its estimation. We estimate the propensity score using a logit model where the binary treatment indicator dschool is regressed on the full set of covariates (which will be included in models below), along with the instrument. Specifically, the specification includes the distance to the nearest secondary school (kmsmp), all interaction terms of covariates with the instrument (variables starting with INT_), and the baseline set of individual and household characteristics contained in the vector X defined above (age and age squared, parental education indicators, religious affiliation, rural status, distance to the closest health post, and province fixed effects). Below we estimate the logit, and then generate each indiviudal’s predicted probability of schooling, ie the propensity score.

form_ps <- as.formula(
  paste("dschool ~", paste(c("kmsmp", INT_vars, X), collapse = " + "))
)

logit_ps <- glm(form_ps, data = data_sample, family = binomial(link = "logit"))

data_sample$ps_manual <- predict(logit_ps, type = "response")  

Given that MTEs can only be identified over the support of the propensity score, hereafter \(P\), it is important to understand its distribution, and in particular the overlap across individuals who did and did not attend upper secondary school.

psvar <- "ps_manual"

p0 <- data_sample %>% 
  filter(dschool == 0, !is.na(.data[[psvar]])) %>%
  ggplot(aes(x = .data[[psvar]], y = after_stat(..count..) / sum(..count..))) +
  geom_histogram(bins = 30, fill = "skyblue3", color = "black") +
  labs(x = "less than upper secondary", y = "") +
  scale_y_continuous(limits = c(0, 0.10)) +
  scale_x_continuous(breaks = seq(0, 1, 0.2))

p1 <- data_sample %>% 
  filter(dschool == 1, !is.na(.data[[psvar]])) %>%
  ggplot(aes(x = .data[[psvar]], y = after_stat(..count..) / sum(..count..))) +
  geom_histogram(bins = 30, fill = "skyblue3", color = "black") +
  labs(x = "upper secondary", y = "") +
  scale_y_continuous(limits = c(0, 0.10)) +
  scale_x_continuous(breaks = seq(0, 1, 0.2))

p0 | p1

Figure 2: Propensity score by treatment status

To further assess the support conditions, it is useful to examine not only the overall distribution of the propensity score, but also how this distribution varies over the support of all covariates \(X\). This relates to a point raised by Carneiro, Heckman, and Vytlacil (2011) whote note:

“If all we are willing to assume is that (\(U_0\), \(U_1\),V) is independent of \(Z\) given \(X\), then it is only possible to estimate the MTE over the support of \(P\) conditional on \(X\).” (Carneiro, Heckman, and Vytlacil (2011), p. 2768)

To consider the relevance of this, we can follow Carneiro, Heckman, and Vytlacil (2011) in plotting the marginal distribution of \(P|X\). However, given that \(X\) is multidimensional, this is considered over an index \(X(\delta_1 - \delta_0)\). To do this, we estimate separate outcome equations for individuals with and without upper secondary schooling, obtain the coefficient vectors \(\hat{\delta}_1\) and \(\hat{\delta}_0\), and take their difference. Multiplying this difference by each individual’s covariates yields the index \(X(\hat{\delta}_1 - \hat{\delta}_0)\), which summarizes the observable component of the return to schooling. We generate this index below:

form_outcome <- as.formula(
  paste("learnhr00", "~", paste(X, collapse = " + "))
)

fit1 <- lm(form_outcome, data = data_sample, subset = (dschool == 1))
fit0 <- lm(form_outcome, data = data_sample, subset = (dschool == 0))

delta1_hat <- coef(fit1)
delta0_hat <- coef(fit0)

delta_diff <- delta1_hat - delta0_hat

X_mat <- model.matrix(form_outcome, data = data_sample)

data_sample$index_X <- as.numeric(X_mat %*% delta_diff)

With this index, we estimate a conditional density \(f(P \mid X)\), which is the density of the propensity score \(P\) at each point of the index. We do this nonparametrically using kde2d to estimate a joint density \(f(X|P)\). The resulting density values are then normalised within each value of \(X\) (each row of the resulting matrix) so that they sum to one, giving an estimate of \(f(P \mid X = x)\) for each bin’s representative value of \(x\). In essence, this simply builds densities locally within small sections of the data. This yields a dataset of triplets \((x, P, \hat{f}(P \mid x))\) across the joint support of \(X\) and \(P\), which can be visualised as a surface showing how the distribution of the propensity scores shifts as the covariate index changes.

library(MASS)    
library(plotly)  

kde <- kde2d(
  x = data_sample$index_X ,
  y = data_sample$ps_manual,
  n = 50,              
  lims = c(range(data_sample$index_X), 0, 1)  
)

z_joint <- kde$z   
row_sums <- rowSums(z_joint)
row_sums[row_sums == 0] <- NA        
z_cond <- sweep(z_joint, 1, row_sums, FUN = "/")

Once we have this triplet of points of the X grid, the propensity score, and the density, we can plot these in a 3-d surface. We do this below using the plotly library, allowing us to observe the density interactively.

library(plotly)

fig <- plot_ly(
  x = kde$x,              
  y = kde$y,             
  z = ~t(z_cond),         #
  type = "surface"
) %>%
  layout(
    title = "Support of P conditional on X index",
    scene = list(
      xaxis = list(title = "X(δ1 - δ0) index"),
      yaxis = list(title = "P (propensity score)"),
      zaxis = list(title = "f(P | X)")
    )
  )

fig

Inspecting the resulting surface, we can see that while there is a positive density at many points of support of the observable index \(X\), there are also areas with essentially no covarage, namely areas with quite high values of the index and high propensity scores. Given this lack of full support over all values of \(X\), typically assumptions are invoked such that all we require is common support of the propensity score across treatment and untreated units. Assumptions such as additive separability allow for this, which implies assuming \(E(U_D|V,X)=E(U_D|V)\), or that the slope of the MTE is independent of \(X\).

If invoking this assumption, all we need to consider is the common support of \(P\) across treatment regimes. In Figure Figure 2 we observe quite broad common support, however below we generate a variable which indicates whether observations are in the region of common support considering maximum and minimum propensity scores in each group. We will limit our analysis by removing the relatively small subset of observations (around 1%) for which there is no overlap.

CS_min <- max(
  min(data_sample$ps_manual[data_sample$dschool == 1]),
  min(data_sample$ps_manual[data_sample$dschool == 0])
)

CS_max <- min(
  max(data_sample$ps_manual[data_sample$dschool == 1]),
  max(data_sample$ps_manual[data_sample$dschool == 0])
)

data_sample$CS_dummy <- 
  ifelse(data_sample$ps_manual >= CS_min & data_sample$ps_manual <= CS_max,  1, 0)

table(data_sample$CS_dummy)

   0    1 
  32 2576 

print(prop.table(table(data_sample$CS_dummy)))

         0          1 
0.01226994 0.98773006 

data_sample <- data_sample %>% filter(CS_dummy == 1)

Estimating Marginal Treatment Effects

A Parametric Approach

There are multiple ways in which we can implement MTEs, and at times these can be quite computationally challenging. A simple and very illustrative way to estimate MTEs is through a parametric approach. In essence, all this requires is for us to model \(E[Y \mid P]\) as a flexible polynomial or spline (in terms of \(P\)), along with all relevant covariates. Then MTE is then the derivative of \(E[Y \mid P]\) with respect to \(P\) across all points of support of the propensity score. Consider the below “manual” implementation of such a parametric approach. Here we include the propensity score in a linear way along with 3 higher polynomial terms. We also include all controls both in levels, as well as interacted with the propensity score. This latter term allows us to consider whether returns to individual characteristics themselves depend on the likelihood of being treated.

# Generate interactions of covariates with propensity score
for (v in X) {
  data_sample[[paste0("PX_", v)]] <- data_sample$ps_manual * data_sample[[v]]
}
PX_vars <- paste0("PX_", X)

# Quartic polynomial in ps_manual plus covariate interactions
fml_mte <- as.formula(paste(
  "learnhr00 ~",
  paste(X, collapse = " + "), "+",
  paste(PX_vars, collapse = " + "), "+",
  "ps_manual + I(ps_manual^2) + I(ps_manual^3) + I(ps_manual^4)"
))

Y_ps <- lm(fml_mte, data = data_sample)
summary(Y_ps)

Call:
lm(formula = fml_mte, data = data_sample)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.7273 -0.5270 -0.0200  0.5163  5.0402 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)      6.255668   0.674591   9.273   <2e-16 ***
age              0.034079   0.033921   1.005   0.3152    
age2            -0.038681   0.041483  -0.932   0.3512    
r_protest        0.282987   0.295971   0.956   0.3391    
r_cathol        -0.755055   0.617426  -1.223   0.2215    
r_other          0.365061   0.334774   1.090   0.2756    
elem_f           0.023521   0.157371   0.149   0.8812    
jsec_f          -0.193658   0.544955  -0.355   0.7223    
edumiss_f        0.162263   0.265804   0.610   0.5416    
elem_m          -0.259223   0.127891  -2.027   0.0428 *  
jsec_m          -1.863968   0.771796  -2.415   0.0158 *  
edumiss_m       -0.236282   0.139476  -1.694   0.0904 .  
rural            0.234101   0.127567   1.835   0.0666 .  
kmsd            -0.001982   0.027458  -0.072   0.9425    
prov_NSUM        0.398837   0.170959   2.333   0.0197 *  
prov_WSUM        0.422549   0.175028   2.414   0.0158 *  
prov_SSUM        0.379108   0.241494   1.570   0.1166    
prov_LAMP        0.073261   0.238319   0.307   0.7586    
prov_JAKA       -0.219936   0.173003  -1.271   0.2037    
prov_CJAV        0.156518   0.120944   1.294   0.1957    
prov_YOGI       -0.038891   0.234756  -0.166   0.8684    
prov_EJAV       -0.056810   0.115098  -0.494   0.6216    
prov_BALI       -0.858596   0.397076  -2.162   0.0307 *  
prov_WNUSA      -0.420611   0.259486  -1.621   0.1052    
prov_SKALI      -0.158946   0.357194  -0.445   0.6564    
prov_SSUL       -0.163502   0.215845  -0.757   0.4488    
PX_age          -0.030083   0.077171  -0.390   0.6967    
PX_age2          0.093111   0.096719   0.963   0.3358    
PX_r_protest    -0.041885   0.497861  -0.084   0.9330    
PX_r_cathol      1.341191   0.887600   1.511   0.1309    
PX_r_other      -0.448937   0.576546  -0.779   0.4362    
PX_elem_f        0.071818   0.499877   0.144   0.8858    
PX_jsec_f        0.679267   1.031570   0.658   0.5103    
PX_edumiss_f    -0.722018   0.824206  -0.876   0.3811    
PX_elem_m        0.503776   0.319951   1.575   0.1155    
PX_jsec_m        2.997030   1.172905   2.555   0.0107 *  
PX_edumiss_m     0.090448   0.371373   0.244   0.8076    
PX_rural        -0.222445   0.293097  -0.759   0.4480    
PX_kmsd         -0.008282   0.071198  -0.116   0.9074    
PX_prov_NSUM    -0.860309   0.366704  -2.346   0.0191 *  
PX_prov_WSUM    -0.393342   0.416727  -0.944   0.3453    
PX_prov_SSUM    -0.348304   0.440687  -0.790   0.4294    
PX_prov_LAMP    -0.399451   0.667934  -0.598   0.5499    
PX_prov_JAKA     0.403414   0.307886   1.310   0.1902    
PX_prov_CJAV    -0.524205   0.325716  -1.609   0.1077    
PX_prov_YOGI    -0.297723   0.444189  -0.670   0.5028    
PX_prov_EJAV     0.113120   0.271512   0.417   0.6770    
PX_prov_BALI     1.358337   0.728812   1.864   0.0625 .  
PX_prov_WNUSA    0.534949   0.565298   0.946   0.3441    
PX_prov_SKALI    1.026031   0.642340   1.597   0.1103    
PX_prov_SSUL     0.357926   0.462994   0.773   0.4396    
ps_manual        5.321029   3.211007   1.657   0.0976 .  
I(ps_manual^2) -16.769728  11.208084  -1.496   0.1347    
I(ps_manual^3)  24.698999  15.951148   1.548   0.1216    
I(ps_manual^4) -14.146071   8.560011  -1.653   0.0985 .  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8785 on 2521 degrees of freedom
Multiple R-squared:  0.1216,    Adjusted R-squared:  0.1028 
F-statistic: 6.461 on 54 and 2521 DF,  p-value: < 2.2e-16

Once we have this parametric implementation, we simply can calculate a marginal treatment effect by exploring how our outcome \(E[Y|X]\) varies with the propensity score—or the resistance to treatment. This quantity is the marginal treatment effect, defined as in (8.42) of the book. Fortunately this is relatively easily calculated at a range of propensity scores using Stata’s margins command. We do this below, calculating the marginal effect of a change in the propensity score, saving the resulting effects, and their default standard errors. Finally, we can plot these marginal treatment effects as they vary with \(P(Z)\).

library(marginaleffects)

plot_slopes(Y_ps, variables = "ps_manual",
            condition = list(ps_manual = seq(0, 1, by = 0.01))) +
  labs(x = "Propensity score, P(Z, X)", y = "MTE") +
  theme_minimal()

Figure 3: Marginal Treatment Effects Under a Polynomial Specification

Figure 3 plots the MTE as a function of the propensity score. The MTE measures the marginal return to upper secondary schooling for individuals with a given probability of enrolling. The figure reveals individuals with a higher propensity to attend upper secondary school tend to experience higher marginal returns, while individuals with a lower propensity exhibit lower returns. While our approach here is quite manual it allows us to see the ideas behind parametric approaches to MTE estimation. Fully-fledged implementations of this which exist as formal routines in other language (namely Stata’s mtefe approach we discuss in the Stata section of this code call-out), produce very similar point estimates.

Semi-parametric (Local IV) methods

Policy Relevant Treatment Effects

As discussed in the Book, given the availability of estimated marginal treatment effects, we can use these as building blocks to estimate a large number of quantities of interest. To see a rough idea of this we can consider what our ATEs, ATTs, ATUs as well as a PRTE under a specific alternative policy might look like. We will do this manually to have an idea of the mechanics, though note that more formal ways about how to do this with weighting are used in practice; see Andresen (2018) for a computational discussion. To gain a rough idea of how we can use these MTEs to calculate marginal PRTEs, we essentially can consider the effect of policies which shift certain individuals into treatment. To build intuition for this, consider two stylised policies applied to currently untreated individuals. The first targets those with low propensity scores (\(P \in [0.05, 0.20]\)); i.e. individuals who are relatively willing to enrol (low propensity score, or low aversion to treatment) but currently not enrolled. These are the “easy shifters”: a modest policy change would be enough to tip them into treatment. The second targets individuals with high propensity scores among the untreated (\(P \in [0.55, 0.70]\)) these are “harder shifters”, or more resistant individuals who would only enrol under a strong intervention.

Based on the MTEs we have estimated above, we can consider what this would imply for the individuals in this setting. To illustrate this, we will use the MTEs we have previously calculated based on the parametric approach. If you refer above, you will remember that we generated a file called mte_grid with the marginal treatment effect for each binned propensity score from 0 to 1 in increments of 0.01. Below, we will merge these marginal treatment effects back into the original estimated propensity scores, allowing us to observe, for each individual (and hence propensity score), our estimate of their MTE. Of course, this will be a simple approximation as we are using a coarse grid of propensity scores.

# Build MTE grid from parametric estimates
b <- coef(Y_ps)
ps_grid_vals <- seq(0, 1, by = 0.01)
mte_grid_vals <- (b["ps_manual"]
                  + 2 * b["I(ps_manual^2)"] * ps_grid_vals
                  + 3 * b["I(ps_manual^3)"] * ps_grid_vals^2
                  + 4 * b["I(ps_manual^4)"] * ps_grid_vals^3)

mte_grid_df <- data.frame(
  at_idx  = seq_along(ps_grid_vals),
  ps_grid = ps_grid_vals,
  margin  = mte_grid_vals
)

data_mte <- data_sample %>%
  mutate(at_idx = round(ps_manual * 100) + 1) %>%
  left_join(mte_grid_df, by = "at_idx")

We can now consider the two movements of individuals into treatment we discussed above.

# Low resistance untreated
data_mte %>% filter(ps_manual >= 0.05, ps_manual < 0.20, dschool == 0) %>%
  summarise(mean_mte = mean(margin), n = n())

# High resistance untreated
data_mte %>% filter(ps_manual >= 0.55, ps_manual < 0.70, dschool == 0) %>%
  summarise(mean_mte = mean(margin), n = n())

We can see that in these two cases, the mean treatment effect is very different, with very high positive returns in the first case, and negative returns in the latter case. The contrast between these two quantities illustrates on of the central message of the MTE framework: the return to a policy depends not just on whether it shifts people into treatment, but on who it shifts. Given positive selection on gains such as those estimated here, policies reaching more willing individuals, i.e. those with low resistance, yield higher average returns than those requiring a more intensive intervention.

However, less abstractly, we can consider specific policy interventions, and how they would map into treatment effects. The idea of a PRTE is to consider an alternative policy \(P(Z')\). In this case, imagine an alternative policy in which distances to schools (the instrument considered above) are reduced, presumably via some sort of school construction program. Using the ideas of PRTEs, we can ask what such a movement in terms of the instrument implies, and for which individuals will such a movement be sufficient to shift them into treatment. Concretely, imagine a policy which reduced distance to secondary schooling for each individual by 0.1km. Because the instrument enters the propensity score via the estimated logit, we can translate this distance reduction directly into a counterfactual propensity score \(P'(Z)\) for each individual, and identify those untreated individuals whose counterfactual propensity score exceeds their original one. These are the compliers of the policy, and the PRTE is their average MTE. We do this below, re-estimating our original logit and “mapping” our new policy into a \(P(Z')\).

# Re-estimate logit (same spec as above)
logit_ps2 <- glm(form_ps, data = data_mte, family = binomial(link = "logit"))

# Predict original linear index
data_mte$xb_orig_idx <- predict(logit_ps2, type = "link")

# Counterfactual: move everyone 0.1km closer, floor at zero
data_mte$kmsmp_policy <- pmax(data_mte$kmsmp - 0.1, 0)

# Swap out kmsmp contribution only (INT_* left unchanged)
kmsmp_coef <- coef(logit_ps2)["kmsmp"]
data_mte$xb_policy_idx <- (data_mte$xb_orig_idx
                            - kmsmp_coef * data_mte$kmsmp
                            + kmsmp_coef * data_mte$kmsmp_policy)

# Counterfactual propensity score
data_mte$ps_policy <- plogis(data_mte$xb_policy_idx)

# Identify shifted individuals
data_mte$shifted <- as.integer(
  data_mte$dschool == 0 & data_mte$ps_policy > data_mte$ps_manual
)
table(data_mte$shifted)

   0    1 
2088  488 

# PRTE: average MTE over shifted individuals
data_mte %>% filter(shifted == 1) %>%
  summarise(PRTE = mean(margin), n = n())

The PRTE for the distance instrument identifies the return to schooling for individuals who would change their enrolment decision in response to a modest reduction in distance to their nearest secondary school. Despite the small size of the policy shift, the strong relationship between distance and enrolment in this setting means that a substantial number of currently untreated individuals are moved into treatment, and these compliers exhibit large and positive marginal returns. This suggests that the policy is reaching individuals who, while not currently enrolled, sit in a region of relatively low resistance and high returns, closer to our “low resistance” individuals considered previously than the higher resistence group.

References

Andresen, Martin Eckhoff. 2018. “Exploring Marginal Treatment Effects: Flexible Estimation Using Stata.” The Stata Journal 18 (1): 118–58.

Carneiro, Pedro, James J. Heckman, and Edward J. Vytlacil. 2011. “Estimating Marginal Returns to Education.” American Economic Review 101 (6): 2754–81. https://doi.org/10.1257/aer.101.6.2754.

Carneiro, Pedro, Michael Lokshin, and Nithin Umapathi. 2017. “Average and Marginal Returns to Upper Secondary Schooling in Indonesia.” Journal of Applied Econometrics 32 (1): 16–36. https://doi.org/https://doi.org/10.1002/jae.2523.

Cattaneo, Matias D., Luke Keele, Rocío Titiunik, and Gonzalo Vazquez-Bare. 2021. “Extrapolating Treatment Effects in Multi-Cutoff Regression Discontinuity Designs.” Journal of the American Statistical Association 116 (536): 1941–52.

Dehejia, Rajeev H., and Sadek Wahba. 1999. “Causal Effects in Nonexperimental Studies: Reevaluating the Evaluation of Training Programs.” Journal of the American Statistical Association 94 (448): 1053–62. http://www.jstor.org/stable/2669919.

Firpo, Sergio. 2007. “Efficient Semiparametric Estimation of Quantile Treatment Effects.” Econometrica 75 (1): 259–76.

LaLonde, Robert J. 1986. “Evaluating the Econometric Evaluations of Training Programs with Experimental Data.” The American Economic Review 76 (4): 604–20.

Londoño-Vélez, Juliana, Catherine Rodríguez, and and Fabio Sánchez. 2020. “Upstream and Downstream Impacts of College Merit-Based Financial Aid for Low-Income Students: Ser Pilo Paga in Colombia.” American Economic Journal: Economic Policy 12 (2): 193–227.