DATOS EN ÁREA

Datos en áreas.

Datos espaciales

Cronología del curso estadística espacial año 2023-1S.

La estadística espacial estudia procesos que ocurren en el espacio y de los que, de una manera u otra, se conoce su ubicación. Desde un punto de vista estadístico, esta dependencia basada en la distancia puede modelizarse de distintas maneras utilizando las herramientas proporcionadas por la estadística espacial.

En el presente curso se implementaron distintas herramientas de acuerdo a cada estudio, comenzando con la geoestadística en donde se tenia datos regulares o irregularmente espaciados continuo y fijo una de sus aplicaciones es la presencia o ausencia de una plaga en un cultivo, contenido de nitrógeno en una parcela entre otros ,para su estudio se hizo uso de técnicas del kriging , seguido se estudió los patrones puntuales (discreto y continuo)para ello se implementaron técnicas de pruebas de aleatoriedad , estimaciones de la densidad con la función Kripley, su mayor área de aplicación se centra en los estudios forestales y la ecología. Para finalizar la ramificación de la estadística espacial se abordara el tema de Lattices o también conocido como datos en áreas(discreto y fijo)su área de aplicación se centra en el campo de la sociología, economía, epidemiología , un ejemplo de estas aplicaciones es el estudio de las tasas de morbilidad por departamento, tasas de muertes violentas, número de accidentes en una ciudad.

Datos en área.

En esta unidad 9 del libro Modeling Areal Data of the first edition of Bivand, R. S., Pebesma, E. and Gómez-Rubio V. Applied Spatial Data Analysis with R se abordarán temas relacionados a los datos en áreas o tambien conocidos como lattices.

Introducción.

Los datos espaciales se visualizan a menudo en ventanas poligonales con limites definidos, estos son definidos por el investigador , dependiendo el campo de estudio y la información recopilada estos pueden ser arbitrarios y en otros casos, las entidades regionales pueden constituir en mi misma las unidades de observación .Los datos observados son frecuentemente agregaciones dentro de los límites, como recuentos de población. Un ejemplo bastante interesante es estudiar el comportamiento de los gobiernos locales, en donde los mandatos o leyes se toman a nivel de entidades , por ejemplo el estableciendo de tasas impositivas locales, Sin embargo, en general, las entidades son agregados, que se utilizan para contabilizar las mediciones, como los resultados de las votaciones en los comicios electorales. Muy a menudo, las entidades de área son una teselación exhaustiva del área de estudio, sin dejar ninguna parte del área total sin asignar a una entidad.

Los límites arbitrarios de las unidades de área presentan problemas si su modificación genera resultados distintos, lo cual obstaculiza un análisis adecuado. Es necesario considerar la escala espacial o el alcance del proceso de generación de datos subyacente, ya que estos deben concordar con los límites previamente establecidos. Si es posible rediseñar la recopilación de datos de manera que las entidades de área coincidan con los datos, se reducirá la influencia de la elección de los agregados. Sin embargo, si no se logra esta concordancia, es necesario recurrir a límites diferentes. Por tanto, es de vital importancia tener en cuenta los posibles desafíos que pueden surgir en este contexto.

Si la recopilación de datos puede diseñarse para que las entidades de área coincidan con los datos la influencia de la elección de los agregados se reducirá. En su intercambio con Tyler en 1889, Galton se preguntaba si las observaciones de las leyes matrimoniales en las entidades de área constituían observaciones independientes, ya que podían reflejar simplemente un patrón general del que todos dependen. En el capítulo 8, se mostró cómo las distancias en una superficie continua pueden utilizarse para estructurar la autocorrelación espacial, por ejemplo con el variograma.

La siguiente sección cubrirá la construcción de vecinos y pesos que se puede aplicar a los barrios. Una vez que este importante y a menudo requisito previo está en su lugar, vamos a buscar formas de medir la autocorrelación espacial. Si bien las pruebas se basan en modelos de procesos espaciales, observamos primero las pruebas y solo después pasar al modelado.

Modelando datos en área.

Después de conocer algunas terminologias y origen del Modelado de áreas se puede concluir lo siguiente: Aplicar datos en área es muy útil en datos espaciales en los cuales, partiendo en un dominio continuo, los pesos de ciertas áreas pueden no tener la representación esperada, lo que conyeba a generar otras dependecias entre ciertos vecinos mientras que en otros no. Por ello, se hace necesario delimitar áreas sobre esa región de estudio. Aunque cabe notar que, definir una ventana de observación suele ser en ocasiones difícil y costoso y sino se tiene claridad de la constitución de la región, la definición de una ventana de observación y de otra con los mismos datos pueden generar resultados completamente distintos.

Lo que dificultaria al análisis , puesto que este tendrá muchas interpretaciones y no se tendrá una correcta validación de nuestro modelo.

La creación de ponderaciones espaciales.

El enfoque de esta sección es comprobar que no queda ningun patrón espacial en los residuos.

La creación de ponderaciones espaciales es un paso necesario en el uso de datos de área, quizás sólo para comprobar que no existen patrones espaciales en los residuos. El primer paso consiste en definir las relaciones entre observaciones a las que se va a dar un peso distinto de cero, es decir, elegir el criterio de vecindad.

Para poder crear estos objetos vecinos se requiere instalar algunas librerias , propias para el estudio de datos en área.

Observaciones de los paquetes:

En el paquete spdep , las relaciones de vencidad entre n observaciones se representa mediante un objeto nb.

Si una observación no tiene vecinos el componente contiene un cero entero.

rm(list = ls())
library(rgdal)

## Loading required package: sp

## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-4, (SVN revision 1196)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/57311/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/57311/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.6-0
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.

library(spdep)

## Warning: package 'spdep' was built under R version 4.2.3

## Loading required package: spData

## Warning: package 'spData' was built under R version 4.2.3

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

## Loading required package: sf

## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(spatialreg)

## Warning: package 'spatialreg' was built under R version 4.2.3

## Loading required package: Matrix

## 
## Attaching package: 'spatialreg'

## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library(rgdal)
library(rgeos)

## Warning: package 'rgeos' was built under R version 4.2.3

## rgeos version: 0.6-3, (SVN revision 696)
##  GEOS runtime version: 3.9.3-CAPI-1.14.3 
##  Please note that rgeos will be retired during October 2023,
## plan transition to sf or terra functions using GEOS at your earliest convenience.
## See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
##  GEOS using OverlayNG
##  Linking to sp version: 1.6-0 
##  Polygon checking: TRUE

library(sphet)

## Warning: package 'sphet' was built under R version 4.2.3

## 
## Attaching package: 'sphet'

## The following object is masked from 'package:spatialreg':
## 
##     impacts

library(coda)
library(hierfstat)

## Warning: package 'hierfstat' was built under R version 4.2.3

library(spData)
library(devtools)

## Warning: package 'devtools' was built under R version 4.2.3

## Loading required package: usethis

## Warning: package 'usethis' was built under R version 4.2.3

library(McSpatial)

Cargando SpatialPointsDataframe - SpatialPolygonsDataFrame - SpatialPointsDataFrame Cargando SpatialPointsDataframe - SpatialPolygonsDataFrame - SpatialPointsDataFrame

 NY8 <- readOGR(".", "NY8_utm18")

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\57311\Downloads\DATOS EN AREA LIBRO", layer: "NY8_utm18"
## with 281 features
## It has 17 fields

 NY_nb <- read.gal("NY_nb.gal", region.id =row.names(NY8))
 summary(NY_nb)

## Neighbour list object:
## Number of regions: 281 
## Number of nonzero links: 1522 
## Percentage nonzero weights: 1.927534 
## Average number of links: 5.41637 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 
##  6 11 28 45 59 49 45 23 10  3  2 
## 6 least connected regions:
## 55 97 100 101 244 245 with 1 link
## 2 most connected regions:
## 34 82 with 11 links

TCE <- readOGR(".", "TCE")

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\57311\Downloads\DATOS EN AREA LIBRO", layer: "TCE"
## with 11 features
## It has 5 fields

cities <- readOGR(".", "NY8cities")

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## Warning: OGR support is provided by the sf and terra packages among others

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\57311\Downloads\DATOS EN AREA LIBRO", layer: "NY8cities"
## with 6 features
## It has 1 fields

Para saber como se almacenan los datos realizamos el siguiente comando , este nos arrojará la clasificación de ellos.

class(NY8)

## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"

