# Cheat sheet
# 1. 

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

# Visualize
pal = colorNumeric("Oranges", GambiaPrev$Prev)
leaflet(GambiaPrev) %>% addTiles() %>% addCircleMarkers(~x, ~y, fillOpacity=1,
                                                             fillColor= ~pal(Prev), radius=~Prev*10, stroke=TRUE, weight=1) %>% 
  addLegend(pal = pal, values = ~Prev)

# Calc log_odds
GambiaPrev$log_odds <- logit(GambiaPrev$Prev+0.0001)

# Generate a distance matrix
GM.dists <- as.matrix(dist(cbind(GambiaPrev$x, GambiaPrev$y)))
dim(GM.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
GM.dists.inv <- 1/GM.dists
diag(GM.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$log_odds, GM.dists.inv) 


# 2 

# Visualize observed data : can use this to calculate moran's i but does not account for 
# population distribution (which we do not have)
pal <- colorBin("YlGn", domain = scot_LL$Observed, bins = 5)
leaflet(scot_LL) %>% addTiles() %>% 
  addPolygons(fillColor = ~pal(Observed), 
              fillOpacity = 0.7, 
              weight = 1) %>%
  addLegend(pal = pal, 
            values = ~Observed,
            labFormat = labelFormat(prefix = ""))


#Standardize residuals (o-E)/sqrt(E) : We do have expected values based on age standardization
scot_LL@data$StdCt<-(scot_LL@data$Observed-scot_LL$Expected)/sqrt(scot_LL@data$Expected)

# Visualize standardized residuals 
pal <- colorBin("YlGn", domain = scot_LL$StdCt, bins = 5)
leaflet(scot_LL) %>% addTiles() %>% 
addPolygons(fillColor = ~pal(StdCt), 
            fillOpacity = 0.7, 
            weight = 1) %>%
  addLegend(pal = pal, 
            values = ~StdCt,
            labFormat = labelFormat(prefix = ""))

# estimate centers of polygons
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,1))
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

# 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(scot_LL$StdCt, scot_W)                                                     # "spdep" package


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

# Plot the spatial data against its spatially lagged values (weighted mean of neighbors)                         
nci<-moran.plot(scot_LL$StdCt, listw=scot_W, 
                xlab="StdCt", ylab="Spatially lagged", labels=T, pch=16, col="grey")


# Map points that are local outliers in the plot
infl<-apply(nci$is.inf,1,any) # which are statistically significant
sum(infl==T)#6 true (11% - more than would expect by chance)
x<-scot_LL$Observed
#Cut data into"Low" and "High"
lhx<-cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=T)
# Calculate spatial lag
wx<-lag(scot_W,scot_LL$Observed)
#Cut lag into "Low" and "High"
lhwx<-cut(wx, breaks=c(min(wx), mean(wx), max(wx)), labels=c("L", "H"), include.lowest=T)
#Interact (LH) and (LH) to form (LL, LH, HL, HH)
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") ##6 outliers
scot_LL$names<-names


# Map localities
factpal <- colorFactor(c( "cyan4","coral4","coral","cyan","lightgrey"), names)
leaflet(scot_LL) %>% addTiles() %>% 
  addPolygons(fillColor = ~factpal(names), 
              fillOpacity = 0.7, 
              weight = 1) %>%
  addLegend(pal = factpal, 
            values = ~names, title="Class")


# 3 
NairobiPPP<-ppp(nairobi_cases$lng, nairobi_cases$lat, range(nairobi_cases$lng), range(nairobi_cases$lat), marks = as.factor(nairobi_cases$case))

kdest = kdest(NairobiPPP, case = 2,nsim=999, level=0.95, correction=c("isotropic", "Ripley"))                                        #"smacpod" package
kdest
plot(kdest)

kdplus.test(kdest)