Metadata-Version: 2.4
Name: rd2dnet
Version: 0.1.1
Summary: Boundary discontinuity designs on street networks: a network-native extension of rd2d (Cattaneo, Titiunik and Yu)
Author: Zach Winship
License-Expression: GPL-3.0-only
Project-URL: Homepage, https://github.com/zwinship/rd2dnet
Project-URL: Tutorial, https://github.com/zwinship/rd2dnet/blob/main/examples/tutorial.ipynb
Project-URL: Changelog, https://github.com/zwinship/rd2dnet/blob/main/CHANGELOG.md
Project-URL: Issues, https://github.com/zwinship/rd2dnet/issues
Project-URL: rd2d (upstream method), https://rdpackages.github.io/rd2d/
Keywords: regression discontinuity,boundary discontinuity,rd2d,spatial,networks,causal inference
Classifier: Development Status :: 3 - Alpha
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Operating System :: OS Independent
Classifier: Intended Audience :: Science/Research
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Python: >=3.10
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Requires-Dist: pandas>=2.0
Requires-Dist: networkx>=3.2
Requires-Dist: shapely>=2.0
Provides-Extra: osm
Requires-Dist: osmnx>=1.9; extra == "osm"
Requires-Dist: pyproj>=3.5; extra == "osm"
Provides-Extra: igraph
Requires-Dist: igraph>=0.11; extra == "igraph"
Provides-Extra: plot
Requires-Dist: matplotlib>=3.7; extra == "plot"
Provides-Extra: dev
Requires-Dist: pytest>=7; extra == "dev"
Requires-Dist: matplotlib>=3.7; extra == "dev"
Requires-Dist: nbformat>=5.7; extra == "dev"
Requires-Dist: nbclient>=0.7; extra == "dev"
Requires-Dist: ipykernel>=6; extra == "dev"
Dynamic: license-file

# rd2dnet

