Introduction

This analysis examines redistricting in New Mexico through computational simulations of the state’s 42-member senate. Using the redist package, we generate 2,000 alternative redistricting plans that comply with New Mexico’s legal requirements: districts must maintain population balance within 5%, consist of contiguous precincts, and respect county boundaries where possible. These requirements are detailed in the legislative redistricting guidelines.

By applying 2024 election results to these simulated maps, we establish a baseline distribution of possible partisan outcomes under legally compliant alternatives. Comparing the existing map to this distribution reveals where it falls within the range of feasible configurations.

The precinct-level data and shapefiles come from the newmexico-political-data repository, and all code is included below for full reproducibility.

Data

This section loads the geographic boundaries and election data needed for redistricting analysis. All data are loaded directly from GitHub. The dataset includes precinct-level results and boundary files for redistricting analysis. Boundary files are available in the boundaries directory of the repository.

if (!requireNamespace("pacman", quietly = TRUE)) {
  install.packages("pacman")
}

suppressPackageStartupMessages(
  pacman::p_load(
    dplyr,
    ggplot2,
    tidyr,
    patchwork,
    data.table,
    redist,
    sf,
    tidycensus,
    tigris,
    purrr,
    DT
  )
)

Load Election Results

# Load precinct-level data from GitHub (2004-2024)
precinct_data <- read.csv("https://raw.githubusercontent.com/jaytimm/newmexico-political-data/main/data/elections/nm_precinct_results_2004-24.csv", 
                          stringsAsFactors = FALSE)

# Load final results data from GitHub
results <- read.csv("https://raw.githubusercontent.com/jaytimm/newmexico-political-data/main/data/elections/nm_election_results_2000-24.csv", 
                    stringsAsFactors = FALSE)

# Load county-level statewide results from GitHub
county_results <- read.csv("https://raw.githubusercontent.com/jaytimm/newmexico-political-data/main/data/elections/nm_county_statewide_results_2000-24.csv", 
                           stringsAsFactors = FALSE)

Load Boundaries

# Load joined precincts with district assignments GeoJSON from GitHub
vtd_with_districts <- st_read("https://raw.githubusercontent.com/jaytimm/newmexico-political-data/main/data/boundaries/nm_vtd_with_districts_2021.geojson") |>
  mutate(ID = as.character(ID))
## Reading layer `nm_vtd_with_districts_2021' from data source 
##   `https://raw.githubusercontent.com/jaytimm/newmexico-political-data/main/data/boundaries/nm_vtd_with_districts_2021.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 2163 features and 11 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -109.0504 ymin: 31.33216 xmax: -103.002 ymax: 37.00023
## Geodetic CRS:  WGS 84
# Process precinct data for redistricting analysis
# Standardize column names if needed
if("COUNTY_NAM" %in% names(precinct_data)) {
  precinct_data$county_name <- precinct_data$COUNTY_NAM
}
if("VTD_NUM" %in% names(precinct_data)) {
  precinct_data$precinct_name <- as.character(precinct_data$VTD_NUM)
}

# Filter to 2024 for current redistricting analysis
precinct_2024 <- precinct_data |>
  filter(election_date == 2024)

# Prepare election results for redistricting analysis
# Extract 2024 results by office
elections_2024 <- precinct_2024 |>
  filter(office_name %in% c("United States Senator", 
                             "President of the United States",
                             "United States Representative"))

Generate Redistricting Plans

We use the redist package to generate 2,000 simulated redistricting plans for the state senate. The simulations employ New Mexico’s legal requirements as constraints: population balance within 5%, contiguous precincts, and respect for county boundaries. The sequential Monte Carlo (SMC) algorithm explores the space of valid district configurations.

# Set up output folder
post_name <- "2025-11-28-new-mexico-redistricting-state-senate"
post_folder <- file.path("content", "posts", post_name)
plans_file <- file.path(post_folder, "redist_plans.rds")

# Check if plans already exist - try both relative and absolute paths
if (!file.exists(plans_file)) {
  plans_file <- file.path(getwd(), post_folder, "redist_plans.rds")
}

if (file.exists(plans_file)) {
  # Load existing plans
  plans <- readRDS(plans_file)
  # Recreate map object for plotting (needed for redist.plot.plans)
  adj <- redist::redist.adjacency(vtd_with_districts)
  map <- redist::redist_map(
    vtd_with_districts,
    adj = adj,
    ndists = 42,
    total_pop = POP,
    county = COUNTY_NAM,
    pop_tol = 0.05
  )
} else {
  # Create folder if it doesn't exist
  if (!dir.exists(post_folder)) {
    dir.create(post_folder, recursive = TRUE)
  }
  
  # Generate adjacency from precinct boundaries
  adj <- redist::redist.adjacency(vtd_with_districts)
  
  # Set up redistricting map object
  map <- redist::redist_map(
    vtd_with_districts,
    adj         = adj,
    ndists      = 42,          # State Senate seats
    total_pop   = POP,
    county = COUNTY_NAM,
    pop_tol     = 0.05
  )
  
  # Run SMC sampler
  plans <- redist::redist_smc(
    map,
    nsims = 2000
  )
  
  # Save plans for future runs
  saveRDS(plans, plans_file)
}

Example Simulated Maps

Six example simulated senate mappings from the 2,000 generated plans. Each map represents a different valid redistricting configuration that satisfies New Mexico’s legal requirements.

redist.plot.plans(plans, draws=c("1", "2", "3", "4", "5", "6"), shp = map)

Process Plans

The generated plans are stored as a matrix where each column represents a different redistricting plan and each row represents a precinct. We convert this to a data frame format for analysis, merging it with precinct identifiers.

# Convert plans matrix to data.frame
df <- attr(plans, "plans") |> data.frame()
df1 <- vtd_with_districts |> 
  st_drop_geometry() |>
  as.data.frame() |>
  select(GEOCODE, COUNTY_NAM, VTD_NUM) |>
  bind_cols(df) 

df2 <- df1 |> tidyr::pivot_longer(cols = 4:103)

Compute Partisan Splits

This function aggregates precinct-level election results to the district level for each redistricting plan, then determines which party wins each district based on vote totals.

# Function to compute D/R splits for one election
# Assumes elect_df and df1 have COUNTY_NAM and VTD_NUM as character columns
get_split <- function(elect_df, df1) {
    elect_df <- as.data.table(elect_df)
    
    # Merge precinct results with district assignments
    m <- merge(elect_df, df1, by = c("COUNTY_NAM", "VTD_NUM"), all.x = TRUE)
    
    plan_cols <- setdiff(colnames(df1), c("COUNTY_NAM", "VTD_NUM", "GEOCODE"))
    
    # For each plan, aggregate votes to districts and count D/R seats
    res <- lapply(plan_cols, function(p) {
        tmp <- m[!is.na(get(p))]
        if(nrow(tmp) == 0) return(data.table(seats_D = 0L, seats_R = 0L, plan = p))
        
        agg <- tmp[, .(dem = sum(dem, na.rm = TRUE), rep = sum(rep, na.rm = TRUE)), by = get(p)]
        setnames(agg, "get", "district")
        
        agg[, winner := fifelse(dem > rep, "D", "R")]
        agg[, .(seats_D = sum(winner == "D"), seats_R = sum(winner == "R"), plan = p)]
    })
    
    rbindlist(res)
}

Analyze Election Results

We evaluate each redistricting plan using three different election baselines: 2024 U.S. Senate, 2024 Presidential, and aggregated 2024 Congressional District results. For each plan, we calculate how many districts would be won by Democrats versus Republicans under each election baseline.

# Ensure merge columns are character type
df1$COUNTY_NAM <- as.character(df1$COUNTY_NAM)
df1$VTD_NUM <- as.character(df1$VTD_NUM)

