# Attach libraries for visualisation
library(rgdal)
library(raster)
library(ggplot2)
library(spatstat)
library(plotrix)
library(fields)

library(leaflet)
library(plotGoogleMaps)
library(maptools)
library(RColorBrewer)
library(lattice)
library(geoR)
library(plotrix) #  Plotting spatial scanning 

# Data management / point processes
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 libraries for point processes
library(spatstat)
library(splancs) # K-function
library(smacpod) # Spatial scanning statistic
##library(stpp) # Spatio-temporal K function


# Global measures of clustering for point (prevalence etc.) data
# Open BF malaria data
BF_malaria_data <- read.csv("https://www.dropbox.com/s/bfs3pinxe1lvvxr/data_bf2_binomial.csv?dl=1",
                            header=T)

BF_Adm_1 <- raster::getData("GADM", country="BFA", level=1)
proj4string(BF_Adm_1) <- CRS('+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ')

# Calc prevalence
BF_malaria_data$prevalence <- BF_malaria_data$positives / BF_malaria_data$examined

# Remind yourself of what data look like - do you see evidence of spatial clustering?
pal = colorNumeric("Oranges", BF_malaria_data$prevalence)
leaflet(BF_malaria_data) %>% addTiles() %>% addCircleMarkers(~longitude, ~latitude, fillOpacity=1,
                                              fillColor= ~pal(prevalence), radius=~prevalence*10, stroke=TRUE, weight=1) %>% 
                                                addLegend(pal = pal, values = ~prevalence)

###########################################################################################################
# Moran's i and correlograms: often used as a formal statistical test if there is presence of spatial
# autocorrelation across the study area
# NB can do this with multiple packages
###########################################################################################################

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

hist(BF_malaria_data$prevalence)
# Use the logistic transformation (log odds) to make more normal
library(car)
BF_malaria_data$log_odds <- logit(BF_malaria_data$prevalence)

# Generate a distance matrix
BF.dists <- as.matrix(dist(cbind(BF_malaria_data$longitude, BF_malaria_data$latitude)))
dim(BF.dists) # 109 x 109 matrix of distance between all sets of points

# Take the inverse of the matrix values so that closer values have a larger weight and vs vs
BF.dists.inv <- 1/BF.dists
diag(BF.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(BF_malaria_data$log_odds, BF.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(BF_malaria_data$longitude, BF_malaria_data$latitude)))
maxDist

xy=cbind(BF_malaria_data$longitude, BF_malaria_data$latitude)
pgi.cor <- correlog(coords=xy, z=BF_malaria_data$log_odds, method="Moran", nbclass=10)               # "pgirmess" package
# coords = xy cordinates, z= vector of values at each location and nbclass = the number of bins
dev.off()
plot(pgi.cor) # statistically significant values (p<0.05) are plotted in red
pgi.cor # distclass is midpoint for the bin

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

# Compare the correlogram to the results from the semivariogram in the last class (copied below)
BF_malaria_data_geo<-as.geodata(BF_malaria_data[,c("longitude","latitude","log_odds")])

# Generate and plot a binned variogram (10 bins) NB: have made for full max distance (even though likely inaccurate) for comparison
Vario<-variog(BF_malaria_data_geo,max.dist=7.53,uvec=seq(0.4121237,7.1595572,l=10))
par(mfrow=c(2,1))
plot(Vario)
plot(pgi.cor)


# Calculate moran's I using a binary distance matrix

#Q3: What do you think 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) # set spatial coordinates to create a spatial object
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 one 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) # list of neighbor structures
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(2,1), mar= c(1, 0, 1, 0))
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="W")   # row standardized binary weights, using minimum distance for one neighbor
weights                                   # "B" is simplest binary weights

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

# Simulate the test statistic using random permutations of BF_malaria_data$log_oddsso 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(BF_malaria_data$log_odds , listw=weights,nsim=999)
bperm
#statistic = 0.15, observed rank = 1000, p-value = 0.001

# Plot simulated test statistics
par(mfrow=c(1,1), mar= c(5, 4, 4, 2))
hist(bperm$res, freq=T, breaks=20, xlab="Simulated Moran's I")
abline(v=0.15, 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(BF_malaria_data$log_odds, 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(BF_malaria_data$log_odds, listw=weights, 
                xlab="Prevalence", ylab="Spatially lagged prev", labels=T, pch=16, col="grey")
text(c(3,3, -5,-5),c(0.9, -1.9,0.9,-1.9), c("High-High", "High-Low", "Low-High", "Low-Low"), cex=0.8)

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)#13 true (12% - more than would expect by chance)
x<-BF_malaria_data$log_odds
lhx<-cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=T)