Se visualiza que es de tipo Polígono y es un dataframe.

getClass("SpatialPolygonsDataFrame")

## Class "SpatialPolygonsDataFrame" [package "sp"]
## 
## Slots:
##                                                                   
## Name:         data    polygons   plotOrder        bbox proj4string
## Class:  data.frame        list     integer      matrix         CRS
## 
## Extends: 
## Class "SpatialPolygons", directly
## Class "Spatial", by class "SpatialPolygons", distance 2
## Class "SpatialPolygonsNULL", by class "SpatialPolygons", distance 2, with explicit coerce

Los polígonos espaciales contienen algunas caracteristicas, para conocerla se hace uso del getclass

getClass("SpatialPolygons")

## Class "SpatialPolygons" [package "sp"]
## 
## Slots:
##                                                       
## Name:     polygons   plotOrder        bbox proj4string
## Class:        list     integer      matrix         CRS
## 
## Extends: "Spatial", "SpatialPolygonsNULL"
## 
## Known Subclasses: 
## Class "SpatialPolygonsDataFrame", directly, with explicit coerce

Un polígono es una secuencia de líneas cerradas;puntos coordenadas donde el primer punto es igual al último.

getClass("Polygon")

## Class "Polygon" [package "sp"]
## 
## Slots:
##                                               
## Name:    labpt    area    hole ringDir  coords
## Class: numeric numeric logical integer  matrix
## 
## Extends: "Line"

labpt - punto de etiqueta, centroide del polígono es una línea es una lista ordenada de coordenadas

getClass("Line")

## Class "Line" [package "sp"]
## 
## Slots:
##              
## Name:  coords
## Class: matrix
## 
## Known Subclasses: "Polygon"

Spatial es la clase madre de todas las clases Spatial usando o en el paquete sp

getClass("Spatial")

## Class "Spatial" [package "sp"]
## 
## Slots:
##                               
## Name:         bbox proj4string
## Class:      matrix         CRS
## 
## Known Subclasses: 
## Class "SpatialPoints", directly
## Class "SpatialMultiPoints", directly
## Class "SpatialGrid", directly
## Class "SpatialLines", directly
## Class "SpatialPolygons", directly
## Class "SpatialPointsDataFrame", by class "SpatialPoints", distance 2
## Class "SpatialPixels", by class "SpatialPoints", distance 2
## Class "SpatialMultiPointsDataFrame", by class "SpatialMultiPoints", distance 2
## Class "SpatialGridDataFrame", by class "SpatialGrid", distance 2
## Class "SpatialLinesDataFrame", by class "SpatialLines", distance 2
## Class "SpatialPixelsDataFrame", by class "SpatialPoints", distance 3
## Class "SpatialPolygonsDataFrame", by class "SpatialPolygons", distance 2

Resumen de los datos para su debida interpretación.

summary(NY8)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##         min       max
## x  358241.9  480393.1
## y 4649755.4 4808545.2
## Is projected: TRUE 
## proj4string :
## [+proj=utm +zone=18 +ellps=WGS84 +units=m +no_defs]
## Data attributes:
##    AREANAME           AREAKEY                X                 Y         
##  Length:281         Length:281         Min.   :-55.482   Min.   :-75.29  
##  Class :character   Class :character   1st Qu.:-19.460   1st Qu.:-30.60  
##  Mode  :character   Mode  :character   Median :-12.469   Median : 31.97  
##                                        Mean   :-11.309   Mean   :  4.98  
##                                        3rd Qu.: -1.213   3rd Qu.: 39.12  
##                                        Max.   : 53.509   Max.   : 56.41  
##       POP8          TRACTCAS        PROPCAS            PCTOWNHOME       
##  Min.   :    9   Min.   :0.000   Min.   :0.0000000   Min.   :0.0008224  
##  1st Qu.: 2510   1st Qu.:0.310   1st Qu.:0.0000930   1st Qu.:0.4588745  
##  Median : 3433   Median :1.890   Median :0.0004130   Median :0.6508585  
##  Mean   : 3764   Mean   :2.107   Mean   :0.0005947   Mean   :0.5872621  
##  3rd Qu.: 4889   3rd Qu.:3.080   3rd Qu.:0.0009170   3rd Qu.:0.7560976  
##  Max.   :13015   Max.   :9.290   Max.   :0.0069930   Max.   :1.0000000  
##    PCTAGE65P              Z              AVGIDIST         PEXPOSURE     
##  Min.   :0.004044   Min.   :-1.9206   Min.   :0.01847   Min.   :0.6134  
##  1st Qu.:0.099926   1st Qu.:-0.7168   1st Qu.:0.02703   1st Qu.:0.9942  
##  Median :0.126415   Median :-0.2876   Median :0.03238   Median :1.1749  
##  Mean   :0.137262   Mean   :-0.2157   Mean   :0.14919   Mean   :1.8042  
##  3rd Qu.:0.160963   3rd Qu.: 0.2498   3rd Qu.:0.13008   3rd Qu.:2.5656  
##  Max.   :0.505050   Max.   : 4.7105   Max.   :3.52637   Max.   :5.8654  
##      Cases               Xm               Ym             Xshift      
##  Min.   :0.00014   Min.   :-55482   Min.   :-75291   Min.   :363839  
##  1st Qu.:0.30928   1st Qu.:-19460   1st Qu.:-30601   1st Qu.:399862  
##  Median :1.88876   Median :-12469   Median : 31970   Median :406852  
##  Mean   :2.10676   Mean   :-11309   Mean   :  4980   Mean   :408013  
##  3rd Qu.:3.08284   3rd Qu.: -1213   3rd Qu.: 39123   3rd Qu.:418108  
##  Max.   :9.28601   Max.   : 53509   Max.   : 56410   Max.   :472830  
##      Yshift       
##  Min.   :4653564  
##  1st Qu.:4698254  
##  Median :4760825  
##  Mean   :4733835  
##  3rd Qu.:4767978  
##  Max.   :4785265

Datos de ciudades.

summary(cities)

## Object of class SpatialPointsDataFrame
## Coordinates:
##                 min       max
## coords.x1  372236.8  445727.9
## coords.x2 4662140.5 4771698.3
## Is projected: TRUE 
## proj4string :
## [+proj=utm +zone=18 +ellps=WGS84 +units=m +no_defs]
## Number of points: 6
## Data attributes:
##     names          
##  Length:6          
##  Class :character  
##  Mode  :character

par(mfrow=c(1,2))
plot(NY8, border="grey60", axes=TRUE)
text(coordinates(cities), labels=as.character(cities$names), font=2, cex=0.9)
text(bbox(NY8)[1,1], bbox(NY8)[2,2], labels="a)", cex=0.8)
plot(NY8, border="grey60", axes=TRUE)
points(TCE, pch=1, cex=0.7)
points(TCE, pch=3, cex=0.7)
text(coordinates(TCE), labels=as.character(TCE$name), cex=0.7,
 font=1, pos=c(4,1,4,1,4,4,4,2,3,4,2), offset=0.3)
text(bbox(NY8)[1,1], bbox(NY8)[2,2], labels="b)", cex=0.8)

Fig. 9.1. (a) Principales ciudades de la zona de estudio de ocho condados de la parte alta del Estado de Nueva York (b) Localización de 11 vertederos de residuos peligrosos inactivos en el área de estudio.

El conjunto de datos de 281 secciones censales de ocho condados del condados del centro del Estado de Nueva York que aparecen de forma destacada en Waller y Gotway (2004) se utilizará en muchos de los ejemplos,2 complementado con los límites de los tramos derivados de TIGER 1992 y distribuidos por SEDAC/CIESIN. las principales ciudades de la zona de estudio y la ubicación de 11 vertederos de residuos peligrosos. Las cifras de Waller y Gotway (2004) incluyen masas de agua, que no están presentes en esta versión de los límites de los tractos,en esta versión de los límites de los tractos, en la que los límites de los tractos siguen los límites de los tractos siguen las líneas centrales de los lagos, en lugar de sus orillas.

9.2. Vecinos Espaciales y Pesos Espaciales.

La creación de pesos espaciales es un paso necesario en el uso de datos en área, tal vez solo para comprobar que no quedan patrones espaciales en los residuos. Dado que el análisis de datos en área depende de las opciones realizadas en la construcción de los pesos espaciales, hemos optado por mover los detalles de discusión sobre la creación de objetos vecinos incluidos en la primera edición de este libro a una viñeta incluida en el paquete spdep.

