library(maptools)
library(spdep)

nc_file = system.file("shapes/sids.shp", package='maptools')[1]
llCRS = CRS("+proj=longlat +datum=NAD27")
nc = readShapePoly(nc_file, ID="FIPSNO", proj4string=llCRS)
rn = sapply(slot(nc, "polygons"), function(x) slot(x, "ID"))
gal_file = system.file("etc/weights/ncCR85.gal", package='spdep')[1]
ncCR85 = read.gal(gal_file, region.id=rn)

# Figure of standardized mortatlity ratio -- using
# independent Poisson with rate independent of county

nc$Observed = nc$SID74
nc$Population = nc$BIR74
r = sum(nc$Observed) / sum(nc$Population)
nc$Expected = r * nc$Population
nc$SMR = nc$Observed / nc$Expected
spplot(nc, "SMR")

# Or, with glm

exp2 = fitted(glm(nc$Observed ~ nc$Population - 1, family=poisson("identity")))
print(exp2)
print(nc$Expected)

# Empirical Bayes estimate of Clayton and Kaldor

library(DCluster)
eb = empbaysmooth(nc$Observed, nc$Expected)
print(c(eb$nu, eb$alpha))
nc$EBPG = exp(eb$smthrr)
spplot(nc, 'EBPG')

# Lognormal model

library(DCluster)
ebln = lognormalEB(nc$Observed, nc$Expected)
print(c(ebln$phi, ebln$sigma2))
nc$EBLN = exp(ebln$smthrr)
spplot(nc, 'EBLN')