---
title: "Boundary discontinuity designs on networks with rd2dnet"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Boundary discontinuity designs on networks with rd2dnet}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4)
```

`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.

```{r setup}
library(rd2dnet)
city <- load_mock_midtown()
obs <- city$obs
y <- obs$outcome
head(obs)
```

## 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.

```{r map, fig.height = 6}
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
```

## 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.

```{r metrics}
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)
```

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.

```{r precomputed}
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)
```

## 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.

```{r coords}
net <- city$network("length")
coords <- network_coordinates(net, city$boundary, obs, n_eval = 21)
coords
head(as.data.frame(coords, j = 11))
head(eval_table(coords), 8)
```

## 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.

```{r bandwidths}
bw <- rdbwnet(y, coords)
head(bw[, c("arc", "h0", "h1", "n0", "n1", "status")])
```

```{r estimate}
res <- rdnet(y, coords, seed = 0)
res
```

```{r plot}
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.

```{r fixed}
res_fixed <- rdnet(y, coords, h = 600, weights = obs$far, seed = 0)
head(res_fixed$table[, c("arc", "tau_bc", "se_rb", "status")])
```

## 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).

```{r aggregate}
ok <- res$table$status == "ok"
summary(res, wbate = as.numeric(ok), lbate = TRUE)$aggregate
# Weight by frontage instead: here, by the number of parcels near each point.
wbate(res, ifelse(ok, res$table$n0 + res$table$n1, 0))
```

Averaging correlated neighbours helps less than independence would suggest:

```{r correlated}
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)
```

## 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.

```{r covariates}
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))
```

**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).

```{r fuzzy}
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))
```

**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.

```{r kink}
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))
```

## 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.

```{r vce}
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))
```

## 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.

```{r density}
density(coords)
density(coords, normalize = "none")
```

## 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.

```{r donut}
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))
```

### 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.

```{r balance}
bal <- rdnet(log(obs$far), coords, seed = 0)
summary(bal, wbate = as.numeric(bal$table$status == "ok"), lbate = TRUE)$aggregate
```

### 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.

```{r placebo}
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
```

### 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)`.

```{r did}
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))
```

## 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.

```{r panel}
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)
```