9.2.1. Objetos Vecinos.

En el paquete spdep, las relaciones vecinas entre n observaciones están representadas por un objeto de clase nb; la clase es una clase de estilo antiguo.

Si alguna observación no tiene vecinos, el componente contiene un número entero cero. También contiene atributos, típicamente un vector de identificadores de región de caracteres y un valor lógico que indica si las relaciones son simétricas. La tarjeta de función auxiliar devuelve la cardinalidad del conjunto vecino para cada objeto, es decir, el número de vecinos; difiere de la aplicación de la longitud a los componentes de la lista porque las entradas no vecinas están codificadas como un vector entero de un solo elemento con el valor de cero. En aras de la simplicidad en la demostración de cómo crear objetos vecinos, trabajamos en un subconjunto del mapa que consiste de las secciones censales dentro de Syracuse, aunque los mismos principios se aplican a el conjunto completo de datos.

Para empezar comenzaremos con la lectura de GAL en una losta de vecinos.

NY_nb <- read.gal("NY_nb.gal", region.id=row.names(NY8))

Se realiza un summy para ver que estados son vecinos.

summary(NY_nb)

## Neighbour list object:
## Number of regions: 281 
## Number of nonzero links: 1522 
## Percentage nonzero weights: 1.927534 
## Average number of links: 5.41637 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 10 11 
##  6 11 28 45 59 49 45 23 10  3  2 
## 6 least connected regions:
## 55 97 100 101 244 245 with 1 link
## 2 most connected regions:
## 34 82 with 11 links

Contigüidades de distritos censales, distritos censales de ocho condados de Nueva York

par(mfrow=c(1,1))
plot(NY8, border="grey60", axes=TRUE)
plot(NY_nb, coordinates(NY8), pch=19, cex=0.6, add=TRUE)

Proximamente creamos una serie de objetos vecinos en otras maneras. Tres son k objetos vecinos más cercanos, tomando k puntos más cercanos como vecinos, para k = 1, 2, 4 – esto se adapta a través del área de estudio, teniendo en cuenta las diferencias en las densidades de las entidades areales. Naturalmente, en la abrumadora mayoría de los casos, conduce a vecinos asimétricos, pero asegurará que todas las zonas tienen k vecinos.

Syracuse <- NY8[NY8$AREANAME == "Syracuse city",]
Sy0_nb <- subset(NY_nb, NY8$AREANAME == "Syracuse city")
summary(Sy0_nb)

## Neighbour list object:
## Number of regions: 63 
## Number of nonzero links: 346 
## Percentage nonzero weights: 8.717561 
## Average number of links: 5.492063 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 
##  1  1  5  9 14 17  9  6  1 
## 1 least connected region:
## 164 with 1 link
## 1 most connected region:
## 136 with 9 links

Los objetos k = 1 también son útiles para encontrar la distancia mínima a la que todas las regiones tienen vecinos en función de la distancia. Usando la función nbdists, podemos calcular una lista de vectores de distancia correspondientes a los objetos vecinos, aquí primero los vecinos más cercanos. El valor máximo será la distancia mínima requerida para garantizar que todas las regiones estén vinculadas al menos a un vecino, por lo que podemos crear fácilmente objetos con observaciones no adyacentes al establecer el umbral superior más bajo que eso.

9.2.2 Objetos de pesos espaciales

Los pesos espaciales se pueden ver como una lista de pesos indexados por una lista de vecinos, donde el peso del enlace entre iyj es el k-ésimo elemento del i-ésimo componente de la lista de pesos, yk nos dice qué valor en la i -ésima lista de vecinos El componente de es igual a j. j no es vecino de i si j no existe en el i-ésimo componente de la lista de vecinos. En esta sección revisaremos los métodos para ponderar objetos: Listar objetos - Construir ; esta clase es una clase de estilo antiguo . La función nb2listw toma un objeto de lista de vecinos y lo convierte en un objeto de peso . El método de impresión del objeto listw muestra las características de los vecinos subyacentes, el estilo de las ponderaciones espaciales es y las constantes de ponderación espacial utilizadas para calcular la prueba de autocorrelación espacial .

Sy0_lw_W <- nb2listw(Sy0_nb)
Sy0_lw_W

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 63 
## Number of nonzero links: 346 
## Percentage nonzero weights: 8.717561 
## Average number of links: 5.492063 
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1      S2
## W 63 3969 63 24.78291 258.564

names(Sy0_lw_W)

## [1] "style"      "neighbours" "weights"

names(attributes(Sy0_lw_W))

## [1] "names"     "class"     "region.id" "call"      "GeoDa"

El estilo de ponderación espacial(style=“W”) puede interpretarse como que permite el cálculo de promedios entre vecinos. Los enlaces que se originan en áreas con menos vecinos tienen más ponderación que los enlaces que se originan en áreas con más vecinos, lo que puede conducir a entidades de área sin darse cuenta al borde del área de estudio.

1/rev(range(card(Sy0_lw_W$neighbours)))

## [1] 0.1111111 1.0000000

summary(unlist(Sy0_lw_W$weights))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.1111  0.1429  0.1667  0.1821  0.2000  1.0000

summary(sapply(Sy0_lw_W$weights, sum))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1

Con el estilo”B” retiene un peso de unidad por cada vecino.

dsts <- nbdists(Sy0_nb, coordinates(Syracuse))
idw <- lapply(dsts, function(x) 1/(x/1000))
Sy0_lw_B <- nb2listw(Sy0_nb, style="B")
summary(unlist(Sy0_lw_B$weights))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1       1       1       1       1       1

coords <- coordinates(Syracuse)
IDs <- row.names(Syracuse)
Sy8_nb <- knn2nb(knearneigh(coords, k=1), row.names=IDs)
Sy9_nb <- knn2nb(knearneigh(coords, k=2), row.names=IDs)
Sy10_nb <- knn2nb(knearneigh(coords, k=4), row.names=IDs)
dsts <- unlist(nbdists(Sy8_nb, coords))
Sy11_nb <- dnearneigh(coords, d1=0, d2=0.75*max(dsts), row.names=IDs)

summary(sapply(Sy0_lw_B$weights, sum))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   4.500   6.000   5.492   6.500   9.000

Sy0_lw_idwB <- nb2listw(Sy0_nb, glist=idw, style="B")
summary(unlist(Sy0_lw_idwB$weights))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.3886  0.7374  0.9259  0.9963  1.1910  2.5274

summary(sapply(Sy0_lw_idwB$weights, sum))

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.304   3.986   5.869   5.471   6.737   9.435

El argumento “glist” para pasar una lista de vectores de pesos generales correspondientes a las relaciones vecinas a nb2listw.

Se propone estableces que los pesos sean proporcioanles a la inversa de la distancia entre los puntos que representan las areas , utilizando “nbdists” para calcular las distancias.**

Presentaciones de pesos espaciales para Syracuse.

library(RColorBrewer)
pal <- brewer.pal(9, "Reds")
oopar <- par(mfrow=c(1,3), mar=c(1,1,3,1)+0.1)
z <- t(listw2mat(Sy0_lw_W))
brks <- c(0,0.1,0.143,0.167,0.2,0.5,1)
nbr3 <- length(brks)-3
image(1:63, 1:63, z[,ncol(z):1], breaks=brks, col=pal[c(1,(9-nbr3):9)],
      main="W style", axes=FALSE)
box()
z <- t(listw2mat(Sy0_lw_B))
image(1:63, 1:63, z[,ncol(z):1], col=pal[c(1,9)], main="B style", axes=FALSE)
box()
z <- t(listw2mat(Sy0_lw_idwB))
brks <- c(0,0.35,0.73,0.93,1.2,2.6)
nbr3 <- length(brks)-3
image(1:63, 1:63, z[,ncol(z):1], breaks=brks, col=pal[c(1,(9-nbr3):9)],
      main="IDW B style", axes=FALSE)
box()

La siguiente figura muestra tres representaciones de los pesos espaciales de Syracuse que se muestran en forma de matriz . La imagen style=“W” de la izquierda es asimétrica , con colores más oscuros que muestran pesos más altos para las regiones con menos vecinos. Los otros dos paneles son simétricos pero expresan diferentes supuestos sobre la fuerza de los lazos vecinales.

try(Sy0_lw_D1 <- nb2listw(Sy11_nb, style="B"))