# Define elections to evaluate
election_config <- list(
  list(name = "Senate2024", data = elections_2024 |> filter(office_name == "United States Senator")),
  list(name = "President2024", data = elections_2024 |> filter(office_name == "President of the United States")),
  list(name = "AggCD2024", data = elections_2024 |> 
       filter(office_name == "United States Representative" & 
              (district_name == "1" | district_name == 1 | 
               district_name == "2" | district_name == 2 |
               district_name == "3" | district_name == 3)))
)

# Ensure election data has character merge columns
for(i in seq_along(election_config)) {
  d <- election_config[[i]]$data
  if(!"COUNTY_NAM" %in% names(d) && "county_name" %in% names(d)) {
    d$COUNTY_NAM <- as.character(d$county_name)
  }
  if(!"VTD_NUM" %in% names(d) && "precinct_name" %in% names(d)) {
    d$VTD_NUM <- as.character(d$precinct_name)
  }
  if("COUNTY_NAM" %in% names(d)) d$COUNTY_NAM <- as.character(d$COUNTY_NAM)
  if("VTD_NUM" %in% names(d)) d$VTD_NUM <- as.character(d$VTD_NUM)
  election_config[[i]]$data <- d
}

# Evaluate all elections
final_results <- map_dfr(election_config, function(config) {
  result <- get_split(config$data, df1)
  result[, election := config$name]
  result
})

Interpreting the Results

Distribution of Democratic seats across 2,000 simulated redistricting plans using 2024 election results. The orange vertical line represents the existing map (26 Democratic seats). The existing map aligns with the simulated distribution for the presidential race, but for down-ballot races (Aggregate CD and Senate), it falls toward the lower end—most alternative plans would have yielded 27-30 Democratic seats.

plot_dt <- final_results[
    election %in% c("Senate2024","President2024","AggCD2024"),
    .(seats_D, election)
]

plot_dt[election == "Senate2024",    election := "Senate"]
plot_dt[election == "President2024", election := "President"]
plot_dt[election == "AggCD2024",     election := "Agg CD"]

enacted <- 26L

# Calculate overall x-axis range across all facets
xmin <- floor(min(plot_dt$seats_D))
xmax <- ceiling(max(plot_dt$seats_D))
xbreaks <- seq(xmin, xmax, by = 1)

ggplot(plot_dt, aes(x = seats_D)) +
    geom_histogram(
        binwidth = 1,
        color = "black",
        fill = "grey80",
        linewidth = 0.3,
        boundary = 0.5
    ) +
    geom_vline(
        xintercept = enacted,
        color = "#fc8d62",
        linewidth = 1
    ) +
    facet_wrap(~ election, nrow = 1) +
    scale_x_continuous(
        breaks = xbreaks,
        limits = c(xmin - 0.5, xmax + 0.5),
        minor_breaks = NULL
    ) +
    labs(
        x = "Democratic seats",
        y = "Count"
    ) +
    theme_minimal()

Analyzing Split Precincts and Senate District Effects

Most simulated redistricting plans would yield 28-30 Democratic state senate seats based on 2024 election results, while the existing map produces 26. To identify where this gap originates, we examine ticket-splitting: comparing state senate election results to federal election results (president, US Senate, US House) at the precinct level, then aggregating to see how these splits affect senate district outcomes.

Precinct-Level Ticket Splits

Winning party for each of the four elections at every precinct. The table shows all combinations of winners. While 84% of precincts voted straight-ticket, 16% showed splitting patterns. The two largest are mirror images: 103 precincts voted Republican for state senate but Democratic for federal races, while 80 precincts did the opposite.

# Get senate district assignments
vtd_senate <- vtd_with_districts |>
  st_drop_geometry() |>
  select(COUNTY_NAM, VTD_NUM, state_senate_district) |>
  mutate(COUNTY_NAM = as.character(COUNTY_NAM),
         VTD_NUM = as.character(VTD_NUM),
         state_senate_district = as.character(state_senate_district))

