ACTIVIDAD 3 – ANÁLISIS DE REGRESIÓN ESPACIAL LOCAL

Preparación previa para realizar la actividad

#INSTALAR LIBRERÍAS
#install.packages("stringi")
#install.packages("spatialreg")
#install.packages("randomForest")
#install.packages(GWmodel)
library(GWmodel)
library(stringi)
library(dplyr)     
library(ggplot2)  
library(readxl)
library(scales)  
library(sf)
library(tmap)
library(tmaptools)
library(spdep)
library(spatialreg)
library(randomForest)
#LEER LOS ARCHIVOS
df_turismo<- read_excel("C:\\Users\\Diego Pérez\\Downloads\\tourism_state_data.xlsx")
df_hoteles<- read_excel("C:\\Users\\Diego Pérez\\Downloads\\turismo_llegadas_extranjeros.xlsx")

# Filtrar solo datos del año 2022
df_turismo <- df_turismo %>%
  filter(year == 2022)

# Renombrar columnas si es necesario
names(df_turismo)[which(names(df_turismo) == "region...17")] <- "region"
names(df_turismo)[which(names(df_turismo) == "region...18")] <- "region_numero"

# Hacer la unión sin modificar los nombres de estado
df_turismo <- df_turismo %>%
  left_join(df_hoteles %>% select(estado, total_visitantes),
            by = c("state" = "estado"))

#Vemos nuestra base de datos
df_turismo <- df_turismo %>%
  mutate(state = ifelse(state == "estado de mexico", "mexico", state))
shapefile <- st_read("C:\\Users\\Diego Pérez\\Downloads\\mx_spatial_data\\mx_spatial_data\\mx_maps\\mx_states\\mexlatlong.shp")
## Reading layer `mexlatlong' from data source 
##   `C:\Users\Diego Pérez\Downloads\mx_spatial_data\mx_spatial_data\mx_maps\mx_states\mexlatlong.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -118.4042 ymin: 14.55055 xmax: -86.73862 ymax: 32.71846
## Geodetic CRS:  WGS 84
names(shapefile)
##  [1] "OBJECTID"   "FIPS_ADMIN" "GMI_ADMIN"  "ADMIN_NAME" "FIPS_CNTRY"
##  [6] "GMI_CNTRY"  "CNTRY_NAME" "POP_ADMIN"  "TYPE_ENG"   "TYPE_LOC"  
## [11] "SQKM"       "SQMI"       "COLOR_MAP"  "Shape_Leng" "Shape_Area"
## [16] "OBJECTID_1" "OBJECTID_2" "longitude"  "latitude"   "geometry"
head(shapefile)
## Simple feature collection with 6 features and 19 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -109.4428 ymin: 22.21486 xmax: -97.1375 ymax: 31.78389
## Geodetic CRS:  WGS 84
##   OBJECTID FIPS_ADMIN GMI_ADMIN ADMIN_NAME FIPS_CNTRY GMI_CNTRY CNTRY_NAME
## 1      888       MX06   MEX-CHH  Chihuahua         MX       MEX     Mexico
## 2      933       MX07   MEX-CDZ   Coahuila         MX       MEX     Mexico
## 3      976       MX19   MEX-NLE Nuevo Leon         MX       MEX     Mexico
## 4      978       MX28   MEX-TML Tamaulipas         MX       MEX     Mexico
## 5      998       MX25   MEX-SIN    Sinaloa         MX       MEX     Mexico
## 6     1004       MX10   MEX-DRN    Durango         MX       MEX     Mexico
##   POP_ADMIN TYPE_ENG TYPE_LOC      SQKM     SQMI COLOR_MAP Shape_Leng
## 1   2656214    State   Estado 247935.02 95727.70        12   22.60928
## 2   2145539    State   Estado 150843.95 58240.87         2   18.99309
## 3   3370912    State   Estado  65173.05 25163.31         3   15.42617
## 4   2272724    State   Estado  79502.24 30695.81        11   18.02314
## 5   2397706    State   Estado  57638.85 22254.36         5   16.46605
## 6   1467826    State   Estado 120674.60 46592.46         4   17.51290
##   Shape_Area OBJECTID_1 OBJECTID_2  longitude latitude
## 1  22.890985          1        888 -106.44431 28.80303
## 2  13.733655          2        933 -102.03356 27.30662
## 3   5.844668          3        976  -99.83125 25.60105
## 4   7.056563          4        978  -98.62181 24.29819
## 5   5.145524          5        998 -107.48280 25.02179
## 6  10.764853          6       1004 -104.92001 24.93519
##                         geometry
## 1 MULTIPOLYGON (((-103.6309 2...
## 2 MULTIPOLYGON (((-102.6669 2...
## 3 MULTIPOLYGON (((-99.7139 27...
## 4 MULTIPOLYGON (((-98.61609 2...
## 5 MULTIPOLYGON (((-108.3942 2...
## 6 MULTIPOLYGON (((-104.3114 2...
# Primero, asegúrate que ambos sean del mismo tipo
shapefile$OBJECTID <- as.numeric(shapefile$OBJECTID)
df_turismo$state_id <- as.numeric(df_turismo$state_id)

