Introduction to Rvoterdistance

library(sf)

# Simplified AZ-NM border: a vertical line at longitude -109.05
border <- st_sf(
  geometry = st_sfc(
    st_linestring(matrix(
      c(
        -109.05, 31.33,
        -109.05, 37.00
      ),
      ncol = 2, byrow = TRUE
    )),
    crs = 4326
  )
)

# Two voters: one in Albuquerque, one near the border
voters <- data.frame(
  name = c("Albuquerque voter", "Border voter"),
  lat  = c(35.08, 35.0),
  lon  = c(-106.65, -109.0)
)

d <- dist_to_boundary(voters, border,
  voter_coords = c("lat", "lon"),
  units = "km", progress = FALSE
)

data.frame(voter = voters$name, dist_km = round(d, 1))
#>               voter dist_km
#> 1 Albuquerque voter   218.6
#> 2      Border voter     4.6

# A rectangular district
district <- st_sf(
  geometry = st_sfc(
    st_polygon(list(matrix(c(
      -110, 35,
      -108, 35,
      -108, 37,
      -110, 37,
      -110, 35
    ), ncol = 2, byrow = TRUE))),
    crs = 4326
  )
)

# One voter inside, one outside
voters2 <- data.frame(
  name = c("Inside district", "Outside district"),
  lat  = c(36.0, 36.0),
  lon  = c(-109.0, -107.0)
)

d2 <- dist_to_boundary(voters2, district,
  voter_coords = c("lat", "lon"),
  units = "miles", progress = FALSE
)

data.frame(voter = voters2$name, dist_miles = round(d2, 1))
#>              voter dist_miles
#> 1  Inside district         56
#> 2 Outside district         56

voters_sf <- st_sf(
  name = c("Voter A", "Voter B"),
  geometry = st_sfc(
    st_point(c(-106.65, 35.08)),
    st_point(c(-109.00, 35.00)),
    crs = 4326
  )
)

d3 <- dist_to_boundary(voters_sf, border,
  units = "miles", progress = FALSE
)

data.frame(voter = voters_sf$name, dist_miles = round(d3, 1))
#>     voter dist_miles
#> 1 Voter A      135.8
#> 2 Voter B        2.8

voters_u <- data.frame(lat = 36.0, lon = -108.0)

d_km <- dist_to_boundary(voters_u, border,
  voter_coords = c("lat", "lon"), units = "km", progress = FALSE
)
d_mi <- dist_to_boundary(voters_u, border,
  voter_coords = c("lat", "lon"), units = "miles", progress = FALSE
)
d_m <- dist_to_boundary(voters_u, border,
  voter_coords = c("lat", "lon"), units = "meters", progress = FALSE
)

data.frame(km = round(d_km, 2), miles = round(d_mi, 2), meters = round(d_m, 1))
#>      km miles  meters
#> 1 94.56 58.76 94560.5

library(sf)
library(ggplot2)
library(rdrobust)

# Simplified Sandia Pueblo reservation boundary
sandia_coords <- matrix(c(
  -106.6140, 35.1850,
  -106.5400, 35.1800,
  -106.4800, 35.2000,
  -106.4500, 35.2350,
  -106.4450, 35.2700,
  -106.4600, 35.3050,
  -106.4900, 35.3200,
  -106.5500, 35.3250,
  -106.5900, 35.3100,
  -106.6100, 35.2750,
  -106.6200, 35.2350,
  -106.6140, 35.1850 # close the ring
), ncol = 2, byrow = TRUE)

sandia <- st_sf(
  name = "Sandia Pueblo",
  geometry = st_sfc(st_polygon(list(sandia_coords)), crs = 4326)
)

set.seed(2024)
n <- 50000

# Population centers: lat, lon, mixture weight, spatial spread
#   1. NE Albuquerque / Sandia Heights  (south of reservation, dense)
#   2. Bernalillo                        (northwest, medium)
#   3. On/near reservation               (sparse)
#   4. Rio Rancho                        (west, medium)
#   5. Placitas                          (northeast, small)
#   6. Corrales                          (west along the river)
centers <- data.frame(
  lat    = c(35.160, 35.310, 35.250, 35.275, 35.340, 35.235),
  lon    = c(-106.520, -106.560, -106.530, -106.680, -106.440, -106.630),
  weight = c(0.35, 0.18, 0.12, 0.17, 0.06, 0.12),
  sd_lat = c(0.035, 0.025, 0.035, 0.025, 0.018, 0.018),
  sd_lon = c(0.035, 0.025, 0.040, 0.025, 0.018, 0.018)
)

