Modelo de regresión espacial para la Tasa de Homicidios en USA en 1990
Parcial 3 Teoria de Estadistica Espacial
Camila Alejandra Pinzon Caro
2025-12-05
Planteamiento del problema
La base de datos US Homicides + Socio-Economics (1960-90) es importante por ser un recurso de referencia que permite una comprensión profunda de las raíces complejas de la violencia gracias a su total cobertura de datos en todos los condados de Estados Unidos durante 4 décadas. La relevancia de estudiar las tasas de homicidios junto con variables socioeconómicas radica en la necesidad de ir más allá de la simple estadística criminal para entender el contexto social, económico e histórico que moldea la violencia.
La tasa de homicidios por cada 100.000 habitantes es la variable central de esta base, y el presente estudio va a buscar la forma de explicar su comportamiento espacial mediante modelos de regresión que involucren a las variables socioeconómicas disponibles.
Descripción de los datos
Esta base de datos de 3.085 observaciones correspondientes a condados en Estados Unidos medidos en el año 1990 cuenta con 69 variables, de las cuales 13 de ellas se usarán como variables explicativas en el modelo de regresión de la Tasa de homicidios. La descripción de las variables consideradas se muestra a continuación:
| Variable | Descripción |
|---|---|
| FIPS | Código fips combinado para el estado y el condado |
| HR | Tasa de homicidios por cada 100.000 habitantes |
| PO | Población del condado |
| RD | Escasez de recursos (componente principal) |
| PS | Estructura de la población (componente principal) |
| UE | Tasa de desempleo |
| DV | Tasa de divorcios |
| MA | Mediana de la edad |
| POL | Logaritmo de la población |
| DNL | Logaritmo de la densidad de población |
| MFIL | Logaritmo de la mediana de los ingresos por familia |
| FP | Porcentaje de familias en condición de pobreza |
| BLK | Porcentaje de población negra |
| GI | Indice de Gini para la desigualdad en ingresos de las familias |
| FH | Porcentaje de hogares encabezados por mujeres |
En el siguiente gráfico se muestra la variable HR90. Los colores de los condados están dados por los cuantiles de variable, donde los colores más oscuros representan una mayor tasa de homicidios (11.52-71.38). Note que la imagen nos deja ver como hacia el sur este del país, la tasa de homicidios es alta y homogénea en el espacio.
Carga de librerías
library(sf)## Warning: package 'sf' was built under R version 4.5.2## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUElibrary(ggplot2)
library(ggspatial)## Warning: package 'ggspatial' was built under R version 4.5.2library(spdep)## Warning: package 'spdep' was built under R version 4.5.2## Cargando paquete requerido: spData## Warning: package 'spData' was built under R version 4.5.2## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(purrr)
library(patchwork)## Warning: package 'patchwork' was built under R version 4.5.2library(deldir)## Warning: package 'deldir' was built under R version 4.5.2## deldir 2.0-4 Nickname: "Idol Comparison"##
## The syntax of deldir() has changed since version
## 0.0-10. Read the help!!!.library(spatialreg)## Warning: package 'spatialreg' was built under R version 4.5.2## Cargando paquete requerido: Matrix##
## Adjuntando el paquete: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption#library(sppois)
library(devtools)## Warning: package 'devtools' was built under R version 4.5.2## Cargando paquete requerido: usethis## Warning: package 'usethis' was built under R version 4.5.2library(pscl)## Warning: package 'pscl' was built under R version 4.5.2## Classes and Methods for R originally developed in the
## Political Science Computational Laboratory
## Department of Political Science
## Stanford University (2002-2015),
## by and under the direction of Simon Jackman.
## hurdle and zeroinfl functions by Achim Zeileis.library(MASS)##
## Adjuntando el paquete: 'MASS'## The following object is masked from 'package:patchwork':
##
## arealibrary(CARBayes)## Warning: package 'CARBayes' was built under R version 4.5.2## Cargando paquete requerido: RcppCarga de datos y arreglo de coordenadas
US_homi <- st_read("ncovr/ncovr/NAT.shp")## Reading layer `NAT' from data source
## `C:\Users\CAMILA\OneDrive - Universidad Nacional de Colombia\Escritorio\Universidad\Maestría Estadística\2do semestre\Espacial\ncovr\ncovr\NAT.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 3085 features and 69 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -124.7314 ymin: 24.95597 xmax: -66.96985 ymax: 49.37173
## Geodetic CRS: WGS 84US_homi_reg <- US_homi[,c("HR90", "PO90", "RD90", "PS90", "UE90","DV90", "FH90",
"MA90", "POL90", "DNL90", "MFIL89", "FP89", "BLK90", "GI89")]
US_homi_utm <- st_transform(US_homi_reg, 5070)
#st_crs(US_homi_utm)Mapa de valores observados
ggplot(US_homi_reg) +
geom_sf(aes(fill = HR90)) +
scale_fill_viridis_c(direction = -1, option='B') + # Escala de colores para la variable
labs(fill = "Tasa de homicidio") +
theme_minimal() +
ggtitle("Tasas de homicidio por condados en EEUU, año 1990") +
# Añadir la leyenda personalizada
theme(legend.position = "right") +
# Añadir escala
annotation_scale(location = "bl", width_hint = 0.15) +
annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"))## Scale on map varies by more than 10%, scale bar may be inaccurate
Este mapa nos deja ver como se comporta la variable respuesta en las diferentes áreas. En general, la tasa de homicidios en los condados se mantiene en valores bajos cercanos a 0, son pocos los condados con tasas de homicidio considerablemente altas.
También se presentan los mapas de las variables explicativas consideradas a continuación:
# Variables a graficar
vars <- c("PO90", "RD90", "PS90", "UE90", "DV90", "FH90", "MA90", "POL90", "DNL90", "MFIL89",
"FP89", "BLK90", "GI89")
names <- c("Población", "Escasez de recursos", "Estructura de la población", "Tasa de desempleo", "Tasa de divorcios ", "Edad mediana", "Logaritmo de la población", "Logaritmo de la densidad de población", "Logaritmo de la mediana de ingresos por familia", "Porcentaje de familias en condición de pobreza ", "Porcentaje de población negra", "Indice de Gini", "Porcentaje de hogares encabezados por mujeres")
# Crear lista de gráficos
plots <- map2(vars, names, function(v, name_label) {
ggplot(US_homi_reg) +
geom_sf(aes(fill = .data[[v]])) +
scale_fill_viridis_c(direction = -1, option='B') +
ggtitle(name_label) +
theme_minimal() +
theme(legend.position = "right")
})
ggsave("covariables_maps.png",
wrap_plots(plots, ncol = 2),
width = 10, height = 20, dpi = 300)
BoxMap
A continuación, y como medida exploratoria gráfica se presenta el BoxMap de la Tasa de Homicidios generado en la herramienta GeoDa
El BoxMap nos deja ver que las áreas consideradas atípicas son aquellas con una tasa de homicidios mayor a 20.34, resultando en 129 áreas en total. Así mismo se observa que estas áreas se concentran mayormente en la parte sur este del país. Por otro lado, en la parte central norte del país se presentan las menores tasas de homicidios.
Finalmente, a excepción de unos cuantos casos aislados, parece ser que la tasa de homicidios presenta una dependencia de tipo espacial, ya que condados “seguros” tienden a estar rodeados de otros condados con las mismas características.
Correlación con las variables explicativas
cor(as.data.frame(US_homi_reg)[,1:14])[,1]## HR90 PO90 RD90 PS90 UE90 DV90 FH90
## 1.0000000 0.1537506 0.5678323 0.1626505 0.3063221 0.1813698 0.6284786
## MA90 POL90 DNL90 MFIL89 FP89 BLK90 GI89
## -0.2013303 0.1430758 0.1706063 -0.2243756 0.4289281 0.5929101 0.4708659Las correlaciones nos dejan ver que la variable Porcentaje de hogares encabezados por mujeres es la más correlacionada linealmente con la variable respuesta (62.84%), seguida de el Porcentaje de población negra (59.29%) y el componente principal de Escasez de recursos (56.78%).
Matriz de vecindades
Como estamos trabajando con variables que no son continuas en el espacio para las cuales los individuos ahora son áreas, es de interés conocer las estructuras de conectividad espacial entre estas unidades. Esto se logra mediante las matrices de vecindad. Para el presente estudio, cada condado de EEUU tiene una “vecindad” que se define por su proximidad con otros condados.
Para identificar estas subregiones se establece una coordenada representativa para cada poligono, conocida como centroide.
# Centroides
centroides <- st_centroid(US_homi_utm)## Warning: st_centroid assumes attributes are constant over geometries# Matriz de distancias entre los centroides
Wdist <- st_distance(centroides)La cuantificación de estas relaciones de acuerdo a algún criterio se resumen en una matriz cuadrada \(\mathbf{W}_{N\times N}\) llamada Matriz de pesos espacial. Entradas distintas de cero indican unidades relacionadas. Como esta matriz no es conocida debe estimarse, y hay varias formas de llevar esto a cabo: Matrices de peso binarias que asignan 1 a las unidades contiguas y 0 e.o.c., medidas gráficas como la construcción de grafos y matrices de distancias.
En el presente trabajo se considerarán los siguientes criterios para estimar la matriz \(\mathbf{W}\):
- Criterio Torre: Dos áreas son vecinas si comparten un borde común.
- Criterio Reina: Dos áreas son vecinas si comparten un borde o un vértice común.
- Triangulación de Delaunay: Considera como vecinas de un área \(A_i\) a aquellas que pueden ser conectadas mediante triángulos, t.q. el centroide \(\mathbf{s}_i\) sea al menos un vértice.
- Gráfica de Gabriel: Dos áreas son vecinas si \[d(\mathbf{s}_i,\mathbf{s}_j)\leq \min\left\{\sqrt{d(\mathbf{s}_i,\mathbf{s}_k)^2+d(\mathbf{s}_j,\mathbf{s}_k)^2}| \forall \mathbf{s}_k\in\mathbf{S}\right\}\]
- Vecinos relativos: Dos áreas son vecinas si \[d(\mathbf{s}_i,\mathbf{s}_j)\leq \min \left\{\max d(\mathbf{s}_i,\mathbf{s}_j), d(y,\mathbf{s}_k) \right\}^{1/2}\]
- k-nn: Estas relaciones se establecen con distancias no simétricas. El criterio de vecindario elige, para cada \(\mathbf{s}_i\), las k unidades más cercanas. Así, el resto de los datos \(k+1, k+2, ...\) son considerados no vecinos.
Ahora bien, dentro de cada uno de estos criterios se pueden considerar ponderaciones para las matrices de vecindad, reflejando distintos grados de proximidad: Ponderaciones binarias, por contigüidad, por pesos y sin pesos.
A continuación se definen todas las matrices de pesos a evaluar
## Criterio Torre ##
w_torre_b <- nb2listw(poly2nb(US_homi_utm, queen=FALSE), style = "B", zero.policy = TRUE)
w_torre_c <- nb2listw(poly2nb(US_homi_utm, queen=FALSE), style = "C", zero.policy = TRUE)
w_torre_w <- nb2listw(poly2nb(US_homi_utm, queen=FALSE), style = "W", zero.policy = TRUE)
w_torre_u <- nb2listw(poly2nb(US_homi_utm, queen=FALSE), style = "U", zero.policy = TRUE)
## Criterio Reina ##
w_reina_b <- nb2listw(poly2nb(US_homi_utm, queen=TRUE), style = "B", zero.policy = TRUE)
w_reina_c <- nb2listw(poly2nb(US_homi_utm, queen=TRUE), style = "C", zero.policy = TRUE)
w_reina_w <- nb2listw(poly2nb(US_homi_utm, queen=TRUE), style = "W", zero.policy = TRUE)
w_reina_u <- nb2listw(poly2nb(US_homi_utm, queen=TRUE), style = "U", zero.policy = TRUE)
## Triangulación de Delaunay ##
xy <- st_coordinates(centroides)
#del <- deldir(xy[,1], xy[,2]) # triangulación que no me daña el R :b
#w_trian_b <- nb2listw(spdep::tri2nb(del$delsgs), style = "B", zero.policy = TRUE)
#w_trian_c <- nb2listw(tri2nb(del$delsgs), style = "C", zero.policy = TRUE)
#w_trian_w <- nb2listw(tri2nb(del$delsgs), style = "W", zero.policy = TRUE)
#w_trian_u <- nb2listw(tri2nb(del$delsgs), style = "U", zero.policy = TRUE)
## Gráfica de Gabriel ##
w_gabriel_b <- nb2listw(graph2nb(gabrielneigh(xy), sym = TRUE), style = "B", zero.policy = TRUE)
w_gabriel_c <- nb2listw(graph2nb(gabrielneigh(xy), sym = TRUE), style = "C", zero.policy = TRUE)
w_gabriel_w <- nb2listw(graph2nb(gabrielneigh(xy), sym = TRUE), style = "W", zero.policy = TRUE)
w_gabriel_u <- nb2listw(graph2nb(gabrielneigh(xy), sym = TRUE), style = "U", zero.policy = TRUE)
## Vecinos relativos ##
w_relativ_b <- nb2listw(graph2nb(relativeneigh(xy), sym = TRUE), style = "B", zero.policy = TRUE)
w_relativ_c <- nb2listw(graph2nb(relativeneigh(xy), sym = TRUE), style = "C", zero.policy = TRUE)
w_relativ_w <- nb2listw(graph2nb(relativeneigh(xy), sym = TRUE), style = "W", zero.policy = TRUE)
w_relativ_u <- nb2listw(graph2nb(relativeneigh(xy), sym = TRUE), style = "U", zero.policy = TRUE)
## K-nn ##
w_knn1_b <- nb2listw(knn2nb(knearneigh(xy, k = 1)), style = "B", zero.policy = TRUE)## Warning in knn2nb(knearneigh(xy, k = 1)): neighbour object has 824 sub-graphsw_knn1_c <- nb2listw(knn2nb(knearneigh(xy, k = 1)), style = "C", zero.policy = TRUE)## Warning in knn2nb(knearneigh(xy, k = 1)): neighbour object has 824 sub-graphsw_knn1_w <- nb2listw(knn2nb(knearneigh(xy, k = 1)), style = "W", zero.policy = TRUE)## Warning in knn2nb(knearneigh(xy, k = 1)): neighbour object has 824 sub-graphsw_knn1_u <- nb2listw(knn2nb(knearneigh(xy, k = 1)), style = "U", zero.policy = TRUE)## Warning in knn2nb(knearneigh(xy, k = 1)): neighbour object has 824 sub-graphsw_knn2_b <- nb2listw(knn2nb(knearneigh(xy, k = 2)), style = "B", zero.policy = TRUE)
w_knn2_c <- nb2listw(knn2nb(knearneigh(xy, k = 2)), style = "C", zero.policy = TRUE)
w_knn2_w <- nb2listw(knn2nb(knearneigh(xy, k = 2)), style = "W", zero.policy = TRUE)
w_knn2_u <- nb2listw(knn2nb(knearneigh(xy, k = 2)), style = "U", zero.policy = TRUE)
w_knn3_b <- nb2listw(knn2nb(knearneigh(xy, k = 3)), style = "B", zero.policy = TRUE)
w_knn3_c <- nb2listw(knn2nb(knearneigh(xy, k = 3)), style = "C", zero.policy = TRUE)
w_knn3_w <- nb2listw(knn2nb(knearneigh(xy, k = 3)), style = "W", zero.policy = TRUE)
w_knn3_u <- nb2listw(knn2nb(knearneigh(xy, k = 3)), style = "U", zero.policy = TRUE)
w_knn4_b <- nb2listw(knn2nb(knearneigh(xy, k = 4)), style = "B", zero.policy = TRUE)
w_knn4_c <- nb2listw(knn2nb(knearneigh(xy, k = 4)), style = "C", zero.policy = TRUE)
w_knn4_w <- nb2listw(knn2nb(knearneigh(xy, k = 4)), style = "W", zero.policy = TRUE)
w_knn4_u <- nb2listw(knn2nb(knearneigh(xy, k = 4)), style = "U", zero.policy = TRUE)Indice de Moran
Este índice mide la correlación espacial entre cada área y sus vecinos, extendiendo la idea del coeficiente de correlación de Pearson al caso espacial. Evalúa si los valores de una variable en una región tienden a estar agrupados en áreas vecinas o, por el contrario, dispersos. Se calcula comparando el valor de la variable en cada unidad espacial con los valores en sus vecinas, ponderados por la distancia o vecindad.
Aunque este índice por si solo no nos dice mucho, su valor p si lo hace. Es por esto que este valor es el que se compara para elegir cual matriz de pesos se ajusta mejor a los datos
matrices <- list(w_torre_b, w_torre_c, w_torre_w, w_torre_u,
w_reina_b, w_reina_c, w_reina_w, w_reina_u,
#w_trian_b, w_trian_c, w_trian_w, w_trian_u,
w_gabriel_b, w_gabriel_c, w_gabriel_w, w_gabriel_u,
w_relativ_b, w_relativ_c, w_relativ_w, w_relativ_u,
w_knn1_b, w_knn1_c, w_knn1_w, w_knn1_u,
w_knn2_b, w_knn2_c, w_knn2_w, w_knn2_u,
w_knn3_b, w_knn3_c, w_knn3_w, w_knn3_u,
w_knn4_b, w_knn4_c, w_knn4_w, w_knn4_u)
aux <- numeric(0)
for(i in 1:length(matrices)){
aux[i] <- moran.test(US_homi_reg$HR90, matrices[[i]], alternative = "two.sided")$p.value
}
which.min(aux) # Criterio Reina con ponderación C## [1] 6matrices[[which.min(aux)]]## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 3085
## Number of nonzero links: 18168
## Percentage nonzero weights: 0.190896
## Average number of links: 5.889141
##
## Weights style: C
## Weights constants summary:
## n nn S0 S1 S2
## C 3085 9517225 3085 1047.691 13043.78moran.test(US_homi_reg$HR90, matrices[[which.min(aux)]], alternative="two.sided")##
## Moran I test under randomisation
##
## data: US_homi_reg$HR90
## weights: matrices[[which.min(aux)]]
##
## Moran I statistic standard deviate = 37.101, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.3879660699 -0.0003242542 0.0001095317Obtenemos que la matriz de vecindad que minimiza el valor p del Índice de Moran es la generada por el Criterio Reina con ponderación por Contigüidad.
La estadística obtenida en el Test \[H_0:\text{Ausencia de correlación espacial (IM = 0)}\\H_1:\text{Presencia de correlación espacial}(\text{IM}\neq0)\] nos sugiere una autocorrelación espacial positiva, mostrando que los datos de HR90 tienden a agruparse. Como el valor p es 0, se rechaza la hipótesis nula, afirmando que hay autocorrelación espacial significativa
moran.plot(US_homi_reg$HR90,
w_reina_c,
xlab="Tasas de homicidio",
ylab="Tasas de homicidio rezagadas",
las=1,
pch=16,
cex=0.5)
legend("bottomright",
legend=c("I de Moran: 0.387", "Valor P: 0"),
cex=1,
bg='lightblue')
title("Dispersograma de Moran para las tasas de homicidio por condado en\n Estados Unidos, año 1990", cex.main=1)
Simulación Monte carlo para el Índice de Moran
A continuación se realiza una simulación de Monte Carlo para obtener la distribución empírica del Índice de Moran mediante permutaciones de la variable respuesta.
moran_obs <- moran.test(US_homi_reg$HR90, w_reina_c, alternative = "two.sided")
# Simulacion montecarlo
set.seed(123)
B <- 999 # número de permutaciones
moran_sim <- replicate(999,
moran(sample(US_homi_reg$HR90), w_reina_c, n = length(US_homi_reg$HR90), S0 = Szero(w_reina_c))$I
)
I_obs <- moran_obs$estimate[["Moran I statistic"]]
# valor p empírico
p_mc <- (sum(abs(moran_sim) >= abs(I_obs)) + 1) / (B + 1)
p_mc## [1] 0.001hist(moran_sim, breaks = 30,
main = "Distribución Monte Carlo del Índice de Moran",
xlab = "IM simulado")
abline(v = I_obs, col = "red", lwd = 3)
legend("topright",
legend=c("I de Moran: 0.387", "P-valor empírico: 0.001"),
cex=1,
bg='lightblue')
En el test de Moran, bajo la hipótesis nula de ausencia de autocorrelación la distribución del estadístico está centrada alrededor de 0 y con un rango pequeño. Este es el comportamiento que vemos en la gráfica porque las permutaciones rompen con la estructura espacial y las diferencias con los vecinos son netamente aleatorias. El Indice de Moran observado (0.387) está muy lejos de la distribución nula, mostrando que hay una autocorrelación espacial significativa.
Esto se traduce en que condados con tasas de homicidio altas tienden a estar rodeados por condados con tasas altas también, igualmente para los que presentan tasas bajas.
Local G
El Local G es un índice de autocorrelación espacial que complementa el índice de Moran al permitir la identificación de patrones locales de agrupamiento o dispersión en los datos. Mientras que el índice de Moran proporciona una medida global de la autocorrelación espacial, el Local G permite detectar si ciertas áreas específicas presentan una concentración significativa de valores similares (o diferentes)
localG <- localG(US_homi_reg$HR90, w_reina_c)Para evaluar la significancia de los resultados obtenidos con el índice Local G se realizan 1000 simulaciones MonteCarlo como se muestra a continuación:
# Simulacion montecarlo
sim.G <- matrix(0,1000,3085)
for(i in 1:1000) sim.G[i,] <- localG(sample(US_homi_reg$HR90), w_reina_c)
mc.pvalor.G <- (colSums(sweep(sim.G,2,localG,">="))+1)/(nrow(sim.G)+1)
#mc.pvalor.G
table(mc.pvalor.G<0.05)##
## FALSE TRUE
## 2599 486Para 486 condados se sugiere que su autocorrelación espacial observada es significativamente diferente de la aleatoriedad, indicando un cluster local
localG <- as.numeric(localG)
# Mapa de G
ggplot(US_homi_utm) +
geom_sf(aes(fill = localG)) +
scale_fill_viridis_c() +
ggtitle("G Getis Ord Local para las tasas de homicidio en EEUU, año 1990") +
theme_minimal()
# Mapa de p-valor
ggplot(US_homi_utm) +
geom_sf(aes(fill = mc.pvalor.G)) +
scale_fill_viridis_c() +
ggtitle("P-Valor de G Getis Ord Local para las tasas de homicidio en EEUU, año 1990") +
theme_minimal()
Los condados con p-valores bajos son los ubicados en el sur este del país. Como se ha visto en análisis previos, este resultado corrobora que esta es una zona con un patrón espacial significativo ya que en esta se agrupan las tasas de homicidios más altas. Podrían considerarse como un cluster de condados donde el riesgo de homicidio es más elevado que en el resto de condados del país. Esta zona es la compuesta por los estados de Alabama, Florida, Georgia, Kentucky, Misisipi, Carolina del Norte, Carolina del Sur, Tennessee, Virginia y Virginia Occidental.
Matrices de vecindades para las covariables
Como se verá más adelante, los modelos de regresión sugieren que las variables explicativas están guardando también una estructura espacial que debe explicarse mediante matrices de pesos. Es por esto que, siguiendo el procedimiento que se hizo para la variable respuesta, se procede a elegir cual matriz de pesos es mejro para dada covariable o ver si hay alguna que no requiera ser rezagada espacialmente.
covs <- c("RD90","PS90","UE90","DV90","FH90","MA90",
"DNL90","MFIL89","FP89","BLK90","GI89")
matrices_nombres <- c("w_torre_b", "w_torre_c", "w_torre_w", "w_torre_u",
"w_reina_b", "w_reina_c", "w_reina_w", "w_reina_u",
#w_trian_b, w_trian_c, w_trian_w, w_trian_u,
"w_gabriel_b", "w_gabriel_c", "w_gabriel_w", "w_gabriel_u",
"w_relativ_b", "w_relativ_c", "w_relativ_w", "w_relativ_u",
"w_knn1_b", "w_knn1_c", "w_knn1_w", "w_knn1_u",
"w_knn2_b", "w_knn2_c", "w_knn2_w", "w_knn2_u",
"w_knn3_b", "w_knn3_c", "w_knn3_w", "w_knn3_u",
"w_knn4_b", "w_knn4_c", "w_knn4_w", "w_knn4_u")
names(matrices) <- matrices_nombres
resultados <- list()
for (var in covs) {
aux <- numeric(length(matrices)) # p-values
names(aux) <- names(matrices) # nombres de matrices
# Calcular p-value para cada matriz
for (i in seq_along(matrices)) {
aux[i] <- moran.test(
US_homi_reg[[var]],
matrices[[i]],
alternative = "two.sided"
)$p.value
}
# Identificar la matriz con menor p-value
mejor_matriz_nombre <- names(aux)[which.min(aux)]
# Guardar todo
resultados[[var]] <- list(
pvalues = aux,
mejor_nombre = mejor_matriz_nombre,
mejor_matriz = matrices[[mejor_matriz_nombre]],
moran_test = moran.test(
US_homi_reg[[var]],
matrices[[mejor_matriz_nombre]],
alternative = "two.sided"
)
)
}
# 6. Mostrar el resultado legible
for (var in covs) {
cat("\n=======================\n")
cat("Variable:", var, "\n")
cat("Mejor matriz:", resultados[[var]]$mejor_nombre, "\n")
print(resultados[[var]]$moran_test)
}##
## =======================
## Variable: RD90
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 64.302, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.6922006545 -0.0003242542 0.0001159912
##
##
## =======================
## Variable: PS90
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 54.107, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.5825880096 -0.0003242542 0.0001160624
##
##
## =======================
## Variable: UE90
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 51.341, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.5524668836 -0.0003242542 0.0001159293
##
##
## =======================
## Variable: DV90
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 48.373, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.5206809840 -0.0003242542 0.0001160053
##
##
## =======================
## Variable: FH90
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 53.568, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.5765299256 -0.0003242542 0.0001159630
##
##
## =======================
## Variable: MA90
## Mejor matriz: w_reina_c
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 37.039, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.3877547975 -0.0003242542 0.0001097776
##
##
## =======================
## Variable: DNL90
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 60.38, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.6501515108 -0.0003242542 0.0001160597
##
##
## =======================
## Variable: MFIL89
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 61.082, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.6576945297 -0.0003242542 0.0001160505
##
##
## =======================
## Variable: FP89
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 61.657, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.6635878411 -0.0003242542 0.0001159455
##
##
## =======================
## Variable: BLK90
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 76.411, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.8223373207 -0.0003242542 0.0001159139
##
##
## =======================
## Variable: GI89
## Mejor matriz: w_torre_b
##
## Moran I test under randomisation
##
## data: US_homi_reg[[var]]
## weights: matrices[[mejor_matriz_nombre]]
##
## Moran I statistic standard deviate = 59.898, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.6449520292 -0.0003242542 0.0001160570Note que, TODAS las covariables rechazan la hipótesis nula de ausencia de autocorrelación espacial. Así mismo, revisando los valores p obtenidos en los Test del Indice de Moran para cada variable se observa que, sin importar la matriz de vecindad que se establezca, este p-valor es igual a 0. Es por esto, y por sencillez del modelo, que todas estas covariables se considerarán rezagadas con la misma matriz de pesos que la seleccionada para la Tasa de homicidios (Criterio Reina con ponderación por contigüidad)
Modelos de regresión
Modelos de regresión Poisson
La variable de interés en el estudio es la Tasa de homicidios. Para poder modelarla de forma más sencilla mediante modelos de regresión de la familia Poisson, transformamos esta variable en Conteo de homicidios por condado.
df_US_homi <- as.data.frame(US_homi_reg)[,1:14]
df_US_homi$conteos_HR <- round(US_homi_reg$HR90*US_homi_reg$PO90/100000)
ecuacion = conteos_HR ~ RD90+PS90+UE90+DV90+FH90+DNL90+MFIL89+FP89+BLK90+GI89+MA90 + offset(log(PO90))Modelo Poisson OLS
Este es el modelo lineal generalizado más sencillo que se puede construir sin asumir estructura espacial en los datos
\[Y_i\sim\text{Poisson}(\exp(\mathbf{X}\beta+\epsilon))\]
glm_pois <- glm(ecuacion, data = df_US_homi, family = "poisson")
summary(glm_pois)##
## Call:
## glm(formula = ecuacion, family = "poisson", data = df_US_homi)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.994e+08 9.984e+07 -4.001 6.32e-05 ***
## RD90 6.713e+07 1.678e+07 4.001 6.32e-05 ***
## PS90 5.946e-01 2.575e-02 23.094 < 2e-16 ***
## UE90 -1.167e-02 5.524e-03 -2.112 0.0347 *
## DV90 9.183e-02 4.649e-03 19.754 < 2e-16 ***
## FH90 -2.844e+06 7.109e+05 -4.001 6.32e-05 ***
## DNL90 -1.765e-01 1.650e-02 -10.699 < 2e-16 ***
## MFIL89 6.642e+07 1.660e+07 4.001 6.32e-05 ***
## FP89 -2.732e+06 6.828e+05 -4.001 6.32e-05 ***
## BLK90 -9.456e+05 2.364e+05 -4.001 6.32e-05 ***
## GI89 -5.324e+08 1.331e+08 -4.001 6.32e-05 ***
## MA90 -3.661e-03 3.048e-03 -1.201 0.2297
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 17535.7 on 3084 degrees of freedom
## Residual deviance: 3736.5 on 3073 degrees of freedom
## AIC: 9958.7
##
## Number of Fisher Scoring iterations: 5res_poi <- residuals(glm_pois, type = "pearson")
moran.test(res_poi, w_reina_c, alternative="two.sided")##
## Moran I test under randomisation
##
## data: res_poi
## weights: w_reina_c
##
## Moran I statistic standard deviate = 22.09, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2310568874 -0.0003242542 0.0001097134El modelo de regresión glm evidentemente presenta autocorrelación espacial en sus residuales.
Modelo Poisson-Moran Eigenvectors
Este es un modelo que captura autocorerlación espacial en regresores no observados mediante los vectores propios de la Matriz de Moran. \[Y_i\sim\text{Poisson}(\mu_i)\\ \log(\mu_i)=\mathbf{X}_i\beta+\sum_{k=1}^K\alpha_kM_{ik}\]
donde \(M_{ik}\) es el k-ésimo eigenvector significativo del operador espacial.
La función ME() busca remover la correlación espacial de los residuales de un modelo lineal generalizado. El calculo de los vectores se lleva a cabo a continuación:
me_pois <- ME(ecuacion, data = df_US_homi, listw = w_reina_c, family = "poisson", verbose = FALSE, alpha = 0.2)df_US_homi <- cbind(df_US_homi, fitted(me_pois))
ecuacion2 = conteos_HR ~ RD90+PS90+UE90+DV90+FH90+DNL90+MFIL89+FP89+BLK90+GI89 + offset(log(PO90)) + fitted(me_pois)
glm_ME <- glm(ecuacion2, data = df_US_homi, family = "poisson")
summary(glm_ME)##
## Call:
## glm(formula = ecuacion2, family = "poisson", data = df_US_homi)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.732e+08 1.023e+08 -2.671 0.00756 **
## RD90 4.592e+07 1.719e+07 2.671 0.00756 **
## PS90 4.644e-01 2.829e-02 16.416 < 2e-16 ***
## UE90 -8.572e-03 5.534e-03 -1.549 0.12140
## DV90 9.230e-02 4.785e-03 19.291 < 2e-16 ***
## FH90 -1.946e+06 7.284e+05 -2.671 0.00756 **
## DNL90 -1.055e-01 1.752e-02 -6.023 1.71e-09 ***
## MFIL89 4.544e+07 1.701e+07 2.671 0.00756 **
## FP89 -1.869e+06 6.996e+05 -2.671 0.00756 **
## BLK90 -6.469e+05 2.422e+05 -2.671 0.00756 **
## GI89 -3.642e+08 1.364e+08 -2.671 0.00756 **
## fitted(me_pois)vec71 2.260e+00 3.642e-01 6.206 5.43e-10 ***
## fitted(me_pois)vec165 -4.274e+00 3.930e-01 -10.875 < 2e-16 ***
## fitted(me_pois)vec2951 -1.736e+00 3.479e-01 -4.992 5.98e-07 ***
## fitted(me_pois)vec967 1.709e+00 3.571e-01 4.785 1.71e-06 ***
## fitted(me_pois)vec2096 1.940e+00 4.069e-01 4.767 1.87e-06 ***
## fitted(me_pois)vec267 2.323e+00 4.535e-01 5.123 3.01e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 17535.7 on 3084 degrees of freedom
## Residual deviance: 3487.5 on 3068 degrees of freedom
## AIC: 9719.6
##
## Number of Fisher Scoring iterations: 6res_ME <- residuals(glm_ME, type = "pearson")
moran.test(res_ME, w_reina_c, alternative="two.sided")##
## Moran I test under randomisation
##
## data: res_ME
## weights: w_reina_c
##
## Moran I statistic standard deviate = 21.778, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2278152679 -0.0003242542 0.0001097365La prueba de moran aplicada sobre los residuales del modelo nos deja ver que aunque la autocorrelación espacial disminuyó en magnitud respecto de un modelo sin estructura espacial, esta sigue siendo altamente significativa. Rechazando la hipótesis de ausencia de correlación espacial, por lo que deben seguirse probando más modelos
Modelo Hurdle - ME
Este modelo es una alternativa al modelo Poisson cuando la variable de conteos a explicar tiene muchos ceros. Este divide las observaciones en dos subgrupos: el primero con todos los ceros y el segundo con los conteos positivos
hurdle_ME <- hurdle(conteos_HR ~ MA90 + MFIL89 + GI89 + FP89 + fitted(me_pois) + offset(log(PO90)) |
RD90 + PS90 + UE90 + DV90 + FH90 + DNL90 + BLK90, data = df_US_homi, dist = "poisson")
summary(hurdle_ME)##
## Call:
## hurdle(formula = conteos_HR ~ MA90 + MFIL89 + GI89 + FP89 + fitted(me_pois) +
## offset(log(PO90)) | RD90 + PS90 + UE90 + DV90 + FH90 + DNL90 + BLK90,
## data = df_US_homi, dist = "poisson")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -9.9799 -0.7728 -0.2830 0.3278 21.5584
##
## Count model coefficients (truncated poisson with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.845e+01 5.943e-01 -64.704 < 2e-16 ***
## MA90 -2.987e-04 2.945e-03 -0.101 0.91921
## MFIL89 1.899e+00 5.399e-02 35.181 < 2e-16 ***
## GI89 2.513e+01 4.562e-01 55.081 < 2e-16 ***
## FP89 -1.137e-02 3.574e-03 -3.181 0.00147 **
## fitted(me_pois)vec71 3.735e+00 3.616e-01 10.329 < 2e-16 ***
## fitted(me_pois)vec165 -1.808e+00 3.564e-01 -5.073 3.92e-07 ***
## fitted(me_pois)vec2951 -1.784e+00 3.316e-01 -5.380 7.43e-08 ***
## fitted(me_pois)vec967 1.674e+00 3.647e-01 4.591 4.42e-06 ***
## fitted(me_pois)vec2096 2.143e+00 4.131e-01 5.187 2.14e-07 ***
## fitted(me_pois)vec267 9.407e-01 4.415e-01 2.131 0.03313 *
## Zero hurdle model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.018380 0.696514 7.205 5.80e-13 ***
## RD90 1.145129 0.135640 8.442 < 2e-16 ***
## PS90 5.048016 0.292858 17.237 < 2e-16 ***
## UE90 -0.044467 0.026370 -1.686 0.0917 .
## DV90 0.271924 0.038049 7.147 8.89e-13 ***
## FH90 -0.009118 0.029796 -0.306 0.7596
## DNL90 -1.464746 0.142870 -10.252 < 2e-16 ***
## BLK90 0.043504 0.008591 5.064 4.10e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of iterations in BFGS optimization: 18
## Log-likelihood: -6246 on 19 Dfres_hurdle_ME_pois <- residuals(hurdle_ME, type = "pearson")
moran.test(res_hurdle_ME_pois, w_reina_c, alternative="two.sided")##
## Moran I test under randomisation
##
## data: res_hurdle_ME_pois
## weights: w_reina_c
##
## Moran I statistic standard deviate = 21.976, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2286337643 -0.0003242542 0.0001085477Este modelo parece ajustarse mejor a los datos dada la separaciónd e las covariables en los dos grupos previamente mencionados, además los vectores propios de Moran son significativos en la regresión. No obstante, el test de Moran deja ver que se sigue teniendo presencia de autocorrelación espacial en los datos. Es por esto que a continuación se explora la autocorrelación de las variables explicativas.
Modelo Poisson SLX
El modelo SLX (Spatial Lag of X) considera que las variables explicativas estén rezagadas espacialmente (\(\mathbf{WX}\theta\)) donde \(\mathbf{W}\) es la matriz de vecindad y \(\theta\) es el vector de coeficientes de las variables rezagadas.
Wmat <- listw2mat(w_reina_c)
df_US_homi$W_conteos_HR <- Wmat %*% df_US_homi$conteos_HR
df_US_homi$W_HR90 <- Wmat %*% df_US_homi$HR90
df_US_homi$W_RD90 <- Wmat %*% df_US_homi$RD90
df_US_homi$W_PS90 <- Wmat %*% df_US_homi$PS90
df_US_homi$W_UE90 <- Wmat %*% df_US_homi$UE90
df_US_homi$W_DV90 <- Wmat %*% df_US_homi$DV90
df_US_homi$W_FH90 <- Wmat %*% df_US_homi$FH90
df_US_homi$W_MA90 <- Wmat %*% df_US_homi$MA90
df_US_homi$W_DNL90 <- Wmat %*% df_US_homi$DNL90
df_US_homi$W_MFIL89 <- Wmat %*% df_US_homi$MFIL89
df_US_homi$W_FP89 <- Wmat %*% df_US_homi$FP89
df_US_homi$W_BLK90 <- Wmat %*% df_US_homi$BLK90
df_US_homi$W_GI89 <- Wmat %*% df_US_homi$GI89ecuacion_lag = conteos_HR ~ RD90 + PS90 + UE90 + DV90 + FH90 + DNL90 + MFIL89 +
FP89 + BLK90 + GI89 + MA90 + offset(log(PO90)) + W_conteos_HR + W_RD90 + W_PS90 + W_UE90 + W_DV90 + W_FH90 + W_DNL90 + W_MFIL89 +
W_FP89 + W_BLK90 + W_GI89 + W_MA90 + fitted(me_pois)
glm_lag_pois <- glm(ecuacion_lag, data = df_US_homi, family = "poisson")
summary(glm_lag_pois)##
## Call:
## glm(formula = ecuacion_lag, family = "poisson", data = df_US_homi)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -9.888e+07 1.061e+08 -0.932 0.351449
## RD90 1.662e+07 1.784e+07 0.932 0.351449
## PS90 4.404e-01 3.924e-02 11.224 < 2e-16 ***
## UE90 1.983e-02 6.844e-03 2.897 0.003772 **
## DV90 8.814e-02 7.301e-03 12.072 < 2e-16 ***
## FH90 -7.041e+05 7.557e+05 -0.932 0.351449
## DNL90 -1.208e-01 2.801e-02 -4.313 1.61e-05 ***
## MFIL89 1.644e+07 1.765e+07 0.932 0.351449
## FP89 -6.763e+05 7.258e+05 -0.932 0.351449
## BLK90 -2.341e+05 2.512e+05 -0.932 0.351449
## GI89 -1.318e+08 1.415e+08 -0.932 0.351449
## MA90 -5.554e-04 3.726e-03 -0.149 0.881495
## W_conteos_HR 1.106e-04 1.471e-04 0.752 0.452226
## W_RD90 -1.218e+00 3.133e-01 -3.889 0.000101 ***
## W_PS90 9.362e-02 6.745e-02 1.388 0.165140
## W_UE90 -8.690e-02 9.629e-03 -9.024 < 2e-16 ***
## W_DV90 3.231e-02 1.040e-02 3.108 0.001883 **
## W_FH90 2.870e-02 1.525e-02 1.882 0.059873 .
## W_DNL90 2.570e-03 4.044e-02 0.064 0.949335
## W_MFIL89 -5.921e-01 1.604e-01 -3.691 0.000223 ***
## W_FP89 9.708e-02 1.639e-02 5.924 3.15e-09 ***
## W_BLK90 2.363e-02 4.660e-03 5.071 3.96e-07 ***
## W_GI89 1.116e+01 3.166e+00 3.525 0.000424 ***
## W_MA90 9.534e-03 5.179e-03 1.841 0.065613 .
## fitted(me_pois)vec71 2.123e+00 4.037e-01 5.258 1.46e-07 ***
## fitted(me_pois)vec165 -4.219e+00 4.378e-01 -9.637 < 2e-16 ***
## fitted(me_pois)vec2951 -1.349e+00 3.914e-01 -3.445 0.000570 ***
## fitted(me_pois)vec967 1.446e+00 3.775e-01 3.831 0.000128 ***
## fitted(me_pois)vec2096 1.298e+00 4.208e-01 3.085 0.002033 **
## fitted(me_pois)vec267 1.901e+00 4.739e-01 4.012 6.03e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 17535.7 on 3084 degrees of freedom
## Residual deviance: 3316.2 on 3055 degrees of freedom
## AIC: 9574.3
##
## Number of Fisher Scoring iterations: 5res_lag_pois <- residuals(glm_lag_pois)
moran.test(res_lag_pois, w_reina_c, alternative="two.sided")##
## Moran I test under randomisation
##
## data: res_lag_pois
## weights: w_reina_c
##
## Moran I statistic standard deviate = 21.125, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2209911381 -0.0003242542 0.0001097596Note que pese a la adición de la estructura de autocorrelación espacial de todas las covariables y los vectores propios de Moran no se logra modelar la dependencia exacta de correlación. Esto quiere decir que poner la variable respuesta rezagada como una covariable más NO corrige la autocorrelación de los errores ni de la Tasa de homicidios. Nos vemos en la necesidad de estimar un modelo SAR-Poisson que, teoricamente, si corrige todos estos inconvenientes.
Modelo SAR Poisson
El modelo SAR (Spatial Autoregressive) incluye el término de rezago en la variable dependiente (\(\mathbf{Wy}\)) donde \(\rho\) es el coeficiente asociado al lag espacial. El término del error también está modelado como un proceso autoregresivo espacial.
Para el presente caso de la distribución Poisson \[Y_i\sim\text{Poisson}(\mu_i)\\ \log(\mu_i)=\mathbf{X}_i\beta+\rho(WY)_i\]
#ecuacion_lag = conteos_HR ~ RD90 + PS90 + UE90 + DV90 + FH90 + DNL90 + MFIL89 +
# FP89 + BLK90 + GI89 + MA90 + offset(log(PO90)) +
# W_RD90 + W_PS90 + W_UE90 + W_DV90 + W_FH90 + W_DNL90 + W_MFIL89 +
# W_FP89 + W_BLK90 + W_GI89 + W_MA90 + fitted(me_pois)
#sar_pois <- sarpoisson(ecuacion_lag, data = df_US_homi, listw = w_reina_c, method = "fiml")#res_SAR <- as.numeric(residuals(sar_pois, type = "pearson"))
#moran.test(res_SAR, w_reina_c, alternative="two.sided")
# Moran I test under randomisation
#data: res_SAR
#weights: w_reina_c
#Moran I statistic standard deviate = 4.986, p-value = 6.163e-07
#alternative hypothesis: two.sided
#sample estimates:
#Moran I statistic Expectation Variance
# 4.655880e-02 -3.242542e-04 8.841405e-05 Se probaron 3 tipos de modelos SAR: modelo SAR sencillo, modelo SAR con los vectores propios de Moran y modelo SAR con los vectores propios de Moran más las variables explicativas rezagadas espacialmente.Tristemente, ninguno de estos 3 mdoelos logró corregir el problema de la autocorrelación espacial en los residuales. Por esto, y por el alto costo computacional que estos modelos representan, no se incluyen en el documento final y solo se deja el resultado del Test de Moran comentado como referencia.
Modelos de regresión Binomiales Negativos
Dentro de la búsqueda de los modelos más adecuados para los datos en estudio se exploran a continuación los modelos análogos a los presentados previamente pero con la salvedad de que ahora se usará la distribución Binomial negativa que también es útil para modelar conteos cuando la varianza no coincide con la media (supuesto de la distribución Poisson).
Modelo BinNeg OLS
Modelo lineal generalizado de la familia bonimal negativa con las covariables originales, sin tener en cuenta la estructura espacial.
glm_negbin <- glm.nb(ecuacion, data = df_US_homi)## Warning in glm.nb(ecuacion, data = df_US_homi): alternation limit reachedsummary(glm_negbin)##
## Call:
## glm.nb(formula = ecuacion, data = df_US_homi, init.theta = 15.83107966,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.512e+08 1.793e+08 -0.844 0.398942
## RD90 2.542e+07 3.013e+07 0.844 0.398942
## PS90 4.358e-01 5.011e-02 8.698 < 2e-16 ***
## UE90 -2.082e-02 8.127e-03 -2.562 0.010405 *
## DV90 1.158e-01 8.309e-03 13.935 < 2e-16 ***
## FH90 -1.077e+06 1.277e+06 -0.844 0.398942
## DNL90 -1.010e-01 2.981e-02 -3.389 0.000702 ***
## MFIL89 2.515e+07 2.981e+07 0.844 0.398942
## FP89 -1.034e+06 1.226e+06 -0.844 0.398942
## BLK90 -3.580e+05 4.244e+05 -0.844 0.398942
## GI89 -2.016e+08 2.390e+08 -0.844 0.398942
## MA90 8.422e-03 4.486e-03 1.877 0.060469 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(15.8311) family taken to be 1)
##
## Null deviance: 5930.6 on 3084 degrees of freedom
## Residual deviance: 2454.6 on 3073 degrees of freedom
## AIC: 9282.1
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 15.83
## Std. Err.: 1.56
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -9256.133res_negbin <- residuals(glm_negbin, type = "pearson")
moran.test(res_negbin, w_reina_c, alternative="two.sided")##
## Moran I test under randomisation
##
## data: res_negbin
## weights: w_reina_c
##
## Moran I statistic standard deviate = 22.403, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2344044179 -0.0003242542 0.0001097801Modelo NegBin- Moran Eigenvectors
Se añaden ahora los vectores propios de la Matriz de Moran como covariables al modelo lineal generalizado binomial negativo. Se hace la claridad de que se agregan los vectroes ME estimados con el modelo poisson porque no hay forma de estimarlos directamente desde la familia binomial negativa.
glm_ME_negbin <- glm.nb(ecuacion2, data = df_US_homi)## Warning in glm.nb(ecuacion2, data = df_US_homi): alternation limit reachedsummary(glm_ME_negbin)##
## Call:
## glm.nb(formula = ecuacion2, data = df_US_homi, init.theta = 18.09448823,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.548e+08 1.751e+08 -0.884 0.376707
## RD90 2.601e+07 2.943e+07 0.884 0.376707
## PS90 3.776e-01 4.924e-02 7.669 1.74e-14 ***
## UE90 -2.217e-02 8.017e-03 -2.765 0.005685 **
## DV90 1.177e-01 8.067e-03 14.592 < 2e-16 ***
## FH90 -1.102e+06 1.247e+06 -0.884 0.376707
## DNL90 -6.431e-02 2.917e-02 -2.204 0.027507 *
## MFIL89 2.574e+07 2.912e+07 0.884 0.376707
## FP89 -1.058e+06 1.197e+06 -0.884 0.376707
## BLK90 -3.664e+05 4.145e+05 -0.884 0.376707
## GI89 -2.063e+08 2.334e+08 -0.884 0.376707
## fitted(me_pois)vec71 1.268e+00 6.589e-01 1.924 0.054311 .
## fitted(me_pois)vec165 -2.582e+00 7.112e-01 -3.631 0.000283 ***
## fitted(me_pois)vec2951 -1.218e+00 6.450e-01 -1.888 0.059030 .
## fitted(me_pois)vec967 1.397e+00 6.227e-01 2.243 0.024893 *
## fitted(me_pois)vec2096 1.199e+00 6.798e-01 1.764 0.077705 .
## fitted(me_pois)vec267 1.037e+00 6.752e-01 1.536 0.124586
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(18.0945) family taken to be 1)
##
## Null deviance: 6195.9 on 3084 degrees of freedom
## Residual deviance: 2483.0 on 3068 degrees of freedom
## AIC: 9266.4
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 18.09
## Std. Err.: 1.95
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -9230.447res_ME_negbin <- residuals(glm_ME_negbin, type = "pearson")
moran.test(res_ME_negbin, w_reina_c, alternative="two.sided")##
## Moran I test under randomisation
##
## data: res_ME_negbin
## weights: w_reina_c
##
## Moran I statistic standard deviate = 22.201, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2322856014 -0.0003242542 0.0001097777Los reaiduales del modelo siguen presentando autocorrelación espacial significativa. Lo mismo que sucedía con el modelo de la familia Poisson.
Modelo Hurdle (NegBin) - ME
Se explroa el modelo Hurdle (explicado enn la sección anterior) pero ahora con el supuesto de que los datos de conteo vienen de una distribución binomial negativa. Se maneja la misma estructura del modelo hurdle poisson para la modelación de los ceros.
hurdle_ME <- hurdle(conteos_HR ~ MA90 + MFIL89 + GI89 + FP89 + fitted(me_pois) + offset(log(PO90)) |
RD90 + PS90 + UE90 + DV90 + FH90 + DNL90 + BLK90, data = df_US_homi, dist = "negbin")
summary(hurdle_ME)##
## Call:
## hurdle(formula = conteos_HR ~ MA90 + MFIL89 + GI89 + FP89 + fitted(me_pois) +
## offset(log(PO90)) | RD90 + PS90 + UE90 + DV90 + FH90 + DNL90 + BLK90,
## data = df_US_homi, dist = "negbin")
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## -1.9487 -0.6091 -0.2294 0.3611 7.6212
##
## Count model coefficients (truncated negbin with log link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -39.179294 1.864020 -21.019 < 2e-16 ***
## MA90 0.003545 0.006464 0.548 0.58339
## MFIL89 2.040950 0.157391 12.967 < 2e-16 ***
## GI89 21.775399 1.189534 18.306 < 2e-16 ***
## FP89 0.009757 0.007831 1.246 0.21281
## fitted(me_pois)vec71 2.287140 0.953997 2.397 0.01651 *
## fitted(me_pois)vec165 -0.857525 1.068269 -0.803 0.42213
## fitted(me_pois)vec2951 -0.872645 0.920667 -0.948 0.34321
## fitted(me_pois)vec967 2.711940 0.887985 3.054 0.00226 **
## fitted(me_pois)vec2096 1.094267 1.019135 1.074 0.28295
## fitted(me_pois)vec267 0.688843 0.968354 0.711 0.47686
## Log(theta) 1.469237 0.075670 19.416 < 2e-16 ***
## Zero hurdle model coefficients (binomial with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 5.018380 0.696514 7.205 5.80e-13 ***
## RD90 1.145129 0.135640 8.442 < 2e-16 ***
## PS90 5.048016 0.292858 17.237 < 2e-16 ***
## UE90 -0.044467 0.026370 -1.686 0.0917 .
## DV90 0.271924 0.038049 7.147 8.89e-13 ***
## FH90 -0.009118 0.029796 -0.306 0.7596
## DNL90 -1.464746 0.142870 -10.252 < 2e-16 ***
## BLK90 0.043504 0.008591 5.064 4.10e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Theta: count = 4.3459
## Number of iterations in BFGS optimization: 16
## Log-likelihood: -4874 on 20 Dfres_hurdle_ME_neg <- residuals(hurdle_ME, type = "pearson")
moran.test(res_hurdle_ME_neg, w_reina_c, alternative="two.sided")##
## Moran I test under randomisation
##
## data: res_hurdle_ME_neg
## weights: w_reina_c
##
## Moran I statistic standard deviate = 24.917, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.2605152584 -0.0003242542 0.0001095894No se observan mejorías en la estructura espacial de los residuales. Sigue sin capturarse su estrutura.
Modelo SLX NegBin
Por último, se explora el modelo glm binomial negativo que recibe como covaribles a las variables explicativas rezagadas espacialmente por la matriz de pesos ‘w_reina_c’ junto con los vectores ME.
SLX_negbin <- glm.nb(ecuacion_lag, data = df_US_homi)## Warning in glm.nb(ecuacion_lag, data = df_US_homi): alternation limit reachedsummary(SLX_negbin)##
## Call:
## glm.nb(formula = ecuacion_lag, data = df_US_homi, init.theta = 19.41663054,
## link = log)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.254e+07 1.751e+08 -0.129 0.897563
## RD90 3.788e+06 2.942e+07 0.129 0.897563
## PS90 3.352e-01 6.841e-02 4.900 9.58e-07 ***
## UE90 1.474e-02 9.932e-03 1.484 0.137704
## DV90 9.862e-02 1.106e-02 8.915 < 2e-16 ***
## FH90 -1.605e+05 1.247e+06 -0.129 0.897563
## DNL90 -6.623e-02 4.699e-02 -1.410 0.158683
## MFIL89 3.748e+06 2.911e+07 0.129 0.897563
## FP89 -1.541e+05 1.197e+06 -0.129 0.897563
## BLK90 -5.335e+04 4.144e+05 -0.129 0.897563
## GI89 -3.004e+07 2.333e+08 -0.129 0.897563
## MA90 8.293e-03 5.123e-03 1.619 0.105518
## W_conteos_HR 7.731e-04 2.878e-04 2.686 0.007230 **
## W_RD90 -9.704e-01 4.725e-01 -2.054 0.039994 *
## W_PS90 1.051e-01 1.031e-01 1.019 0.308189
## W_UE90 -8.993e-02 1.290e-02 -6.971 3.16e-12 ***
## W_DV90 5.048e-02 1.488e-02 3.392 0.000695 ***
## W_FH90 5.683e-03 2.309e-02 0.246 0.805554
## W_DNL90 8.829e-03 6.153e-02 0.143 0.885910
## W_MFIL89 -5.555e-01 2.401e-01 -2.313 0.020706 *
## W_FP89 6.524e-02 2.531e-02 2.578 0.009950 **
## W_BLK90 2.433e-02 7.163e-03 3.396 0.000684 ***
## W_GI89 1.298e+01 4.636e+00 2.799 0.005119 **
## W_MA90 -1.773e-03 8.492e-03 -0.209 0.834572
## fitted(me_pois)vec71 9.031e-01 6.653e-01 1.357 0.174633
## fitted(me_pois)vec165 -2.166e+00 7.310e-01 -2.963 0.003043 **
## fitted(me_pois)vec2951 -1.224e+00 6.354e-01 -1.927 0.053946 .
## fitted(me_pois)vec967 1.564e+00 6.223e-01 2.514 0.011932 *
## fitted(me_pois)vec2096 1.010e+00 6.751e-01 1.496 0.134771
## fitted(me_pois)vec267 7.124e-01 6.711e-01 1.062 0.288455
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Negative Binomial(19.4166) family taken to be 1)
##
## Null deviance: 6338.8 on 3084 degrees of freedom
## Residual deviance: 2381.2 on 3055 degrees of freedom
## AIC: 9163.8
##
## Number of Fisher Scoring iterations: 1
##
##
## Theta: 19.42
## Std. Err.: 2.07
## Warning while fitting theta: alternation limit reached
##
## 2 x log-likelihood: -9101.823res_SLX_negbin <- residuals(SLX_negbin, type = "pearson")
moran.test(res_SLX_negbin, w_reina_c, alternative="two.sided")##
## Moran I test under randomisation
##
## data: res_SLX_negbin
## weights: w_reina_c
##
## Moran I statistic standard deviate = 18.958, p-value < 2.2e-16
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.1983004280 -0.0003242542 0.0001097659La inclusión de las covariables rezagadas espacialmente no mejoró la corelación que los datos están presentando. Esto se debe a que no se está modelando directamente el carácter espacial de la variable respuesta ni de los errores del modelo (componentes que normalmente son los que cargan con la mayoria de la estructura espacial).
Modelo de regresión S-CAR Bayesiano
Debido a que todos los modelos clásicos no lograron modelar correctamente la autocorrelación espacial se prueba ahora un Modelo lineal generalizado mixto espacial, cuyos efectos aleatorios siguen una prori autoregresiva Bayesiana.
Este modelo se puede escribir como \[Y_i|\mu_i\sim \text{Poisson}(\mu_i)\\ \log(\mu_i)=\log(E_i)+\beta_0+\mathbf{X}_i^T\beta+u_i+v_i\]
donde \(Y_i\) es el número de homicidios en el área \(i\), \(\log(E_i)\) es el offset de población por condado que permite trabajar los conteos como tasas, \(u_i\) es el efecto espacial CAR que captura la correlación entre áreas vecinas y \(v_i\) es el efecto no estructurado de variación no espacial.
\[u_i|u_{-i}, \tau_u\sim N\left(\frac{1}{n_i}\sum_{j\sim i}u_j, \frac{1}{\tau_u n_i}\right)\\v_i\sim N(0, \tau_v^{-1})\\ n_i=\text{número de vecinos de i}\]
donde los hiperparámetros tienen previas \[\tau_u\sim\text{Gamma}(a, b)\quad \tau_v\sim\text{Gamma}(a, b)\quad \beta\sim N(0,M)\]
modCAR <- S.CARbym(conteos_HR ~ RD90 + PS90 + UE90 + DV90 + FH90 + DNL90 +
FP89 + BLK90 + MA90 + offset(log(PO90)),
family = "poisson",
data = df_US_homi,
W = Wmat,
burnin = 10000, n.sample = 40000)## Setting up the model.
##
## Markov chain 1 - generating 30000 post burnin and thinned samples.
## | | | 0% | |= | 1% | |= | 2% | |== | 3% | |=== | 4% | |==== | 5% | |==== | 6% | |===== | 7% | |====== | 8% | |====== | 9% | |======= | 10% | |======== | 11% | |======== | 12% | |========= | 13% | |========== | 14% | |========== | 15% | |=========== | 16% | |============ | 17% | |============= | 18% | |============= | 19% | |============== | 20% | |=============== | 21% | |=============== | 22% | |================ | 23% | |================= | 24% | |================== | 25% | |================== | 26% | |=================== | 27% | |==================== | 28% | |==================== | 29% | |===================== | 30% | |====================== | 31% | |====================== | 32% | |======================= | 33% | |======================== | 34% | |======================== | 35% | |========================= | 36% | |========================== | 37% | |=========================== | 38% | |=========================== | 39% | |============================ | 40% | |============================= | 41% | |============================= | 42% | |============================== | 43% | |=============================== | 44% | |================================ | 45% | |================================ | 46% | |================================= | 47% | |================================== | 48% | |================================== | 49% | |=================================== | 50% | |==================================== | 51% | |==================================== | 52% | |===================================== | 53% | |====================================== | 54% | |====================================== | 55% | |======================================= | 56% | |======================================== | 57% | |========================================= | 58% | |========================================= | 59% | |========================================== | 60% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 63% | |============================================= | 64% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 67% | |================================================ | 68% | |================================================ | 69% | |================================================= | 70% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 73% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 76% | |====================================================== | 77% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 80% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 83% | |=========================================================== | 84% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 87% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 90% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 93% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 97% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 100%
##
## Summarising results.
## Finished in 396.6 seconds.res_car <- modCAR$residuals[,2] # residuales de pearson
moran.test(res_car, w_reina_c, alternative = "two.sided")##
## Moran I test under randomisation
##
## data: res_car
## weights: w_reina_c
##
## Moran I statistic standard deviate = 0.67607, p-value = 0.499
## alternative hypothesis: two.sided
## sample estimates:
## Moran I statistic Expectation Variance
## 0.0067589347 -0.0003242542 0.0001097668Como el p-valor del test de Moran es de 0.55, a un nivel de significancia del 5% NO se rechaza la hipótesis nula por lo que se concluye que los residuales del modelo ya no tienen autocorrelación espacial.
Observe que, este resulta ser en el mejor modelo de todos y también el más parsimonioso, debido a que basta con darle como covariables a las variables explicativas originales para que se desempeñe de forma excelente.
Medidas de desempeño del modelo S-CAR
A continuación se presentan unas medidas de desempeño del modelo Espacial Bayesiano CAR.}
Nota: Como ningún otro modelo logró explicar correctamente la correlación espacial además del S-CAR bayesiano, no es justo compararlo con los demás modelos ajustados.
cat("\n========== Resultados Modelo CAR ==========\n")##
## ========== Resultados Modelo CAR ==========cat("MAD: ", mad(abs(fitted(modCAR)-df_US_homi$conteos_HR)), "\n")## MAD: 0.5100683cat("MSE: ", mean((fitted(modCAR)-df_US_homi$conteos_HR)^2), "\n")## MSE: 2.235655PseudoR2_CAR <- 1 - sum((fitted(modCAR) - df_US_homi$conteos_HR)^2) /
(length(df_US_homi$conteos_HR) * var(df_US_homi$conteos_HR))
cat("Pseudo-R2: ", PseudoR2_CAR, "\n")## Pseudo-R2: 0.9993956modCAR$modelfit## DIC p.d WAIC p.w LMPL
## 8643.6357 350.4952 8572.7554 228.8778 -4342.8500
## loglikelihood
## -3971.3227cat("===========================================\n")## ===========================================mean((fitted(hurdle_ME)-df_US_homi$conteos_HR)^2)## [1] 325.0633Solo como comparación (un poco injusta) con el modelo Hurdle ME con distribución Binomia negativa, se observa que el error cuadrático medio (MSE) es muchísimo menor para el modelo bayesiano porque este si está explicando debidamente la estructura de los datos.
modCAR$summary.results## Mean 2.5% 97.5% n.sample % accept n.effective
## (Intercept) -11.2291 -11.8032 -10.6818 30000 47.3 78.3
## RD90 0.2169 0.0956 0.3650 30000 47.3 14.0
## PS90 0.4360 0.3404 0.5622 30000 47.3 15.1
## UE90 0.0230 0.0048 0.0404 30000 47.3 100.0
## DV90 0.0957 0.0733 0.1171 30000 47.3 164.8
## FH90 0.0207 0.0046 0.0354 30000 47.3 25.1
## DNL90 -0.1393 -0.2263 -0.0753 30000 47.3 11.8
## FP89 0.0042 -0.0115 0.0182 30000 47.3 18.3
## BLK90 0.0098 0.0052 0.0146 30000 47.3 38.1
## MA90 0.0121 0.0025 0.0211 30000 47.3 213.1
## tau2 0.0260 0.0210 0.0318 30000 100.0 251.9
## sigma2 0.0046 0.0018 0.0091 30000 100.0 23.9
## Geweke.diag
## (Intercept) 0.1
## RD90 -0.1
## PS90 -0.3
## UE90 1.1
## DV90 0.1
## FH90 -0.2
## DNL90 0.1
## FP89 -0.1
## BLK90 0.4
## MA90 -0.6
## tau2 0.0
## sigma2 1.4Diagnóstico de los residuales del modelo final
hist(modCAR$residuals[,2], main = "Histograma de los residuales de Pearson del modelo SCAR \n Conteos de homicidios en EEUU, 1990",
xlab = "Residuales de Pearson", ylab = "Frecuencia")
plot(fitted(modCAR), df_US_homi$conteos_HR, pch=20, xlab="Conteos de homicidios ajustados", ylab="Conteos de homicidios observados",
main="Valores observados VS Valores Predichos\n Modelo SCAR Bayesiano para los conteos de homicidios en EEUU")
abline(0,1, col="red")
EL histograma de los residuales de pearson (residuales estandarizados) nos deja ver que estos se comportan como ruido blanco centrados alrededor de cero. Así mismo, el segundo gráfico nos deja ver el buen ajuste que tiene el modelo, manteniéndose casi totalmente sobre la diagonal de ajuste ideal.
Mapa de valores predichos
tasa_ajustada <- modCAR$fitted.values / df_US_homi$PO90 * 100000
ggplot(US_homi_reg) +
geom_sf(aes(fill = tasa_ajustada)) +
scale_fill_viridis_c(direction = -1, option='B') + # Escala de colores para la variable
labs(fill = "Tasa ajustadas") +
theme_minimal() +
ggtitle("Valores ajustados mediante el modelo CAR Bayesiano\n Tasas de homicidio por condados en EEUU, año 1990") +
# Añadir la leyenda personalizada
theme(legend.position = "right") +
# Añadir escala
annotation_scale(location = "bl", width_hint = 0.15) +
annotation_north_arrow(location = "tr", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"))## Scale on map varies by more than 10%, scale bar may be inaccurate
Para visualizar mejor las Tasas ajustadas por el modelo CAR Bayesiano se construye el gráfico BoxMap
US_homi_reg$tasa_ajustada <- modCAR$fitted.values / df_US_homi$PO90 * 100000
# 2. Calcular cortes del boxplot
Q1 <- quantile(US_homi_reg$tasa_ajustada, 0.25, na.rm = TRUE)
Q3 <- quantile(US_homi_reg$tasa_ajustada, 0.75, na.rm = TRUE)
IQR <- Q3 - Q1
lower_fence <- Q1 - 1.5 * IQR
upper_fence <- Q3 + 1.5 * IQR
# 3. Clasificación tipo boxmap
US_homi_reg$box_class <- cut(
US_homi_reg$tasa_ajustada,
breaks = c(-Inf, lower_fence, Q1, Q3, upper_fence, Inf),
labels = c("Outlier Bajo", "Bajo", "Medio", "Alto", "Outlier Alto"),
include.lowest = TRUE
)
# 4. Paleta personalizada para categorías
colores_boxmap <- viridis::viridis(4, option = "B", direction = -1)
# 5. Gráfico tipo boxmap
ggplot(US_homi_reg) +
geom_sf(aes(fill = box_class), color = "gray40", size = 0.1) +
scale_fill_manual(values = colores_boxmap) +
labs(fill = "Boxmap (tasa ajustada)",
title = "Boxmap de Tasa de Homicidio Ajustada (1990)") +
theme_minimal() +
theme(legend.position = "right")
Los datos clasificados en el BoxMap como medios corresponden a valores que se encuentran dentro del 50% de valores medios (entre el Q1 y el Q3), mientras que los datos clasificados como Bajo y ALto corresponden a aquellos fuera del “boxplot” y dentro de los bigotes. Se observa que al igual que las Tasas de homicidio originales, el modelo nos muestra que hacia el sur este de Estados Unidos se tienen las tasas de homicidos más altas consideradas valores atípicos. Del mismo modo, si miramos detalladamente, hacia la parte alta central del país los valores predichos para 2 condados se consideran outliers, coincidiendo con outliers marcados por las tasas de homicidio originales.