This function loops over all boundary points and locally estimates a non-parametric RD (using local linear regression)
using the rdrobust function from the rdrobust package from Calonico, Cattaneo, Titiunik (2014).
It takes in the discretized cutoff point file (the RDcutoff, a linestring chopped into parts by the discretise_border function)
and the sf object (which essentially is just a conventional data.frame with a geometry() column) containing all the observations (treated and untreated).
The treated indicator variable has to be assigned before (potentially with assign_treated) and be part of the sf object as a column.
spatialrd(
y,
data,
cutoff.points,
treated,
minobs = 50,
bwfix_m = NA,
sample = FALSE,
samplesize = NA,
sparse.exclusion = FALSE,
store.CIs = FALSE,
spatial.object = TRUE,
...
)The name of the dependent variable in the points frame in the form of a string
sf data.frame with points that describe the observations
sf object of borderpoints (provided by user or obtained with discretise_border)
column that contains the treated dummy (as string)
the minimum amount of observations in each estimation for the point estimate to be included (default is 50)
fixed bandwidth in meters (in case you want to impose one yourself)
draw a random sample of points (default is FALSE)
if random, how many points
in case we want to try to exclude sparse border points before the estimation (should reduce warnings)
set TRUE of confidence intervals should be stored
return a spatial object (deafult is TRUE, needed if you want to plot the point estimates on a map)?
in addition you can use all options in rdrobust
a data.frame or spatial data.frame (sf object) in case spatial.object = TRUE (default)
This function nests rdrobust. All its options (aside from running variable x and cutoff c) are available here as well (e.g. bw selection, cluster level, kernel, weights).
Check the documentation in the rdrobust package for details. (bandwidth selection default in rdrobust is bwselect = 'mserd')
To visualise the output, use plotspatialrd for a graphical representation. You can use printspatialrd (or an R package of your choice) for a table output. .
Calonico, Cattaneo and Titiunik (2014): Robust Nonparametric Confidence Intervals for Regression-Discontinuity Designs, Econometrica 82(6): 2295-2326.
points_samp.sf <- sf::st_sample(polygon_full, 1000) # create points
# make it an sf object bc st_sample only created the geometry list-column (sfc):
points_samp.sf <- sf::st_sf(points_samp.sf)
# add a unique ID to each observation:
points_samp.sf$id <- 1:nrow(points_samp.sf)
# assign treatment:
points_samp.sf$treated <- assign_treated(points_samp.sf, polygon_treated, id = "id")
#> Warning: attribute variables are assumed to be spatially constant throughout all geometries
# first we define a variable for the number of "treated" and control
NTr <- length(points_samp.sf$id[points_samp.sf$treated == 1])
NCo <- length(points_samp.sf$id[points_samp.sf$treated == 0])
# the treated areas get a 10 percentage point higher literacy rate
points_samp.sf$education[points_samp.sf$treated == 1] <- 0.7
points_samp.sf$education[points_samp.sf$treated == 0] <- 0.6
# and we add some noise, otherwise we would obtain regression coeffictions with no standard errors
points_samp.sf$education[points_samp.sf$treated == 1] <- rnorm(NTr, mean = 0, sd = .1) +
points_samp.sf$education[points_samp.sf$treated == 1]
points_samp.sf$education[points_samp.sf$treated == 0] <- rnorm(NCo, mean = 0, sd = .1) +
points_samp.sf$education[points_samp.sf$treated == 0]
# create distance to cutoff
points_samp.sf$dist2cutoff <- as.numeric(sf::st_distance(points_samp.sf, cut_off))
points_samp.sf$distrunning <- points_samp.sf$dist2cutoff
# give the non-treated one's a negative score
points_samp.sf$distrunning[points_samp.sf$treated == 0] <- -1 *
points_samp.sf$distrunning[points_samp.sf$treated == 0]
# create borderpoints
borderpoints.sf <- discretise_border(cutoff = cut_off, n = 10)
borderpoints.sf$id <- 1:nrow(borderpoints.sf)
# finally, carry out estimation alongside the boundary:
results <- spatialrd(y = "education", data = points_samp.sf, cutoff.points = borderpoints.sf,
treated = "treated", minobs = 20, spatial.object = FALSE)