17  Part 3: Regression Discontinuity Design

TipHow to use the code in this section

The code blocks throughout this section contain inline annotations explaining what data structures go in, what the key function arguments control, and what to look for in the output. Click any “▶ Code” triangle to expand a block and read the annotations alongside the code. Treat them as a guided checklist: wherever you see a comment starting with # ──, it marks a decision point you would need to revisit with your own data.

17.1 The Core Idea

AlterEco’s retailer awards a “Certified Sustainable” badge to any coffee product scoring 70 or above on a third-party environmental audit (scored 0–100). Products at 69 miss the badge; products at 71 get it. Because the audit score has measurement noise, which side of the cutoff a product lands on is nearly random near 70. We exploit this near-randomness.

NoteWhy the jump at 70 is causal — the “as-good-as-random” intuition

Imagine two products with audit scores of 69 and 71. They have almost identical green credentials — similar supply chains, similar materials — but the auditor’s ruler is imprecise. Those two points of difference are largely measurement noise. Yet one gets the badge and the other doesn’t. Near the cutoff, badge assignment is as good as a coin flip.

This local near-randomness is what makes the jump causal. It is equivalent to a mini natural experiment: products just above and just below 70 are exchangeable in every way except for badge receipt. The jump in WTP at that threshold is therefore attributable entirely to the badge.

The cost of this elegance: the estimate is a local average treatment effect (LATE). It tells you the causal effect of receiving the badge for brands hovering around a 70-point score — not for the top-scoring 90s or the low-scoring 40s, which may respond to the badge very differently. This LATE is conceptually identical to the LATE from Part 1 of this module — both arise because the causal estimate is anchored to a specific sub-population (compliers near the threshold) rather than the full population.

A Module 2 connection: The “as-good-as-random near the cutoff” logic is local randomisation — a naturally occurring version of the randomised experiment studied in Module 2. The same assumption that makes experiments valid (exchangeability of treated and control) holds here, but only locally. In Module 2, randomisation made exchangeability a design guarantee; here, it is an empirical claim that must be verified through the diagnostics below.

A Module 1 connection: The running variable itself is a measurement. If audit scores are measured with substantial error — auditors apply different standards, or scores reflect lobbying as much as sustainability — then the near-threshold units may not be as similar as assumed. Measurement error in the running variable biases the RDD estimate toward zero (attenuation at the boundary) and inflates the apparent smoothness of pre-cutoff trends. This is the causal-inference analogue of the reliability problem in Module 1: a noisy ruler undermines the precision of every inference drawn from it.

A practical check to run first: If brands can game the audit score to land just above 70, the “as-good-as-random” assumption breaks down. The density test below checks for suspicious bunching right at the cutoff.

17.2 Simulating an RDD

▶ Simulate 1,200 products near the 70-point certification threshold
set.seed(2024)
N_rdd <- 1200; cutoff <- 70; true_jump <- 0.80
# ── YOUR DATA: replace audit_score with your running variable (the continuous
#    score that determines treatment assignment), cutoff with the actual policy
#    threshold, and WTP_rdd with your outcome variable.
#    df_rdd must contain: the running variable, a 0/1 treatment indicator
#    (= 1 if running variable >= cutoff), the outcome, and running = running - cutoff.
audit_score <- runif(N_rdd, 30, 100)
treat_rdd   <- as.integer(audit_score >= cutoff)

WTP_rdd <- 4.0 + 0.02*(audit_score-cutoff) - 0.0003*(audit_score-cutoff)^2 +
           true_jump*treat_rdd + rnorm(N_rdd, 0, 0.60)
WTP_rdd <- pmax(1, pmin(10, WTP_rdd))
df_rdd  <- tibble(audit_score, treat_rdd, WTP_rdd, running=audit_score-cutoff)
# ── CHECK: plot outcome vs. running variable (as below) — a visible jump at
#    the cutoff, with smooth trends on each side, is the visual signature of
#    a valid RDD. Irregular or absent jumps warrant investigating the DGP.

