Introduccion

Este documento presenta una exploracion espacial sencilla del PIB turistico en Mexico.
La idea es visualizar patrones territoriales, medir si existe autocorrelacion espacial y comparar algunos modelos espaciales basicos.

Se utilizan dos niveles de analisis:

  • Corte transversal: para mapas y analisis espacial exploratorio.
  • Panel completo: para los modelos econometricos.

El objetivo es entender si existen patrones geograficos en la actividad turistica entre los estados del pais.

Carga de librerias

Se cargan los paquetes necesarios para trabajar con datos espaciales, mapas y modelos econometricos.

library(dplyr)
library(sf)
library(sp)
library(tmap)
library(spdep)
library(rgeoda)
library(ggplot2)
library(fields)
library(spatialreg)
library(stargazer)
library(plm)
library(splm)
library(pspatreg)

Configuracion general

Se define el modo de visualizacion de mapas y algunas opciones para evitar problemas graficos.

tmap_mode("plot")
graphics.off()
tmap_options(component.autoscale = FALSE)

Importacion de datos

Se carga la base de datos panel con variables economicas y el shapefile con los estados de Mexico.

mx_state <- read.csv("~/Proyecto R/CONCENTRACION/updated_tourism_data_model_2.csv")

mx_state_map <- st_read(
  "~/Proyecto R/CONCENTRACION/mx_maps/mx_maps/mx_states/mexlatlong.shp",
  quiet = TRUE
)

Union de datos espaciales

Se combinan los datos economicos con la geometria de los estados.

state_geodata <- mx_state_map %>%
  left_join(mx_state, by = c("OBJECTID" = "state_id"))

Construccion de la corte transversal

Para el analisis espacial exploratorio se necesita una sola observacion por estado.
Aqui se utiliza el año mas reciente disponible.

latest_year <- max(state_geodata$year, na.rm = TRUE)

state_geodata_cs <- state_geodata %>%
  filter(year == latest_year) %>%
  arrange(OBJECTID)

latest_year
[1] 2022
nrow(state_geodata_cs)
[1] 32

Estilo base para mapas

Se define un estilo visual simple para mantener consistencia en los mapas.

map_layout <- tm_layout(
  frame = FALSE,
  legend.position = c("left", "bottom"),
  legend.title.size = 0.8,
  legend.text.size = 0.7,
  inner.margins = c(0.04, 0.04, 0.08, 0.04)
)

Funcion para construir mapas

Esta funcion permite crear mapas coropleticos de forma mas eficiente.

make_map <- function(data, var, palette, legend_title, map_title, show_legend = TRUE) {
  tm_shape(data) +
    tm_polygons(
      col = var,
      palette = palette,
      style = "quantile",
      n = 8,
      title = legend_title,
      border.col = "gray70",
      lwd = 0.6,
      legend.show = show_legend
    ) +
    tm_title(map_title) +
    map_layout
}

Distribucion del PIB turistico

El siguiente mapa muestra como se distribuye el PIB turistico entre los estados.

tourism_map <- make_map(
  data = state_geodata_cs,
  var = "tourism_gdp",
  palette = "Blues",
  legend_title = "PIB turistico",
  map_title = "Importancia economica del turismo por estado"
)

tourism_map

Variables explicativas

A continuacion se muestran algunas variables que pueden ayudar a explicar diferencias territoriales.

Poder adquisitivo

real_wage_map <- make_map(
  state_geodata_cs,
  "real_wage",
  "BuGn",
  "Salario real",
  "Poder adquisitivo"
)

real_wage_map

Densidad poblacional

pop_density_map <- make_map(
  state_geodata_cs,
  "pop_density",
  "OrRd",
  "Habitantes por km cuadrado",
  "Densidad poblacional"
)

pop_density_map

Educacion superior

college_education_map <- make_map(
  state_geodata_cs,
  "college_education",
  "Purples",
  "% educacion superior",
  "Educacion Superior segun el Estado"
)

college_education_map

Inversion publica

public_investment_map <- make_map(
  state_geodata_cs,
  "ratio_public_investment",
  "YlGn",
  "Inversion publica",
  "Inversion publica relativa"
)

public_investment_map

Panel de mapas

Este panel resume visualmente varias variables relevantes.

tmap_arrange(
  real_wage_map,
  pop_density_map,
  college_education_map,
  public_investment_map,
  ncol = 2,
  nrow = 2
)

Matriz de contigüidad espacial

Se construye una matriz de vecindad tipo Queen, donde los estados vecinos comparten frontera o vertice.

