Table of Contents

• The Flip MCMC Algorithm
• Strengths of Flip
• Initializing Flip
• Redistricting with Flip MCMC
• Using Multiple Chains
• Tidy Flip MCMC using redist_flip()
• Diagnostic Plots
• Tuning Flip Constraints
• Final Thoughts

The Flip MCMC Algorithm

data(fl25)
data(fl25_enum)
plan <- fl25_enum$plans[, 7241]
fl25$plan <- plan
fl_map <- redist_map(fl25, existing_plan = plan, pop_tol = 0.2, total_pop = pop)
#> Projecting to CRS 3857
constr <- redist_constr(fl_map) %>% 
    add_constr_edges_rem(0.02)
set.seed(1)
sims <- redist_flip(map = fl_map, nsims = 6, constraints = constr)
#> 
#> ── redist_flip() ───────────────────────────────────────────────────────────────
#> 
#> ── Automated Redistricting Simulation Using Markov Chain Monte Carlo ──
#> ℹ Preprocessing data.
#> ℹ Starting swMH().
#> 
■■■■■■                            17% | ETA:  0s | MH Acceptance: 1.00

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s | MH Acceptance: 1.00

The flip algorithm is one of the more straightforward redistricting algorithms. Beginning with an initial partition of a graph, it proposes flipping a node from one partition to an adjacent partition. By checking that the proposed flip meets basic constraints, such as keeping partitions contiguous and staying within a certain population parity, it ensures that all proposed new partitions are also valid partitions. The implementation within redist is a bit more advanced that this, as it allows for multiple flips and rejecting valid partitions based on a Metropolis Hastings algorithm. The following walks through the basics of this algorithm to provide an introduction to using flip correctly and efficiently.

Suppose we are redistricting this small map above on the left. To use the flip algorithm, we need to consider the adjacency graph that underlies this map, which is above on the right. Each of the 25 precincts on the left are displayed as a node on the right, connected if they are contiguous on the map. If we use the above district as an initial plan, we can then run flip for a few steps.

While this map is extremely small, the five iterations give a basic idea of what is going on behind the scenes. At each iteration, it searches the boundary for possible swaps, selects one, and then accepts or rejects the proposals. With very weak constraints, like those used to create the above example, almost every swap is accepted. Even then, though, it doesn’t guarantee that some iterations won’t repeat other plans sampled. In fact, in the above, the second iteration is the same plan as the initialization.

This possibility is very important for ensuring that the sampled plans are representative of the desired target distribution, which is controlled by the constraints chosen. The possible constraints are discussed below, as is information on setting up simulations and some advice on ensuring that your simulations are efficient.

Back to top

Strengths of Flip

Flip is incredibly powerful for local exploration. If you can make large changes to a summary statistic of interest without making large changes to the map itself, this may tell an important story of what went into making the map.

Flip is one of the easiest to understand algorithms and has theoretical guarantees behind it. This can make it especially useful when the audience of interest does not have an advanced background in mathematics or statistics.

Flip has the power to make less compact maps than many other algorithms. This can be especially powerful when a blind allegiance to compactness makes otherwise viable plans appear to be outliers.

Our implementation of flip has many more Gibbs constraints than our other implementations. This can allow you to consider different forms of partisan and countysplit constraints among others.

However, with these strengths do come weaknesses. Like most Markov Chain Monte Carlo methods, convergence can’t be shown, it can only be suggested. Diagnostics, like those in the section on diagnostic plots, can help ensure that convergence is likely, but can never show that it has indeed happened. Additionally, flip makes relatively small moves per iteration, so many more iterations are needed to move around the space. If your map is particularly large, you may require several hundred iterations to make the map substantively different, which leads to thinning the chain, which is dropping many sequential iterations. However, thinning doesn’t make the algorithm more efficient, so you still need to work through those plans, which comes with a time cost.

Initializing Flip

