Lee and Lemieux (p. 31, 2009) suggest the researcher to present the graphs while doing Regression discontinuity design analysis (RDD). They suggest the following procedure:

"...for some bandwidth $h$, and for some number of bins $K_0$ and $K_1$ to the left and right of the cutoff value, respectively, the idea is to construct bins ($b_k$,$b_{k+1}$], for $k = 1, . . . ,K = K_0$+$K_1$, where $b_k = c−(K_0−k+1) \cdot h.$"

c=cutoff point or threshold value of assignment variable
h=bandwidth or window width.


...then compare the mean outcomes just to the left and right of the cutoff point..."

..in all cases, we also show the ﬁtted values from a quartic regression model estimated separately on each side of the cutoff point...(p. 34 of the same paper)

My question is how do we program that procedure in Stata or R for plotting the graphs of outcome variable against assignment variable (with confidence intervals) for the sharp RDD.. A sample example in Stata is mentioned here and here (replace rd with rd_obs) and a sample example in R is here. However, I think both of these didn't implement the step 1. Note, that both have the raw data along with the fitted lines in the plots.

Sample graph without confidence variable [Lee and Lemieux,2009] Thank you in advance.

• In response to your flag, a good way to revive your question is to edit it and offer a bounty: This will bump your question and get more people interested in it. If you feel this question might be better served on Stack Overflow, let us know and we can migrate it for you.
– chl
Feb 9, 2013 at 20:06
• I would like this to be migrated to Stack Overflow. Feb 12, 2013 at 12:10
• Unfortunately, this question is too old to be migrated to Stack Overflow. I believe it belongs on Cross Validated but if you want to ask on Stack Overflow (putting emphasis on the programming aspect and providing a minimal reproducible example), let me know and I will close it here.
– chl
Feb 14, 2013 at 10:29
• You should use cmogram. It does everything you need. Apr 7, 2013 at 19:49

Is this much different from doing two local polynomials of degree 2, one for below the threshold and one for above with smooth at $K_i$ points? Here's an example with Stata:

use votex // the election-spending data that comes with rd

tw
(scatter lne d, mcolor(gs10) msize(tiny))
(lpolyci lne d if d<0, bw(0.05) deg(2) n(100) fcolor(none))
(lpolyci lne d if d>=0, bw(0.05) deg(2) n(100) fcolor(none)), xline(0)  legend(off)


Alternatively, you can just save the lpoly smoothed values and standard errors as variables instead of using twoway. Below $x$ is the bin, $s$ is the smoothed mean, $se$ is the standard error, and $ul$ and $ll$ are the upper and lower limits of the 95% Confidence Interval for the smoothed outcome.

lpoly lne d if d<0, bw(0.05) deg(2) n(100) gen(x0 s0) ci se(se0)
lpoly lne d if d>=0, bw(0.05) deg(2) n(100) gen(x1 s1) ci se(se1)

/* Get the 95% CIs */
forvalues v=0/1 {
gen ulv' = sv' + 1.95*sev'
gen llv' = sv' - 1.95*sev'
};

tw
(line ul0 ll0 s0 x0, lcolor(blue blue blue) lpattern(dash dash solid))
(line ul1 ll1 s1 x1, lcolor(red red red) lpattern(dash dash solid)), legend(off)


As you can see, the lines in the first plot are the same as in the second.

• @Dimitry: +1 for the solution. However, I would like to have the mean value for each bin (please run the stata example above) rather than the scatter plot showing raw values. CI is great. Feb 15, 2013 at 2:54
• I am not quite sure what you mean. I added coded showing how you get the smoothed means in each bin by hand. If that's not what you are looking for, please explain what you have in mind in more detail. As far as I can tell, these graphs usually show the raw data and the smoothed means. Feb 15, 2013 at 21:34
• To quote Lee and Lemieux (p. 31, 2009): "A standard way of graphing the data is to divide the assignment variable(d here) into a number of bins, making sure there are two separate bins on each side of the cutoff point (to avoid having treated and untreated observations mixed together in the same bin). Then, the average value of the outcome variable can be computed for each bin and graphed against the mid-points of the bins". So, if there are 50 bins, then we will have only 25 data points on the left and right and not all the raw data (e.g, Graph 6(b) of the reference: updated in question) Feb 15, 2013 at 22:39
• Now it's clear! I agree on the kernel. But are you certain it's now not degree 0? That would correspond to equally-weighted mean smoothing. Feb 15, 2013 at 22:51
• I believe that corresponds to lpoly with a regular kernel and a degree 0 polynomial Feb 19, 2013 at 4:35

Here's a canned algorithm. Calonico, Cattaneo, and Titiunik recently proposed a procedure for robust bandwidth selection. They implemented their theoretical work for both Stata and R, and it also comes with a plot command. Here's an example in R:

# install.packages("rdrobust")
library(rdrobust)
set.seed(26950) # from random.org
x<-runif(1000,-1,1)
y<-5+3*x+2*(x>=0)+rnorm(1000)
rdplot(y,x)


That will give you this graph:

• Hi, how to add CI? Oct 14, 2019 at 1:00
• @Krantz just use rdplot(y,x, ci=95) May 13, 2022 at 19:58