df_rdd |>
  ggplot(aes(x=audit_score, y=WTP_rdd, colour=factor(treat_rdd))) +
  geom_point(alpha=0.30, size=0.9) +
  geom_vline(xintercept=cutoff, linewidth=1, colour="grey30") +
  geom_smooth(method="lm", formula=y~poly(x,2), se=TRUE, alpha=0.2) +
  scale_colour_manual(values=c("0"=clr_ctrl,"1"=clr_eco),
                      labels=c("Below 70 (no badge)","Above 70 (badge)")) +
  annotate("text", x=71, y=9.2, label="Cutoff = 70", hjust=0, size=3.5, fontface="bold") +
  annotate("segment", x=70, xend=70, y=3.85, yend=4.85, colour="firebrick",
           arrow=arrow(ends="both", length=unit(.2,"cm")), linewidth=0.8) +
  annotate("text", x=71, y=4.35,
           label=sprintf("~$%.2f jump\n(badge effect)", true_jump),
           hjust=0, size=3, colour="firebrick") +
  labs(x="Audit Score (running variable)", y="Consumer WTP ($)", colour=NULL,
       title="Regression Discontinuity: the jump at 70 estimates the badge's causal effect") +
  theme_mod3()

17.3 Key Assumptions and How to Test Them

17.3.1 Assumption 1: No Manipulation

Code
# ── YOUR DATA: replace df_rdd$audit_score with your running variable column
#    and c=cutoff with your policy threshold value.
# ── KEY ARGS: rddensity(X, c) tests whether the density of the running variable
#    is continuous at the cutoff; a significant p-value (< .05) means units are
#    sorting — manipulating the running variable to land just above (or below)
#    the threshold, which invalidates the "as-good-as-random" RDD assumption.
dens_test <- rddensity(X=df_rdd$audit_score, c=cutoff)
summary(dens_test)

Manipulation testing using local polynomial density estimation.

Number of obs =       1200
Model =               unrestricted
Kernel =              triangular
BW method =           estimated
VCE method =          jackknife

c = 70                Left of c           Right of c          
Number of obs         677                 523                 
Eff. Number of obs    256                 209                 
Order est. (p)        2                   2                   
Order bias (q)        3                   3                   
BW est. (h)           13.827              12.101              

Method                T                   P > |T|             
Robust                0.6254              0.5317              


P-values of binomial tests (H0: p=0.5).

Window Length              <c     >=c    P>|T|
1.276     + 1.276          28      20    0.3123
2.553     + 2.479          44      41    0.8284
3.829     + 3.682          64      57    0.5856
5.105     + 4.885          89      80    0.5384
6.382     + 6.087         118     101    0.2796
7.658     + 7.290         140     117    0.1698
8.934     + 8.493         165     138    0.1351
10.211    + 9.696         185     168    0.3945
11.487    + 10.898        213     190    0.2731
12.763    + 12.101        235     209    0.2354
Code
# ── CHECK: look at the p-value in the summary output; p > .05 is required for
#    the no-manipulation assumption; also inspect the density plot for visible
#    bunching or a gap immediately around the cutoff.
invisible(rdplotdensity(dens_test, X=df_rdd$audit_score,
              title="Density test: checking for bunching at the 70-point cutoff",
              xlabel="Audit score", ylabel="Density"))

A significant p-value signals that brands are gaming audit scores to land just above 70. With our simulated uniform running variable, the density is continuous.

17.3.2 Assumption 2: Covariate Continuity

Code
env_rdd <- rnorm(N_rdd, 0 + 0.01*(audit_score-cutoff)/10, 1)
rdd_env <- rdrobust(y=env_rdd, x=df_rdd$audit_score, c=cutoff)
cat("Covariate balance — env. concern should NOT jump at the threshold:\n")
Covariate balance — env. concern should NOT jump at the threshold:
Code
summary(rdd_env)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 1200
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                  677          523
Eff. Number of Obs.             136          120
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   7.411        7.411
BW bias (b)                  11.566       11.566
rho (h/b)                     0.641        0.641
Unique Obs.                     677          523

=====================================================================
                   Point    Robust Inference
                Estimate         z     P>|z|      [ 95% C.I. ]       
---------------------------------------------------------------------
     RD Effect    -0.091    -0.122     0.903    [-0.604 , 0.533]     
=====================================================================

The covariate continuity test is the RDD’s randomisation check — conceptually identical to the balance tables you ran in Module 2 before analysing experimental data. In an experiment, balance across pre-treatment covariates is a consequence of random assignment; here it must be demonstrated empirically as evidence that the threshold is not systematically sorting units on unobserved characteristics.

17.4 Estimating the RDD Effect

