.shp o shapefileslibrary(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
#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
poligonos <- poligonos[poligonos$DPA_PROVIN !=20,]
poligonos <- poligonos[poligonos$DPA_PROVIN !=90,]
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")
centroides <- coordinates(poligonos)
plot(poligonos)
text(centroides,AcrProv,cex=0.7)
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 ...
indig <- datos[,3]/apply(datos[,3:10],1,sum)
indig <- as.data.frame(cbind(1:23))
names(indig) <- "indig"
row.names(indig) <- row.names(poligonos)
poligonos.data <- SpatialPolygonsDataFrame(poligonos,indig)
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)
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)
leafletlibrary(leaflet)
leaflet() %>%
addTiles() %>%
addMarkers(lng=-78.4876692, lat=-0.1962967, popup="¡Todo ojo lo verá!")
Potencialidades: (https://rstudio.github.io/leaflet/shiny.html)
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")
Coeficiente de Moran (MC) y Ratio de Geary (GR).
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.
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}\).
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} } \]
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.
\[ 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:
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}\))
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")
(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)
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:
wd10 <- dnearneigh(xy, 0, 10)
wd25 <- dnearneigh(xy, 0, 25, longlat=TRUE)
k3 <- knn2nb(knearneigh(xy, k=3, RANN=FALSE))
k6 <- knn2nb(knearneigh(xy, k=6, RANN=FALSE))
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
}
plotitplotit <- 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')
# Farm Density in 2007
farm.den07 <- pr$nofarms_07/pr$area
A continuación:
quantile.# 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")
moran(farm.den07, pr.listw, length(pr.listw$weights),Szero(pr.listw),zero.policy = FALSE)
## $I
## [1] 0.3342575
##
## $K
## [1] 3.255188
# 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
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)
gstatEl 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)
# Farm Density in 2002
farm.den02 <- pr$nofarms_02/pr$area
# 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.
# 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")
Un curso completo de estadística espacial en R estaría compuesto de los siguientes elementos:
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.
A. Cliff and Ord (1981) desarrollan una fórmula exacta (pero complicada) para la varianza de MC bajo supuestos de normalidad.↩