wx<-lag(weights,BF_malaria_data$log_odds)
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)

names<-rep("none", length(lhlh))
names[lhlh=="L.L.TRUE"]<-"LL"
names[lhlh=="H.L.TRUE"]<-"HL"
names[lhlh=="L.H.TRUE"]<-"LH"
names[lhlh=="H.H.TRUE"]<-"HH"


sum(names!="none") ##7 outliers
BF_malaria_localM<-as.data.frame(cbind(xy,names))
colnames(BF_malaria_localM)<-c("longitude", "latitude", "names")
BF_malaria_localM[c("longitude", "latitude")] <- lapply( BF_malaria_localM[c("longitude", "latitude")], function(x) as.numeric(as.character(x)) )

# Map points
factpal <- colorFactor(c( "cyan4","coral4","coral","cyan","lightgrey"), names)
leaflet(BF_malaria_localM) %>% addTiles() %>% addCircleMarkers(~longitude, ~latitude, fillOpacity=1,
                                                               color= ~factpal(names), radius=4, stroke=TRUE, weight=1) %>% 
  addLegend(pal = factpal, values = ~names, title="Class")


# In class assignment on Moran's i and local moran's i

# 1. Visualize data and test whether there is global spatial autocorrelation apparent using malaria point data from the Gambia
# NB you will need to use some techniques from earlier lectures to calculate prevalence

GM_malaria_data <- read.csv("https://www.dropbox.com/s/dgaldemru009pir/GambiaID.csv?dl=1")


# 2. Map the incidence of lip cancer using the Scottish lip cancer dataset and visually assess for local high or low clusters
#    Calculate LISA's to test for any statistically significant local clusters
#library(gcmr)   # data available in package 'gcmr'
#data("scotland")
library(classInt)

scot <- read.csv("https://www.dropbox.com/s/ngizuxblkazplk8/scotland.csv?dl=1")
# Shapefile available for download online
scot_LL <- readOGR("/Users/smithj1/Box Sync Copy for other comp2/Spatial Epi Course/2016/Lecture 5 Analysis of spatial clustering/Lab 1", "scot")
proj4string(scot_LL) <- CRS("+proj=longlat +datum=WGS84") 
# Link data to shapefile
row.names(scot) <- formatC(scot$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[match(ID, row.names(scot)), ])
plot(scot_LL)

###########################################################################################################
#  Part II:  Spatial point processes

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

# Open Namibia malaria case data
CaseControl<-read.csv("https://www.dropbox.com/s/hwma9q39z2axjcg/CaseControl.csv?dl=1")


# And boundary file
NAM_Adm0<-raster::getData('GADM',country='NAM',level=0)

# Convert to a SPDF
CaseControl_SPDF <- SpatialPointsDataFrame(coords = CaseControl[,c("long", "lat")],
                                           data = CaseControl[,c("household_id", "case")])
cases<-CaseControl_SPDF[CaseControl$case==1,]
controls<-CaseControl_SPDF[CaseControl$case==0,]

# Let's plot and see what we have
case_color_scheme <- colorNumeric(c("blue", "red"), CaseControl_SPDF$case)
leaflet() %>% addTiles() %>% addCircleMarkers(data=CaseControl_SPDF, color = case_color_scheme(CaseControl_SPDF$case))


# As part of the last lecture, you already generated first order kernel density estimates and 
# calculatred the ratio of the density estimate of cases:controls
# Now you will look at second order functions, summarizing the spatial dependence between events

# Promote the case data to a PPP data type
CasesPPP<-as(cases, "ppp")

# Ripley's K function: summarizes the spatial dependence between events at wide range of spatial scales
K<-Kest(CasesPPP,correction=c("isotropic", "Ripley"))                                                                                   # "spatstat" package

par(mfrow=c(1,1)) # Plot the estimate of K(r); note different border-corrected estimates ('iso', 'border' and 'trans')
plot(K, xlab="d (dd)", ylab="K(dd)") # Red dashed line is expected K value computed for a CRS process
E<-envelope(CasesPPP, Kest, nsim=999) # Plot confidence envelope using MC simulation
plot(E)


