Boundary discontinuity designs on networks
with rd2dnet
rd2dnet estimates boundary discontinuity designs on
networks. It carries the location-based framework of
rd2d (Cattaneo, Titiunik and Yu) over from Euclidean
coordinates to network coordinates and user-defined network metrics:
side-specific local polynomials, robust bias-corrected inference,
uniform bands, data-driven bandwidths and aggregate effects. Please cite
the rd2d papers (see citation("rd2dnet")) when you use
it.
This vignette uses the bundled mock city, a small
synthetic street grid with a boundary along “60th Street” (treated side
to the south), a highway with ramps, and a park on the boundary. It is
only an illustration; nothing here depends on that application.
Everything runs offline.
library(rd2dnet)
city <- load_mock_midtown()
obs <- city$obs
y <- obs$outcome
head(obs)
#> x y treated far block outcome
#> 1 1592.404 863.317 TRUE 11.762 6010 0.020990
#> 2 102.434 52.888 TRUE 1.990 0 -0.404383
#> 3 2033.176 2920.818 FALSE 4.258 8036 1.577470
#> 4 1359.062 2992.232 FALSE 4.078 5037 2.234262
#> 5 2039.634 8.763 TRUE 2.826 8000 -0.290319
#> 6 2143.511 107.474 TRUE 2.429 8001 -1.057961
1. Data and network
city$nodes holds the street intersections with projected
x/y coordinates (metres);
city$edges the street segments with length,
highway class, speed_kph and
travel_time. city$boundary is a
boundary(): a line, the treated side, and a
mask for the park frontage, where one side has no
parcels and no effect can be estimated.
nd <- city$nodes
e <- city$edges
plot(nd$x, nd$y, type = "n", asp = 1, xlab = "x (m)", ylab = "y (m)",
main = "mock city: treated (red) / control (blue)")
segments(nd$x[match(e$from, nd$name)], nd$y[match(e$from, nd$name)],
nd$x[match(e$to, nd$name)], nd$y[match(e$to, nd$name)], col = "grey75", lwd = 0.6)
polygon(city$park, col = adjustcolor("forestgreen", 0.3), border = NA)
points(obs$x, obs$y, pch = 16, cex = 0.25, col = ifelse(obs$treated, "firebrick", "steelblue"))
lines(city$boundary$line, lwd = 2)
city$boundary
#> <rd2dnet boundary> length=3000.0, treated=region, masks=1
2. Distance is a metric you choose
A network is a graph plus a metric. Distance means
dissimilarity and does not have to be metres. Coordinates, kernel
distances and bandwidths are all in the metric’s units, while the
treatment side stays geographic.
metrics <- list(
"length (m)" = "length",
"travel time (s)" = travel_time(),
# 50 extra units per unit change in floor-area ratio between adjacent nodes
"length + land-use contrast" = node_contrast("far", 50),
# switching road class costs 500 m; ramps are transparent connectors, so
# highway -> ramp -> highway is cheap but highway -> ramp -> street is not
"length + class switches" = transition_cost("length", class_attr = "highway",
default_penalty = 500, connectors = "motorway_link")
)
pts <- obs[1:5, c("x", "y")]
round(sapply(metrics, function(m) {
network_coordinates(city$network(m), city$boundary, pts, n_eval = 5)$s
}), 1)
#> length (m) travel time (s) length + land-use contrast
#> [1,] 829.1 78.1 1487.1
#> [2,] 1649.5 107.9 2190.1
#> [3,] -1354.0 -118.9 -2327.9
#> [4,] -1783.2 -178.7 -2926.7
#> [5,] 1648.4 149.7 3173.8
#> length + class switches
#> [1,] 829.1
#> [2,] 1649.5
#> [3,] -1354.0
#> [4,] -2283.2
#> [5,] 1648.4
A metric can also be any vectorised function of the edge table, for
example
function(edges) edges$length * ifelse(edges$highway == "motorway", 3, 1).
Any node-by-node dissimilarity matrix works too
(precomputed()). Observations then snap to the nearest
node, because per-edge costs are unknown.
base <- city$network("length")
D <- distances_from(base, seq_len(base$n_nodes)) # or your own matrix, in node order
pre <- city$network(precomputed(sqrt(D))) # e.g. a concave transform
round(network_coordinates(pre, city$boundary, pts, n_eval = 5)$s, 2)
#> Warning: precomputed(450 x 450) cannot split edges; make sure boundary nodes
#> exist in the graph
#> [1] 26.83 38.99 -36.88 -43.13 40.00
3. Network coordinates
network_coordinates() snaps each observation to the
street network on its own side and computes, as rd2d’s
(X1 - b1, X2 - b2):
s: signed network distance to the boundary (positive on
the treated side);- at each evaluation point
b: d(b), the
network distance used for kernel weights, and o(b), the
signed offset along the boundary.
Evaluation points default to an even grid along the boundary; masked
ones are skipped.
net <- city$network("length")
coords <- network_coordinates(net, city$boundary, obs, n_eval = 21)
coords
#> <rd2dnet coordinates> n=3000, treated=1764, eval_points=21, masked=7, locations=3000
head(as.data.frame(coords, j = 11))
#> x y treated s crossing_arc depth d o
#> 1 1592.404 863.317 TRUE 829.087 1900 16.683 1229.087 400
#> 2 102.434 52.888 TRUE 1649.546 400 27.112 2544.678 -1100
#> 3 2033.176 2920.818 FALSE -1353.994 2400 33.176 2253.994 900
#> 4 1359.062 2992.232 FALSE -1783.170 2150 32.232 2351.294 650
#> 5 2039.634 8.763 TRUE 1648.397 2400 8.763 2548.397 900
#> 6 2143.511 107.474 TRUE 1653.963 2650 27.474 2590.985 1150
head(eval_table(coords), 8)
#> arc node node_arc x y masked
#> 1 100 9060 100 -300 1600 FALSE
#> 2 240 9060 100 -300 1600 FALSE
#> 3 380 60 400 0 1600 FALSE
#> 4 520 60 520 120 1600 FALSE
#> 5 660 160 660 260 1600 FALSE
#> 6 800 260 800 400 1600 FALSE
#> 7 940 260 940 540 1600 FALSE
#> 8 1080 360 1080 680 1600 FALSE
4. Estimation
rdbwnet() ports rd2d’s direct plug-in bandwidth selector
(mserd by default; also msetwo,
imserd, imsetwo and the coverage-error-optimal
cer* variants). rdnet() selects bandwidths
itself when h is not given, fits order p
(default 1) for the point estimate and order q = p + 1 for
robust bias-corrected inference, and adds a uniform band.
bw <- rdbwnet(y, coords)
head(bw[, c("arc", "h0", "h1", "n0", "n1", "status")])
#> arc h0 h1 n0 n1 status
#> 1 100 NA NA NA NA singular
#> 2 240 NA NA NA NA singular
#> 3 380 880.3830 880.3830 1236 1764 ok
#> 4 520 931.7889 931.7889 1236 1764 ok
#> 5 660 725.6914 725.6914 1236 1764 ok
#> 6 800 1033.4306 1033.4306 1236 1764 ok
res <- rdnet(y, coords, seed = 0)
res
#> rd2dnet: network boundary discontinuity estimates (sharp)
#> metric=edge_cost("length") kernel=triangular p=1 q=2 vce=hc1
#> evaluation points: 21 estimated: 12
#> uniform band critical value: 2.699
#> arc h0 h1 n0 n1 tau se tau_bc ci_lo ci_hi pvalue
#> 100.000 NA NA 0 0 NA NA NA NA NA NA
#> 240.000 NA NA 0 0 NA NA NA NA NA NA
#> 380.000 645.519 645.519 82 100 2.221 0.597 2.5427 0.711 4.374 6.50e-03
#> 520.000 730.659 730.659 125 146 2.216 0.422 2.4247 1.138 3.712 2.22e-04
#> 660.000 560.790 560.790 108 117 2.224 0.413 2.4153 1.050 3.781 5.27e-04
#> 800.000 821.495 821.495 194 266 2.157 0.290 2.2103 1.097 3.324 9.99e-05
#> 940.000 743.512 743.512 163 233 2.142 0.314 2.0946 1.027 3.163 1.21e-04
#> 1080.000 678.826 678.826 95 191 2.066 0.360 1.7011 0.869 2.533 6.13e-05
#> 1220.000 NA NA 0 0 NA NA NA NA NA NA
#> 1360.000 NA NA 0 0 NA NA NA NA NA NA
#> 1500.000 NA NA 0 0 NA NA NA NA NA NA
#> 1640.000 NA NA 0 0 NA NA NA NA NA NA
#> 1780.000 NA NA 0 0 NA NA NA NA NA NA
#> 1920.000 NA NA 0 0 NA NA NA NA NA NA
#> 2060.000 NA NA 0 0 NA NA NA NA NA NA
#> 2200.000 743.565 743.565 121 213 0.671 0.330 0.1552 -0.982 1.292 7.89e-01
#> 2340.000 552.100 552.100 91 111 0.555 0.353 0.0236 -1.099 1.147 9.67e-01
#> 2480.000 637.297 637.297 120 147 0.591 0.328 0.0819 -0.929 1.093 8.74e-01
#> 2620.000 729.850 729.850 154 179 0.758 0.304 0.3458 -0.621 1.313 4.83e-01
#> 2760.000 866.820 866.820 179 191 0.793 0.322 0.5559 -0.497 1.609 3.01e-01
#> 2900.000 727.187 727.187 102 107 0.651 0.483 0.1603 -1.505 1.825 8.50e-01
#> status
#> no bandwidth
#> no bandwidth
#> ok
#> ok
#> ok
#> ok
#> ok
#> ok
#> masked
#> masked
#> masked
#> masked
#> masked
#> masked
#> masked
#> ok
#> ok
#> ok
#> ok
#> ok
#> ok
plot(res, true_tau = city$true_tau(res$table$arc),
main = "BATEC along the boundary: estimate, robust CIs, uniform band, truth")
Bandwidths can also be given directly: a scalar,
c(h_control, h_treated), one per evaluation point, or a
J x 2 matrix. weights multiplies the kernel
weights (for example floor area), and pooled = TRUE fits a
common-slope model.
res_fixed <- rdnet(y, coords, h = 600, weights = obs$far, seed = 0)
head(res_fixed$table[, c("arc", "tau_bc", "se_rb", "status")])
#> arc tau_bc se_rb status
#> 1 100 NA NA insufficient
#> 2 240 NA NA insufficient
#> 3 380 2.464088 0.8753670 ok
#> 4 520 2.935654 0.5621422 ok
#> 5 660 2.722338 0.5773783 ok
#> 6 800 2.586707 0.6496810 ok
5. Aggregate effects: WBATE and LBATE
As in rd2d’s summary(WBATE = ..., LBATE = TRUE):
- WBATE
w' tau: weights are normalised
to sum to one, and the standard error sqrt(w' Sigma w) uses
the joint covariance res$cov of the estimates across
evaluation points. Points that were not estimated (masked) must get
weight 0.
- LBATE
max_j tau(b_j): its interval
uses the uniform critical value, so it inherits the band’s joint
coverage (conservative for the maximum).
ok <- res$table$status == "ok"
summary(res, wbate = as.numeric(ok), lbate = TRUE)$aggregate
#> tau tau_bc se_rb ci_lo ci_hi pvalue
#> WBATE 1.420357 1.225953 0.3286006 0.5819076 1.869998 0.0001908501
#> LBATE 2.224411 2.542737 NA 0.6768318 5.064925 NA
# Weight by frontage instead: here, by the number of parcels near each point.
wbate(res, ifelse(ok, res$table$n0 + res$table$n1, 0))
#> tau tau_bc se_rb ci_lo ci_hi pvalue
#> 1.4461755880 1.2528537159 0.3265479181 0.6128315571 1.8928758746 0.0001247184
Averaging correlated neighbours helps less than independence would
suggest:
se_naive <- sqrt(sum((ok / sum(ok))^2 * ifelse(ok, res$table$se_rb, 0)^2))
c(with_covariance = wbate(res, as.numeric(ok))[["se_rb"]], ignoring_it = se_naive)
#> with_covariance ignoring_it
#> 0.3286006 0.1820890
6. Covariates, fuzzy designs, derivatives
Covariates enter with coefficients common to both
sides (rd2d’s covs_eff) and are estimated jointly with the
local polynomials. In the mock city the outcome depends only weakly on
log floor-area ratio, so adjusting for it leaves estimates and standard
errors almost unchanged; a covariate that predicts the outcome well
shrinks the standard errors.
res_cov <- rdnet(y, coords, h = 600, covs = log(obs$far), seed = 0)
res_raw <- rdnet(y, coords, h = 600, seed = 0)
cmp <- data.frame(tau_bc = res_raw$table$tau_bc, tau_bc_covs = res_cov$table$tau_bc,
se_rb = res_raw$table$se_rb, se_rb_covs = res_cov$table$se_rb)
head(round(na.omit(cmp), 3))
#> tau_bc tau_bc_covs se_rb se_rb_covs
#> 3 2.563 2.579 0.954 0.973
#> 4 2.671 2.670 0.664 0.669
#> 5 2.339 2.338 0.660 0.662
#> 6 2.314 2.296 0.780 0.773
#> 7 1.849 1.836 0.650 0.648
#> 8 1.567 1.615 0.425 0.430
Fuzzy designs: pass take-up fuzzy = D.
The estimate is the ratio of the outcome jump to the take-up jump at
each point; the first stage is reported as
fs/fs_se. Below, take-up is 75% on the treated
side and 25% on the other (true first stage 0.5) and the effect of
take-up is 2, in a larger simulated city (12,000 parcels).
big <- make_mock_midtown(n_obs = 12000, seed = 11)
big_coords <- network_coordinates(big$network(), big$boundary, big$obs, n_eval = 21)
set.seed(2)
T <- big_coords$treated
D_takeup <- as.numeric(ifelse(T, runif(length(T)) < 0.75, runif(length(T)) < 0.25))
y_fuzzy <- 0.0008 * (big$obs$y - 1600) + 2 * D_takeup + rnorm(length(T), 0, 0.5)
res_fz <- rdnet(y_fuzzy, big_coords, fuzzy = D_takeup, seed = 0)
head(round(na.omit(res_fz$table[, c("arc", "fs", "fs_se", "tau", "ci_lo", "ci_hi")]), 3))
#> arc fs fs_se tau ci_lo ci_hi
#> 3 380 0.467 0.120 2.096 1.477 3.067
#> 4 520 0.461 0.068 1.998 1.316 2.361
#> 5 660 0.447 0.073 1.713 1.118 2.157
#> 6 800 0.427 0.077 1.730 0.699 2.133
#> 7 940 0.371 0.086 1.730 0.054 2.392
#> 8 1080 0.404 0.076 1.869 0.353 2.536
Derivatives: deriv = c(1, 0) is the
jump in the slope across the boundary (a kink),
deriv = c(0, 1) the jump in the slope along it,
and tangvec = c(t_s, t_o) any directional derivative. The
mock outcome has the same slope on both sides, so the true kink is
zero.
kink <- rdnet(y, coords, h = 800, deriv = c(1, 0), seed = 0)
head(signif(na.omit(kink$table[, c("arc", "tau_bc", "ci_lo", "ci_hi")]), 3))
#> arc tau_bc ci_lo ci_hi
#> 3 380 -4.29e-04 -0.00843 0.00757
#> 4 520 9.98e-04 -0.00621 0.00821
#> 5 660 -4.64e-05 -0.00592 0.00582
#> 6 800 4.18e-04 -0.00622 0.00706
#> 7 940 -9.11e-04 -0.00641 0.00459
#> 8 1080 -2.57e-03 -0.00804 0.00290
7. Standard errors
vce takes "hc0"-"hc3";
cluster gives cluster-robust errors;
conley_cutoff gives a Bartlett spatial HAC over
network distance (or
conley_distance = "euclidean"). The same choice drives the
uniform band and WBATE/LBATE through res$cov.
A caution on cluster: as in rd2d, cluster-robust errors
use the residuals of the local fit. When clusters are small relative to
the bandwidth (here, city blocks inside a 600 m window), the local
polynomial absorbs much of the within-cluster variation and the cluster
standard errors come out too small. For spatially correlated shocks, the
Conley option, whose residuals come from a wider auxiliary fit, is the
validated choice.
fits <- list(
HC1 = rdnet(y, coords, h = 600, seed = 0),
cluster_block = rdnet(y, coords, h = 600, cluster = obs$block, seed = 0),
conley_network = rdnet(y, coords, h = 600, conley_cutoff = 150, seed = 0),
conley_euclidean = rdnet(y, coords, h = 600, conley_cutoff = 150, conley_distance = "euclidean", seed = 0)
)
head(round(na.omit(sapply(fits, function(r) r$table$se_rb)), 3))
#> HC1 cluster_block conley_network conley_euclidean
#> [1,] 0.954 0.322 0.898 0.749
#> [2,] 0.664 0.431 0.660 0.602
#> [3,] 0.660 0.455 0.528 0.449
#> [4,] 0.780 0.454 0.621 0.491
#> [5,] 0.650 0.446 0.512 0.404
#> [6,] 0.425 0.257 0.393 0.370
8. Density (sorting) test
The local polynomial density test of Cattaneo, Jansson and Ma
(rddensity), as a density() method for network
coordinates. By default position is measured in cumulative
street length from the boundary
(normalize = "network"), so the test asks whether there are
more parcels per metre of street on one side, instead of reacting to
differences in street supply.
density(coords)
#> rd2dnet density test (network): jump=-1.98e-07 (bc 3.37e-06) se=3.2e-06 z=1.053 p=0.2923 log-ratio=-0.020
#> [n=903/896, in window 241/272, h=12249.6/12249.6]
density(coords, normalize = "none")
#> rd2dnet density test (none): jump=2.53e-05 (bc -5.15e-05) se=9.4e-05 z=-0.548 p=0.5834 log-ratio=0.643
#> [n=903/896, in window 136/144, h=326.8/326.8]
9. Recipes
These need no special API; they combine the pieces above.
Donut
Drop observations within delta (in metric units) of the
boundary. Missing outcomes are skipped.
y_donut <- ifelse(abs(coords$s) < 40, NA, y)
donut <- rdnet(y_donut, coords, h = 700, seed = 0)
head(round(na.omit(data.frame(full = res_fixed$table$tau_bc, donut_40m = donut$table$tau_bc)), 3))
#> full donut_40m
#> 3 2.464 1.340
#> 4 2.936 1.085
#> 5 2.722 1.435
#> 6 2.587 1.169
#> 7 2.050 1.539
#> 8 1.796 1.320
Covariate balance
Pass a pre-determined characteristic as the outcome; there should be
no jump. In the mock city, floor-area ratio is unrelated to
treatment.
bal <- rdnet(log(obs$far), coords, seed = 0)
summary(bal, wbate = as.numeric(bal$table$status == "ok"), lbate = TRUE)$aggregate
#> tau tau_bc se_rb ci_lo ci_hi pvalue
#> WBATE 0.005795486 0.06236909 0.1496254 -0.2308913 0.3556295 0.6767992
#> LBATE 0.105378732 0.24280099 NA -0.4586853 1.3672864 NA
Placebo boundary
Build a boundary() along a fake line inside one side and
use only that side’s observations. Here: a street 8 blocks north, on the
control side. It runs through the park, which is masked.
fake_y <- city$boundary_y + 8 * 80
north <- obs$y > city$boundary_y
fake <- boundary(rbind(c(-400, fake_y), c(2600, fake_y)),
treated_region = rect_polygon(-1e5, -1e5, 1e5, fake_y), mask = list(city$park))
pc <- network_coordinates(city$network(), fake, obs[north, ], n_eval = 15)
placebo <- rdnet(y[north], pc, seed = 0)
summary(placebo, wbate = as.numeric(placebo$table$status == "ok"))$aggregate
#> tau tau_bc se_rb ci_lo ci_hi pvalue
#> WBATE -0.004436142 0.1262539 0.2861236 -0.434538 0.6870458 0.6590273
Difference-in-discontinuities
With a panel (the same units before and after), run
on y_post - y_pre. With independent
samples, estimate each period with the same evaluation points
and bandwidths and difference the estimates, with
se = sqrt(se_pre^2 + se_post^2).
set.seed(5)
f <- y - city$true_tau(coords$crossing_arc) * coords$treated # outcome without the effect
y_pre <- f + rnorm(length(y), 0, 0.5) # same parcels, earlier period
y_post <- y + rnorm(length(y), 0, 0.5)
did_panel <- rdnet(y_post - y_pre, coords, h = 700, seed = 0)
c_pre <- make_mock_midtown(n_obs = 3000, seed = 21)
c_post <- make_mock_midtown(n_obs = 3000, seed = 22)
co_pre <- network_coordinates(c_pre$network(), c_pre$boundary, c_pre$obs, n_eval = 21)
co_post <- network_coordinates(c_post$network(), c_post$boundary, c_post$obs, n_eval = 21)
y0 <- c_pre$obs$outcome - c_pre$true_tau(co_pre$crossing_arc) * co_pre$treated
r_pre <- rdnet(y0, co_pre, h = 700, seed = 0)$table
r_post <- rdnet(c_post$obs$outcome, co_post, h = 700, seed = 0)$table
did <- data.frame(arc = r_post$arc, panel = did_panel$table$tau_bc,
independent = r_post$tau_bc - r_pre$tau_bc,
se_independent = sqrt(r_pre$se_rb^2 + r_post$se_rb^2),
truth = city$true_tau(r_post$arc))
head(round(na.omit(did), 3))
#> arc panel independent se_independent truth
#> 3 380 1.876 1.037 0.910 1.462
#> 4 520 1.577 0.621 0.739 1.723
#> 5 660 1.805 0.963 0.683 1.956
#> 6 800 2.156 1.538 0.651 2.133
#> 7 940 2.445 2.261 0.623 2.233
#> 8 1080 2.869 2.277 0.825 2.243
10. Large data
- Repeated measurements (sensor panels, many readings
per site) are snapped and measured once per distinct location: compare
coords$n_locations and coords$n_obs.
- Shortest paths, snapping and the Conley kernel run in compiled code.
Conley cost grows with the number of observation pairs closer than the
cutoff, which grows with the square of local density. With millions of
nearby points, a smaller cutoff or clustered errors are much
cheaper.
panel_xy <- as.matrix(obs[rep(1:500, each = 20), c("x", "y")]) # 500 sensors x 20 readings
pcoords <- network_coordinates(net, city$boundary, panel_xy, n_eval = 11)
c(n_obs = pcoords$n_obs, n_locations = pcoords$n_locations)
#> n_obs n_locations
#> 10000 500