# Join por ID
estados_datos <- left_join(shapefile, df_turismo, by = c("OBJECTID" = "state_id"))

# Obtener solo las columnas numéricas sin geometría
datos_numericos <- estados_datos %>%
  st_drop_geometry() %>%
  select(tourism_gdp, crime_rate, college_education, unemployment,
         employment, business_activity, real_wage, pop_density,
         good_governance, ratio_public_investment)

# Estandarizar con scale()
datos_estandarizados <- as.data.frame(scale(datos_numericos))

# Combinar con la geometría y otras variables importantes
estados_datos_std <- estados_datos %>%
  select(total_visitantes, geometry) %>%
  bind_cols(datos_estandarizados)

# Crear vecinos y lista de pesos
vecinos <- poly2nb(estados_datos_std)
pesos <- nb2listw(vecinos, style = "W", zero.policy = TRUE)

Pregunta 1. Describir 3 – 5 diferencias entre la estimación de un modelo de regression global (no especial, SAR, SEM, SDM, y/o SARAR) y un modelo de regresión local (GWR)

  1. Los modelos globales como el SAR, SEM, etc, estiman un solo conjunto de coeficientes por un terrirotio, los modelos locales, tienen una estimación por cada ubicación.

  2. Mientras que los métodos globales estiman pocas parámetros (1 por variable), los métodos locales estiman varias por cada una de las variables que estemos tratando.

  3. Mientras que en los métodos globales, el efecto de las variables es constante (Homogeneidad espacial), para los métodos locales, las variables son diferentes para cada región.

# Asegurarte de que los datos están alineados antes de combinar
stopifnot(nrow(estados_datos) == nrow(datos_estandarizados))

# Combinar datos con cbind (más seguro que bind_cols en este caso)
estados_datos_std <- estados_datos %>%
  select(total_visitantes, geometry)

estados_datos_std <- st_as_sf(cbind(estados_datos_std, datos_estandarizados))

Pregunta 2. Calcular 2 matrices de conectividad que incluya uno de los siguientes criterios: contiguedad (queen, rook, bishop), vecino más próximo (k nearest neighbor), y/o distancia.

# Crear vecinos tipo Queen
vecinos_queen <- poly2nb(estados_datos, queen = TRUE)
# Obtener centroides para trazar líneas
coords <- st_coordinates(st_centroid(estados_datos))
# Mapa de conectividad Queen
plot(st_geometry(estados_datos), border = "gray", main = "Conectividad tipo Queen")
plot(vecinos_queen, coords, add = TRUE, col = "darkgreen", lwd = 1)

# Crear vecinos tipo Rook
vecinos_rook <- poly2nb(estados_datos, queen = FALSE)

# Mapa de conectividad Rook
plot(st_geometry(estados_datos), border = "gray", main = "Conectividad tipo Rook")
plot(vecinos_rook, coords, add = TRUE, col = "purple", lwd = 1)

Pregunta 3. Especificar y estimar 1 modelo de regresión global no espacial

