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.695
#>       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.6794174 5.060673           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.4613147 1.3722724        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