One of the keys to ensuring good performance is the choice of initialization. In some cases, a starting point may be obvious, such as when you want to explore the local area around an existing map. If that’s the use case, then it is straightforward to use that plan as the starting point. However, if the goal is to understand the larger space of possibilities, then starting from just one map can be misleading. Why? Since constraint tuning is not a perfect science, you could be setting the constraints too strong and, if that map is very good on some dimension, the flip algorithm may have difficulty getting away from that point without a very large number of iterations.

Our implementation defaults to using the Sequential Monte Carlo (SMC) algorithm via redist_smc() to create an initial partition of the districts, if no district is provided.

While the implementations of Random Seed and Grow (RSG) and Compact Random Seed and Grow (CRSG) via redist.rsg() and redist.crsg do not sample from a defined target distribution, they can serve as useful initializations for flip as they help provide a more diverse set of starting states. SMC is often faster and provides more theoretical guarantees, but tends to sample very compact districts, even when decreasing the compactness constraint. As such, when trying to decide if chains have likely converged or not, it can be misleading to only check chains that start from very compact states.

Redistricting with Flip MCMC

With the basics of what the flip algorithm is doing down, we can proceed into how to use the algorithm.

To begin running the MCMC algorithm, we have to provide some basic information, typically beginning with a shapefile. The below loads an Iowa dataset included within the redist package and plots the actual congressional districts from 2012-2021. (Iowa is a favorite choice for redistricting simulation examples, as it requires keeping counties together in plans which allows us to use the counties as the unit for redistricting, rather than thousands of precincts.)

data(iowa)
redist.plot.map(iowa, plan = cd_2010)

map_ia <- redist_map(iowa, existing_plan = cd_2010, pop_tol = 0.05)

From there, we need to build an adjacency graph which identifies which counties are touching which other counties on a map. If you have an existing plan, it’s generally advised to supply this to the optional plan argument to ensure that the existing plan is a valid, connected plan. If you get a warning, the geomander R package can help solve potential issues.

In addition, we need population for each unit. We’ve included iowa$pop as the total population as of the 2010 Census. From there, we have the basic information that we need to run our first simulation. The below indicates that we are simulating 1000 plans (with nsims) for the state of Iowa that have at most a population parity deviation of 0.05 (with pop_tol).

sims <- redist_flip(map_ia, nsims = 100)
#> 
#> ── redist_flip() ───────────────────────────────────────────────────────────────
#> 
#> ── Automated Redistricting Simulation Using Markov Chain Monte Carlo ──
#> ℹ Preprocessing data.
#> ℹ Starting swMH().
#> 
■                                  1% | ETA: 0s

■■■■■■■■■■■■■■■■■■■■■■            69% | ETA:  0s | MH Acceptance: 0.63

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s | MH Acceptance: 0.65

The printed output can be silenced by setting verbose = FALSE, however it displays very important information. First, it displays when preprocessing begins and when the algorithm actually starts. Each 10% of the way through the flip algorithm, it outputs the current estimated Metropolis acceptance. Here, we’ve specified no Gibbs constraints, so the acceptance will always be near 100%.

The output is an object of class redist.

class(sims)
#> [1] "redist_plans" "tbl_df"       "tbl"          "data.frame"

The sims object includes various pieces of information that were tracked while simulating, but we focus on get_plans_matrix(sims), which is a matrix that contains the plans.

dim(get_plans_matrix(sims))
#> [1]  99 101

Checking the dimensions shows that each plan is saved as a column, where each row is a precinct. From this, we can extract a single plan as we would from a normal matrix, like below, where we plot the final simulated plan.

redist.plot.map(shp = iowa, plan = get_plans_matrix(sims)[, 100])

Now, this plan is incredibly non-compact, which can be an issue. However, we should expect this type of outcome, as we didn’t include a compactness constraint while simulating. Thus, the only things checked were contiguity and that no plan would be outside of the pop_tol set above. Since there are many more non-compact plans than compact plans in the space of all redistricting plans, we end up with highly non-compact districts. We can fix this by specifying a constraint, as below:

constr <- redist_constr(map_ia) %>% add_constr_edges_rem(0.4)

sims_comp <- redist_flip(map_ia, nsims = 100, constraints = constr)
#> 
#> ── redist_flip() ───────────────────────────────────────────────────────────────
#> 
#> ── Automated Redistricting Simulation Using Markov Chain Monte Carlo ──
#> ℹ Preprocessing data.
#> ℹ Starting swMH().
#> 
■                                  1% | ETA: 0s

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s | MH Acceptance: 0.61

The first arguments as the same, but this adds three key arguments. First, setting constraint to any combination of the nine implemented constraints allows us to specify the target distribution. Setting constraintweights = 0.4 means that we want to put a relatively weak weight on the compactness, though a weak constraint still does a lot of work. There are four compact constraints implemented currently. The recommended is to use edges-removed because it can be calculated very quickly.

If we plot the final map sampled from the above code, we can see that it is far more compact.

redist.plot.map(shp = iowa, plan = get_plans_matrix(sims_comp)[, 100])

Back to top

Using Multiple Chains

When running larger redistricting analyses, one important step is to run multiple chains of the MCMC algorithm. This will also allow us to diagnose convergence better, using the Gelman-Rubin plot, as seen in the section on Diagnostic Plots.

On Windows and in smaller capacities, it is useful to run the algorithm within an lapply loop. First, we set up the seed for replicability and decide on the number of chains and simulations.

set.seed(1)
nchains <- 4
nsims <- 100

Here, we opt to initialize using the SMC algorithm. When we want to initialize without providing an initial partition, we need to specify the number of districts, ndists.

constr <- redist_constr(map_ia) %>% add_constr_edges_rem(0.4)
map_ia <- redist_map(iowa, ndists = 4, pop_tol = 0.05)
flip_chains <- lapply(1:nchains, function(x){
  redist_flip(map_ia, nsims = nsims,
              constraints = constr, verbose = FALSE)
})

In Unix-based systems, this can be run considerably faster by running this in parallel.

mcmc_chains <- parallel::mclapply(1:nchains, function(x){
  redist_flip(map_ia, nsims = nsims,
              constraints = constr, verbose = FALSE)
}, mc.set.seed = 1, mc.cores = parallel::detectCores())

Back to top

Tidy Flip MCMC using redist_flip()

The new, tidy interface to functions with redist introduces a pair of key objects, redist_map and redist_plans. The Get Started page goes into depth about these, but this shows the basics of how to work with the flip algorithm within the newer interface.

As in the standard interface, we need a data set to work with. This example will also follow with using the included Iowa data.

data(iowa)

Rather than building the adjacency graph manually, here we can set this up using redist_map which will build it an add it as a column.

iowa_map <- redist_map(iowa, existing_plan = cd_2010, pop_tol=0.01)

We set a population tolerance of 1%. While this is generally a good population parity tolerance for most simulations, be careful when using the default within flip. If your starting partition sits outside of that population deviation, flip may take a very, very long time to find a valid partition to flip.

Now, we can pass the redist_map object to redist_flip to begin simulating.

tidy_sims <- redist_flip(iowa_map, nsims = 100)
#> 
#> ── redist_flip() ───────────────────────────────────────────────────────────────
#> 
#> ── Automated Redistricting Simulation Using Markov Chain Monte Carlo ──
#> ℹ Preprocessing data.
#> ℹ Starting swMH().
#> 
■                                  1% | ETA: 0s

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s | MH Acceptance: 0.80

redist_flip’s constraint includes a relatively weak compactness constraint by default because simulating compact maps is far more efficient and completely non-compact maps are not super useful for most purposes.

You can override this by making a blank redist_constr object

cons <- redist_constr(iowa_map)

Then, you can pass this to redist_flip.