#Creación de formula espacial
formula_espacial <- total_visitantes ~ tourism_gdp + crime_rate + college_education +
                    unemployment + business_activity + real_wage +pop_density

# Seleccionar datos de entrenamiento
set.seed(123)  # para reproducibilidad

# Variables predictoras
X_ml <- estados_datos_std %>%
  st_drop_geometry() %>%
  select(tourism_gdp, crime_rate, college_education, unemployment, 
         employment, business_activity, real_wage, pop_density)

# Variable dependiente
y_ml <- estados_datos_std$total_visitantes

# Entrenar el modelo Random Forest
modelo_rf <- randomForest(x = X_ml, y = y_ml, ntree = 500, importance = TRUE)

# Ver resultados
print(modelo_rf)
## 
## Call:
##  randomForest(x = X_ml, y = y_ml, ntree = 500, importance = TRUE) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 2
## 
##           Mean of squared residuals: 4.613301e+12
##                     % Var explained: 4.31
# Importancia de cada variable
importance(modelo_rf)
##                     %IncMSE IncNodePurity
## tourism_gdp       3.5447938  3.120073e+13
## crime_rate        0.0538395  8.626571e+12
## college_education 6.8567600  3.191269e+13
## unemployment      3.2192822  5.681166e+12
## employment        1.1349750  8.682096e+12
## business_activity 0.3369611  7.193952e+12
## real_wage         1.8161283  9.162462e+12
## pop_density       0.3626430  1.201754e+13
varImpPlot(modelo_rf)

Pregunta 4. Especificar y estimar 1 modelo de regresión de global espacial (e.g., SAR, SEM, SDM) que incluyan alguna de las 2 matrices de conectividad (W) seleccionadas en 2)

# SDM
model_SDM<- lagsarlm(formula_espacial, data = estados_datos_std, listw = pesos, zero.policy = TRUE, type="mixed") 
summary(model_SDM)
## 
## Call:lagsarlm(formula = formula_espacial, data = estados_datos_std, 
##     listw = pesos, type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1957395  -666180   -92008   731111  3644042 
## 
## Type: mixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##                       Estimate Std. Error z value  Pr(>|z|)
## (Intercept)            1495064     317330  4.7114  2.46e-06
## tourism_gdp            1374930     504781  2.7238 0.0064532
## crime_rate             -482576     274910 -1.7554 0.0791921
## college_education       409909     407063  1.0070 0.3139394
## unemployment           -542126     344241 -1.5748 0.1152917
## business_activity      -551471     309103 -1.7841 0.0744068
## real_wage               367577     386761  0.9504 0.3419099
## pop_density           -1168390     544294 -2.1466 0.0318239
## lag.tourism_gdp         162513    1323184  0.1228 0.9022501
## lag.crime_rate        -1981441     520527 -3.8066 0.0001409
## lag.college_education  2578928     971675  2.6541 0.0079519
## lag.unemployment      -3324327    1104626 -3.0095 0.0026171
## lag.business_activity  -732334     484433 -1.5117 0.1306018
## lag.real_wage          1907740     771552  2.4726 0.0134134
## lag.pop_density       -1430235    1401019 -1.0209 0.3073239
## 
## Rho: -0.66286, LR test value: 6.2771, p-value: 0.012231
## Approximate (numerical Hessian) standard error: 0.23625
##     z-value: -2.8058, p-value: 0.0050191
## Wald statistic: 7.8725, p-value: 0.0050191
## 
## Log likelihood: -493.5157 for mixed model
## ML residual variance (sigma squared): 1.3044e+12, (sigma: 1142100)
## Number of observations: 32 
## Number of parameters estimated: 17 
## AIC: 1021, (AIC for lm: 1025.3)

Pregunta 5. Identificar la posible presencia de autocorrelación espacial en los residuales Ɛi para cada uno de los modelos de regresión estimados en 3) y 4)

# Calcular predicciones
pred_RF <- predict(modelo_rf)

# Calcular residuos manualmente (observado - predicho)
residuos_RF <- y_ml - pred_RF

