# # Barry Edwards # Replication file for Formulating Voting Rights Act Remedies to Address Current Conditions # # Simulating Effect of Varying VRA Standards on African American Voter Success # Simulate result of different cut-offs on district race composition # This is a fairly complex script. It uses nested "for" loops to simulate elections # within districts, within states, at varying district composition standards. It may # not reproduce Figure 6(a) exactly as it appears in the text because it is a simulation. # # Remove all objects just to be safe # rm(list=ls(all=TRUE)) # # # ### state-by-state look stateData <- read.csv("http://www.poliscidata.com/replication/vra/AfAmRepSimulationData.csv") attach(stateData) # names(stateData) # fix(stateData) # [1] "state" "stNumber" "allPopSources" # [4] "totPop" "overseasPop" "nReps2010" # [7] "hispanicPop" "blackPop" "avgDistSize" #[10] "section5" "neRegion" "stateFamIncome" #[13] "stateDemSupport" "popOver18" "blackOver18" #[16] "hispOver18" "blackVoterProp" "latinoVoterProp" #[19] "stateVoterProp" "avgVotersDist" "interceptChangeAfAm" #[22] "sdInterceptChangeAfAm" "slopeChangeAfAm" "sdSlopeChangeAfAm" #[25] "interceptChangeLatino" "sdInterceptChangeLatino" "slopeChangeLatino" #[28] "sdSlopeChangeLatino" "afAmBreakEven" "latinoBreakEven" # [2] "stateAbbr" # "afAmReps2010" ####################### based on analysis of 2002-2010 data # fixed effects common to all the states # bugs: node mean sd MC error 2.5% median 97.5% start sample commonIntercept <- -8.414 # 1.564 0.1168 -11.95 -8.235 -5.861 10000 10002 commonSlope <- 34.69 # 5.057 0.3881 25.76 34.21 46.05 10000 10002 openSeatCoef <- -2.888 # 1.095 0.05005 -5.255 -2.822 -0.8924 10000 10002 whiteIncCoef <- -7.662 # 1.528 0.0941 -10.94 -7.544 -5.038 10000 10002 hvapCoef <- -6.68 # 4.179 0.3125 -17.46 -6.103 -0.482 10000 10002 # what if more Af Am Repubs elected? # this are approximations designed to test effects # commonIntercept <- commonIntercept + 5 # commonSlope <- commonSlope - 15 ########################### declare variables used in simulation i <- 1 totalExpectedAfAmReps <- NULL totalExpectedLatinoReps <- NULL targetProportion <- NULL minAfAmProportion <- NULL minPopNeededPerDistrict <- NULL afAmDistricts <- NULL afAmRemainder <- NULL afAmOtherDist <- NULL afAmTargetXB <- NULL probTargetAfAmDist <- NULL expectAfAmRepsTarget <- NULL remainderAfAmXB <- NULL probRemainderAfAmDist <- NULL expectAfAmRepsRemainder <- NULL totalAfAmReps <- NULL totalExpectedAfAmRepsSection5only <- NULL totalExpectedAfAmRepsNEonly <- NULL afAmMeanConstant <- NULL latinoVAPtargetDist <- NULL latinoVAPotherDist <- NULL blackVotersProportion <- NULL otherVotersProportion <- NULL blackVotersNeeded <- NULL targetDistOtherPopulation <- NULL targetDistOtherVoters <- NULL targetDistLatinoVoters <- NULL blackPopTargetDistrict <- NULL simulatedAfAmRepsTarget <- NULL simulatedAfAmRepsRemainder <- NULL totalSimulated <- NULL simulationIterations <- 5000 # change number of iterations seatOpensRandomly <- .1 probAfAmDist <- NULL simulatedMean <- NULL simulatedSD <- NULL lower95quantile <- NULL upper95quantile <- NULL avgDistVoters <- NULL # i <- 400 # j <- 33 # north carolina # i <- 394 # at issue in NC litigation while(i <= 1000) { simulatedAfAmRepsTotal <- rep(0,simulationIterations) simulatedAfAmRepsSection5 <- rep(0,simulationIterations) j <- 1 while (j <= 50) { # determine how many target districts can be made targetProportion[j] <- i/1000 minAfAmProportion[j] <- blackOver18[j]/popOver18[j] if(targetProportion[j] < minAfAmProportion[j]) { targetProportion[j] <- minAfAmProportion[j] } afAmProportion <- blackOver18[j]/popOver18[j] geoDensityParameter <- 1 - (1/(1 - afAmProportion))*(1/(1 - afAmProportion))*(targetProportion[j] - afAmProportion)*(targetProportion[j] - afAmProportion) # alternative specification: linear relationship # geoDensityParameter <- (-1/(1 - afAmProportion))*targetProportion[j] + (1/(1 - afAmProportion)) if (targetProportion[j] < afAmProportion) { geoDensityParameter <- 1 } avgDistSize[j] <- totPop[j]/nReps2010[j] avgDistVoters[j] <- popOver18[j]/nReps2010[j] blackVotersProportion[j] <- blackOver18[j]/blackPop[j] otherVotersProportion[j] <- (popOver18[j] - blackOver18[j]) / popOver18[j] blackVotersNeeded[j] <- (avgDistSize[j]*otherVotersProportion[j]) / ( ((1 - targetProportion[j])/targetProportion[j]) + (1/blackVotersProportion[j])*otherVotersProportion[j]) blackVotersNeeded[j] <- blackVotersNeeded[j] blackPopTargetDistrict[j] <- blackVotersNeeded[j]*(1/blackVotersProportion[j]) # blackPopTargetDistrict[j] <- avgDistSize[j]*targetProportion[j] afAmDistricts[j] <- floor(blackPop[j]*geoDensityParameter/blackPopTargetDistrict[j]) if (afAmDistricts[j] >= nReps2010[j]) # can reach target in all districts { afAmDistricts[j] <- nReps2010[j] afAmRemainder[j] <- 0 afAmOtherDist[j] <- 0 } if (afAmDistricts[j] < nReps2010[j]) # can reach target in less than all districts { afAmRemainder[j] <- blackOver18[j] - (afAmDistricts[j] * blackVotersNeeded[j]) afAmOtherDist[j] <- afAmRemainder[j] / ((nReps2010[j] - afAmDistricts[j])* avgDistVoters[j]) } if (afAmDistricts[j] == 0) # cannot reach target in any district { blackPopTargetDistrict[j] <- 0 afAmOtherDist[j] <- blackOver18[j]/(avgDistVoters[j]*nReps2010[j]) } targetDistOtherPopulation[j] <- avgDistSize[j] - (blackPopTargetDistrict[j]) targetDistOtherVoters[j] <- targetDistOtherPopulation[j]* otherVotersProportion[j] targetDistLatinoVoters[j] <- targetDistOtherVoters[j]*(hispOver18[j]/(popOver18[j] - blackOver18[j])) latinoVAPtargetDist[j] <- targetDistLatinoVoters[j] / (blackVotersNeeded[j] + targetDistOtherVoters[j]) latinoVAPotherDist[j] <- (1 - afAmOtherDist[j])*(hispOver18[j]/(popOver18[j] - blackOver18[j])) # here's where I am specifying an open seat, if the status varies, it would need to change here # I could determine the number of seats in each state, how many af am incumbents they have, allocate # incumbents based on current representation, want to include some random open seat factor (roughly 10%) totalAfAmReps[1:simulationIterations] <- 0 CD <- 1 while(CD <= nReps2010[j]) { if (afAmDistricts[j] >= CD) { BVAP <- targetProportion[j] HVAP <- latinoVAPtargetDist[j] # need to fix this } if (afAmDistricts[j] < CD) { BVAP <- afAmOtherDist[j] HVAP <- latinoVAPotherDist[j] } if (afAmReps2010[j] >= CD) { AfAmInc <- 1 whiteInc <- 0 openSeat <- 0 } if (afAmReps2010[j] < CD) { AfAmInc <- 0 whiteInc <- 1 openSeat <- 0 } openSeat <- rbinom(1, 1, seatOpensRandomly) if (openSeat == 1) { AfAmInc <- 0 whiteInc <- 0 openSeat <- 1 } afAmMeanConstant <- commonIntercept + openSeatCoef*openSeat + whiteIncCoef*whiteInc + hvapCoef*HVAP afAmXB <- (commonSlope + slopeChangeAfAm[j])*BVAP + afAmMeanConstant + interceptChangeAfAm[j] probAfAmDist[CD] <- exp(afAmXB )/(1 + exp(afAmXB )) simulatedAfAmReps <- rbinom(simulationIterations, 1, probAfAmDist[CD]) totalAfAmReps <- totalAfAmReps + simulatedAfAmReps CD <- CD + 1 } # simulatedAfAmRepsTotal <- simulatedAfAmRepsTotal + totalAfAmReps simulatedAfAmRepsSection5 <- simulatedAfAmRepsSection5 + (totalAfAmReps*section5[j]) j <- j + 1 } simulatedMean[i] <- mean(simulatedAfAmRepsTotal) simulatedSD[i] <- sd(simulatedAfAmRepsTotal) lower95quantile[i] <- quantile(simulatedAfAmRepsTotal,probs=0.025) upper95quantile[i] <- quantile(simulatedAfAmRepsTotal,probs=0.975) # totalExpectedAfAmReps[i] <- sum(totalAfAmReps) totalExpectedAfAmReps[i] <- mean(simulatedAfAmRepsTotal) totalExpectedAfAmRepsSection5only[i] <- mean(simulatedAfAmRepsSection5) # totalExpectedAfAmRepsNEonly[i] <- sum(totalAfAmReps*stateData$neRegion) i <- i + 1 } # write.table(stateData,file="H:/voting rights/Myfile70.csv",sep=",",row.names=TRUE) #pnorm((45.021-48.79504)/ 2.698489) #pnorm(1.96) # comparing 46.8% and 50% standard # at 47% mean = 73.7252, sd = 2.452772 # at 50% mean = 67.5198, sd = 1.765398 # simulatedMean[470] # simulatedMean[500] # simulatedSD[470] # simulatedSD[500] simAt47 <- rnorm(1000,simulatedMean[468],simulatedSD[468]) simAt50 <- rnorm(1000,simulatedMean[500],simulatedSD[500]) mean(simAt47 >= simAt50) # how often 47% greater than nationwide proportionality mean(simAt47 >= 0.1121725*435) mean(simAt50 >= 0.1121725*435) ####################################################################### create the plot ########### uncomment to write graphics files to disk rather than display in R # setwd('H:/voting rights/hi res figures/') # png(file="afAmSimulationResultsFigure7panelA.png",width=3,height=3, units="in", pointsize=10, res=800) par(mar=c(3.2,2.8,1.7,0.1), mgp=c(2,.6,.5), family = "serif") par(xpd=FALSE) plot(x="",y="",xlim=c(1,1000),ylim=c(0,70),xlab="Target African American VAP (%)",ylab="Expected Seats in Congress", axes=F,cex.lab=1.0,main="(a) African American Representation",cex.main=1.0) abline(h= 0.1121725*435,lty=2,lwd=1) abline(h=0.168954*95,lty=3,col="gray30") # nationwide polygon(x=c(seq(1,1000,by=1),seq(1000,1,by=-1)),y=c(totalExpectedAfAmReps-qnorm(.975)*simulatedSD,rev(totalExpectedAfAmReps+qnorm(.975)*simulatedSD)),col="gray90",border=F) polygon(x=c(seq(1,1000,by=1),seq(1000,1,by=-1)),y=c(totalExpectedAfAmReps-qnorm(.95)*simulatedSD,rev(totalExpectedAfAmReps+qnorm(.95)*simulatedSD)),col="gray80",border=F) lines(x=seq(1,1000,by=1),y=totalExpectedAfAmReps,lwd=1) # section 5 states lines(x=seq(1,1000,by=1),y=totalExpectedAfAmRepsSection5only,lwd=1,lty=1,col="gray30") axis(side=1,seq(0,1000,by=100),labels=seq(0,100,by=10),cex.axis=.7,line=0) axis(side=2,seq(0,80,by=10),las=2,cex.axis=.7,line=0) box() xAxis <- seq(1,1000,by=1) xAxis[totalExpectedAfAmReps== max(totalExpectedAfAmReps)] ########### uncomment to write graphics files to disk rather than display in R # dev.off() #### just the legend now ########### uncomment to write graphics files to disk rather than display in R # setwd('H:/voting rights/hi res figures/') # png(file="simulationFigure7legendOnly.png",width=5.5,height=.6, units="in", pointsize=10, res=800) par(mar=c(3.2,0.1,0.1,0.1), mgp=c(2,.6,.5), family = "serif") plot(x="",y="",xlim=c(1,1000),ylim=c(0,60),xlab="",ylab="",axes=F,main="") par(xpd=TRUE) legend(x=0, y=0, border=c(F,F,T,T,F,F), fill=c('white','white','gray80','gray90','white','white'), lty=c(1,2,0,0,1,3), lwd=c(1,1,0,0,1,1), ncol=3, cex=.7, col=c('black','black','white','white','gray30','gray30'), legend=c('Expected Seats Nationwide', 'National Proportionality','90% Confidence Interval','95% Confidence Interval','Expected Seats, Section 5 States','Section 5 State Proportionality')) ########### uncomment to write graphics files to disk rather than display in R # dev.off()