# Configurar: borrar todo y configurar el directorio de trabajo
rm(list=ls())

Resumen

Este artículo tiene como objetivo replicar los ejemplos del capítulo 9 del libro Aplied Spatial Data Analysis With R (segunda edición), elaborado por Roger S. Bivand • Edzer Pebesma • Virgilio Gómez-Rubio, el cuál está asociado a los modelos de datos en área.

9.1. Introducción

Los datos espaciales a menudo se observan en entidades poligonales con límites definidos. Los límites de los polígonos son definidos por el investigador en algunos campos de estudio, pueden ser arbitrarios en otros y pueden crearse límites administrativos para propósitos muy diferentes en otros de nuevo. Los datos observados son frecuentemente agregaciones dentro de los límites, como recuentos de población. Las ebtidades areales pueden constituir por sí mismas las unidades de observación, por ejemplo cuando se planea estudiar el comportamiento del gobierno local donde las decisiones se toman a nivel de la entidad, por ejemplo estableciendo tasas impositivas locales. En general, sin embargo, las entidades areales son agregados, contenedores, utilizados para contar medidas, como resultados de votaciones en los colegios electorales. Muy a menudo, las entidades areales son un teselado exhaustivo del área de estudio, sin dejar ninguna parte del área total sin asignar a una entidad. Por supuesto, las entidades areales pueden estar formadas por múltiples entidades geométricas, tales como islas pertenecientes a un mismo condado; también pueden rodear otras entidades completamente, y pueden contener agujeros, como lagos.

Los límites de las entidades de área pueden definirse para algún otro propósito con su uso en el análisis de datos. Las áreas de códigos postales pueden ser útiles para el análisis, pero se crearon para facilitar las entregas postales. Es solo recientemente que Las organizaciones nacionales de censos han aceptado que los cambios frecuentes, aparentemente justificados, en los límites son un problema importante para los análisis longitudinales. Los límites arbitrarios de las unidades de área son un problema si su modificación puede conducir a resultados diferentes, como el caso de la manipulación política. siendo un recordatorio aleccionador de cómo los cambios en la agregación pueden cambiar los resultados, también pueden obstaculizar el análisis si la escala espacial o la huella de un proceso de generación de datos subyacente no coincide con los límites elegidos.

Si la recopilación de datos se puede diseñar para hacer coincidir las entidades de área con los datos, se reducirá la influencia de la elección de los áridos. Un ejemplo podría ser el cotejo de los datos del mercado laboral con los mercados laborales locales, tal vez definido por los viajes al trabajo. Por otra parte, si nos vemos obligados a utilizar arbitrariamente fronteras, a menudo porque no hay otras fuentes viables de datos secundarios, debemos ser conscientes de las posibles dificultades. Tales desajustes se encuentran entre las razones para encontrar autocorrelación espacial en el análisis de agregados de áreas; otras razones incluyen procesos espaciales sustantivos en los que las entidades influyen entre sí por contagio, como la adopción de políticas similares por parte de los vecinos, y la especificación errónea del modelo dejando información con patrones espaciales en los residuos del modelo. Estas causas de la autocorrelación espacial observada pueden ocurrir en combinación, haciendo la identificación correcta de la espacialidad real procesada a una empresa interesante; Dray et al. (2012) discuten muchos de los cuestiones involucradas en la realización de dicha identificación.

Una amplia gama de disciplinas científicas se han encontrado con la autocorrelación espacial entre entidades de área, con el término “problema de Galton” utilizado en muchas ocaciones. El problema es establecer cuántas observaciones efectivamente independientes están presentes, cuando se han utilizado límites arbitrarios para dividir un estudio de área. En su intercambio con Tyler en 1889, Galton cuestionó si las observaciones de las leyes matrimoniales en entidades del área constituían observaciones independientes, ya que podían reflejar un patrón general del cual se habían basado. todos descendidos. De modo que la dependencia espacial positiva tiende a reducir la cantidad de información contenida en las observaciones, porque las observaciones próximas pueden usarse en parte para predecirse entre sí.

En el cap. 8, hemos visto cómo se pueden usar las distancias en una superficie continua para estructurar la autocorrelación espacial, por ejemplo con el variograma. Aquí nos ocuparemos de las entidades de área que se definen como vecinos, por elección. Definiciones de vecinos. En una superficie continua, todos los puntos son vecinos. unos de otros, aunque algunos pueden tener muy poco peso, porque son muy distante En una superficie teselada, podemos elegir definiciones vecinas que dividen el conjunto de todas las entidades (excluyendo la observación i) en miembros o no miembros del conjunto vecino de observación i. También podemos decidir dar a cada relación vecina un peso igual, o variar los pesos en el arco del gráfico dirigido que describen la dependencia espacial.

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. También nos interesará mostrar cómo se puede introducir la autocorrelación espacial en datos, de manera que se puedan realizar simulaciones. Vale la pena tener en cuenta que el patrón espacial que encontramos solo puede indicar que nuestro modelo actual de los datos no son apropiados. Esto se aplica a las unidades de área que no se ajustan a los datos. proceso de generación, a las variables que faltan, incluidas las variables con la mal forma funcional, y las diferencias entre nuestras suposiciones acerca de los datos y sus distribuciones reales, a menudo mostradas como una dispersión excesiva en los datos de conteo.

# Paquetes...
library(rgdal)
library(spdep)
library(RColorBrewer)
library(boot)
library(spatialreg)

# Cargando SpatialPointsDataframe - SpatialPolygonsDataFrame - SpatialPointsDataFrame