## Error in nb2listw(Sy11_nb, style = "B") : Empty neighbour sets found

Sy0_lw_D1 <- nb2listw(Sy11_nb, style="B", zero.policy=TRUE)

print(Sy0_lw_D1, zero.policy = TRUE)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 63 
## Number of nonzero links: 230 
## Percentage nonzero weights: 5.794911 
## Average number of links: 3.650794 
## 2 regions with no links:
## 154 168
## 
## Weights style: B 
## Weights constants summary:
##    n   nn  S0  S1   S2
## B 61 3721 230 460 4496

9.3 Pruebas de Autocorrelación Espacial.

Aqui se aborda la importancia de probar la presencia de autocorrelación espacial en los datos. El autor destaca que estas pruebas asumen que el modelo de media elimina cualquier patrón espacial sistemático de los datos. Sin embargo, también señala que la autocorrelación espacial aparente puede surgir debido a una mala especificación del modelo de media, lo que significa que es crucial asegurarse de que el modelo esté correctamente especificado. se menciona que una de las pruebas más comunes para la autocorrelación espacial es el índice de Moran, que se utiliza para variables continuas. También se mencionan otras pruebas globales, como el índice de Geary y la prueba de Mantel, que se utilizan para diferentes tipos de variables. El autor destaca la importancia de comprender que las pruebas de autocorrelación espacial también pueden reaccionar a un modelo de media mal especificado, lo que puede hacer que la falta de inclusión de una variable con patrón espacial en la función de media parezca autocorrelación espacial en las pruebas.Aqui se enfatiza la importancia de realizar pruebas de autocorrelación espacial para comprender mejor la estructura espacial de los datos, pero también advierte sobre la necesidad de asegurarse de que el modelo de media esté correctamente especificado para evitar interpretaciones erróneas de la autocorrelación espacial.

La prueba que se utilizará en esta discusión introductoria es el índice de Moran, que se calcula como una relación entre el producto de la variable de interés y su retraso espacial, con el producto cruzado de la variable de interés, y se ajusta en función de los pesos espaciales utilizados asi:

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

donde \(y_i\) representa la i-ésima observación, \(\bar{y}\) es la media de la variable de interés, y \(W_{ij}\) es el peso espacial del enlace entre \(i\) y \(j\).

set.seed(987654)
n <- length(Sy0_nb)
uncorr_x <- rnorm(n)
rho <- 0.5
autocorr_x <- spatialreg::invIrW(Sy0_lw_W, rho) %*% uncorr_x

library(rgdal)
library(spdep)
moran.test(uncorr_x, listw=Sy0_lw_W)

## 
##  Moran I test under randomisation
## 
## data:  uncorr_x  
## weights: Sy0_lw_W    
## 
## Moran I statistic standard deviate = -0.22698, p-value = 0.5898
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.033287124      -0.016129032       0.005714128

La salida previa corresponde al escenario en el que no hay correlación. Evidentemente, el valor p de 0.5898 obtenido respalda la hipótesis de que no existe dependencia espacial con los pesos proporcionados por el objeto Sy0_lw_W. Esto contradice la hipótesis de que sí existe dependencia, recordando que los pesos mencionados se obtienen al aplicar la función nb2listw a los pesos originales asignados a los enlaces de la ciudad de Syracuse.

moran.test(autocorr_x, listw=Sy0_lw_W)

## 
##  Moran I test under randomisation
## 
## data:  autocorr_x  
## weights: Sy0_lw_W    
## 
## Moran I statistic standard deviate = 3.1027, p-value = 0.0009588
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.218178182      -0.016129032       0.005702793

Con base en la salida anterior, donde se utilizan las ponderaciones espaciales para una variable que muestra autocorrelación espacial, como era de esperarse, se obtiene un resultado significativo.

moran.test(autocorr_x, listw=nb2listw(Sy9_nb, style="W"))

## 
##  Moran I test under randomisation
## 
## data:  autocorr_x  
## weights: nb2listw(Sy9_nb, style = "W")    
## 
## Moran I statistic standard deviate = 1.8619, p-value = 0.03131
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.19207401       -0.01612903        0.01250486

En este caso, al utilizar pesos que garantizan que todas las zonas tengan un vecino y una lista de pesos codificados en el estilo W, donde se realiza una normalización por filas (sumando sobre todos los enlaces hasta n), se encuentra evidencia de autocorrelación espacial.

Los siguientes códigos y sus resultados tienen como objetivo demostrar lo que puede suceder cuando la suposición de una media constante es incorrecta, resaltando la importancia de especificar modelos adecuados al probar la autocorrelación espacial en la construcción de modelos.

et <- coords[,1] - min(coords[,1])
trend_x <- uncorr_x + 0.00025 * et
moran.test(trend_x, listw=Sy0_lw_W)

## 
##  Moran I test under randomisation
## 
## data:  trend_x  
## weights: Sy0_lw_W    
## 
## Moran I statistic standard deviate = 3.3438, p-value = 0.0004131
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.237474257      -0.016129032       0.005751989

lm.morantest(lm(trend_x ~ et), listw=Sy0_lw_W)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = trend_x ~ et)
## weights: Sy0_lw_W
## 
## Moran I statistic standard deviate = -0.31166, p-value = 0.6223
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     -0.053830651     -0.030937453      0.005395814

Para calcular estos índices, es necesario establecer la proximidad entre las áreas. Esto se logra mediante el cálculo de una matriz de proximidad espacial. Dado un conjunto de n áreas (\(A_1\),…,\(A_n\)), se crea una matriz \(W\) de tamaño \(n×n\), donde cada elemento \(W_{ij}\) representa una medida de proximidad entre \(A_i\) y \(A_j\).

Si \(Y(s)\) es el atributo \(Y\) observado en el sitio \(s\), entonces el término autocorrelación espacial se refiere a la correlación entre \(Y_{si}\) y \(Z_{sj}\). Se trata de una correlación entre el mismo atributo en dos sitios diferentes. Si existe una autocorrelación espacial positiva, la cercanía en el espacio debe estar relacionada con la similitud en los valores de los atributos.

9.3.1 Test Globales

Aplicado en R con la función moran.test el I de Moran es el test más común para probar autocorrelación espacial, aunque no es el único test implementado en el paquete spdep para variables continuas:

C de Geary (geary.test()) **Global GetisOrd (globalG.test())* Prueba de producto cruzado espacial general de Mantel, que incluye la I de Moran, la C de Geary y la variante Sokal de la C de Geary (sp.mantel.mc()) Cuando la variable de interés es categórica, representada como factor, pueden implementarse las funciones joincount.test() para uniones del mismo color, joincount.multi() para uniones del mismo y de distinto color.

Nuevamente, se debe tener en cuenta que los resultados pueden depender de las elecciones del tipo de peso espacial. Para ello, las pruebas de Monte Carlo o las basadas en permutaciones bootstrap, en las que los valores de la variable de interés se asignan aleatoriamente a entidades espaciales, pueden ser de ayuda contra los errores de inferencia. De hecho, dado que las pruebas de autocorrelación espacial son sensibles a los patrones espaciales en la variable de interés, los métodos mencionados no son necesariamente buenas guías para decidir lo que ocurre en el proceso de generación de datos. A veces se necesita un bootstrapping paramétrico o pruebas específicamente adaptadas al entorno, o una mejor especificación de la variable de interés.

Cuando se desarrollaron las medidas de autocorrelación se supone que todas las entidades tenían vecinos, por lo que no es obvio lo que hay que hacer cuando algunas no los tienen. Algo similar a lo que pasaba con la elección de la anchura de los variogramas en geoestadística.

Inicialmente, el interés radica en busca de autocorrelación, probando el número de casos de leucemia por sección censal utilizando el estilo predeterminado de pesos espaciales de estandarización de filas, y utilizando el supuesto de aleatoriedad analítica en el cálculo de la varianza.

La noción de autocorrelación espacial de estas variables está asociada con la idea de que valores observados en áreas geográficas adyacentes serán más similares que los esperados bajo el supuesto de independencia espacial.

Indice de Moran Para probar autocorrelación espacial, el test más usado es el I de Moran, que se calcula como un cociente del producto de la variable de interés y su “desfase” espacial, con el producto cruzado de la variable de interés, y ajustado por las ponderaciones espaciales utilizadas, considerando la información de los vecinos más cercanos, el índice de autocorrelación de Moran es definido como:

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

