I have pasted my full code below.
I have been using verbose=TRUE in inla, but can't locate the output. So I am attaching my full code here, which reads the data.
The shape file from US Census (cb_2016_us_county_500k.shp) is too big to attach. But can provide directly if needed.
Any advice is appreciated.
library(sp)
library(rgeos)
library(rgdal)
library(maptools)
library(dplyr)
library(leaflet)
library(scales)
require(maptools)
# Read Data
dat <- read.csv(url, stringsAsFactors = FALSE)
# Colnames tolower
names(dat) <- tolower(names(dat))
dat$countyname <- tolower(dat$countyname)
# Wide data set, subset only what we need.
county_dat <- subset(dat, measureid == "296",
select = c("reportyear","countyfips","statename", "countyname", "value", "unitname")) %>%
subset(reportyear==2011, select = c("countyfips", "value"))
# Rename columns to make for a clean df merge later.
colnames(county_dat) <- c("GEOID", "airqlty")
# Have to add leading zeos to any FIPS code that's less than 5 digits long to get a good match.
# I'm cheating by using C code. sprintf will work as well.
county_dat$GEOID <- formatC(county_dat$GEOID, width = 5, format = "d", flag = "0")
### End data prep
## read shape file
##us.poly<-readShapePoly(system.file("cb_2016_us_county_500k.shp", package="spdep"))
us.poly<-readOGR("cb_2016_us_county_500k.shp")
# Remove Alaska(2), Hawaii(15), Puerto Rico (72), Guam (66), Virgin Islands (78), American Samoa (60)
# Mariana Islands (69), Micronesia (64), Marshall Islands (68), Palau (70), Minor Islands (74)
#us.poly <- us.poly[!us.poly$STATEFP %in% c("02", "15", "72", "66", "78", "60", "69",
# "64", "68", "70", "74"),]
# Make sure other outling islands are removed.
#us.poly <- us.poly[!us.poly$STATEFP %in% c("81", "84", "86", "87", "89", "71", "76",
# "95", "79"),]
# create inla weight matrices
list.queen<-poly2nb(us.poly, queen=TRUE)
cards<-card(list.queen)
tr<- which(cards==0)
# Merge spatial df with downloade ddata.
sub.poly<- us.poly[-tr,]
leafmap2 <- merge(sub.poly, county_dat, by=c("GEOID"), all=FALSE)
list.queen2<-poly2nb(leafmap2, queen=TRUE)
set.ZeroPolicyOption(TRUE)
W<-nb2listw(list.queen2, style="W", zero.policy=TRUE)
#W1<-nb2mat(list.queen2, style="W", zero.policy=TRUE)
WW<-listw2mat(W)
W2<-as(WW, "Matrix")
# graph data
popup_dat2 <- paste0("<strong>County: </strong>",
leafmap2$NAME,
"<br><strong>Value: </strong>",
leafmap2$airqlty)
pal2 <- colorQuantile("YlOrRd", NULL, n = 20)
leaflet(data = leafmap2) %>% addTiles() %>%
addPolygons(fillColor = ~pal(airqlty),
fillOpacity = 0.8,
color = "#BDBDC3",
weight = 1,
popup = popup_dat2)
plot(W, coordinates(leafmap2))
## INLA
library(INLA)
# setup data
D <- data.frame(leafmap2)
leafdata<-transform(D, GEOID = as.numeric(GEOID))
leafdata <- leafdata[ , c("GEOID", "airqlty")]
ID<-seq(1:NROW(leafdata))
leafdata <- cbind(leafdata, ID)
E <-rep(100, NROW(leafdata))
formula <- airqlty ~ 1 + f(ID, model="bym", graph = W2, hyper =
list(prec.unstruct = list(prior="loggamma", param=c(1,0.01)),
prec.spatial = list(prior="loggamma", param=c(1,0.001 ))))
mod<-inla(formula, family="poisson", E=E, data=leafdata, verbose=TRUE)