# 8 Condados de Nueva York Superior y ubicaciones de 11 sitios de desechos peligrosos inactivos.
NY8 <- readOGR(".", "NY8_utm18")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\jdbul\OneDrive\Escritorio\UNIVERSIDAD\Estadistica Espacial\Datos-En-Area", layer: "NY8_utm18"
## with 281 features
## It has 17 fields

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

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\jdbul\OneDrive\Escritorio\UNIVERSIDAD\Estadistica Espacial\Datos-En-Area", layer: "TCE"
## with 11 features
## It has 5 fields

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

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\jdbul\OneDrive\Escritorio\UNIVERSIDAD\Estadistica Espacial\Datos-En-Area", layer: "NY8cities"
## with 6 features
## It has 1 fields

# Como se almacenan los datos?
class(NY8)

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

# Es un SpatialPolygonsDataFrame - un data.frame en polígono
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 "SpatialVector", by class "SpatialPolygons", distance 2, with explicit coerce

# Polígonos espaciales contiene polígonos y características espaciales
getClass("SpatialPolygons")

## Class "SpatialPolygons" [package "sp"]
## 
## Slots:
##                                                       
## Name:     polygons   plotOrder        bbox proj4string
## Class:        list     integer      matrix         CRS
## 
## Extends: "Spatial", "SpatialVector"
## 
## 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
# una línea es una lista ordenada de coordenadas
getClass("Line")

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

# finalmente... Spatial es la clase madre de todas las clases Spatial
# usado 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 cargados

NY8

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

Algunos gráficos

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)

El conjunto de datos de 281 distritos censales para ocho condados centrales del estado de Nueva York destacado en Waller y Gotway (2004) se utilizará en muchos de los ejemplos, complementado con límites de tracto derivados de TIGER 1992 y distribuido por SEDAC/CIESIN. Este archivo no es idéntico a los límites usados en la fuente original, pero es muy cercano y puede ser redistribuido, a diferencia de la versión usada en el libro. El área tiene una extensión de unos 160 km de norte a sur y 120 km de este a oeste; La gráfica anterior muestra los principales ciudades del área de estudio y la ubicación de 11 sitios de desechos peligrosos. Las cifras de Waller y Gotway (2004) incluyen cuerpos de agua, que no son presente en esta versión de los límites del tracto, en la que los límites del tracto siguen las líneas centrales de los lagos, en lugar de sus orillas.

Tracemos una de las características - porcentaje de edad > 65

spplot(NY8, c("PCTAGE65P"))

col=“transparent”

spplot(NY8, c("PCTAGE65P"), col="transparent")

Hagamos una trama diferente, con una nueva paleta de colores.

library("RColorBrewer")
#función de creador de paleta de colores
rds <- colorRampPalette(brewer.pal(8, "RdBu"))
#obtener un rango para los valores
tr_at <- seq(min(NY8$PCTAGE65P), max(NY8$PCTAGE65P), length.out=20)
#crear una función de interpolación de color tomando el requerido
#número de sombras como argumento
tr_rds <- rds(20)
#parámetros
# en - en qué valores cambian los colores
# col.regions - especificar colores de relleno
tr_pl <- spplot(NY8, c("PCTAGE65P"), at=tr_at, col="transparent",
                col.regions=tr_rds, main=list(label="Age>65", cex=0.8))
plot(tr_pl)

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. El primer paso es definir a qué relaciones entre las observaciones se les va a dar un peso distinto de cero, es decir, elegir el criterio vecino a utilizar; el segundo paso es asignar pesos a los enlaces vecinos identificados. tratando de detectar patrón en mapas de residuales visualmente no es una opción aceptable, aunque a veces se escuchan comentarios que explican la falta de análisis formal como “parecían aleatorios” o, alternativamente, “puedo ver los grupos”. Haciendo los vecinos y pesos no es, sin embargo, fácil de hacer, por lo que se incluyen varias funciones en el paquete spdep para ayudar. Se encuentran otras funciones en algunos paquetes ecológicos, como el paquete ade4, este paquete también proporciona convertidores nb2neig y neig2nb para la interoperabilidad. La construcción de los pesos espaciales es mencionado por Cressie (1993, pp. 384-385), Schabenberger y Gotway (2005, pág. 18), Waller y Gotway (2004, págs. 223–225), Fortin y Dale (2005, págs. 113 a 118), O’Sullivan y Unwin (2010, págs. 200 a 204), Ward y Gleditsch (2008, págs. 14 a 22), Chun y Griffith (2013, págs. 9, 57 a 59) y Banerjee et al. (2004, págs. 70 y 71). La escasez de tratamientos en la literatura contrasta con la fuerza de la información previa que está introduciendo el analista en esta etapa. 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. Las viñetas: “nb”, “CO69” y “sids” en spdep incluyen discusiones sobre la creación y el uso de pesos espaciales, y se puede acceder a través de:

# vignette
vignette("nb", package = "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. Es una lista de longitud n con los números índice de vecinos de cada componente registrado como un vector entero. 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. Los identificadores de región pueden ser utilizados para verificar la integridad entre los datos mismos y el objeto vecino. 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.

A partir de las contigüidades de la reina del distrito censal, donde todos los polígonos que se tocan son vecinos, utilizada en Waller y Gotway (2004) y proporcionada como DBF en su sitio web, se ha creado y leído un archivo en formato GAL R.

Como ahora tenemos un objeto nb para examinar, podemos presentar el método estándar para estos objetos. Existen métodos de impresión, resumen, trama y otros; el método de resumen presenta una tabla de la distribución del número de enlace, y tanto el método de impresión como el de resumen informan asimetría y la presencia de observaciones no vecinas; la asimetría está presente cuando yo soy un vecino de j pero j no es vecino de i. El siguiente gráfico muestra el vecino completo, gráfico para el área de estudio de los ocho condados. 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. Recuperamos la parte de la lista de vecinos en Syracuse usando el método de subconjuntos.

# lee un archivo de celosía GAL en una lista de vecinos
NY_nb <- read.gal("NY_nb.gal", region.id=row.names(NY8))

¿Qué estados son vecinos?

# ¿Qué 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)

Para usar más adelante en esta discusión, 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

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)

