Regresión Lineal Espacial

Actividad 03

Planeación estratégica basada en analítica prescriptiva

Grupo 503

Profesor Rodolfo Miguel Gameros

Equipo 7:

A00833113 - Avril Lobato

A01771127 - Lesly Darian Romero Vázquez

A00831105 - Jazmín Cortez Mendoza

A01284611 - Lisset Hernández

Análisis Exploratorio de Datos Tradicional

Importación de base de datos y geocercas

El conjunto de datos boston.c contiene 506 filas y 20 columnas. Incluye los datos de Harrison y Rubinfeld (1978), corregidos por pequeños errores y complementados con la latitud y longitud de las observaciones. Gilley y Pace también señalan que la variable MEDV está censurada, es decir, los valores medianos iguales o superiores a USD 50,000 se establecen en USD 50,000.

Este data frame contiene las siguientes columnas:

• TOWN: Factor con los nombres de las localidades.

• TOWNNO: Vector numérico correspondiente a TOWN.

• TRACT: Vector numérico con los identificadores de los tramos censales.

• LON: Longitud del punto central del tramo censal en grados decimales.

• LAT: Latitud del punto central del tramo censal en grados decimales.

• MEDV: Valor mediano de viviendas ocupadas por sus propietarios en miles de USD.

• CMEDV: Valor mediano corregido de viviendas ocupadas por sus propietarios en miles de USD.

• CRIM: Tasa de criminalidad per cápita.

• ZN: Proporción de terreno residencial zonificado para lotes mayores a 25,000 pies² (constante para todos los tramos de Boston).

• INDUS: Proporción de acres de negocios no minoristas por localidad (constante para todos los tramos de Boston).

• CHAS: Factor con niveles: 1 si el tramo colinda con el río Charles; 0 en caso contrario.

• NOX: Concentración de óxidos nítricos (partes por 10 millones) por localidad.

• RM: Promedio de habitaciones por vivienda.

• AGE: Proporción de unidades ocupadas por propietarios construidas antes de 1940.

• DIS: Distancias ponderadas a cinco centros de empleo en Boston.

• RAD: Índice de accesibilidad a autopistas radiales por localidad (constante para todos los tramos de Boston).

• TAX: Tasa de impuesto predial por cada 10,000 USD de valor de propiedad (constante para todos los tramos de Boston).

• PTRATIO: Relación alumno/profesor por localidad (constante para todos los tramos de Boston).

• B: Cálculo de 1000*(Bk - 0.63)^2, donde Bk es la proporción de personas negras.

• LSTAT: Porcentaje de población de menor nivel socioeconómico.

data(boston, package="spData")
head(boston.c)

##         TOWN TOWNNO TRACT      LON     LAT MEDV CMEDV    CRIM ZN INDUS CHAS
## 1     Nahant      0  2011 -70.9550 42.2550 24.0  24.0 0.00632 18  2.31    0
## 2 Swampscott      1  2021 -70.9500 42.2875 21.6  21.6 0.02731  0  7.07    0
## 3 Swampscott      1  2022 -70.9360 42.2830 34.7  34.7 0.02729  0  7.07    0
## 4 Marblehead      2  2031 -70.9280 42.2930 33.4  33.4 0.03237  0  2.18    0
## 5 Marblehead      2  2032 -70.9220 42.2980 36.2  36.2 0.06905  0  2.18    0
## 6 Marblehead      2  2033 -70.9165 42.3040 28.7  28.7 0.02985  0  2.18    0
##     NOX    RM  AGE    DIS RAD TAX PTRATIO      B LSTAT
## 1 0.538 6.575 65.2 4.0900   1 296    15.3 396.90  4.98
## 2 0.469 6.421 78.9 4.9671   2 242    17.8 396.90  9.14
## 3 0.469 7.185 61.1 4.9671   2 242    17.8 392.83  4.03
## 4 0.458 6.998 45.8 6.0622   3 222    18.7 394.63  2.94
## 5 0.458 7.147 54.2 6.0622   3 222    18.7 396.90  5.33
## 6 0.458 6.430 58.7 6.0622   3 222    18.7 394.12  5.21

glimpse(boston.c)

## Rows: 506
## Columns: 20
## $ TOWN    <fct> Nahant, Swampscott, Swampscott, Marblehead, Marblehead, Marble…
## $ TOWNNO  <int> 0, 1, 1, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, 4, 4,…
## $ TRACT   <int> 2011, 2021, 2022, 2031, 2032, 2033, 2041, 2042, 2043, 2044, 20…
## $ LON     <dbl> -70.9550, -70.9500, -70.9360, -70.9280, -70.9220, -70.9165, -7…
## $ LAT     <dbl> 42.2550, 42.2875, 42.2830, 42.2930, 42.2980, 42.3040, 42.2970,…
## $ MEDV    <dbl> 24.0, 21.6, 34.7, 33.4, 36.2, 28.7, 22.9, 27.1, 16.5, 18.9, 15…
## $ CMEDV   <dbl> 24.0, 21.6, 34.7, 33.4, 36.2, 28.7, 22.9, 22.1, 16.5, 18.9, 15…
## $ CRIM    <dbl> 0.00632, 0.02731, 0.02729, 0.03237, 0.06905, 0.02985, 0.08829,…
## $ ZN      <dbl> 18.0, 0.0, 0.0, 0.0, 0.0, 0.0, 12.5, 12.5, 12.5, 12.5, 12.5, 1…
## $ INDUS   <dbl> 2.31, 7.07, 7.07, 2.18, 2.18, 2.18, 7.87, 7.87, 7.87, 7.87, 7.…
## $ CHAS    <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
## $ NOX     <dbl> 0.538, 0.469, 0.469, 0.458, 0.458, 0.458, 0.524, 0.524, 0.524,…
## $ RM      <dbl> 6.575, 6.421, 7.185, 6.998, 7.147, 6.430, 6.012, 6.172, 5.631,…
## $ AGE     <dbl> 65.2, 78.9, 61.1, 45.8, 54.2, 58.7, 66.6, 96.1, 100.0, 85.9, 9…
## $ DIS     <dbl> 4.0900, 4.9671, 4.9671, 6.0622, 6.0622, 6.0622, 5.5605, 5.9505…
## $ RAD     <int> 1, 2, 2, 3, 3, 3, 5, 5, 5, 5, 5, 5, 5, 4, 4, 4, 4, 4, 4, 4, 4,…
## $ TAX     <int> 296, 242, 242, 222, 222, 222, 311, 311, 311, 311, 311, 311, 31…
## $ PTRATIO <dbl> 15.3, 17.8, 17.8, 18.7, 18.7, 18.7, 15.2, 15.2, 15.2, 15.2, 15…
## $ B       <dbl> 396.90, 396.90, 392.83, 394.63, 396.90, 394.12, 395.60, 396.90…
## $ LSTAT   <dbl> 4.98, 9.14, 4.03, 2.94, 5.33, 5.21, 12.43, 19.15, 29.93, 17.10…

Seguidamente, se cargan las geocercas asociadas a la data transversal previamente cargada

boston.tr<-st_read(system.file("shapes/boston_tracts.gpkg", package="spData")[1])

## Reading layer `boston_tracts' from data source 
##   `C:\Users\AVRIL\AppData\Local\R\win-library\4.3\spData\shapes\boston_tracts.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 506 features and 36 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -71.52311 ymin: 42.00305 xmax: -70.63823 ymax: 42.67307
## Geodetic CRS:  NAD27

Summary Boston

