Introducción

¿Cómo hacer mapas en R?

  • Una de las formas de hacer mapas es R es usando arhivos .shp o shapefiles
library(rgdal) #readOGR
library(RColorBrewer)
library(classInt)
poligonos <- readOGR("nxprovincias.shp",layer="nxprovincias")
## OGR data source with driver: ESRI Shapefile 
## Source: "nxprovincias.shp", layer: "nxprovincias"
## with 25 features
## It has 7 fields
  • ¿Con qué tipo de objeto estamos tratando?
#clase:
class(poligonos)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
#tipo:
typeof(poligonos)
## [1] "S4"
#estructura:
str(poligonos@polygons[1])
## List of 1
##  $ :Formal class 'Polygons' [package "sp"] with 5 slots
##   .. ..@ Polygons :List of 1
##   .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
##   .. .. .. .. ..@ labpt  : num [1:2] 704263 9667236
##   .. .. .. .. ..@ area   : num 8.33e+09
##   .. .. .. .. ..@ hole   : logi FALSE
##   .. .. .. .. ..@ ringDir: int 1
##   .. .. .. .. ..@ coords : num [1:14239, 1:2] 770262 770301 770358 770416 770423 ...
##   .. ..@ plotOrder: int 1
##   .. ..@ labpt    : num [1:2] 704263 9667236
##   .. ..@ ID       : chr "0"
##   .. ..@ area     : num 8.33e+09
  • Puedo indexar el objeto, por ejeplo, removemos galápagos y no delimitadas:
poligonos <- poligonos[poligonos$DPA_PROVIN !=20,]
poligonos <- poligonos[poligonos$DPA_PROVIN !=90,]
  • Creo acrónimos de provincias:
AcrProv <- c("AZ",  "BO",   "C?",   "CA",   "CO",   "CH",   
             "EL",  "ES",   "GU",   "IM",   "LO",   "LR",   
             "MA",  "MO",   "NA",   "PA",   "PI",   "TU",   
             "ZA",  "SU",   "OR",   "SD",   "SE")
  • Calculo los centroides y grafico:
centroides <- coordinates(poligonos)
plot(poligonos)
text(centroides,AcrProv,cex=0.7)

  • Leo los datos que deseo graficar y los adjunto a polígonos:
datos <- read.csv("datosCPV10.csv",header=TRUE,sep=";")
datos$Cod <- poligonos@data$DPA_PROVIN
str(datos)
## 'data.frame':    23 obs. of  10 variables:
##  $ Cod                   : Factor w/ 25 levels "01","02","03",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ Nombre.de.la.Provincia: Factor w/ 23 levels "AZUAY","BOLIVAR",..: 1 2 3 4 6 5 7 8 9 10 ...
##  $ Indig                 : int  17638 46719 34213 5649 90437 174211 4060 15022 46241 102640 ...
##  $ Afro                  : int  10838 1206 4162 6767 4833 3586 24152 123076 204271 12154 ...
##  $ Negrx                 : int  890 160 513 1675 375 212 4676 56571 36434 4051 ...
##  $ Mulatx                : int  3924 581 1277 2120 1605 1162 12613 54864 111372 5221 ...
##  $ Montubix              : int  2941 2067 2399 445 7266 1182 16858 13017 410991 1196 ...
##  $ Mestizx               : int  637912 127795 172616 142933 294840 267880 489843 238619 2461749 261684 ...
##  $ Blancx                : int  36672 4921 9602 4711 9349 9975 46801 31333 355284 10776 ...
##  $ Otrx                  : int  1312 192 402 224 500 373 1656 1590 19141 522 ...
  • Calculo la proporción de indígenas:
indig <- datos[,3]/apply(datos[,3:10],1,sum)
  • Preparo los datos a ser graficados:
indig <- as.data.frame(cbind(1:23))
names(indig) <- "indig"
row.names(indig) <- row.names(poligonos)
poligonos.data <- SpatialPolygonsDataFrame(poligonos,indig)
  • Fijo la paleta de colores:
plotvar <- poligonos.data$indig
nclr <- 5 # Numero de colores
plotclr <- brewer.pal(nclr,"Blues")
class <- classIntervals(round(plotvar,3)*100,nclr,style="equal")
colcode <- findColours(class,plotclr) 
  • Genero el mapa final:
plot (poligonos.data, col=colcode, border="grey", axes=T)
title(main = "Proporción de Indígenas",cex=3)
legend("bottomright",legend = names(attr(colcode,"table")),
       fill= attr(colcode,"palette"),cex=1.25)
text(centroides,AcrProv,cex=1)

leaflet