El objeto k = 1 también es útil para encontrar la distancia mínima a la que todas las áreas tienen un vecino basado en la distancia. Usando la función nbdists, podemos calcular una lista de vectores de distancias correspondientes al objeto vecino, aquí para los primeros vecinos más cercanos. El mayor valor será la mínima distancia necesaria para asegurarse de que todas las áreas estén vinculadas al menos a un vecino, por lo que podemos crear fácilmente un objeto con observaciones no vecinas estableciendo el umbral superior por debajo de este valor.

9.2.2 Objetos de pesos espaciales

La literatura sobre pesos espaciales es sorprendentemente pequeña, dada su importancia para medir y modelar la dependencia espacial en datos en área. Griffith (1995) proporciona buenos consejos prácticos, mientras que Bavaud (1998) busca insertar fundamentos conceptuales bajo pesos espaciales ad hoc. Los pesos espaciales pueden verse como una lista de pesos indexados por una lista de vecinos, donde el peso del enlace entre i y j es el k-ésimo elemento del i-ésimo componente de la lista de pesos, y k nos dice cuál de los i-ésimos valores de los componentes de la lista de vecinos es igual a j. Si j no está presente en el i-ésimo componente de la lista de vecinos, j no es vecino de i. En consecuencia, algunos pesos wij en la representación de la matriz de pesos W establecido en cero, donde j no es un vecino de i. Aquí, seguimos a Tiefelsdorf et al. (1999) en nuestro tratamiento, usando sus estilos de abstracción de pesos espaciales.

Una vez establecida la lista de conjuntos de vecinos para nuestra área de estudio, se procede a asignar pesos espaciales a cada relación. Si sabemos poco sobre el proceso espacial asumido, tratamos de evitar alejarnos de la representación binaria de un peso de unidad para los vecinos (Bavaud, 1998), y cero de lo contrario. En esta sección, revisamos las formas en que pondera los objetos: listw objetos – se construyen; la clase es una clase de estilo antiguo. A continuación, se mostrará la conversión de estos objetos en representaciones matriciales densas y dispersas, concluyendo con las funciones de importación y exportación.

La función nb2listw toma un objeto de lista de vecinos y lo convierte en un objeto de pesos. El estilo de conversión predeterminado es W, donde los pesos para cada entidad areal están estandarizados para sumar a la unidad; esto también se llama a menudo fila de Estandarización. El método de impresión para objetos listw muestra las características de los vecinos subyacentes, el estilo de los pesos espaciales y la constante de ponderación espacial utilizada en el cálculo de pruebas de autocorrelación espacial. El componente objeto vecino es el objeto nb subyacente, que da la indexación del componente de pesos.

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"

Para style=“W”, los pesos varían entre la unidad dividida por el mayor y los números más pequeños de vecinos, y las sumas de los pesos para cada entidad de área son unidad. Este estilo de ponderaciones espaciales puede interpretarse como que permite el cálculo de valores promedio entre vecinos. Los pesos de los enlaces que se originan en áreas con pocos vecinos son más grandes que los que se originan en áreas con muchos vecinos, tal vez impulsando entidades areales en el borde del estudio de área sin querer. Esta representación ya no es simétrica, pero es similar a la simétrica.

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

Establecer style = “B” - ‘binary’ - retiene un peso de unidad para cada vecino relación, pero en este caso, las sumas de pesos para las áreas difieren según a la cantidad de áreas vecinas que tienen.

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

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 se puede usar para pasar una lista de vectores de pesos generales correspondiente a las relaciones vecinas a nb2listw. Decir que creemos que la fuerza de las relaciones de vecindad se atenúa con la distancia, uno de los casos considerados por Cliff y Ord (1981, pp. 17-18) y O’Sullivan y Unwin (2010, pp. 200–204) brindan una discusión similar. Podríamos establecer qu los pesos sean proporcionales a la inversa de la distancia entre los puntos que representan las áreas, usando nbdists para calcular las distancias para el objeto nb dado. Usando lapply para invertir las distancias, podemos obtener una estructura diferente de pesos espaciales de los anteriores. Si no tenemos motivos para suponer más conocimiento sobre las relaciones vecinales que su existencia o ausencia, este paso es potencialmente engañoso. Si sabemos, en cambio, que la migración o los flujos de desplazamiento describen la estructura de los pesos espaciales mejor que los de alternativa binaria, puede valer la pena usarlos como pesos generales; sin embargo, habrá problemas de simetría, porque tales flujos, a diferencia de las distancias inversas – rara vez son simétricos

El siguiente gráfico muestra tres representaciones de ponderaciones espaciales para Syracuse mostradas como matrices. La imagen style=“W” de la izquierda es evidentemente asimétrica, con colores más oscuros mostrando mayores pesos para zonas con pocos vecinos. los otros dos paneles son simétricos, pero expresan diferentes supuestos sobre la fortalezas de las relaciones vecinales.

