###########################################################################################################
#  Set up your workplace

###########################################################################################################

# Clean workspace
rm(list = ls())

# Load paths to file directories (adapt as needed)
root<-"C://"
path<-"Users//SmithJ1//Box Sync//MEI//Teaching//Spatial Epi Course//Lecture 5 Analysis of spatial clustering//"


# Attach libraries for visualisation
library(rgdal)
library(raster)
library(ggplot2)

# Data management
library(sp)

# Moran's I and spatial dependencies
library(spdep) # Spatial Dependence: Weighting Schemes, Statistics and Models
library(ape) # Analyses of Phylogenetics and Evolution
library(pgirmess) # Data Analysis in Ecology

# Attach library to convert between data types used by different packages (ppp and spatial points)
library(maptools) # Must have rgeos installed as well

###########################################################################################################
#  Part I:  Global measures of clustering for point-level data

###########################################################################################################

# Data input and visualisation
###########################################################################################################
setwd(paste(root, path, "Lab 1", sep=""))

# Read in prevalence data from the Gambia
Gambia <- readRDS("GambiaID.RData")

# Read in shapefile for the Gambia
GMB_Adm2<-getData('GADM', country='GMB', level=2)                                          # "raster" package
proj4string(GMB_Adm2) <- CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ')

# Calculate prevalence from individual malaria data by cluster ID
GambiaPrev<-as.data.frame(cbind(unique(Gambia$x), unique(Gambia$y), tapply(Gambia$pos, Gambia$ID, mean)))
colnames(GambiaPrev)<-c("x", "y", "Prev")

# Visualise data
map <- fortify(GMB_Adm2, region="NAME_2")                                                  #"ggplot2" package
ggplot() +
  geom_path(data = map, aes(x = long, y = lat, group = group)) +
  geom_point(data = GambiaPrev, aes(x = x, y = y, size = Prev), color = "red") 

# Convert to sp
coordinates(GambiaPrev)<-~x+y # convert to SPDF                                            # "sp" package

# adding Coordinate Referent Sys.
proj4string(GambiaPrev) <- CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ')


###           ADD CODE    #1                     ##
# Plot the data in google maps by prevalence as you learned in the earlier practical

# Q1: Visually assess whether the prevalence of malaria clusters within the study area.

###########################################################################################################
# Moran's i and correlograms
# NB can do this with multiple packages
###########################################################################################################

# Approach 1: Calculate moran's I using a distance based matrix in "Ape"

# First generate a distance matrix
gambia.dists <- as.matrix(dist(cbind(GambiaPrev$x, GambiaPrev$y)))

# Take the inverse of the matrix values so that closer values have a larger weight and vs vs
gambia.dists.inv <- 1/gambia.dists
diag(gambia.dists.inv) <- 0 # replace the diagonal values with zero

# Computes Moran's I autocorrelation coefficient of x giving a matrix of weights (here based on distance) 
Moran.I(GambiaPrev$Prev, gambia.dists.inv)                                                   # "ape" package


# Create a correlogram to explore moran's i over different spatial lags
       # "pgirmess" requires spdep (which also has correlogram options) but is much simplier and user-friendly. 

# Calculate the maximum distance between points
maxDist<-max(dist(cbind(GambiaPrev$x, GambiaPrev$y)))
maxDist

xy=cbind(GambiaPrev$x, GambiaPrev$y)
pgi.cor <- correlog(coords=xy, z=GambiaPrev$Prev, method="Moran", nbclass=10)               # "pgirmess" package
# 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

#Q2: Based on the correlogram, over what spatial lags are there evidence for spatial autocorrelation? 
#    Is this clustering positive or negative?


# Calculate moran's I using a binary distance matrix

#Q3: What is a sensible distance to classify points as neighbors/not? 
#Considerations might include the scale of analysis and the distribution of points

coords<-coordinates(xy)
class(coords)
IDs<-row.names(as.data.frame(coords))

# In this approach, we chose a distance d such that pairs of points with distances less than 
# d are neighbors and those further apart are not. 

Neigh_nb<-knn2nb(knearneigh(coords, k=1, longlat = TRUE), row.names=IDs)                 # "spdep" package
#assigns at least on neighbor to each and calculates the distances between
dsts<-unlist(nbdists(Neigh_nb,coords)) 
# returns the distance between nearest neighbors for each point
summary(dsts)