library(leaflet)
leaflet() %>%
  addTiles() %>%  
  addMarkers(lng=-78.4876692, lat=-0.1962967, popup="¡Todo ojo lo verá!")

Potencialidades: (https://rstudio.github.io/leaflet/shiny.html)

Autocorrelación espacial

Vecino: Dos áreas serán vecinas si dos áreas de una partición comparten un borde de longitud no nula.

La siguiente figura ilustra algunos conceptos:

library(spdep)
library(RColorBrewer)
library(classInt)
library(maptools)

setwd("~/Documents/Consultorias&Cursos/DataLectures/SpatialData")

pr <- readShapePoly("PuertoRico_SPCS.shp")
pr.nb <- read.gal("PuertoRico.txt")
pr.listw <- nb2listw(pr.nb, style="B")

# Plotting Spatal connectivity 
cent <- coordinates(pr)
plot(pr, lwd=1.5,sub = "Partición de las municipalidades de Puerto Rico. 
     \nLas líneas en rojo denotan la red topológica que \n
     conecta los centroides de las municipalidades con bordes comunes.")
plot(pr.nb, cent, add=T, col="red")

Índices para medir dependencia espacial

Coeficiente de Moran (MC) y Ratio de Geary (GR).

Coeficiente de Moran (MC)

  • Es una medida de covarianza que se relaciona directamente con el coeficiente de correlación de Pearson
  • Puede trabajar con todas las escalas de medición: nominal, ordinal, inervalo y ratio.

Consideremos el coeficiente de correlación \(r\) para dos variables \(X\) y \(Y\):

\[ r = \frac{\sum_{i = 1}^{n}1(x_i-\bar{x})(y_i-\bar{y})/n}{\sqrt{\sum_{i = 1}^{n}(x_i-\bar{x})^2/n}\sqrt{\sum_{i = 1}^{n}(y_i-\bar{y})^2/n}}\label{eq:001} \]

donde \(1\) es la frecuencia de la observación \(i\); \(x_i\) y \(y_i\) son los valores de las parejas de dos variables para la observación \(i\). \(\bar{x}\) y \(\bar{y}\) son los promedios de las variables. El numerador es un término de covarianza; el denominador es el producto de las desviaciones estándar de dos variables; y la división para \(n\) se puede reemplazar por \(n-1\) para eliminar el sesgo en muestras pequeñas.

Consideremos ahora el MC muestral para una variable \(Y\):

\[ MC = \frac{\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}(y_i-\bar{y})(y_j-\bar{y})/\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}}{\sqrt{\sum_{i = 1}^{n}(y_i-\bar{y})^2/n}\sqrt{\sum_{i = 1}^{n}(y_i-\bar{y})^2/n}}\label{eq:002} \]

donde \(c_{ij}\) es el valor correspondiente (0 o 1) en la matriz \(\mathbf{C}\). La frecuencia \(c_{ij}\) en \(\eqref{eq:002}\) reemplaza la frecuencia denotada por \(1\) en \(\eqref{eq:001}\).

Justo como \(n\) cuenta el número de 1s en el numerador de \(\eqref{eq:001}\), \(\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}\) cuenta el número de 1s en \(\eqref{eq:002}\).

Recuerda, la autocorrelación se refiere a la correlación dentro de una misma variable, por lo que se aprecia el cambio de \(x_i\) y \(\bar{x}\) por \(y_i\) y \(\bar{y}\) entre las ecuaciones \(\eqref{eq:001}\) y \(\eqref{eq:002}\) respectivamente.

La ecuación \(\eqref{eq:002}\) suele presentarse en la siguiente forma simplificada:

\[ MC = \frac{n}{\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}}\frac{\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}(y_i-\bar{y})(y_j-\bar{y})}{\sum_{i = 1}^{n}(y_i-\bar{y})^2}\label{eq:003} \]

pero su relación con \(\eqref{eq:001}\) es menos notoria.

Propiedades de MC

  • El rango del MC no es \([-1,1]\) (este intervalo puede contraerse, pero usualmente se expande), y el valor central, que denota cero autocorrelación espacial es \(-1/(n-1)\) en lugar de \(0\).

  • El MC para un teselado regualar (por ejemplo, los pixeles de una imagen) convergen asintóticamente al intervalo \([-1,1]\) cuando hay un número creciente de cuadrados.

  • Una partición irregular (por ejemplo, las provincias de Ecuador) no necesariamente tienen la propiedad de convergencia anterior.

  • El valor central converge asintóticamente a \(0\) a mayor número de posiciones.

  • La varianza (para variables normales o no normales simétricas) de MC (para valores \(n\geq25\) cuando es simétrica y \(n\geq100\) cuando no) es aproximadamente1 (Griffith (2010)):