Code
# ── YOUR DATA: replace y= with your outcome variable, x= with your running
#    variable, and c= with your cutoff value. Both y and x must be numeric vectors
#    of the same length with no missing values at the cutoff.
# ── KEY ARGS: rdrobust() selects bandwidth automatically (MSE-optimal rule);
#    you can override with h= to specify a fixed bandwidth, or kernel= to change
#    from the default triangular kernel (options: "epanechnikov", "uniform").
#    The default uses local linear regression on each side; p=2 switches to
#    local quadratic (more flexible, higher variance).
rdd_est <- rdrobust(y=df_rdd$WTP_rdd, x=df_rdd$audit_score, c=cutoff)
summary(rdd_est)
Sharp RD estimates using local polynomial regression.

Number of Obs.                 1200
BW type                       mserd
Kernel                   Triangular
VCE method                       NN

Number of Obs.                  677          523
Eff. Number of Obs.             181          172
Order est. (p)                    1            1
Order bias  (q)                   2            2
BW est. (h)                   9.821        9.821
BW bias (b)                  17.119       17.119
rho (h/b)                     0.574        0.574
Unique Obs.                     677          523

=====================================================================
                   Point    Robust Inference
                Estimate         z     P>|z|      [ 95% C.I. ]       
---------------------------------------------------------------------
     RD Effect     1.042     6.435     0.000     [0.762 , 1.429]     
=====================================================================
Code
# ── CHECK: in the summary output, look at "Conventional" coef for the point
#    estimate and "Robust" CI for the bias-corrected confidence interval —
#    use the Robust CI for inference. Also note the automatically chosen
#    bandwidth (h) and effective sample size on each side of the cutoff.
rdplot(y=df_rdd$WTP_rdd, x=df_rdd$audit_score, c=cutoff,
       title="RDD estimate: causal effect of sustainability badge on WTP",
       x.label="Audit Score", y.label="WTP ($)")

Code
# ── YOUR DATA: adjust the bandwidth grid c(5,8,10,15,20,25,30) to be
#    meaningful in your running variable's units — e.g., if your score runs
#    0–1000, sensible bandwidths might be c(25,50,75,100,150,200).
#    Keep the grid symmetric around the MSE-optimal h chosen by rdrobust() above.
# ── KEY ARGS: h=bw fixes the bandwidth; est$coef["Conventional",1] is the
#    point estimate; est$se["Robust",1] is the bias-corrected robust SE.
bw_results <- map_dfr(c(5,8,10,15,20,25,30), function(bw) {
  est <- rdrobust(y=df_rdd$WTP_rdd, x=df_rdd$audit_score, c=cutoff, h=bw)
  tibble(bandwidth=bw, estimate=est$coef["Conventional",1], se=est$se["Robust",1])
})
# ── CHECK: estimates should be roughly flat across bandwidths — large swings
#    indicate sensitivity to bandwidth choice (a red flag). If estimates collapse
#    at very narrow bandwidths, the effective N is too small for precision.

bw_results |>
  mutate(lo=estimate-1.96*se, hi=estimate+1.96*se) |>
  ggplot(aes(x=bandwidth, y=estimate, ymin=lo, ymax=hi)) +
  geom_hline(yintercept=true_jump, linetype="dashed", colour="grey40", linewidth=1) +
  geom_ribbon(alpha=0.15, fill=clr_eco) +
  geom_line(colour=clr_eco, linewidth=1) + geom_point(colour=clr_eco, size=2.5) +
  annotate("text", x=5.5, y=true_jump+0.05, label=sprintf("True = $%.2f", true_jump),
           hjust=0, size=3.5) +
  labs(x="Bandwidth (audit score units)", y="Estimated badge effect ($)",
       title="RDD estimates stable across bandwidths h = 5 to 30",
       subtitle="Instability would indicate sensitivity to bandwidth — a red flag") +
  theme_mod3()

NoteThe bandwidth trade-off: bias vs. variance

The bandwidth \(h\) controls which products are included in the RDD comparison — only those with audit scores in \([70-h,\ 70+h]\). This creates a fundamental tension:

Narrow \(h\) Wide \(h\)
Products are very similar to each other near the cutoff ✔ More observations → tighter CI ✔
Few observations → wide CI ✖ Products far from 70 differ systematically from each other ✖

rdrobust selects the bandwidth automatically using a mean-squared-error optimal rule — it finds the \(h\) that best balances bias and variance. The bandwidth-sensitivity plot above is a sanity check: stable estimates across a wide range of \(h\) suggest the result is not an artefact of a specific choice.

When estimates collapse or swing wildly as \(h\) changes, it usually means the polynomial fit on one side of the cutoff is extrapolating too aggressively. That’s a red flag that warrants narrowing \(h\) or using a local linear fit.