# Q: The K-function computed for cases assumes that H0 is complete spatial randomness. What are the limitations of this assumption?


###########################################################################################################
# Difference in Ripley's K function between cases and controls
# Two approaches below that do essentially the same thing; #2 with hypothesis testing
###########################################################################################################

# Approach #1 
# K function vignette from Bradley et simply calculates the K function for cases and controls, and evaluates the difference

# Create a marked point process
CaseControlPPP<-ppp(CaseControl$long, CaseControl$lat, range(CaseControl$long), range(CaseControl$lat), marks = as.factor(CaseControl$case))

# Calculate the K-function for cases
KX <- Kest(CaseControlPPP[CaseControlPPP$marks==1],correction=c("isotropic", "Ripley"))
plot(KX, sqrt(iso/pi) ~ r)

# Calculate the K-function for controls
KY <- Kest(CaseControlPPP[CaseControlPPP$marks==0],correction=c("isotropic", "Ripley"))
plot(KY, sqrt(iso/pi) ~ r)

# Calulate the difference in the two functions
Kdiff <- eval.fv(KX - KY)
plot(Kdiff, legendpos="float")


# Approach #2 "Smacpod" package includes a function to estimate the difference in K function and plot simulated CI
# Also includes a function to the test the sifnificance based on these simulations

kdest = kdest(CaseControlPPP, case = 2,nsim=999, level=0.95, correction=c("isotropic", "Ripley"))                                        #"smacpod" package
# Note that the case = is position of the marks, not the value!  levels(CaseControlPPP$marks)
kdest
plot(kdest)
# grey is min/max; blue is confidence envelope (can change these with options)

kdplus.test(kdest)
# Performs test of significance based on simulated confidence envelope and observed statistic

###########################################################################################################
# In class assignment

# 3. Estimate the scale and type of clustering using the K functions in the simulated case-control data
nairobi_cases <- read.csv("https://www.dropbox.com/s/ds839x22jmyyihy/cases_nairobi.csv?dl=1")


# 4. Open the SatScan_Practical word document and follow the instructions to identify clustering of cases using 
#   Kulldorff's spatial scan statistic. This corresponds to the Namibia data.

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

# After running the SatScan analysis, use "smacpod" library for spatial scan test of Kulldorf (1997)
# Convert CaseControl to a "PPP" object for spatial scan
CaseControlPPP<-ppp(CaseControl$long, CaseControl$lat, range(CaseControl$long), range(CaseControl$lat), marks = as.factor(CaseControl$case))

out<-spscan.test(CaseControlPPP, nsim = 999, case = 2, maxd=0.15, alpha = 0.05)           # "smacpod" library
out
plot(CaseControlPPP)

case_color_scheme <- colorNumeric(c("blue", "red"), CaseControl_SPDF$case)
leaflet() %>% addTiles() %>% addCircleMarkers(data=CaseControl_SPDF, color = case_color_scheme(CaseControl_SPDF$case), 
                                              stroke = FALSE, radius=2, fillOpacity=1)%>% 
  addCircles(lng = out$clusters[[1]]$coords[,1], lat = out$clusters[[1]]$coords[,2], weight = 2,
             radius = out$clusters[[1]]$r*112*1000, color="grey")
  
# Read in shapefiles from SatScan and visualize for comparison
SatScan <- readOGR("/Users/smithj1/Box Sync Copy for other comp2/Spatial Epi Course/2017/Cluster Analysis", "SatScan_output.col")
SatScan@data


leaflet() %>% addTiles() %>% addCircleMarkers(data=CaseControl_SPDF, color = case_color_scheme(CaseControl_SPDF$case), 
                                              stroke = FALSE, radius=2, fillOpacity=1)%>% 
  addCircles(lng = out$clusters[[1]]$coords[,1], lat = out$clusters[[1]]$coords[,2], weight = 2,
             radius = out$clusters[[1]]$r*112*1000, color="grey")%>% 
  addCircles(lng = SatScan$LONGITUDE[1], lat = SatScan$LATITUDE[1], weight = 2,
             radius = SatScan$RADIUS[1]*1000, color="green")
  

## Assignment
# Use SatScan to carry out an analysis of the count data from the scottish lip cancer
# dataset. What model would you use for these data? Export and map the results.