redist: Simulation Methods for Legislative Redistricting

Documented in redist.ipw

###########################################
## Author: Ben Fifield
## Institution: Princeton University
## Date Created: 2015/02/04
## Date Modified: 2015/03/09
## Purpose: R wrapper to run swMH() code (non-mpi)
###########################################

#' Inverse probability reweighting for MCMC Redistricting
#'
#' \code{redist.ipw} properly weights and resamples simulated redistricting plans
#' so that the set of simulated plans resemble a random sample from the
#' underlying distribution. \code{redist.ipw} is used to correct the sample when
#' population parity, geographic compactness, or other constraints are
#' implemented.
#'
#' @param plans An object of class `redist_plans` from `redist_flip()`.
#' @param resampleconstraint The constraint implemented in the simulations: one
#' of "pop", "compact", "segregation", or "similar".
#' @param targetbeta The target value of the constraint.
#' @param targetpop The desired level of population parity. \code{targetpop} =
#' 0.01 means that the desired distance from population parity is 1%. The
#' default is \code{NULL}.
#' @param temper A flag for whether simulated tempering was used to improve the
#' mixing of the Markov Chain. The default is \code{1}.
#'
#' @details This function allows users to resample redistricting plans using
#' inverse probability weighting techniques described in Rubin (1987). This
#' techniques reweights and resamples redistricting plans so that the resulting
#' sample is representative of a random sample from the uniform distribution.
#'
#' @return \code{redist.ipw} returns an object of class "redist". The object
#' \code{redist} is a list that contains the following components (the
#' inclusion of some components is dependent on whether tempering
#' techniques are used):
#' \item{plans}{Matrix of congressional district assignments generated by the
#' algorithm. Each row corresponds to a geographic unit, and each column
#' corresponds to a simulation.}
#' \item{distance_parity}{Vector containing the maximum distance from parity for
#' a particular simulated redistricting plan.}
#' \item{mhdecisions}{A vector specifying whether a proposed redistricting plan
#' was accepted (1) or rejected (0) in a given iteration.}
#' \item{mhprob}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm.}
#' \item{pparam}{A vector containing the draw of the \code{p} parameter for each
#' simulation, which dictates the number of swaps attempted.}
#' \item{constraint_pop}{A vector containing the value of the population
#' constraint for each accepted redistricting plan.}
#' \item{constraint_compact}{A vector containing the value of the compactness
#' constraint for each accepted redistricting plan.}
#' \item{constraint_segregation}{A vector containing the value of the
#' segregation constraint for each accepted redistricting plan.}
#' \item{constraint_similar}{A vector containing the value of the similarity
#' constraint for each accepted redistricting plan.}
#' \item{constraint_vra}{A vector containing the value of the
#' vra constraint for each accepted redistricting plan.}
#' \item{constraint_partisan}{A vector containing the value of the
#' partisan constraint for each accepted redistricting plan.}
#' \item{constraint_minority}{A vector containing the value of the
#' minority constraint for each accepted redistricting plan.}
#' \item{constraint_hinge}{A vector containing the value of the
#' hinge constraint for each accepted redistricting plan.}
#' \item{constraint_qps}{A vector containing the value of the
#' QPS constraint for each accepted redistricting plan.}
#' \item{beta_sequence}{A vector containing the value of beta for each iteration
#' of the algorithm. Returned when tempering is being used.}
#' \item{mhdecisions_beta}{A vector specifying whether a proposed beta value was
#' accepted (1) or rejected (0) in a given iteration of the algorithm. Returned
#' when tempering is being used.}
#' \item{mhprob_beta}{A vector containing the Metropolis-Hastings acceptance
#' probability for each iteration of the algorithm. Returned when tempering
#' is being used.}
#'
#' @references Fifield, Benjamin, Michael Higgins, Kosuke Imai and Alexander
#' Tarr. (2016) "A New Automated Redistricting Simulator Using Markov Chain
#' Monte Carlo." Working Paper.
#' Available at \url{http://imai.princeton.edu/research/files/redist.pdf}.
#'
#' Rubin, Donald. (1987) "Comment: A Noniterative Sampling/Importance Resampling
#' Alternative to the Data Augmentation Algorithm for Creating a Few Imputations
#' when Fractions of Missing Information are Modest: the SIR Algorithm."
#' Journal of the American Statistical Association.
#'
#' @examples
#' \donttest{
#' data(iowa)
#' map_ia <- redist_map(iowa, existing_plan = cd_2010, pop_tol = 0.01)
#' cons <- redist_constr(map_ia)
#' cons <- add_constr_pop_dev(cons, strength = 5.4)
#' alg <- redist_flip(map_ia, nsims = 500, constraints = cons)
#'
#' alg_ipw <- redist.ipw(plans = alg,
#'     resampleconstraint = "pop_dev",
#'     targetbeta = 1,
#'     targetpop = 0.05)
#' }
#'
#' @concept post
#' @export
redist.ipw <- function(plans,
                       resampleconstraint = c("pop_dev", "edges_removed",
                           "segregation", "status_quo"),
                       targetbeta,
                       targetpop = NULL,
                       temper = 0) {

    ## Warnings:
    if (missing(plans) | !inherits(plans, "redist_plans")) {
        cli_abort("Please provide {.arg plans} as a {.cls redist_plans}.")
    }

    plans_ref <- subset_ref(plans)
    plans <- subset_sampled(plans)

    if (length(resampleconstraint) != 1) {
        cli_abort("We currently only support one resamplingconstraint at a time.")
    }
    if (!(resampleconstraint %in% c("pop_dev", "edges_removed", "segregation", "status_quo"))) {
        cli_abort("We do not provide support for that constraint at this time")
    }
    if (missing(targetbeta)) {
        cli_abort("Please specify the target beta value")
    }

    ## Get indices drawn under target beta if tempering
    if (temper == 1) {
        indbeta <- which(plans$beta_sequence == targetbeta)
    } else {
        indbeta <- seq_len(ncol(get_plans_matrix(plans)))
    }

    ## Get indices of draws that meet target population
    if (!is.null(targetpop)) {
        indpop <- which(plans$distance_parity <= targetpop)
    } else {
        indpop <- seq_len(ncol(get_plans_matrix(plans)))
    }

    ## Get intersection of indices
    inds <- intersect(indpop, indbeta)
    ## Construct weights
    psi <- plans[[paste0("constraint_", resampleconstraint)]][inds]
    weights <- 1/exp(targetbeta*psi)

    ## Resample indices
    inds <- sample(inds, length(inds), replace = TRUE, prob = weights)
    ndists <- max(plans$district)
    indx <- unlist(lapply(inds, function(x) {seq(ndists*(x - 1) + 1, ndists*x, by = 1)}))

    ## Subset the entire list
    plans %>% slice(indx)
}