\[ Var(MC) = \frac{2}{\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}} \]

La cual es análoga a la varianza del coeficiente de Pearson transormado \(1/(n-3)\). El \(2\) aparece en el numerador porque el cálculo involucra a \(c_{ij}\) y \(c_{ji}\).

Ratio de Geary (GC)

  • Es una medida de comparación de parejas que se relaciona directamente al gráfico de semivariograma utilizado en geoestadística (tranqui, veremos qué es un semivariograma más adenate).

  • Puede ser escrito en términos del MC:

\[ GR = \frac{n-1}{2\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}}\frac{\sum_{i = 1}^{n}(x_i-\bar{x})^2\left(\sum_{j = 1}^{n}c_{ij}\right)}{\sum_{i = 1}^{n}(x_i-\bar{x})^2}-\frac{n-1}{n}MC \] lo que enfatiza la relación inversa entre los dos índices y subraya que GR es más sensible a la presencia de outliers. ¿Por qué?

  • Si \(MC+GR\approx 1\) entonces los datos tienen buen comportamiento (no hay outliers extremos).

  • El error estándar aproximado del GR es:

\[ GR_{sd} =\sqrt{ \frac{2}{\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}}+2\frac{\sum_{i=1}^n\left( \sum_{j=1}^nc_{ij}\right)^2}{\left(\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}\right)^2} } \]

  • Lo que indica que la distribución de GR es más variable que la de MC. Este es uno de los motivos por los que MC es preferido a GR.

Gráfico de Moran y algunas herramientas:

El gráfico de Moran es un diagrama de dos dimensiones que usa coordenadas cartesianas para mostrar parejas de valores de modo que se resuma la relación entre obervaciones que comprenden un conjunto de datos georreferenciados.

  • La construcción de este gráfico se puede apreciar re escribiendo la ecuación \(\eqref{eq:003}\):

\[ MC = \frac{1}{\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}}\sum_{i=1}^n\left(\frac{(y_i-\bar{y})}{\sqrt{\sum_{i = 1}^{n}(y_i-\bar{y})^2/n}}\right)\left(\sum_{i = 1}^{n}c_{ij}\frac{(y_j-\bar{y})}{\sqrt{\sum_{i = 1}^{n}(y_i-\bar{y})^2/n}}\right) \] el cual se refiere a las parejas ordenadas \((z_i,\sum_{j = 1}^{n}c_{ij}z_j)\). En otras palabras:

  • Paso 1: Convertir cada \(y_i\) en un puntaje \(z\).
  • Paso 2: Graficar cada puntaje \(z\) (eje horizontal) contra su correspondiente suma de puntaje \(z\) (eje vertical).

La línea de tendencia de este gráfico es el MC no estandarizado (esto es, la pendiente necesita ser dividida para \(\sum_{i = 1}^{n}\sum_{j = 1}^{n}c_{ij}\))

En R:

  • Creamos un teselado regular:
x <- c(1,2,2,1,1)
y <- c(1,1,2,2,1)

vl1 <- cbind(x,y);vl2 <- cbind(x-1,y+1);vl3 <- cbind(x,y+1);vl4 <- cbind(x+1,y+1);vl5 <- cbind(x-1,y);vl6 <- cbind(x+1,y);vl7 <- cbind(x-1,y-1);vl8 <- cbind(x,y-1);vl9 <- cbind(x+1,y-1)

Sr1 = Polygon(vl1)
Sr2 = Polygon(vl2)
Sr3 = Polygon(vl3)
Sr4 = Polygon(vl4)
Sr5 = Polygon(vl5)
Sr6 = Polygon(vl6)
Sr7 = Polygon(vl7)
Sr8 = Polygon(vl8)
Sr9 = Polygon(vl9)

Srs1 = Polygons(list(Sr1), "s1")
Srs2 = Polygons(list(Sr2), "s2")
Srs3 = Polygons(list(Sr3), "s3")
Srs4 = Polygons(list(Sr4), "s4")
Srs5 = Polygons(list(Sr5), "s5")
Srs6 = Polygons(list(Sr6), "s6")
Srs7 = Polygons(list(Sr7), "s7")
Srs8 = Polygons(list(Sr8), "s8")
Srs9 = Polygons(list(Sr9), "s9")
SpP = SpatialPolygons(list(Srs1,Srs2,Srs3,Srs4,Srs5,Srs6,Srs7,Srs8,Srs9))
plot(SpP, col = 1:9, pbg="white")

  • Calculamos los vecinos del teselado:
(wr <- poly2nb(SpP, queen=FALSE))
## Neighbour list object:
## Number of regions: 9 
## Number of nonzero links: 24 
## Percentage nonzero weights: 29.62963 
## Average number of links: 2.666667
(wm <- nb2mat(wr, style='B', zero.policy = TRUE))
##    [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
## s1    0    0    1    0    1    1    0    1    0
## s2    0    0    1    0    1    0    0    0    0
## s3    1    1    0    1    0    0    0    0    0
## s4    0    0    1    0    0    1    0    0    0
## s5    1    1    0    0    0    0    1    0    0
## s6    1    0    0    1    0    0    0    0    1
## s7    0    0    0    0    1    0    0    1    0
## s8    1    0    0    0    0    0    1    0    1
## s9    0    0    0    0    0    1    0    1    0
## attr(,"call")
## nb2mat(neighbours = wr, style = "B", zero.policy = TRUE)
  • Graficamos los enlaces entre los polígonos:
plot(SpP, col='gray', border='blue')
xy <- coordinates(SpP)
plot(wr, xy, col='red', lwd=2, add=TRUE)

# Creo lista de vecinos
wrl <- nb2listw(wr, style="B") 
# Número de vecinos por área
i <- rowSums(wm)
vals <- i/sum(i)
moran.plot(vals,wrl)

Otras formas de influencia entre polígonos:

  • Basado en la distancia:
wd10 <- dnearneigh(xy, 0, 10)
wd25 <- dnearneigh(xy, 0, 25, longlat=TRUE)
  • Vecino más cercano:
k3 <- knn2nb(knearneigh(xy, k=3, RANN=FALSE))
k6 <- knn2nb(knearneigh(xy, k=6, RANN=FALSE))
  • Rezago 2:
wr2 <- wr
for (i in 1:length(wr)) {
    lag1 <- wr[[i]]
    lag2 <- wr[lag1]
    lag2 <- sort(unique(unlist(lag2)))
    lag2 <- lag2[!(lag2 %in% c(wr[[i]], i))]
    wr2[[i]] <- lag2
}
  • Graficamos todas las influencias con la función plotit
plotit <- function(nb, lab='') {
  plot(SpP, col='gray', border='white')
  plot(nb, xy, add=TRUE, pch=20)
  text(2.5, 2.8, paste0('(', lab, ')'), cex=1.25)
}
par(mfrow=c(2, 3), mai=c(0,0,0,0))
plotit(wr, 'adjacency')
plotit(wr2, 'lag-2 adj.')
plotit(wd10, '10 km')
plotit(wd25, '25 km')
plotit(k3, 'k=3')
plotit(k6, 'k=6')

Un ejemplo real:

  • Analizamos el número de granjas por área geográfica (municipalidades) en Puerto Rico:
# Farm Density in 2007
farm.den07 <- pr$nofarms_07/pr$area

A continuación:

  • Creamos una paleta con 5 colores.
  • Definimos 5 clases para farm.den07 con la opción quantile.
  • Emparejamos los colores con cada valor de farm.den07
  • Graficamos las municipalidades con los colores emparejados:
# Mapping Farm Density in 2007
#pal.gray <- gray.colors(5)
pal.red <- brewer.pal(5,"Reds")
q5.den <- classIntervals(farm.den07,5,style="quantile") 
cols.den <- findColours(q5.den, pal.red)
plot(pr, col=cols.den)
brks.den <- round(q5.den$brks,3)
leg.txt  <- paste(brks.den[-6], brks.den[-1], sep=" - ")
legend("bottomright", fill=attr(cols.den,"palette"), legend=leg.txt ,bty="n")

  • Calculamos el coeficiente de Moran según la ecuación y el de Geary:
moran(farm.den07, pr.listw, length(pr.listw$weights),Szero(pr.listw),zero.policy = FALSE)
## $I
## [1] 0.3342575
## 
## $K
## [1] 3.255188
  • Realizamos las pruebas de autocorrelación espacial para la densidad de las granjas con el coeficiente de Moran y el ratio de Geary
# Spatial autocorrelation tests
moran.test(farm.den07, pr.listw)
## 
##  Moran I test under randomisation
## 
## data:  farm.den07  
## weights: pr.listw  
## 
## Moran I statistic standard deviate = 4.7987, p-value = 7.986e-07
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.334257491      -0.013888889       0.005263585
geary.test(farm.den07, pr.listw)
## 
##  Geary C test under randomisation
## 
## data:  farm.den07 
## weights: pr.listw 
## 
## Geary C statistic standard deviate = 2.4637, p-value = 0.006876
## alternative hypothesis: Expectation greater than statistic
## sample estimates:
## Geary C statistic       Expectation          Variance 
##       0.773159440       1.000000000       0.008477634
  • Realizamos el gráfico de Moran
      1. La concentración de datos en el primer y cuarto cuadrante muestran autocorrelación espacial positiva. Es decir: Valores altos son cercanos a valores altos, valores intermedios a intermedios y bajos a bajos.
      1. Concentración en el segundo y cuarto cuadrante indican autocorrelación espacial negativa. Esto es: valores altos son cercanos a valores bajos, intermedios a intermedios y bajos a altos.
    • Sea en el caso 1 o 2, valores que quedan en los otros cuadrantes indican heterogeneidad que puede ser detectada con estadísticos locales de Moran.