mx_cs_map <- state_geodata_cs

swm <- poly2nb(mx_cs_map, queen = TRUE)
sswm <- nb2listw(swm, style = "W", zero.policy = TRUE)

summary(swm)
Neighbour list object:
Number of regions: 32 
Number of nonzero links: 138 
Percentage nonzero weights: 13.47656 
Average number of links: 4.3125 
Link number distribution:

1 2 3 4 5 6 7 8 9 
1 6 6 6 5 2 3 2 1 
1 least connected region:
31 with 1 link
1 most connected region:
8 with 9 links

Contigüidad tipo Queen

La matriz de contigüidad tipo Queen considera vecinos a los estados que comparten frontera o vertice.
Esto permite representar las relaciones espaciales directas entre entidades federativas.

mx_cs_map_sp <- as(mx_cs_map, "Spatial")
mx_cs_centroid <- coordinates(mx_cs_map_sp)

plot(mx_cs_map_sp, border = "gray60", axes = FALSE,
     main = "Conectividad espacial: matriz Queen")

plot(mx_cs_map_sp, col = "gray90", border = "white", add = TRUE)

plot(sswm,
     coords = mx_cs_centroid,
     pch = 19,
     cex = 0.2,
     col = "red",
     add = TRUE)

Vecinos mas cercanos (KNN)

En este enfoque cada estado se conecta con sus k vecinos mas cercanos en terminos de distancia geografica.

knn5 <- knn2nb(knearneigh(mx_cs_centroid, k = 5))
summary(knn5)
Neighbour list object:
Number of regions: 32 
Number of nonzero links: 160 
Percentage nonzero weights: 15.625 
Average number of links: 5 
Non-symmetric neighbours list
Link number distribution:

 5 
32 
32 least connected regions:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 with 5 links
32 most connected regions:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 with 5 links
plot(mx_cs_map_sp,
     border = "gray70",
     main = "Vecinos mas cercanos (k = 5)")

plot(knn5,
     mx_cs_centroid,
     pch = 19,
     cex = 0.4,
     add = TRUE,
     col = "red")

Clusteres espaciales locales (LISA)

El analisis LISA permite identificar agrupamientos locales de valores altos o bajos en el espacio.

sswm_geoda <- queen_weights(mx_cs_map)

lisa_tourism <- local_moran(
  sswm_geoda,
  state_geodata_cs["tourism_gdp"]
)

state_geodata_cs$cluster_tourism <- as.factor(
  lisa_tourism$GetClusterIndicators()
)

levels(state_geodata_cs$cluster_tourism) <- lisa_tourism$GetLabels()

ggplot(state_geodata_cs) +
  geom_sf(aes(fill = cluster_tourism), color = "white", linewidth = 0.2) +
  labs(
    title = "Clusters espaciales del PIB turistico",
    subtitle = paste("Year", latest_year),
    fill = "Tipo de cluster"
  ) +
  theme_minimal()

Modelos espaciales de corte transversal

Despues del analisis exploratorio, se estiman algunos modelos econometricos para evaluar la relacion entre el PIB turistico y las variables que mostraron mayor correlacion durante el EDA previo.

Las variables utilizadas en los modelos son:

  • salario real (real_wage)
  • densidad poblacional (pop_density)
  • educacion superior (college_education)
  • inversion publica relativa (ratio_public_investment)
  • PIB turistico rezagado (tourism_gdp_lag1)

Preparacion de los datos

Primero se prepara una base limpia para la estimacion de los modelos.

model_vars <- c(
  "tourism_gdp",
  "real_wage",
  "pop_density",
  "college_education",
  "ratio_public_investment",
  "tourism_gdp_lag1"
)

cs_model <- state_geodata_cs %>%
  st_drop_geometry() %>%
  select(all_of(model_vars)) %>%
  na.omit()

nrow(cs_model)
[1] 32

Modelo OLS

Se estima primero una regresion lineal tradicional como punto de referencia.

tourism_formula <- tourism_gdp ~ real_wage +
  pop_density +
  college_education +
  ratio_public_investment +
  tourism_gdp_lag1

model_ols <- lm(tourism_formula, data = cs_model)

summary(model_ols)

Call:
lm(formula = tourism_formula, data = cs_model)

Residuals:
   Min     1Q Median     3Q    Max 
