library(maptools)
library(spdep)

data(eire)
# The 'eire.df' data frame has 26 rows and 9 columns.      #
# The data set includes:                                   #
# eire.polys.utm: a 'polylist' object,                     #
# eire.coords.utm: polygon centroids,                      #
# eire.nb: the original Cliff and Ord binary contiguities. #

### 1. Create nearest neighbors ###
## 1a. k nearest neighbors ##
eire.coords = cbind(eire.coords.utm$V1, eire.coords.utm$V2)

col.knn <- knearneigh(eire.coords, k=4)
# knearneigh(spdep): returns a matrix with the indices of regions belonging to the set of the k nearest neighbors of each other.#

plot(eire.polys.utm, border="grey", forcefill=FALSE) 

k.nb <- knn2nb(col.knn) 
# knn2nb: converts a 'knn' object returned by 'knearneigh' into a neighbors list of class 'nb' with a list of integer vectors #  
# containing neighbor region number ids. #

plot(knn2nb(col.knn), eire.coords, add=TRUE)
title(main="K nearest neighbors, k = 4")

## 1b. nearest distance based neighbors ##

all.linked <- max(unlist(nbdists(k.nb, eire.coords)))

col.nb.0.all <- dnearneigh(eire.coords, 0, all.linked)
# dnearneigh: identifies neighbors of region points by Euclidean distance between lower (greater than) and upper (less than or #
# equal to) bounds, or with lonlat = TRUE, by Great Circle distance in kilometers.                                              #

# another example of dnearneigh, distance <= 3.0: #
# dnb = dnearneigh(eire.coords, 0, 100)                #
# max(unlist(nbdists(dnb,eire.coords)))                #
# [1] 99.96239                                    #

plot(eire.polys.utm, border="grey", forcefill=FALSE)
plot(col.nb.0.all, eire.coords, add=TRUE)
title(main=paste("Distance based neighbors 0-",  format(all.linked), " distance units", sep=""))

### 2. Fit CAR & SAR model ###
## 2a. CAR model ##
# ROADACC: arterial road network accessibility in 1961 #
# OWNCONS: percentage in value terms of gross agricultural output of each county consumed by itself #

mcar = spautolm(OWNCONS ~ ROADACC, data=eire.df, listw=nb2listw(k.nb, style="W"), family="CAR", method="full")
summary(mcar)

# nb2listw:  supplements a neighbors list with spatial weights for the chosen coding scheme. #
# options of style in nb2listw: B is the basic binary coding, W is row standardized (sums over all links to n), #
# C is globally standardized (sums over all links to n) ...                                                     #

# spautolm: CAR & SAR models #

## 2b. SAR model ##
msar1 = spautolm(OWNCONS ~ ROADACC, data=eire.df, listw=nb2listw(k.nb, style="W"), family="SAR", method="full")
summary(msar1)

msar2 = errorsarlm(OWNCONS ~ ROADACC, data=eire.df, nb2listw(k.nb, style="W"), method="eigen")
summary(msar2, correlation=TRUE)
   
### 3. R sample code for plotting choropleth polygons using "maps" package ###
source("http://www.biostat.umn.edu/~brad/data/MinnesotaChoropleth.R")


### 4. Construct and plot objects for the neighborhoods for CH3 Hw Prob 9 ###

data(state)

sat<-read.table("http://www.biostat.umn.edu/~brad/data/state-sat.dat",header=F,skip=15)
# Discard DC and overall average #
sat<-sat[-c(9,52),]

sat[,1]<-state.abb
# 2 and 11 correspond to Alaska and Hawaii, discard them #
sat<-sat[-c(2,11),]
names(sat)<-c("state","verb","math","p")
nb<-read.csv("http://www.biostat.umn.edu/~brad/data/contig-lower48.dat",header=F,skip=16,sep=" ")

# Generate the neighborhood structure for R #
n.state<-dim(nb)[1]

W.nb<-matrix(0,n.state,n.state)

for( i in 1:n.state)
{
    a<- ! is.na(nb[i,])
      W.nb[i,unlist(nb[i,a])]<-1
}

listw<-mat2listw(W.nb)
# Now listw is a neighborhood structure that can be used in R #