Some computational redistricting methods: or, how to sniff out a gerrymander in a pinch

• Statewide historical election data
• New Mexico State Senate
• Building an ensemble
• Partition ensemble & historical votes
• Results
• Summary
• References

Statewide historical election data

The MGGG States project provides an open collection of precinct shapefiles for US states. Here, we consider the state of New Mexico; the shapefile for which includes precinct-level results for six statewide elections.

library(tidyverse)
nm_sf <- sf::st_read(dsn = paste0(gitdir, # local folder 
                                  fname), # shapefile name from MGGG --
                     layer = fname, 
                     quiet = TRUE) %>% 
  sf::st_transform(crs = sf::st_crs('NAD83'))

## as simple df for viewing -- 
nm_df <- nm_sf %>%
  data.table::data.table() %>%  
  select(-geometry) 

Column names below. So, lots of information compiled/included in the shapefile; in addition to district details and election returns, precinct-level race and ethnicity counts are also made available.

colnames(nm_sf)
##  [1] "NAME10"    "prec_num"  "County"    "AG18D"     "AG18R"     "AG18L"    
##  [7] "SOS18D"    "SOS18L"    "SEN18R"    "SEN18D"    "SEN18L"    "PRES16L"  
## [13] "PRES16R"   "PRES16D"   "SOS16R"    "SOS16D"    "GOV18D"    "GOV18R"   
## [19] "CDDIST"    "HDIST"     "SDIST"     "TOTPOP"    "NH_WHITE"  "NH_BLACK" 
## [25] "NH_AMIN"   "NH_ASIAN"  "NH_NHPI"   "NH_OTHER"  "NH_2MORE"  "HISP"     
## [31] "H_WHITE"   "H_BLACK"   "H_AMIN"    "H_ASIAN"   "H_NHPI"    "H_OTHER"  
## [37] "H_2MORE"   "VAP"       "HVAP"      "WVAP"      "BVAP"      "AMINVAP"  
## [43] "ASIANVAP"  "NHPIVAP"   "OTHERVAP"  "X2MOREVAP" "Area"      "Perimeter"
## [49] "SOS18R"    "geometry"

In the following analyses, we focus on voting results from four historical statewide elections:

geos <- c('NAME10', 'prec_num', 'County', 'SDIST')

short <- c('SEN18', 'GOV18', 'PRES16', 'SOS16')

election <- c('2018 Senate',
              '2018 Governor',
              '2016 President',
              '2016 Secretary of State')

Lastly, some data re-structuring. The resulting returns data frame contains precinct-level vote returns disaggregated by party and election, in long format.

ecs <- grep(paste0(esum$short, collapse = '|'), colnames(nm_df), value = T)

returns <- data.table::melt(nm_df, geos, ecs) %>%
  mutate(party = gsub('^.*[0-9]', '', variable),
         election = gsub('.$', '', variable)) %>%
  select(-variable)

New Mexico State Senate

The existing legislative map for the New Mexico State Senate was drawn during the 2010 cycle and, after a hiccup or two, went into effect in time for the general election in 2012. This partition will be referred to as NM2012.

#NAME10, partition, district
partition2012 <- nm_df %>%
  mutate(partition = 'NM2012') %>%
  mutate(district = as.character(as.numeric(SDIST))) %>%
  select(NAME10, partition, district)

Election results from 2016 are available via the nmelectiondatr package – an election in which all 42 seats were up for election. Dems won 26/42 seats. Note that Senate elections were most recently held in 2020; however, for our purposes here, 2016 is fine.

nm16 <- nmelectiondatr::nmel_results_summary %>%
  filter(Type == 'State Senator',
         Year == '2016') %>%
  rename(district = Type_Sub,
         value = Votes,
         party = Party) %>%
  mutate(party = substr(party, 1,1),
         partition = 'NM2016')

Building an ensemble

The redist package includes multiple approaches to randomly sampling graph partitions; here, we implement the Sequential Monte Carlo (SMC) method (McCartan and Imai 2020) via the redist.smc function. The SMC method requires two intermediary spatial data structures to build our ensemble of redistricting plans: (1) a distance matrix, which is derived from centroids of the precinct-level shapefile, and (2) an adjacency matrix.

centers <- sf::st_centroid(nm_sf) # %>% select(NAME10))
distmat <- sf::st_distance(centers)
attr(distmat, "units") <- NULL
distmat <- sqrt(distmat)

adj_list <- redist::redist.adjacency(nm_sf)

Parameters of relevance: (1) popvec: a vector specifying precinct-level voting age population counts; (2) ndists: the number of legislative boundaries to create; (3) nsim: the number of partitions to create; (4) popcons: the weight of population constraint; and (5) compactness: the weight of the compactness constraint.

Building the ensemble, then, amounts to piecing together 42 districts from 1,483 precincts – 10,000 times over – while abiding by the population & compactness constraints.

set.seed(9999)
smch <- redist::redist.smc(adjobj = adj_list,
                          popvec = nm_sf$VAP,
                          nsims = 10000,
                          ndists = 42,
                          popcons = 0.2, 
                          compactness = 1) 