redist.warmup.chain <- function(algout, warmup = 1) {
    if (warmup <= 0) {
        return(algout)
    }
    inds <- seq_len(warmup)
    algout_new <- vector(mode = "list", length = length(algout))
    for (i in seq_along(algout)) {

        ## Subset the matrix first, then the vectors
        if (i == 1) {
            algout_new[[i]] <- algout[[i]][, -inds]
        } else if (length(algout[[i]]) == 1) {
            algout_new[[i]] <- algout[[i]]
        } else if (all(names(algout[[i]]) == "adj")) {
            algout_new[[i]] <- algout[[i]]
        } else {
            algout_new[[i]] <- algout[[i]][-inds]
        }

    }
    names(algout_new) <- names(algout)
    class(algout_new) <- "redist"
    algout_new
}


redist.thin.chain <- function(algout, thin = 100) {
    if (thin <= 1) {
        return(algout)
    }

    inds <- seq(1, ncol(algout$plans), by = thin)
    algout_new <- vector(mode = "list", length = length(algout))
    for (i in seq_along(algout)) {

        ## Subset the matrix first, then the vectors
        if (is.matrix(algout[[i]])) {
            algout_new[[i]] <- algout[[i]][, inds]
        } else if (length(algout[[i]]) == 1) {
            algout_new[[i]] <- algout[[i]]
        } else if (!is.null(names(algout[[i]])) & all(names(algout[[i]]) == "adj")) {
            algout_new[[i]] <- algout[[i]]
        } else {
            algout_new[[i]] <- algout[[i]][inds]
        }

    }
    names(algout_new) <- names(algout)
    class(algout_new) <- "redist"
    algout_new
}