# Verificar tipo
is.numeric(residuos_RF)  # debe ser TRUE
## [1] TRUE
# Ahora sí puedes aplicar Moran's I
moran.test(residuos_RF, pesos, zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  residuos_RF  
## weights: pesos    
## 
## Moran I statistic standard deviate = 0.25289, p-value = 0.4002
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.020926088      -0.032258065       0.002007999
residuos_SDM <- residuals(model_SDM)

moran.test(residuos_SDM, pesos, zero.policy = TRUE)
## 
##  Moran I test under randomisation
## 
## data:  residuos_SDM  
## weights: pesos    
## 
## Moran I statistic standard deviate = -0.48264, p-value = 0.6853
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.08988856       -0.03225806        0.01425798

En este caso, vemos que al extraer el residuo de nuestros modelos, en niguno de ellos tenemos un valor significante (p<.1) Por lo que se considera que, no existe autocorrelación especial. Mi modelo SDM absorbió la dependencia espacial, lo que confirma que es un buen modelo.

Pregunta 6. Especificar y 1 modelo de regresión local espacial (e.g, GWR)

# Convertir a objeto Spatial si es sf
estados_sp <- as(estados_datos_std, "Spatial")

# Definir fórmula
formula_espacial1 <- total_visitantes ~ tourism_gdp + crime_rate + college_education +
                     unemployment + business_activity + real_wage + pop_density

# Calcular bandwidth óptimo y ajustar modelo GWR
modelo_gwr <- gwr.basic(
  formula = formula_espacial1,
  data = estados_sp,
  bw = bw.gwr(formula_espacial1, data = estados_sp, approach = "AICc", 
              kernel = "bisquare", adaptive = TRUE),
  kernel = "bisquare",
  adaptive = TRUE
)
## Adaptive bandwidth (number of nearest neighbours): 27 AICc value: 1027.464 
## Adaptive bandwidth (number of nearest neighbours): 25 AICc value: 1035.885 
## Adaptive bandwidth (number of nearest neighbours): 30 AICc value: 1025.743 
## Adaptive bandwidth (number of nearest neighbours): 30 AICc value: 1025.743
# Ver resumen 
summary(modelo_gwr)
##               Length Class                    Mode
## GW.arguments  11     -none-                   list
## GW.diagnostic  8     -none-                   list
## lm            14     lm                       list
## SDF           32     SpatialPolygonsDataFrame S4  
## timings        5     -none-                   list
## this.call      6     -none-                   call
## Ftests         0     -none-                   list

Pregunta 7. Apartir de la estimación de los modelos de regresión 4) - 6) seleccionar 1 - 2 modelos para la interpretación de los principales resultados

Para la actividad se seleccionaron los modelos Spatial Durbin Model y Geographically Weighted Regression:

Modelo SDM:

  • La variable tourism_gdp fue estadísticamente significativa (p = 0.006), con un efecto positivo sobre el total de visitantes, lo que indica que un mayor PIB turístico local se asocia a mayor número de visitantes.

  • El lag espacial de tourism_gdp también fue significativo (p = 0.002), lo cual sugiere que los estados vecinos con altos niveles de PIB turístico también influyen positivamente sobre el número de visitantes.

  • El valor del parámetro espacial ρ (rho) fue 0.66 con un p-valor = 0.012, lo cual indica una fuerte dependencia espacial, y justifica el uso de un modelo espacial.

Modelo GWR:

  • Al analizar los coeficientes estimados localmente para tourism_gdp, se observa que el efecto varía notablemente entre estados.
# modelo SDM
AIC(model_SDM)
## [1] 1021.031
# Para GWR
modelo_gwr$GW.diagnostic$AICc
## [1] 1025.743

Comparando los valores de AIC, el modelo SDM obtiene un AIC menor (1021.031) en comparación con el modelo GWR (AICc = 1025.743), lo cual sugiere que el SDM logra un mejor equilibrio dentro del modelo.