wr <- poly2nb(pr, queen=FALSE)
wm <- nb2mat(wr, style='B', zero.policy = TRUE)
plot(farm.den07,wm%*%farm.den07) # A mano

# Moran Scatterplot
moran.plot(farm.den07, pr.listw, pch=20)

  • Cargamos el paquete gstat
  • Calculamos el variograma empírico
  • Ajustamos el variograma empírico con el teórico (asumiendo un modelo esférico)
  • Graficamos ambos variogramas (en puntos el empírico y en línea el teórico)

El semi-variogama es un diagrama bi-dimensional donde en el eje x ss muestra la distancia y en el eje y se tiene

\[ \hat\gamma(\mathbf{h}) = \frac{1}{2\#N(\mathbf{h})}\sum(Z(s_i+\mathbf{h})-Z(s_i))^2 \]

# Variogram
library(gstat)
pr.v <- variogram(farm.den07 ~ 1,  pr)
pr.v.fit <- fit.variogram(pr.v, vgm(11,"Sph", "30000", 1))
plot(pr.v, pr.v.fit)

¿Por qué tener en cuenta la autocorrelación espacial?

  • Calculamos la densidad de las granjas en el 2002
# Farm Density in 2002
farm.den02 <- pr$nofarms_02/pr$area
  • Realizamos una regresión lineal de la densidad 2007 y la 2002. Reportamos los resultados y realizamos un test sobre los residuos:
# Run a regression model
lm.farm <- lm(farm.den07 ~ farm.den02)
summary(lm.farm)
## 
## Call:
## lm(formula = farm.den07 ~ farm.den02)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.5157 -1.1769 -0.5443  0.5473  7.0563 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.56014    0.31353   4.976 4.37e-06 ***
## farm.den02   0.61161    0.04692  13.035  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.849 on 71 degrees of freedom
## Multiple R-squared:  0.7053, Adjusted R-squared:  0.7011 
## F-statistic: 169.9 on 1 and 71 DF,  p-value: < 2.2e-16
lm.morantest(lm.farm,pr.listw)
## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = farm.den07 ~ farm.den02)
## weights: pr.listw
## 
## Moran I statistic standard deviate = 2.0715, p-value = 0.01916
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##       0.12765978      -0.02096606       0.00514772

Test de Moran: la hipótesis nula establece que el atributo que se analiza está distribuido en forma aleatoria entre las entidades del área de estudio.

  • Obtenermos los residuos de la regresión
  • Graficamos el mapa
# Mapping residuals of the regression model
res <- resid(lm.farm)
q5.res <- classIntervals(res,5,style="quantile")
cols.res <- findColours(q5.res, pal.red)
plot(pr, col=cols.res)
brks.res <- round(q5.res$brks,3)
leg.txt  <- paste(brks.res[-6], brks.res[-1], sep=" - ")
legend("bottomright", fill=attr(cols.res,"palette"), legend=leg.txt ,bty="n")

Outline

Un curso completo de estadística espacial en R estaría compuesto de los siguientes elementos:

Referencias

Cliff, AD, and J Keith Ord. 1981. Spatial Processes: Models and Applications. London: Pion.

Cliff, Andrew D, and J Keith Ord. 1973. “Spatial Autocorrelation, Monographs in Spatial Environmental Systems Analysis.” London: Pion Limited.

Geary, Robert C. 1954. “The Contiguity Ratio and Statistical Mapping.” The Incorporated Statistician 5 (3). JSTOR: 115–46.

Griffith, Daniel A. 2010. “The Moran Coefficient for Non-Normal Data.” Journal of Statistical Planning and Inference 140 (11). Elsevier: 2980–90.

Moran, Patrick AP. 1948. “The Interpretation of Statistical Maps.” Journal of the Royal Statistical Society. Series B (Methodological) 10 (2). JSTOR: 243–51.


  1. A. Cliff and Ord (1981) desarrollan una fórmula exacta (pero complicada) para la varianza de MC bajo supuestos de normalidad.