boston.trSP<-as(boston.tr, "Spatial")
boston_nb<-poly2nb(boston.trSP)
summary(boston.trSP)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##         min       max
## x -71.52311 -70.63823
## y  42.00305  42.67307
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=NAD27 +no_defs]
## Data attributes:
##    poltract             TOWN               TOWNNO          TRACT     
##  Length:506         Length:506         Min.   : 0.00   Min.   :   1  
##  Class :character   Class :character   1st Qu.:26.25   1st Qu.:1303  
##  Mode  :character   Mode  :character   Median :42.00   Median :3394  
##                                        Mean   :47.53   Mean   :2700  
##                                        3rd Qu.:78.00   3rd Qu.:3740  
##                                        Max.   :91.00   Max.   :5082  
##                                                                      
##       LON              LAT             MEDV           CMEDV      
##  Min.   :-71.29   Min.   :42.03   Min.   : 5.00   Min.   : 5.00  
##  1st Qu.:-71.09   1st Qu.:42.18   1st Qu.:17.02   1st Qu.:17.02  
##  Median :-71.05   Median :42.22   Median :21.20   Median :21.20  
##  Mean   :-71.06   Mean   :42.22   Mean   :22.53   Mean   :22.53  
##  3rd Qu.:-71.02   3rd Qu.:42.25   3rd Qu.:25.00   3rd Qu.:25.00  
##  Max.   :-70.81   Max.   :42.38   Max.   :50.00   Max.   :50.00  
##                                                                  
##       CRIM                ZN             INDUS           CHAS          
##  Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46   Length:506        
##  1st Qu.: 0.08205   1st Qu.:  0.00   1st Qu.: 5.19   Class :character  
##  Median : 0.25651   Median :  0.00   Median : 9.69   Mode  :character  
##  Mean   : 3.61352   Mean   : 11.36   Mean   :11.14                     
##  3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10                     
##  Max.   :88.97620   Max.   :100.00   Max.   :27.74                     
##                                                                        
##       NOX               RM             AGE              DIS        
##  Min.   :0.3850   Min.   :3.561   Min.   :  2.90   Min.   : 1.130  
##  1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02   1st Qu.: 2.100  
##  Median :0.5380   Median :6.208   Median : 77.50   Median : 3.207  
##  Mean   :0.5547   Mean   :6.285   Mean   : 68.57   Mean   : 3.795  
##  3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08   3rd Qu.: 5.188  
##  Max.   :0.8710   Max.   :8.780   Max.   :100.00   Max.   :12.127  
##                                                                    
##       RAD              TAX           PTRATIO            B         
##  Min.   : 1.000   Min.   :187.0   Min.   :12.60   Min.   :  0.32  
##  1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40   1st Qu.:375.38  
##  Median : 5.000   Median :330.0   Median :19.05   Median :391.44  
##  Mean   : 9.549   Mean   :408.2   Mean   :18.46   Mean   :356.67  
##  3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20   3rd Qu.:396.23  
##  Max.   :24.000   Max.   :711.0   Max.   :22.00   Max.   :396.90  
##                                                                   
##      LSTAT           units             cu5k            c5_7_5      
##  Min.   : 1.73   Min.   :   5.0   Min.   : 0.000   Min.   : 0.000  
##  1st Qu.: 6.95   1st Qu.: 115.0   1st Qu.: 0.000   1st Qu.: 1.000  
##  Median :11.36   Median : 511.5   Median : 1.000   Median : 3.000  
##  Mean   :12.65   Mean   : 680.8   Mean   : 2.921   Mean   : 5.534  
##  3rd Qu.:16.95   3rd Qu.:1152.0   3rd Qu.: 4.000   3rd Qu.: 7.000  
##  Max.   :37.97   Max.   :3031.0   Max.   :35.000   Max.   :70.000  
##                                                                    
##     C7_5_10            C10_15           C15_20          C20_25     
##  Min.   :  0.000   Min.   :  0.00   Min.   :  0.0   Min.   :  0.0  
##  1st Qu.:  2.000   1st Qu.: 14.00   1st Qu.: 19.0   1st Qu.: 13.0  
##  Median :  6.000   Median : 33.00   Median : 85.0   Median :101.5  
##  Mean   :  9.984   Mean   : 55.41   Mean   :141.8   Mean   :166.2  
##  3rd Qu.: 13.000   3rd Qu.: 75.50   3rd Qu.:210.0   3rd Qu.:281.8  
##  Max.   :121.000   Max.   :520.00   Max.   :937.0   Max.   :723.0  
##                                                                    
##      C25_35           C35_50           co50k            median     
##  Min.   :   0.0   Min.   :  0.00   Min.   :  0.00   Min.   : 5600  
##  1st Qu.:   7.0   1st Qu.:  1.00   1st Qu.:  0.00   1st Qu.:16800  
##  Median :  95.0   Median : 17.00   Median :  3.00   Median :21000  
##  Mean   : 170.8   Mean   : 82.74   Mean   : 45.44   Mean   :21749  
##  3rd Qu.: 292.0   3rd Qu.: 97.00   3rd Qu.: 20.00   3rd Qu.:24700  
##  Max.   :1189.0   Max.   :769.00   Max.   :980.00   Max.   :50000  
##                                                     NA's   :17     
##        BB           censored             NOX_ID           POP       
##  Min.   : 0.000   Length:506         Min.   : 1.00   Min.   :  434  
##  1st Qu.: 0.200   Class :character   1st Qu.:18.00   1st Qu.: 3697  
##  Median : 0.500   Mode  :character   Median :44.00   Median : 5105  
##  Mean   : 6.082                      Mean   :41.77   Mean   : 5340  
##  3rd Qu.: 1.600                      3rd Qu.:62.00   3rd Qu.: 6825  
##  Max.   :96.400                      Max.   :96.00   Max.   :15976  
## 

library(tmap)

tmap_mode("view") 

tm_shape(boston.tr) +
  tm_polygons(col = "gray90", border.col = "black") +
  tm_text("TOWN", size = 0.5) +
  tm_compass(position = c("left", "bottom")) +
  tm_title("Distritos de Boston (TOWN)")

Histogramas

El gráfico presenta los histogramas de las variables numéricas del conjunto de datos de Boston. A continuación se destacan algunos hallazgos clave:

• Sesgo a la derecha: Variables como CRIM (tasa de criminalidad), DIS (distancia a centros laborales), y TAX (tasa de impuesto predial) presentan una distribución altamente asimétrica, con la mayoría de las observaciones concentradas en valores bajos y pocos valores extremos altos.

• Variables censuradas: CMEDV y MEDV muestran un claro corte en el valor máximo de 50, lo cual refleja la censura superior impuesta en el conjunto de datos.

• Distribuciones bimodales o discretas: Variables como RAD (accesibilidad a autopistas) y PTRATIO presentan distribuciones con valores repetidos o agrupamientos discretos.

• Normalidad aparente: Algunas variables como RM (número promedio de habitaciones por vivienda) y LAT/LON (coordenadas geográficas) muestran distribuciones más cercanas a una forma normal o simétrica.

• Concentración extrema: ZN (zonificación para terrenos residenciales grandes) y B (proporción transformada de población negra) presentan una alta concentración de observaciones en un valor específico, lo que indica poca variabilidad en muchas zonas.

Dicho análisis preliminar sugiere que será necesario transformar algunas variables previo a usarlas en modelos estadísticos.

boston.c.num <- boston.c %>% select(where(is.numeric))

boston.c.num %>% 
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor") %>%
  ggplot(aes(x = valor)) +
  geom_histogram(bins = 30, fill = "steelblue", color = "white") +
  facet_wrap(~variable, scales = "free") +
  labs(title = "Histogramas de Variables Numéricas")

Matriz de Correlación

Relaciones con el valor de la vivienda (MEDV y CMEDV):

• RM (habitaciones por vivienda) tiene una fuerte correlación positiva con MEDV (0.70), lo que sugiere que más habitaciones se asocian con mayor valor de vivienda.

• LSTAT (población de menor nivel socioeconómico) tiene una fuerte correlación negativa con MEDV (-0.74), indicando que zonas con mayor proporción de personas de bajos ingresos tienden a tener viviendas más baratas.

También hay una correlación negativa notable con NOX (-0.43) y PTRATIO (-0.51), lo que sugiere que altos niveles de contaminación y una mayor proporción alumno/profesor se relacionan con menores valores de vivienda.

Contaminación y urbanización:

• NOX está fuertemente correlacionado con INDUS (0.76) y AGE (0.73), lo cual tiene sentido, ya que las zonas más industriales y con construcciones más antiguas suelen tener más contaminación.

También está negativamente correlacionado con DIS (-0.77), lo que indica que zonas más cercanas a los centros urbanos (menor distancia) presentan mayor contaminación.

Infraestructura y zonificación:

• RAD y TAX tienen una correlación muy alta (0.91), lo que sugiere que las áreas con mejor acceso a autopistas tienden a tener impuestos prediales más altos.

• ZN tiene correlaciones negativas con variables como INDUS, NOX y CRIM, lo que sugiere que las áreas con más zonificación residencial son menos industriales, menos contaminadas y más seguras.

library(corrplot)

mat_cor <- cor(boston.c.num, use = "complete.obs")

corrplot(mat_cor, method = "color", type = "upper", 
         col = colorRampPalette(c("blue", "white", "red"))(200),
         tl.col = "black", tl.srt = 45, 
         title = "Matriz de Correlaciones", mar = c(0,0,1,0))

