All data is spatial data because everithing happen somewhere

0.1 Introducción

Todo sucede en alguna parte, y cada vez más el lugar donde todas estas cosas suceden está siendo registrado en una base de datos. Hay algo de verdad detrás de la afirmación, a menudo repetida, de que el 80% de los datos tienen un componente espacial. Entonces, ¿qué podemos hacer con estos datos espaciales? Estadísticas espaciales, por supuesto! La ubicación es una variable explicativa importante en tantas cosas - ya sea un brote de enfermedad, la elección del hábitat de un animal, un choque de tráfico o una veta de oro en las montañas - que sería prudente incluirla siempre que sea posible.

Es parte de la naturaleza humana tratar de descubrir patrones de un conjunto de eventos aparentemente arbitrarios. Se nos enseña desde una edad temprana a “conectar los puntos”, aprendiendo que si conectamos los puntos correctos de la manera correcta, surgirá una imagen significativa.

En los estudios cientificos, los métodos formalizados para “conectar los puntos” proporcionan herramientas poderosas para identificar asociaciones y patrones entre los resultados y sus supuestas causas. En salud pública, identi???cation y quanti???cation de los patrones de ocurrencia de la enfermedad proporcionan los primeros pasos hacia una mayor comprensión y posiblemente, el control de esa enfermedad en particular.

Como componente del patrón observado, el lugar donde ocurre un evento puede proporcionar alguna indicación de por qué ocurre ese evento en particular. Los métodos estadísticos espaciales ofrecen un medio para que utilicemos dicha información de localización para detectar y cuantificar patrones en los datos de salud pública e investigar el grado de asociación entre los factores de riesgo potenciales y la enfermedad.

0.1.1 Primera ley de la geografía

La primera ley de la geografía, o principio de autocorrelación espacial, fue formulado por el geógrafo Waldo Tobler de esta forma:

Todas las cosas están relacionadas entre sí, pero las cosas más próximas en el espacio tienen una relación mayor que las distantes.

Esta observación se encaja en el modelo de la gravedad de la distribución de viajes o en la ley de la demanda, que indica que las interacciones entre los lugares son inversamente proporcionales al coste del recorrido entre ellos. También se relaciona con las ideas de la ley de Isaac Newton sobre la gravitación universal y viene a ser un concepto similar al de dependencia espacial que constituye la base de la geoestadística.

  • GeoEstadistica: Desarrollada originalmente para la industria minera, la geoestadística cubre el análisis de datos de medición basados en la ubicación. Permite la interpolación basada en modelos de mediciones con estimación de la incertidumbre. Medidas tomadas en unas localizaciones, monitores de polución.

0.2 R y estadística espacial

160+ paquetes en CRAN Task View: Analysis of Spatial Data https://cran.r-project.org/web/views/Spatial.html

  • Clases para datos espaciales (y espacio-temporales)
  • Importación y exportación de datos
  • Análisis exploratorio de datos espaciales
  • Soporte para operaciones en vector y raster
  • Estadística espacial
  • Visualización en gráficos estáticos o interactivos (web)
  • Integración con software SIG

Hay una gran cantidad de paquetes de R útiles para trabajar con datos espaciales. Para los ejemplos de este documento se requiere tener instalados los siguientes:

# install.packages(c("rgdal", "raster", "mapview", "ggplot2", "rgl", "spdep", "caret", "tmap", "geospt", "twitteR", "MODISTools"))
# install.packages("rgdal")
# install.packages("raster")

0.2.1 Cargar los datos

a. Importar capas vector

En R podemos importar archivos vector con el paquete [rgdal] como objetos de clase Spatial___DataFrame:

setwd("D:\\DriveW7\\2019\\Sabana\\Estadistica aplicada Salud Publica\\10. Intro Estadistica Espacial\\data")
library(rgdal)
miShp <- readOGR("llanos.shp", layer = "llanos")
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\DriveW7\2019\Sabana\Estadistica aplicada Salud Publica\10. Intro Estadistica Espacial\data\llanos.shp", layer: "llanos"
## with 1 features
## It has 5 fields
## Integer64 fields read as strings:  cat
summary(miShp)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
##          min        max
## x -73.890084 -67.409803
## y   2.581074   7.104381
## Is projected: FALSE 
## proj4string :
## [+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0]
## Data attributes:
##  cat       Country   ecoregion      AREA             PERIMETER      
##  3:1   Colombia:1   Llanos:1   Min.   :1.533e+11   Min.   :2792892  
##                                1st Qu.:1.533e+11   1st Qu.:2792892  
##                                Median :1.533e+11   Median :2792892  
##                                Mean   :1.533e+11   Mean   :2792892  
##                                3rd Qu.:1.533e+11   3rd Qu.:2792892  
##                                Max.   :1.533e+11   Max.   :2792892
# Ejemplo KML
# myKML <- readOGR("data/llanos_kml.kml", layer = "llanos_kml")

Shapefiles también pueden ser importados con el comando shapefile del paquete [raster] de manera más simple:

library(raster)
shp <- shapefile("deptR2.shp")
shp
## class       : SpatialPolygonsDataFrame 
## features    : 34 
## extent      : 449993.2, 1802275, 27685.35, 1872411  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 3
## names       : DANE, NOM_DPTO,      Otros 
## min values  :   11, AMAZONAS, 0.00000000 
## max values  :   99,  VICHADA, 0.50000000
spplot(shp)

Podemos generar mapas web para visualizar información geográfica usando paquetes como [leaflet] o [mapview]. Con [mapview] solo necesitamos una línea de código para crear un mapa web con los elementos básicos:

# install.packages("mapview")
shp <- shapefile("llanos.shp")
library(mapview)
mapView(shp)


b. Importar capas raster

La importación de archivos raster en R se puede realizar con el paquete [raster]. Archivos de una sola capa se importan con el comando raster, mientras que archivos con múltiples capas se pueden importar con los comandos stack o brick:

miTIFF <- raster("cropimg.tif")
miTIFF
## class       : RasterLayer 
## dimensions  : 6949, 7215, 50137035  (nrow, ncol, ncell)
## resolution  : 0.0002694946, 0.0002694946  (x, y)
## extent      : -73.89029, -71.94589, 3.403178, 5.275895  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## data source : D:\DriveW7\2019\Sabana\Estadistica aplicada Salud Publica\10. Intro Estadistica Espacial\data\cropimg.tif 
## names       : cropimg 
## values      : 0, 8  (min, max)
plot(miTIFF)

El paquete mapview ofrece algunos métodos avanzados de visualización, por ejemplo la vista en slide (correr en la consola):

data(poppendorf)
stck1 <- subset(poppendorf, c(3, 4, 5))
stck2 <- subset(poppendorf, c(2, 3, 4))
slideView(stck1, stck2)

También podemos generar una vista de los datos en cubo 3D (correr en la consola y abrir en una nueva ventana):

kili_data <- system.file("extdata", "kiliNDVI.tif", package = "mapview")
kiliNDVI <- stack(kili_data)
cubeView(kiliNDVI)


c. Obtener datos georreferenciados usando APIs

Actualmente es posible acceder a datos desde aplicaciones web como Facebook o Twitter, entre otras, haciendo uso de sus APIs. Puedes ver un ejemplo de una app web para recuperar tweets geocodificados

https://amsantac.shinyapps.io/twitter-r/

0.3 Explorar

a. Estadísticas descriptivas

Después de haber importado los objetos espaciales podemos llevar a cabo análisis típicos de estadística descriptiva. Con una capa vector, por ejemplo, accedemos a los datos de la tabla de atributos de la siguiente forma:

vis <- shapefile("ba_LC80940792015255LGN00_sr.shp")@data
summary(vis[, 13:20])
##       red              nir             swir1            swir2       
##  Min.   :0.0889   Min.   :0.1693   Min.   :0.2504   Min.   :0.1931  
##  1st Qu.:0.1580   1st Qu.:0.2298   1st Qu.:0.3598   1st Qu.:0.2923  
##  Median :0.1756   Median :0.2510   Median :0.3831   Median :0.3179  
##  Mean   :0.1794   Mean   :0.2545   Mean   :0.3853   Mean   :0.3217  
##  3rd Qu.:0.1983   3rd Qu.:0.2768   3rd Qu.:0.4157   3rd Qu.:0.3495  
##  Max.   :0.2956   Max.   :0.3764   Max.   :0.5063   Max.   :0.4425  
##       dvi               evi               gdvi             ndvi       
##  Min.   :0.03440   Min.   :0.04770   Min.   :0.0937   Min.   :0.0837  
##  1st Qu.:0.06557   1st Qu.:0.08317   1st Qu.:0.1406   1st Qu.:0.1428  
##  Median :0.07760   Median :0.10071   Median :0.1582   Median :0.1794  
##  Mean   :0.07507   Mean   :0.09942   Mean   :0.1558   Mean   :0.1763  
##  3rd Qu.:0.08485   3rd Qu.:0.11388   3rd Qu.:0.1717   3rd Qu.:0.2035  
##  Max.   :0.11980   Max.   :0.16793   Max.   :0.2114   Max.   :0.3415

b. Gráficos exploratorios

Gráficos como diagramas de dispersión, boxplots, o gráficos de barras, pueden ser muy útiles para entender las características de nuestros datos:

library(ggplot2)
p <- qplot(vis$evi, vis$ndvi, color = vis$Landtype, shape = vis$Landtype, 
             main = "EVI vs. NDVI", xlab = "EVI", ylab = "NDVI") + geom_point(size = 4) +
    guides(colour = guide_legend("Type"), shape = guide_legend("Type"))
p

c. Autocorrelación espacial

Cuando trabajamos con datos espaciales es fundamental determinar si la autocorrelación espacial es estadísticamente significativa, como lo muestra el ejemplo siguiente donde se calcula el test I de Moran para autocorrelación espacial global:

# install.packages("spdep")
library(spdep)
data(oldcol)
# spplot(oldcol)
coords.OLD <- cbind(COL.OLD$X, COL.OLD$Y)
moran.test(COL.OLD$CRIME, nb2listw(COL.nb, style="B"))
## 
##  Moran I test under randomisation
## 
## data:  COL.OLD$CRIME  
## weights: nb2listw(COL.nb, style = "B")    
## 
## Moran I statistic standard deviate = 6.2116, p-value = 2.622e-10
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.52063815       -0.02083333        0.00759872

De manera similar podemos evaluar la autocorrelación espacial local calculando estadísticos como Gi:

data(getisord)
xycoords <- cbind(xyz$x, xyz$y)
nb30 <- dnearneigh(xycoords, 0, 30)
G30 <- localG(xyz$val, nb2listw(nb30, style="B"))
brks <- seq(-5,5,1)
cm.col <- cm.colors(length(brks)-1)
image(x, y, t(matrix(G30, nrow=16, ncol=16, byrow=TRUE)), breaks=brks, col=cm.col, asp=1)
text(xyz$x, xyz$y, round(G30, digits=1), cex=0.7)
polygon(c(195,225,225,195), c(195,195,225,225), lwd=2)
title(main=expression(paste("Valores del estadístico ", G[i])))


0.4 Modelar

a. Regresión

Cuando detectamos que la autocorrelación espacial de nuestros datos es estadísticamente significativa entonces debemos aplicar modelos de regresión apropiados como los modelos espaciales autorregresivos:

library(spdep)
data(oldcol)
COL.lag.eig <- lagsarlm(CRIME ~ INC + HOVAL, data = COL.OLD, listw = nb2listw(COL.nb, style="W"), method = "eigen")
summary(COL.lag.eig, correlation=TRUE)
## 
## Call:
## lagsarlm(formula = CRIME ~ INC + HOVAL, data = COL.OLD, listw = nb2listw(COL.nb, 
##     style = "W"), method = "eigen")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -37.68585  -5.35636   0.05421   6.02013  23.20555 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 45.079249   7.177346  6.2808 3.369e-10
## INC         -1.031616   0.305143 -3.3808 0.0007229
## HOVAL       -0.265926   0.088499 -3.0049 0.0026570
## 
## Rho: 0.43102, LR test value: 9.9736, p-value: 0.001588
## Asymptotic standard error: 0.11768
##     z-value: 3.6626, p-value: 0.00024962
## Wald statistic: 13.415, p-value: 0.00024962
## 
## Log likelihood: -182.3904 for lag model
## ML residual variance (sigma squared): 95.494, (sigma: 9.7721)
## Number of observations: 49 
## Number of parameters estimated: 5 
## AIC: 374.78, (AIC for lm: 382.75)
## LM test for residual autocorrelation
## test value: 0.31954, p-value: 0.57188
## 
##  Correlation of coefficients 
##             sigma rho   (Intercept) INC  
## rho         -0.14                        
## (Intercept)  0.12 -0.83                  
## INC         -0.05  0.35 -0.61            
## HOVAL       -0.01  0.08 -0.25       -0.44

0.5 Comunicar

Este es uno de los pasos claves en un ciclo de trabajo de Data Science. Si trabajamos con datos espaciales R nos permite en la actualidad generar mapas temáticos, por ejemplo haciendo uso del paquete tmap:

# install.packages("tmap")
library(tmap)
data(land, World)
tm_shape(land, projection="eck4") +
    tm_raster("elevation", breaks=c(-Inf, 250, 500, 1000, 1500, 2000, 2500, 3000, 4000, Inf),
        palette = terrain.colors(9), title="Elevation", auto.palette.mapping=FALSE) +
tm_shape(World) +
    tm_borders("grey20") +
    tm_grid(projection="longlat", labels.size = .5) +
    tm_text("name", size="AREA") +
tm_compass(position = c(.65, .15), color.light = "grey90") +
tm_credits("Eckert IV projection", position = c(.85, 0)) +
tm_style_classic(inner.margins=c(.04,.03, .02, .01), legend.position = c("left", "bottom"),
    legend.frame = TRUE, bg.color="lightblue", legend.bg.color="lightblue",
    earth.boundary = TRUE, space.color="grey90")

Con paquetes como shiny y flexdashboard podemos ahora desarrollar fácilmente aplicaciones web que permitan a los usuarios manipular, visualizar e interactuar con los datos y los resultados de nuestros modelos.

En este ejemplo, una aplicación desarrollada con Shiny, el usuario puede conocer el número de vehículos de buses por rutas para un área urbana: http://shiny.rstudio.com/gallery/bus-dashboard.html

En esta segunda aplicación, desarrollada con Shiny/flexdashboard, el usuario puede explorar la diversidad de los vecindarios para un área metropolitana: [https://walkerke.shinyapps.io/neighborhood_diversity/].

0.5.1 Cartogramas

Un cartograma es un mapa o diagrama que muestra datos de cantidad asociados a sus respectivas áreas, mediante la modificación de los tamaños de las unidades de enumeración. La información es aportada mediante la distorsión de las superficies reales, utilizando cada superficie de enumeración como un símbolo proporcional, el cual aumenta o disminuye en función de los valores correspondientes.