# Get precinct winners for each election
precinct_winners <- precinct_2024 |>
  filter(office_name %in% c("State Senate", "President of the United States", 
                            "United States Senator", "United States Representative")) |>
  group_by(COUNTY_NAM, VTD_NUM, office_name) |>
  summarise(dem = sum(dem, na.rm = TRUE), rep = sum(rep, na.rm = TRUE), .groups = "drop") |>
  #filter(dem > 0 | rep > 0) |>
  mutate(winner = ifelse(dem > rep, "D", "R")) |>
  mutate(office_key = case_when(
    office_name == "State Senate" ~ "winner_state_sen",
    office_name == "President of the United States" ~ "winner_pres",
    office_name == "United States Senator" ~ "winner_us_sen",
    office_name == "United States Representative" ~ "winner_house"
  )) |>
  select(COUNTY_NAM, VTD_NUM, office_key, winner) |>
  pivot_wider(names_from = office_key, values_from = winner)

# Merge with senate districts - ensure all columns exist
precinct_comparison <- vtd_senate |>
  left_join(precinct_winners, by = c("COUNTY_NAM", "VTD_NUM"))

# Count precinct ticket splits showing all race combinations
# Filter out rows with ANY NA before counting (ignore NAs in sums/counts)
precinct_splits_full <- precinct_comparison |>
  filter(!is.na(winner_state_sen) & !is.na(winner_pres) & !is.na(winner_us_sen) & !is.na(winner_house)) |>
  count(winner_state_sen, winner_pres, winner_us_sen, winner_house, name = "Count") |>
  mutate(
    StateSen = as.character(winner_state_sen),
    Pres = as.character(winner_pres),
    USSen = as.character(winner_us_sen),
    House = as.character(winner_house)
  ) |>
  select(StateSen, Pres, USSen, House, Count) |>
  arrange(desc(Count))

# Format table with red/blue coloring
DT::datatable(
  precinct_splits_full,
  options = list(
    dom = 't',
    paging = FALSE,
    searching = FALSE,
    ordering = TRUE,
    info = FALSE
  ),
  rownames = FALSE
) |>
  DT::formatStyle(
    columns = c("StateSen", "Pres", "USSen", "House"),
    backgroundColor = DT::styleEqual(
      levels = c("D", "R"),
      values = c("rgba(0, 102, 204, 0.3)", "rgba(204, 0, 0, 0.3)")
    )
  )

Senate District Splits

Precinct-level results aggregated to the state senate district level, showing how ticket-splitting aggregates into district outcomes and where state senate winners differ from federal election winners.

# Get vote totals by district for each election - aggregate using GeoJSON districts
district_votes <- precinct_2024 |>
  filter(office_name %in% c("State Senate", "President of the United States", 
                            "United States Senator", "United States Representative")) |>
  left_join(vtd_senate, by = c("COUNTY_NAM", "VTD_NUM")) |>
  filter(!is.na(state_senate_district)) |>
  group_by(state_senate_district, office_name) |>
  summarise(dem = sum(dem, na.rm = TRUE), rep = sum(rep, na.rm = TRUE), .groups = "drop") |>
  filter(dem > 0 | rep > 0) |>
  mutate(winner = ifelse(dem > rep, "D", ifelse(rep > dem, "R", NA))) |>
  mutate(office_key = case_when(
    office_name == "State Senate" ~ "winner_state_sen_dist",
    office_name == "President of the United States" ~ "winner_pres_dist",
    office_name == "United States Senator" ~ "winner_us_sen_dist",
    office_name == "United States Representative" ~ "winner_house_dist"
  )) |>
  select(state_senate_district, office_key, winner) |>
  pivot_wider(names_from = office_key, values_from = winner)

senate_district_results <- district_votes