donde \(y_i\) representa la i-ésima observación, \(\bar{y}\) es la media de la variable de interés, y \(W_{ij}\) es el peso espacial del enlace entre \(i\) y \(j\).

moran.test(NY8$Cases, listw=nb2listw(NY_nb))

## 
##  Moran I test under randomisation
## 
## data:  NY8$Cases  
## weights: nb2listw(NY_nb)    
## 
## Moran I statistic standard deviate = 3.9778, p-value = 3.477e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.146882990      -0.003571429       0.001430595

lw_B <- nb2listw(NY_nb, style="B")
moran.test(NY8$Cases, listw=lw_B)

## 
##  Moran I test under randomisation
## 
## data:  NY8$Cases  
## weights: lw_B    
## 
## Moran I statistic standard deviate = 3.1862, p-value = 0.0007207
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.110387402      -0.003571429       0.001279217

moran.test(NY8$Cases, listw=lw_B, randomisation=FALSE)

## 
##  Moran I test under normality
## 
## data:  NY8$Cases  
## weights: lw_B    
## 
## Moran I statistic standard deviate = 3.1825, p-value = 0.0007301
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.110387402      -0.003571429       0.001282228

lm.morantest(lm(Cases ~ 1, NY8), listw=lw_B)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = Cases ~ 1, data = NY8)
## weights: lw_B
## 
## Moran I statistic standard deviate = 3.1825, p-value = 0.0007301
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.110387402     -0.003571429      0.001282228

lm.morantest.sad(lm(Cases ~ 1, NY8), listw=lw_B)

## 
##  Saddlepoint approximation for global Moran's I (Barndorff-Nielsen
##  formula)
## 
## data:  
## model:lm(formula = Cases ~ 1, data = NY8)
## weights: lw_B
## 
## Saddlepoint approximation = 2.9929, p-value = 0.001382
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I 
##        0.1103874

lm.morantest.exact(lm(Cases ~ 1, NY8), listw=lw_B)

## 
##  Global Moran I statistic with exact p-value
## 
## data:  
## model:lm(formula = Cases ~ 1, data = NY8)
## weights: lw_B
## 
## Exact standard deviate = 2.9923, p-value = 0.001384
## alternative hypothesis: greater
## sample estimates:
## [1] 0.1103874

En los anteriores comandos podemos comprobar que una explicación detallada sobre los diferentes enfoques y supuestos utilizados en el análisis de la asociación espacial. Se discuten los efectos de cambiar los pesos espaciales en los resultados y se prueba los supuestos de normalidad y aleatorización utilizados en las pruebas de Moran. Además, se destaca la conexión entre la prueba de Moran y la regresión de residuos, lo que permite la inclusión de variables adicionales en el modelo. También se mencionan métodos alternativos, como la aproximación de Saddlepoint y la prueba exacta, aunque se señala que su aplicación computacional puede ser más demandante. En general, la prueba proporciona una visión amplia y reflexiva sobre las diferentes consideraciones y opciones disponibles en el análisis de asociación espacial.

utilizando simulacion montecarlo

set.seed(1234)
bperm <- moran.mc(NY8$Cases, listw=lw_B, nsim=999)
bperm

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  NY8$Cases 
## weights: lw_B  
## number of simulations + 1: 1000 
## 
## statistic = 0.11039, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greater

También podemos utilizar una prueba de Monte Carlo, una prueba de bootstrap por permutación, en la cual los valores observados se asignan aleatoriamente a regiones y se calcula la estadística de interés nsim veces. Dado que, en el caso global en contraste con el local, tenemos suficientes observaciones y podemos repetir esta permutación potencialmente muchas veces sin repetición.

r <- sum(NY8$Cases)/sum(NY8$POP8)
rni <- r*NY8$POP8
CR <- function(var, mle) rpois(length(var), lambda=mle)
MoranI.pboot <- function(var, i, listw, n, S0, ...) {
  return(moran(x=var, listw=listw, n=n, S0=S0)$I)
}
set.seed(1234)

library(boot)
boot2 <- boot(NY8$Cases, statistic=MoranI.pboot, R=999, sim="parametric",
  ran.gen=CR, listw=lw_B, n=length(NY8$Cases), S0=Szero(lw_B), mle=rni)

pnorm((boot2$t0 - mean(boot2$t))/sd(boot2$t[,1]), lower.tail=FALSE)

## [1] 0.1472112

Se emplean recuentos esperados RNI como parámetros lambda en la función Rpois para generar conjuntos de datos sintéticos basados en la distribución de Poisson. La probabilidad de salida se calcula restando la media de los valores simulados de Moran’s I del valor observado, y se divide por su desviación estándar. En la Figura 9.5, se muestra la representación de Waller y Gotway (2004, p. 232, Fig. 7.8), con simulaciones paramétricas que desplazan la distribución de Moran’s I hacia la derecha para tener en cuenta la influencia de las poblaciones heterogéneas en las áreas.

Asimismo, existe una versión modificada de INdice de MOran que utiliza una tasa de Bayes Empírico. A diferencia de los resultados previos, esta versión reduce las tasas extremas en áreas con poblaciones pequeñas en riesgo y las ajusta hacia la tasa general del área. Además, se emplean métodos de Monte Carlo para realizar inferencias.

oopar <- par(mfrow=c(1,2))
xlim <- range(c(bperm$res, boot2$t[,1]))
hist(bperm$res[-length(bperm$res)], main="Permutation bootstrap", xlab=expression(I[std]), xlim=xlim, density=15, angle=45, ylim=c(0,260))
abline(v=bperm$statistic, lty=2)
hist(boot2$t, col=rgb(0.4,0.4,0.4), main="Parametric bootstrap", xlab=expression(I[CR]), xlim=xlim, ylim=c(0,260))
hist(bperm$res[-length(bperm$res)], density=15, angle=45, add=TRUE)
abline(v=boot2$t0, lty=2)

par(oopar)

En la figura se ve que los histogramas de valores simulados de indice de Moran bajo permutaciones aleatorias de los datos y muestras paramétricas de valores esperados de riesgo constante; el valor observado de indice de Moran está marcado por una línea vertical. La presencia de valores inusuales y extremos en áreas con poblaciones pequeñas podría explicar la falta de significancia en los resultados de la bootstrap paramétrica para las tasas de riesgo observadas y esperadas. Después de suavizar las tasas, se observa cierta autocorrelación global.

set.seed(1234)
EBImoran.mc(n=NY8$Cases, x=NY8$POP8, listw=nb2listw(NY_nb, style="B"), nsim=999)

## The default for subtract_mean_in_numerator set TRUE from February 2016

## 
##  Monte-Carlo simulation of Empirical Bayes Index (mean subtracted)
## 
## data:  cases: NY8$Cases, risk population: NY8$POP8
## weights: nb2listw(NY_nb, style = "B")
## number of simulations + 1: 1000
## 
## statistic = 0.072497, observed rank = 983, p-value = 0.017
## alternative hypothesis: greater

cor8 <- sp.correlogram(neighbours=NY_nb, var=NY8$Cases, order=8, method="I", style="C")

library(pgirmess)

## Warning: package 'pgirmess' was built under R version 4.2.3

corD <- correlog(coordinates(NY8), NY8$Cases, method="Moran")

oopar <- par(mfrow=c(1,2))
plot(cor8, main="Contiguity lag orders")
plot(corD, main="Distance bands")

par(oopar)

En la figura muestra las Correlogramas: (izquierda) valores de Indice Moran para ocho órdenes de rezago sucesivos de vecinos contiguos; (derecha) valores de Indice de Moran para una secuencia de pares de vecinos en bandas de distancia.

Aqui La función selecciona automáticamente bandas de distancia de aproximadamente 10 km, abarcando toda el área de estudio. En este caso, las dos primeras bandas de 0-10 y 10-20 km.

Test Locales.

Las pruebas globales de autocorrelación espacial se calculan a partir de las relaciones locales entre los valores observados en una entidad espacial y sus vecinos, según la definición de vecindad elegida.