El argumento final de nb2listw nos permite manejar listas de vecinos con zonas de no vecinos. No es obvio que la representación de peso del conjunto vacío es cero; tal vez debería ser NA, lo que generaría problemas luego

Por este motivo, el valor predeterminado del argumento es zero.policy=FALSE, dando lugar a un error cuando se da un argumento nb con áreas sin vecinos. Establecer el argumento en VERDADERO permite la creación del espacio objeto de pesos, con pesos cero. El argumento zero.policy posteriormente deberá usarse en cada función llamada, a menos que se establezca en TRUE para la función actual. La Sesión R con set.ZeroPolicyOption(TRUE). El contraste entre la comprensión basada en conjuntos de los vecinos y la conversión a una representación matricial es discutido por Bivand y Portnov (2004), y se reduce a si el producto de los pesos de un área no vecina y un vector n arbitrario debe ser un valor faltante o un cero numérico. Como veremos más adelante, manteniendo las entidades de área no vecinas, plantea preguntas sobre el tamaño relevante de n teniendo en cuenta pruebas de autocorrelación, entre otras cuestiones

Tres representaciones de pesos espaciales para Syracuse

# Tres representaciones de pesos espaciales para Syracuse
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()

par(oopar)

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

El problema paralelo de los conjuntos de datos con valores faltantes en las variables, pero con pesos espaciales completamente especificados, se aborda a través de subset.listw, método que vuelve a generar los pesos para el subconjunto dado de áreas, para ejemplo dado por complete.cases. Sabiendo qué observaciones están incompletas, los vecinos y pesos subyacentes se pueden dividir en subconjuntos en algunos casos, con el fin de evitar la propagación de valores NA al calcular valores retrasados espacialmente. Muchas pruebas y funciones de ajuste de modelos pueden llevar esto afuera internamente si se establece el indicador de argumento apropiado, aunque el analista preferirá dividir en subconjuntos los datos de entrada y los pesos antes de probar o modelar.

Comentarios secc 9.1 - 9.2.2.

Considerando una superficie continua, todos los puntos son vecinos, unos de otros, a pesar de que algunos puedan tener muy poco peso debido a la distancia en una superficie teselada.
La creación de los objetos vecinos en el análisis de datos en área, es un paso fundamental para modelar dichos datos, pues a éstos, son asignados enlaces entre sí, con el fin de retribuirles pesos, que serán de gran utilidad para aplicar test de autocorrelación espacial y proceder a la construcción de modelos autorregresivos simultáneos, condicionales o de regresión espacial.
Los pesos espaciales juegan un papel fundamental para medir y modelar la dependencia espacial en datos en área.
Es importante tener construidos los objetos vecinos y los pesos a sus respectivos enlaces antes de aplicar métodos para medir la autocorrelación espacial.

9.2.3. Manejo de objetos de pesos espaciales

Hay varios paquetes contribuidos que brindan soporte para matrices dispersas, entre los cuales Matrix es un paquete recomendado. El as_dgRMatrix_listw wrapper convierte un objeto listw en una matriz dispersa de matriz ordenada, comprimida y orientada a filas, como un objeto dgRMatrix, una subclase de la matriz virtual Clase RsparseMatrix. Es más fácil hacer una matriz dispersa orientada por filas a partir de un objeto de ponderaciones espaciales, ya que las ponderaciones están orientadas por filas. Una función que se usa mucho en las funciones de prueba y ajuste de modelos es listw2U, que devuelve un objeto listw simétrico que representa el 1/2 (W +WT) espacial de la matriz de pesos

Sy14_nb <- read.gal("Sy_GeoDa1.GAL")
isTRUE(all.equal(Sy0_nb, Sy14_nb, check.attributes=FALSE))

## [1] TRUE

Sy16_nb <- read.gwt2nb("Sy_GeoDa4.GWT")
isTRUE(all.equal(Sy10_nb, Sy16_nb, check.attributes=FALSE))

## [1] TRUE

9.2.4 Uso de pesos para simular Autocorrelación Espacial

# Using Weights to Simulate Spatial Autocorrelation

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

# Simulación de autocorrelación espacial: diagramas de retraso espacial,
# mostrando una línea más suave ponderada localmente
oopar <- par(mfrow=c(1,2), mar=c(4,4,3,2)+0.1)
plot(uncorr_x, lag(Sy0_lw_W, uncorr_x), xlab="", cex.lab=0.8,
 ylab="spatial lag", main="Uncorrelated random variable", cex.main=0.8)
lines(lowess(uncorr_x, lag(Sy0_lw_W, uncorr_x)), lty=2, lwd=2)
plot(autocorr_x, lag(Sy0_lw_W, autocorr_x),
 xlab="", ylab="",
 main="Autocorrelated random variable", cex.main=0.8, cex.lab=0.8)
lines(lowess(autocorr_x, lag(Sy0_lw_W, autocorr_x)), lty=2, lwd=2)

par(oopar)

9.3. Prueba de autocorrelación espacial

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

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

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

et <- coords[,1] - min(coords[,1])
trend_x <- uncorr_x + 0.00025 * et
moran_t <- moran.test(trend_x, listw=Sy0_lw_W)
moran_t1 <- lm.morantest(lm(trend_x ~ et), listw=Sy0_lw_W)
# K <- moran(NY8$Cases, listw=nb2listw(NY_nb, style="B"), n=length(NY8$Cases), S0=Szero(nb2listw(NY_nb, style="B")))$K

9.3.1 Pruebas globales