# Count senate district splits showing all race combinations
# Filter out rows with ANY NA before counting (ignore NAs in sums/counts)
district_splits_full <- senate_district_results |>
  filter(!is.na(winner_state_sen_dist) & !is.na(winner_pres_dist) & !is.na(winner_us_sen_dist) & !is.na(winner_house_dist)) |>
  count(winner_state_sen_dist, winner_pres_dist, winner_us_sen_dist, winner_house_dist, name = "Count") |>
  mutate(
    StateSen = as.character(winner_state_sen_dist),
    Pres = as.character(winner_pres_dist),
    USSen = as.character(winner_us_sen_dist),
    House = as.character(winner_house_dist)
  ) |>
  select(StateSen, Pres, USSen, House, Count) |>
  arrange(desc(Count))

# Format table with red/blue coloring
DT::datatable(
  district_splits_full,
  options = list(
    dom = 't',
    paging = FALSE,
    searching = FALSE,
    ordering = TRUE,
    info = FALSE
  ),
  rownames = FALSE
) |>
  DT::formatStyle(
    columns = c("StateSen", "Pres", "USSen", "House"),
    backgroundColor = DT::styleEqual(
      levels = c("D", "R"),
      values = c("rgba(0, 102, 204, 0.3)", "rgba(204, 0, 0, 0.3)")
    )
  )

Districts Where State Senate Winners Differ from Federal Races

Districts where the state senate winner differs from federal election winners—the seats that “flip” relative to federal patterns. Three districts (12, 21, 28) elected Republican state senators while voting Democratic in federal races, accounting for the 2-4 seat difference between the existing map’s 26 Democratic seats and the 28-30 seats that most alternative maps would produce.

# Get counties for each district from GeoJSON data
district_counties <- vtd_with_districts |>
  st_drop_geometry() |>
  group_by(state_senate_district) |>
  summarise(Counties = paste(sort(unique(COUNTY_NAM)), collapse = ", "), .groups = "drop") |>
  mutate(state_senate_district = as.character(state_senate_district))

# Identify districts where state senate winner differs from federal winners
flip_districts <- senate_district_results |>
  filter(!is.na(winner_state_sen_dist) & !is.na(winner_pres_dist) & !is.na(winner_us_sen_dist) & !is.na(winner_house_dist)) |>
  mutate(
    # Check if state senate differs from federal races
    differs_from_pres = winner_state_sen_dist != winner_pres_dist,
    differs_from_us_sen = winner_state_sen_dist != winner_us_sen_dist,
    differs_from_house = winner_state_sen_dist != winner_house_dist,
    is_flip = differs_from_pres | differs_from_us_sen | differs_from_house
  ) |>
  filter(is_flip) |>
  mutate(state_senate_district = as.character(state_senate_district)) |>
  left_join(district_counties, by = "state_senate_district") |>
  mutate(
    StateSen = as.character(winner_state_sen_dist),
    Pres = as.character(winner_pres_dist),
    USSen = as.character(winner_us_sen_dist),
    House = as.character(winner_house_dist)
  ) |>
  select(District = state_senate_district, Counties, StateSen, Pres, USSen, House) |>
  arrange(District)

# Format table with red/blue coloring
DT::datatable(
  flip_districts,
  options = list(
    dom = 't',
    paging = FALSE,
    searching = FALSE,
    ordering = TRUE,
    info = FALSE
  ),
  rownames = FALSE
) |>
  DT::formatStyle(
    columns = c("StateSen", "Pres", "USSen", "House"),
    backgroundColor = DT::styleEqual(
      levels = c("D", "R"),
      values = c("rgba(0, 102, 204, 0.3)", "rgba(204, 0, 0, 0.3)")
    )
  )

Summary

This post demonstrates the redist package for computational redistricting analysis, using New Mexico’s 42-member state senate as a case study. We generated 2,000 simulated maps meeting New Mexico’s legal requirements. Most would produce 28-30 Democratic seats using 2024 down-ballot results; the existing map produces 26. Three districts (12, 21, 28) account for the gap—each elected Republican state senators while voting Democratic in federal races.