Regresión Lineal Espacial
Equipo 1
2026-04-21
Ejercicio Regresión Lineal Espacial
El propósito de esta actividad es poner en práctica lo estudiado hasta ahora en los temas de Análisis Exploratorio de Datos Espaciales, así como Modelos de Regresión Lineal Espaciales.
Primeramente aquí se proporcionan las librerías principales que deberás considerar. De todas maneras si deseas agregar cualquier otra librería lo puedes hacer y recomendaría hacerlo de una vez dentro del siguiente chunk de código:
library(foreign)
library(dplyr)
library(spdep)
library(tigris)
library(rgeoda)
library(RColorBrewer)
library(viridis)
library(ggplot2)
library(tmap)
library(sf)
library(sp)
library(spatialreg)
library(stargazer)
library(plm)
library(splm)
library(pspatreg)
library(regclass)
library(mctest)
library(lmtest)
library(spData)
library(mapview)
library(naniar)
library(dlookr)
library(caret)
library(e1071)
library(SparseM)
library(Metrics)
library(randomForest)
library(rpart.plot)
library(insight)
library(jtools)
library(xgboost)
library(DiagrammeR)
library(effects)
library(leaflet)
library(corrplot)
library(spatialreg)
library(stargazer)
library(spdep)A continuación se carga la base de datos a ser analizada: ‘boston’ de la librería ‘spData’
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.21glimpse(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…#?bostonSeguidamente, 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\Alberto JIménez\AppData\Local\R\win-library\4.5\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: NAD27Ejemplo de primer visualización en Mapas utilizando los datos ya cargados y visualizando los valores utilizando una variedad de paletas de colores.
boston.trSP<-as(boston.tr, "Spatial")
boston_nb<-poly2nb(boston.trSP)
mapview(boston.trSP, zcol="MEDV", col.regions = colorRampPalette(c('lightgrey', 'darkblue')))mapview(boston.trSP, zcol="MEDV", col.regions = brewer.pal(7,'Blues'))mapview(boston.trSP, zcol="MEDV", col.regions = viridisLite::magma(20))#?mapview
#?colorRampPalette
#?brewer.pal
#?viridisLa meta esencial del análisis que llevarás a cabo será explicar la variable ‘CRIM’ o ‘CMEDV’ de acuerdo a la paridad del número de su equipo: - Equipos Impares 1, 3, 5 explicarán ‘CRIM’ - Equipos Pares 2 y 4 explicarán ‘CMEDV’ (en este caso deberán de omitir la variable ‘MEDV’ de la base de datos)
Realiza un breve análisis exploratorio de datos tradicional. Cuestiones que se deberán de incluir en esta sección incluyen: - Histogramas - Matriz de Correlaciones - Boxplots
hist(boston.c$CRIM,
main = "Histograma de CRIM",
xlab = "CRIM (Tasa de criminalidad)",
col = "lightblue",
border = "black",
breaks = 30)
boxplot(boston.c$CRIM,
main = "Boxplot de CRIM",
ylab = "CRIM",
col = "lightgreen",
horizontal = TRUE)
# Filtrar solo columnas numéricas
boston_num <- boston.c[, sapply(boston.c, is.numeric)]
# Calcular matriz de correlación
cor_matrix <- cor(boston_num)
# Ver correlaciones de CRIM
cor_matrix["CRIM", ]## TOWNNO TRACT LON LAT MEDV CMEDV
## 0.44791970 -0.54716534 0.06510061 -0.08429296 -0.38830461 -0.38958244
## CRIM ZN INDUS NOX RM AGE
## 1.00000000 -0.20046922 0.40658341 0.42097171 -0.21924670 0.35273425
## DIS RAD TAX PTRATIO B LSTAT
## -0.37967009 0.62550515 0.58276431 0.28994558 -0.38506394 0.45562148corrplot(cor_matrix["CRIM", , drop = FALSE],
method = "color",
tl.cex = 0.8)
boston_num <- boston.c[, sapply(boston.c, is.numeric)]
corrplot(cor_matrix,
method = "color",
type = "upper",
tl.cex = 0.6,
tl.col = "black",
addCoef.col = "black")
#?bostonRealiza un análisis exploratorio de datos espaciales. Los puntos esenciales a generar en esta sección son los siguientes: - Definir claramente la matriz de conectividad espacial - Graficar dicha matriz, utilizar la estructura ‘Reina’ para determinar los vecinos - Realizar mapa coloreado de acuerdo al valor de la variable que se desea explicar y al menos 2 de las variables independientes que se utilizarán - Realizar análisis de clusters: + Incluir mapa de un lag espacial de la variable que se desea explicar junto con el mapa original previamente realizado + Visualización Espacial de Clusters (HotSpots, ColdSpots, NonSignificant y Atípicos: HighLows, LowHighs) - Índice Global de Moran + Calcular el Estadístico del Índice Global de Moran para las variables previamente graficadas (al menos son 3: 1 dependiente, 2 independientes) + Generar sus respectivos gráficos de dispersión
# Crear vecinos tipo Reina
nb_queen <- poly2nb(boston.tr, queen = TRUE)
# Convertir a lista de pesos espaciales
lw_queen <- nb2listw(nb_queen, style = "W")
# Ver resumen
summary(nb_queen)## 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 linksplot(st_geometry(boston.tr), border = "gray")
plot(nb_queen, st_coordinates(st_centroid(boston.tr)),
add = TRUE, col = "red")
tmap_mode("plot")
# CRIM
mapview(boston.trSP,
zcol = "CRIM",
col.regions = colorRampPalette(c("navy", "cyan", "yellow", "red")))# RAD
mapview(boston.trSP,
zcol = "RAD",
col.regions = colorRampPalette(c("purple", "pink", "orange")))# TAX
mapview(boston.trSP,
zcol = "TAX",
col.regions = viridisLite::plasma(20))# Lag espacial de CRIM
boston.tr$lag_CRIM <- lag.listw(lw_queen, boston.tr$CRIM)
# Mapa comparativo
tm_shape(boston.tr) +
tm_fill(c("CRIM", "lag_CRIM"),
palette = "Reds",
title = c("CRIM", "Lag CRIM")) +
tm_facets(ncol = 2) +
tm_borders()
# Moran local
local_moran <- localmoran(boston.tr$CRIM, lw_queen)
# Agregar resultados
boston.tr$Ii <- local_moran[,1]
boston.tr$pvalue <- local_moran[,5]
# Clasificación de clusters
boston.tr$cluster <- "No significativo"
mean_crim <- mean(boston.tr$CRIM)
for(i in 1:nrow(boston.tr)){
if(boston.tr$pvalue[i] < 0.05){
if(boston.tr$CRIM[i] > mean_crim & boston.tr$lag_CRIM[i] > mean_crim){
boston.tr$cluster[i] <- "High-High"
} else if(boston.tr$CRIM[i] < mean_crim & boston.tr$lag_CRIM[i] < mean_crim){
boston.tr$cluster[i] <- "Low-Low"
} else if(boston.tr$CRIM[i] > mean_crim & boston.tr$lag_CRIM[i] < mean_crim){
boston.tr$cluster[i] <- "High-Low"
} else {
boston.tr$cluster[i] <- "Low-High"
}
}
}
# Mapa de clusters
tm_shape(boston.tr) +
tm_fill("cluster",
palette = c("red", "blue", "lightblue", "pink", "grey"),
title = "Clusters LISA") +
tm_borders()
# Moran global
moran.test(boston.tr$CRIM, lw_queen)##
## Moran I test under randomisation
##
## data: boston.tr$CRIM
## weights: lw_queen
##
## 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.0006736983moran.test(boston.tr$RAD, lw_queen)##
## Moran I test under randomisation
##
## data: boston.tr$RAD
## weights: lw_queen
##
## Moran I statistic standard deviate = 31.894, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.8587313616 -0.0019801980 0.0007282609moran.test(boston.tr$TAX, lw_queen)##
## Moran I test under randomisation
##
## data: boston.tr$TAX
## weights: lw_queen
##
## Moran I statistic standard deviate = 30.389, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.818327988 -0.001980198 0.000728656moran.plot(boston.tr$CRIM, lw_queen,
main = "Moran Plot - CRIM")
moran.plot(boston.tr$RAD, lw_queen,
main = "Moran Plot - RAD")
moran.plot(boston.tr$TAX, lw_queen,
main = "Moran Plot - TAX")
Por último deberás de ajustar los Modelos de Regresión estudiados en este Módulo: - Modelo de Regresión Lineal Tradicional - Modelo de Regresión Espacial AutoRegresivo (SAR) - Modelo de Regresión Espacial de Errores (SEM) - Modelo de Regresión Espacial Durbin
# Seleccionar solo variables numéricas
vars_model <- boston.tr[, c(
"CRIM",
"ZN", "INDUS", "NOX", "RM", "AGE",
"DIS", "RAD", "TAX", "PTRATIO", "B", "LSTAT", "MEDV"
)]
# Fórmula automática (CRIM contra todo lo demás)
formula_full <- CRIM ~ ZN + INDUS + NOX + RM + AGE +
DIS + RAD + TAX + PTRATIO +
B + LSTAT + MEDV
formula_full## CRIM ~ ZN + INDUS + NOX + RM + AGE + DIS + RAD + TAX + PTRATIO +
## B + LSTAT + MEDV#Modelo de OLS
model_a <- lm(formula_full, data = vars_model)
summary(model_a)##
## Call:
## lm(formula = formula_full, data = vars_model)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.906 -2.130 -0.314 1.048 75.061
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.713e+01 7.229e+00 2.369 0.018218 *
## ZN 4.497e-02 1.872e-02 2.402 0.016674 *
## INDUS -6.920e-02 8.293e-02 -0.834 0.404433
## NOX -1.053e+01 5.262e+00 -2.001 0.045960 *
## RM 4.400e-01 6.123e-01 0.719 0.472729
## AGE 8.842e-04 1.789e-02 0.049 0.960604
## DIS -9.933e-01 2.815e-01 -3.529 0.000457 ***
## RAD 5.843e-01 8.778e-02 6.656 7.5e-11 ***
## TAX -3.461e-03 5.128e-03 -0.675 0.500046
## PTRATIO -2.660e-01 1.862e-01 -1.429 0.153686
## B -7.611e-03 3.669e-03 -2.074 0.038576 *
## LSTAT 1.260e-01 7.568e-02 1.666 0.096432 .
## MEDV -2.045e-01 5.984e-02 -3.417 0.000686 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.435 on 493 degrees of freedom
## Multiple R-squared: 0.4536, Adjusted R-squared: 0.4403
## F-statistic: 34.1 on 12 and 493 DF, p-value: < 2.2e-16AIC(model_a)## [1] 3334.936#Modelos espaciales
# SAR
model_b <- lagsarlm(formula_full,
data = vars_model,
listw = lw_queen)# SEM
model_c <- errorsarlm(formula_full,
data = vars_model,
listw = lw_queen)# Durbin
model_d <- lagsarlm(formula_full,
data = vars_model,
listw = lw_queen,
type = "mixed")#Comparación de Modelos
stargazer(model_a, model_b, model_c, model_d,
type = "text",
title = "Modelos espaciales - CRIM")##
## Modelos espaciales - CRIM
## =====================================================================================
## Dependent variable:
## -----------------------------------------------------------------
## CRIM
## OLS spatial spatial spatial
## autoregressive error autoregressive
## (1) (2) (3) (4)
## -------------------------------------------------------------------------------------
## ZN 0.045** 0.026 0.024 0.002
## (0.019) (0.017) (0.020) (0.022)
##
## INDUS -0.069 -0.042 -0.075 -0.117
## (0.083) (0.076) (0.101) (0.125)
##
## NOX -10.527** -8.345* -12.025* -19.008*
## (5.262) (4.803) (7.147) (10.550)
##
## RM 0.440 0.115 0.427 0.124
## (0.612) (0.559) (0.654) (0.693)
##
## AGE 0.001 0.009 0.015 0.017
## (0.018) (0.016) (0.020) (0.022)
##
## DIS -0.993*** -0.593** -1.020*** -0.902
## (0.281) (0.260) (0.384) (0.895)
##
## RAD 0.584*** 0.400*** 0.667*** 0.745***
## (0.088) (0.085) (0.108) (0.140)
##
## TAX -0.003 -0.003 -0.005 -0.006
## (0.005) (0.005) (0.006) (0.006)
##
## PTRATIO -0.266 -0.239 -0.374* -0.469*
## (0.186) (0.170) (0.223) (0.270)
##
## B -0.008** -0.004 -0.009* -0.009
## (0.004) (0.003) (0.005) (0.006)
##
## LSTAT 0.126* -0.015 -0.083 -0.183**
## (0.076) (0.069) (0.079) (0.085)
##
## MEDV -0.204*** -0.174*** -0.237*** -0.206***
## (0.060) (0.055) (0.065) (0.070)
##
## lag.ZN 0.042
## (0.034)
##
## lag.INDUS 0.080
## (0.180)
##
## lag.NOX 16.090
## (13.024)
##
## lag.RM -0.063
## (1.108)
##
## lag.AGE -0.035
## (0.033)
##
## lag.DIS 0.373
## (0.977)
##
## lag.RAD -0.512***
## (0.191)
##
## lag.TAX 0.006
## (0.010)
##
## lag.PTRATIO 0.465
## (0.388)
##
## lag.B 0.003
## (0.008)
##
## lag.LSTAT 0.494***
## (0.138)
##
## lag.MEDV 0.193*
## (0.107)
##
## Constant 17.126** 15.343** 23.079*** 2.957
## (7.229) (6.603) (8.297) (11.013)
##
## -------------------------------------------------------------------------------------
## Observations 506 506 506 506
## R2 0.454
## Adjusted R2 0.440
## Log Likelihood -1,623.569 -1,620.062 -1,608.948
## sigma2 34.355 33.304 32.392
## Akaike Inf. Crit. 3,277.137 3,270.124 3,271.897
## Residual Std. Error 6.435 (df = 493)
## F Statistic 34.101*** (df = 12; 493)
## Wald Test (df = 1) 76.800*** 107.953*** 71.581***
## LR Test (df = 1) 59.799*** 66.812*** 53.965***
## =====================================================================================
## Note: *p<0.1; **p<0.05; ***p<0.01Generar breve comparativa de estos modelos y elegir ¿cuál consideran qué es el mejor modelo? y ¿por qué?
Tras evaluar la estructura de los datos, se determinó que el modelo de Error Espacial (SEM) es la herramienta más robusta para explicar la criminalidad, superando al modelo lineal tradicional y a otras variantes espaciales.
- ¿Por qué SEM? Inadecuación del OLS: Las pruebas Wald y Likelihood Ratio (LR) confirmaron una dependencia espacial significativa. Ignorar la ubicación geográfica invalidaría los resultados de una regresión tradicional.
Eficiencia y Precisión (AIC): El modelo SEM obtuvo el menor valor de AIC (3270.124). Esto significa que ofrece el mejor ajuste con la menor complejidad posible.
Rechazo del modelo Durbin (SDM): Aunque el modelo SDM es más complejo, presentó sobreajuste y variables no significativas. Bajo el principio de parsimonia, el SEM resultó superior por ser más interpretable y eficiente.
Similar a las exposiciones de las Actividades 01 y 02, esta se llevará a cabo al iniciar la clase del Martes 21 de Abril y para esto pueden utilizar únicamente el archivo HTML kniteado a partir de su código RMD.
Sin embargo, la entrega final deberá incluir tanto tu código completamente ejecutado en archivo HTML, así como una breve presentación en la que se mencionen claramente los principales patrones o insights detectados en su análisis, así como sus conclusiones. Los archivos se deberán subir en Canvas por equipo a más tardar el día Viernes 24 de Abril en la sección habilitada por el profesor para esto, junto con los archivos correspondientes a las actividades 01 y 02.