-11366  -2975  -1397   3595  10565 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)             -5.241e+03  9.197e+03  -0.570 0.573681    
real_wage                1.538e+01  2.655e+01   0.579 0.567276    
pop_density              6.498e+00  1.598e+00   4.067 0.000393 ***
college_education       -3.703e+03  2.209e+04  -0.168 0.868198    
ratio_public_investment  2.203e+05  8.009e+04   2.751 0.010683 *  
tourism_gdp_lag1         1.085e+00  2.774e-02  39.129  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5250 on 26 degrees of freedom
Multiple R-squared:  0.9955,    Adjusted R-squared:  0.9946 
F-statistic:  1144 on 5 and 26 DF,  p-value: < 2.2e-16
AIC(model_ols)
[1] 646.3853

Modelo autorregresivo espacial (SAR)

Este modelo incorpora dependencia espacial directamente en la variable dependiente.

model_sar <- lagsarlm(
  tourism_formula,
  data = cs_model,
  listw = sswm,
  zero.policy = TRUE
)

summary(model_sar)

Call:lagsarlm(formula = tourism_formula, data = cs_model, listw = sswm, 
    zero.policy = TRUE)

Residuals:
     Min       1Q   Median       3Q      Max 
-11876.4  -3155.7  -1291.7   3104.6  11203.1 

Type: lag 
Coefficients: (numerical Hessian approximate standard errors) 
                           Estimate  Std. Error z value  Pr(>|z|)
(Intercept)             -3.2152e+03  8.5978e+03 -0.3740   0.70844
real_wage                1.4209e+01  2.3780e+01  0.5975   0.55015
pop_density              6.8804e+00  1.5175e+00  4.5341 5.786e-06
college_education       -3.5047e+03  2.0019e+04 -0.1751   0.86102
ratio_public_investment  2.1192e+05  7.2531e+04  2.9218   0.00348
tourism_gdp_lag1         1.0809e+00  2.5562e-02 42.2852 < 2.2e-16

Rho: -0.024572, LR test value: 0.54477, p-value: 0.46046
Approximate (numerical Hessian) standard error: 0.033146
    z-value: -0.74132, p-value: 0.4585
Wald statistic: 0.54956, p-value: 0.4585

Log likelihood: -315.9203 for lag model
ML residual variance (sigma squared): 22010000, (sigma: 4691.4)
Number of observations: 32 
Number of parameters estimated: 8 
AIC: 647.84, (AIC for lm: 646.39)

Modelo de error espacial (SEM)

En este caso la dependencia espacial se modela en el termino de error.

model_sem <- errorsarlm(
  tourism_formula,
  data = cs_model,
  listw = sswm,
  zero.policy = TRUE
)

summary(model_sem)

Call:errorsarlm(formula = tourism_formula, data = cs_model, listw = sswm, 
    zero.policy = TRUE)

Residuals:
     Min       1Q   Median       3Q      Max 
-11380.3  -3022.4  -1372.5   3475.8   9827.5 

Type: error 
Coefficients: (asymptotic standard errors) 
                           Estimate  Std. Error z value  Pr(>|z|)
(Intercept)             -5.4475e+03  7.8915e+03 -0.6903  0.490004
real_wage                1.1627e+01  2.3219e+01  0.5008  0.616539
pop_density              6.0754e+00  1.4454e+00  4.2033 2.631e-05
college_education        1.5415e+03  1.8842e+04  0.0818  0.934795
ratio_public_investment  2.2611e+05  7.1793e+04  3.1495  0.001636
tourism_gdp_lag1         1.0881e+00  2.5558e-02 42.5758 < 2.2e-16

Lambda: -0.16374, LR test value: 0.1588, p-value: 0.69026
Approximate (numerical Hessian) standard error: 0.43539
    z-value: -0.37607, p-value: 0.70686
Wald statistic: 0.14143, p-value: 0.70686

Log likelihood: -316.1132 for error model
ML residual variance (sigma squared): 22133000, (sigma: 4704.6)
Number of observations: 32 
Number of parameters estimated: 8 
AIC: 648.23, (AIC for lm: 646.39)

Rezago espacial

Se calcula el promedio ponderado de los valores de los estados vecinos.

state_geodata_cs$sp_lag_tourism <- lag.listw(
  sswm,
  state_geodata_cs$tourism_gdp,
  zero.policy = TRUE
)

Comparacion de modelos

Finalmente se comparan los resultados de los diferentes modelos.

stargazer(
  model_ols,
  model_sar,
  model_sem,
  type = "text",
  title = "Comparacion de modelos espaciales"
)

