# Area-level data
library(RColorBrewer)

library(rgdal) 
library(spdep)
library(classInt)

# Import in the shapefile
setwd("C://Users//SmithJ1//Box Sync//MEI//Teaching//Spatial Epi Course//Lecture 5 Analysis of spatial clustering//Lab 1//")
scot_LL <- readOGR(".", "scot")
proj4string(scot_LL) <- CRS("+proj=longlat +datum=WGS84") 

# Import in the data set
setwd("C://Users//SmithJ1//Box Sync//MEI//Teaching//Spatial Epi Course//JLS - Point processes//Data//Scottish lip cancer//")

library(maptools)
scot_dat <- read.csv("scotland 2.csv", header=T, , stringsAsFactors=FALSE)


names(scot_dat) <- c("District", "Observed", "Expected", "PcAFF",
                       "Latitude", "Longitude")
row.names(scot_dat) <- formatC(scot_dat$District, width = 2,
                                 flag = "0")
ID <- formatC(scot_LL$ID, width = 2, flag = "0")
scot_LL1 <- spChFIDs(scot_LL, ID)
scot_LL <- spCbind(scot_LL1, scot_dat[match(ID, row.names(scot_dat)), ])
plot(scot_LL)

# Look at maps of observed and expected counts
nclr <- 6
plotclr <- brewer.pal(nclr, "BuPu")
brks <- c(0, 5, 10, 25, 50, 100)

O_CI <- classIntervals(scot_LL$Observed, style = "fixed", fixedBreaks = brks)
E_CI <- classIntervals(scot_LL$Expected, style = "fixed", fixedBreaks = brks)
O_colcode <- findColours(O_CI, plotclr)
E_colcode <- findColours(E_CI, plotclr)

par(mfrow=c(1,2))
plot(scot_LL, col=O_colcode, main="Observed")
plot(scot_LL, col=E_colcode, main="Expected")

coords<-coordinates(scot_LL)

# Contiguity neighbors - all that share a boundary point
Scot_nb_Cont <- poly2nb(scot_LL)
par(mfrow=c(1,1))
plot(scot_LL)
plot(Scot_nb_Cont, coords, col="green",add=T)


# Distance-based neighbors - returns the k nearest points as neighbors
IDs<-row.names(as(scot_LL, "data.frame"))
Scot_nb_Dist <- knn2nb(knearneigh(coords, k=1, longlat=TRUE),row.names=IDs)
par(mfrow=c(1,1))
plot(scot_LL)
plot(Scot_nb_Dist, coords, col="green",add=T)


# Find the minimum distance at which all areas have a distance-based neighbor
dsts<-unlist(nbdists(Scot_nb_Dist, coords))
summary(dsts)

max_1nn<-max(dsts)
max_1nn

# use dnearneigh to find neighbors with an interpoint distance, with distances settings
# upper and lower distance bounds
Scot_nb_Interpoint0.5 <- dnearneigh(coords, d1=0, d2=0.5*max_1nn, row.names=IDs)
Scot_nb_Interpoint1 <- dnearneigh(coords, d1=0, d2=1*max_1nn, row.names=IDs)
Scot_nb_Interpoint1.5 <- dnearneigh(coords, d1=0, d2=1.5*max_1nn, row.names=IDs)

nb<-list(d1=Scot_nb_Interpoint0.5, d2=Scot_nb_Interpoint1, d3=Scot_nb_Interpoint1.5)

sapply(nb, function(x) is.symmetric.nb(x, verbose=F, force=T))
sapply(nb, function(x) n.comp.nb(x)$nc)

par(mfrow=c(1,3))
plot(scot_LL)
plot(Scot_nb_Interpoint0.5, coords, col="green",add=T)
plot(scot_LL)
plot(Scot_nb_Interpoint1, coords, col="green",add=T)
plot(scot_LL)
plot(Scot_nb_Interpoint1.5, coords, col="green",add=T)

# What do you notice about the number of distance-based neighbors in the maps? 
# NB: due to the large nearest neighbor distances in the north of the study area, 
# the distance criterion will have to be set so that many areas have many neighbors

# Assign weights
scot_W<-nb2listw(Scot_nb_Cont)
scot_W

# Test for spatial autorcorrelation using moran's i
moran<-moran.test(scot_LL$Observed, listw=scot_W)
 
library(pgirmess) # Requires spdep but is much simplier and user-friendly. It does not test the significance
# by a randomization test but uses a normal approximetly. \


pgi.cor <- correlog(coords=coordinates(scot_LL), z=scot_LL$Observed, method="Moran", nbclass=10)
# coordys = xy corredinates, z= vector of values at each location and nbclass = the number of bins
plot(pgi.cor) # statistically significant values (p<0.05) are plotted in red
pgi.cor

# what are the limitations of this approach? How could this approach be improved?