## Sampling 10000 1483-unit maps with 42 districts and population tolerance 20%.

The ensemble of 10K partitions is included in model output as a matrix. Partitions are not shapefiles or geometries; instead, each partition is represented as a vector in which each precinct is assigned a district number.

ensembles <- smc$cdvec %>% 
  data.table::data.table() %>%
  mutate(NAME10 = nm_df$NAME10) %>%
  data.table::melt(., 'NAME10', c(1:ncol(smc$cdvec))) %>%
  rename(partition = variable,
         district = value) %>%
  mutate(partition = as.character(partition),
         district = as.character(district)) %>%
  bind_rows(partition2012) ###

## sample --
set.seed(99)
samp4 <- c('NM2012', 
           sample(unique(ensembles$partition), size = 3))

Sample of partition matrix re-structured as a data frame:

lc <- tmaptools::geocode_OSM (q = 'Albuquerque, NM', as.sf = T)
lc$bbox <- sf::st_set_crs(lc$bbox, sf::st_crs(nm_sf))
cropped <- sf::st_crop(nm_sf, lc$bbox)

The map below illustrates an example partition, zoomed into the Albuquerque Metro Area for a better look at things. Color reflects precinct-level district assignment per the SMC algorithm.

dists <- 42
nc <- cropped %>%
  select(-prec_num:-SOS18R) %>%
  left_join(ensembles %>% filter(partition %in% samp4), 
            by = 'NAME10') 

nc %>%
  filter(partition == samp4[[2]]) %>%
  ggplot() + 
  geom_sf(aes(fill = district),
          color = 'white',
          alpha = .65,
          lwd = .1) +
  scale_fill_manual(
    values = colorRampPalette(ggthemes::stata_pal()(8))(dists)) +
  theme_minimal() + map_theme() +
  labs(title = paste0('Partition ', samp4[[2]]), 
       subtitle = 'ABQ Metro Area')

We then aggregate precincts by district assignment to build out formal legislative boundaries for partition V5552.

nc %>% 
  filter(partition == samp4[[2]]) %>%
  group_by(partition, district) %>%
  summarise(geometry = sf::st_union(geometry)) %>%
  ungroup() %>%
  
  ggplot() + 
  geom_sf(aes(fill = district),
          color = 'white',
          alpha = .65,
          lwd = .3) +
  geom_sf_text(aes(label = district),
               color = 'black',
               size = 2.5, 
               check_overlap = TRUE) +
  scale_fill_manual(
    values = colorRampPalette(ggthemes::stata_pal()(8))(dists)) +
  theme_minimal() + map_theme() +
  labs(title = paste0('Partition ', samp4[[2]]), 
       subtitle = 'ABQ Metro Area')

Partition ensemble & historical votes

returns1 <- returns %>% 
  select(NAME10, election, party, value) %>%
  filter(party %in% c('R', 'D'))

From redistricting plan to election result, then, is fairly straightforward. The first step is to assign each partition in our ensemble voting data – again, using election results from four historical elections in New Mexico.

ps <- unique(ensembles$partition)
fens <- list()

for (q in 1:length(ps)) {
  x <- subset(ensembles, partition == ps[q])
  x1 <- x[returns1, on = 'NAME10']
  fens[[q]] <- x1[ , list(value = sum(value)),
      by = list(election, partition, district, party)] 
  print(q)
}

x2 <- fens %>% data.table::rbindlist()

As a simple starting point, the map below details precinct-level election results in the ABQ Metro Area for the 2018 statewide race for Governor. The darker the shade of blue, the larger the vote margin in favor of the Democratic candidate (here MLG).

c1 <- cropped %>% mutate(g18 = GOV18R - GOV18D)

c1 %>% 
  ggplot() + 
  geom_sf(aes(fill = g18),
          color = 'black',
          alpha = .85,
          lwd = .05) +
   scale_fill_distiller(palette = "RdYlBu",
                        limit = max(abs(c1$g18)) * c(-1, 1)) +
  theme_minimal() + map_theme() +
  labs(title = '2018 Gubernatorial election results', 
       subtitle = 'ABQ Metro Area')

nc1 <- nc %>% 
  group_by(partition, district) %>%
  summarise(geometry = sf::st_union(geometry)) %>%
  ungroup() 

Then we overlay four partitions on top of these results, three from our randomly generated ensemble, along with the actual NM2012 partition.

ggplot(data = c1) +
    geom_sf(aes(fill = g18),
            #color = 'black',
            alpha = .85,
            lwd = .05) +
   scale_fill_distiller(palette = "RdYlBu",
                        limit = max(abs(c1$g18)) * c(-1, 1)) +
  geom_sf(data = nc1,
          fill = NA,
          color = 'black',
          lwd = .25) +
  facet_wrap(~partition) +
  theme_minimal() + map_theme() +
  labs(title = 'Redistricting plans as overlay', 
       subtitle = 'ABQ Metro Area')