tidy_sims_no_comp <- redist_flip(iowa_map, nsims = 100, constraints = cons)
#> 
#> ── redist_flip() ───────────────────────────────────────────────────────────────
#> 
#> ── Automated Redistricting Simulation Using Markov Chain Monte Carlo ──
#> ℹ Preprocessing data.
#> ℹ Starting swMH().
#> 
■                                  1% | ETA: 0s

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s | MH Acceptance: 0.97

redist_flip outputs a redist_plans object.

class(tidy_sims)
#> [1] "redist_plans" "tbl_df"       "tbl"          "data.frame"

To extract the plans, use get_plans_matrix().

plans <- get_plans_matrix(tidy_sims)

Alternatively, you can directly use functions on the redist_plans object. For example, if we want to measure the competitiveness of each plan:

tidy_sims <- tidy_sims %>% 
  mutate(competitiveness = compet_talisman(pl(), iowa_map, rvote = rep_08, dvote = dem_08))

tidy_sims %>% 
  ggplot(aes(x = competitiveness)) +
  geom_density() + 
  theme_bw()

For more information on using redist_plans objects, see the Get Started page.

Back to top

Diagnostic Plots

When using the MCMC algorithms, there are various useful diagnostic plots. The redist.diagplot function creates familiar plots by converting numeric entries into mcmc objects to use with coda.

We use the dissimilarity index in Massey and Denton 1988 as a summary statistic for the following examples. This can be computed with seg_dissim. In this case, we create a Republican dissimilarity index. We can work with two examples, the first is a single vector of the segregation index, while the second is a list of vectors, with one vector for each chain.

seg <- by_plan(seg_dissim(tidy_sims, iowa_map, rep_08, pop))

The first three plots only need a single index.

• Autocorrelation Plot

redist.diagplot(seg, plot = "autocorr")

• Density Plot

redist.diagplot(seg, plot = "densplot")

• Mean Plot

redist.diagplot(seg, plot = "mean")

As examples for the next two plots, we can use the example above which ran 4 chains. This is the same index, but computed for each chain.

seg_chains <- lapply(1:nchains, function(i) { 
    seg_dissim(flip_chains[[i]], iowa_map, rep_08, pop)
})

• Trace Plot

redist.diagplot(sumstat = seg_chains, plot = "trace")

• Gelman Rubin Plot

redist.diagplot(sumstat = seg_chains, plot = 'gelmanrubin')

Back to top

Tuning Flip Constraints

When using the flip algorithm, the most important and difficult step is setting the right constraint weights. While there may be some general pieces of advice for doing so, no advice can replace working with your data. The bottom line is that every data set is a bit different. What works for one state’s redistricting process, with the data specific to that state at that time may not transfer to another state or municipality or school district. The general process of finding what works might be very similar, but getting the right set of constraint weights and other parameters will vary immensely. Even starting from a different plan within the same time and place can change the weights that perform best. Like most things, the key to tuning flip is patience. Going for a full scale simulation without testing some parameter configurations is likely an inefficient use of time and computing power.

The following highlights some advice on how to tune flip to make it work for your particular redistricting problem. For the advice, we’ll use the following example:

data(iowa)
iowa_map <- redist_map(iowa, existing_plan = cd_2010, pop_tol = 0.02, total_pop = pop)

cons <- redist_constr(iowa_map) %>% 
    add_constr_edges_rem(0.5) %>% 
    add_constr_pop_dev(100)

sims <- redist_flip(map = iowa_map,  nsims = 100)
#> 
#> ── redist_flip() ───────────────────────────────────────────────────────────────
#> 
#> ── Automated Redistricting Simulation Using Markov Chain Monte Carlo ──
#> ℹ Preprocessing data.
#> ℹ Starting swMH().
#> 
■                                  1% | ETA: 0s

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s | MH Acceptance: 0.55

mean(sims$mhdecisions, na.rm = TRUE)
#> [1] 0.55

sims_new <- redist_flip(map = iowa_map, nsims = 100, constraints = cons, 
                        eprob = 0.10, lambda = 2, verbose = FALSE)