corrplot(mat_cor, method = "number", type = "upper", 
         tl.col = "black", tl.srt = 45,
         title = "Matriz de Correlaciones (con valores)", mar = c(0,0,1,0))

Correlación con variable CRIME

En este caso, la variable CRIME será nuestra variable y, ante lo cual, se determina que tiene una correlación positiva moderada con TAX, RAD y PTRATIO, y negativa con ZN y DIS, lo que indica que áreas con mayor infraestructura o menor planificación residencial tienden a tener mayor crimen.

mat_cor <- cor(boston.c.num, use = "complete.obs")

cor_crim <- mat_cor["CRIM", ]

sort(cor_crim, decreasing = TRUE)

##        CRIM         RAD         TAX       LSTAT      TOWNNO         NOX 
##  1.00000000  0.62550515  0.58276431  0.45562148  0.44791970  0.42097171 
##       INDUS         AGE     PTRATIO         LON         LAT          ZN 
##  0.40658341  0.35273425  0.28994558  0.06510061 -0.08429296 -0.20046922 
##          RM         DIS           B        MEDV       CMEDV       TRACT 
## -0.21924670 -0.37967009 -0.38506394 -0.38830461 -0.38958244 -0.54716534

corrplot(mat_cor["CRIM", , drop = FALSE], is.corr = FALSE,
         method = "number", tl.col = "black", title = "Correlación de CRIM con otras variables",
         mar = c(0,0,2,0))

Boxplots

Variables con alta dispersión y outliers notables:

• CRIM (criminalidad), DIS (distancia a centros laborales) y B (proporción de personas negras) muestran muchos valores atípicos (puntos rojos), lo que sugiere una gran variabilidad entre zonas.

• ZN, RAD, y INDUS también presentan concentración de valores extremos, lo cual es esperable dado que son variables relacionadas con zonificación o infraestructura urbana.

Variables fuertemente sesgadas:

• ZN y CRIM parecen tener una asimetría positiva (cola larga a la derecha), indicando que la mayoría de los tramos censales tienen valores bajos, pero unos pocos tienen valores muy altos.

• B está fuertemente sesgada con valores máximos agrupados cerca del tope.

Variables más concentradas (baja dispersión y pocos outliers):

• RM (habitaciones) y PTRATIO (alumnos/profesor) tienen una distribución bastante estrecha, lo que indica que esas variables son más homogéneas entre localidades.

Impuestos (TAX) y tramos censales (TRACT):

• TAX tiene una distribución bastante amplia, lo cual sugiere diferencias importantes en política fiscal entre localidades.

• TRACT muestra amplísima variación, pero como es un identificador, esta dispersión no implica un patrón útil directamente, sino que debe tratarse como variable categórica o de localización.

Valor de vivienda (MEDV y CMEDV):

• Ambas presentan outliers por arriba, lo que sugiere presencia de zonas con viviendas significativamente más caras que el promedio, aunque la mayoría de los valores están agrupados en un rango estrecho.

library(dplyr)
library(tidyr)
library(ggplot2)

boston.c.num %>%
  pivot_longer(cols = everything(), names_to = "variable", values_to = "valor") %>%
  ggplot(aes(x = "", y = valor)) +  # x vacío porque no hay categorías internas
  geom_boxplot(fill = "lightblue", color = "darkblue", outlier.color = "red", outlier.shape = 16) +
  facet_wrap(~ variable, scales = "free_y") +  # facetas por variable con escala libre en Y
  labs(title = "Boxplots por Variable Numérica",
       x = NULL,
       y = "Valor") +
  theme_minimal() +
  theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())

Análisis Exploratorio de Datos Espaciales (ESDA: Exploratory Spatial Data Analysis)

Matriz de Conectividad Espacial

Se crea matriz de conectividad espacial o lista de vecinos entre polígonos para lo cual se usa el parámetro de QUEEN para definir la vecindad, en dónde se tomarán en cuaenta si dos polígonos son vecinos o, si comparten un borde o un punto (hace la cercanía más inclusiva).

Con base en el summary obtenemos lo siguiente:

• Hay 506 polígonos (506 localidades en Boston).

• Hay 2910 conexiones de vecindad entre ellos.

• Proporción de pares vecinos con relación al total posible = 1.14

• Cada estado tiene, en promedio, 5.75 vecinos

• Sólo tres localidades tienen 1 vecino: Boston North End, Boston Beacon Hill y Bedford

• Sólo 1 localidad está conectada con 15 vecinos: Boston Mattapan

boston_nb <- poly2nb(boston.trSP, queen = TRUE)
summary(boston_nb)

## Neighbour list object:
## Number of regions: 506 
## Number of nonzero links: 2910 
## Percentage nonzero weights: 1.136559 
## Average number of links: 5.750988 
## Link number distribution:
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  15 
##   3   9  28  81 107 120  87  40  22   5   2   1   1 
## 3 least connected regions:
## 18 51 345 with 1 link
## 1 most connected region:
## 112 with 15 links

Conexiones entre areas

boston.trSP$TOWN[112]   # Area más conectada (con 15 vecinos)

## [1] "Boston Mattapan"

boston.trSP$TOWN[c(18, 15, 345)]  # Areas menos conectadas (con 1 vecino)

## [1] "Boston North End"   "Boston Beacon Hill" "Bedford"

# Se muestran las 9 conexiones del area con mayor número de conexiones, es decir, de Boston Mattapan
boston_nb[[112]]

##  [1]  69  75  77  78 100 101 110 111 113 114 115 117 118 119 120

# Listado de localidades conectados con Boston Mattapan
boston.trSP$TOWN[boston_nb[[112]]]

##  [1] "Boston Roxbury"      "Boston Roxbury"      "Boston Roxbury"     
##  [4] "Boston Savin Hill"   "Boston Savin Hill"   "Boston Dorchester"  
##  [7] "Boston Dorchester"   "Boston Dorchester"   "Boston Mattapan"    
## [10] "Boston Mattapan"     "Boston Mattapan"     "Boston Mattapan"    
## [13] "Boston Forest Hills" "Boston Forest Hills" "Boston Forest Hills"

Visualización Gráfica de la matriz de conectividad (Vecindades estilo ‘Reina’)

boston_listw <- nb2listw(boston_nb, style = "W", zero.policy = TRUE)

centroids <- coordinates(boston.trSP)
plot(boston.trSP, border = "blue", axes = FALSE, las = 1, main = "Boston Tracts - Queen Contiguity")
plot(boston.trSP, col = "lightgrey", border = grey(0.9), add = TRUE)
plot(boston_listw, coords = centroids, pch = 19, cex = 0.1, col = "red", add = TRUE)

CRIME

boston.tr_clean <- na.omit(boston.tr)
boston.trSP<-as(boston.tr, "Spatial")
boston_nb<-poly2nb(boston.trSP)
mapview(boston.trSP, zcol="CRIM", col.regions = viridisLite::magma(20))

Elección de variables explicativas

Con base en la combinación de variables realizada por la regresión stepwise, la matriz de correlación con CRIM y la significancia de cada variable según en random forest, se determina que las principales variables explicativas a utilizar y que conviene rezagar espacialmente son:

MEDV:

• Tiene una correlación con CRIM de -0.39, refleja que el valor de vivienda más alto reduce crimen (relación negativa).

• A mayor valor de la vivienda, menor crimen. Resultado esperado y significativo. Relación clara y altamente significativa (p < 0.001).

• Su alta importancia como predictor indica que el valor mediano de las viviendas en un área está fuertemente asociado con las otras características consideradas.

DIS

• Tiene correlación de -0.38, lo que puede implicar que a mayor distancia a centros laborales = menos crimen.

• A mayor distancia a centros laborales, menor crimen; además es muy significativo (p < 0.001).

• Su alta importancia en ambas métricas (%IncMSE e IncNodePurity) indica que la cercanía a los centros de empleo es un factor crucial para predecir CRIM.

VARImportance

En resumen, considerando ambas gráficas, algunas de las variables que consistentemente aparecen como importantes son:

• LSTAT: Porcentaje de estatus bajo de la población.

• RM: Número promedio de habitaciones por vivienda.

• DIS: Distancias ponderadas a cinco centros de empleo de Boston.

• RAD: Índice de accesibilidad a carreteras radiales.

• CMEDV / MEDV: Valor mediano de las viviendas ocupadas por sus dueños (podrían ser la misma variable o una transformación).

