#####################
### Random fields ###
#####################

###############
# Ising model #
###############

# Wolfram Demonstration:
# http://msekce.karlin.mff.cuni.cz/~dvorak/teaching/2019_2020/teaching_2019_2020.html


# The following code works fine with R version 3.1.2
# (but not the newest one). See comments below.

# We will need a couple of external libraries,
# please either run the following command or install the libraries manually
install.packages(c("spdep","maptools","RANN"))


###################
# library 'spdep' #
###################

library(spdep)

# examples of regional data
data(package="spdep")


library(maptools)
# crime data from 49 regions ('neighbourhoods') in Columbus (Ohio)
# https://en.wikipedia.org/wiki/Columbus,_Ohio
example(columbus)
plot(columbus,border="grey")
coords <- coordinates(columbus)
points(coords,pch=19)

# If a nice graph with vertices and edges is plotted,
# with some background lines depicting the neighbourhoods,
# the code works fine with your R version.

# If instead you see a strange kind of matrix with
# very small plots in each position, the code does not work well.


### neigbours
col.gal.nb                         # neigbourhood relation
plot(col.gal.nb,coords,add=TRUE)
summary(col.gal.nb)
is.symmetric.nb(col.gal.nb)
print(cards <- card(col.gal.nb))   # counts of neigbours


### weights for assessing the spatial autocorrelation
col.b <- nb2listw(col.gal.nb,style="B")   # binary weights
col.b$weights
col.w <- nb2listw(col.gal.nb,style="W")   # normalized binary weights (sum up to 1 for each point)
col.w$weights


### second-order neigbours
col.lags <- nblag(col.gal.nb, 2)   # neighbours of neighbours
plot(col.lags[[2]], coords, add=TRUE, col="red", lty=2)
summary(col.lags[[2]], coords)


### k nearest neighbours (in Euclidean space, not necessarily in the neigbourhood relation)
library(RANN)
col.k4.nb <- knn2nb(knearneigh(coords, 4))  # four nearest neighbours
is.symmetric.nb(col.k4.nb)
plot(columbus,border="grey")
plot(col.k4.nb,coords,arrows=TRUE,add=TRUE)


### neighbourhood relation based on the Euclidean distance
col.d.nb <- dnearneigh(coords,0,0.5)     # neighbours if closer than 0.5
is.symmetric.nb(col.d.nb)
plot(columbus,border="grey")
plot(col.d.nb,coords,add=TRUE)



### continuous data
# residential burglaries and vehicle thefts per thousand households
columbus$CRIME
summary(columbus$CRIME)


# Tests for "no spatial autocorrelation" hypothesis:

# We assume a constant mean value and constant variance!
# Moran index for normalized weights, test based on the randomization assumption
# (all permutations of assigning the observed values to the lattice points are equally probable)
moran.test(columbus$CRIME,listw=col.w)

# Moran index for normalized weights, test based on the normality assumption
# (uncorrelated random variables with Gaussian distribution)
moran.test(columbus$CRIME,listw=col.w,randomisation=FALSE)

# Geary index for normalized weights, test based on the randomization assumption:
geary.test(columbus$CRIME,listw=col.w)



### binary data
HICRIME <- cut(columbus$CRIME, breaks=c(0,35,80), labels=c("low","high"))
plot(columbus)
points(coords[HICRIME=="high",],col="red",pch=19)
points(coords[HICRIME=="low",],col="blue",pch=19)

# BB statistic = 1/2 sum_{i,j} ( w_ij * Z_i * Z_j )
# (here we assume that the two possible values are 0 and 1)
joincount.test(HICRIME, listw=col.w)
joincount.test(HICRIME, listw=col.b)
# the tests are performed under the assumption of 'nonfree sampling', i.e. the hypergeometric sampling;
# we have a fixed number of red and blue regions