Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Calculate local indices of disproportionality
## Libraries
library(maptools)
library(rgdal)
library(plyr)
library(rgeos)
library(RColorBrewer)
## Read in the shape
x <- readShapeSpatial("boundaries/westminster_const_region.shp",IDvar="NAME")
## Get centroids
cents <- gCentroid(x,byid=TRUE)@coords
## Get all combinations
combs <- expand.grid(ConstA = rownames(cents),
ConstB = rownames(cents))
combs$dist <- apply(combs,1,function(row){
centA <- cents[row["ConstA"],]
centB <- cents[row["ConstB"],]
xdist <- centA[1] - centB[1]
ydist <- centA[2] - centB[2]
dist <- sqrt(xdist^2 + ydist^2)
return(dist)
})
## For each constituency, select the twenty closest
combs <- ddply(combs,.(ConstA),function(df){
df <- df[order(df$dist),][1:20,]
df
})
## Merge on the right hand column with seat shares and vote totals
## Must convert to Press Association number
name2pa <- read.csv("name2pa.csv",header=T)
combs <- merge(combs,name2pa,
by.x="ConstB",by.y="ShapeFileName")
pnorris2010 <- read.csv("pnorris_const_2010.csv",header=T)
combs <- merge(combs,pnorris2010,
by.x="refno",by.y="RefNo")
gallagher <- function(votes,seats) {
seats <- seats/sum(seats,na.rm=T) * 100
votes <- votes/sum(votes,na.rm=T) * 100
res <- votes - seats
res <- res ^ 2
res <- sum(res)
res <- res / 2
res <- sqrt(res)
res
}
### Calculate the disproportionality by ConstA
ldi <- ddply(combs,.(ConstA),function(df){
voteslice <- df[,c("Convt10","Labvt10","LDvt10","SNPvt10","PCvt10","Greenvt10","BNPvt10","UKIPvt10")]
votes <- colSums(voteslice,na.rm=T)
seats <- apply(voteslice,1,function(x)x==max(x,na.rm=T))
seats <- rowSums(seats,na.rm=T)
gallagher(votes,seats)
})
### Now plot this
names(ldi) <- c("ConstA","LDI")
rownames(ldi) <- ldi$ConstA
ldi.spdf <- SpatialPolygonsDataFrame(x, ldi)
my.ylims <- bbox(ldi.spdf)[2,]
my.ylims[2] <- my.ylims[2] * 10/11
my.xlims <- bbox(ldi.spdf)[1,]
my.xlims[1] <- my.xlims[1] + (my.xlims[2] - my.xlims[1]) * .15
pdf(file="ldi.pdf",width=6,height=8)
trellis.par.set(list(regions = list(col = colorRampPalette(brewer.pal(9, "BuPu"))(100))))
spplot(ldi.spdf,c("LDI"),col="transparent",
ylim = my.ylims,
xlim = my.xlims
)
dev.off()
## Write out the data
ldi <- ldi[order(ldi$LDI),]
write.csv(ldi,"ldi.csv",row.names=FALSE)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.