R/simulate_election.R

In i3pack: Incentives for Inter- And Intra-Party Electoral Competition

#' Function to simulate a full election in a single district
#'
#' The function runs a complete election in a single district, 
#' using the simulation framework described in detail in
#' Chapter 4 of Crisp et al. 2024.  
#'
#' @param voters Optional vector of voter positions in 1d ideological space.
#' @param parties Optional vector of party positions in 1d ideological space. Maximum of 10 parties allowed. 
#' @param cands Optional matrix with three columns: candidate 1d ideological position, unique numerical candidate ID, and positive numerical candidate valence  
#' @param nominated Optional data.frame with five variables: `rank` (candidate ranking in the party list); `candidate` (numeric candidate ID); `pos` (1d ideological position of candidate); `list` (numeric list ID; equal to 1, unless parties are allowed to have multiple lists); `party` (numeric party ID). 
#' @param nvoters Number of voters; defaults to 3,000.
#' @param nparties Number of parties; defaults to 5; maximum allowable: 10.
#' @param nvotes Number of votes per voter; defaults to 1. Can also take on special values `0` (which then is internally replaced by the district magnitude) and `-1` (which is then internally replaced by 1 fewer vote than the district magnitude). 
#' @param M District magnitude; defaults to 5.
#' @param rank_cand Boolean: should candidates be ranked on the party list? Defaults to `TRUE`.
#' @param strategic Boolean: do parties and voters behave strategically? Defaults to `TRUE`.
#' @param strategic_error Numeric probability with which strategic actors fail to choose the optimal alternative.
#' @param who_ranks Character actor who arranges party lists, one of `parties`, `voters`,`none`; defaults to `parties`. 
#' @param gamma_val Numeric weight assigned to the valence component of voters' utility function.
#' @param gamma_rank Numeric weight assigned to the candidate ranking on the party list when computing the voter's utility. 
#' @param elec_fun_name Name of function implementing electoral system formula. 
#' @param ballot_type Character string indicating type of ballot, one of `open`,`closed`, or `flexible`; defaults to `open`. 
#' @param primary Boolean: should a primary election be conducted? Defaults to `FALSE`.
#' @param two_round Boolean: should a second election round be conducted? Defaults to `FALSE`.
#' @param pool_level Character level at which votes are pooled, one of `party_list` (or sub-party list),`party`, or `candidate`. Defaults to `party_list`
#' @param ranked_vote Boolean: Do voters cast a ranked vote? Defaults to `FALSE`.
#' @param free_vote Boolean: If voters can cast multiple votes, can the be for candidates in different parties? Defaults to `FALSE`. 
#' @param max_cand Numeric maximum number of candidates running in a party list; defaults to 0, which is internally interpreted as the district magnitude.
#' @param threshold Numerical legal electoral threshold; defaults to 0 (i.e., no threshold).
#' @param lists_per_party Integer allowed number of lists per party; defaults to 1.
#' @param seed Random number generator seed; defaults to 123.
#' @param elec_results_only Boolean: Should function return ancillary information on election, or just election results? Defaults to `FALSE`.
#' @param multiplier Numeric factor by which to multiply the votes cast by voters with the same ideological position; defaults to 1. 
#' @param system_name Character name of electoral system used, one of 'AV', 'BC', 'STV', 'MNTV', 'LV', 'PR', or 'SMDP'
#' @param ... Additional arguments passed to `elec_fun_name`.
#' 
#'
#' @return data.frame with the following variables (if `elec_results_only=FALSE`,
#' otherwise, data.frame with candidate id's, positions, valences, votes obtained,
#' and whether they won a seat or not):
#' \describe{
#' \item{gamma_val}{See `Usage` above}
#' \item{epsilon}{Maximum acceptable ideological distance used in voters' utility function}
#' \item{hetero}{Measure of elected candidate heterogeneity}
#' \item{pers}{Average valence of elected candidates}
#' \item{lsq}{Least Squares measure of disproportionality}
#' \item{enp_v}{Effective number of electoral parties}
#' \item{enp_s}{Effective number of legislative parties}
#' \item{avg_dist}{Average distance between elected candidates and voters}
#' \item{var_elect}{Variance of ideological positions of elected candidates}
#' \item{avg_vote_util}{Average utility of voters w.r.t. candidates they voted for}
#' \item{avg_elect_util}{Average utility of voters w.r.t. elected candidates}
#' \item{sample_parties}{Parties that initially could have entered the election}
#' \item{ran_parties}{Parties that decided to enter the election}
#' }
#' 
#' @examples
#' # Simulate a PR (D'Hondt) election with 3 parties, 5 candidates per party,
#' # 100 voters, and a district magnitude of 2, allowing for strategic voting
#' 
#' simulate_election(parties = c(-1, 0, 1),
#'                   nvoters = 100,
#'                   M = 2,
#'                   strategic = TRUE,
#'                   elec_fun_name = "dhondt",
#'                   system_name = "PR")
#' @export
#' @import dplyr
#' @importFrom magrittr "%>%"
#' @importFrom rlang "!!" sym .data
#' @importFrom stats rnorm rbeta var
#' @importFrom tidyr replace_na
#' @importFrom utils head
#'
simulate_election <- function(voters = NULL,
                     parties = NULL,
                     cands = NULL,
                     nominated = NULL,
                     nvoters = 3000,
                     nparties = 5,
                     nvotes = 1,
                     M = 5,
                     rank_cand = TRUE,
                     strategic = TRUE,
                     strategic_error = 0.05, 
                     who_ranks = c("parties", "voters","none"),
                     gamma_val = NULL,
                     gamma_rank = 1, # How much weight does list rank have?
                     elec_fun_name = "dhondt",
                     ballot_type = "open",
                     primary = FALSE,
                     two_round = FALSE,
                     pool_level = c("party_list","party", "candidate"),
                     ranked_vote = FALSE,
                     free_vote = FALSE,
                     max_cand = 0,
                     threshold = 0,
                     lists_per_party = 1,
                     seed = 123,
                     elec_results_only = FALSE,
                     multiplier = 1, 
                     system_name,
                     ...)
{
  set.seed(seed)
  ## Check allowed values
  who_ranks <- match.arg(who_ranks)
  pool_level <- match.arg(pool_level)
  if(strategic & elec_fun_name %in% c("stv","a_v")){
    stop("Strategic voting not defined for STV or AV")
  }
  if((rank_cand == FALSE) & (who_ranks != "none")){
    stop("If rank_cand is FALSE, who_ranks must be 'none'.")
  }
  ## Set max_cand equal to M
  if(max_cand == 0){
    max_cand <- M
  }
  ## Ditto with nvotes
  if(nvotes == 0){
    nvotes <- M
  } else if (nvotes == -1){
    nvotes <- M-1
  }
  
  ## More system checks
  if((elec_fun_name == "stv") & ((ballot_type != "open") | (nvotes == 1) | (ranked_vote == FALSE) | (rank_cand == TRUE))){
    stop("STV only implemented for `open' ballot type, nvotes > 1, ranked_vote TRUE, and rank_cand FALSE.")
  }
  if((elec_fun_name == "a_v") & ((free_vote != TRUE) | (nvotes == 1) | (ranked_vote == FALSE) | (M > 1) | (lists_per_party > 1) | (max_cand > 1))){
    stop("AV only implemented for `free' votes, nvotes > 1, ranked_vote TRUE, M = 1, lists_per_party = 1, and max_cand = 1.")
  }
  if(nparties > 10 | length(parties) > 10){
    stop("Maximum number of parties is 10.")
  }
  
  ## Set electoral system formula
  elec_fun <- get(elec_fun_name)
  
  ## Define pool level as the default pool level
  pool_level <- pool_level[1]
  
  ## Sample voters, and party leaders 
  if(is.null(voters)){
    voters <- rnorm(nvoters)
  }
  if(is.null(parties)){
    party_locs <- seq(-2, 2, length.out = nparties)
    sel_parties <- sample(length(party_locs), nparties)
    parties <- party_locs[sel_parties]
    sel_parties <- sel_parties[order(parties)]
    parties <- sort(parties) 
    names(parties) <- 1:nparties
  } else {
    nparties <- length(parties)
    sel_parties <- order(parties)
    parties <- sort(parties) 
    if(is.null(names(parties))){
      names(parties) <- 1:nparties
    }
  }
  
  ## Select preliminary pool of candidates by party/list
  if(is.null(nominated)){
    nominated <- nominating(parties, lists_per_party, rank_cand, max_cand)
  }
  ## Create cand matrix and add valence
  if(is.null(cands)){
    cands <- nominated[,c("pos","candidate")]
    ncands <- nrow(cands)
    cands <- cbind(cands, val = rbeta(ncands, 2, 4))
  }
  ## Sample epsilon tolerance for squared distance
  spat_comp <- outer(cands[,1],
                     voters,
                     FUN = function(x,y){(x - y)^2})
  epsilon <- sample(spat_comp, 1) 
  
  ## Sample ≠ gammas values for all voters
  if(is.null(gamma_val)){
    gamma_val <- sample(c(1, 5, 10), nvoters, replace=TRUE) 
  }
  
  
  ## Actual nomination of candidates and form lists
  if(!primary){ ##Without primaries
    
    if(strategic){
      ## Run Poll for strategic entry and voting
      polled_votes <- sample(1:nvoters, nvoters * 0.1, replace = FALSE)
      poll_res <- simulate_election(voters = voters[polled_votes],
                           parties = parties,
                           nparties = length(parties),
                           cands = cands,
                           nominated = nominated,
                           nvotes = nvotes,
                           M = M,
                           rank_cand = rank_cand,
                           strategic = FALSE,
                           strategic_error = 0,
                           who_ranks = who_ranks,
                           gamma_val = gamma_val, 
                           gamma_rank = 0, 
                           elec_fun_name = elec_fun_name,
                           ballot_type = ballot_type,
                           primary = FALSE,
                           two_round = FALSE,
                           pool_level = pool_level,
                           ranked_vote = ranked_vote,
                           free_vote = free_vote,
                           max_cand = max_cand,
                           threshold = 0,
                           lists_per_party = lists_per_party,
                           seed = seed,
                           elec_results_only = TRUE,
                           multiplier = 1, 
                           system_name = system_name,
                           ...)
      party_info <- poll_res %>%
        mutate(party = factor(.data$party, levels = 1:nparties)) %>%
        group_by(.data$party, .drop=FALSE) %>%
        summarize(tot_votes = sum(votes),
                  tot_seats = sum(.data$winner), .groups='drop_last') %>%
        mutate(loc = parties) %>%
        arrange(desc(.data$tot_votes), desc(.data$tot_seats), abs(.data$loc)) %>%
        mutate(higher_seats = c(1, 1, head(.data$tot_seats, -2))) %>%
        mutate(innit = ifelse(.data$tot_seats >= 1 | .data$higher_seats != 0, 1, 0))
      party_info$innit_fail <- party_info$innit
      party_info$innit_fail[party_info$innit_fail == 0] <- sample(c(0, 1), sum(party_info$innit_fail == 0), replace=TRUE, prob = c(1-strategic_error, strategic_error))
      
      entry_parties <- party_info$party[party_info$innit_fail != 0]
      safe_parties <- party_info$party[party_info$innit != 0]
    } else {
      entry_parties <- safe_parties <-  1:nparties
    }
    nominated <- nominated[order(nominated[,"candidate"]),]
    if (who_ranks == "none"){
      nominated[,"rank"] <- ifelse(nominated[,"rank"] > 0, 1, 0)
    }
    nominated <- cbind(nominated, "safe_party"= nominated[,"party"] %in% safe_parties)
    
  } else { ## With primaries
    if(lists_per_party > 1){
      stop("Primaries not implemented when lists_per_party > 1.")
    }
    
    if(strategic){
      ## Run Poll for strategic entry and voting
      polled_votes <- sample(1:nvoters, nvoters * 0.1, replace = FALSE)
      poll_res <- simulate_election(voters = voters[polled_votes],
                           parties = parties,
                           cands = cands,
                           nominated = nominated,
                           nvotes = nvotes,
                           nparties = nparties,
                           M = M,
                           rank_cand = rank_cand,
                           strategic = FALSE,
                           strategic_error = 0,
                           who_ranks = who_ranks,
                           gamma_val = gamma_val[polled_votes], 
                           gamma_rank = 0, 
                           elec_fun_name = elec_fun_name,
                           ballot_type = ballot_type,
                           primary = FALSE,
                           two_round = FALSE,
                           pool_level = pool_level,
                           ranked_vote = ranked_vote,
                           free_vote = free_vote,
                           max_cand = max_cand,
                           threshold = 0,
                           lists_per_party = lists_per_party,
                           seed = seed,
                           elec_results_only = TRUE,
                           multiplier = multiplier,
                           system_name = system_name,
                           ...)
      
      party_info <- poll_res %>%
        mutate(party = factor(.data$party, levels = 1:nparties)) %>%
        group_by(.data$party, .drop=FALSE) %>%
        summarize(tot_votes = sum(votes),
                  tot_seats = sum(.data$winner), .groups='drop_last') %>%
        mutate(loc = parties) %>%
        arrange(desc(.data$tot_votes), desc(.data$tot_seats), abs(.data$loc)) %>%
        mutate(higher_seats = c(1, 1, head(.data$tot_seats, -2))) %>%
        mutate(innit = ifelse(.data$tot_seats >= 1 | .data$higher_seats != 0, 1, 0))
      party_info$innit_fail <- party_info$innit
      party_info$innit_fail[party_info$innit_fail == 0] <- sample(c(0, 1), sum(party_info$innit_fail == 0), replace=TRUE, prob = c(1-strategic_error, strategic_error))
      
      entry_parties <- party_info$party[party_info$innit_fail != 0]
      safe_parties <- party_info$party[party_info$innit != 0]
    } else {
      entry_parties <- safe_parties <-  1:nparties
    }
    
    nominated_tmp <- nominated %>% 
      mutate(rank = 1,
             safe_party = .data$party %in% safe_parties, 
             val = cands[,"val"])
    df_cands_tmp <- nominated_tmp
    
    
    ## Run primary with entrants  
    votes_primary_res <- voting(voters,
                                df_cands_tmp,
                                1,
                                gamma_val,
                                gamma_rank,
                                epsilon,
                                free = TRUE,
                                closed_primary = FALSE,
                                strategic = strategic,
                                strategic_error = strategic_error,
                                party_pos = parties) 
    votes_primary <- votes_primary_res$votes
    
    ##Rank candidates according to primary votes
    nominated_tmp <- do.call(rbind,by(nominated_tmp,
                                      list(nominated_tmp[,"party"]),
                                      function(x, votes, n, rank_cand){
                                        cand_ind <- x$candidate
                                        rank_tmp <- rank(-votes[cand_ind], ties.method = "random")
                                        x$rank <- ifelse(rank_tmp <= n, rank_tmp, 0)
                                        if(!rank_cand){
                                          x$rank <- ifelse(x$rank > 0, 1, 0)
                                        }
                                        return(x)
                                      }, votes = table(factor(votes_primary, levels = 1:ncands)), n = max_cand, rank_cand = rank_cand))
    if(who_ranks == "parties"){
      for(i in seq_along(parties)){
        nominated_tmp[nominated_tmp[,"party"]==i, "rank"] <- rank(abs(nominated_tmp[nominated_tmp[,"party"]==i, "pos"] - parties[i]), ties.method = "random")
      }
    } else if (who_ranks == "none"){
      nominated_tmp[,"rank"] <- ifelse(nominated_tmp[,"rank"] > 0, 1, 0)
    }
    nominated <- nominated_tmp[order(nominated_tmp[,"candidate"]),]
  } ## end with primaries
  
  df_cands <- as.data.frame(nominated)
  df_cands$party_list <- with(df_cands, factor(paste(party,list, sep="_")))
  df_cands$val <- cands[,3]
  
  df_cands <- df_cands %>%
    filter(rank > 0)
  
  if(nvotes > 1){
    nvotes <- min(nvotes, sum((df_cands$rank)>0))
  }
  
  ## Vote (votes is nvotes x voters)
  votes_res <- voting(voters,
                      df_cands,
                      nvotes,
                      gamma_val,
                      gamma_rank,
                      epsilon,
                      free_vote,
                      strategic = strategic,
                      strategic_error = strategic_error)
  
  votes <- votes_res$votes
  vote_utils <- votes_res$max_utils
  
  # vector to keep winners in the first round if M > 1 
  winners_fr <- c()
  if(two_round){
    cand_votes_fr <- table(votes) # first round
    fraction = (1/2)/M
    fraction_out = 0.025
    
    if(M == 1){
      if(all(cand_votes_fr < (sum(cand_votes_fr) * 0.5))){
        top_two_cand <- names(cand_votes_fr)[max_n(cand_votes_fr, 2)]
        
        df_cands[!(df_cands[,"candidate"] %in% top_two_cand), "rank"] <- 0
        votes_res <- voting(voters,
                            df_cands,
                            nvotes,
                            gamma_val,
                            gamma_rank,
                            epsilon,
                            free_vote,
                            strategic = FALSE)
        votes <- votes_res$votes
      }
      
    }else {
      
      if(all(cand_votes_fr < (sum(cand_votes_fr) * fraction))){
        winners_fr <- cand_votes_fr[cand_votes_fr >= (sum(cand_votes_fr) * fraction)]
        # Update M
        M <- M - length(winners_fr)
        # Candidates for the second round
        cand_votes_keep <- cand_votes_fr[cand_votes_fr < (sum(cand_votes_fr) * fraction)]
        top_cand <- names(cand_votes_keep)[cand_votes_keep > (sum(cand_votes_keep) * fraction_out)]
        df_cands[!(df_cands[,"candidate"] %in% top_cand), "rank"] <- 0
        votes_res <- voting(voters,
                            df_cands,
                            nvotes,
                            gamma_val,
                            gamma_rank,
                            epsilon,
                            free_vote,
                            strategic = FALSE)
        votes <- votes_res$votes      
      }
    }
  }  
  
  
  df_cands <- df_cands %>%
    dplyr::arrange(.data$party, .data$list, .data$rank)
  
  ## Get final distance to choice
  if(!elec_results_only){
    distances <- array(NA, nvoters)
    for(i in 1:nvoters){
      distances[i] <- ifelse(votes[1,i] > 0, abs(df_cands[df_cands$candidate==votes[1,i],"pos"] - voters[i]), NA)
    }
  }
  
  ## Tally the votes by candidate
  if(elec_fun_name %in% c('a_v', 'stv', 'bc')){
    cand_votes <- table(votes[1, ])  
  }else if(two_round){
    cand_votes <- cand_votes_fr
  }else{
    cand_votes <- table(c(votes))  
  }
  cand_votes_df <- data.frame(candidate = as.numeric(names(cand_votes)),
                              votes = c(cand_votes)) 
  suppressMessages(df_cands <- df_cands %>% 
                     left_join(cand_votes_df) %>%
                     tidyr::replace_na(list(votes = 0)))
  
  
  ##Re-sort lists for systems in which personal votes matter
  if(ballot_type %in% c("flexible", "open") & (ranked_vote == FALSE) & (pool_level != "candidate")){
    by_party_list <- split.data.frame(df_cands, df_cands$party_list)
    ranked_by_pl <- do.call(rbind,lapply(by_party_list,
                                         function(x, nvoters, nvotes, M){
                                           if(ballot_type == "flexible"){
                                             quota <- nvoters * nvotes * 0.03
                                             if(any(x$votes >= quota)){
                                               cand_over <- (x$votes >= quota)
                                               n_over <- sum(cand_over)
                                               x$rank[cand_over][order(x[cand_over,]$votes, decreasing = TRUE)] <- 1:n_over
                                               x$rank[!cand_over] <- x$rank[!cand_over] + n_over
                                             }
                                           }
                                           if((ballot_type == "open") & (elec_fun_name != "stv")){
                                             x$rank <- rank(-x$votes, ties.method = "random")#
                                           } 
                                           return(x)
                                         }, nvoters = nvoters, nvotes = nvotes, M = M))
  } else {
    ranked_by_pl <- df_cands
  }
  ## Pool votes to desired level
  if(ranked_vote == FALSE){
    tmp <- ranked_by_pl %>%
      group_by(!!sym(pool_level)) %>%
      summarize(tot_votes = sum(.data$votes), .groups='drop_last')
    pooled_votes <- array(unlist(tmp[,2]), c(1, nrow(tmp)), dimnames = list(NULL,unlist(tmp[,1])))
  } else {
    if(pool_level != "candidate"){
      stop("Ranked vote is only defined when pool_level is 'candidate'")
    }
    pooled_votes <- votes
  }
  
  
  ## Distribute seats at vote-pooling level
  if(!(system_name %in% c('AV', 'BC', "STV", 'MNTV'))){
    pooled_votes <- (pooled_votes + 1) * multiplier
  }
  
  if((elec_fun_name == "plurality") & (M > 1) & (ballot_type == "closed") & two_round == FALSE){ ##PBV
    tmp_seats <- elec_fun(v = pooled_votes, m = 1, threshold = threshold, ...)
    tmp_seats <- tmp_seats * M
  } else {
    tmp_seats <- elec_fun(v = pooled_votes, m = M,  threshold = threshold, ...)
  }
  if(length(winners_fr) > 0){ # only for TR with M>1
    tmp_seats <- c(tmp_seats, winners_fr)
  }
  
  pooled_seats <- data.frame(seats = c(tmp_seats),
                             key = as.factor(names(tmp_seats)))
  names(pooled_seats)[2] <- pool_level
  
  if(pool_level %in% c("candidate","party")){
    pooled_seats[,pool_level] <- as.numeric(levels(pooled_seats[,pool_level]))[pooled_seats[,pool_level]]
  }
  
  # If two_round == TRUE and M > 1
  if(two_round == TRUE & M > 1){
    ranked_by_pl[,pool_level] <- as.numeric(as.character(ranked_by_pl[,pool_level]))
  }
  suppressMessages(ranked_by_pl <- ranked_by_pl %>%
                     left_join(pooled_seats) %>%
                     tidyr::replace_na(list(seats = 0)))
  
  
  ## Elect candidates
  ranked_by_pl <- ranked_by_pl %>%
    group_by(!!sym(pool_level)) %>%
    mutate(winner = ifelse((.data$rank <= .data$seats)&(.data$rank > 0), 1, 0)) %>%
    ungroup() 
  # If two_round == TRUE and M > 1
  if(two_round == TRUE &  M > 1){
    ranked_by_pl[,pool_level] <- as.factor(as.character(ranked_by_pl[,pool_level]))  
  }
  
  if(elec_results_only){
    return(ranked_by_pl)
  } else {
    if(lists_per_party > 1){
      res_by_party_list <- ranked_by_pl %>%
        group_by(.data$party, .data$party_list) %>%
        summarise(
          across(c(.data$votes, .data$seats), ~ max(.x, na.rm = TRUE)),
          .groups = "drop")
      res_by_party <- res_by_party_list %>% 
        group_by(.data$party) %>%
        summarise(
          across(c(.data$all_votes, .data$all_seats), ~ sum(.x)),
          .groups = "drop")      
      res_by_party$all_votes[!(res_by_party$party %in% as.numeric(as.character(safe_parties)))] <- NA 
      res_by_party$all_seats[!(res_by_party$party %in% as.numeric(as.character(safe_parties)))] <- NA 
    } else if(system_name %in% c('STV', 'MNTV', 'SNTV', 'LV') | (two_round == TRUE & M > 1)){
      res_by_party <- ranked_by_pl %>% 
        group_by(.data$party) %>%
        summarise(
          across(c(.data$ll_votes, .data$all_seats), ~ sum(.x)),
          .groups = "drop")         
      res_by_party$all_votes[!(res_by_party$party %in% as.numeric(as.character(safe_parties)))] <- NA 
      res_by_party$all_seats[!(res_by_party$party %in% as.numeric(as.character(safe_parties)))] <- NA 
      
    } else{
      res_by_party <- ranked_by_pl %>% 
        group_by(.data$party) %>%
        summarise(all_votes = sum(.data$votes),
                  all_seats = max(.data$seats),
                  .groups = "drop") 
      res_by_party$all_votes[!(res_by_party$party %in% as.numeric(as.character(safe_parties)))] <- NA 
      res_by_party$all_seats[!(res_by_party$party %in% as.numeric(as.character(safe_parties)))] <- NA 
    }
    
    res_by_party$all_votes <- (res_by_party$all_votes + 1) * multiplier
    
    lsq_res <- lsq(res_by_party$all_votes[!is.na(res_by_party$all_votes)],
                   res_by_party$all_seats[!is.na(res_by_party$all_seats)])
    enp_v_res <- 1/sum(((res_by_party$all_votes[!is.na(res_by_party$all_votes)])/sum(res_by_party$all_votes[!is.na(res_by_party$all_votes)]))^2)
    enp_s_res <- 1/sum(((res_by_party$all_seats[!is.na(res_by_party$all_seats)])/sum(res_by_party$all_seats[!is.na(res_by_party$all_seats)]))^2)
    
    
    just_win <- dplyr::filter(ranked_by_pl, .data$winner == 1)
    
    ## Compute utilities w.r.t elected cands
    spat_comp <- outer(just_win$pos,
                       voters,
                       FUN = function(x,y){(x - y)^2})
    elec_utils <- apply( array(gamma_val, dim(spat_comp)) * just_win$val - spat_comp, 2, max, na.rm=TRUE)
    
    
    ## Compute and return QIs
    cohesion <- sum(mapply(
      function(party.pos, cand.pos)
      {
        mean((party.pos - unlist(cand.pos))^2)# * length(cand.pos)/M
      },
      party.pos=as.list(parties[unique(just_win$party)]),
      cand.pos=split(just_win$pos, just_win$party)))
    pers <- mean(just_win$val, na.rm=TRUE)
    
    ## Form df of party information
    names(parties) <- paste0("Loc",1:nparties)
    party_vote <- array(NA, nparties, dimnames = list(paste0("Votes", 1:nparties)))
    party_vote[sel_parties[res_by_party$party]] <- res_by_party$all_votes
    party_seat <- array(NA, nparties, dimnames = list(paste0("Seats", 1:nparties)))
    party_seat[sel_parties[res_by_party$party]] <- res_by_party$all_seats
    party_df <- as.data.frame(c(as.list(parties),
                                as.list(party_vote),
                                as.list(party_seat)))
    
    # Get parties that were sampled
    sampledParties <- paste(sort(c(1:nparties)[sel_parties[res_by_party$party]]), collapse = ', ')
    
    # Get entry parties
    entryParties <- paste(sort(sel_parties[as.numeric(as.character(safe_parties))]), collapse = ', ')
    
    
    out <- cbind(data.frame(gamma_val = stat_mode(gamma_val),
                            epsilon = epsilon,
                            hetero = cohesion,
                            pers = pers, 
                            lsq = lsq_res,
                            enp_v = enp_v_res,
                            enp_s = enp_s_res,
                            avg_dist = mean(distances, na.rm = TRUE),
                            var_elect = ifelse(M == 1, 0, var(just_win$pos, na.rm = TRUE)),
                            avg_vote_util = mean(vote_utils[is.finite(vote_utils)], na.rm = TRUE),
                            avg_elect_util = mean(elec_utils[is.finite(elec_utils)], na.rm = TRUE),
                            sample_parties = sampledParties,
                            ran_parties = entryParties),
                 party_df)
    return(out)
  }
}