Boundary discontinuity designs on street networks. `rd2dnet` adapts the location-based
boundary discontinuity framework of [`rd2d`](https://github.com/rdpackages/rd2d)
(Cattaneo, Titiunik & Yu) from Euclidean coordinates to **network-native coordinates**
and **user-defined network metrics**.

> **Built on rd2d.** `rd2dnet` is a direct extension of the boundary discontinuity
> methodology and software of **Matias D. Cattaneo, Rocío Titiunik and Ruiqi Rae Yu**
> ([rd2d](https://github.com/rdpackages/rd2d), part of the
> [rdpackages](https://rdpackages.github.io/) project). The estimands, side-specific local
> polynomial estimation, robust bias-corrected inference, uniform bands and the
> bandwidth selector (a port of `rdbw2d`'s direct plug-in) follow their work. The
> density test follows Cattaneo, Jansson and Ma's `rddensity`. What is new here is
> carrying that framework over to network coordinates and user-defined network
> metrics. If you use `rd2dnet`, please cite the rd2d papers listed under
> [References](#references), as well as this package. `rd2dnet` is an independent
> project and is not affiliated with or endorsed by the rd2d authors.

## Installation

```bash
pip install rd2dnet              # core
pip install "rd2dnet[osm]"       # plus OpenStreetMap download (osmnx, pyproj)
pip install "rd2dnet[plot]"      # plus plotting (matplotlib)
```

Requires Python 3.10 or newer. Releases are on [PyPI](https://pypi.org/project/rd2dnet/).

**R package.** The same methods are available in R, in [`R-package/`](R-package/)
(same functions in snake case, same numbers on the bundled data):

```r
remotes::install_github("zwinship/rd2dnet", subdir = "R-package")
```

See [`R-package/README.md`](R-package/README.md) and `vignette("rd2dnet")`.

## Method

`rd2d` fits a local polynomial in `(X1 - b1, X2 - b2)` on each side of the boundary at
every evaluation point `b`, using kernel weights based on Euclidean distance. `rd2dnet`
replaces both pieces:

| | rd2d | rd2dnet |
|---|---|---|
| coordinates | `(X1 - b1, X2 - b2)` | `s_i`: signed network distance to the boundary (negative on the control side)<br>`o_i(b)`: signed offset along the boundary from `b` to where `i` reaches it |
| kernel distance | `‖X_i - b‖` | `d_i(b)`: shortest-path distance from `b` to `i` |
| estimate | side-specific intercept difference | same, with polynomial order `p` in `(s, o)` |
| inference | robust bias correction with `q = p + 1` | same, plus uniform bands |

Points snap to the nearest street segment on their own side. The off-network
"depth" is charged through `snap_scale`.

## Distance is a metric you choose

Each `NetworkGraph` carries the metric that defines distance on it. Distance here means
**dissimilarity**. It does not have to be metres.

```python
import rd2dnet
from rd2dnet import EdgeCost, TransitionCost, Precomputed, travel_time, node_contrast

rd2dnet.NetworkGraph(G, "length")                      # metres along streets
rd2dnet.NetworkGraph(G, travel_time())                 # seconds, from speed_kph
rd2dnet.NetworkGraph(G, node_contrast("far", 50.0))    # length + 50 * |land-use change|
rd2dnet.NetworkGraph(G, lambda u, v, d: ...)           # any non-negative edge cost
rd2dnet.NetworkGraph(G, Precomputed(D, nodes))         # your own node-by-node matrix

# Switching road class costs extra. Ramps are transparent connectors, so
# highway -> ramp -> highway is close while highway -> ramp -> local street is far.
rd2dnet.NetworkGraph(G, TransitionCost(
    "length", class_attr="highway",
    penalty={("motorway", "residential"): 2000.0},   # or default_penalty=..., or a callable
    connectors={"motorway_link", "trunk_link"},
))
```

The coordinates `s` and `o`, the kernel distances and the bandwidths are all in the
metric's units. Treatment side is always geographic.

## Quick start

```python
import rd2dnet
from rd2dnet.datasets import load_mock_midtown   # bundled, works offline

city = load_mock_midtown()
net = city.network("length")
coords = rd2dnet.network_coordinates(net, city.boundary, city.obs[["x", "y"]], n_eval=21)

res = rd2dnet.rdnet(city.obs["outcome"], coords,
                  weights=city.obs["far"],          # intensity / floor-area scaling
                  cluster=city.obs["block"])        # or conley_cutoff=150.0
print(res)          # pointwise BATEC table
res.plot()

rd2dnet.rdbwnet(city.obs["outcome"], coords)          # bandwidths
print(rd2dnet.density(coords))                         # sorting test per metre of street
rd2dnet.rdnet(city.obs["outcome"], coords, fuzzy=take_up)   # fuzzy design (take_up = your 0/1 or share)
```

The tutorial notebook [`examples/tutorial.ipynb`](examples/tutorial.ipynb) walks through
every feature on the mock city (metrics, coordinates, bandwidths, aggregate effects,
covariates, fuzzy designs, derivatives, cluster and Conley errors, the density test) and
the recipes for donut designs, covariate balance, placebo boundaries and
difference-in-discontinuities.

For real data, use `rd2dnet.NetworkGraph.from_place("Manhattan, New York")`, which needs
`pip install "rd2dnet[osm]"`. The query's top geocoding match is used: its outline when it
is an area of at least `min_area` (250,000 m² by default), otherwise (a neighbourhood
mapped as a point, a single building) the streets within `dist` metres of it (1,000 by
default). It falls back to the bundled mock network with a warning when offline. You can also pass your own graph:
`NetworkGraph(G)`, `NetworkGraph.from_adjacency(A, coords)` or `NetworkGraph.from_igraph(g)`.

### Features

- `Boundary(line, treated_region=... | treated_side="left"/"right", mask=[park_polygon, (a, b)])`:
  masked stretches with no mass, such as Central Park, are skipped by both the
  estimator and the density test.
- Edges that cross the boundary are split automatically, so every cross-side path
  passes through a boundary node. Points on a street segment that runs along the
  boundary count as already being on it.
- **Bandwidths** (`rdbwnet`): a port of rd2d's two-stage direct plug-in. It uses a
  rule-of-thumb pilot, a separate pilot bandwidth for the bias, regularisation
  (`scaleregul`) and a higher-order bias term (`scalebiascrct`). `bwselect` takes
  `mserd`, `msetwo` (one bandwidth per side), `imserd` or `imsetwo`, plus the
  coverage-error-optimal `cerrd`, `certwo`, `icerrd` and `icertwo`. You can also
  pass `h` as a scalar, `(h_control, h_treated)`, or one value per evaluation point.
- **Covariates** (`covs=`): the coefficients are common to both sides, as in rd2d's
  `covs_eff`, and estimated in one joint local regression whose influence function
  accounts for estimating them. The point estimate matches rd2d's plug-in exactly.
  Covariates that are constant in a local sample are dropped there.
- **Derivative targets**: `deriv=(1, 0)` gives the jump in the slope *across* the
  boundary (a kink) and `deriv=(0, 1)` the jump in the slope *along* it.
  `tangvec=(t_s, t_o)` gives any directional derivative. Bandwidth selection adapts,
  as in rd2d.
- **Fuzzy designs**: `rdnet(y, coords, fuzzy=D)` estimates the ratio of the outcome
  jump to the take-up jump at each boundary point. It reports the first stage `fs`
  and `fs_se`, and selects the bandwidth from the linearised outcome, as rd2d does.
- **Standard errors** (`vce`): `hc0` to `hc3`, `cluster` (CR1 with rd2d's
  `(n - 1)/(n - k) * g/(g - 1)` factor), or `conley` (Bartlett
  spatial HAC over **network** or Euclidean distance). Conley is computed block by
  block, so memory stays bounded, and it also drives the uniform band's joint
  covariance across boundary points. Conley residuals come from a smoother
  auxiliary fit (`conley_resid_scale`, default 3 × h, and `conley_resid_order`,
  default q). Residuals from the estimation fit itself absorb spatially correlated
  shocks and make the HAC variance far too small. The default was chosen by
  simulation across five error designs (see Validation).
- `weights`: parcel intensity (floor-area ratio (FAR), building height) multiplies the
  kernel weights.
- `pooled=True` gives the common-slope model `y ~ 1 + T + s + o`. The default is
  side-specific fits, as in rd2d.
- **Aggregate effects**: `res.wbate(weights)` gives the weighted boundary average
  treatment effect `w' tau` (weights normalised to sum to one) with
  `se = sqrt(w' Sigma w)`, and `res.lbate()` the largest effect `max_j tau(b_j)` with
  the interval `[max_j(tau_bc - c se), max_j(tau_bc + c se)]`, `c` the uniform critical
  value. Both follow rd2d's `summary(WBATE=, LBATE=)`; `res.summary(wbate=w, lbate=True)`
  prints them. `res.cov` always holds the J x J covariance `Sigma` of the bias-corrected
  estimates across evaluation points (NaN for points not estimated; its diagonal is
  `se_rb**2`), so aggregates account for the correlation between nearby points.
  Masked or unestimated points need weight 0.
- **Density test** (`density`): the local polynomial density estimator of
  Cattaneo, Jansson & Ma (as in `rddensity`), with jackknife standard errors and
  bias-corrected inference. The default `normalize="network"` measures position in
  *cumulative street length* from the boundary, so the test asks whether there are
  more parcels per metre of street on one side. On a street grid the raw density
  of `s` is jagged, and street supply differs across sides (for example a park
  frontage), so the unnormalised test (`normalize="none"`) can flag geometry as
  sorting. The default `bwselect="common"` uses one bandwidth for both sides.

## Validation (Monte Carlo on the bundled mock city)

| check | result |
|---|---|
| robust 95% CI, pointwise coverage | 0.925 (n = 4,000), 0.940 (n = 16,000) |
| uniform 95% band coverage | 0.873 (n = 4,000), 0.940 (n = 16,000) |
| first stage, fuzzy design | unbiased; HC standard error matches Monte Carlo sd within 5% |
| kink (`deriv=(1, 0)`) | exact without noise; unbiased, standard error matches Monte Carlo sd, coverage 0.94–0.955 |
| WBATE (equal weights), 95% CI coverage | 0.910 (n = 4,000), 0.947 (n = 16,000); ignoring the cross-point covariance gives 0.717 / 0.774 |
| WBATE standard error vs Monte Carlo sd | 0.128 vs 0.129 (n = 16,000); median 0.237 vs robust sd 0.262 (n = 4,000) |
| LBATE 95% CI coverage | 0.982 (n = 4,000), 0.997 (n = 16,000): conservative by construction |
| density test size at 5%, `network` / `none` (no-park city) | 0.030 / 0.065 |
| density test size at 5%, `network` / `none` (park city, masked) | 0.025 / 0.085 |

Spatially correlated errors (n = 8,000, h = 500 m, Conley cutoff = 2 × correlation scale;
each cell is pointwise / uniform-band coverage):

| error design | HC1 | Conley (default) |
|---|---|---|
| independent | 0.917 / 0.842 | 0.962 / 0.958 |
| smooth field, 100 m | 0.825 / 0.617 | 0.927 / 0.808 |
| smooth field, 200 m | 0.895 / 0.825 | 0.955 / 0.933 |
| block shocks, 150 m | 0.900 / 0.800 | 0.958 / 0.950 |
| block shocks, 300 m | 0.885 / 0.858 | 0.958 / 0.883 |

**Known limitations.**
- Coverage is below nominal when local samples are small (about 150 observations
  per side) and improves with n. rd2d notes the same.
- Conley standard errors are conservative (up to about 1.5 × the true sd) when the
  spatial correlation scale is close to the bandwidth. A residual-based HAC cannot
  fully separate correlated noise from the regression function at that scale, and
  the default errs toward wider intervals.
- The uniform band is weakest for short-range smooth correlation (0.81). Pointwise
  coverage there is 0.93.
- At n = 4,000 the WBATE error is heavy-tailed (a few replications with tiny local
  samples, max error about 8 with matching large standard errors). The plain Monte Carlo
  sd (0.42) then exceeds the mean SE (0.25); on the robust scale the SE is about 10%
  too small, and coverage is 0.91.
- Cluster-robust errors use the local fit's own residuals, as rd2d does. When clusters
  are small relative to the bandwidth they are too small: clustering by city block
  (n = 4,000, selected bandwidths, 400 replications) gives 95% coverage of 0.78 with
  independent errors (HC1: 0.93), and neither HC1 nor cluster errors capture block-level
  shocks (coverage 0.69 for both; their SEs are about half the Monte Carlo sd). The
  local polynomial absorbs within-cluster variation, the same effect that motivates
  the Conley default's wider auxiliary fit. For spatially correlated shocks use Conley.
- Observations at exactly the same location are at network distance 0 in the Conley
  kernel (as with `conley_distance="euclidean"`). Before the sparse kernel they were
  treated as `2 x snapping depth` apart, which understated the correlation between
  repeated readings at one site.
- The LBATE interval inherits the uniform band's joint coverage and is conservative
  for the maximum (0.98-1.00). Its point estimate `max_j tau_bc` is biased upward in
  small samples (+0.42 at n = 4,000, +0.05 at n = 16,000), as any maximum of noisy
  estimates is.

## Scale

Timings in seconds on one desktop (AMD Ryzen 7 3700X, 8 cores, 32 GB), 25 evaluation
points, bandwidths selected by `rdbwnet`, Conley cutoff 150 m (city), 200 m (suburb) or
1 km (river). The networks are synthetic (`benchmarks/synthetic.py`): a jittered grid
city with arterials, a motorway with ramps and a park on the boundary; an irregular
Delaunay suburb with cul-de-sacs and a curved boundary that cuts edges; and a dendritic
river basin (a tree with edges directed downstream). Observations are uniform
parcels, GPS-like pings near the streets, or a sensor panel (many readings at a few
thousand fixed sites).

| network (nodes) | observations | coords | `rdbwnet` | `rdnet` HC1 | cluster | Conley | `density` |
|---|---|---|---|---|---|---|---|
| city (50k) | 100k parcels | 2.1 | 0.9 | 0.2 | 0.2 | 1.1 | 0.7 |
| city (50k) | 1M parcels | 4.5 | 7.1 | 2.9 | 3.3 | 11.8 | 3.9 |
| city (50k) | 1M GPS pings | 4.0 | 6.5 | 2.7 | 2.7 | 12.9 | 4.1 |
| city (50k) | 5M GPS pings | 18.1 | 35.9 | 12.7 | 14.0 | 280 | 26.2 |
| city (50k) | 1M readings, 5k sensors | 1.0 | 5.4 | 1.2 | 1.4 | 4.2 | 3.0 |
| metro (200k) | 1M parcels | 7.3 | 9.0 | 2.6 | 3.2 | 8.7 | 6.8 |
| suburb (50k) | 200k parcels | 5.2 | 1.8 | 0.6 | 0.6 | 4.0 | 1.3 |
| river (50k) | 200k GPS pings | 2.1 | 1.3 | 0.5 | 0.5 | 2.8 | 0.5 |
| river (50k) | 10M readings, 10k sensors | 2.4 | 62.0 | 10.3 | 10.9 | 24.9 | 8.9 |

Building a 50k-node `TransitionCost` engine takes about 0.6 s (2.4 s for 200k nodes).

How the cost grows:
- **Coordinates**: snapping and boundary distances are computed once per distinct
  location (x, y, side), so a panel costs as much as its sites. Snapping and the
  geometry tests run on threads.
- **`rdbwnet` and `rdnet`**: linear in n per evaluation point. rd2d's bias pilot uses
  the median distance, so those fits see about half of each side's sample; this is
  the largest cost at 10M observations.
- **Conley**: proportional to the number of location pairs closer than the cutoff.
  That grows with the square of local density: 5x the observations in the same city
  is about 25x the pairs (5M pings, roughly 10^9 pairs, takes under 5 minutes). Node
  pairs come from Dijkstra searches bounded at the cutoff on the local subgraph,
  and the kernel is assembled in sparse row blocks of bounded memory. With very dense
  data, a smaller cutoff or clustered standard errors are far cheaper.
- **Density test**: O(n log n) plus a pass over the edges.

`RD2DNET_NUM_THREADS` caps the number of threads (default: up to 16 for geometry,
8 for the Conley kernel, 4 for evaluation points in `rdbwnet`; `1` disables
threading). Threaded results are identical to serial ones.

Reproduce with `python benchmarks/profile_scale.py --net grid_city --n 1000000
--design gps` (add `--profile STAGE` for a cProfile breakdown).

## Development

```bash
git clone https://github.com/zwinship/rd2dnet.git && cd rd2dnet
python -m venv .venv
.venv\Scripts\activate          # Windows   (source .venv/bin/activate on macOS/Linux)
pip install -e ".[dev]"
pytest                          # offline, includes doctests; `-m online` runs the live OSM test
pytest -m notebook              # executes examples/tutorial.ipynb
```

On Windows, check that `python` is your real interpreter and not one bundled with
another app (such as Inkscape). `py -0p` lists the installed interpreters.

GitHub Actions (`.github/workflows/ci.yml`) runs the tests on Python 3.10-3.13 (Linux) and
3.12 (Windows), executes the tutorial, and builds and checks the distribution.
`examples/make_tutorial.py` regenerates the notebook.

The R package is developed in `R-package/` (C++ core via Rcpp, roxygen2 documentation,
testthat tests including parity checks against this package's numbers):

```bash
cd R-package
Rscript -e 'Rcpp::compileAttributes(); roxygen2::roxygenise()'
R CMD build . && R CMD check --as-cran rd2dnet_*.tar.gz
```

`.github/workflows/R-CMD-check.yml` runs `R CMD check --as-cran` on Linux (R release,
devel and oldrel), macOS and Windows. `R-package/data-raw/` regenerates the bundled data
from this package's files.

## References

`rd2dnet` builds directly on the following work. Please cite it.

**rd2d (Cattaneo, Titiunik and Yu)**, see https://rdpackages.github.io/rd2d/

- Cattaneo, M. D., Titiunik, R. and Yu, R. R. (2026). rd2d: Causal Inference in Boundary
  Discontinuity Designs. Working paper. arXiv:2505.07989.
- Cattaneo, M. D., Titiunik, R. and Yu, R. R. (2026). Boundary Discontinuity Designs:
  Theory and Practice. *Advances in Economics and Econometrics: Thirteenth World
  Congress*, Cambridge University Press, Vol. 1, Ch. 2, to appear.
- Cattaneo, M. D., Titiunik, R. and Yu, R. R. (2026). Estimation and Inference in Boundary
  Discontinuity Designs: Location-Based Methods. *Journal of the American Statistical
  Association*, revise and resubmit.
- Cattaneo, M. D., Titiunik, R. and Yu, R. R. (2026). Estimation and Inference in Boundary
  Discontinuity Designs: Distance-Based Methods. *Journal of Econometrics* 256(A): 106266.
- Cattaneo, M. D., Titiunik, R. and Yu, R. R. (2026). Estimation and Inference in Boundary
  Discontinuity Designs: Pooling-Based Methods. Working paper.

**Density test (rddensity)**

- Cattaneo, M. D., Jansson, M. and Ma, X. (2020). Simple Local Polynomial Density
  Estimators. *Journal of the American Statistical Association* 115(531): 1449–1455.
- Cattaneo, M. D., Jansson, M. and Ma, X. (2018). Manipulation Testing Based on Density
  Discontinuity. *Stata Journal* 18(1): 234–261.

**Robust bias-corrected inference and spatial HAC**

- Calonico, S., Cattaneo, M. D. and Titiunik, R. (2014). Robust Nonparametric Confidence
  Intervals for Regression-Discontinuity Designs. *Econometrica* 82(6): 2295–2326.
- Conley, T. G. (1999). GMM Estimation with Cross Sectional Dependence. *Journal of
  Econometrics* 92(1): 1–45.

## License

GPL-3.0-only, the same license as rd2d and rddensity, from which parts of the
implementation are ported. See [LICENSE](LICENSE).
