Part 1: Basic Autocorrelation

1.1 Autocorrelation Data Generation

set.seed(0) 
d <- sample(100, 10) 
d 
##  [1]  14  68  39   1  34  87  43 100  82  59

set.seed(0): Ensures reproducibility by setting the random number generator to a fixed state. sample(100, 10): Generates a random sample of 10 numbers from 1 to 100.

1.2 Remove First and Last Values

a <- d[-length(d)] 
a
## [1]  14  68  39   1  34  87  43 100  82
b <- d[-1] 
b
## [1]  68  39   1  34  87  43 100  82  59

1.3 Plotting Autocorrelation Pairs

plot(a, b, xlab='t', ylab='t-1')

1.4 Calculating Correlation

cor(a, b) 
## [1] 0.1227634

1.5 Autocorrelation Function (ACF)

acf(d)

# Part 2: Spatial Data Preparation ## 2.1 Loading the Terra Library

 library(terra) 
## terra 1.7.83

2.2 Importing Spatial Data

p <- vect(system.file("ex/lux.shp", package="terra"))
head(p)

2.3 Filtering Spatial Data

p <- p[p$NAME_1=="Diekirch", ] 
head(p)

2.4 Adding a New Column

p$value <- c(10, 6, 4, 11, 6)
head(p)

Part 3: Visualization & Neighbor Analysis

3.1 Plotting Spatial Data

par(mai=c(0,0,0,0)) 
plot(p, col=2:7) 
xy <- centroids(p) 
points(xy, cex=6, pch=20, col='white') 
text(p, 'ID_2', cex=1.5)

## 3.2 Identifying Adjacent Polygons

library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
w <- adjacent(p, symmetrical=TRUE) 
class(w) 
## [1] "matrix" "array"
head(w) 
##      from to
## [1,]    1  2
## [2,]    1  4
## [3,]    1  5
## [4,]    2  3
## [5,]    2  4
## [6,]    2  5

3.3 Plotting Adjacency Relationships

plot(p, col='gray', border='blue', lwd=2) 
p1 <- xy[w[,1], ] 
p1
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 7, 7  (geometries, attributes)
##  extent      : 5.886502, 6.127425, 49.80014, 50.07064  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP value
##  type        : <num>    <chr> <num>    <chr> <num> <int> <num>
##  values      :     1 Diekirch     1 Clervaux   312 18081    10
##                    1 Diekirch     1 Clervaux   312 18081    10
##                    1 Diekirch     1 Clervaux   312 18081    10
p2 <- xy[w[,2], ]
p2
##  class       : SpatVector 
##  geometry    : points 
##  dimensions  : 7, 7  (geometries, attributes)
##  extent      : 5.886502, 6.165081, 49.80014, 49.93892  (xmin, xmax, ymin, ymax)
##  coord. ref. : lon/lat WGS 84 (EPSG:4326) 
##  names       :  ID_1   NAME_1  ID_2   NAME_2  AREA   POP value
##  type        : <num>    <chr> <num>    <chr> <num> <int> <num>
##  values      :     1 Diekirch     2 Diekirch   218 32543     6
##                    1 Diekirch     4  Vianden    76  5163    11
##                    1 Diekirch     5    Wiltz   263 16735     6
lines(p1, p2, col='red', lwd=2)

## 3.4 Creating Adjacency Matrix

wm <- adjacent(p, pairs=FALSE) 
wm 
##   1 2 3 4 5
## 1 0 1 0 1 1
## 2 1 0 1 1 1
## 3 0 1 0 0 1
## 4 1 1 0 0 0
## 5 1 1 1 0 0

Part 4: Moran’s I and Spatial Correlation

4.1 Computing Moran’s I

n <- length(p) 
y <- p$value 
ybar <- mean(y) 
dy <- y - ybar 
yi <- rep(dy, each=n) 
yj <- rep(dy) 
yiyj <- yi * yj 
pm <- matrix(yiyj, ncol=n) 
pmw <- pm * wm 
pmw 
##       1     2    3     4     5
## 1  0.00 -3.64 0.00  9.36 -3.64
## 2 -3.64  0.00 4.76 -5.04  1.96
## 3  0.00  4.76 0.00  0.00  4.76
## 4  9.36 -5.04 0.00  0.00  0.00
## 5 -3.64  1.96 4.76  0.00  0.00
spmw <- sum(pmw) 
spmw 
## [1] 17.04
smw <- sum(wm) 
sw <- spmw / smw 
vr <- n / sum(dy^2) 
MI <- vr * sw 
MI 
## [1] 0.1728896
EI <- -1/(n-1) 
EI 
## [1] -0.25
ww <- adjacent(p, "queen", pairs=FALSE) 
ww 
##   1 2 3 4 5
## 1 0 1 0 1 1
## 2 1 0 1 1 1
## 3 0 1 0 0 1
## 4 1 1 0 0 0
## 5 1 1 1 0 0
ac <- autocor(p$value, ww, "moran") 
ac
## [1] 0.1728896

n <- length(p): Determines the number of spatial units. y <- p\(value: Extracts the value column for analysis. ybar <- mean(y): Calculates the mean of y. dy <- y - ybar: Centers the data by subtracting the mean. yi & yj: Create replicated vectors for pairwise multiplication. yiyj <- yi * yj: Computes the product of deviations for each pair. pm <- matrix(yiyj, ncol=n): Reshapes the products into a matrix. pmw <- pm * wm: Applies the adjacency weights to the products. spmw <- sum(pmw): Sums the weighted products. smw <- sum(wm): Sums the adjacency weights. sw <- spmw / smw: Calculates the spatially weighted sum. vr <- n / sum(dy^2): Computes the variance ratio. MI <- vr * sw: Calculates Moran’s I. EI <- -1/(n-1): Computes the Expected Index under spatial randomness. ww <- adjacent(p, "queen", pairs=FALSE): Generates adjacency weights using the Queen criterion. ac <- autocor(p\)value, ww, “moran”): Final Moran’s I calculation using autocor function.

4.2 Monte Carlo Simulation for Significance Testing

# test for significance using Monte Carlo simulatio
m <- sapply(1:99, function(i) { autocor(sample(p$value), ww, "moran") }) 
hist(m)

4.3 Significance level of the spatial autocorrelation

pval <- sum(m >= ac) / 100 
pval
## [1] 0.04

Part 5: Moran Scatter Plot

5.1 Preparing Neighboring Values

n <- length(p) 
ms <- cbind(id=rep(1:n, each=n), y=rep(y, each=n), value=as.vector(wm * y))

5.2 Removing Zero Values

ms <- ms[ms[,3] > 0, ]

5.3 Aggregating Neighboring Values

ams <- aggregate(ms[,2:3], list(ms[,1]), FUN=mean) 
ams <- ams[,-1] 
colnames(ams) <- c('y', 'spatially lagged y') 
head(ams)

5.4 Plotting the Moran Scatter Plot

plot(ams, pch=20, col="red", cex=2) 
reg <- lm(ams[,2] ~ ams[,1]) 
abline(reg, lwd=2) 
abline(h=mean(ams[,2]), lt=2) 
abline(v=ybar, lt=2)