Comparacion de modelos espaciales
===============================================================================
                                          Dependent variable:                  
                        -------------------------------------------------------
                                              tourism_gdp                      
                                   OLS               spatial        spatial    
                                                  autoregressive     error     
                                   (1)                 (2)            (3)      
-------------------------------------------------------------------------------
real_wage                        15.384               14.209         11.627    
                                (26.550)             (23.780)       (23.219)   
                                                                               
pop_density                     6.498***             6.880***       6.075***   
                                 (1.598)             (1.518)        (1.445)    
                                                                               
college_education              -3,702.823           -3,504.738     1,541.527   
                              (22,093.800)         (20,018.630)   (18,841.860) 
                                                                               
ratio_public_investment       220,298.300**       211,920.400*** 226,110.000***
                              (80,088.990)         (72,530.850)   (71,793.280) 
                                                                               
tourism_gdp_lag1                1.086***             1.081***       1.088***   
                                 (0.028)             (0.026)        (0.026)    
                                                                               
Constant                       -5,240.656           -3,215.200     -5,447.496  
                               (9,196.790)         (8,597.814)    (7,891.456)  
                                                                               
-------------------------------------------------------------------------------
Observations                       32                   32             32      
R2                                0.995                                        
Adjusted R2                       0.995                                        
Log Likelihood                                       -315.920       -316.113   
sigma2                                            22,009,558.000 22,132,962.000
Akaike Inf. Crit.                                    647.841        648.226    
Residual Std. Error        5,249.569 (df = 26)                                 
F Statistic             1,144.489*** (df = 5; 26)                              
Wald Test (df = 1)                                    0.550          0.141     
LR Test (df = 1)                                      0.545          0.159     
===============================================================================
Note:                                               *p<0.1; **p<0.05; ***p<0.01

Autocorrelacion espacial global

El indice de Moran permite evaluar si valores similares tienden a agruparse en el espacio.

moran.test(state_geodata_cs$tourism_gdp, sswm)

    Moran I test under randomisation

data:  state_geodata_cs$tourism_gdp  
weights: sswm    

Moran I statistic standard deviate = 0.38237, p-value = 0.3511
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.006512234      -0.032258065       0.010281098

Cálculo del RMSE para comparar modelos

rmse <- function(model){
  sqrt(mean(residuals(model)^2))
}

rmse_results <- data.frame(
  Modelo = c("OLS","SAR","SEM"),
  RMSE = c(
    rmse(model_ols),
    rmse(model_sar),
    rmse(model_sem)
  )
)

rmse_results

Comparacion visual

tmap_arrange(
  make_map(state_geodata_cs,"tourism_gdp","Blues","PIB turistico","PIB turistico"),
  make_map(state_geodata_cs,"sp_lag_tourism","Blues","Rezago espacial","Rezago espacial del PIB turistico"),
  ncol = 2
)

Conclusion

El analisis espacial permite explorar si la actividad turistica presenta patrones segun su ubicacion.

Los mapas muestran diferencias claras en la distribucion del PIB turistico, con algunos estados concentrando una mayor actividad economica asociada al turismo.

Los resultados del analisis indican que tres variables presentan una relacion significativa con el PIB turistico: la densidad poblacional, la inversion publica relativa y el nivel de actividad turistica del periodo anterior. En particular, la variable rezagada del PIB turistico muestra un efecto muy fuerte, lo que sugiere que la actividad turistica tiende a persistir en el tiempo dentro de cada estado.

En contraste, variables como el salario real y el nivel de educacion superior no resultan estadisticamente significativas en los modelos estimados.

La comparacion entre el modelo OLS y los modelos espaciales (SAR y SEM) muestra que la inclusion de efectos espaciales no necesariamente mejora de forma sustancial el ajuste del modelo. Ademas, la prueba global de autocorrelacion espacial de Moran no resulta significativa, lo que sugiere que no existe una dependencia espacial fuerte en la distribucion del PIB turistico entre estados.

Dicho esto, el modelo SAR si mostró un un error medio más bajo (RMSE) que el resto de modelos estimados en este sprint. Comparados con el primer sprint, los modelos de corte transversal mostraron una mejor generalización, pues el RMSE de los modelos de la primera fase superaron en cifra los siete mil, con excepción del modelo de XGBOOST, que sigue teniendo el mejor desempeño general en cuanto a predicción.

En conjunto, los resultados sugieren que las diferencias en la actividad turistica entre estados estan mas asociadas a caracteristicas economicas y estructurales internas que a efectos de interaccion espacial entre entidades vecinas.