cluster <- sample(nrow(centers), n, replace = TRUE, prob = centers$weight)
voters <- data.frame(
  voter_id = seq_len(n),
  lat = rnorm(n, mean = centers$lat[cluster], sd = centers$sd_lat[cluster]),
  lon = rnorm(n, mean = centers$lon[cluster], sd = centers$sd_lon[cluster])
)

# Distance to the reservation boundary (unsigned, in miles)
dist_miles <- dist_to_boundary(
  voters, sandia,
  voter_coords = c("lat", "lon"),
  units = "miles",
  progress = FALSE
)

# Determine which voters fall inside the reservation
voters_sf <- st_as_sf(voters, coords = c("lon", "lat"), crs = 4326)
inside <- lengths(st_intersects(voters_sf, sandia)) > 0

# Signed score: positive inside, negative outside
voters$score <- ifelse(inside, dist_miles, -dist_miles)
voters$inside <- inside

cat("Voters inside reservation:", sum(inside), "\n")
#> Voters inside reservation: 15869
cat("Voters outside reservation:", sum(!inside), "\n")
#> Voters outside reservation: 34131
summary(voters$score)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
#> -10.3429  -3.0934  -1.2068  -1.2683   0.5282   4.5446

# Turnout probability: ~62% far outside, dropping to ~45% just inside
turnout_prob <- 0.62 + 0.005 * voters$score # gentle gradient
turnout_prob[inside] <- turnout_prob[inside] - 0.15 # discontinuity
turnout_prob <- pmin(pmax(turnout_prob, 0.10), 0.90) # clamp

voters$voted <- rbinom(n, 1, turnout_prob)

cat("Overall turnout rate:", round(mean(voters$voted), 3), "\n")
#> Overall turnout rate: 0.565
cat("Turnout inside:     ", round(mean(voters$voted[inside]), 3), "\n")
#> Turnout inside:      0.481
cat("Turnout outside:    ", round(mean(voters$voted[!inside]), 3), "\n")
#> Turnout outside:     0.604

ggplot(voters, aes(x = score)) +
  geom_histogram(
    aes(fill = inside),
    bins = 100, alpha = 0.8, boundary = 0
  ) +
  geom_vline(xintercept = 0, linetype = "dashed", linewidth = 0.8) +
  scale_fill_manual(
    values = c("FALSE" = "steelblue", "TRUE" = "firebrick"),
    labels = c("Outside reservation", "Inside reservation"),
    name = NULL
  ) +
  labs(
    x = "Distance to reservation boundary (miles)",
    y = "Number of voters",
    title = "Score variable: signed distance to Sandia Pueblo boundary"
  ) +
  theme_minimal() +
  theme(legend.position = "top")

# Plot a random sample of 5,000 voters for readability
set.seed(99)
voter_sample <- voters[sample(n, 5000), ]

ggplot() +
  geom_point(
    data = voter_sample,
    aes(x = lon, y = lat, color = score),
    size = 0.4, alpha = 0.6
  ) +
  geom_sf(
    data = sandia,
    fill = NA, color = "black", linewidth = 1
  ) +
  scale_color_gradient2(
    low = "steelblue", mid = "grey90", high = "firebrick",
    midpoint = 0, name = "Score\n(miles)"
  ) +
  labs(
    x = "Longitude", y = "Latitude",
    title = "Simulated voters around Sandia Pueblo reservation",
    subtitle = "Red = inside reservation, Blue = outside"
  ) +
  coord_sf(
    xlim = c(-106.78, -106.35),
    ylim = c(35.08, 35.42)
  ) +
  theme_minimal()

rd <- rdrobust(y = voters$voted, x = voters$score, c = 0)
summary(rd)
#> Call: rdrobust
#> 
#> Sharp RD estimates using local polynomial regression.
#> 
#> Number of Obs.                50000
#> BW type                       mserd
#> Kernel                   Triangular
#> VCE method                       NN
#> 
#>                                Left        Right
#> Number of Obs.                34131        15869
#> Eff. Number of Obs.           10951         8250
#> Order est. (p)                    1            1
#> Order bias (q)                    2            2
#> BW est. (h)                   1.457        1.457
#> BW bias (b)                   2.275        2.275
#> rho (h/b)                     0.641        0.641
#> Unique Obs.                   34131        15869
#> 
#> =====================================================================
#>                    Point    Robust Inference
#>                 Estimate         z     P>|z|      [ 95% C.I. ]       
#> ---------------------------------------------------------------------
#>      RD Effect    -0.140    -7.732     0.000    [-0.178 , -0.106]    
#> =====================================================================

rdplot(
  y = voters$voted,
  x = voters$score,
  c = 0,
  title = "Geographic RD: Voter turnout at Sandia Pueblo boundary",
  x.label = "Distance to reservation boundary (miles)",
  y.label = "Voter turnout"
)