Autocorrelación espacial

Análisis espacial con R

Dr. Francisco Zambrano

Introducción

La autocorrelación espacial es un concepto clave en estadística espacial.
Es tanto una molestía, por que complica los test estadísticos, y una ventaja ya que permite la interpolación espacial.
La autocorrelación (sea espacial o no) es una medida de similaridad (correlación) entre observaciones cercanas.
Para entender autocorrelación espacial, ayuda primero considerar la autocorrelación temporal.

Temporal autocorrelation

Si mides algo sobre un objeto en el tiempo, por ejemplo el peso de una persona, es probable que dos observaciones que están cerca tiene valores similares.
Digamos que tu peso en algunos años varió de 50kg a 80kg.
Es poco probable que fue de 50kg un día y al día siguiente 80kg.
El proceso debe haber sido gradual.
Para medir el grado de asociación en el tiempo, podemos calcular la correlación de cada observación con respecto a la siguiente.

Autocorrelación temporal

Un ejemplo

Sea d un vector de observaciones diarias

set.seed(0)
(d <- sample(100, 10))

 [1]  14  68  39   1  34  87  43 100  82  59

Autocorrelación temporal

Un ejemplo

Calculamos la autocorrelación con desfase de un día

a <- d[-length(d)]
b <- d[-1]
plot(a, b, xlab='t', ylab='t-1')

Autocorrelación temporal

Un ejemplo

La autocorrelación es muy baja.

cor(a,b)

[1] 0.1227634

Pero es una muestra aleatoria
Calculamos la autocorrelación con un desfase de 1, para los valores inmediatamente próximos.

Autocorrelación temporal

Un ejemplo

Si ordenamos los datos mejoramos la autocorrelación.

(d <- sort(d))

 [1]   1  14  34  39  43  59  68  82  87 100

a <- d[-length(d)]
b <- d[-1]

Autocorrelación temporal

Un ejemplo

Si ordenamos los datos mejoramos la autocorrelación.

cor(a,b)

[1] 0.9819258

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

Autocorrelación temporal

Un ejemplo

La función acf muestra la autocorrelación para diferentes desfases.

acf(d)

Autocorrelación espacial

El concepto de autocorrelación espacial es una extensión del de autocorrelación temporal.
Es un poco más complicado
- El tiempo tiene una dimensión
- Objetos espaciales tiene al menos dos dimensiones
Las mediciones de autocorrelación espacial describe el grado al cual dos observaciones en ubicaciones espaciales (sea puntos, áreas, raster) son similares.
Necesitamos observaciones y sus ubicaciones

Autocorrelación espacial

Puede ser exógena: causada por algun otra variable autocorrelacionada
Puede ser endogena: causada por el mismo proceso (ej, covid-19)
Un estadístico usado para describir la autocorrelación espacial es el índice de Moran.
Otros son Geary y semivariograma

Autocorrelación espacial

Un ejemplo

library(terra)
p <- vect(system.file("external/lux.shp", package="raster"))
p <- p[p$NAME_1=="Diekirch", ]
#p$value <- c(10, 6, 4, 11, 6)
data.frame(p)

  ID_1   NAME_1 ID_2   NAME_2 AREA
1    1 Diekirch    1 Clervaux  312
2    1 Diekirch    2 Diekirch  218
3    1 Diekirch    3  Redange  259
4    1 Diekirch    4  Vianden   76
5    1 Diekirch    5    Wiltz  263

Autocorrelación espacial

Un ejemplo con polígonos

Calcular la autocorrelación de la variable AREA

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)

Autocorrelación espacial

Polígonos adyacentes

Debemos determinar que polígonos están cerca.
Utilizaremos adyacencia con el paquete {spdep}.

library(spdep)
w <- adjacent(p, symmetrical = TRUE)
class(w)

[1] "matrix" "array"

## [1] "nb"

Autocorrelación espacial

Polígonos adyacentes

summary(w)

      from             to     
 Min.   :1.000   Min.   :2.0  
 1st Qu.:1.000   1st Qu.:3.5  
 Median :2.000   Median :4.0  
 Mean   :1.714   Mean   :4.0  
 3rd Qu.:2.000   3rd Qu.:5.0  
 Max.   :3.000   Max.   :5.0

Autocorrelación espacial

Polígonos adyacentes

Conección entre polígonos

plot(p, col='gray', border='blue', lwd=2)
p1 <- xy[w[,1],]
p2 <- xy[w[,2],]
lines(p1,p2,col = 'red', lwd = 2)

Autocorrelación espacial

Polígonos adyacentes

Ahora podemos transforma w a una matriz de pesos espaciales.
La matriz de pesos refleja la intensidad de la relación entre las observaciones.

(wm <- adjacent(p, pairs=FALSE))

  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

Autocorrelación espacial

Indice de Moran