oopar <- par(mfrow=c(1,2))
msp <- moran.plot(NY8$Cases, listw=nb2listw(NY_nb, style="C"), quiet=TRUE)
title("Moran scatterplot")
infl <- ifelse(packageVersion("spdep") > "1.1.4", msp$is_inf, apply(msp$is.inf, 1, any))
x <- NY8$Cases
lhx <- cut(x, breaks=c(min(x), mean(x), max(x)), labels=c("L", "H"), include.lowest=TRUE)
wx <- lag(nb2listw(NY_nb, style="C"), NY8$Cases)
lhwx <- cut(wx, breaks=c(min(wx), mean(wx), max(wx)), labels=c("L", "H"), include.lowest=TRUE)
lhlh <- interaction(lhx, lhwx, infl, drop=TRUE)
cols <- rep(1, length(lhlh))
cols[lhlh == "H.L.TRUE"] <- 2
cols[lhlh == "L.H.TRUE"] <- 3
cols[lhlh == "H.H.TRUE"] <- 4
plot(NY8, col=brewer.pal(4, "Accent")[cols])
legend("topright", legend=c("None", "HL", "LH", "HH"), fill=brewer.pal(4, "Accent"), bty="n", cex=0.8, y.intersp=0.8)
title("Tracts with influence")

par(oopar)

Aqui, vemos un diagrama de dispersión de Moran de la variable de recuento de casos de leucemia; para una discusión sobre los diagramas de dispersión de Moran. El gráfico (mostrado en la Figura) coloca convencionalmente la variable de interés en el eje x y la suma ponderada espacialmente de los valores de los vecinos (los valores retrasados espacialmente) en el eje y. Moran’s I global es una relación lineal entre ellos y se representa como una pendiente. Dado que el Indice de Moran global es, al igual que los coeficientes de correlación similares, una relación lineal, también podemos aplicar técnicas estándar para detectar observaciones con una influencia inusualmente fuerte en la pendiente.

lm1 <- localmoran(NY8$Cases, listw=nb2listw(NY_nb, style="C"))
lm2 <- as.data.frame(localmoran.sad(lm(Cases ~ 1, NY8), nb=NY_nb, style="C"))
lm3 <- as.data.frame(localmoran.exact(lm(Cases ~ 1, NY8), nb=NY_nb, style="C"))

r <- sum(NY8$Cases)/sum(NY8$POP8)
rni <- r*NY8$POP8
lw <- nb2listw(NY_nb, style="C")
sdCR <- (NY8$Cases - rni)/sqrt(rni)
wsdCR <- lag(lw, sdCR)
I_CR <- sdCR * wsdCR

library(RColorBrewer)
gry <- c(rev(brewer.pal(8, "Reds")[1:6]), brewer.pal(6, "Blues"))
NY8$Standard <- lm1[,1]
NY8$"Constant_risk" <- I_CR
#nms <- match(c("Standard", "Constant_risk"), names(NY8))
spplot(NY8, c("Standard", "Constant_risk"), at=c(-2.5,-1.4,-0.6,-0.2,0,0.2,0.6,4,7), col.regions=colorRampPalette(gry)(8))

La Figura muestra los dos conjuntos de valores del Moran’s Ii local, calculados de manera estándar y utilizando la suposición de Poisson para la hipótesis de riesgo constante. Ya sabemos que el Moran’s Ii global puede variar en valor y en inferencia dependiendo de nuestras suposiciones, por ejemplo, que la inferencia debe tener en cuenta las desviaciones de nuestras suposiciones de distribución.

set.seed(1234)
nsim <- 999
N <- length(rni)
sims <- matrix(0, ncol=nsim, nrow=N)
for (i in 1:nsim) {
  y <- rpois(N, lambda=rni)
  sdCRi <- (y - rni)/sqrt(rni)
  wsdCRi <- lag(lw, sdCRi)
  sims[,i] <- sdCRi * wsdCRi 
}
xrank <- apply(cbind(I_CR, sims), 1, function(x) rank(x)[1])
diff <- nsim - xrank
diff <- ifelse(diff > 0, diff, 0)
pval <- punif((diff + 1)/(nsim + 1))

NY8$Normal <- lm2[,3]
NY8$Randomisation <- lm1[,5]
NY8$Saddlepoint <- lm2[,5]
NY8$Exact <- lm3[,5]
NY8$Constant_risk <- pval
gry <- c(rev(brewer.pal(6, "Reds")), brewer.pal(6, "Blues"))
spplot(NY8, c("Normal", "Randomisation", "Saddlepoint", "Exact", "Constant_risk"), at=c(0,0.01,0.05,0.1,0.9,0.95,0.99,1), col.regions=colorRampPalette(gry)(7))

La figura muestra losValores de probabilidad para todos los distritos censales, Valores de probabilidad para todos los distritos censales, INdice de Moran locales: suposiciones de normalidad y randomización, aproximación de Saddlepoint, valores exactos y hipótesis de riesgo constante INdice de Moran locales: suposiciones de normalidad y randomización, aproximación de Saddlepoint, valores exactos y hipótesis de riesgo constante.

spplot(NY8, c("Normal", "Exact", "Constant_risk"), xlim=c(405200, 432200), ylim=c(4652700, 4672000), at=c(0,0.01,0.05,0.1,0.9,0.95,0.99,1), col.regions=colorRampPalette(gry)(7))

En la figura Parece que el uso del enfoque de riesgo constante maneja mejor la heterogeneidad en los recuentos que las alternativas, pero debemos tener en cuenta que nuestro modelo subyacente es muy simplista. ahora vamos a probar los entrenamientos de algoritmos para modelar datos de areas.

9.4 Ajuste de modelos de datos de Area

Como hemos visto anteriormente, la falta de independencia entre las observaciones en datos espaciales, es decir, la autocorrelación espacial, es común y existen pruebas disponibles para detectarla. En un mundo ideal, preferiríamos recopilar datos en los que las observaciones fueran mutuamente independientes, para evitar problemas en la inferencia a partir de los resultados analíticos. ##9.4.1 Enfoques de la Estadistica Espacial Existen diversas formas de abordar la dependencia espacial mediante el uso de modelos estadísticos. En muchos casos, es habitual asumir que las observaciones son independientes y tienen la misma distribución, pero esto puede no ser cierto cuando se trabaja con datos espaciales. Las observaciones no son independientes debido a la posible correlación entre áreas adyacentes. Además, puede resultar difícil distinguir entre el impacto de la autocorrelación espacial y las disparidades espaciales en la distribución de las observaciones.

Desde una perspectiva estadística, es posible tomar en cuenta la correlación entre las observaciones al considerar una estructura específica en el modelo. Si el vector de variables de respuesta sigue una distribución normal multivariada, podemos expresar el modelo de la siguiente manera: \[Y=μ+e\] En el contexto del modelo, el vector de medias de área \(μ\) puede ser modelado de varias formas, y el vector de errores aleatorios e se asume que sigue una distribución normal con media cero y varianza genérica V. En muchos casos, se supone que la media depende linealmente de ciertas covariables X, por lo que reemplazamos la media por \(Xtβ\) en el modelo. Por otro lado, consideramos la correlación entre áreas utilizando una forma específica de la matriz de varianza V. Si las variables no siguen una distribución normal, podríamos transformar los datos originales para lograr la normalidad deseada. Por lo tanto, las técnicas descritas a continuación pueden aplicarse también a los datos transformados. Existen varias estructuras de correlación posibles para tener en cuenta la correlación espacial, pero nos centraremos en dos enfoques comúnmente utilizados en la práctica: los Modelos SAR (Autoregresivos Simultáneos) y los Modelos CAR (Condicionalmente Autoregresivos).

nylm <- lm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8)
summary(nylm)

## 
## Call:
## lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -1.7417 -0.3957 -0.0326  0.3353  4.1398 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.51728    0.15856  -3.262  0.00124 ** 
## PEXPOSURE    0.04884    0.03506   1.393  0.16480    
## PCTAGE65P    3.95089    0.60550   6.525 3.22e-10 ***
## PCTOWNHOME  -0.56004    0.17031  -3.288  0.00114 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6571 on 277 degrees of freedom
## Multiple R-squared:  0.1932, Adjusted R-squared:  0.1844 
## F-statistic:  22.1 on 3 and 277 DF,  p-value: 7.306e-13

NY8$lmresid <- residuals(nylm)