max_1nn<-max(dsts)
max_1nn # maximum distance to provide at least one neighbor to each point

# Create neighbor structures
Neigh_kd1<-dnearneigh(coords,d1=0, d2=max_1nn, row.names=IDs)   # neighbors within maximum distance
Neigh_kd2<-dnearneigh(coords,d1=0, d2=2*max_1nn, row.names=IDs) # neighbors within 2X maximum distance

nb_1<-list(d1=Neigh_kd1, d2=Neigh_kd2)
sapply(nb_1, function(x) is.symmetric.nb(x, verbose=F, force=T))
# Checks for symmetry (i.e. if i is a neighbor of j, then j is a neighbor of i). Does not always hold for k-nearest neighbours
sapply(nb_1, function(x) n.comp.nb(x)$nc)
# Number of disjoint connected subgraphs

# Plot neighbors comparing the two distances
par(mfrow=c(1,1))
plot(xy, pch=16)
plot(Neigh_kd1, coords, col="green",add=T)
plot(xy, pch=16)
plot(Neigh_kd2, coords,col="green", add=T)

#STEP 2
#assign weights; 
weights<-nb2listw(Neigh_kd1, style="B")   # simplest binary weights, using minimum distance for one neighbor
weights

##STEP 3
moran.test(GambiaPrev$Prev, listw=weights)  #using row standardised weights

# Simulate the test statistic using random permutations of GambiaPrev$Prev so that the values are randomlhy
# assigned to locations and the statistic is computed nsim times.
# NB: if you have a lot of observations can repeat the simulation many times without repetition
set.seed(1234)
bperm<-moran.mc(GambiaPrev$Prev, listw=weights,nsim=999)
bperm
#statistic = 0.582, observed rank = 1000, p-value = 0.001

# Plot simulated test statistics
par(mfrow=c(1,1))
hist(bperm$res, freq=T, breaks=20, xlab="Simulated Moran's I")
abline(v=0.58284, col="red")

# Q4: Assess the evidence for spatial autocorrelation from the Moran's I test


###########################################################################################################
#  Part II:  Local measures of clustering for point-level data

###########################################################################################################
# Local Moran's i

# Calculate the local spatial statistic Moran's I around each point based on the spatial weights object (binary based on at least one neighbor)
I <-localmoran(GambiaPrev$Prev, weights)                                                     # "spdep" package

# Print LISA for each point
Coef<-printCoefmat(data.frame(I[IDs,], row.names=row.names(coords),
                              check.names=FALSE))

# Plot the spatial data against its spatially lagged values (weighted mean of neighbors)                         
nci<-moran.plot(GambiaPrev$Prev, listw=weights, 
                xlab="Prevalence", ylab="Spatially lagged prev", labels=F, pch=16, col="grey")


require(GeoXp)

# Map points that are local outliers in the plot
infl<-apply(nci$is.inf,1,any) # which are statistically significant
sum(infl==T)#7 true (11% - more than would expect by chance)
x<-GambiaPrev$Prev
lhx<-cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=T)

wx<-lag(weights,GambiaPrev$Prev)
lhwx<-cut(wx, breaks=c(min(wx), mean(wx), max(wx)), labels=c("L", "H"), include.lowest=T)
lhlh<-interaction(lhx,lhwx,infl,drop=T)
cols<-rep(1, length(lhlh))
cols[lhlh=="L.L.TRUE"]<-2
cols[lhlh=="H.L.TRUE"]<-3
cols[lhlh=="L.H.TRUE"]<-4
cols[lhlh=="H.H.TRUE"]<-5

sum(cols>1) ##7 outliers

# Map points
plot(GMB_Adm2,border="darkgrey")
points(xy, col=c("lightgrey", "cyan", "coral","coral4","cyan4")[cols], pch=16, cex=0.6, add=T)
legend("topright", legend=c("None", "LL","HL", "LH", "HH"),
       fill=c("lightgrey", "cyan", "coral","coral4","cyan4"), bty="n", cex=0.8,
       y.intersp=0.8)


# Q5: In the Moran scatter plot, the centred weighted average is plotted against the observations. 
# Depending on their position on the plot, the Moran plot data points express the level of spatial 
# association of each observation with its neighbouring ones. Are the results plotted on the map as 
# expected based on the scatter plot?