\[I = \frac{n}{\sum_{i=1}^n (y_i - \bar{y})^2} \frac{\sum_{i=1}^n \sum_{j=1}^n w_{ij}(y_i - \bar{y})(y_j - \bar{y})}{\sum_{i=1}^n \sum_{j=1}^n w_{ij}}\]

Parecida a la ecuación de calculo de correlación lineal
La mayor diferencia es que agrega la matriz de pesos.

El número de observaciones

(n <- length(p))

[1] 5

Obtenemos \(y\) e \(\bar y\)

y <- p$AREA
ybar <- mean(y)

Autocorrelación espacial

Indice de Moran

Ahora necesitamos \((y_i - \bar{y})(y_j - \bar{y})\)

dy <- y - ybar
g <- expand.grid(dy, dy)
yiyj <- g[,1] * g[,2]

Creamos la matriz con pares multiplicados

pm <- matrix(yiyj, ncol=n)

Multiplicamos la matríz por la matriz de pesos

pmw <- pm * wm

Autocorrelación espacial

Indice de Moran

Ahora sumamos los valores para obtener

\[\sum_{i=1}^n \sum_{j=1}^n w_{ij}(y_i - \bar{y})(y_j - \bar{y})\]

spmw <- sum(pmw)

El siguiente paso es dividir este valor por la suma de pesos

smw <- sum(wm)
sw  <- spmw / smw

Autocorrelación espacial

Indice de Moran

Y el inverso de la varianza

(vr <- n / sum(dy^2))

[1] 0.0001542391

El valor del índice de Moran es

(MI <- vr*sw)

[1] -0.1873494

Autocorrelación espacial

Indice de Moran

El valor esperado en ausencia de autocorrelación espacial.

\[E(I) = \frac{-1}{n-1}\]

Valores \(< E(I)\) indican autocorrelación espacial negativa
Valores \(> E(I)\) indican autocorrelación espacial positiva

En nuestro ejemplo tenemos

(EI <- -1/(n-1))

[1] -0.25

Y el valor del indice de Moran fue de -0.19

Autocorrelación espacial

Indice de Moran

Pero, tenemos la función autocor en el paquete terra
Creamos una matriz de pesos conadjacent.

(ww <-  adjacent(p, "queen", pairs = FALSE))

  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

Autocorrelación espacial

Indice de Moran

Pero, tenemos la función autocor en el paquete terra

autocor(p$AREA, ww, "moran")

[1] -0.1873494

Autocorrelación espacial

Indice de Moran

Podemos correr un análisis de significancia utilizando montecarlo/bootstraping

m <- sapply(1:99, function(i) {
    autocor(sample(p$AREA), ww, "moran")
})
hist(m)

Autocorrelación espacial

Indice de Moran

Graficamos para ver la autocorrelación espacial

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

      id   y value
 [1,]  1 312     0
 [2,]  1 312   218
 [3,]  1 312     0
 [4,]  1 312    76
 [5,]  1 312   263
 [6,]  2 218   312
 [7,]  2 218     0
 [8,]  2 218   259
 [9,]  2 218    76
[10,]  2 218   263
[11,]  3 259     0
[12,]  3 259   218
[13,]  3 259     0
[14,]  3 259     0
[15,]  3 259   263
[16,]  4  76   312
[17,]  4  76   218
[18,]  4  76     0
[19,]  4  76     0
[20,]  4  76     0
[21,]  5 263   312
[22,]  5 263   218
[23,]  5 263   259
[24,]  5 263     0
[25,]  5 263     0

Autocorrelación espacial

Indice de Moran

Graficamos para ver la autocorrelación espacial

El promedio en cada vecindario

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

    y spatially lagged y
1 312           185.6667
2 218           227.5000
3 259           240.5000
4  76           265.0000
5 263           263.0000

Autocorrelación espacial

Indice de Moran

El promedio en cada vecindario

Finalmente, el gráfico.

plot(ams)
reg <- lm(ams[,2] ~ ams[,1])
abline(reg, lwd=2)
abline(h=mean(ams[,2]), lt=2)
abline(v=ybar, lt=2)

Autocorrelación espacial

Indice de Moran

El promedio en cada vecindario

Finalmente, el gráfico.

Autocorrelación espacial

Indice de Moran

La pendiente

  ams[, 1] 
-0.2302019

Autocorrelación espacial

Indice de Moran

Otro ejemplo para clarificar

Autocorrelación espacial

¿Cómo será el índice de Moran en cada caso?

Autocorrelación espacial

Si miramos en los puntos de algunos de los cluster

Autocorrelación espacial

Si miramos en los puntos de algunos de los cluster

Autocorrelación espacial

Indice de Moran

Si miramos en los datos meuse en la variable de concentración de copper

