Calculate local indices of disproportionality

local_disproportionality.R

## 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)