El primer ejemplo es probar el número de casos por tramo censal (siguiendo a Waller y Gotway, 2004, p. 231) para la autocorrelación utilizando el valor predeterminado estilo de ponderaciones espaciales de estandarización de filas, y uso de la suposición de aleatorización analítica al calcular la varianza de la estadística. El resultado, como vemos, es que el patrón espacial de la variable de interés es significativo, y es muy probable que los tramos vecinos tengan valores similares.

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

Cambiando el estilo de los pesos espaciales para que todos los pesos sean iguales y sumando el número de observaciones, vemos que la probabilidad resultante, el valor se reduce unas 20 veces; recordamos que la estandarización por filas favorece observaciones con pocos vecinos, y que los estilos ‘B’, ‘C’ y ‘U’ ‘pesoup’ observaciones con muchos vecinos. En este caso, el estilo ‘S’ desciende entre ‘C’ y ‘W’.

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

De forma predeterminada, moran.test utiliza la suposición de aleatorización, que difiere de la suposición de normalidad más simple al introducir un término de corrección basado en la curtosis de la variable de interés. Cuando la curtosis corresponde al de una variable distribuida normalmente, los dos supuestos arrojan la misma varianza, pero a medida que la variable se aleja de la normalidad, la suposición de aleatorización compensa aumentando la varianza y disminuyendo la desviación estándar. En este caso, hay poca diferencia y los dos devuelven resultados similares.

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

Es útil mostrar aquí que la prueba estándar bajo normalidad está de hecho, la misma prueba que la prueba de Moran para los residuos de regresión del modelo, incluyendo sólo el intercepto. Hacer esta conexión aquí muestra que podríamos introducir variables adicionales en el lado derecho de nuestro modelo, una y otra vez por encima del intercepto, y potencialmente otras formas de manejar errores de especificación.

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

Usando la misma construcción, también podemos usar una aproximación de Saddlepoint en lugar del supuesto analítico normal (Tiefelsdorf, 2002), y una prueba exacta (Tiefelsdorf, 1998, 2000; Hepple, 1998; Bivand et al., 2009). Estos métodos son sustancialmente más exigentes desde el punto de vista computacional y originalmente se consideraron poco prácticos. Para conjuntos de datos de tamaño moderado como el que estamos usando, sin embargo, necesitamos menos del doble del tiempo requerido para alcanzar un resultado. En general, los métodos exactos y Saddlepoint hacen poca diferencia para resultados, para pruebas globales cuando el número de entidades espaciales no es pequeño, como aquí, con el valor de probabilidad cambiando solo por un factor de 2. Vemos más tarde que el impacto de las diferencias entre el supuesto de normalidad y la aproximación de Saddlepoint y la prueba exacta es más fuerte para los indicadores locales de asociación espacial.

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

También podemos usar una prueba de Monte Carlo, una prueba de arranque de permutación, en que los valores observados se asignan aleatoriamente a las zonas, y la estadística de interés calculado nsim tiempos. Dado que, en el caso global en contraste con el local, tenemos suficientes observaciones y podemos repetir esta permutación potencialmente muchas veces sin repetir.

# Monte Carlo test
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

Waller y Gotway (2004, p. 231) también incluyen una constante de riesgo de Poisson de evaluación de arranque paramétrico de la importancia de la autocorrelación en el caso cuenta. La tasa global constante r se calcula primero y se usa para crear recuentos esperados para cada sector censal multiplicándolos por la població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)


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

Los recuentos esperados también se pueden expresar como los valores ajustados de un nulo. Regresión de Poisson con una compensación establecida en el logaritmo de la población del tracto con un enlace de registro, esto muestra la relación con los modelos lineales generalizados (debido a que los Casos no son todos enteros, se generan advertencias):

Estos recuentos esperados rni se transmiten al argumento lambda a rpois para generar los conjuntos de datos sintéticos mediante el muestreo de la distribución de Poisson. El valor de probabilidad de salida se calcula a partir de la misma que observó el I de Moran menos la media de los valores I simulados, y lo dividió por su desviación estándar. La siguiente gráfica corresponde a Waller y Gotway (2004,pags. 232, Fig. 7.8), con las simulaciones paramétricas cambiando la distribución de Moran’s I hacia la derecha, porque está tomando en cuenta el impacto de la heterogénea, las poblaciones del tracto.

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)

fitted glm

rni <- fitted(glm(Cases ~ 1 + offset(log(POP8)), data=NY8, family="poisson"))

Hay una versión del I de Moran adaptada para usar una tasa Empirica de Bayes por Assun¸c˜ao y Reis (1999) que, a diferencia de los resultados de la tasa anteriores, se contrae de forma extrema en tasas para tramos con pequeñas poblaciones en riesgo hacia la tasa para el área como un todo, también utiliza métodos de Monte Carlo para la inferencia:

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

## 
##  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

El panel derecho en el siguiente gráfico presenta la salida de la función correlog en el paquete pgirmess de Patrick Giraudoux; la función es un envoltorio para dnearneigh y moran.test. La función selecciona automáticamente la distancia. bandas de casi 10 km, que abarcan toda el área de estudio. En este caso, la primera dos bandas de 0-10 y 10-20 km tienen valores significativos.

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

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

Correlogramas: (izquierda) valores del I de Moran para ocho órdenes de retraso sucesivos de vecinos contiguos; (derecha) valores del I de Moran para una secuencia de banda de distancia pares de vecinos

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

par(oopar)

El panel derecho en la gráfica anterior presenta la salida de la función correlog en el paquete pgirmess de Patrick Giraudoux; la función es un envoltorio para dnearneigh y moran.test. La función selecciona automáticamente la distancia. bandas de casi 10 km, que abarcan toda el área de estudio. En este caso, las primeras dos bandas de 0-10 y 10-20 km tienen valores significativos