winners <- x2 %>%
  filter(election == 'GOV18') %>%
  select(-election) %>%
  filter(partition %in% samp4) %>%
  
  ## append Actual partition --
  bind_rows(nm16 %>% select(partition, district, party, value)) %>%
  group_by(partition, district) %>%
  filter(value == max(value))%>% 
  ungroup() 

Lastly, we tabulate election results for our example redistricting plans by aggregating precinct-level voting outcomes from the GOV18 election – in comparison to existing plan NM2012. The basic gist, then, is to (1) fix precinct-level vote counts, (2) adjust partition, and (3) tabulate results.

nc1 %>%
  left_join(winners) %>%
  ggplot() +
  
  geom_sf(aes(fill = party),
          color = 'white',
          alpha = .85,
          lwd = .25) +
  geom_sf_text(aes(label = district),
               color = 'black',
               size = 2.5,
               check_overlap = TRUE) +
  
  scale_fill_manual(values = c('#678fc3', '#e76a53')) +
  theme_minimal() + 
  map_theme() +
  facet_wrap(~partition) +
  labs(title = 'Redistricting plans: party winners', 
       subtitle = 'ABQ Metro Area')

Results

So, 10K redistricting plans, four historical elections. In theory, if an existing plan – or some proposed plan – is on the level, it should result in a legislature similar in partisan composition to those generated by the SMC algorithm.

comp <- x2 %>%
  group_by(election, partition, district) %>%
  filter(party == party[which.max(value)]) %>%
  group_by(election, partition, party) %>%
  summarize(seats = n()) %>%
  ungroup()

ds <- comp %>% filter(party == 'D')

ds1 <- ds %>%
  group_by(election, seats) %>%
  count() %>%
  group_by(election) %>%
  mutate(pertot = round(n/sum(n)*100, 2))

ds %>%
  ggplot() +
  geom_histogram(aes(seats, fill = election), 
                 binwidth = .5) +
  
  geom_vline(xintercept = 26, 
             linetype = 3, 
             color = '#e37e00',
             size = 1) +
  scale_x_continuous(breaks=seq(min(ds$seats), max(ds$seats),1)) +
  facet_wrap(~election) + #, scales = 'free') +
  ggthemes::scale_fill_economist()+
  theme_minimal()+
  theme(legend.position = 'none') +
  xlab('Number of Democrats Elected') +
  ggtitle('Distributions of elected Democrats')

x3 <- x2 %>%
  group_by(election, partition, district) %>%
  mutate(per = round(value/sum(value)*100, 3)) %>%
  #ungroup() %>%
  select(-value) %>%
  spread(party, per) %>%
  group_by(election, partition) %>%
  mutate(rank = rank(D)) %>%
  ungroup()

vote_dist <- nm16 %>%
  filter(Winner == 'Winner')  %>% 
  select(district, party, Percent) %>%
  mutate(D = ifelse(party == 'R', 1 - Percent, Percent),
         # D = ifelse(D == 1, 0.75, D),
         # D = ifelse(D == 0, 0.25, D),
         rank = rank(D, ties.method = 'first'))

x3 %>%
  ggplot() +
  geom_boxplot(aes(factor(rank), 
                   D, 
                   color = election), 
               size = 0.5, 
               outlier.size = 0.2) + 
  
  geom_point(data = vote_dist,
             aes(x = factor(rank), 
                 y = D * 100),
             color = '#e37e00',
             size = 1) +
  
  ggthemes::scale_color_economist() +
  geom_hline(yintercept = 50, lty = 2, color = 'steelblue') +
  facet_wrap(~election, ncol = 2) +
  theme_minimal() +
  theme(axis.text.x=element_blank(),
        legend.position = 'none') +
  xlab('Districts ordered from least to most Democratic') +
  ylab('Percentage of votes for a Democrat') +
  
  ggtitle('Marginal distributions')

x4 <- x3 %>%
  group_by(election, rank) %>%
  
  mutate(med = median(D),
         diff = (D - med) ^ 2) %>%
  group_by(election, partition) %>%
  
  summarize(marg_dev = sqrt(sum(diff))) 

x5 <- x4 %>% filter(partition %in% c(samp4))

x4 %>%
  ggplot() +
  geom_density(aes(marg_dev, fill = election),
               alpha = .5, color = 'gray') +
  geom_vline(data = x5,
             aes(xintercept = marg_dev), 
             color = '#e37e00', 
             linetype = 2,    
             size = .5) +
  
  geom_text(data = x5,
            aes(x = marg_dev, y = 0.3, label = partition), 
            color = '#e37e00', size = 3,
            check_overlap = TRUE) +
  
  theme_minimal() + 
  theme(legend.position = 'none')+
  ggthemes::scale_fill_economist() +
  facet_wrap(~election) +
  ggtitle('Ranked-Marginal Deviation')

Summary

So, hopefully a useful resource. Again, best utilized as an R-based appendix for the work highlighted in the reference section. The end of this three-part post. opefully timely, with the decennial census complete and redistricting in progress.

References