Datos en áreas
Cristina Álvarez Múnera
2023-06-28
Modelación de datos en áreas
Introducción
Los datos espaciales suelen ser más visibles en elementos poligonales con límites definidos. Los límites de estos los polígonos los definen los investigadores, sin embargo, en algunos campos de estudio pueden ser arbitrarios y en otros pueden ser límites administrativos creados para propósitos muy diferentes. Los datos observados a menudo se agrupan en límites, como por ejemplo la población. Las unidades territoriales pueden formar por sí mismas unidades de control, por ejemplo, examinando las acciones de las autoridades locales cuando se toman decisiones a nivel de unidad, como el establecimiento de tasas impositivas locales. Sin embargo, en general, las funciones de región son agregados, contenedores que se utilizan para calcular métricas como los resultados de encuestas. En la mayoría de los casos, un elemento de área es una subdivisión completa del área de estudio en la que cualquier parte del área total se asigna a una función. Desde luego, las características del área pueden constar de varios elementos geométricos, como islas que pertenecen a un solo condado; también, pueden cubrir completamente otras características del área y pueden contener agujeros como lagos.
Los límites de las áreas geográficas pueden ser establecidos con fines distintos a su utilización en el análisis de datos. Aunque las divisiones por códigos postales pueden resultar útiles en el análisis, fueron creadas originalmente para facilitar los envíos postales. Solo recientemente, las instituciones censales nacionales han reconocido que los cambios frecuentes e aparentemente justificados en los límites representan un problema importante para el análisis longitudinal. La naturaleza arbitraria de los límites de las unidades geográficas se convierte en un problema cuando su modificación puede dar lugar a resultados diferentes, como se ejemplifica en el caso del gerrymandering político, que nos sirve como un recordatorio revelador de cómo los cambios en la agrupación pueden alterar los resultados. Además, estos límites pueden dificultar el análisis si la escala espacial o la extensión de un proceso de generación de datos subyacente no se corresponden con los límites elegidos.
Si se diseña la recopilación de datos de manera que las entidades geográficas se alineen con los datos, se puede minimizar la influencia de la elección de los agregados. Un ejemplo de esto podría ser vincular los datos del mercado laboral con las áreas de mercado laboral local, que podrían definirse según los desplazamientos al trabajo. Sin embargo, cuando nos vemos obligados a utilizar límites arbitrarios, generalmente debido a la falta de otras fuentes de datos secundarios, debemos ser conscientes de las dificultades potenciales que esto conlleva. Estos desajustes son una de las razones por las cuales se encuentra autocorrelación espacial en el análisis de los agregados geográficos. Otras razones incluyen la influencia mutua entre las entidades debido a procesos espaciales contagiosos, como la adopción de políticas similares por parte de los vecinos, así como la especificación incorrecta del modelo, lo que puede dejar residuos con patrones espaciales. Identificar correctamente los procesos espaciales reales es una tarea interesante, ya que pueden ser causantes de la autocorrelación espacial observada.
Diversas disciplinas científicas han encontrado evidencia de autocorrelación espacial entre entidades geográficas, y este fenómeno se conoce como el “problema de Galton”, que ha sido abordado en varios contextos. Este problema se centra en determinar cuántas observaciones son verdaderamente independientes cuando se utilizan límites arbitrarios para dividir un área de estudio. En un intercambio de ideas con Tyler en 1889, Galton planteó la pregunta de si las observaciones de las leyes matrimoniales entre entidades geográficas constituían observaciones independientes, dado que podrían reflejar únicamente un patrón general de descendencia en todos los casos. Por lo tanto, la dependencia espacial positiva tiende a reducir la cantidad de información contenida en las observaciones, ya que las observaciones cercanas pueden ser utilizadas en parte para predecir unas a otras.
Se ha observado cómo las distancias en una superficie continua pueden ser utilizadas para establecer la autocorrelación espacial, como lo es el caso del variograma. En esta ocasión, nos centraremos en las entidades geográficas que se definen como vecinas. En una superficie continua, todos los puntos son vecinos entre sí, aunque algunos pueden tener un peso insignificante debido a su gran distancia. En una superficie discretizada, podemos elegir diferentes definiciones de vecindad que dividen el conjunto de todas las entidades (excluyendo la observación i) en miembros y no miembros del conjunto de vecinos de la observación i. Además, podemos decidir asignar a cada relación de vecindad un peso igual o variar los pesos en los arcos del grafo dirigido que describe la dependencia espacial.
En la siguiente sección, abordaremos la construcción de vecindarios y los pesos que se pueden aplicar a estas vecindades. Una vez que hayamos establecido este prerrequisito importante y a menudo exigente, procederemos a estudiar cómo medir la autocorrelación espacial. Si bien las pruebas se basan en modelos de procesos espaciales, primero examinaremos las pruebas y luego pasaremos a la modelización. También nos interesa explorar cómo podemos introducir la autocorrelación espacial en datos que inicialmente eran independientes, para así poder realizar simulaciones.
Es importante tener en cuenta que el patrón espacial que identifiquemos puede indicar que nuestro modelo actual de los datos no es adecuado. Esto puede deberse a unidades de área que no se ajustan al proceso de generación de datos, a variables faltantes, incluyendo aquellas con una forma funcional incorrecta, así como a las diferencias entre nuestras suposiciones sobre los datos y las distribuciones reales, las cuales a menudo se manifiestan como una dispersión excesiva en los datos de recuento.
Dataset
Se empleará los datos de censo para 8 condados centrales de New York.
NY8 <- readOGR("data", "NY8_utm18")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\Usuario\Videos\trabajo3\data", layer: "NY8_utm18"
## with 281 features
## It has 17 fieldsny <- st_read("data/NY8_utm18.shp")## Reading layer `NY8_utm18' from data source
## `C:\Users\Usuario\Videos\trabajo3\data\NY8_utm18.shp' using driver `ESRI Shapefile'
## Simple feature collection with 281 features and 17 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 358241.9 ymin: 4649755 xmax: 480393.1 ymax: 4808545
## Projected CRS: WGS 84 / UTM zone 18Nplot(st_geometry(ny), axes = T)
Vecinos y pesos espaciales
La creación de pesos espaciales es un paso necesario en el análisis de datos en área. Para ello, primeramente se debe definir cuales relaciones entre las observaciones tendrán pesos no nulos, es decir, se debe escoger el criterio para que dos enlaces sean considerados como “vecinos”; el segundo paso es asignar los pesos a los enlaces identificados como “vecinos” en el paso anterior.
Este prceso no es facil y por ello existen varios métodos y funciones
en el paquete spdep que ayudan a facilitar este trabajo.
Además, el paquete ade4 tiene las funciones
nb2neig() y neig2nb(), las cuales ayudan a
transformar objetos entre las clases nb y
neig.
Objetos vecinos
En el paquete spdep, las observaciones vinculadas o
vecinas se representan mediante objetos de clase nb. Estos
objetos se pueden visualizar como una lista de tamaño \(n\), donde cada elemento contiene un vector
con los índices de los vecinos correspondientes a esa observación. Si
alguna observación no tiene vecinos, la entrada correspondiente contiene
el valor cero (0). Además de estos índices de vecinos, los objetos nb
también pueden contener otros atributos, como un vector de caracteres
para identificar cada una de las regiones consideradas, y un vector
lógico (booleano) que indica si la relación entre observaciones es
simétrica o no (este aspecto se discute más adelante). La función
card() se utiliza para obtener el número de vecinos para
cada elemento de la lista (objeto nb). A continuación, se
presenta de manera concisa la estructura de estos objetos
espaciales.
NY_nb <- read.gal("data/NY_nb.gal", # read.gal carga el archivo de clase nb
region.id = as.numeric(row.names(ny)) - 1
) # Debe empezar en 0
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 linksplot(st_geometry(ny), border = "grey60", axes = T)
plot(NY_nb, st_geometry(ny), pch = 19, cex = 0.6, add = T)
Cabe resaltar, que un archivo gal es una matriz de pesos
o vecinos producida por el software GeoDa. Para esto, todas las regiones
contiguas son consideradas como “regiones vecinas”.
Una vez ilustrados los datos espaciales que van a ser utilizados, se
usará un subconjunto de ellos para mostrar el uso de las principales
funciones y métodos para los objetos nb, es decir, se usará
“Syracuse city”.
Syracuse <- filter(ny, AREANAME == "Syracuse city")
Sy0_nb <- subset(NY_nb, ny$AREANAME == "Syracuse city")
plot(st_geometry(Syracuse), border = "grey60", axes = T)
plot(Sy0_nb, st_geometry(Syracuse), pch = 19, cex = 0.6, add = T)
Se tiene que uno de los principales métodos de los objetos
nb son print(), summary() y
plot(), sin embargo, estos no son los únicos que existen;
tanto print como summary presentan el número
de regiones, enlaces no nulos, porcentaje de pesos no nulos y promedio
de enlaces entre regiones, sin embargo, summary presenta
adicionalmente la distribución del número de enlaces del objeto
nb, además, de las regiones con menor y mayor frecuencia de
enlaces dentro del objeto.
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 linksEl método summary también reporta la presencia o no de
asimetría en el objeto nb
La asimetria está presente cuando i es un vecino de j, pero j no es un vecino de i.
El método summary proporciona la siguiente
información:
- El número de regiones (polígonos) en la base de datos: 63
- El número de relaciones diferentes de cero: 346
- El porcentaje de pesos diferentes de cero: 8.717561%
- El número promedio de relaciones entre los polígonos: 5.492063
Además, se muestra la distribución de las relaciones, lo que indica que el número máximo de vecinos encontrados fue 9, y esto ocurrió en relación con un polígono específico. Además, se observa que 17 de las 63 regiones presentes tienen 6 conexiones, y así sucesivamente.
La función card permite obtener rápidamente el número de
vecinos de cada región.
card(Sy0_nb)## [1] 5 5 2 7 7 5 4 5 7 5 5 7 5 7 6 6 6 6 6 4 4 7 6 8 4 3 3 9 6 8 6 6 6 8 6 4 4 8
## [39] 6 6 7 5 7 6 4 5 3 5 6 4 6 7 4 8 5 1 5 8 5 5 6 3 3En el documento, se generarán objetos de lista de vecinos para su uso posterior. Para ello, se utiliza la función knearneigh que nos permite obtener los índices de los puntos que son los \(k\) vecinos más cercanos. Esta función recibe como parámetros la lista de coordenadas de los puntos y el valor de \(k\), que indica el número de puntos vecinos cercanos que se desean obtener.
Una vez obtenido el resultado anterior, se pasa como entrada a la función knn2nb, que convierte el objeto resultante en una lista de vecinos de la clase nb. Esta transformación nos permite obtener un objeto nb que contiene la información sobre los vecinos de cada punto en forma de una lista.
coords <- st_point_on_surface(st_geometry(Syracuse))
ids <- row.names(Syracuse)
Sy8_nb <- knn2nb(knearneigh(coords, k = 1), row.names = ids)
Sy8_nb## Neighbour list object:
## Number of regions: 63
## Number of nonzero links: 63
## Percentage nonzero weights: 1.587302
## Average number of links: 1
## Non-symmetric neighbours listLa función nbdists se utiliza para calcular las distancias entre los puntos de los objetos vecinos obtenidos en el paso anterior. Esta función es útil porque nos permite determinar la distancia mínima necesaria para garantizar que todas las áreas tengan al menos un vecino.
Después de calcular las distancias con nbdists, se emplea la función dnearneigh para obtener la lista de objetos vecinos (nb) que se encuentran a una distancia específica (d1, d2). Esto nos permite seleccionar los vecinos que se encuentran dentro de un rango determinado de distancias y obtener la lista correspondiente de vecinos para cada punto.
dsts <- unlist(nbdists(Sy8_nb, coords))
Sy11_nb <- dnearneigh(x = coords, d1 = 0, d2 = 0.75 * max(dsts))
Sy11_nb## Neighbour list object:
## Number of regions: 63
## Number of nonzero links: 156
## Percentage nonzero weights: 3.930461
## Average number of links: 2.47619
## 8 regions with no links:
## 1 10 20 36 46 57 60 61Como resultado, se observa que hay 8 regiones que no poseen vecinos (1, 10, 20…).
Objetos espaciales ponderados
Resulta sorprendente que la literatura sobre pesos espaciales sea limitada, considerando la gran importancia que estos tienen al medir y modelar la dependencia espacial en datos de áreas. Los pesos espaciales desempeñan un papel fundamental en el análisis espacial al cuantificar la influencia de las observaciones vecinas en un determinado fenómeno. Estos pesos permiten capturar patrones de dependencia espacial y son utilizados en diversas técnicas, como la autocorrelación espacial, la regresión espacial y la interpolación espacial.
Aunque la importancia de los pesos espaciales es ampliamente reconocida en la comunidad científica, la disponibilidad de métodos y herramientas para su aplicación y la exploración de sus propiedades sigue siendo limitada. Es posible que esta falta de atención en la literatura se deba a la complejidad y la naturaleza multidimensional de los datos espaciales, así como a los desafíos metodológicos asociados con la selección y el cálculo de pesos espaciales adecuados.
A medida que el análisis espacial y el interés en comprender los patrones espaciales continúan creciendo en diferentes disciplinas, es necesario fomentar una mayor investigación y desarrollo en el campo de los pesos espaciales. Esto permitirá mejorar nuestra comprensión de la dependencia espacial y promover el uso adecuado de los pesos espaciales en el análisis y modelado de datos de áreas.
Los pesos espaciales pueden ser vistos como una lista de pesos indexada por una lista de vecinos. En este enfoque, el peso de la relación entre dos regiones, \(i\) y \(j\), se encuentra en el \(k\)-ésimo elemento del componente \(i\) de la lista de pesos.
Una vez que se ha definido la lista de vecinos espaciales, el siguiente paso consiste en asignar pesos a los enlaces entre las regiones. Si se cuenta con información relevante sobre el proceso espacial que se está analizando, esta información puede ser utilizada para asignar pesos a los enlaces. Sin embargo, si no se dispone de esta información, es recomendable evitar hacer suposiciones sobre los pesos para evitar problemas de mala especificación del modelo.
La función nb2listw se utiliza para convertir un objeto de lista de vecinos (nb) en un objeto de pesos. La conversión por defecto es un objeto de pesos denominado W, donde los pesos varían entre la unidad dividida por el número máximo y mínimo de vecinos, y la suma de los pesos para cada región es igual a la unidad. Esta conversión permite establecer una medida adecuada de la influencia espacial entre las regiones.
Sy0_lw_w <- nb2listw(neighbours = 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.564Existe otro tipo de conversión que se puede obtener al configurar el parámetro \(style = B\). En este caso, se asigna un peso de 1 a cada relación vecina y un peso de 0 para las relaciones que no son vecinas. Sin embargo, es importante tener en cuenta que en esta configuración, la suma de los pesos para cada área puede variar según el número de vecinos que tenga cada área.
Al establecer style = B, se crea un objeto de pesos que refleja únicamente la presencia o ausencia de una relación vecina entre las áreas. Cada relación vecina se representa con un peso de 1, lo que indica que existe una conexión directa entre las áreas correspondientes. Por otro lado, las relaciones que no son vecinas tienen asignado un peso de 0, indicando que no existe una conexión directa entre esas áreas.
Es recomendable utilizar pesos binarios cuando se tiene poco conocimiento acerca del proceso que se está estudiando.
El argumento glist proporciona la opción de asignar manualmente los pesos a cada enlace entre regiones en el objeto nb. Esta función resulta útil en situaciones en las que se dispone de información específica sobre las posibles relaciones entre dos regiones, lo cual puede contribuir a obtener resultados más precisos y adecuados a la situación particular. A continuación se presenta un ejemplo que ilustra el uso del argumento glist, suponiendo que los pesos siguen la regla IDW binaria (inverso de las distancias entre las regiones).
dsts <- nbdists(Sy0_nb, coords)
idw <- lapply(dsts, function(x) 1 / (x / 1000))
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.3062 0.7481 0.9547 1.0154 1.2232 2.7686summary(sapply(Sy0_lw_idwB$weights, sum))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.321 4.129 5.896 5.577 6.967 9.860Se presentan 3 tipos de pesos espaciales:
Sy0_lw_B <- nb2listw(Sy0_nb, style = "B")
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)El último argumento a tener en cuenta en la función nb2listw() es zero.policy. Cuando hay regiones que no tienen vecinos, es necesario especificar el argumento zero.policy = TRUE para evitar que se produzca un error. Esto se debe a que, por defecto, la función nb2listw() asume que cada región tiene al menos un vecino.
Al establecer zero.policy = TRUE, se le indica a la función que permita regiones sin vecinos y que no genere un error en tales casos. Esto es útil cuando se trabaja con conjuntos de datos en los que algunas regiones pueden no tener vecinos debido a su ubicación geográfica o a otras características particulares.
try(
expr = {
nb2listw(neighbours = Sy11_nb, style = "B", zero.policy = F)
}
)## Error in nb2listw(neighbours = Sy11_nb, style = "B", zero.policy = F) :
## Empty neighbour sets foundLuego de especificar como verdadero dicho argumento, se obtiene el resultado, sin fallas computacionales:
Sy0_lw_D1 <- nb2listw(neighbours = Sy11_nb, style = "B", zero.policy = T)
print(Sy0_lw_D1, zero.policy = T)## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 63
## Number of nonzero links: 156
## Percentage nonzero weights: 3.930461
## Average number of links: 2.47619
## 8 regions with no links:
## 1 10 20 36 46 57 60 61
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 55 3025 156 312 2240Evaluando la autocorrelación espacial
Una vez que hemos construido diferentes tipos de pesos espaciales, podemos comenzar a utilizarlos para analizar la presencia de autocorrelación espacial. Sin embargo, antes de adentrarnos en un análisis más detallado, es importante revisar los supuestos subyacentes a las pruebas de autocorrelación espacial; utilizaremos la I de Moran como ejemplo, pero las consecuencias se aplican también a otras pruebas.
Existen dos tipos de test de autocorrelación espacial para este tipo de datos:
Test globales: Se utilizan para verificar si existe autocorrelación espacial entre dos o más regiones en algún lugar del área de estudio (no se especifica donde).
Test locales: Detectan clusters de vecinos con observaciones muy diferentes de los vecinos en general. Encuentran hotspots.
Test globales
Test de Moran I
Basicamete, es la forma más común de abordar la autocorrelación espacial. Se calcula como la razón del producto de la variable de interés y su rezago espacial, con el producto cruzado de la variable de interés y ajustada por los pesos espaciales usados.
\[I = \frac{n}{\sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij}}\ \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} (y_i - \bar{y})(y_j - \bar{y})}{\sum_{i=1}^{n} (y_i - \bar{y})^2}\] donde:
- \(y_i\): \(i\)-ésima observación
- \(\bar{y}\): media de la variable de interés
- \(w_{ij}\): ponderación espacial de las relaciones entre \(i\) e \(j\).
Al enfocarnos en la media, estamos asumiendo que el modelo adecuado tiene una media constante, y cualquier patrón residual que persista después de centrar los datos se debe a las relaciones espaciales codificadas en los pesos espaciales.
Al centrar los datos, eliminamos la influencia de la media en el análisis, lo que nos permite examinar y evaluar específicamente los patrones espaciales presentes en los datos. Si después de centrar los datos aún se observa un patrón espacial significativo, esto sugiere que hay una estructura espacial más allá de las variaciones en la media.
Se calcula la esperanza como:
\[E(I) = \frac{-1}{n-1}\]
En el siguiente paso, procedemos restando el valor observado al valor esperado analíticamente (I - E(I)), y luego dividimos esta diferencia por la raíz cuadrada de la varianza analítica. El resultado obtenido se compara con la distribución normal para determinar el valor de probabilidad del estadístico observado, bajo la suposición nula de ausencia de dependencia espacial para los pesos espaciales seleccionados. Es importante tener en cuenta que los resultados pueden variar según las elecciones de las ponderaciones, los vecindarios y el cumplimiento de los supuestos requeridos.
\[H_0: Ausencia \ de \ dependencia \ espacial\] \[H_1: El \ estadístico \ observado \ es \ significativamente \ más \ grande \ que \ el \ valor \ esperado\]
Es importante tener en ceunta que el dominio del I de Moran es [-1, 1]. Además, se evidencia autocorrelación espacial positiva cuando I > 0. Esto permite concluir que las unidades de análisis vecinas tienden a ser similares.
moran.test(x = ny$Cases, listw = nb2listw(NY_nb))##
## Moran I test under randomisation
##
## data: ny$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.001430595Por defecto, al utilizar moran.test, se aplica el supuesto de aleatorización, lo cual difiere del supuesto más simple de normalidad al introducir un término de corrección basado en la curtosis de la variable de interés. Este enfoque considera que la distribución de la variable puede desviarse de la normalidad, y ajusta la varianza en consecuencia.
Cuando la curtosis de la variable se asemeja a la de una distribución normal, ambas suposiciones resultan en la misma varianza. Sin embargo, a medida que la variable se aleja de la normalidad, la suposición de aleatorización compensa este hecho aumentando la varianza y disminuyendo la desviación estándar.
moran.test(x = ny$Cases, listw = nb2listw(NY_nb), randomisation = F)##
## Moran I test under normality
##
## data: ny$Cases
## weights: nb2listw(NY_nb)
##
## Moran I statistic standard deviate = 3.9731, p-value = 3.547e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.146882990 -0.003571429 0.001433975Se observa que la varianza no cambia, por lo tanto, el supuesto de normalidad se cumple.
Es posible también realizar test de Monte Carlo para evaluar la hipótesis de dependencia espacial con el test de Moran I.
moran.mc(x = ny$Cases, listw = nb2listw(NY_nb), nsim = 999)##
## Monte-Carlo simulation of Moran I
##
## data: ny$Cases
## weights: nb2listw(NY_nb)
## number of simulations + 1: 1000
##
## statistic = 0.14688, observed rank = 1000, p-value = 0.001
## alternative hypothesis: greaterWaller y Gotway (2004, p. 231) proponen un enfoque adicional para evaluar la autocorrelación en los recuentos de casos, utilizando un modelo de riesgo constante de Poisson. Este enfoque se basa en calcular una tasa global constante, denotada como \(r\), que se utiliza posteriormente para generar los recuentos esperados para cada unidad geográfica. Estos recuentos esperados se obtienen multiplicando la tasa global por la población correspondiente a cada área censal.
lw_B <- nb2listw(NY_nb, style = "B")
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.1472112Los recuentos esperados también pueden ser obtenidos como los valores ajustados de una regresión de Poisson, teniendo en cuenta un desplazamiento ajustado al logaritmo de la población de cada unidad geográfica. Este enfoque muestra una conexión con los modelos lineales generalizados. Sin embargo, es importante tener en cuenta que como los casos no siempre son números enteros, pueden generarse advertencias al aplicar este método.
Para generar los conjuntos de datos sintéticos, se utilizan los recuentos esperados rni como parámetros en la función rpois, que simula la distribución de Poisson. A partir de estos datos sintéticos, se calcula la probabilidad de salida, que se obtiene a partir de la diferencia entre el I de Moran observado y la media de los valores simulados de I, dividida por la desviación estándar de estos valores simulados.
# 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: greaterr <- 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.1472112oopar <- 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)Los resultados del índice de Bayes empírico indican que una posible explicación para la falta de significancia en el bootstrap paramétrico del riesgo constante observado y los valores esperados podría ser la presencia de valores inusuales y extremos en segmentos con poblaciones reducidas. Después de aplicar técnicas de suavizado a las tasas, se observa cierta autocorrelación a nivel global.
Test locales
Al utilizar las medidas globales de asociación espacial, se puede descomponer dichas medidas para crear pruebas locales que buscan identificar clusters (observaciones con vecinos altamente similares) y hotspots (observaciones con vecinos altamente diferentes). Estas pruebas también pueden utilizarse para evaluar el impacto de los estadísticos individuales en los indicadores globales correspondientes.
Indicadores locales de asociación espacial (LISA)
(Anselin, 1995)
Cualquier estadístico que cumpla las siguientes condiciones:
El LISA (Local Indicators of Spatial Association) para cada observación proporciona una medida del grado de agrupación espacial significativa de valores similares alrededor de dicha observación.
La suma de los LISA para todas las observaciones es proporcional al indicador global de asociación espacial.
Estos estadísticos permiten verificar si la contribución local a la agrupación alrededor de la observación “i” es significativamente diferente de cero, lo que permite identificar áreas donde existe un grado significativo de agrupación espacial.
Partiendo del índice global de Moran:
\[ I = \frac{n}{\sum_{i=1}^{n} \sum_{j=1}^{n} w_{ij}}\ \frac{\sum_{i=1}^{n} \sum_{j=1}^{n} (y_i - \bar{y})(y_j - \bar{y})}{\sum_{i=1}^{n} (y_i - \bar{y})^2} \]
Se reescribe como:
\[ I = n \left[ \sum_k \sum_j w_{kj} \right]^{-1}\left[ \sum_k (y_k-\bar{y})^2 \right]^{-1} \sum_i I_i \]
Donde:
\[ I_i = (y_i-\bar{y}) \sum_j w_{ij} (y_j - \bar{y}) \]
De manera que:
\(n \left[ \sum_k \sum_j w_{kj} \right]^{-1}\left[ \sum_k (y_k-\bar{y})^2 \right]^{-1}\) no depende de i, entonces para un conjunto dado de \(z_i\) esta expresión se comporta como una constante.
\[ I = cte \cdot \sum_i I_i \]
Reemplazando \(z6_i = y_i - \bar{y}\), se obtiene la expresión original propuesta por Anselin, 1995.
\[ I_i = z_i \sum_j w_{ij}z_j \]
\(I_i\) es el producto de \(z_i\) y la media de los valores \(z_j\) de los vecinos del polígono \(i\)
Si tanto el valor \(z_i\) como el valor medio de los vecinos \(z_j\) están por encima de la media, la cantidad será alta, lo que indica la presencia de un cluster de valores por encima de la media, centrado en el polígono “i”.
Si los signos de \(z_i\) y \(\sum_j w_{ij}z_j\) son diferentes y el valor de \(I_i\) es ampliamente negativo, esto podría sugerir la existencia de repulsión local (outliers).
Si la magnitud de \(I_i\) no es especialmente grande, esto indica que hay poca evidencia para concluir sobre la presencia de clusters o repulsión.
Finalmente, se puede aplicar un test de significancia estadística a cada valor de \(I_i\), utilizando una hipótesis nula que postula la ausencia de asociación espacial.
Cuadrantes de asociación espacial. Modificado de Siabato et al., 2018
¿Qué es un rezago espacial?
Antes de abordar la creación del diagrama de dispersión de Moran, se proporciona una breve explicación sobre los rezagos espaciales, los cuales son fundamentales para la construcción de dicho gráfico.
Una vez que se ha definido la estructura de pesos en la matriz \(W\), los rezagos espaciales se definen como el resultado de multiplicar la matriz de pesos por el valor de la variable de interés \(y\). En otras palabras, un rezago espacial es una suma de los valores de los vecinos, ponderados de acuerdo a su proximidad espacial.
\[ Wy \] Esta operación está bien definida ya que \(W \in \mathbb{R}^{n\times n}\) y \(y \in \mathbb{R}^{n \times 1}\).
En la práctica, suele hacerse como:
\[ [Wy]_i = \sum_{j=1}^n w_{ij}y_j \] Esto se puede expresar de manera equivalente a la operación matricial mencionada anteriormente. Sin embargo, en la práctica, al trabajar con muestras de gran tamaño, almacenar la matriz \(W\) como una “matriz pura” puede ser ineficiente en términos de uso de memoria. Por lo tanto, es común utilizar una estructura de datos conocida como “sparse matrix” para almacenar los valores de \(W\), ya que estos tienden a ser dispersos.
En el gráfico de dispersión de Moran, siguiendo una convención establecida, se coloca la variable de interés en el eje x y los valores de los rezagos espaciales en el eje y. El índice global de Moran se basa en una relación lineal entre estos valores y se representa gráficamente como la pendiente de la línea de tendencia.
(ggplot(moran_plot, aes(x, wx)) +
geom_point(aes(shape = is_inf, color = is_inf), size = 0.75) +
scale_color_manual(values = c("black", "red")) +
geom_hline(
yintercept = moran_plot$wx %>% mean(),
linetype = "dashed", alpha = 0.5
) +
geom_vline(
xintercept = moran_plot$x %>% mean(),
linetype = "dashed", alpha = 0.5
) +
geom_text(aes(x, wx, label = labels),
data = filter(moran_plot, is_inf),
size = 2.5
) +
geom_smooth(method = "lm", se = F, size = .6, col = "black") +
labs(
x = "Cases", y = "Spatial lagged cases",
title = "Moran Scatterplot",
shape = "Significativo", color = "Significativo"
) +
theme_bw() +
theme(plot.title = element_text(hjust = 0.5))) %>%
plotly::ggplotly()El gráfico se divide en cuatro cuadrantes, cada uno con su significado específico. En el cuadrante superior derecho se encuentran las áreas con una alta incidencia de casos de leucemia, rodeadas por áreas que también presentan niveles superiores al promedio de casos (high-high). Estas áreas se conocen como hotspots. Por otro lado, en el cuadrante inferior izquierdo se encuentran las áreas con una baja incidencia de casos de leucemia, rodeadas por áreas con una incidencia por debajo del promedio (low-low). Estas áreas se conocen como coldspots.
En la diagonal opuesta se encuentran los valores atípicos espaciales, lo que indica que estas áreas están rodeadas por vecinos que difieren significativamente de ellas.
Al igual que en el caso de la estadística global, se puede evaluar la divergencia de los estadísticos locales con respecto a los valores esperados, asumiendo la normalidad y aleatoriedad, utilizando métodos analíticos como el de punto de silla y exacto. Estos métodos permiten comprobar la significancia estadística de los resultados obtenidos.
local <- as.data.frame(localmoran(ny$Cases,
listw = nb2listw(NY_nb, style = "C")
))
local2 <- as.data.frame(localmoran.sad(lm(Cases ~ 1, ny),
nb = NY_nb, style = "C"
))
local3 <- as.data.frame(localmoran.exact(lm(Cases ~ 1, ny), nb = NY_nb, style = "C"))
ny %>%
select(Cases) %>%
mutate(
p_value = moran_plot$is_inf,
rezago = lag.listw(nb2listw(NY_nb, style = "B"), Cases),
cuadrante = case_when(
(Cases > mean(Cases)) & (rezago > mean(rezago)) &
(p_value == T) ~ "High-High",
(Cases <= mean(Cases)) & (rezago <= mean(rezago)) &
(p_value == T) ~ "Low-Low",
(Cases > mean(Cases)) & (rezago <= mean(rezago)) &
(p_value == T) ~ "High-Low",
(Cases <= mean(Cases)) & (rezago > mean(rezago)) &
(p_value == T) ~ "Low-High",
TRUE ~ "NS"
),
index = moran_plot$labels
) -> resultados
resultados_coords <- filter(resultados, p_value == T) %>%
mutate(
x = (st_centroid(resultados$geometry) %>% st_coordinates())[1],
y = (st_centroid(resultados$geometry) %>% st_coordinates())[2]
)
(ggplot(resultados) +
geom_sf(aes(fill = cuadrante), alpha = 0.8) +
scale_fill_manual(values = viridisLite::viridis(5)) +
labs(
title = "Lugares con influencia de Leucemia \n",
x = "Longitud", y = "Latitud"
) +
geom_sf_text(
data = resultados_coords,
mapping = aes(x, y, label = index),
color = "white"
) +
theme_bw() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank()
)) %>%
plotly::ggplotly()Es posible realizar una evaluación de riesgo constante al asumir una distribución Poisson con un valor de \(\lambda\) que representa la probabilidad de ocurrencia de la enfermedad. Para estandarizar la variable, se utiliza la propiedad de que tanto la media como la desviación estándar son iguales en esta distribución. Luego, se calculan los rezagos espaciales y se compara el índice de Moran actual con el índice de riesgo constante.
ny %>%
transmute(
p_afectada = sum(Cases) / sum(POP8) * POP8,
casos_std = (Cases - p_afectada) / sqrt(p_afectada),
rezago = stats::lag(nb2listw(NY_nb, style = "C"), casos_std),
indice_riesgo = casos_std * rezago,
indice = local$Ii
) -> riesgo
m_riesgo <- ggplot(riesgo) +
geom_sf(aes(fill = indice_riesgo)) +
labs(title = "Mapa de riesgo Constante") +
scale_fill_gradientn(
colours = viridisLite::viridis(20,
begin = 0.17,
end = 0.98
),
limits = c(-2.2, 3.7)
) +
theme_bw() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank()
)
m_actual <- ggplot(riesgo) +
geom_sf(aes(fill = indice)) +
labs(title = "Mapa de riesgo Standard") +
scale_fill_gradientn(
colours = viridisLite::viridis(20,
begin = 0.17,
end = 0.98
),
limits = c(-2.2, 3.7)
) +
theme_bw() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank()
)
cowplot::plot_grid(m_actual, m_riesgo, align = "hv", ncol = 2)
# 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))
monte_carlo_sim <- function(nsim, prop, Ii, nb) {
N <- length(prop)
sims <- matrix(0, ncol = nsim, nrow = N)
for (i in 1:nsim) {
y <- rpois(N, lambda = prop)
z <- (y - prop) / sqrt(prop)
rezago <- stats::lag(nb2listw(nb, style = "C"), z)
sims[, i] <- z * rezago
}
xrank <- apply(cbind(Ii, sims), 1, function(x) rank(x)[1])
diff <- nsim - xrank
diff <- ifelse(diff > 0, diff, 0)
pval <- punif((diff + 1) / (nsim + 1))
return(pval)
}
mc_sim <- riesgo %>%
select(p_afectada, indice_riesgo) %>%
mutate(
pvalue_cr = monte_carlo_sim(
999, p_afectada,
indice_riesgo, NY_nb
),
pvalue_sad = local2$`Pr. (Sad)`,
pvalue_exact = local3$`Pr. (exact)`,
pvalue_rand = local$`Pr(z != E(Ii))`,
pvalue_norm = local2$`Pr. (N)`
) %>%
pivot_longer(
cols = c(
"pvalue_cr", "pvalue_sad",
"pvalue_exact", "pvalue_rand", "pvalue_norm"
),
names_to = "type", values_to = "p_value"
) %>%
mutate(discrete = cut(p_value,
breaks = c(
0, 0.05, 0.1,
0.9, 0.95, 1, max(p_value)
),
include.lowest = T
))aux_labs <- c(
"pvalue_cr" = "Riesgo constante",
"pvalue_exact" = "Exacto",
"pvalue_norm" = "Normal",
"pvalue_rand" = "Aleatorio",
"pvalue_sad" = "Saddlepoint"
)
ggplot(mc_sim) +
geom_sf() +
geom_sf(aes(fill = discrete),
data = filter(mc_sim, p_value < 0.1 | p_value > 0.9 &
p_value < 1)
) +
facet_wrap(~type, labeller = labeller(type = aux_labs)) +
labs(title = "") +
scale_fill_manual(values = viridisLite::viridis(5), name = "P-value") +
theme_bw() +
theme(
plot.title = element_text(size = 14, face = "bold", hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank()
)
En el gráfico presentado, podemos realizar una verificación considerando
diferentes distribuciones con el fin de determinar si un área es inusual
o presenta agrupamiento. Esta información nos permite evaluar si los
índices observados tienen respaldo estadístico, lo cual refuerza su
validez y relevancia.
Ajuste de modelos para datos en área
Los modelos estadísticos nos ofrecen diversas formas de modelar la dependencia espacial. Aunque normalmente se asume independencia entre las observaciones o distribuciones, esto no se aplica a los datos espaciales. Sin embargo, podemos explicar la autocorrelación espacial entre áreas vecinas. En el caso de que el vector de variables respuesta siga una distribución normal multivariada, se puede expresar de la siguiente manera:
\[ Y=\mu + e, \] Una de las alternativas recomendadas es utilizar una transformación logarítmica que incluya los pesos, en nuestro caso la cantidad de población, para tener en cuenta la heterogeneidad del riesgo. Esta transformación logarítmica puede ser una opción viable para abordar este aspecto.
\[ Z_i=log \left( \frac{1000(Y_i+1)}{n_i}\right) \]
TCE <- readOGR("data", "TCE")## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\Usuario\Videos\trabajo3\data", layer: "TCE"
## with 11 features
## It has 5 fieldsny <- readOGR("data",'NY8_utm18')## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\Usuario\Videos\trabajo3\data", layer: "NY8_utm18"
## with 281 features
## It has 17 fieldsnylm <- lm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=ny)
summary(nylm)##
## Call:
## lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = ny)
##
## 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-13Para analizar los datos de un área determinada, se utiliza el centroide de cada polígono y se asignan diversas covariables con el objetivo de construir un modelo de regresión lineal que explique la variable Z. A continuación, se presenta un ejemplo que considera tres covariables en el modelo.
ny$lmresid <- residuals(nylm) # residuales modelo sin pesos
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 = ny)
## 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.001242320Si consideramos que la tendencia de la variable de interés ha sido eliminada y que las ponderaciones espaciales están estandarizadas por filas, el estimador de un paso del parámetro de regresión espacial puede proporcionar más información comparativa que el Índice de Moran. Para realizar una comparación, presentamos la estimación ajustada del parámetro de máxima verosimilitud utilizando los mismos datos:
NYlistwW <- nb2listw(NY_nb, style = "W")
aple(residuals(nylm), listw=NYlistwW)## [1] 0.2051536spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=ny, listw=NYlistwW)$lambda## lambda
## 0.2169261nysar <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=ny, listw=NYlistw)
summary(nysar)##
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = ny,
## 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.017201
##
## 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 los resultados obtenidos, se observa una correlación en los residuos, dado que el valor estimado de \(\lambda\) es 0.0405 y el valor p asociado es 0.022 en la prueba de máxima verosimilitud. Aunque el parámetro resulte significativo a un nivel de confianza del 95%, es recomendable realizar una verificación adicional con otro tipo de modelo para no descartar la posibilidad de una asociación. De esta manera, se evita basar las conclusiones únicamente en un único enfoque estadístico y se asegura una evaluación más completa y precisa.
Por otro lado, se observa que la covariable exposición no es significativa.
nylam1 <- c(nysar$lambda)
nylam2 <- c(LR1.Spautolm(nysar)$p.value)
ny$sar_trend <- nysar$fit$signal_trend
ny$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(ny, 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(ny, 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)
Es importante tener en cuenta que este modelo no considera la distribución heterogénea de la población, es decir, asume que la población es uniforme en todas las zonas. Por lo tanto, se incluye el argumento “weights = POP8” en el modelo LM.
nylmw <- lm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=ny, weights=POP8)
summary(nylmw)##
## Call:
## lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = ny,
## 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-13ny$lmwresid <- residuals(nylmw)Cuando se incorpora la ponderación de pesos, se observa que la variable covariable “exposición” se vuelve significativa, lo cual difiere del modelo anterior. Esto tiene sentido, ya que es plausible que exista una relación entre los casos de leucemia y el nivel de exposición.
lm.morantest(nylmw, NYlistw)##
## Global Moran I for regression residuals
##
## data:
## model: lm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = ny,
## 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.001245248Los hallazgos son intrigantes, ya que indican que la mala especificación identificada por el índice de Moran está más relacionada con la heteroscedasticidad que con la autocorrelación espacial.
nysarw <- spautolm(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME,
data = ny, listw = NYlistw, weights = POP8)
summary(nysarw)##
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = ny,
## 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.016522
##
## 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)Observando esto, se nota que no hay rastros de espacio.
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(ny, c("lmresid", "lmwresid"), sp.layout=list(TCEpts), col.regions=gry, col="transparent", lwd=0.5, at=seq(-2,4.5,0.5))
Se observa que el modelo con ponderación obtiene resudiales mas pequeños que con el modelo sin los pesos, esto dado en pequeños poligonos.
ny$sarw_trend <- nysarw$fit$signal_trend
ny$sarw_stochastic <- nysarw$fit$signal_stochastic
tr_pl <- spplot(ny, c("sarw_trend"), at=tr_at, col="transparent", col.regions=tr_rds, main=list(label="Trend", cex=0.8))
st_pl <- spplot(ny, 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)
Se observa que la variación espacial en la población entre tractos es responsable de la correlación espacial residual después del ajuste por covariables.
Modelos autorregresivos condicionales
La definición de CAR se basa en la distribución condicional de la distribución de los errores, es decir:
\[ e_i | e_{j \sim \: i } \sim N \left( \sum_{j \sim i } \frac{c_{ij}e_j }{\sum_{j\sim i }c_{ij} }, \frac{\sim^{2}_{e_i} }{\sum_{j\sim i }c_{ij} } \right ) \]
Donde \(e_{j\sim i }\) representa los vecinos del área i y \(c_{ij}\) son parámetros dependientes de \(b_{ij}\).
Se obtiene que:
\[ E[Y] =X^T\beta \\ Var(Y) =\sigma^2(1-\lambda W )^{-1}(1-\lambda W^T )^{-1} \]
nycar <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME , data=ny, family="CAR",listw=NYlistw)
summary(nycar)##
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = ny,
## 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.030872
##
## 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)Aunque la estimación de los parámetros es similar a la del modelo SAR, se observa que los valores de p para dos covariables, la distancia al TCE más cercano y el porcentaje de personas propietarias de una vivienda, superan el umbral de 0.05 con un nivel de confianza del 95%. Además, la prueba de verosimilitud revela la presencia de una autocorrelación espacial significativa. Es importante destacar que esta conclusión se realiza sin tener en cuenta los pesos ponderados de POP8.
nylam1 <- c(nycar$lambda)
nycarw <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=ny, family="CAR",listw=NYlistw, weights=POP8)
summary(nycarw)##
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = ny,
## 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.038977
##
## 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)Aplicando los pesos, se observa que las tres covariables son significativas ya que el p valor obtenido es muy pequeño.
Ajuste de modelos de regresión espacial
Mediante la función “spautolm”, se realiza el ajuste de modelos de regresión espacial utilizando el método de máxima verosimilitud. Este proceso implica encontrar el coeficiente autoregresivo que maximiza la función de verosimilitud para la familia de modelos seleccionada y luego ajustar los demás coeficientes mediante mínimos cuadrados generalizados en ese punto.
En otras palabras, el coeficiente autoregresivo espacial se obtiene mediante la búsqueda de líneas utilizando el método de optimización, en lugar de optimizar simultáneamente todos los parámetros del modelo. Se tiene lo siguiente:
\[ log(|I-\lambda W|)=log \left( \prod_{i=1}^n (1-\lambda \zeta_i) \right) \] Donde $ _i$ son los valores propios de \(W\).
nysarwM <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=ny,family="SAR",listw=NYlistw, weights=POP8, method="Matrix")
summary(nysarwM)##
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = ny,
## 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)nysar_ll <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=ny, family="SAR",listw=NYlistw, llprof=100)
nysarw_ll <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=ny, family="SAR",listw=NYlistw, weights=POP8, llprof=100)
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")
Se observa la función de log verosimilitud en función de \(\lambda\) y donde esta alcanza su máximo.
1/range(eigenw(NYlistw))## [1] -0.3029200 0.1549552nysmaw <- spautolm(Z~PEXPOSURE+PCTAGE65P+PCTOWNHOME, data=ny, family="SMA",listw=NYlistw, weights=POP8)
summary(nysmaw)##
## Call: spautolm(formula = Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME, data = ny,
## 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.01657
##
## Log likelihood: -251.6112
## ML residual variance (sigma squared): 1105.3, (sigma: 33.245)
## Number of observations: 281
## Number of parameters estimated: 6
## AIC: NA (not available for weighted model)Enfoques de econometría espacial
Un punto a favor acerca del análisis de datos espaciales es la amplia gama de disciplinas involucradas, lo cual permite dar multiples enfoques en esto cado para terminos ecónomicos, uno de ellos es la econometría espacial.La econometría espacial es descrita por Anselin (1988, 2002), avanzado por LeSage y Pace (2009) y con comentarios adicionales de Bivand (2002, 2006) con respecto a hacer econometría espacial en el sotfware R.
Cabe resaltar que en el caso de econometria, proponer diferentes supuestos en el área de econometria es algo complejo y no se cuenta con buenos conomientos en esa área.
Otros modelos
Se pueden utilizar otros metodos para explicar la correlación espacial como los GAM donde se asume una distribución, no necesariamente normal.
ny$x <- coordinates(ny)[, 1]/1000
ny$y <- coordinates(ny)[, 2]/1000
nyGAM1 <- gam(Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME +
s(x, y), weights = POP8, data = ny)
anova(nylmw, nyGAM1, test = "Chisq")## Analysis of Variance Table
##
## Model 1: Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME
## Model 2: Z ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME + s(x, y)
## Res.Df RSS Df Sum of Sq Pr(>Chi)
## 1 277.00 310778
## 2 273.19 305229 3.8096 5549.7 0.2671En este caso, no agrega valor a los modelos ya vistos anteriormente.
nyGLMp <- glm(Cases ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME +
offset(log(POP8)), data = NY8, family = "poisson")
summary(nyGLMp)##
## Call:
## glm(formula = Cases ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME + offset(log(POP8)),
## family = "poisson", data = NY8)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.13442 0.18261 -44.545 < 2e-16 ***
## PEXPOSURE 0.14893 0.03121 4.773 1.82e-06 ***
## PCTAGE65P 3.99816 0.59783 6.688 2.27e-11 ***
## PCTOWNHOME -0.35710 0.19027 -1.877 0.0605 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 428.25 on 280 degrees of freedom
## Residual deviance: 353.35 on 277 degrees of freedom
## AIC: Inf
##
## Number of Fisher Scoring iterations: 5Se observa que hay poca sobredispersión, además, el ajuste con la familia poisson en el que el parámetro de dispersión no se fija en la unidad, por lo que pueden modelar una dispersión excesiva que no da como resultado grandes cambios.
nyGAMp <- gam(Cases ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME + offset(log(POP8)) + s(x, y), data = ny, family = "poisson")
summary(nyGAMp)##
## Family: poisson
## Link function: log
##
## Formula:
## Cases ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME + offset(log(POP8)) +
## s(x, y)
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -8.13659 0.21475 -37.889 < 2e-16 ***
## PEXPOSURE 0.16807 0.05985 2.808 0.00498 **
## PCTAGE65P 3.71993 0.64308 5.785 7.27e-09 ***
## PCTOWNHOME -0.36019 0.19940 -1.806 0.07087 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(x,y) 7.708 10.71 9.989 0.504
##
## R-sq.(adj) = 0.394 Deviance explained = 21.4%
## UBRE = 0.2815 Scale est. = 1 n = 281anova(nyGLMp, nyGAMp, test = "Chisq")## Analysis of Deviance Table
##
## Model 1: Cases ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME + offset(log(POP8))
## Model 2: Cases ~ PEXPOSURE + PCTAGE65P + PCTOWNHOME + offset(log(POP8)) +
## s(x, y)
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 277.00 353.35
## 2 269.29 336.68 7.7079 16.666 0.02907 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1