print(cor8, p.adj.method="holm")

## Spatial correlogram for NY8$Cases 
## method: Moran's I
##            estimate expectation    variance standard deviate Pr(I) two sided
## 1 (281)  0.11038740 -0.00357143  0.00127922           3.1862        0.010090
## 2 (281)  0.09511336 -0.00357143  0.00056381           4.1561        0.000259
## 3 (281)  0.01671070 -0.00357143  0.00034815           1.0870        0.831113
## 4 (281)  0.03750596 -0.00357143  0.00025549           2.5699        0.061036
## 5 (281)  0.02691999 -0.00357143  0.00020312           2.1394        0.129603
## 6 (281)  0.02642797 -0.00357143  0.00017498           2.2679        0.116677
## 7 (281)  0.00934133 -0.00357143  0.00017186           0.9850        0.831113
## 8 (281)  0.00211893 -0.00357143  0.00019703           0.4054        0.831113
##            
## 1 (281) *  
## 2 (281) ***
## 3 (281)    
## 4 (281) .  
## 5 (281)    
## 6 (281)    
## 7 (281)    
## 8 (281)    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

9.3.2 Pruebas 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. para la definición de vecino elegida. Debido a esto, podemos romper de forma global medidas en sus componentes y, por extensión, construir pruebas localizadas destinadas a detectar ‘clusters’ (observaciones con vecinos muy similares) y ‘hotspots’ (observaciones con vecinos muy diferentes). Estas son discutidos brevemente por Schabenberger y Gotway (2005, pp. 23-25), y en mayor extensión por O’Sullivan y Unwin (2010, pp. 216-226), Waller y Got way (2004, pp. 236-242), Chun y Griffith (2013, pp. 94-98) y Fortin y Dale (2005, págs. 153–159). Están cubiertos con cierto detalle por Lloyd (2007, pp. 65–70) en un libro que se concentra en modelos locales; otros ejemplos son dado por Bivand (2010, pp. 236–244).

Primero, examinemos un diagrama de dispersión de Moran de la variable de conteo de casos de leucemia; para una discusión de los diagramas de dispersión de Moran, ver Ward y Gleditsch (2008,págs. 24–27). La siguiente gráfica por convención coloca 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. El I de Moran global es de relación lineal entre estos y se dibuja como una pendiente. La gráfica se divide aún más en cuadrantes en los valores medios de la variable y sus valores rezagados: bajo-bajo, bajo-alto, alto-bajo y alto-alto.

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=grey.colors(4, 0.95, 0.55, 2.2), bty="n", cex=0.8, y.intersp=0.8)
title("Tracts with influence")

par(oopar)

Tratar de detectar patrones locales residuales en presencia de espacial global la autocorrelación es difícil. Por esta razón, los resultados para la dependencia local no deben ser vistos como ‘absolutos’, sino que están condicionados al menos por el espacio espacial global. autocorrelación, y más generalmente por la posible influencia de los datos espaciales generando procesos en un rango de escalas desde global pasando por local hasta dependencia no detectada en la escala de las observaciones.

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"))

Waller y Gotway (2004, p. 239) amplían su hipótesis de riesgo de constante tratamiento a Moran’s Ii local, y podemos seguir su ejemplo:

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

La siguiente gráfica muestra los dos conjuntos de valores del Ii local de Moran, calculados en la forma estándar y utilizando el supuesto de Poisson para el riesgo de constante hipótesis. Ya sabemos que los I globales de Moran pueden variar en valor y en inferencia dependiendo de nuestras suposiciones - por ejemplo, que la inferencia debe tener en cuenta las desviaciones de nuestros supuestos distributivos. Lo mismo se aplica aquí a la suposición de la distribución de Poisson de que su media y desviación estándar son iguales, mientras que la dispersión excesiva parece ser un problema en datos que también muestran autocorrelación. Hay algunos cambios de signo entre los mapas, con los valores de la hipótesis de riesgo constante algo más lejos de cero.

Valores locales de Ii de Moran calculados directamente y usando la constante de riesgo (hipótesis)

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))

También podemos construir una prueba de Monte Carlo simple del riesgo constante (hipótesis) de los valores locales de Ii de Moran, simulando mucho como en el caso global, pero ahora conservando todos los resultados locales. Una vez finalizada la simulación, extraemos el rango del valor Ii de Moran local de riesgo constante observado para cada tramo, y calcular su valor de probabilidad para el número de simulaciones hecha. Usamos un enfoque paramétrico para simular los conteos locales usando el conteo esperado local como el parámetro para rpois, porque en el vecino los recuentos son muy bajos y hacen que la permutación sea imprudente. La Realización de la permutación para probar, usando todo el conjunto de datos también parece imprudente, porque entonces se estaría comparando iguales con diferentes (Siguiente gráfico).

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))

Valores de probabilidad para todas las secciones censales, Ii local de Moran: normalidad y suposiciones de aleatorización, aproximación de Saddlepoint, valores exactos y constante hipótesis de riesgo

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))

Valores de probabilidad para distritos censales en y cerca de la ciudad de Binghampton, Moran’s Ii local: supuesto de normalidad, valores exactos e 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))

9.4. Ajuste de modelos de datos de área

