###########################################################################################################
#  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//JLS - Point processes//"


# Attach libraries for visualisation
library(rgdal)

library(geoR)
library(gcmr)


# Attach libraries for point processes
library(sp)
library(spatstat)
library(splancs)
library(smacpod) # Spatial scanning statistic

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

###########################################################################################################
#  Point processes

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

# 1. Simulating complete spatial randomness (CSR)
# Waller, L.A. and Gotway, C.A. (2004) Applied Spatial Statistics
#   for Public Health Data.  John Wiley and Sons: New York.

par(mfrow=c(2,3),pty="s")
for (i in 1:6) {
  x <- runif(30) #generates 30 random deviates from the uniform distribution
  y <- runif(30)
  plot(x,y,xlim=c(0,1),ylim=c(0,1),xlab="u",ylab="v",cex.lab=1.5,cex.axis=1.1)
}

###########################################################################################################
#  Part II:  Global cluster detection for event data 

###########################################################################################################
# Data input and visualisation
###########################################################################################################

# 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//"

setwd(paste(root, path, "Lab 2", sep=""))

# Read in case-control data
CaseControl<-read.csv("CaseControl_Date.csv", header=T, stringsAsFactors=FALSE)

project2<-"+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 " #
StudyArea<-readShapePoly("StudySiteDistricts.shp",IDvar=NULL, proj4string=CRS(project2))

###           ADD CODE    #2                     ##
# Plot the data using techniques you have already learned
# Q1: Visually assess the presence of clustering across the study area.

# Promote data frame to spatial
coordinates(CaseControl) = ~long + lat
class(CaseControl)  # class is now SpatialPointsDataFrame
CaseControl@data # View attribute data

# Take cases only 
cases<-CaseControl[CaseControl@data$case==1,]
controls<-CaseControl[CaseControl@data$case==0,]
coordinates(cases)

# back to data
#CaseControl<-as.data.frame(CaseControl)
###########################################################################################################
#  Ripley's K function

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

# Ripley's K function
K<-Kest(CasesPPP)                                                                                   # "spatstat" package
par(mfrow=c(1,1))
plot(K)
E<-envelope(CasesPPP, Kest, nsim=999)
plot(E)

# Q2: 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

# 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])
plot(KX, sqrt(iso/pi) ~ r)

# Calculate the K-function for controls
KY <- Kest(CaseControlPPP[CaseControlPPP$marks==0])
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)                                        #"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

# Approach #3 - Bivand et al 2008. Goes through the code and process in which above smacpod functions are based on
# Create a grid over the study area where the risk ratio will be estimated
sG<-Sobj_SpatialGrid(StudyArea)$SG
gt<-slot(sG,"grid")

# compute the risk ratio by estimating the intensity of cases and controls and taking the ratio
grid<-slot(slot(slot(StudyArea, "polygons")[[1]], "Polygons")[[1]], "coords")

mean(dist(coordinates(cases)))

s<-seq(0,0.12, by=0.0001) #A vector of distances at which to calculate the K function

# Calculate the K function for cases and controls over the grid
khcases<-khat(coordinates(cases), grid, s)
khcontrols<-khat(coordinates(controls), grid, s)

# Calculate the covariance matrix for the difference between two K-functions under random labelling of the corresponding two sets of points.
khcov<-khvmat(coordinates(cases), coordinates(controls), grid, s)

# Test statistic
T0<-sum(((khcases-khcontrols))/sqrt(diag(khcov+0.00000000001)))

niter<-999  # Number of iterations to simulate
T<-rep(NA,niter) # Empty vector in which to store simulated values of T

khcasesrel<-matrix(NA, nrow=length(s), ncol=niter)
khcontrolsrel<-matrix(NA, nrow=length(s), ncol=niter)

# Simulate values
DIFF<-NULL
for(i in 1:niter){
  idxrel<-sample(CaseControl$case)==0
  casesrel<-coordinates(CaseControl[idxrel,])
  controlsrel<-coordinates(CaseControl[!idxrel,])
  khcasesrel[,i]<-khat(casesrel, grid, s)
  khcontrolsrel[,i]<-khat(controlsrel, grid, s)
  khdiff<-khcasesrel[,i]-khcontrolsrel[,i]
  DIFF<-cbind(DIFF,khdiff)
  length(khdiff); length(sqrt(diag(khcov)))
  T[i]<-sum(khdiff/sqrt(diag(khcov+0.00000000001)))  # check this line, did not work as entered
}

# Calculate the P-value (proportion of times the simulated value is more extreme than observed)
pvalue<-(sum(T>T0)+1)/(niter+1)
pvalue

# Estimate 95% CI by taking the maximum and minimum simulated value for each point estimated
UppEnv<-apply(DIFF,1, quantile, probs = 0.95)
LowerEnv<-apply(DIFF,1,quantile, probs = 0.05)

# Plot the data
plot(s, khcases-khcontrols, type="l", ylim=c(-0.01,0.01))
polygon(c(s,rev(s)),c(UppEnv, rev(LowerEnv)),col=gray(0.8))
lines(s, khcases-khcontrols, type="l", col="red")


###########################################################################################################
#//  Part IV:  Local cluster detection for point-process data 
# We will mainly use satscan for this part of the practical, but a quick intro to packages available in R
#for you to look into on your own
###########################################################################################################

###           ADD CODE    #3                     ##
# Export data as instructed in part IV of Lab 2.


# Approach 1: Use SatScan to locate clusters (use lab)


# 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))


# Approach 2: Using "smacpod" library for spatial scan test of Kulldorf (1997)

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

draw.circle(out$clusters[[1]]$coords[,1], out$clusters[[1]]$coords[,2], out$clusters[[1]]$r, border="red")