• B: 1000(Bk−0.63) donde Bk es la proporción de afroamericanos por ciudad.

library(randomForest)

set.seed(123)
modelo_rf <- randomForest(CRIM ~ ., data = boston.c.num, importance = TRUE, ntree = 500)

varImpPlot(modelo_rf)

Regresión Stepwise

El modelo incluye 8 variables explicativas de entre las 12 posibles que se integraron inicialmente. Se eligieron por balancear calidad del ajuste y simplicidad del modelo (criterio AIC).

• RAD (+0.53): A mayor acceso a autopistas, mayor crimen. Relación clara y altamente significativa (p < 0.001).

• LSTAT (+0.11): A mayor % de población vulnerable, mayor crimen. No es estadísticamente significativo (p = 0.11).

• ZN (+0.043): Mayor proporción de zonas residenciales parece asociarse a más crimen.

• DIS (–0.91): Mayor distancia a centros laborales, menor crimen. Muy significativo.

• MEDV (–0.17): A mayor valor de la vivienda, menor crimen. Resultado esperado y significativo.

• NOX (–12.75): Menor calidad del aire (mayor NOX), menor crimen (algo contraintuitivo, aunque es posible colinealidad).

• PTRATIO (–0.31): A mayor ratio alumno-profesor, menor crimen. Marginalmente significativo (p ~ 0.09).

• B (–0.0079): A mayor población de personas negras (B), menor crimen. Sin embargo, esta variable suele estar correlacionada con otras condiciones estructurales.

boston_sf <- st_as_sf(boston.trSP)

modelo_full <- lm(CRIM ~ RAD + TAX + LSTAT + ZN + INDUS + DIS + AGE + MEDV+ NOX + RM + PTRATIO + B, data = boston_sf)

# Stepwise regresión usando AIC
modelo_step <- step(modelo_full, direction = "both")

## Start:  AIC=1896.97
## CRIM ~ RAD + TAX + LSTAT + ZN + INDUS + DIS + AGE + MEDV + NOX + 
##     RM + PTRATIO + B
## 
##           Df Sum of Sq   RSS    AIC
## - AGE      1      0.10 20417 1895.0
## - TAX      1     18.86 20436 1895.4
## - RM       1     21.39 20438 1895.5
## - INDUS    1     28.84 20446 1895.7
## <none>                 20417 1897.0
## - PTRATIO  1     84.55 20501 1897.1
## - LSTAT    1    114.89 20532 1897.8
## - NOX      1    165.79 20582 1899.1
## - B        1    178.18 20595 1899.4
## - ZN       1    238.94 20656 1900.9
## - MEDV     1    483.47 20900 1906.8
## - DIS      1    515.68 20932 1907.6
## - RAD      1   1834.83 22252 1938.5
## 
## Step:  AIC=1894.97
## CRIM ~ RAD + TAX + LSTAT + ZN + INDUS + DIS + MEDV + NOX + RM + 
##     PTRATIO + B
## 
##           Df Sum of Sq   RSS    AIC
## - TAX      1     18.81 20436 1893.4
## - RM       1     22.76 20440 1893.5
## - INDUS    1     28.82 20446 1893.7
## <none>                 20417 1895.0
## - PTRATIO  1     84.57 20501 1895.1
## - LSTAT    1    129.63 20546 1896.2
## + AGE      1      0.10 20417 1897.0
## - NOX      1    175.96 20593 1897.3
## - B        1    178.37 20595 1897.4
## - ZN       1    241.26 20658 1898.9
## - MEDV     1    483.38 20900 1904.8
## - DIS      1    563.37 20980 1906.8
## - RAD      1   1842.82 22260 1936.7
## 
## Step:  AIC=1893.44
## CRIM ~ RAD + LSTAT + ZN + INDUS + DIS + MEDV + NOX + RM + PTRATIO + 
##     B
## 
##           Df Sum of Sq   RSS    AIC
## - RM       1      23.0 20459 1892.0
## - INDUS    1      64.4 20500 1893.0
## <none>                 20436 1893.4
## - PTRATIO  1      87.4 20523 1893.6
## - LSTAT    1     137.9 20574 1894.8
## + TAX      1      18.8 20417 1895.0
## + AGE      1       0.0 20436 1895.4
## - B        1     178.1 20614 1895.8
## - NOX      1     181.9 20617 1895.9
## - ZN       1     222.9 20658 1896.9
## - MEDV     1     465.3 20901 1902.8
## - DIS      1     556.9 20992 1905.0
## - RAD      1    4693.4 25129 1996.0
## 
## Step:  AIC=1892.01
## CRIM ~ RAD + LSTAT + ZN + INDUS + DIS + MEDV + NOX + PTRATIO + 
##     B
## 
##           Df Sum of Sq   RSS    AIC
## - INDUS    1      74.0 20533 1891.8
## <none>                 20459 1892.0
## - PTRATIO  1      88.2 20547 1892.2
## - LSTAT    1     118.9 20577 1892.9
## + RM       1      23.0 20436 1893.4
## + TAX      1      19.1 20440 1893.5
## + AGE      1       1.2 20457 1894.0
## - NOX      1     176.9 20636 1894.4
## - B        1     202.4 20661 1895.0
## - ZN       1     233.9 20693 1895.8
## - MEDV     1     458.7 20917 1901.2
## - DIS      1     572.2 21031 1904.0
## - RAD      1    4811.3 25270 1996.9
## 
## Step:  AIC=1891.83
## CRIM ~ RAD + LSTAT + ZN + DIS + MEDV + NOX + PTRATIO + B
## 
##           Df Sum of Sq   RSS    AIC
## <none>                 20533 1891.8
## + INDUS    1      74.0 20459 1892.0
## + TAX      1      58.8 20474 1892.4
## - LSTAT    1     104.7 20637 1892.4
## - PTRATIO  1     119.0 20652 1892.8
## + RM       1      32.6 20500 1893.0
## + AGE      1       1.3 20531 1893.8
## - B        1     198.4 20731 1894.7
## - ZN       1     239.6 20772 1895.7
## - NOX      1     296.6 20829 1897.1
## - MEDV     1     430.2 20963 1900.3
## - DIS      1     507.8 21040 1902.2
## - RAD      1    4739.5 25272 1994.9

summary(modelo_step)

## 
## Call:
## lm(formula = CRIM ~ RAD + LSTAT + ZN + DIS + MEDV + NOX + PTRATIO + 
##     B, data = boston_sf)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -9.860 -2.102 -0.363  0.895 75.702 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  19.683128   6.086010   3.234 0.001301 ** 
## RAD           0.532617   0.049727  10.711  < 2e-16 ***
## LSTAT         0.110173   0.069219   1.592 0.112097    
## ZN            0.043293   0.017977   2.408 0.016394 *  
## DIS          -0.918318   0.261932  -3.506 0.000496 ***
## MEDV         -0.174207   0.053988  -3.227 0.001334 ** 
## NOX         -12.753708   4.760157  -2.679 0.007623 ** 
## PTRATIO      -0.310541   0.182941  -1.697 0.090229 .  
## B            -0.007922   0.003615  -2.191 0.028897 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.428 on 497 degrees of freedom
## Multiple R-squared:  0.4505, Adjusted R-squared:  0.4416 
## F-statistic: 50.92 on 8 and 497 DF,  p-value: < 2.2e-16

MEDV y DIS

library(tmap)
# Calcular rezago espacial de MEDV
boston_sf$sp_lag_MEDV <- lag.listw(boston_listw, boston_sf$MEDV, zero.policy = TRUE)

# Calcular rangos comunes para una escala continua compartida
medv_range <- range(c(boston_sf$MEDV, boston_sf$sp_lag_MEDV), na.rm = TRUE)
medv_breaks <- seq(medv_range[1], medv_range[2], length.out = 10)