Hemos visto anteriormente que la falta de independencia entre las observaciones en datos espaciales - autocorrelación espacial - es un lugar común, y que las pruebas son disponibles. En un mundo ideal, uno preferiría recopilar datos en los que las observaciones fueran mutuamente independientes, y así evitar problemas en la inferencia. de los resultados analíticos. La mayoría de los analistas de datos aplicados, sin embargo, no tienen esta opción, y deben trabajar con los datos que están disponibles, o que pueden recolectarse con las tecnologías disponibles. Muy a menudo sucede que las observaciones sobre las covariables relevantes no están disponibles en absoluto, y que la detección de autocorrelación espacial en datos o residuos del modelo constituye de hecho el único camino que queda para modelar la variación restante

9.4.1 Enfoques de estadísticas espaciales

La dependencia espacial se puede modelar de diferentes maneras usando modelos estadísticos. En muchos casos, es común suponer que las observaciones son independientes e idénticamente distribuidas, pero esto puede no ser el caso cuando se trabaja con datos espaciales. Las observaciones no son independientes porque puede existir alguna correlación entre áreas vecinas. También puede ser difícil de separar el impacto de la autocorrelación espacial y las diferencias espaciales en la distribución de la observación. Cressie (1993, págs. 402–448, 458–477, 548–568) proporciona una discusión muy amplia de estos enfoques, incluyendo revisiones de los antecedentes para su desarrollo y completos ejemplos trabajados. Schabenberger y Gotway (2005, págs. 335–348) y Waller y Gotway (2004, pp. 362–380) se concentran en los modelos autorregresivos espaciales que se han utilizado en esta sección. Wall (2004) proporciona una revisión comparativa útil de la forma en que se modelan los procesos espaciales para los datos de área. Banerjee et al. (2004, pp. 79–87) también se centran en estos modelos, porque las características clave llevan hasta modelos jerárquicos. Fortin y Dale (2005, pp. 229–233) indican que los modelos autorregresivos espaciales pueden desempeñar un papel diferente en ecología, aunque revisiones como la de Dormann et al. (2007) sugieren que pueden ser de utilidad. En esta sección, hemos seguido a Waller y Gotway (2004, cap. 9) bastante de cerca, ya que sus ejemplos resaltan cuestiones como la transformación de la respuesta variable y usando pesos para tratar de manejar la heteroscedasticidad.

Desde un punto de vista estadístico, es posible dar cuenta de la correlación entre observaciones considerando una estructura del siguiente tipo en el modelo. Si el vector de variables de respuesta es normal multivariado, podemos expresar el modelo de la siguiente manera:

Y = μ + e,

donde μ es el vector de medias de área, que se puede modelar de diferentes maneras y e es el vector de errores aleatorios, que suponemos que se distribuye normalmente con media cero y varianza genérica V . A menudo se supone que la media depende en un término lineal en algunas covariables X, por lo que sustituiremos la media por XTβ en el modelo. Por otro lado, se toma la correlación entre áreas en cuenta considerando una forma específica de la matriz de varianza V .

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

La siguiente gráfica muestra la distribución espacial de los valores residuales para el estudio tramos censales del área. Las dos variables censales parecen contribuir a explicar la varianza en la variable respuesta, pero la exposición a CT no contribuye.

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

Antes de ver la estimación del modelo con más detalle, vale la pena mencionar el estimador de probabilidad de perfil aproximado introducido por Li et al. (2007) y perfeccionado en trabajos posteriores (Li et al., 2012). Suponiendo que se ha eliminado la tendencia de la variable de interés y que las ponderaciones espaciales están estandarizadas por filas, el estimador de un paso del parámetro de regresión espacial puede producir más información comparativa que la I de Moran. A modo de comparación, mostramos la estimación del parámetro de máxima verosimilitud ajustada utilizando los mismos datos:

NYlistwW <- nb2listw(NY_nb, style = "W")
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

Estos modelos se pueden estimar eficientemente por máxima verosimilitud. En R esto se puede hacer usando la función spautolm en el paquete spdep. El modelo se puede especificar usando una fórmula para el predictor lineal, mientras que la matriz W debe pasarse como un objeto listw. Para crear este objeto de la lista de vecinos podemos usar la función nb2listw, que tomará un objeto de clase nb, como se explicó anteriormente.

El siguiente código muestra cómo ajustar una autorregresión simultánea al modelo elegido. Hemos ajustado el modelo estándar y el modelo ponderado utilizando el tamaño de la población en 1980 (según el censo de EE. UU.) en las áreas como pesos. Esto reproduce el ejemplo desarrollado en Waller y Gotway (2004, cap. 9, pp. 375–379), y se remite al lector a su discusión para más información. En la llamada a nb2listw, especificamos style = “B” para construir W usando un indicador binario de vecindad

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: 564.21

De acuerdo con los resultados obtenidos parece que existe una importante correlación en los residuos porque el valor estimado de λ es 0.0405 y el valor p de la prueba de razón de verosimilitud es 0,022. En la prueba de razón de verosimilitud comparamos el modelo sin autocorrelación espacial (es decir, λ = 0) con uno que lo permita (es decir, el modelo ajustado con autocorrelación distinta de cero parámetro).

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

Izquierda): componente de tendencia de los valores ajustados del modelo SAR; (derecha): espacial componente estocástico de los valores ajustados del modelo SAR

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)

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 de proporciones de incidencia. La versión ponderada de estos modelos se puede instalar para que los distritos se ponderen proporcionalmente a la inversa del tamaño de su población. Para ello incluimos el parámetro weights=POP8 en la llamada a la función lm

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)

Comenzando con el modelo lineal ponderado, podemos ver que la variable de exposición a TCE se ha vuelto significativa con el signo esperado, lo que indica que las extensiones más cercanas a los sitios de TCE tienen proporciones de incidencia transformada ligeramente más altas. Las otras dos covariables ahora también tienen coeficientes más significativos. El panel derecho de la siguiente gráfica muestra que la información se ha desplazado del residuos del modelo al propio modelo, con poca estructura espacial restante visible en el mapa.