Se muestra la distribución espacial de los valores residuales para el estudio de tramos censales del área. Las dos variables censales parecen contribuir a explicar la varianza en la respuesta. Además, aunque hay menos autocorrelación espacial en los residuos del modelo con covariables que en el modelo nulo, es claro que hay información en los residuos que deberíamos tratar de usar. Una prueba exacta para la autocorrelación espacial en los residuos lleva a conclusiones similares. Dado que la prueba de Moran está destinada a detectar la autocorrelación espacial, podemos trate de ajustar un modelo teniendo esto en cuenta. Sin embargo, no debemos olvidar que las especificaciones erróneas detectadas por el I de Moran pueden tener una variedad de causas . También se da el caso de que si el modelo ajustado exhibe multicolinealidad, los resultados de la prueba pueden verse afectados debido a las consecuencias numéricas de la matriz del modelo que no es de rango completo para la expectativa y la varianza del estadistico.

NYlistw<-nb2listw(NY_nb, style = "B")
lm.morantest(nylm, NYlistw)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8)
## weights: NYlistw
## 
## Moran I statistic standard deviate = 2.638, p-value = 0.004169
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.083090278     -0.009891282      0.001242320

Se muestra el calculo del I de Moran de la regresion global de los datos de NY8, donde podemos ver que el estadistico es significativo a una significancia del 0.05 y el I de Moran observado es I=0.083090278

NYlistwW <- nb2listw(NY_nb, style = "W")
library(spatialreg)
aple(residuals(nylm), listw=NYlistwW)

## [1] 0.2051536

spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, listw=NYlistwW)$lambda

##    lambda 
## 0.2169261

9.4.1.1 Modelos autorregresivos simultáneos

La especificación SAR utiliza una regresión sobre los valores de las otras áreas para dar cuenta de la dependencia espacial. Esto significa que los términos de error ei son modelado para que dependan unos de otros de la siguiente manera:

\[e_i = \sum_{i=1}^{m} b_{ij}e_i + \xi_i, \] Aquí, \(\xi_i\) se utiliza para representar errores residuales, que se supone que se distribuyen de forma independiente de acuerdo con una distribución normal con media cero y matriz de covarianza diagonal \(\Sigma_{\xi}\) con elementos \(\sigma^2_{\xi_i}\), i = 1,…,m. Los valores \(b_{ij}\) se utilizan para representar dependencia espacial entre áreas. \(b_{ii}\) debe establecerse en cero para que cada área no retroceda sobre sí misma. Tenga en cuenta que si expresamos los términos de error como: \[Y = \beta X^T + B(Y - X^T \beta) + e, \] el modelo también se puede expresar como \[(I-B)(Y-X^T \beta)=e,\]

donde B es una matriz que contiene los parámetros de dependencia \(b_{ij}\) e I es la matriz identidad de la dimensión requerida. Es importante señalar que para que este modelo SAR esté bien definido, la matriz I − B debe ser no singular. Bajo este modelo, \(Y\) se distribuye de acuerdo a una normal multivariante con media \(E[Y] = X^T \beta\) y una matriz de covarianzas: \[Var[Y]= (I-B)^{-1} \Sigma_{\xi}(I-B^T)^{-1}\] Usualmente, \(\Sigma_\xi\) depende de un solo parametro \(\sigma^2\), tal que \(\Sigma_\xi = \sigma^2 I\) y \(Var[Y]\) se simplifica a: \[Var[Y] = \sigma^2(I-B)^{-1}(I-B^T)^{-1}\] Tambien es posible especificar a \(\Sigma_\xi\) como una matriz diagonal de pesos asociada con la hetérogeneidad entre las observaciones. Se puede obtener una reparametrización útil de este modelo escribiendo \(B = \lambda W\), donde \(\lambda\) es un parámetro de autocorrelación espacial y W es una matriz que representa la dependencia espacial, con esta especificación la varianza de Y se convierte en: \[Var[Y] = \sigma^2(I-\lambda W)^{-1}(I-\lambda W^T)^{-1}\]

Sin embargo, este modelo no tiene en cuenta la distribución heterogénea de la población por tramos más allá de la corrección introducida en la transformación proporciones de incidencia. La versión ponderada de estos modelos se puede instalar para que los distritos se ponderan proporcionalmente a la inversa del tamaño de su población. Para ello incluimos el parámetro weights=POP8 en la llamada de la funcion.

nysar<-spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, listw=NYlistw)
summary(nysar)

## 
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8, 
##     listw = NYlistw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.56754 -0.38239 -0.02643  0.33109  4.01219 
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -0.618193   0.176784 -3.4969 0.0004707
## PEXPOSURE    0.071014   0.042051  1.6888 0.0912635
## PCTAGE65P    3.754200   0.624722  6.0094 1.862e-09
## PCTOWNHOME  -0.419890   0.191329 -2.1946 0.0281930
## 
## Lambda: 0.040487 LR test value: 5.2438 p-value: 0.022026 
## Numerical Hessian standard error of lambda: 0.017199 
## 
## Log likelihood: -276.1069 
## ML residual variance (sigma squared): 0.41388, (sigma: 0.64333)
## Number of observations: 281 
## Number of parameters estimated: 6 
## AIC: NA (not available for weighted model)

Al analizar el modelo lineal ponderado, observamos que la variable de exposición al TCE ha adquirido relevancia con la orientación esperada, lo que sugiere que las extensiones más próximas a los lugares con TCE presentan tasas de incidencia transformadas ligeramente más elevadas. Además, las otras dos variables independientes también presentan coeficientes más significativos en este momento.

nylam1 <- c(nysar$lambda)
nylam2 <- c(LR1.Spautolm(nysar)$p.value)

NY8$sar_trend <- nysar$fit$signal_trend
NY8$sar_stochastic <- nysar$fit$signal_stochastic
rds <- colorRampPalette(brewer.pal(8, "RdBu"))
tr_at <- seq(-1, 1.3, length.out=21)
tr_rds <- rds(sum(tr_at >= 0)*2)[-(1:(sum(tr_at >= 0)-sum(tr_at < 0)))]
tr_pl <- spplot(NY8, c("sar_trend"), at=tr_at, col="transparent", col.regions=tr_rds, main=list(label="Trend", cex=0.8))
st_at <- seq(-0.16, 0.39, length.out=21)
st_rds <- rds(sum(st_at >= 0)*2)[-(1:(sum(st_at >= 0)-sum(st_at < 0)))]
st_pl <- spplot(NY8, c("sar_stochastic"), at=st_at, col="transparent", col.regions=st_rds, main=list(label="Stochastic", cex=0.8))
plot(tr_pl, split=c(1,1,2,1), more=TRUE)
plot(st_pl, split=c(2,1,2,1), more=FALSE)

En la figura muestra Componente de tendencia de los valores ajustados del modelo SAR; (derecha): Componente estocástica espacial de los valores ajustados del modelo SAR. Este modelo no tiene en cuenta la distribucion heterogénea de la poblacion por distritos, por eso incluimos el parametro weights=POP8.

nylmw <- lm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, weights=POP8)
summary(nylmw)

## 
## Call:
## lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8, 
##     weights = POP8)
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -129.067  -14.714    5.817   25.624   70.723 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.77837    0.14116  -5.514 8.03e-08 ***
## PEXPOSURE    0.07626    0.02731   2.792  0.00560 ** 
## PCTAGE65P    3.85656    0.57126   6.751 8.60e-11 ***
## PCTOWNHOME  -0.39869    0.15305  -2.605  0.00968 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 33.5 on 277 degrees of freedom
## Multiple R-squared:  0.1977, Adjusted R-squared:  0.189 
## F-statistic: 22.75 on 3 and 277 DF,  p-value: 3.382e-13

NY8$lmwresid <- residuals(nylmw)

library(RColorBrewer)
gry <- c(rev(brewer.pal(6, "Reds")[1:4]), colorRampPalette(brewer.pal(5, "Blues"))(9))
TCEpts <- list("sp.points", TCE, pch=16, col="grey5")
spplot(NY8, c("lmresid", "lmwresid"), sp.layout=list(TCEpts), col.regions=gry, col="transparent", lwd=0.5, at=seq(-2,4.5,0.5))

En la figura (Izquierda): Residuos del modelo lineal de proporciones de incidencia transformadas; (derecha): Residuos del modelo lineal ponderado de proporciones de incidencia transformadas; se muestran las ubicaciones de los sitios TCE con fines comparativos.

lm.morantest(nylmw, NYlistw)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8,
## weights = POP8)
## weights: NYlistw
## 
## Moran I statistic standard deviate = 0.4773, p-value = 0.3166
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.007533246     -0.009309771      0.001245248

