I've got a timeseries model I am fitting to a city dataset with about 45 zones. The data are daily, stratified by Zone, Age and Sex. Following on from learning about spatiallly correlated errors I want to see if the Standard Error on the estimated \(\beta_{1}\) from the timeseries model is affected.

I think the simplest option is to use the spatial lag model, which can be fitted with just adding a term that is the average of the set of each Zone's neighbours outcome level on each day. For this I need to find the list of each region's neighbours. Then I'll use this to assign each zone/day/age group their neighbours values and then collapse that to get their daily means.

# we have access to a classic dataset for studying spatial dependence
# in the spdep package
if(!require(spdep))    install.packages(spdep); require(spdep)     
if(!require(rgdal))    install.packages(rgdal); require(rgdal) 
if(!require(maptools)) install.packages(maptools); require(maptools) 
if(!require(maps))     install.packages(maps); require(maps) 
fn <- system.file("etc/shapes/eire.shp", package="spdep")[1]
prj <- CRS("+proj=utm +zone=30 +units=km")
eire <- readShapeSpatial(fn, ID="names", proj4string=prj)
str(eire)
# reproject into a better coordinate system
eire <- spTransform(eire, CRS("+proj=longlat +datum=WGS84"))
# check out the attributes
head(eire@data)

nb <- poly2nb(eire)
str(nb)
#List of 26
nb[[1]]
#[1]  9 10 11 25 26
# So this returns the set of index values for each area's neighbours
# I'd prefer to read their names
eire[['names']][1]
# > [1] Carlow
# so therefore the neighbours of area 1 "Carlow" are in the first
# element of the list
eire[['names']][nb[[1]]]
# > [1] Kildare  Kilkenny Laoghis  Wexford  Wicklow

################################################################
# name:plot these
png("images/Fig1.png")
plot(eire)
plot(nb, coordinates(eire), add=TRUE, pch=".", lwd=2)
map.scale(ratio = F)
box()
dev.off()

I THINK I actually want this as a dataframe so I can merge it with the master table of outcome data.

################################################################
# name:adjacency_df
adjacency_df <- function(NB, shp, zone_id)
  {
    adjacencydf <- as.data.frame(matrix(NA, nrow = 0, ncol = 2))
    for(i in 1:length(NB))
    {
      if(length(shp[[zone_id]][NB[[i]]]) == 0) next
      adjacencydf <- rbind(
                           adjacencydf,
                           cbind(
                                 as.character(shp[[zone_id]][i]),
                                 as.character(shp[[zone_id]][NB[[i]]])
                                 )
                           )
    }
    return(adjacencydf)
  }

################################################################
# name:adjacency_df
adj <- adjacency_df(NB = nb, shp = eire, zone_id = 'names')
adj