lmwresid, Residuos del modelo lineal ponderado de incidencia transformada dimensiones; Las ubicaciones de los sitios de TCE se muestran con fines comparativos

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( "lmwresid"), sp.layout=list(TCEpts), col.regions=gry, col="transparent", lwd=0.5, at=seq(-2,4.5,0.5))

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 <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME , data=NY8, listw=NYlistw, weights=POP8)
summary(nysarw)

## 
## Call: 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: 515.2

(Izquierda): componente de tendencia de los valores ajustados del modelo SAR ponderado; (derecha): Componente estocástico espacial de los valores ajustados del modelo SAR ponderado

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)

9.4.1.2 Modelos autorregresivos condicionales

Se remite al lector a Schabenberger y Gotway (2005, pp. 338–339) para obtener una descripción detallada de las especificaciones CAR. 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 establecemos 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: 563.66

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 que son propietarios 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 del valor de λ es 0.0841.

nylam1 <- c(nycar$lambda)
nycarw <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, family="CAR",listw=NYlistw, weights=POP8)
summary(nycarw)

## 
## Call: 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: 515.14

9.4.1.3 Ajuste de modelos de regresión espacial

La función spautolm ajusta los modelos de regresión espacial por máxima verosimilitud, encontrando primero el valor del coeficiente autorregresivo espacial, que maximiza la función de verosimilitud logarítmica para la familia de modelos elegida, y luego ajustando los otros coeficientes por mínimos cuadrados generalizados en ese punto. Esto significa que el coeficiente autorregresivo espacial se puede encontrar mediante la búsqueda de líneas utilizando optimizar, en lugar de optimizar sobre todos los parámetros del modelo en el Mismo tiempo.

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

## 
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8, 
##     listw = NYlistw, weights = POP8, family = "SAR", method = "Matrix")
## 
## 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.816730   0.576037  6.6258 3.453e-11
## PCTOWNHOME  -0.380777   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.013835 
## 
## Log likelihood: -251.6017 
## ML residual variance (sigma squared): 1104.1, (sigma: 33.229)
## Number of observations: 281 
## Number of parameters estimated: 6 
## AIC: 515.2

El resultado de ajustar el modelo SAR ponderado usando funciones del paquete de matrices es idéntico al que se obtiene al usar los valores propios de W. Si es de interés examinar los valores de la función logarítmica de verosimilitud para un rango de valores de λ, el argumento llprof puede usarse para dar el número de valores de λ igualmente espaciados que se elegirán entre el inverso del más pequeño y valores propios más grandes para method=“eigen”, o una secuencia de tales valores más en general.

1/range(eigenw(NYlistw))

## [1] -0.3029200  0.1549552

nysar_ll <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, family="SAR",listw=NYlistw, llprof=100)
nysarw_ll <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, family="SAR",listw=NYlistw, weights=POP8, llprof=100)

Valores de probabilidad de registro para un rango de valores de λ, ponderados y no ponderados modelos SAR; se muestran los valores máximos y los coeficientes espaciales ajustados

ylim <- range(c(nysarw_ll$llprof$ll, nysar_ll$llprof$ll), na.rm=TRUE)
plot(nysarw_ll$llprof$lambda, nysarw_ll$llprof$ll, type="l", xlab=expression(lambda), ylab="log likelihood", ylim=ylim, lwd=2)
abline(v=nysarw_ll$lambda)
abline(h=nysarw_ll$LL)
lines(nysar_ll$llprof$lambda, nysar_ll$llprof$ll, lty=2, lwd=2)
abline(v=nysar_ll$lambda, lty=2)
abline(h=nysar_ll$LL, lty=2)
legend("bottom", legend=c("weighted SAR", "SAR"), lty=c(1,2), lwd=2, bty="n")

La gráfica anterior muestra la forma de los valores de la función logarítmica de verosimilitud a lo largo del rango factible de λ para los modelos SAR ponderados y no ponderados. Podemos ver fácilmente que las curvas son muy planas en los máximos, lo que significa que podríamos desplazar λ bastante sin impactar mucho el valor de la función. La gráfica también muestra la fuerte caída en los valores de la función como la gran negativa para valores de la patada jacobiana cerca de los extremos del rango factible. Finalmente, family=“SMA” para modelos de promedios móviles simultáneos también esta disponible dentro del mismo marco general, pero siempre implica el manejo de matrices densas para el ajuste.

nysmaw <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=NY8, family="SMA",listw=NYlistw, weights=POP8)
summary(nysmaw)

## 
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = NY8, 
##     listw = NYlistw, weights = POP8, family = "SMA")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.487080 -0.268990  0.093956  0.466055  4.284087 
## 
## Coefficients: 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -0.795243   0.143749 -5.5321 3.163e-08
## PEXPOSURE    0.080153   0.028237  2.8386  0.004531
## PCTAGE65P    3.820316   0.575463  6.6387 3.165e-11
## PCTOWNHOME  -0.382529   0.156160 -2.4496  0.014302
## 
## Lambda: 0.0091843 LR test value: 0.30771 p-value: 0.57909 
## Numerical Hessian standard error of lambda: 0.016571 
## 
## Log likelihood: -251.6112 
## ML residual variance (sigma squared): 1105.3, (sigma: 33.245)
## Number of observations: 281 
## Number of parameters estimated: 6 
## AIC: 515.22

Modelado de datos en área