R codes

This document presents the R codes used to obtain the computational results included in the paper ``Exploring the Effects on the Electoral College of National and Regional Popular Vote Interstate Compact: An Electoral Engineering Perspective’’.

Presentation of SimulElec

We have conceived a program, called SimuElect which simulates the electoral outcomes of mechanisms M1 and M2. The program has so far been conceived for the case of equipopulated districts but the integers K and r can be arbitrary. The main arguments of this program are the following :

K : an integer (an odd) which gives the number of states.
n.simu : the number of simulations.
r : an integer, which gives as intra-size of one state 2r +1,
type : a integer 0 for IC, 1 for IAC and 2 for IAC*.

To make our scripts run faster, we code the main algorithm into C++ which is called from SimuElec() thanks to the R package Rcpp. The C++ code of that function is available at this link and may be copied into a folder whose path is given in the argument path.

SimuElect <- function(x, K, n.simu, r, type, path) {
library("Rcpp")  
sourceCpp(paste(path,"/codesCplusplus.cpp", sep = ""))
res <- f(K, n.simu, r, type)
return(res)
}

To run parallel computing, we write the function SimuElect which allows to send several simulations on several cores via a simple modification of the parameter x.

We present here a small example of \(SimuElec()\) in a non parallel computing framework.

We consider \(K=51\) states :

K <- 51

The intra-size of one state is equal to 210 +1 = 21

r <- 10

The number of replications is equal to 50000 :

n.simu <- 50000

We consider the IC probabilistic model to simulate the choice of the electors :

type <- 0 # IC model
prob.model <- "IC" # IC model

To execute the pogram, we simply use this command :

path <- getwd()
res_IC_r10_K51 <- SimuElect(1, K = K, n.simu = n.simu, r = r, 
                            type = type, path = path)

Values returned by the function

The result of SimuElect is a matrix of size \((K+1) \times 19\), where each row corresponds to a value of \(L\) (the size of the coalition \(L=0,...,K\)). The columns correspond to the following informations :

column 1 : number of times that an election inversion occured in M1 among the n.simu simulations,
column 2 : number of people satisfied inside the coalition by the election result in M1 among the n.simu simulations,
column 3 : number of people satisfied outside the coalition by the election result in M1 among the n.simu simulations,
column 4 : number of total people satisfied by the election result in M1 among the n.simu simulations,
column 5 : number of pivotal people inside the coalition in M1 among the n.simu simulations,
column 6 : number of pivotal people outside the coalition in M1 among the n.simu simulations,
column 7 : number of pivotal people in M1 among the n.simu simulations,
column 8 : number of times that an election inversion occured in M2,
column 9 : number of people satisfied inside the coalition by the election result in M2 among the n.simu simulations,
column 10 : number of people satisfied outside the coalition by the election result in M2 among the n.simu simulations,
column 11 : number of total people satisfied by the election result in M2 among the n.simu simulations,
column 12 : number of pivotal people inside the coalition in M2 among the n.simu simulations,
column 13 : number of pivotal people outside the coalition in M2 among the n.simu simulations,
column 14 : number of pivotal people in M2 among the n.simu simulations,
column 15 : group DO among the pivotal people outside the coalition in M2,
column 16 : group IO among the pivotal people outside the coalition in M2,
column 17 : group BD among the pivotal people outside the coalition in M2,
column 18 : group BI among the pivotal people outside the coalition in M2,
column 19 : group B among the pivotal people outside the coalition in M2.

Normalize the results

The values returned by the function are not given in percent. To obtain the different figures plotted in the article, one needs to normalize the results by dividing it by the number of simulations and eventually the size of the population (either inside, outside or total).

Defining the size of the population

We compute the size of the population (total, outside and inside) depending on \(L\). We will then divide the appropriated columns by one of these sizes to get the probabilities.

  size_pop <- K * (2*r + 1) 
  size_in <- (0:K)*(2*r + 1)
  size_out <- size_pop - size_in

Election inversions

To get the probability of an election inversion, we only divide by the number of simulations:

pct.inversion.M1 <- res_IC_r10_K51[,1]/n.simu
pct.inversion.M2 <- res_IC_r10_K51[,8]/n.simu

Utilities

To get the utilities, we divide by both the number of simulations and the size of the appropriate population (inside, outside or total) :

mean.Ut.in.M1 <- res_IC_r10_K51[,2]/n.simu/size_in
mean.Ut.out.M1 <- res_IC_r10_K51[,3]/n.simu/size_out
mean.Ut.M1 <- res_IC_r10_K51[,4]/n.simu/size_pop

mean.Ut.in.M2 <- res_IC_r10_K51[,9]/n.simu/size_in
mean.Ut.out.M2 <- res_IC_r10_K51[,10]/n.simu/size_out
mean.Ut.M2 <- res_IC_r10_K51[,11]/n.simu/size_pop

Pivotalities

