La contigüidad espacial se registra mediante una matriz de unos y ceros mediante la cual es posible evaluar si dos unidades terriroriales se encuentran cercana este sí. Una de las razones para revisar los aspectos relacionados con las contigüidades es la construcción de la matriz de pesos espaciales la cual, como se vio, está en función de la distancia que dos unidades registran entre sí.
Si se trata de patrones puntuales las contogüidades se registran normalmente, pero si se están evaluando datos de áreas (o a partir de polígonos) contruir de manera adecauda la matriz de contigüidades permite evaluar mejor las distancias y, por lo tanto, la matriz de pesos espaciales. Dependerá del tipo de ente territorial que se esté analizando. En el caso que se ha venido estudiando (Colombia) cuando se evalúan unidades territoriales como los departamentos no sucede lo mismo cuando el centroide señala el centro del mismo que cuando señala la capital del departamento.
En esta sesión se verá cómo es posible cambiar este aspecto para calcular la matriz de pesos espaciales.
El primer paso es cargar el mapa de departamentos que se tiene en la base de datos:
library(rgdal) # Se cargan los mapasLoading required package: sp
rgdal: version: 1.3-6, (SVN revision 773)
Geospatial Data Abstraction Library extensions to R successfully loaded
Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
Path to GDAL shared files: D:/Documentos/R/win-library/3.5/rgdal/gdal
GDAL binary built with GEOS: TRUE
Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
Path to PROJ.4 shared files: D:/Documentos/R/win-library/3.5/rgdal/proj
Linking to sp version: 1.3-1
departamentos <- readOGR(dsn = 'Mapas', layer = 'DepartamentosVeredas')OGR data source with driver: ESRI Shapefile
Source: "D:\Dropbox (Personal)\UIS\0000 SEMINARIO ESPACIAL\Bases\CENSO AGROPECUARIO\Total_nacional(csv)\Mapas", layer: "DepartamentosVeredas"
with 33 features
It has 1 fields
Warning message in readOGR(dsn = "Mapas", layer = "DepartamentosVeredas"):
"Z-dimension discarded"
names(departamentos)'DPTO_CCDGO'
data <- departamentos
IDs <- departamentos$DPTO_CCDGO
names(data)'DPTO_CCDGO'
Antes de continuar, es útil crear un par de funciones que permitan hacer cambios de escala facilmente:
## Función: Convertir km a grados
km2d <- function(km){
out <- (km/1.852)/60
return(out)
}km2d(100) ## 100 km0.899928005759539
## Función: Convertir grados a km
d2km <- function(d){
out <- d*60*1.852
return(out)
}d2km(1) ## 1 grado111.12
Un primer análisis de las contigüidades:
library(spdep)Loading required package: Matrix
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')`
## Vecindades sin establecer la distancia
W_cont <- poly2nb(data, queen=T)
W_cont_mat <- nb2listw(W_cont, style="W", zero.policy=TRUE)
## Vecindades estableciendo distancia (100 km)
W_cont_s <- poly2nb(data, queen=T, snap=km2d(100))
W_cont_s_mat <- nb2listw(W_cont_s, style="W", zero.policy=TRUE)map_crd <- coordinates(departamentos) # Se extraen las coordenadasY las gráficas quedan:
par(mfrow=c(1,2),mar=c(0,0,1,0))
plot(data,border="grey")
plot(W_cont_mat,coords=map_crd,pch=19, cex=0.1, col="blue", add=T)
title("Vecindad Directa")
plot(data,border="grey")
plot(W_cont_s_mat,coords=map_crd,pch=19, cex=0.1, col="blue", add=T)
title("Vecindad + 100 km")
Como se observa, ampliar la distancia en 100 Km (0.9 grados) cambia las relaciones. Ahora, ¿cómo cambiarán estas relaciones si se considera la capital del departamento y no el centro del polígono?
Para respoderlo, es necesario incorporar esa información al mapa:
library(readr)
capitales <- read_delim("capitales.txt", "\t", escape_double = FALSE, trim_ws = TRUE)
head(capitales)Parsed with column specification:
cols(
MPIO = col_double(),
MPIO_CCDGO = col_double(),
MUNICIPIO = col_character(),
LONGITUD = col_double(),
LATITUD = col_double(),
DPTO_CCDGO = col_character()
)
| MPIO | MPIO_CCDGO | MUNICIPIO | LONGITUD | LATITUD | DPTO_CCDGO |
|---|---|---|---|---|---|
| 1110 | 5001 | MEDELLIN | -75.61103 | 6.257590 | 05 |
| 21 | 8001 | BARRANQUILLA | -74.82772 | 10.981520 | 08 |
| 1100 | 11001 | BOGOTA | -74.18107 | 4.316108 | 11 |
| 986 | 13001 | CARTAGENA DE INDIAS | -75.45890 | 10.463430 | 13 |
| 1066 | 15001 | TUNJA | -73.37802 | 5.518473 | 15 |
| 153 | 17001 | MANIZALES | -75.50726 | 5.083429 | 17 |
Se transforma la variable de interés (DPTO_CCDGO) en factor:
capitales<-as.data.frame(capitales)
capitales$DPTO_CCDGO<-factor(capitales$DPTO_CCDGO)
str(capitales)'data.frame': 33 obs. of 6 variables:
$ MPIO : num 1110 21 1100 986 1066 ...
$ MPIO_CCDGO: num 5001 8001 11001 13001 15001 ...
$ MUNICIPIO : chr "MEDELLIN" "BARRANQUILLA" "BOGOTA" "CARTAGENA DE INDIAS" ...
$ LONGITUD : num -75.6 -74.8 -74.2 -75.5 -73.4 ...
$ LATITUD : num 6.26 10.98 4.32 10.46 5.52 ...
$ DPTO_CCDGO: Factor w/ 33 levels "05","08","11",..: 1 2 3 4 5 6 7 8 9 10 ...
- attr(*, "spec")=
.. cols(
.. MPIO = col_double(),
.. MPIO_CCDGO = col_double(),
.. MUNICIPIO = col_character(),
.. LONGITUD = col_double(),
.. LATITUD = col_double(),
.. DPTO_CCDGO = col_character()
.. )
Ahora, se puede unir todo en un solo data frame:
data <- as.data.frame(departamentos)
merged <- merge(data, capitales, by.x="DPTO_CCDGO", by.y="DPTO_CCDGO", all.x=T, all.y=F)
head(merged)| DPTO_CCDGO | MPIO | MPIO_CCDGO | MUNICIPIO | LONGITUD | LATITUD |
|---|---|---|---|---|---|
| 05 | 1110 | 5001 | MEDELLIN | -75.61103 | 6.257590 |
| 08 | 21 | 8001 | BARRANQUILLA | -74.82772 | 10.981520 |
| 11 | 1100 | 11001 | BOGOTA | -74.18107 | 4.316108 |
| 13 | 986 | 13001 | CARTAGENA DE INDIAS | -75.45890 | 10.463430 |
| 15 | 1066 | 15001 | TUNJA | -73.37802 | 5.518473 |
| 17 | 153 | 17001 | MANIZALES | -75.50726 | 5.083429 |
Organizar un poco los fragmentos de datos:
merged <- merged[order(merged$DPTO_CCDGO),]
map2<-departamentos
rownames(merged) <- map2$DPTO_CCDGO
rownames(map2@data) <- map2$DPTO_CCDGOE in corporporar la data al mapa:
library(dplyr)Attaching package: 'dplyr'
The following objects are masked from 'package:stats':
filter, lag
The following objects are masked from 'package:base':
intersect, setdiff, setequal, union
map2@data <- left_join(map2@data, merged)Joining, by = "DPTO_CCDGO"
head(map2@data)| DPTO_CCDGO | MPIO | MPIO_CCDGO | MUNICIPIO | LONGITUD | LATITUD |
|---|---|---|---|---|---|
| 05 | 1110 | 5001 | MEDELLIN | -75.61103 | 6.257590 |
| 08 | 21 | 8001 | BARRANQUILLA | -74.82772 | 10.981520 |
| 11 | 1100 | 11001 | BOGOTA | -74.18107 | 4.316108 |
| 13 | 986 | 13001 | CARTAGENA DE INDIAS | -75.45890 | 10.463430 |
| 15 | 1066 | 15001 | TUNJA | -73.37802 | 5.518473 |
| 17 | 153 | 17001 | MANIZALES | -75.50726 | 5.083429 |
Extraer las coordenadas para la capitales
map_crd2 <- cbind(map2$LONGITUD, map2$LATITUD) # Coordenadas de las Capitales
head(map_crd2)| -75.61103 | 6.257590 |
| -74.82772 | 10.981520 |
| -74.18107 | 4.316108 |
| -75.45890 | 10.463430 |
| -73.37802 | 5.518473 |
| -75.50726 | 5.083429 |
Crear un ID para poder construir las matrices de pesos espaciales:
data <- map2
IDs <- map2$DPTO_CCDGO
names(data)- 'DPTO_CCDGO'
- 'MPIO'
- 'MPIO_CCDGO'
- 'MUNICIPIO'
- 'LONGITUD'
- 'LATITUD'
Ahora, se construyen las matrices (una para los centroides y una para las capitales). En este ejemplo, se calculan las matrices con una condicion adicional: sólo se tiene en cuenta el primer vecino (k = 1):
library(spdep)#########
## k = 1
#########
## Centroides
W_knn1 <- knn2nb(knearneigh(map_crd, k=1), row.names=IDs)
W_knn1_mat <- nb2listw(W_knn1)
## Capitales
W_knn1_2 <- knn2nb(knearneigh(map_crd2, k=1), row.names=IDs)
W_knn1_mat_2 <- nb2listw(W_knn1_2)
Y el respectivo mapa:
par(mfrow=c(1,2),mar=c(0,0,1,0))
plot(data,border="grey")
plot(W_knn1_mat,coords=map_crd,pch=19, cex=0.1, col="blue", add=T)
title("k=1 (Centroides)")
plot(data,border="grey")
plot(W_knn1_mat_2,coords=map_crd2,pch=19, cex=0.1, col="blue", add=T)
title("k=1 (Capitales)")
¿Que pasa si se cambia el número de vecinos a evaluar?
#########
## k = 5
#########
## Centroides
W_knn1.5 <- knn2nb(knearneigh(map_crd, k=5), row.names=IDs)
W_knn1_mat.5 <- nb2listw(W_knn1.5)
## Capitales
W_knn1_2.5 <- knn2nb(knearneigh(map_crd2, k=5), row.names=IDs)
W_knn1_mat_2.5 <- nb2listw(W_knn1_2.5)par(mfrow=c(1,2),mar=c(0,0,1,0))
plot(data,border="grey")
plot(W_knn1_mat.5,coords=map_crd,pch=19, cex=0.1, col="blue", add=T)
title("k=5 (Centroides)")
plot(data,border="grey")
plot(W_knn1_mat_2.5,coords=map_crd2,pch=19, cex=0.1, col="blue", add=T)
title("k=5 (Capitales)")
¿Qué ocurre, entonces, con la autocorrelación espacial? Hay que recordar que la autocorrelación espacial se evalúa en función de la matriz de pesos y de las distancias. Como se observa en los mapas, las relaciones cambian.
En las sesiones anteriores se ha empleado el cultivo de la Yuca para ilustrar los ejemplos. Se evaluará la autocorrelación espacial con esta nueva configuración. En el ejercicio anterior se usaron los municipios como unidades territoriales; para apreciar si hay cambios o no se debe cambiar el análisis al nivel departamental.
Lo primero es cargar los datos:
load('Total_nacional(csv)/Cultivos.RData')Y seleccionar el cultivo:
Cultivos$P_S6P46<-factor(Cultivos$P_S6P46) # Esto ya se discutió
cultivo1<-Cultivos[Cultivos$P_S6P46=='00159201001',]Organizar algunas variables:
cultivo1$P_DEPTO<-factor(cultivo1$P_DEPTO)
cultivo1$P_MUNIC<-factor(cultivo1$P_MUNIC)
cultivo1$COD_VEREDA<-factor(cultivo1$COD_VEREDA)Se va a emplear el área cosechada para ilustar:
areas<-aggregate(AREA_COSECHADA~P_DEPTO, FUN = sum, data = cultivo1)head(areas)| P_DEPTO | AREA_COSECHADA |
|---|---|
| 05 | 16195.838 |
| 08 | 1083.152 |
| 13 | 44495.876 |
| 15 | 14380.682 |
| 17 | 12544.785 |
| 18 | 14663.987 |
Se cambia el nombre de la variable con la que se van a unir las tablas
areas$DPTO_CCDGO<-areas$P_DEPTO
areas$P_DEPTO<-NULL
head(areas)| AREA_COSECHADA | DPTO_CCDGO |
|---|---|
| 16195.838 | 05 |
| 1083.152 | 08 |
| 44495.876 | 13 |
| 14380.682 | 15 |
| 12544.785 | 17 |
| 14663.987 | 18 |
Se unen las tablas:
data@data <- left_join(data@data, areas)Joining, by = "DPTO_CCDGO"
Warning message:
"Column `DPTO_CCDGO` joining factors with different levels, coercing to character vector"
head(data@data)| DPTO_CCDGO | MPIO | MPIO_CCDGO | MUNICIPIO | LONGITUD | LATITUD | AREA_COSECHADA |
|---|---|---|---|---|---|---|
| 05 | 1110 | 5001 | MEDELLIN | -75.61103 | 6.257590 | 16195.838 |
| 08 | 21 | 8001 | BARRANQUILLA | -74.82772 | 10.981520 | 1083.152 |
| 11 | 1100 | 11001 | BOGOTA | -74.18107 | 4.316108 | NA |
| 13 | 986 | 13001 | CARTAGENA DE INDIAS | -75.45890 | 10.463430 | 44495.876 |
| 15 | 1066 | 15001 | TUNJA | -73.37802 | 5.518473 | 14380.682 |
| 17 | 153 | 17001 | MANIZALES | -75.50726 | 5.083429 | 12544.785 |
Se cambian los datos perdidos por cero (sólo para efectos prácticos)
data$AREA_COSECHADA[is.na(data$AREA_COSECHADA)] <- 0summary(data@data) DPTO_CCDGO MPIO MPIO_CCDGO MUNICIPIO
Length:33 Min. : 21.0 Min. : 5001 Length:33
Class :character 1st Qu.: 745.0 1st Qu.:20001 Class :character
Mode :character Median : 991.0 Median :52001 Mode :character
Mean : 852.7 Mean :52153
3rd Qu.:1066.0 3rd Qu.:81001
Max. :1110.0 Max. :99001
LONGITUD LATITUD AREA_COSECHADA
Min. :-81.72 Min. :-3.530 Min. : 0
1st Qu.:-75.71 1st Qu.: 3.319 1st Qu.: 6239
Median :-74.83 Median : 5.083 Median :10181
Mean :-74.22 Mean : 5.555 Mean :13452
3rd Qu.:-72.96 3rd Qu.: 8.112 3rd Qu.:17260
Max. :-68.14 Max. :12.543 Max. :44496
Ahora, se procede a evaluar la autocorrelación (usando la I de Morán) tomando en cuenta las distintas matrices de pesos espaciales que se construyero
# Normal
moran.test(data$AREA_COSECHADA, listw=W_cont_mat, zero.policy=T) Moran I test under randomisation
data: data$AREA_COSECHADA
weights: W_cont_mat n reduced by no-neighbour observations
Moran I statistic standard deviate = 0.84942, p-value = 0.1978
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.06220558 -0.03225806 0.01236765
# 100 Km
moran.test(data$AREA_COSECHADA, listw=W_cont_s_mat, zero.policy=T) Moran I test under randomisation
data: data$AREA_COSECHADA
weights: W_cont_s_mat n reduced by no-neighbour observations
Moran I statistic standard deviate = 0.84343, p-value = 0.1995
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.035784845 -0.032258065 0.006508253
# Centroides
moran.test(data$AREA_COSECHADA, listw=W_knn1_mat, zero.policy=T) Moran I test under randomisation
data: data$AREA_COSECHADA
weights: W_knn1_mat
Moran I statistic standard deviate = 1.0489, p-value = 0.1471
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.19212952 -0.03125000 0.04535755
# Capitales
moran.test(data$AREA_COSECHADA, listw=W_knn1_mat_2, zero.policy=T) Moran I test under randomisation
data: data$AREA_COSECHADA
weights: W_knn1_mat_2
Moran I statistic standard deviate = -0.59215, p-value = 0.7231
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
-0.1624198 -0.0312500 0.0490682
# centroides k=5
moran.test(data$AREA_COSECHADA, listw=W_knn1_mat.5, zero.policy=T) Moran I test under randomisation
data: data$AREA_COSECHADA
weights: W_knn1_mat.5
Moran I statistic standard deviate = 1.0076, p-value = 0.1568
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.058484843 -0.031250000 0.007931913
# Capitales k=5
moran.test(data$AREA_COSECHADA, listw=W_knn1_mat_2.5, zero.policy=T) Moran I test under randomisation
data: data$AREA_COSECHADA
weights: W_knn1_mat_2.5
Moran I statistic standard deviate = 1.3186, p-value = 0.09365
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
0.086779727 -0.031250000 0.008012538
La autocorrelación local se puede evaluar, ahora, a la luz de un par de matrices:
lm1 <- localmoran(data$AREA_COSECHADA, listw=W_cont_mat, zero.policy=T)
data$lm1 <- abs(lm1[,4]) ## Extract z-scoreslm.palette <- colorRampPalette(c("white","orange", "red"), space = "rgb")spplot(data, zcol="lm1", col.regions=lm.palette(20), main="Local Moran's I (|z| scores)", pretty=T, col="transparent",
par.settings = list(panel.background=list(col="lightblue")))
lm2 <- localmoran(data$AREA_COSECHADA, listw=W_knn1_mat_2.5, zero.policy=T)
data$lm2 <- abs(lm2[,4]) ## Extract z-scoresspplot(data, zcol="lm2", col.regions=lm.palette(20), main="Local Moran's I (|z| scores)", pretty=T,col="transparent",
par.settings = list(panel.background=list(col="lightblue")))
data@data| DPTO_CCDGO | MPIO | MPIO_CCDGO | MUNICIPIO | LONGITUD | LATITUD | AREA_COSECHADA | lm1 | lm2 |
|---|---|---|---|---|---|---|---|---|
| 05 | 1110 | 5001 | MEDELLIN | -75.61103 | 6.257590 | 16195.83837 | 0.503372733 | 0.04427836 |
| 08 | 21 | 8001 | BARRANQUILLA | -74.82772 | 10.981520 | 1083.15221 | 3.712280076 | 2.40682148 |
| 11 | 1100 | 11001 | BOGOTA | -74.18107 | 4.316108 | 0.00000 | 1.079767408 | 0.78675210 |
| 13 | 986 | 13001 | CARTAGENA DE INDIAS | -75.45890 | 10.463430 | 44495.87649 | 2.326270179 | 2.83712870 |
| 15 | 1066 | 15001 | TUNJA | -73.37802 | 5.518473 | 14380.68195 | 0.009676157 | 0.05299484 |
| 17 | 153 | 17001 | MANIZALES | -75.50726 | 5.083429 | 12544.78491 | 0.158177310 | 0.20466535 |
| 18 | 179 | 18001 | FLORENCIA | -75.55824 | 1.749139 | 14663.98651 | 0.045192763 | 0.20051551 |
| 19 | 436 | 19001 | POPAYAN | -76.59194 | 2.471704 | 29607.02923 | 0.512821344 | 0.74905559 |
| 20 | 988 | 20001 | VALLEDUPAR | -73.46326 | 10.382036 | 9039.56282 | 0.580503495 | 0.66780768 |
| 23 | 495 | 23001 | MONTERIA | -75.95055 | 8.584698 | 24211.03392 | 2.368927573 | 2.12289387 |
| 25 | 523 | 25001 | AGUA DE DIOS | -74.67111 | 4.372745 | 2851.22349 | 1.154988956 | 1.84073174 |
| 27 | 935 | 27001 | QUIBDO | -76.66368 | 5.938953 | 30361.32792 | 0.954680075 | 1.63488305 |
| 41 | 380 | 41001 | NEIVA | -75.27236 | 2.993360 | 8198.29907 | 0.267222264 | 0.06114394 |
| 44 | 1102 | 44001 | RIOHACHA | -72.95876 | 11.242972 | 8661.29960 | 0.380907125 | 0.73896740 |
| 47 | 1101 | 47001 | SANTA MARTA | -73.88528 | 11.121894 | 32606.15550 | 0.869482144 | 1.49850836 |
| 50 | 1096 | 50001 | VILLAVICENCIO | -73.49290 | 4.091666 | 14179.13937 | 0.001547429 | 0.02072196 |
| 52 | 1035 | 52001 | PASTO | -77.20610 | 1.083605 | 36607.87119 | 1.464441103 | 0.38885263 |
| 54 | 991 | 54001 | CUCUTA | -72.48863 | 8.112098 | 9582.26717 | 0.211739659 | 0.41440174 |
| 63 | 745 | 63001 | ARMENIA | -75.72490 | 4.499504 | 1631.72042 | 1.122941878 | 1.47614010 |
| 66 | 574 | 66001 | PEREIRA | -75.71236 | 4.781292 | 1502.51927 | 0.044820730 | 1.22276949 |
| 68 | 958 | 68001 | BUCARAMANGA | -73.11157 | 7.155836 | 9026.31188 | 0.392608805 | 0.32567312 |
| 70 | 1109 | 70001 | SINCELEJO | -75.43175 | 9.316674 | 22449.87174 | 2.323640192 | 1.67417727 |
| 81 | 955 | 81001 | ARAUCA | -70.50921 | 6.796281 | 6902.68753 | 0.212522504 | 0.41933695 |
| 85 | 1064 | 85001 | YOPAL | -72.25803 | 5.242745 | 6238.92233 | 0.480490084 | 0.73738437 |
| 86 | 1039 | 86001 | MOCOA | -76.68436 | 1.228575 | 7411.60239 | 1.114968662 | 0.71241032 |
| 88 | 1005 | 88001 | SAN ANDRES | -81.71762 | 12.543117 | 88.75956 | NaN | 3.01421072 |
| 91 | 1007 | 91001 | LETICIA | -70.04513 | -3.530142 | 17259.81546 | 0.240456472 | 0.21304551 |
| 94 | 1054 | 94001 | INIRIDA | -68.45724 | 3.319405 | 20341.50395 | 0.615043984 | 0.84214651 |
| 95 | 1081 | 95001 | SAN JOSE DEL GUAVIARE | -71.91917 | 2.484286 | 3190.93267 | 0.092754181 | 1.78114891 |
| 97 | 1047 | 97001 | MITU | -70.41755 | 1.015572 | 2665.65668 | 0.010370971 | 0.33966620 |
| 99 | 934 | 99001 | PUERTO CARRENO | -68.14122 | 5.836530 | 15039.50987 | 0.027219171 | 0.05652976 |
| 76 | 771 | 76001 | CALI | -76.57649 | 3.399044 | 10729.02856 | 0.011521202 | 0.25401276 |
| 73 | 1098 | 73001 | IBAGUE | -75.25259 | 4.451922 | 10181.25309 | 0.406534370 | 0.71783903 |
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent