spatialAutocorrelation.R

library(gstat)
library(latticeExtra)
library(raster)
library(sp)
library(spdep)
library(tidyverse)

set.seed(345)
size <- 100

xy <- expand.grid(1:size, 1:size)%>%
    setNames(c("x", "y"))
g <- gstat::gstat(formula = z~1, 
                  locations = ~x+y, 
                  dummy = T, 
                  beta = 0, 
                  model = vgm(psill = 0.1, 
                              range = 10, 
                              model = "Exp"), 
                  nmax = 20)
yy <- predict(g, newdata = xy, nsim = 2)

# Normalize Data from 0 to 1 
# This is strictly unnecessary but given the example is an interval that makes sense
yy <- yy %>% 
  mutate(sim1 = (sim1 - min(sim1))/(max(sim1)- min(sim1)),
         sim2 =  (sim2 - min(sim2))/(max(sim2) - min(sim2)))

sp::gridded(yy) <- ~x + y # will return True to indicate it's a grid 

# Sample Points
# set.seed(123)
samples <- 50 
x <- runif(samples, 0, size) %>%
  floor + 1 # get us to nearest integer 
y <- runif(samples, 0, size)%>%
  floor + 1
dataPoints <- data.frame(x,y)
y1 <- raster::raster(yy["sim1"])
y2 <- raster::raster(yy["sim2"])

plot_data <- data.frame(sim1 = raster::extract(y1, dataPoints), 
                        sim2 = raster::extract(y2, dataPoints))

# plot grids 
spplot(obj = yy, 
           names.attr = c("Hopscotch Enthusiasts in 1900", 
                          "Berkeley Admit Rate"),
           main = list("An obviously spurious correlation between two Spatial Variables"))+
  layer(panel.points(x,y, col = "white", pch = 19), data = dataPoints)

# Estimate the spatial dependence 
out <- summary(lm(sim1 ~ sim2, data = plot_data))

# Calculate Moran's Index via Monte Carlo 
plot_data_sp <- st_as_sf(plot_data %>% cbind(dataPoints), coords = c("x", "y"), crs=NA, agr= "constant")

# Get sum of the weights out to a distance of 20 away 
# This makes sense given our construction. Consult a variogram for a different dataset
Wij <- spdep::dnearneigh(plot_data_sp, 0, 20)
dlist <- nbdists(Wij, st_coordinates(plot_data_sp))
ilist <- lapply(dlist, function(x) exp(-x))
                
# Get weight matrix 
lw <- nb2listw(Wij, glist = ilist, style = "W", zero.policy = TRUE)

# Moran test MC 
MC <- moran.mc(plot_data$sim1, lw, nsim = 10000)
plot(MC, main = "", las=1)