To get the pivotalities, we also divide by both the number of simulations and the size of the appropriate population (inside, outside or total) :

mean.pivot.in.M1 <- res_IC_r10_K51[,5]/n.simu/size_in
mean.pivot.out.M1 <- res_IC_r10_K51[,6]/n.simu/size_out
mean.pivot.M1 <- res_IC_r10_K51[,7]/n.simu/size_pop

mean.pivot.in.M2 <- res_IC_r10_K51[,12]/n.simu/size_in
mean.pivot.out.M2 <- res_IC_r10_K51[,13]/n.simu/size_out
mean.pivot.M2 <- res_IC_r10_K51[,14]/n.simu/size_pop

To get the sharing of pivotality under M2, we divide the different cases by the total number of observed pivotalities under M2 :

mean.pivot.out.DO.M2 <- res_IC_r10_K51[,15]/res_IC_r10_K51[,13]
mean.pivot.out.IO.M2 <- res_IC_r10_K51[,16]/res_IC_r10_K51[,13]
mean.pivot.out.BD.M2 <- res_IC_r10_K51[,17]/res_IC_r10_K51[,13]
mean.pivot.out.BI.M2 <- res_IC_r10_K51[,18]/res_IC_r10_K51[,13]
mean.pivot.out.B.M2 <- res_IC_r10_K51[,19]/res_IC_r10_K51[,13]

Plotting

Finally, we use the following commands to obtain the figures presented in the article :

Election inversions

Under M1

  par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
  plot(0:K, pct.inversion.M1, type = "l", xlab = "Number of states inside the coalition", 
       ylab = paste("Probability of an election inversions under ", prob.model, " (Mechanism 1)", sep=""), 
       lty = 1, lwd = 3, col = "magenta", xaxt="n", las=2, 
       cex.axis = 0.8)
  axis(1)
  par(tcl =  -0.2)
  axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
  grid()
  legend("topright", legend = c("Maximum", "Back to initial probability"), lwd = 3, 
         pch = c(2, 6))
  
  # the max
  points(which.max(pct.inversion.M1)-1, max(pct.inversion.M1), col = "magenta", pch=2, cex=2)
  segments(which.max(pct.inversion.M1)-1, -0.5, which.max(pct.inversion.M1)-1, 
           max(pct.inversion.M1), lty = 2, col="grey")
  # the return to the beginning horizontal
  segments(which(pct.inversion.M1<pct.inversion.M1[1])[2] - 1, -0.5, 
           which(pct.inversion.M1<pct.inversion.M1[1])[2] - 1, 
           pct.inversion.M1[which(pct.inversion.M1<pct.inversion.M1[1])[2]], 
           lty = 2, col="grey")
  points(which(pct.inversion.M1<pct.inversion.M1[1])[2] - 1, 
         pct.inversion.M1[which(pct.inversion.M1<pct.inversion.M1[1])[2]], 
         col = "magenta", pch=6, cex=2)

Under M2

  op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
  plot(0:K, pct.inversion.M2, type = "l", xlab = "Number of states inside the coalition", 
       ylab = paste("Probability of an election inversions under ", prob.model, " (Mechanism 2)", sep=""), 
       lty = 1, lwd = 3, col = "magenta", xaxt="n", las=2, 
       cex.axis = 0.8)
  axis(1)
  par(tcl = -0.2)
  axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
  grid()

Utilities

Under M1

p <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
  plot(0:K, mean.Ut.M1, type = "l", xlab = "Number of states inside the coalition", 
       ylab = paste("Utility under ", prob.model, " (Mechanism 1)", sep=""), 
       lty = 1, lwd = 3, col = "cyan", xaxt = "n", las = 2, cex.axis = 0.8,
       ylim = range(c(mean.Ut.in.M1, mean.Ut.out.M1), na.rm = T))
  axis(1)
  par(tcl =  -0.2)
  axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
  grid()
  lines(0:K, mean.Ut.out.M1, type = "l", lty = 3, lwd = 3, col = "cyan")
  lines(0:K, mean.Ut.in.M1, type = "l", lty = 2, lwd = 3, col = "cyan")
  
  legend("topright", legend = c("Minimum", "Back to initial level", "Total", "Inside", "Outside"), 
         lwd = 3, lty = c(NA, NA, 1:3),
         pch = c(6, 2, NA, NA, NA))
  
  # the min of the total 
  points(which.min(mean.Ut.M1)-1, min(mean.Ut.M1, na.rm=T), col = "cyan", pch=6, cex=2)
  segments(which.min(mean.Ut.M1)-1, -0.5, which.min(mean.Ut.M1)-1, min(mean.Ut.M1, na.rm=T), 
           lty = 2, col="grey")
  # the return to the beginning horizontal
  segments(which(mean.Ut.M1>mean.Ut.M1[1])[1] - 1, -0.5, which(mean.Ut.M1>mean.Ut.M1[1])[1] - 1, 
           mean.Ut.M1[which(mean.Ut.M1>mean.Ut.M1[1])[1]], lty = 2, col="grey")
  points(which(mean.Ut.M1>mean.Ut.M1[1])[1] - 1, 
         mean.Ut.M1[which(mean.Ut.M1>mean.Ut.M1[1])[1]], col = "cyan", pch=2, cex=2)