Pregunta 8. Mediante el uso de 2 - 3 mapas, visualizar los principales resultados estimados (e.g, prediccion de principal variable de análisis, significancia estadística de principal variable expectativa, precencia de autocorrelación espacial) de al menos 1 de los modelos de regresión seleccionados.

Mapa 1 Coeficiente local para tourism_gdp

# Agregar coeficientes de GWR al sf
estados_datos_std$coef_tourism_gdp <- modelo_gwr$SDF$tourism_gdp

# Visualizar mapa
library(tmap)
tm_shape(estados_datos_std) +
  tm_polygons("coef_tourism_gdp", palette = "RdBu", style = "quantile",
              title = "Coef. Local - tourism_gdp") +
  tm_layout(main.title = "Modelo GWR: Efecto local del PIB Turístico",
            legend.outside = TRUE)

El efecto del PIB turístico es más fuerte en estados del sureste.

Mapa 2 Significancia estadística local para tourism_gdp

# Calcular p-values locales
coef <- modelo_gwr$SDF$tourism_gdp
se <- modelo_gwr$SDF$"tourism_gdp_SE"
p_values <- 2 * (1 - pnorm(abs(coef / se)))
estados_datos_std$p_tourism_gdp <- p_values

# Visualizar mapa de p-values
tm_shape(estados_datos_std) +
  tm_polygons("p_tourism_gdp", palette = "-Blues", style = "quantile",
              title = "p-value de tourism_gdp") +
  tm_layout(main.title = "Significancia local del efecto tourism_gdp (GWR)",
            legend.outside = TRUE)

  • Las zonas con tonos azul oscuro indican mayor significancia estadística (p < 0.05).

  • Esto se concentra principalmente en el sur y centro del país, donde el modelo estima con mayor confianza que el PIB turístico influye en la llegada de visitantes.

Mapa 3 Predicción vs residuos

# Predicción de GWR
predicciones_gwr <- modelo_gwr$SDF$yhat
estados_datos_std$pred_gwr <- predicciones_gwr
estados_datos_std$resid_gwr <- estados_datos_std$total_visitantes - predicciones_gwr

# Mapa de residuos
tm_shape(estados_datos_std) +
  tm_polygons("resid_gwr", palette = "-RdBu", style = "quantile",
              title = "Residuales GWR") +
  tm_layout(main.title = "Errores de predicción - Modelo GWR",
            legend.outside = TRUE)
## Variable(s) "resid_gwr" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

### Mapa 4 t value de tourism_gdp

# Calcular t-value local
estados_datos_std$tval_tourism_gdp <- modelo_gwr$SDF$tourism_gdp / modelo_gwr$SDF$"tourism_gdp_SE"

# Graficar mapa del t-value
library(tmap)
tm_shape(estados_datos_std) +
  tm_polygons("tval_tourism_gdp", palette = "RdBu", style = "cont", midpoint = 0,
              title = "t-value - tourism_gdp") +
  tm_layout(main.title = "T-Value Local del PIB Turístico (GWR)",
            legend.outside = TRUE)

Conclusión

El modelo SDM mostró un buen ajuste general, reflejado en un menor valor de AIC y en la ausencia de autocorrelación espacial en los residuales. Además, permitió identificar efectos directos e indirectos significativos del PIB turístico, lo que confirma la influencia de las características de los estados vecinos en la dinámica del turismo.

Por otro lado, el modelo GWR ofreció una visión más detallada del territorio, revelando que el impacto del PIB turístico no es uniforme en todo el país. En estados del sur y sureste, se observó un efecto mayor y estadísticamente significativo. .

Referencias

  1. Esri. (2025). Conceptos básicos del análisis de regresión. ArcGIS Pro. Recuperado de https://pro.arcgis.com/es/pro-app/latest/tool-reference/spatial-statistics/regression-analysis-basics.html

  2. FasterCapital. (2025). Regresión ponderada geográficamente (GWR): local versus global – una inmersión profunda en la regresión ponderada geográficamente. Recuperado de https://fastercapital.com/es/contenido/Regresion-ponderada-geograficamente--GWR--local-versus-global--una-inmersion-profunda-en-la-regresion-ponderada-geograficamente.html