Las pruebas de Moran para residuos de regresión también se pueden utilizar con un objeto de modelo lineal. Los resultados son interesantes y sugieren que la especificación errónea detectada por el I de Moran está relacionada con la heteroscedasticidad. más que a la autocorrelación espacial. Podemos verificar esto para el modelo SAR también, ya que spautolm también toma un argumento de pesos.

nysarw <- spatialreg::spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, listw=NYlistw, weights=POP8)
summary(nysarw)

## 
## Call: spatialreg::spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, 
##     data = NY8, listw = NYlistw, weights = POP8)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.48488 -0.26823  0.09489  0.46552  4.28343 
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -0.797063   0.144054 -5.5331 3.146e-08
## PEXPOSURE    0.080545   0.028334  2.8428  0.004473
## PCTAGE65P    3.816731   0.576037  6.6258 3.453e-11
## PCTOWNHOME  -0.380778   0.156507 -2.4330  0.014975
## 
## Lambda: 0.0095636 LR test value: 0.32665 p-value: 0.56764 
## Numerical Hessian standard error of lambda: 0.016466 
## 
## Log likelihood: -251.6017 
## ML residual variance (sigma squared): 1104.1, (sigma: 33.229)
## Number of observations: 281 
## Number of parameters estimated: 6 
## AIC: NA (not available for weighted model)

NY8$sarw_trend <- nysarw$fit$signal_trend
NY8$sarw_stochastic <- nysarw$fit$signal_stochastic
tr_pl <- spplot(NY8, c("sarw_trend"), at=tr_at, col="transparent", col.regions=tr_rds, main=list(label="Trend", cex=0.8))
st_pl <- spplot(NY8, c("sarw_stochastic"), at=st_at, col="transparent", col.regions=st_rds, main=list(label="Stochastic", cex=0.8))
plot(tr_pl, split=c(1,1,2,1), more=TRUE)
plot(st_pl, split=c(2,1,2,1), more=FALSE)

En la figura Izquierda: Componente de tendencia de los valores ajustados del modelo SAR ponderado. Derecha: Componente espacial estocástica de los valores ajustados del modelo SAR ponderado. Para comparar ambos modelos y elegir el mejor, utilizamos el Criterio de Información de Akaike (AIC) reportado en los resúmenes del modelo. El AIC es una suma ponderada deEl logaritmo de verosimilitud del modelo y el número de coeficientes ajustados; según el criterio, los mejores modelos son aquellos con valores más bajos de AIC. Por lo tanto, el modelo ponderado proporciona un mejor ajuste ya que su AIC es considerablemente más bajo. Esto indica la importancia de tener en cuenta las poblaciones heterogéneas en el análisis de este tipo de datos de malla.

9.4.1.2 Modelos Condicionales Autoregresivos

La especificación CAR se basa en la distribución condicional de la distribución espacial en términos de error. En este caso, la distribución de \(e_i\) condicionada a \(e_{−i}\) (el vector de todos los términos de error aleatorio menos ei mismo). En lugar de todo el e−i vector, sólo se utilizan los vecinos del área i, definidos de una manera elegida. Los representamos por \(e_{j∼i}\). Entonces, una forma simple de poner el condicional distribución de \(e_i\) es:

\[e_i |e_{j∼i} \sim N ( \sum_ {j∼i} \frac{c_{ij}e_j}{\sum_ {j∼i} c_{}ij},\frac{\sigma^2_{}e_i}{\sum_ {j∼i} c_{}ij})\]

donde \(c_{ij}\) son parámetros de dependencia similares a \(b_{ij}\). Sin embargo, especificar el distribuciones condicionales de los términos de error no implica que la combinación existe. Para tener una distribución adecuada, se deben establecer algunas restricciones en los parámetros del modelo. Para nuestros propósitos de modelado, las pautas anteriores serán suficientes para obtener una especificación CAR adecuada en la mayoría de los casos. Para adaptarse a un modelo de COCHE, podemos usar la función spautolm nuevamente. esta vez nosotros establezca el argumento family=“CAR” para especificar que estamos ajustando este tipo de modelos

nycar<-spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME , data=NY8, family="CAR",
   listw=NYlistw)
summary(nycar)

## 
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8, 
##     listw = NYlistw, family = "CAR")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.539732 -0.384311 -0.030646  0.335126  3.808848 
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -0.648362   0.181129 -3.5796 0.0003442
## PEXPOSURE    0.077899   0.043692  1.7829 0.0745986
## PCTAGE65P    3.703830   0.627185  5.9055 3.516e-09
## PCTOWNHOME  -0.382789   0.195564 -1.9574 0.0503053
## 
## Lambda: 0.084123 LR test value: 5.8009 p-value: 0.016018 
## Numerical Hessian standard error of lambda: 0.030868 
## 
## Log likelihood: -275.8283 
## ML residual variance (sigma squared): 0.40758, (sigma: 0.63842)
## Number of observations: 281 
## Number of parameters estimated: 6 
## AIC: NA (not available for weighted model)

nylam1 <- c(nycar$lambda)

Los coeficientes estimados de las covariables en el modelo son muy similares a los obtenidos con los modelos SAR. Sin embargo, los valores de p de dos covariables, la distancia al TCE más cercano y el porcentaje de personas ser propietario de una vivienda, están ligeramente por encima del umbral de 0,05. La prueba de la razón de verosimilitud indica que existe una autocorrelación espacial significativa y la estimación el valor de \(λ\) es 0.0841. Considerando una regresión ponderada, usando el tamaño de la población como pesos, para que el mismo modelo dé cuenta de la distribución heterogénea de la la población elimina por completo la autocorrelación espacial en los datos. Los coeficientes de las covariables no cambian mucho y todos ellos se vuelven importante. Por lo tanto, modelar la autocorrelación espacial mediante especificaciones SAR o CAR no cambia los resultados obtenidos.

nycarw <- spatialreg::spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, family="CAR",
   listw=NYlistw, weights=POP8)
summary(nycarw)

## 
## Call: spatialreg::spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, 
##     data = NY8, listw = NYlistw, weights = POP8, family = "CAR")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.491042 -0.270906  0.081435  0.451556  4.198134 
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -0.790154   0.144862 -5.4545 4.910e-08
## PEXPOSURE    0.081922   0.028593  2.8651  0.004169
## PCTAGE65P    3.825858   0.577720  6.6223 3.536e-11
## PCTOWNHOME  -0.386820   0.157436 -2.4570  0.014010
## 
## Lambda: 0.022419 LR test value: 0.38785 p-value: 0.53343 
## Numerical Hessian standard error of lambda: 0.038543 
## 
## Log likelihood: -251.5711 
## ML residual variance (sigma squared): 1102.9, (sigma: 33.21)
## Number of observations: 281 
## Number of parameters estimated: 6 
## AIC: NA (not available for weighted model)

nysarwM<-spatialreg::spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, family="SAR",
   listw=NYlistw, weights=POP8, method="Matrix")

summary(nysarwM)

## 
## Call: spatialreg::spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, 
##     data = NY8, listw = NYlistw, weights = POP8, family = "SAR", 
##     method = "Matrix")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.561100 -0.374524  0.057405  0.591094  5.700142 
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -0.578546   0.090006 -6.4279 1.294e-10
## PEXPOSURE    0.035402   0.013959  2.5361   0.01121
## PCTAGE65P    4.651137   0.421285 11.0404 < 2.2e-16
## PCTOWNHOME  -0.666898   0.091443 -7.2931 3.029e-13
## 
## Lambda: -0.34423 LR test value: 264.24 p-value: < 2.22e-16 
## Numerical Hessian standard error of lambda: 0.036776 
## 
## Log likelihood: -119.6468 
## ML residual variance (sigma squared): 2129.1, (sigma: 46.142)
## Number of observations: 281 
## Number of parameters estimated: 6 
## AIC: NA (not available for weighted model)

Comentarios a destacar.

• La creación de los objetos vecinos en el análisis en datos en área juegan un papel muy importante puesto que todo parte desde aqui para su posterior modelación.

• Antes de realizar un método para medir la autocorrelación espacial hay que definir correctamente los enlaces entre cada vecino.

• En este modelamiento de datos en áreas se observa que la medición de las medidas de las distancias (longitud) no es algo tan relevante ,la conversión a UTM pasa a un segundo plano, algo muy particular, puesto que desde el inicio esto era algo bastante delicado, y se tenía como unos de los ítems más importantes antes de conocer un estudio, ya que dependía en gran medida la interpretación del modelado