Under M2

  op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
  plot(0:K, mean.Ut.M2, type = "l", xlab = "Number of states inside the coalition", 
       ylab = paste("Utility under ", prob.model, " (Mechanism 2)", sep=""), 
       lty = 1, lwd = 3, col = "cyan", xaxt = "n", las = 2, cex.axis = 0.8,
       ylim = range(c(mean.Ut.in.M2, mean.Ut.out.M2), na.rm = T))
  axis(1)
  par(tcl =  -0.2)
  axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
  grid()
  lines(0:K, mean.Ut.out.M2, type = "l", lty = 3, lwd = 3, col = "cyan")
  lines(0:K, mean.Ut.in.M2, type = "l", lty = 2, lwd = 3, col = "cyan")
  
  legend("topright", legend = c("Total", "Inside", "Outside"), lwd = 3, lty = 1:3)

Pivotality

Under M1

  op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
  plot(0:K, mean.pivot.in.M1, type = "l", lty = 2, lwd = 3, col = "green", 
       ylim = range(c(mean.pivot.M1, mean.pivot.in.M1, mean.pivot.out.M1), na.rm = T),
       ylab = paste("Probability to be pivot ", " (Mechanism 1)", sep=""), 
       xlab = "Number of states inside the coalition", las=2, 
       xaxt = "n", cex.axis = 0.8)
  axis(1)
  par(tcl =  -0.2)
  axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
  grid()
  lines(0:K, mean.pivot.out.M1, type = "l", lty = 3, lwd = 3, col = "green")
  lines(0:K, mean.pivot.M1, type = "l", lty = 1, lwd = 3, col = "green")
  legend("topright", legend = c("Total", "Inside", "Outside"), lwd = 3, lty = 1:3)

Under M2

  op <- par(mar = c(4, 4, 0.5, 0.5), mgp = c(3, 0.5, 0))
  plot(0:K, mean.pivot.in.M2, type = "l", lty = 2, lwd = 3, col = "green", 
       ylim = range(c(mean.pivot.M2, mean.pivot.in.M2, mean.pivot.out.M2), na.rm = T),
       ylab = paste("Probability to be pivot ", " (Mechanism 2)", sep=""), 
       xlab = "Number of states inside the coalition", las=2, 
       xaxt = "n", cex.axis = 0.8)
  axis(1)
  par(tcl =  -0.2)
  axis(1, at = 1:K, labels = F, lwd.ticks = 0.5)
  grid()
  lines(0:K, mean.pivot.out.M2, type = "l", lty = 3, lwd = 3, col = "green")
  lines(0:K, mean.pivot.M2, type = "l", lty = 1, lwd = 3, col = "green")
  legend("topright", legend = c("Total", "Inside", "Outside"), lwd = 3, lty = 1:3)

Decomposition of the outside pivotality under M2

require("ggplot2")

## Loading required package: ggplot2

data.ggplot <- data.frame(Sharing = c(mean.pivot.out.DO.M2[-(K+1)], 
                                      mean.pivot.out.IO.M2[-(K+1)], 
                                      mean.pivot.out.BD.M2[-(K+1)], 
                                      mean.pivot.out.BI.M2[-(K+1)], 
                                      mean.pivot.out.B.M2[-(K+1)]),
                            Sector = c(rep("DO", K), rep("IO", K), rep("BD", K), 
                                       rep("BI", K), rep("B", K)),
                            L = rep(0:(K-1), 5))
  
  data.ggplot$Sharing[is.na(data.ggplot$Sharing)] <- 0
  myplot <- ggplot(data.ggplot, aes(x = L, y = Sharing, fill = Sector)) + 
    geom_area(alpha = 0.6 , size = 1, colour = "black")
  print(myplot)

Final programs

The figures given in the article have been done for the three probalistic models (“IC”, “IAC” and “IAC*”), for different values of r (1, 10, 100, 1000, 10000, 100000, 1000000) and with a number of replications equal to \(1~000~000~000\). It has been done in a parallel framework by using 40 cores. The computationnal time needed was nearly equivalent to 3 days. The codes are available here : full program.

R codes

Olivier de Mouzon, Thibault Laurent, Michel Lebreton and Dominique Lepelley

Presentation of SimulElec

Values returned by the function

Normalize the results

Defining the size of the population

Election inversions

Utilities

Pivotalities

Plotting

Election inversions

Under M1

Under M2

Utilities

Under M1

Under M2

Pivotality

Under M1

Under M2

Decomposition of the outside pivotality under M2

Final programs