library(spdep)
library(sp)
data(meuse)
meuse.sf <- st_as_sf(meuse,coords=c('x','y'))
w <- 1/as.matrix(st_distance(meuse.sf))
diag(w) <- 0
w <- mat2listw(w)
moran(meuse$copper,w,n=length(w$neighbours),S0=Szero(w))

$I
[1] 0.08963843

$K
[1] 4.277607

Autocorrelación espacial

Indice de Moran

Si miramos en los datos meuse en la variable de concentración de copper

moran.mc(meuse$copper,w,nsim=10000)


    Monte-Carlo simulation of Moran I

data:  meuse$copper 
weights: w  
number of simulations + 1: 10001 

statistic = 0.089638, observed rank = 10001, p-value = 9.999e-05
alternative hypothesis: greater

Autocorrelación espacial

Indice de Moran

Si miramos en los datos meuse en la variable de concentración de copper

moran.plot(meuse$copper,w)

Autocorrelación espacial

Indice de Moran

Si miramos en los datos meuse en la variable de concentración de zinc

moran.mc(meuse$zinc,w,nsim=10000)


    Monte-Carlo simulation of Moran I

data:  meuse$zinc 
weights: w  
number of simulations + 1: 10001 

statistic = 0.10101, observed rank = 10001, p-value = 9.999e-05
alternative hypothesis: greater

Autocorrelación espacial

Indice de Moran

Si miramos en los datos meuse en la variable de concentración de zinc

moran.plot(meuse$zinc,w)

Autocorrelación espacial

Indice de Moran

Si miramos en los datos meuse en la variable de concentración de cadmium

moran.mc(meuse$cadmium,w,nsim=10000)


    Monte-Carlo simulation of Moran I

data:  meuse$cadmium 
weights: w  
number of simulations + 1: 10001 

statistic = 0.079836, observed rank = 10001, p-value = 9.999e-05
alternative hypothesis: greater

Autocorrelación espacial

Indice de Moran

Si miramos en los datos meuse en la variable de concentración de cadmium

moran.plot(meuse$cadmium,w)

Autocorrelación espacial

Indice de Moran local

Generalmente hay interés en proporcionar medidas de autocorrelación local. Ej. para ver como se comporta una observación respecto a su entorno.

Local Indicators of Spatial Association (LISA)

lmoran <- localmoran(meuse$zinc, w, alternative = "greater")
head(lmoran)

            Ii          E.Ii       Var.Ii       Z.Ii Pr(z > E(Ii))
1  0.035611908 -0.0019483794 7.220805e-04  1.3977712    0.08109089
2  0.035940322 -0.0030238443 1.072375e-03  1.1898495    0.11705278
3  0.006105341 -0.0001921442 5.414280e-05  0.8558484    0.19604083
4 -0.001309196 -0.0002833898 6.285942e-05 -0.1293839    0.55147305
5  0.006527050 -0.0002858155 8.101434e-05  0.7569181    0.22454947
6  0.007118978 -0.0002329552 5.216030e-05  1.0179618    0.15434805

Autocorrelación espacial

Indice de Moran local

meuse_sf <- st_as_sf(meuse[c('x','y','zinc')],coords = c('x','y'),crs = 28992)
meuse_sf$Ii <- lmoran[,'Ii']
meuse_sf$Z <- lmoran[,'Z.Ii']
meuse_sf$Pr <- lmoran[, "Pr(z > E(Ii))"]

library(tmap)
p1 <- tm_shape(meuse_sf) +
  tm_dots(col = "zinc", title = "Zinc", style = "quantile") +
  tm_layout(legend.outside = TRUE)

p2 <- tm_shape(meuse_sf) +
  tm_dots(col = "Ii", title = "Índice de Moran Local",
              style = "quantile") +
  tm_layout(legend.outside = TRUE)

p3 <- tm_shape(meuse_sf) +
  tm_dots(col = "Z", title = "Z-score",
              breaks = c(-Inf, 1.65, Inf)) +
  tm_layout(legend.outside = TRUE)

p4 <- tm_shape(meuse_sf) +
  tm_dots(col = "Pr", title = "p-value",
              breaks = c(-Inf, 0.05, Inf)) +
  tm_layout(legend.outside = TRUE)

Autocorrelación espacial

Indice de Moran local

tmap_mode('view')
tmap_arrange(p1,p2,p3,p4)

Autocorrelación espacial

Indice de Moran local

quad <- attr(lmoran,'quadr') 
head(quad)

      mean   median     pysal
1 High-Low High-Low High-High
2 High-Low High-Low High-High
3 High-Low High-Low High-High
4  Low-Low  Low-Low  Low-High
5  Low-Low  Low-Low   Low-Low
6  Low-Low  Low-Low   Low-Low

Autocorrelación espacial

Indice de Moran local

library(dplyr)
meuse_sf <- meuse_sf |> 
  bind_cols(quad = quad$median) |>
  mutate(quad = case_when(
    Pr < 0.05 ~ quad,
    .default = 'No Significativo'))