# Mapa MEDV original
medv_map <- tm_shape(boston_sf) +  
  tm_polygons(
    col = "MEDV",
    palette = "YlGnBu",
    style = "fixed",
    breaks = medv_breaks,
    title = "MEDV",
    textNA = "",     # Ocultar "Missing"
    showNA = FALSE   # Excluir NA de la leyenda
  ) +
  tm_title("Valor Vivienda (MEDV)") +
  tm_layout(
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Calcular rangos comunes para una escala continua compartida
dis_range <- range(c(boston_sf$DIS, boston_sf$sp_lag_DIS), na.rm = TRUE)
dis_breaks <- seq(dis_range[1], dis_range[2], length.out = 10)

# Mapa DIS original
dis_map <- tm_shape(boston_sf) +  
  tm_polygons(
    col = "DIS",
    palette = "YlOrBr",  # Paleta naranjosa
    style = "fixed",
    breaks = dis_breaks,
    title = "DIS",
    textNA = "",        # Oculta etiqueta "Missing"
    showNA = FALSE      # Excluye NA de la leyenda
  ) +
  tm_title("Distancia a Centros (DIS)") +
  tm_layout(
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

tmap_arrange(medv_map, dis_map, ncol = 1)

Índice Global de Moran

Teóricamente el rango de Moran se muestra entre –1 y +1, en dónde: * +1 Alta autocorrelación espacial positiva (valores similares están cerca entre sí).

• 0 Distribución aleatoria (no hay patrón espacial).

• -1 Alta autocorrelación espacial negativa (valores diferentes están cerca entre sí).

Acorde a lo anterior, se genera codigo para determinar si CRIME tiene autocorrelación espacial, es decir, si las areas ya sea con alto o bajo crimen tienden a estar cerca unos de otros. También se analiza si MEDV y/o DIS, tienem o no autocorrelación espacial, es decir, si las localidades ya sea con alto o bajo nivel tienden a estar cerca unos de otros.

moran.test(boston.trSP$CRIM, boston_listw) # Global Moran's I is 0.527 and statistically significant (p-value < 10%).

## 
##  Moran I test under randomisation
## 
## data:  boston.trSP$CRIM  
## weights: boston_listw    
## 
## Moran I statistic standard deviate = 20.398, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.5274633075     -0.0019801980      0.0006736983

moran.test(boston.trSP$MEDV, boston_listw) # Global Moran's I is 0.627 and statistically significant (p-value < 10%).

## 
##  Moran I test under randomisation
## 
## data:  boston.trSP$MEDV  
## weights: boston_listw    
## 
## Moran I statistic standard deviate = 23.35, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.6266753872     -0.0019801980      0.0007248686

moran.test(boston.trSP$DIS, boston_listw) # Global Moran's I is 0.956 and statistically significant (p-value < 10%).

## 
##  Moran I test under randomisation
## 
## data:  boston.trSP$DIS  
## weights: boston_listw    
## 
## Moran I statistic standard deviate = 35.56, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.956380200      -0.001980198       0.000726315

Presentación de resultados

Crimen (CRIM)
Se determina que sí existe un patrón espacial significativo en la distribución del crimen (CRIM) en Boston. Las áreas con tasas similares de crimen están espacialmente agrupadas. El crimen no está distribuida aleatoriamente; hay agrupamientos claros.

• Moran’s I = 0.527 → Hay una autocorrelación espacial positiva moderada en los valores de criminalidad, es decir, las zonas con alta (o baja) criminalidad tienden a estar cerca unas de otras, en lugar de distribuirse aleatoriamente en el espacio.

• P-value = 2.2e-16 → El valor es extremadamente significativo, es decir, la probabilidad de observar un valor de Moran’s I tan alto por puro azar es virtualmente cero.

Valor de vivienda (MEDV)
Existe un patrón espacial fuerte. El valor de vivienda no está distribuida aleatoriamente ya que, existen hay agrupamientos claros.

• Moran’s I = 0.627 → Existe una fuerte autocorrelación espacial positiva en los valores de la vivienda. Por ende, los vecindarios con viviendas de alto o bajo valor se agrupan geográficamente, lo cual es un patrón esperado en la realidad urbana (zonas caras vs. zonas más accesibles).

• P-value = 2.2e-16 → Sí es estadísticamente significativo (p<0.10).

Distancia a centros (DIS)
Las distancias se comportan como una variable espacialmente estructurada. La distancia a centros no está distribuida aleatoriamente ya que, existen hay agrupamientos claros.

• Moran’s I = 0.956 → Hay una autocorrelación espacial extremadamente fuerte. Las distancias a los centros de empleo están muy agrupadas espacialmente, es decir, hay zonas muy cercanas o muy lejanas a los centros laborales que tienden a estar juntas.

• P-value = 2.2e-16 → Sí es estadísticamente significativo (p<0.10).

table <- data.frame(Variable = c("CRIM", "MEDV", "DIS"), GM = c(0.527, 0.627, 0.956), Significance = c("***", "***", "***"))
table

##   Variable    GM Significance
## 1     CRIM 0.527          ***
## 2     MEDV 0.627          ***
## 3      DIS 0.956          ***

Análisis de Autocorrelación Espacial

Mapa de Lag Espacial CRIM

Se realiza cálculo de lags de acuerdo con la matriz de conectividad anteriormente creada

# Cargar librerías necesarias
library(spdep)
library(tmap)
library(sf)

# Si boston.trSP es de tipo SpatialPolygonsDataFrame, conviértelo a sf
boston_sf <- st_as_sf(boston.trSP)

# Crear objeto de vecinos (Queen contiguity)
boston_nb <- poly2nb(boston_sf, queen = TRUE)
boston_listw <- nb2listw(boston_nb, style = "W", zero.policy = TRUE)

# Calcular rezago espacial de CRIM
boston_sf$sp_lag_CRIM <- lag.listw(boston_listw, boston_sf$CRIM, zero.policy = TRUE)

Se realiza cálculo de rangos comunes y grafico de mapa de CRIM original con mapa de CRIM-lag

# Calcular rangos comunes y ajustar breaks
crim_range <- range(c(boston_sf$CRIM, boston_sf$sp_lag_CRIM), na.rm = TRUE)
crim_breaks <- seq(crim_range[1], crim_range[2], length.out = 10)

# Mapa CRIM original
crim_map <- tm_shape(boston_sf) +  
  tm_polygons(
    col = "CRIM",
    palette = "Purples",
    style = "fixed",
    breaks = crim_breaks,
    title = "CRIM",
    textNA = "",     # Oculta etiqueta "Missing"
    showNA = FALSE   # Excluye NA de la leyenda
  ) +
  tm_title("Tasa de Crimen (CRIM)") +
  tm_layout(
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mapa CRIM con rezago espacial
crim_lag_map <- tm_shape(boston_sf) + 
  tm_polygons(
    col = "sp_lag_CRIM",
    palette = "Purples",
    style = "fixed",
    breaks = crim_breaks,
    title = "CRIM (Lag-1)",
    textNA = "",
    showNA = FALSE
  ) +
  tm_title("Tasa de Crimen Rezagada (Lag-1)") +
  tm_layout(
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mostrar los mapas
tmap_arrange(crim_map, crim_lag_map, ncol = 1)

Mapa de Lag Espacial MEDV

# Calcular rezago espacial de MEDV
boston_sf$sp_lag_MEDV <- lag.listw(boston_listw, boston_sf$MEDV, zero.policy = TRUE)

# Calcular rangos comunes para una escala continua compartida
medv_range <- range(c(boston_sf$MEDV, boston_sf$sp_lag_MEDV), na.rm = TRUE)
medv_breaks <- seq(medv_range[1], medv_range[2], length.out = 10)

# Mapa MEDV original
medv_map <- tm_shape(boston_sf) +  
  tm_polygons(
    col = "MEDV",
    palette = "YlGnBu",
    style = "fixed",
    breaks = medv_breaks,
    title = "MEDV",
    textNA = "",     # Ocultar "Missing"
    showNA = FALSE   # Excluir NA de la leyenda
  ) +
  tm_title("Valor Vivienda (MEDV)") +
  tm_layout(
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mapa MEDV con rezago espacial
medv_lag_map <- tm_shape(boston_sf) + 
  tm_polygons(
    col = "sp_lag_MEDV",
    palette = "YlGnBu",
    style = "fixed",
    breaks = medv_breaks,
    title = "MEDV (Lag-1)",
    textNA = "",
    showNA = FALSE
  ) +
  tm_title("Valor Vivienda Rezagado (Lag-1)") +
  tm_layout(
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mostrar mapas lado a lado
tmap_arrange(medv_map, medv_lag_map, ncol = 1)

Mapa de Lag Espacial DIS

# Calcular rezago espacial de DIS
boston_sf$sp_lag_DIS <- lag.listw(boston_listw, boston_sf$DIS, zero.policy = TRUE)

# Calcular rangos comunes para una escala continua compartida
dis_range <- range(c(boston_sf$DIS, boston_sf$sp_lag_DIS), na.rm = TRUE)
dis_breaks <- seq(dis_range[1], dis_range[2], length.out = 10)

# Mapa DIS original
dis_map <- tm_shape(boston_sf) +  
  tm_polygons(
    col = "DIS",
    palette = "YlOrBr",  # Paleta naranjosa
    style = "fixed",
    breaks = dis_breaks,
    title = "DIS",
    textNA = "",        # Oculta etiqueta "Missing"
    showNA = FALSE      # Excluye NA de la leyenda
  ) +
  tm_title("Distancia a Centros (DIS)") +
  tm_layout(
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mapa DIS con rezago espacial
dis_lag_map <- tm_shape(boston_sf) + 
  tm_polygons(
    col = "sp_lag_DIS",
    palette = "YlOrBr",
    style = "fixed",
    breaks = dis_breaks,
    title = "DIS (Lag-1)",
    textNA = "",
    showNA = FALSE
  ) +
  tm_title("Distancia a Centros Rezagada (Lag-1)") +
  tm_layout(
    title.position = c("right", "top"),
    legend.position = c("left", "bottom"),
    title.size = 1
  )

# Mostrar mapas lado a lado
tmap_arrange(dis_map, dis_lag_map, ncol = 1)

Scatterplot Moran

Los tres gráficos muestran que CRIM, MEDV y DIS presentan patrones espaciales estructurados, lo que tiene implicaciones clave para el modelado espacial:

• La criminalidad no es aleatoria en el espacio. Está influida por lo que ocurre en las zonas vecinas. Por lo tanto, cualquier modelo que explique el crimen debe incorporar dependencia espacial.

• Variables como MEDV y DIS están también espacialmente autocorrelacionadas, lo cual podría introducir efectos indirectos si no se controla adecuadamente.

• MEDV podría tener un efecto negativo sobre CRIM, dado que barrios con viviendas más caras suelen ser más seguros.

• Asimismo, DIS puede estar relacionado con menos crimen en áreas más lejanas (menos densas o más suburbanas), aunque esto debe verificarse en el modelo.

CRIM

• Hay una clara autocorrelación espacial positiva, donde barrios con alta criminalidad están rodeados de barrios similares.

• Se destacan zonas como Boston Roxbury, South Boston y Charlestown, que muestran valores extremos.

• Lo anterior apoya el valor de Moran’s I = 0.527, evidenciando agrupamientos espaciales no aleatorios.

# Ajuste de regresión
M1 <- lm(sp_lag_CRIM ~ CRIM, boston_sf)

# Calcular los residuos absolutos
boston_sf$residuals <- abs(resid(M1))

# Identificar outliers
n_outliers <- 5
top_outliers <- order(boston_sf$residuals, decreasing = TRUE)[1:n_outliers]

# Promedios (para líneas divisorias y cuadrantes)
mean_x <- mean(boston_sf$CRIM)
mean_y <- mean(boston_sf$sp_lag_CRIM)

# Graficar scatterplot de Moran
plot(sp_lag_CRIM ~ CRIM, boston_sf, pch=21, asp=1, las=1, 
     col = "grey40", bg="grey80", main="Tasa de Crimen vs Rezago Espacial",
     xlab = "CRIM", ylab = "Rezago Espacial de CRIM")

abline(M1, col="blue")  # Línea de regresión
abline(v = mean_x, lty=3, col = "grey80")  # Línea vertical media CRIM
abline(h = mean_y, lty=3, col = "grey80")  # Línea horizontal media lag

# Resaltar los outliers
points(boston_sf$CRIM[top_outliers], 
       boston_sf$sp_lag_CRIM[top_outliers], 
       pch=21, bg="red", col="black", cex=1.4)

# Etiquetas de outliers
text(boston_sf$CRIM[top_outliers], 
     boston_sf$sp_lag_CRIM[top_outliers], 
     labels = boston_sf$TOWN[top_outliers], 
     pos = 3, cex = 0.8, col = "black")

MEDV

• Existe una fuerte correlación positiva entre el valor de una vivienda y el de sus vecinas.

• Zonas como Brookline y Beacon Hill muestran altos valores, mientras que Boston North End y Waltham destacan con valores más bajos.

• Confirma el Moran’s I = 0.627, lo que indica una estructura espacial clara en el valor de las viviendas.

M2 <- lm(sp_lag_MEDV ~ MEDV, boston_sf)

# Calcular los residuos absolutos (distancia vertical a la línea de regresión)
boston_sf$residuals <- abs(resid(M2))

# Identificar los n puntos más lejanos (por ejemplo, los 5 mayores residuos)
n_outliers <- 5
top_outliers <- order(boston_sf$residuals, decreasing = TRUE)[1:n_outliers]

# Graficar los datos (Moran scatterplot) + línea de regresión
plot(sp_lag_MEDV ~ MEDV, boston_sf, pch=21, asp=1, las=1, 
     col = "grey40", bg="grey80", main="Tasa de Valor Vivienda vs Rezago Espacial")
abline(M1, col="blue") # Línea de regresión
abline(v = mean(boston_sf$MEDV), lty=3, col = "grey80")
abline(h = mean(boston_sf$sp_lag_MEDV), lty=3, col = "grey80")

# Resaltar los puntos más lejanos (outliers)
points(boston_sf$MEDV[top_outliers], 
       boston_sf$sp_lag_MEDV[top_outliers], 
       pch=21, bg="red", col="black", cex=1.4)

# Agregar etiquetas con el nombre de la localidad
text(boston_sf$MEDV[top_outliers], 
     boston_sf$sp_lag_MEDV[top_outliers], 
     labels = boston_sf$TOWN[top_outliers], 
     pos = 3, cex = 0.8, col = "black")

DIS

• Hay una autocorrelación espacial extremadamente alta (Moran’s I = 0.956). Las zonas con similares distancias a centros laborales están muy agrupadas.

• Ejemplos: Lexington, Duxbury y Garland se agrupan por sus mayores distancias.

M3 <- lm(sp_lag_DIS ~ DIS, boston_sf)

# Calcular los residuos absolutos (distancia vertical a la línea de regresión)
boston_sf$residuals <- abs(resid(M3))

# Identificar los n puntos más lejanos (por ejemplo, los 5 mayores residuos)
n_outliers <- 5
top_outliers <- order(boston_sf$residuals, decreasing = TRUE)[1:n_outliers]

# Graficar los datos (Moran scatterplot) + línea de regresión
plot(sp_lag_DIS ~ DIS, boston_sf, pch=21, asp=1, las=1, 
     col = "grey40", bg="grey80", main="Tasa de Distancia de Centros vs Rezago Espacial")
abline(M1, col="blue") # Línea de regresión
abline(v = mean(boston_sf$DIS), lty=3, col = "grey80")
abline(h = mean(boston_sf$sp_lag_DIS), lty=3, col = "grey80")

# Resaltar los puntos más lejanos (outliers)
points(boston_sf$DIS[top_outliers], 
       boston_sf$sp_lag_DIS[top_outliers], 
       pch=21, bg="red", col="black", cex=1.4)

# Agregar etiquetas con el nombre de la localidad
text(boston_sf$DIS[top_outliers], 
     boston_sf$sp_lag_DIS[top_outliers], 
     labels = boston_sf$TOWN[top_outliers], 
     pos = 3, cex = 0.8, col = "black")

Análisis de clusters

Visualización Espacial de Clusters de CRIM

• High-High: Alta criminalidad rodeada de alta criminalidad. Indica zonas vulnerables, usualmente urbanas, con posibles problemas estructurales.

• Low-High: Baja criminalidad en zonas rodeadas de alta criminalidad. Puede reflejar zonas de resistencia o enclaves seguros en entornos conflictivos.

• Low-Low: Baja criminalidad rodeada de baja criminalidad. Áreas seguras, probablemente con buena infraestructura y calidad de vida.

• Not significant: Sin patrón espacial claro. Áreas heterogéneas o en transición.

swm_a <- queen_weights(boston_sf)  # Matriz de pesos espaciales Queen

lisa_crim <- local_moran(swm_a, boston_sf["CRIM"])

boston_sf$cluster_crim <- as.factor(lisa_crim$GetClusterIndicators())
levels(boston_sf$cluster_crim) <- lisa_crim$GetLabels()

library(ggplot2)

ggplot(data = boston_sf) +
  geom_sf(aes(fill = cluster_crim), color = "white") +
  ggtitle(label = "Tasa de Crimen", subtitle = "Clústeres Espaciales en Boston") +
  theme_minimal() +
  scale_fill_manual(values = c("salmon", "turquoise", "grey", "purple", "green")) + # Colores para los clusters
  theme(legend.title = element_blank(), 
        legend.position = "bottom")

Visualización Espacial de Clusters de MEDV

• High-High: Alto valor de vivienda en zonas rodeadas también de alto valor. Barrios acomodados o suburbanos consolidados.

• Low-High: Bajo valor en zonas de alto valor. Pueden ser focos de gentrificación o rezago en zonas desarrolladas.

• Low-Low: Bajo valor rodeado de bajo valor. Zonas marginadas o de menor desarrollo económico.

• Not significant: Sin patrón claro en el valor de vivienda. Pueden ser zonas con mezcla de clases sociales o desarrollo desigual.

lisa_medv <- local_moran(swm_a, boston_sf["MEDV"])

boston_sf$cluster_medv <- as.factor(lisa_medv$GetClusterIndicators())
levels(boston_sf$cluster_medv) <- lisa_medv$GetLabels()

library(ggplot2)

ggplot(data = boston_sf) +
  geom_sf(aes(fill = cluster_medv), color = "white") +
  ggtitle(label = "Tasa de Valor Vivienda", subtitle = "Clústeres Espaciales en Boston") +
  theme_minimal() +
  scale_fill_manual(values = c("salmon", "turquoise", "grey", "purple", "green")) + # Colores para los clusters
  theme(legend.title = element_blank(), 
        legend.position = "bottom")

Visualización Espacial de Clusters de DIS

• High-High: Alta distancia a centros de empleo en zonas también alejadas. Zonas periféricas con mala accesibilidad.

• Low-Low: Baja distancia a centros de empleo en zonas también bien conectadas. Zonas centrales o privilegiadas para movilidad laboral.

• Not significant: Sin patrón espacial evidente. Pueden combinar accesibilidad variable o estar en zonas en transformación urbana.

lisa_dis <- local_moran(swm_a, boston_sf["DIS"])

boston_sf$cluster_dis <- as.factor(lisa_dis$GetClusterIndicators())
levels(boston_sf$cluster_dis) <- lisa_dis$GetLabels()

library(ggplot2)

ggplot(data = boston_sf) +
  geom_sf(aes(fill = cluster_dis), color = "white") +
  ggtitle(label = "Tasa de Distancia a Centros", subtitle = "Clústeres Espaciales en Boston") +
  theme_minimal() +
  scale_fill_manual(values = c("salmon", "turquoise", "grey", "purple", "green")) + # Colores para los clusters
  theme(legend.title = element_blank(), 
        legend.position = "bottom")

Modelos de Regresión

Modelo de Regresión Lineal Tradicional

model_a <- lm(CRIM ~ RAD + ZN + DIS + MEDV + NOX + B, data = boston_sf)
summary(model_a)

## 
## Call:
## lm(formula = CRIM ~ RAD + ZN + DIS + MEDV + NOX + B, data = boston_sf)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.240  -1.915  -0.376   0.852  75.438 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 14.642639   3.709443   3.947 9.04e-05 ***
## RAD          0.499838   0.044036  11.351  < 2e-16 ***
## ZN           0.053963   0.017305   3.118 0.001923 ** 
## DIS         -0.992811   0.255075  -3.892 0.000113 ***
## MEDV        -0.195990   0.037685  -5.201 2.90e-07 ***
## NOX         -9.238768   4.477580  -2.063 0.039597 *  
## B           -0.008711   0.003612  -2.412 0.016237 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.452 on 499 degrees of freedom
## Multiple R-squared:  0.444,  Adjusted R-squared:  0.4373 
## F-statistic: 66.42 on 6 and 499 DF,  p-value: < 2.2e-16

AIC(model_a)

## [1] 3331.699

Modelo de Regresión Espacial AutoRegresivo (SAR)

model_b <- lagsarlm(CRIM ~ RAD + ZN + DIS + MEDV + NOX + B, data = boston_sf, listw = boston_listw) 
summary(model_b)

## 
## Call:
## lagsarlm(formula = CRIM ~ RAD + ZN + DIS + MEDV + NOX + B, data = boston_sf, 
##     listw = boston_listw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -15.43339  -1.23757  -0.13543   0.60312  77.51277 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  9.1439555  3.4397417  2.6583 0.0078530
## RAD          0.3151860  0.0469326  6.7157 1.872e-11
## ZN           0.0300711  0.0158958  1.8918 0.0585232
## DIS         -0.5420924  0.2372669 -2.2847 0.0223283
## MEDV        -0.1254791  0.0349231 -3.5930 0.0003269
## NOX         -6.9910289  4.0975511 -1.7061 0.0879805
## B           -0.0048688  0.0033311 -1.4616 0.1438538
## 
## Rho: 0.45994, LR test value: 64.289, p-value: 1.1102e-15
## Asymptotic standard error: 0.051889
##     z-value: 8.864, p-value: < 2.22e-16
## Wald statistic: 78.57, p-value: < 2.22e-16
## 
## Log likelihood: -1625.705 for lag model
## ML residual variance (sigma squared): 34.621, (sigma: 5.884)
## Number of observations: 506 
## Number of parameters estimated: 9 
## AIC: 3269.4, (AIC for lm: 3331.7)
## LM test for residual autocorrelation
## test value: 3.161, p-value: 0.075418

AIC(model_b)

## [1] 3269.41

Modelo de Regresión Espacial de Errores (SEM)

model_c <- errorsarlm(CRIM ~ RAD + ZN + DIS + MEDV + NOX + B, data = boston_sf, listw = boston_listw)
summary(model_c)

## 
## Call:errorsarlm(formula = CRIM ~ RAD + ZN + DIS + MEDV + NOX + B, 
##     data = boston_sf, listw = boston_listw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -15.71377  -1.39281  -0.28216   0.67829  74.72508 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  14.4531897   4.9291680  2.9322  0.003366
## RAD           0.5281683   0.0619590  8.5245 < 2.2e-16
## ZN            0.0338710   0.0190154  1.7812  0.074873
## DIS          -0.9547810   0.3515253 -2.7161  0.006605
## MEDV         -0.1443547   0.0443807 -3.2526  0.001143
## NOX         -11.1216131   6.2787627 -1.7713  0.076510
## B            -0.0086773   0.0044375 -1.9555  0.050529
## 
## Lambda: 0.50686, LR test value: 67.679, p-value: 2.2204e-16
## Asymptotic standard error: 0.052503
##     z-value: 9.654, p-value: < 2.22e-16
## Wald statistic: 93.199, p-value: < 2.22e-16
## 
## Log likelihood: -1624.01 for error model
## ML residual variance (sigma squared): 34.029, (sigma: 5.8335)
## Number of observations: 506 
## Number of parameters estimated: 9 
## AIC: 3266, (AIC for lm: 3331.7)

AIC(model_c)

## [1] 3266.019

Modelo de Regresión Espacial Durbin

model_d <- lagsarlm(CRIM ~ RAD + ZN + DIS + MEDV + NOX + B, data = boston_sf, listw = boston_listw, type="mixed") 
summary(model_d)

## 
## Call:
## lagsarlm(formula = CRIM ~ RAD + ZN + DIS + MEDV + NOX + B, data = boston_sf, 
##     listw = boston_listw, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -15.90472  -1.40501  -0.17087   0.76848  73.53893 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)  8.8905e+00  4.3786e+00  2.0304  0.042312
## RAD          5.6078e-01  9.6201e-02  5.8292 5.568e-09
## ZN           1.7655e-02  2.0960e-02  0.8423  0.399608
## DIS         -9.3051e-01  9.0294e-01 -1.0305  0.302755
## MEDV        -1.1644e-01  5.0754e-02 -2.2943  0.021776
## NOX         -2.1876e+01  9.7217e+00 -2.2502  0.024433
## B           -5.7419e-03  5.6323e-03 -1.0195  0.307986
## lag.RAD     -3.3954e-01  1.1709e-01 -2.8998  0.003734
## lag.ZN       3.6048e-02  3.1789e-02  1.1340  0.256810
## lag.DIS      2.7964e-01  9.6850e-01  0.2887  0.772788
## lag.MEDV    -2.7766e-02  7.0705e-02 -0.3927  0.694537
## lag.NOX      1.8741e+01  1.1436e+01  1.6387  0.101270
## lag.B       -5.6835e-04  7.6990e-03 -0.0738  0.941152
## 
## Rho: 0.48978, LR test value: 64.365, p-value: 9.992e-16
## Asymptotic standard error: 0.053295
##     z-value: 9.1899, p-value: < 2.22e-16
## Wald statistic: 84.455, p-value: < 2.22e-16
## 
## Log likelihood: -1620.707 for mixed model
## ML residual variance (sigma squared): 33.723, (sigma: 5.8071)
## Number of observations: 506 
## Number of parameters estimated: 15 
## AIC: 3271.4, (AIC for lm: 3333.8)
## LM test for residual autocorrelation
## test value: 1.026, p-value: 0.31109

AIC(model_d)

## [1] 3271.415

Comparativa de modelos

Con base en los resultados de los cuatro modelos de regresión estimados para explicar la tasa de crimen (CRIM), se determina que:

• Todos los modelos espaciales (2 a 4) tienen mejor desempeño que el OLS (Modelo 1), lo cual indica que existe autocorrelación espacial en los datos y debe ser tomada en cuenta.

• Modelo 3 (SEM) tiene el AIC más bajo, lo que le permite explicar bien con menos complejidad.

• Modelo 4 (SDM) incluye rezagos de las variables explicativas. Tiene la mejor log-likelihood y menor sigma², lo cual sugiere que ajusta mejor los datos, aunque su AIC es ligeramente más alto.

Acorde a lo anterior, se determina que el mejor modelo es el de Regresión Espacial Durbin debido a que:

• Mayor log-likelihood (−1,620.707): El Modelo 4 tiene el valor menos negativo, lo que implica un mejor ajuste general, por ende, es el que “explica” de forma más probable los datos observados dados sus parámetros estimados.

• Menor error residual (σ² = 33.723): El Modelo 4 muestra la menor varianza del residuo entre los modelos, lo cual es indicio de una mayor precisión en las predicciones (comete menos errores al predecir la dependiente).

• Aunque el AIC de Modelo 3 tiene una diferencia de 5.396 con el AIC del Modelo 4. Este último modelo incorpora simultáneamente la dependencia espacial en la variable dependiente (lag de CRIM) y la dependencia espacial en el error (estructura de autocorrelación del error).

• Resultados significativos en Wald y Likelihood Ratio Tests Ambas pruebas (Wald y LR) evalúan la relevancia de incluir componentes espaciales. En el Modelo 4, estos tests son estadísticamente significativos (p < 0.05), lo que indica que los efectos espaciales no solo están presentes, sino que su inclusión mejora significativamente la explicación del modelo.

stargazer(model_a, model_b, model_c, model_d, type = "text", title="Estimated Regression Results")

## 
## Estimated Regression Results
## ====================================================================================
##                                           Dependent variable:                       
##                     ----------------------------------------------------------------
##                                                   CRIM                              
##                               OLS              spatial      spatial      spatial    
##                                             autoregressive   error    autoregressive
##                               (1)                (2)          (3)          (4)      
## ------------------------------------------------------------------------------------
## RAD                        0.500***            0.315***     0.528***     0.561***   
##                             (0.044)            (0.047)      (0.062)      (0.096)    
##                                                                                     
## ZN                         0.054***             0.030*       0.034*       0.018     
##                             (0.017)            (0.016)      (0.019)      (0.021)    
##                                                                                     
## DIS                        -0.993***           -0.542**    -0.955***      -0.931    
##                             (0.255)            (0.237)      (0.352)      (0.903)    
##                                                                                     
## MEDV                       -0.196***          -0.125***    -0.144***     -0.116**   
##                             (0.038)            (0.035)      (0.044)      (0.051)    
##                                                                                     
## NOX                        -9.239**            -6.991*      -11.122*    -21.876**   
##                             (4.478)            (4.098)      (6.279)      (9.722)    
##                                                                                     
## B                          -0.009**             -0.005      -0.009*       -0.006    
##                             (0.004)            (0.003)      (0.004)      (0.006)    
##                                                                                     
## lag.RAD                                                                 -0.340***   
##                                                                          (0.117)    
##                                                                                     
## lag.ZN                                                                    0.036     
##                                                                          (0.032)    
##                                                                                     
## lag.DIS                                                                   0.280     
##                                                                          (0.969)    
##                                                                                     
## lag.MEDV                                                                  -0.028    
##                                                                          (0.071)    
##                                                                                     
## lag.NOX                                                                   18.741    
##                                                                          (11.436)   
##                                                                                     
## lag.B                                                                     -0.001    
##                                                                          (0.008)    
##                                                                                     
## Constant                   14.643***           9.144***    14.453***     8.891**    
##                             (3.709)            (3.440)      (4.929)      (4.379)    
##                                                                                     
## ------------------------------------------------------------------------------------
## Observations                  506                506          506          506      
## R2                           0.444                                                  
## Adjusted R2                  0.437                                                  
## Log Likelihood                                -1,625.705   -1,624.010   -1,620.707  
## sigma2                                          34.621       34.029       33.723    
## Akaike Inf. Crit.                             3,269.410    3,266.019    3,271.415   
## Residual Std. Error    6.452 (df = 499)                                             
## F Statistic         66.418*** (df = 6; 499)                                         
## Wald Test (df = 1)                            78.570***    93.199***    84.455***   
## LR Test (df = 1)                              64.289***    67.679***    64.365***   
## ====================================================================================
## Note:                                                    *p<0.1; **p<0.05; ***p<0.01

Conclusión

El análisis exploratorio de la variable dependiente CRIM (tasa de crimen por distrito) muestra una clara concentración espacial del crimen en la zona central de Boston. Esta área presenta los valores más altos de criminalidad, representados por tonalidades oscuras en los mapas temáticos. Al comparar la distribución espacial de CRIM con su rezago espacial (Lag-1), se observa una gran similitud, lo cual sugiere la presencia de dependencia espacial: las zonas con alto nivel de crimen tienden a estar rodeadas por otras zonas con niveles también elevados.

El análisis LISA (Indicadores Locales de Asociación Espacial) permitió identificar agrupamientos significativos de valores altos y bajos de crimen:

• High-High: zonas con alta criminalidad rodeadas de zonas igualmente altas. Se concentran principalmente en el centro-norte de la ciudad y representan áreas críticas desde la perspectiva de seguridad pública.

• Low-High: zonas con baja criminalidad pero rodeadas de vecindarios con alta criminalidad. Estas pueden representar zonas de riesgo por contagio espacial del crimen.

• Low-Low: zonas con baja criminalidad rodeadas por otras zonas similares, posiblemente asociadas a condiciones socioeconómicas favorables.

• Not Significant: zonas sin patrones espaciales significativos.

El modelo seleccionado fue la Regresión Espacial Durbin, ya que incorpora tanto los efectos locales como los rezagos espaciales de las variables explicativas. Asimismo, mostró un buen desempeño estadístico, con un log-likelihood de –1,620.7 y un criterio de información de Akaike (AIC) relativamente bajo, confirmando su ajuste óptimo. Los resultados más relevantes demuestran que:

• RAD (acceso a autopistas): tiene un efecto directo positivo y significativo (coef. = 0.561), lo cual indica que a mayor accesibilidad a autopistas, mayor es la tasa de crimen. Sin embargo, el efecto espacial de esta variable es negativo y también significativo (lag.RAD = –0.340***), lo cual sugiere que estar rodeado de zonas bien conectadas puede disminuir la criminalidad local, posiblemente por la redistribución de actividades o vigilancia indirecta.

• MEDV (valor medio de las viviendas): tiene un efecto negativo (–0.116**), lo que implica que en vecindarios con mayor valor inmobiliario la criminalidad tiende a ser menor.

• NOX (contaminación por óxidos de nitrógeno): presenta un efecto negativo y significativo (–21.876**), lo cual podría estar relacionado con un entorno físico poco habitable que desincentiva la actividad delictiva o concentra menor densidad poblacional.

• ZN (porcentaje de terrenos residenciales grandes), DIS (distancia a centros laborales) y B (índice racial) pierden significancia estadística en este modelo, lo que indica que sus efectos son absorbidos por otras variables o por sus componentes espaciales.

Con base en lo anterior, se determina que las estrategias de prevención del crimen deben considerar no solo las condiciones internas de cada distrito, sino también su contexto espacial. En otras palabras, las políticas públicas deben diseñarse con una perspectiva regional, reconociendo la interdependencia entre zonas urbanas para lograr intervenciones más efectivas y sostenibles.

Referencias

What is exploratory data analysis? https://www.ibm.com/topics/exploratory-data-analysis

Exploratory Spatial Data Analysis https://cran.r-project.org/web/packages/geostan/vignettes/measuring-sa.html