## Note for the NYleukemia example, 4 census tracts were completely surrounded by another
## unique census tract; when applying the Bayesian cluster detection model in
## \code{\link{bayes.cluster}}, we merge them with the surrounding census tracts yielding
## \code{n=277} areas.
## Load data and convert coordinate system from latitude/longitude to grid
data(NYleukemia)
map <- NYleukemia$spatial.polygon
population <- NYleukemia$data$population
cases <- NYleukemia$data$cases
centroids <- latlong2grid(NYleukemia$geo[, 2:3])
## Identify the 4 census tract to be merged into their surrounding census tracts.
remove <- NYleukemia$surrounded
add <- NYleukemia$surrounding
## Merge population and case counts and geographical objects accordingly
population[add] <- population[add] + population[remove]
population <- population[-remove]
cases[add] <- cases[add] + cases[remove]
cases <- cases[-remove]
map <- SpatialPolygons(map@polygons[-remove], proj4string=CRS("+proj=longlat"))
centroids <- centroids[-remove, ]
## Expected numbers
E <- expected(population, cases, 1)
## Set Parameters
max.prop <- 0.15
k <- 0.00005
shape <- c(2976.3, 2.31); rate <- c(2977.3, 1.31)
J <- 7
pi0 <- 0.95
n.sim.imp <- 0.5*10^4
n.sim.prior <- 0.5*10^4
n.sim.post <- 0.5*10^5
## (Uncomment first) Compute output
# output <- bayes.cluster(E, cases, population, centroids, map, max.prop, k,
# shape, rate, J, pi0, n.sim.imp, n.sim.prior, n.sim.post)
# plotmap(output$prior.map$high.area, map)
# plotmap(output$post.map$high.area, map)
# plotmap(output$post.map$RR.est.area, map, log=TRUE)
# barplot(output$pj.y, names.arg=0:J, xlab="j", ylab="P(j|y)")
Run the code above in your browser using DataLab