mean(sims_new$mhdecisions, na.rm = TRUE)
#> [1] 0.46

dists <- redist.distances(plans = get_plans_matrix(sims))$Hamming
dists_new <- redist.distances(plans = get_plans_matrix(sims_new))$Hamming
adj_dists <- rep(NA_integer_, 100)
adj_dists_new <- rep(NA_integer_, 100)
for(i in 1:100){
  adj_dists[i] <- dists[i, i + 1]
  adj_dists_new[i] <- dists_new[i, i + 1]
}
tibble(Hamming = c(adj_dists, adj_dists_new), 
       `lambda/eprob` = c(rep('0/0.05', 100), rep('2/0.10', 100))) %>% 
  ggplot() + 
  geom_density(aes(x = Hamming, color = `lambda/eprob`)) + 
  theme_bw()

sims <- sims %>% mutate(par = plan_parity(map = iowa_map))

sims <- sims %>% filter(par <= 0.01)

cons <- redist_constr(iowa_map) %>% 
    add_constr_edges_rem(0.25) %>% 
    add_constr_pop_dev(50) %>% 
    add_constr_compet(10, rvote = rep_08, dvote = dem_08) %>% 
    add_constr_splits(10, admin = region)

sims <- redist_flip(iowa_map, 100, constraints = cons)
#> 
#> ── redist_flip() ───────────────────────────────────────────────────────────────
#> 
#> ── Automated Redistricting Simulation Using Markov Chain Monte Carlo ──
#> ℹ Preprocessing data.
#> ℹ Starting swMH().
#> 
■                                  1% | ETA: 0s

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s | MH Acceptance: 0.75

cons <- cons <- redist_constr(iowa_map) %>% 
    add_constr_edges_rem(1.5) %>% 
    add_constr_pop_dev(100) %>% 
    add_constr_compet(40, rvote = rep_08, dvote = dem_08) %>% 
    add_constr_splits(20, admin = region)

sims <- redist_flip(iowa_map, 100, constraints = cons)
#> 
#> ── redist_flip() ───────────────────────────────────────────────────────────────
#> 
#> ── Automated Redistricting Simulation Using Markov Chain Monte Carlo ──
#> ℹ Preprocessing data.
#> ℹ Starting swMH().
#> 
■■■■■■■■■■■■■■■■■                 55% | ETA:  0s | MH Acceptance: 0.34

■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■■  100% | ETA:  0s | MH Acceptance: 0.19

summary(sims$constraint_edges_removed, na.rm = TRUE)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
#>   33.00   36.00   36.50   38.18   41.00   46.00       4

summary(sims$constraint_pop_dev, na.rm = TRUE)
#>     Min.  1st Qu.   Median     Mean  3rd Qu.     Max.     NA's 
#> 0.000007 0.000064 0.000129 0.000168 0.000215 0.000633        4

Some Final Thoughts

Redistricting simulation is very much statistics rather than hard science. When working with flip, or any redistricting sampler, there will be a component that resembles art. Each important variable needs to be included, but getting every variable to the correct target space is not necessarily easy. In general, it may be best to start with one or two constraints and slowly add them to the model. This can help ensure that one single constraint doesn’t dominate the entire process.

When starting off, it’s never a bad idea to run a single simulation to make sure that everything works. If it doesn’t do what you’re expecting, that’s much better than waiting for 1,000,000 iterations to run. If that works, try 100 or 1000. Only once you’ve seen that it’s moving and appears to be moving in reasonable directions should you try for those large numbers of simulations. Remember that running 1,000,000 steps of flip with completely useless parameters is not a very good use of time or computing power.

And finally, when in doubt, it never hurts to run a few extra simulations. Once you know that the code is working, it shouldn’t cost much at all to run just a few extra iterations or a few simulations from new starting points. If the results agree with your prior findings, that’s more support for them. If they disagree, then you know what could be wrong and can run even more additional simulations to figure out what’s right!

Back to top