Actividad 1

Autores

Equipo 3 - Gamma.

Integrantes:
* César Romeo Vega
* Gamaliel Ostos
* Luis Carlos Borbón
* Rodrigo Arroyo

Introducción a la Actividad

1) Brevemente, describir con sus propias palabras qué es un ESDA y cuál es su principal propósito en el proceso de analítica de datos.

Un ESDA (Explanatory Spacial Data Analysis) por sus siglas en inglés es un conjunto de herramientas estadísticas que buscan comprender la estructura y distribución espacial de los datos antes de realizar un análisis más profundo o una modelación de datos,identificando patrones y relaciones espaciales en la información.

En comparación de un EDA (Explanatory Data Analysis) por sus siglas en inglés, el ESDA toma en cuenta el factor geoespacial para datos con ubicaciones geográfica y con posibilidad de tener autocorrelación espacial alta.

Algunos propósitos del ESDA son:

• Identifica patrones y anomalías espaciales en los datos de un periodo definido
• Visualiza la distribución espacial de los datos en mapas temáticos o de calor
• Detecta autocorrelación espacial entre locaciones vecinas
• Evalua supuestos para modelos espaciales mediante dependencias y variabilidad entre regiones

2) Brevemente, describir con sus propias palabras el concepto de autocorrelación espacial así como 1-2 ejemplos relacionados con dicho concepto.

La autocorrelación espacial es una medida que mide la relación entre el valor de una variable y las ubicaciones geográficas cercanas. Es decir, el valor de una región es influenciado por las características o valores de los lugares cercanos.

El Indice de Moran y el Indice de Geary son los indicadores más usados para medir la correlación espacial global y local de los datos.

Algunos ejemplos donde podemos observar autocorrelación espacial son:

• Propagación de enfermedades contagiosas en lugares cercanos (Pandemia de COVID-19)
• Distribución de poderes adquisitivos y nivel de pobreza en una ciudad (División entre San Pedro y el centro de Monterrey)
• Renta de vivienda temporal en distintas regiones urbanas (Zona TEC y Zonas Industriales de Nuevo León)

Situación Problema

El turismo ha sido uno de las mayores actividades económicas de México, donde destaca por sus amplias culturas y tradicionaes, gastronomía, celebraciones y recintos naturales/urbanos propicios para viajes vacacinonales o de negocios.

Con el Mundial de Futbol organizado por la FIFA en 2026, México será una de las cedes por lo que el turismo se verá beneficiado durante los 40 días que dura este macroevento. Sin embargo, ¿existirá alguna relación de beneficio entre estados que potencie el turismo? o ¿habrá algun patrón de los estados que generan más turismo con los que no?

Análisis Exploratorio Espacial de los Datos

Instalar Librerías, Paquetes y Bases de Datos

#install.packages("readxl")
#install.packages("ggplot2")
#install.packages("leaflet")
#install.packages("sf")
#install.packages("tmap")
#install.packages("tmaptools")
#install.packages("htmlwidgets")
#install.packages("naniar")
#install.packages("sp")
#install.packages("spdep")

library(readxl)
library(ggplot2)
library(leaflet)
library(sf)
library(tmap)
library(tmaptools)
library(dplyr)
library(htmlwidgets)
library(naniar)
library(sp)
library(spdep)

# Read the Excel file
data <- read_excel("tourism_state_data.xlsx")

Identificando Variables

Dentro de la Base de Datos se encuentra la siguientes variables por analizar: * tourism_gdp: Viajes y turismo que contribuye directamente al PIB
* crime_rate: Ratio de criminalidad por cada 100,000 personas
* unemployment: Porcentaje de población desempleada
* employement Porcentaje de población empleada
* business_activity: Índice económico ponderado por la distancia al puerto estadounidense más cercano
* real_wage: Salario real tomando de base el INPC del 2018 ($100 MXN)
* pop_density: Población por cada km^2 * good_governance: Ratio entre la inversión y deuda pública del estado
* ratio_public_investment: Ratio entre la inversión pública y el PIB del estado
* exchange_rate: Cambio de 1 USD por MXN * inpc: Índice Nacional de Precio al Consumidor tomando de base el 2018 ($100 MX)
* college_education: Porcentaje de población con por lo menos educación avanzada

str(data)

## tibble [544 × 18] (S3: tbl_df/tbl/data.frame)
##  $ state                  : chr [1:544] "Aguascalientes" "Baja California" "Baja California Sur" "Campeche" ...
##  $ year                   : num [1:544] 2006 2006 2006 2006 2006 ...
##  $ state_id               : num [1:544] 1057 2304 2327 1086 1182 ...
##  $ tourism_gdp            : num [1:544] 11995 54372 24238 13032 22730 ...
##  $ crime_rate             : num [1:544] 2.32 15.49 4.54 4.19 11.59 ...
##  $ college_education      : num [1:544] 0.1642 0.1808 0.1913 0.1544 0.0875 ...
##  $ unemployment           : num [1:544] 0.05 0.03 0.02 0.01 0.04 0.04 0.06 0.09 0.04 0.04 ...
##  $ employment             : num [1:544] 0.96 0.98 0.98 0.98 0.98 0.97 0.96 0.96 0.97 0.97 ...
##  $ business_activity      : num [1:544] -2.29 2.31 -2.43 -1.79 -2.49 -1.53 -2.19 -1.51 -2.57 -2.14 ...
##  $ real_wage              : num [1:544] 301 333 313 348 259 ...
##  $ pop_density            : num [1:544] 199.63 42.01 7.74 13.69 63.33 ...
##  $ good_governance        : num [1:544] 0.16 0.23 0.35 0.02 0.5 0.05 5.19 0 0.29 0.48 ...
##  $ ratio_public_investment: num [1:544] 0.01 0 0.01 0 0.01 0.02 0 0 0.01 0.01 ...
##  $ exchange_rate          : num [1:544] 10.8 10.8 10.8 10.8 10.8 ...
##  $ inpc                   : num [1:544] 62.7 62.7 62.7 62.7 62.7 ...
##  $ border_distance        : num [1:544] 625.59 8.83 800.32 978.33 1111.82 ...
##  $ region...17            : chr [1:544] "Bajio" "Norte" "Occidente" "Sur" ...
##  $ region...18            : num [1:544] 2 3 4 5 5 3 1 3 4 4 ...

# Formatear valores numéricos (2 decimales)
data$tourism_gdp <- round(data$tourism_gdp, 2)
data$crime_rate <- round(data$crime_rate, 2)
data$business_activity <- round(data$business_activity, 2)
data$good_governance <- round(data$good_governance, 2)
data$real_wage <- round(data$real_wage, 2)
data$pop_density <- round(data$pop_density, 2)

# Formatear valores porcentages
data$unemployment <- round(data$unemployment * 100, 2)
data$employment <- round(data$employment * 100, 2)  
data$college_education <- round(data$college_education * 100, 2)

Análisis Exploratorio

summary(data)

##     state                year         state_id     tourism_gdp    
##  Length:544         Min.   :2006   Min.   : 888   Min.   :  6240  
##  Class :character   1st Qu.:2010   1st Qu.:1047   1st Qu.: 22685  
##  Mode  :character   Median :2014   Median :1081   Median : 32482  
##                     Mean   :2014   Mean   :1219   Mean   : 56520  
##                     3rd Qu.:2018   3rd Qu.:1118   3rd Qu.: 59014  
##                     Max.   :2022   Max.   :2357   Max.   :472642  
##    crime_rate      college_education  unemployment      employment   
##  Min.   :  1.710   Min.   : 8.75     Min.   : 1.000   Min.   :89.00  
##  1st Qu.:  8.107   1st Qu.:16.70     1st Qu.: 3.000   1st Qu.:95.00  
##  Median : 13.880   Median :20.30     Median : 4.000   Median :97.00  
##  Mean   : 22.163   Mean   :21.11     Mean   : 4.251   Mean   :96.39  
##  3rd Qu.: 26.315   3rd Qu.:25.09     3rd Qu.: 5.000   3rd Qu.:97.54  
##  Max.   :181.510   Max.   :43.76     Max.   :10.000   Max.   :99.28  
##  business_activity   real_wage      pop_density      good_governance  
##  Min.   :-2.980    Min.   :239.3   Min.   :   7.74   Min.   :  0.000  
##  1st Qu.:-2.260    1st Qu.:282.5   1st Qu.:  39.56   1st Qu.:  0.180  
##  Median :-2.070    Median :306.2   Median :  61.77   Median :  0.500  
##  Mean   :-1.757    Mean   :314.9   Mean   : 299.46   Mean   :  2.362  
##  3rd Qu.:-1.768    3rd Qu.:335.4   3rd Qu.: 150.46   3rd Qu.:  1.350  
##  Max.   : 2.470    Max.   :481.7   Max.   :6211.45   Max.   :200.020  
##  ratio_public_investment exchange_rate        inpc        border_distance  
##  Min.   :0.000000        Min.   :10.85   Min.   : 62.69   Min.   :   8.83  
##  1st Qu.:0.000000        1st Qu.:12.87   1st Qu.: 74.93   1st Qu.: 613.26  
##  Median :0.000000        Median :14.51   Median : 87.19   Median : 751.64  
##  Mean   :0.005736        Mean   :15.91   Mean   : 89.08   Mean   : 704.92  
##  3rd Qu.:0.010000        3rd Qu.:19.47   3rd Qu.:103.02   3rd Qu.: 875.76  
##  Max.   :0.067644        Max.   :20.52   Max.   :126.48   Max.   :1252.66  
##  region...17         region...18   
##  Length:544         Min.   :1.000  
##  Class :character   1st Qu.:2.000  
##  Mode  :character   Median :3.000  
##                     Mean   :3.188  
##                     3rd Qu.:4.250  
##                     Max.   :5.000

Serie de Tiempo de Variable Dependiente

# Serie de tiempo de tourism_gdp
ggplot(data, aes(x = year, y = as.numeric(gsub(",", "", tourism_gdp)), group = state, color = state)) +  # Convert tourism_gdp to numeric
  geom_line() +
  labs(title = "GDP de Turismo por Estado a través del Tiempo",
       x = "Año",
       y = "GDP de Turismo") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

A lo largo de los años se aprecia que el PIB directo del Turismo no es estacionario durante el periodo de 2006 y 2022, teniendo una caida en 2020 debido a la pandemia. En caso de querer realizar un analisis temporal se necesitaría contar con registros mensuales o semanales que indiquen un patrón estacionario o en su defecto, hacer una transformación logarítmica.

Para este ejercicio, no se considerará el factor temporal y se enfocará en el análisis espacial.

Identificación de Valores Nulos

gg_miss_var(data, show_pct=TRUE)

Como se puede apreciar, no existen valores nulos en ninguna de las variables de la base de datos

Análisis descriptivo y de dispersión

# Histograma y Box Plot de tourism_gdp
options(repr.plot.width=12, repr.plot.height=8)

ggplot(data, aes(x = state, y = tourism_gdp)) +
  geom_col() +
  labs(x = "Estado", y = "Turismo PIB", title = "Histograma de Turismo PIB por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(data, aes(x = state, y = tourism_gdp)) +
  geom_boxplot() +
  labs(x = "Estado", y = "Turismo PIB", title = "Box Plot de Turismo PIB por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

La Ciudad de México supera en PIB de turismo a todos los demás estados. Los que menos tienen son Tlaxcala y Colima.

# Histograma de crime_rate
options(repr.plot.width=12, repr.plot.height=8)
ggplot(data, aes(x = state, y = crime_rate)) +
  geom_col() +
  labs(x = "Estado", y = "Crime Rate", title = "Histograma de Crime Rate por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(data, aes(x = state, y = crime_rate)) +
  geom_boxplot() +
  labs(x = "Estado", y = "Crime Rate", title = "Box Plot de Crime Rate por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Chihuahua es el que tiene mayor índice de crimen, seguido de Baja California y Sinaloa. Los que menos tienen son Yucatán y Aguascalientes.

# Histograma de college_education
options(repr.plot.width=12, repr.plot.height=8)

ggplot(data, aes(x = state, y = college_education)) +
  geom_col() +
  labs(x = "Estado", y = "College Education", title = "Histograma de College Education por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(data, aes(x = state, y = college_education)) +
  geom_boxplot() +
  labs(x = "Estado", y = "College Education", title = "Box Plot de College Education por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

La Ciudad de México es el estado con mayor educación. Chiapas, Michoacán, Oaxaca y Zacatecas son los que tienen menor educación.

# Histograma de unemployment
options(repr.plot.width=12, repr.plot.height=8)

ggplot(data, aes(x = state, y = unemployment)) +
  geom_col() +
  labs(x = "Estado", y = "Unemployment", title = "Histograma de Unemployment por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(data, aes(x = state, y = unemployment)) +
  geom_boxplot() +
  labs(x = "Estado", y = "Unemployment", title = "Box Plot de Unemployment por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

El estado con más desempleo es Coahuila, junto con Ciudad de México, Tabasco y Tlaxcala. Los que tienen menos son Morelos, Yucatán y Campeche.

# Histograma de real_wage
options(repr.plot.width=12, repr.plot.height=8)

ggplot(data, aes(x = state, y = real_wage)) +
  geom_col() +
  labs(x = "Estado", y = "Real Wage", title = "Histograma de Real wage por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(data, aes(x = state, y = real_wage)) +
  geom_boxplot() +
  labs(x = "Estado", y = "Real Wage", title = "Box Plot de Real wage por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Hay mucha variabilidad entre los estados, siendo Durango y Sinaloa los que menor valor de “Real Wage” poseen. Mientras que Ciudad de México y Campeche son las que más alto salario real tienen.

# Histograma de real_wage
options(repr.plot.width=12, repr.plot.height=8)

ggplot(data, aes(x = state, y = ratio_public_investment)) +
  geom_col() +
  labs(x = "Estado", y = "Public Investment (%)", title = "Histograma de Public Investment (%) por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(data, aes(x = state, y = ratio_public_investment)) +
  geom_boxplot() +
  labs(x = "Estado", y = "Public Investment (%)", title = "Box Plot de Public Investment (%) por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Si bien no existe mucha variabilidad en la mayoría de los estados, Baja California, Jalisco y Ciudad de México cuentan con un menor ratio entre la inversión publica estatal y el PIB nacional. Por otro lado, Oaxaca cuenta con el mayor ratio de inversión pública con respecto al PIB nacional.

# Histograma de business_activity
options(repr.plot.width=12, repr.plot.height=8)

ggplot(data, aes(x = state, y = business_activity)) +
  geom_col() +
  labs(x = "Estado", y = "Business Activity", title = "Histograma de Business Activity por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(data, aes(x = state, y = business_activity)) +
  geom_boxplot() +
  labs(x = "Estado", y = "Business Activity", title = "Box Plot de Business Activity por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

El único estado con actividad de negocios positiva es Baja California. Los más bajos son Quintana Roo y Yucatán.

# Histograma de business_activity
options(repr.plot.width=12, repr.plot.height=8)

ggplot(data, aes(x = state, y = good_governance)) +
  geom_col() +
  labs(x = "Estado", y = "Ratio Investment/Debt", title = "Histograma de Ratio Investment/Debt por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

ggplot(data, aes(x = state, y = good_governance)) +
  geom_boxplot() +
  labs(x = "Estado", y = "Ratio Investment/Debt", title = "Box Plot de Ratio Investment/Debt por Estado") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

Por otro ladoel estado con mayor ratio entre la inversión y deuda pública es Nuevo León, alcanzando el valor de 200.

Estas variables cuentan con algún estado como outlier, rompiendo el patrón de los valores de otros estados y marcando una dispersión más grande. Así mismo, las variables exchange_rate e inpc no cuentan con variabilidad geoespacial a lo largo de los años, por lo que no existen patrones visibles para incluirlos en el análisis espacial de los datos.

Por esta razón se eligirán las siguientes variables para medir el impacto espacial: * Tourism_gdp * Crime_rate * Unemployment * Real_wage * College_eductaion

Visualización de Variables en Mapas

mexico_states <- st_read('mexlatlong.shp')

## Reading layer `mexlatlong' from data source 
##   `C:\Users\rodio\OneDrive\Escritorio\8semestre\Bloque Final\Actividad 1\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

#Cambiar Distrito Federal a Ciudad de México
mexico_states$ADMIN_NAME <- gsub("Distrito Federal", "Ciudad de Mexico", mexico_states$ADMIN_NAME)

#Crear base de datos con datos geográficos
states <- merge(data, mexico_states, by.x = "state", by.y = "ADMIN_NAME")

# Crear el mapa base
GDP <- leaflet(states) %>% addTiles()

# Añadir capas para cada año
for (year in unique(states$year)) {
  # Subset the data and ensure it retains its spatial attributes
  states_subset <- states[states$year == year, ] 
  states_subset <- st_as_sf(states_subset)
  
  GDP <- GDP %>%
    addPolygons(
      data = states_subset,  # Use the subset here
      fillColor = ~colorNumeric("GnBu", domain = tourism_gdp)(tourism_gdp),
      fillOpacity = 0.7,
      weight = 2,
      color = "black",
      dashArray = "3",
      group = as.character(year), # Asignar un grupo para cada año
      popup = ~paste("Estado:", state, "<br>", "PIB Turismo:", tourism_gdp)
    )
}

# Añadir control de capas
GDP <- GDP %>%
  addLayersControl(
    baseGroups = unique(states$year), # Usar los años como grupos base
    options = layersControlOptions(collapsed = FALSE) # Mostrar el control expandido
  )

# Mostrar el mapa
GDP

Ciudad de México es quien posee más PIB del Turismo, sin embargo Quintana Roo, Estado de México, Jalisco, Guanajuato y Nuevo León han incrementado este valor en los últimos 15 años.

# Crear el mapa base
CRIM <- leaflet(states) %>% addTiles()

# Añadir capas para cada año
for (year in unique(states$year)) {
  # Subset the data and ensure it retains its spatial attributes
  states_subset <- states[states$year == year, ] 
  states_subset <- st_as_sf(states_subset)
  
  CRIM <- CRIM %>%
    addPolygons(
      data = states_subset,  # Use the subset here
      fillColor = ~colorNumeric("YlOrRd", domain = crime_rate)(crime_rate),
      fillOpacity = 0.7,
      weight = 2,
      color = "black",
      dashArray = "3",
      group = as.character(year), # Asignar un grupo para cada año
      popup = ~paste("Estado:", state, "<br>", "Tasa criminalidad (%):", crime_rate)
    )
}

# Añadir control de capas
CRIM <- CRIM %>%
  addLayersControl(
    baseGroups = unique(states$year), # Usar los años como grupos base
    options = layersControlOptions(collapsed = FALSE) # Mostrar el control expandido
  )

# Mostrar el mapa
CRIM

La variable de la tasa de criminalidad ha variado entre estados del 2006 al 2022 siendo Guerrero, Michoacan, Chihuahua, Colima, Baja California, Baja California Sur, Guanajuato y Zacatecas epicentros de una gran tasa de criminalidad por arriba del 70% en uno o varios años entre este periodo.

# Crear el mapa base
EMP <- leaflet(states) %>% addTiles()

# Añadir capas para cada año
for (year in unique(states$year)) {
  # Subset the data and ensure it retains its spatial attributes
  states_subset <- states[states$year == year, ] 
  states_subset <- st_as_sf(states_subset)
  
  EMP <- EMP %>%
    addPolygons(
      data = states_subset,  # Use the subset here
      fillColor = ~colorNumeric("YlOrRd", domain = unemployment)(unemployment),
      fillOpacity = 0.7,
      weight = 2,
      color = "black",
      dashArray = "3",
      group = as.character(year), # Asignar un grupo para cada año
      popup = ~paste("Estado:", state, "<br>", "Tasa de Desempleo:", unemployment)
    )
}

# Añadir control de capas
EMP <- EMP %>%
  addLayersControl(
    baseGroups = unique(states$year), # Usar los años como grupos base
    options = layersControlOptions(collapsed = FALSE) # Mostrar el control expandido
  )

# Mostrar el mapa
EMP

La tasa de desempleo ha afectado a la mayoría de los estados de México siendo los años 2012, 2015, 2017 y 2021 donde más estados han tenido tasas mayores a 5 puntos. Por otro lado en años como 2022, 2019, 2016 y 2008, solo algunos estados de la región del norte y sur cuentan con una alta tasa de desempleo, sin embargo la mayoria son menores a los 5 puntos.

# Crear el mapa base
WAGE <- leaflet(states) %>% addTiles()

# Añadir capas para cada año
for (year in unique(states$year)) {
  # Subset the data and ensure it retains its spatial attributes
  states_subset <- states[states$year == year, ] 
  states_subset <- st_as_sf(states_subset)
  
  WAGE <- WAGE %>%
    addPolygons(
      data = states_subset,  # Use the subset here
      fillColor = ~colorNumeric("GnBu", domain = real_wage)(real_wage),
      fillOpacity = 0.7,
      weight = 2,
      color = "black",
      dashArray = "3",
      group = as.character(year), # Asignar un grupo para cada año
      popup = ~paste("Estado:", state, "<br>", "Salario Real:", real_wage)
    )
}

# Añadir control de capas
WAGE <- WAGE %>%
  addLayersControl(
    baseGroups = unique(states$year), # Usar los años como grupos base
    options = layersControlOptions(collapsed = FALSE) # Mostrar el control expandido
  )

# Mostrar el mapa
WAGE

Los estados con más salario real de los últimos 5 años se encuentran en la zona noroeste, suroeste y el centro del país, resaltando Campeche, Ciudad de México, Nuevo León y Querétaro. Mientras qye la zona centro-este y sur tienen las concentraciones más bajas entorno al salario real

# Crear el mapa base
EDU <- leaflet(states) %>% addTiles()

# Añadir capas para cada año
for (year in unique(states$year)) {
  # Subset the data and ensure it retains its spatial attributes
  states_subset <- states[states$year == year, ] 
  states_subset <- st_as_sf(states_subset)
  
  EDU <- EDU %>%
    addPolygons(
      data = states_subset,  # Use the subset here
      fillColor = ~colorNumeric("GnBu", domain = college_education)(college_education),
      fillOpacity = 0.7,
      weight = 2,
      color = "black",
      dashArray = "3",
      group = as.character(year), # Asignar un grupo para cada año
      popup = ~paste("Estado:", state, "<br>", "Educación Avanzada (%):", college_education,
                     "<br>", "Densidad Poblacional:", pop_density)
    )
}

# Añadir control de capas
EDU <- EDU %>%
  addLayersControl(
    baseGroups = unique(states$year), # Usar los años como grupos base
    options = layersControlOptions(collapsed = FALSE) # Mostrar el control expandido
  )

# Mostrar el mapa
EDU

En términos generales, la mayoría de los estados ha tenido un incremento en la tasa de educación avanzada de sus habitantes siendo muy alto en la zona noreste, noroeste, el centro de méxico y el suroeste del país. Mientras que en la zona sur y el centro-norte cuentan con las menores tasas de educación avanzada.

Matrices de Conectividad Rook y Queen

# Load the datasets
mexlatlong <- st_read("mexlatlong.shp")  # Read shapefile

## Reading layer `mexlatlong' from data source 
##   `C:\Users\rodio\OneDrive\Escritorio\8semestre\Bloque Final\Actividad 1\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

data <- read_excel("tourism_state_data.xlsx")  # Read Excel file

# Merge the spatial data with the attribute data based on the state column
states <- mexlatlong %>%
  left_join(data, by = c("ADMIN_NAME" = "state"))  # Ensure column names match

# Define the variables to analyze
#variables <- c("tourism_gdp", "crime_rate", "college_education",
#               "unemployment", "real_wage")

variables <- c("tourism_gdp", "crime_rate", "college_education",
               "unemployment", "employment", "business_activity",
               "real_wage", "pop_density", "good_governance",
               "ratio_public_investment")

# Convertir 'states' a un objeto espacial
data_sp <- na.omit(states[, c(variables, "geometry")])

mexlatlong_sp <- as(data_sp, "Spatial")

# Crear la matriz de vecindad Queen
nb_queen <- poly2nb(mexlatlong_sp, queen = TRUE)

# Crear la matriz de pesos espaciales
lw_queen <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)

# Obtener las coordenadas de los centroides
centroids <- coordinates(mexlatlong_sp)

plot(mexlatlong_sp,border="blue",axes=FALSE,las=1, main="Mexico Spatial Connectivity Matrix - Queen")
plot(mexlatlong_sp,col="grey",border=grey(0.9),axes=T,add=T) 
plot(lw_queen, centroids, add=TRUE, col='dark green')

# Crear la matriz de vecindad Rook
nb_rook <- poly2nb(mexlatlong_sp, queen = FALSE)

# Crear la matriz de pesos espaciales
lw_rook_subset <- nb2listw(nb_rook, style = "W", zero.policy = TRUE)

plot(mexlatlong_sp,border="red",axes=FALSE,las=1, main="Mexico Spatial Connectivity Matrix - Rook")
plot(mexlatlong_sp,col="grey",border=grey(0.9),axes=T,add=T) 
plot(lw_queen, centroids, add=TRUE, col='red')

Global Moran’s I Statistic - Autocorrelación Global y Local

# Convert variables to numeric to avoid errors
for (var in variables) {
  states[[var]] <- as.numeric(states[[var]])
}

# Loop through each variable for Global & Local Moran's I
for (var in variables) {
  cat("\n--- Global Moran's I for", var, "---\n")

  # Check for missing values and remove them for this variable
  data_for_moran <- na.omit(states[, c(var, "geometry")])

  if (nrow(data_for_moran) > 1) { # Ensure enough data after removing NAs

    # Convert the subset to a spatial object
    mexlatlong_sp <- as(data_for_moran, "Spatial")

    # Recalculate neighbors and spatial weights for the subset data
    # Rook contiguity
    nb_rook <- poly2nb(mexlatlong_sp, queen = FALSE)
    lw_rook_subset <- nb2listw(nb_rook, style = "W", zero.policy = TRUE)

    # Queen contiguity
    nb_queen <- poly2nb(mexlatlong_sp, queen = TRUE)
    lw_queen_subset <- nb2listw(nb_queen, style = "W", zero.policy = TRUE)

    # Global Moran’s I - Use the Rook contiguity weights
    global_moran_rook <- moran.test(data_for_moran[[var]],
                                    listw = lw_rook_subset, zero.policy = TRUE,
                                    na.action = na.omit)
    cat("\n--- Global Moran's I with Rook contiguity ---\n")
    print(global_moran_rook)

    # Global Moran’s I - Use the Queen contiguity weights
    global_moran_queen <- moran.test(data_for_moran[[var]],
                                     listw = lw_queen_subset, zero.policy = TRUE,
                                     na.action = na.omit)
    cat("\n--- Global Moran's I with Queen contiguity ---\n")
    print(global_moran_queen)

    # Moran Scatterplot - Rook contiguity
    moran.plot(data_for_moran[[var]], lw_rook_subset, main = paste("Moran’s I (Rook) -", var))

    # Moran Scatterplot - Queen contiguity
    moran.plot(data_for_moran[[var]], lw_queen_subset, main = paste("Moran’s I (Queen) -", var))

    # Local Moran’s I - Use the Rook contiguity weights
    local_moran_res_rook <- localmoran(data_for_moran[[var]], lw_rook_subset, zero.policy = TRUE)

    # Local Moran’s I - Use the Queen contiguity weights
    local_moran_res_queen <- localmoran(data_for_moran[[var]], lw_queen_subset, zero.policy = TRUE)

    # Append Local Moran’s I values and p-values to the dataset for both Rook and Queen
    pval_name_rook <- paste0("pval_", var, "_rook")
    cluster_name_rook <- paste0("cluster_", var, "_rook")
    pval_name_queen <- paste0("pval_", var, "_queen")
    cluster_name_queen <- paste0("cluster_", var, "_queen")

    states[[pval_name_rook]] <- NA  # Initialize as NA
    states[[cluster_name_rook]] <- "Not Significant"  # Default value
    states[[pval_name_queen]] <- NA  # Initialize as NA
    states[[cluster_name_queen]] <- "Not Significant"  # Default value

    # Transfer Local Moran's results for Rook and Queen into the merged_data
    states[rownames(local_moran_res_rook), pval_name_rook] <- local_moran_res_rook[, 5]
    states[rownames(local_moran_res_rook), cluster_name_rook] <- ifelse(local_moran_res_rook[, 5] < 0.05, "Significant", "Not Significant")

    states[rownames(local_moran_res_queen), pval_name_queen] <- local_moran_res_queen[, 5]
    states[rownames(local_moran_res_queen), cluster_name_queen] <- ifelse(local_moran_res_queen[, 5] < 0.05, "Significant", "Not Significant")

    # Convert to factor
    states[[cluster_name_rook]] <- as.factor(states[[cluster_name_rook]])
    states[[cluster_name_queen]] <- as.factor(states[[cluster_name_queen]])

    # Local cluster map using ggplot for Rook contiguity
    print(
      ggplot(states) +
        geom_sf(aes(fill = !!sym(cluster_name_rook))) +
        scale_fill_manual(values = c("Significant" = "red", "Not Significant" = "grey")) +
        labs(title = paste("Local Spatial Autocorrelation -", var, "(Rook)"), fill = "Cluster")
    )

    # Local cluster map using ggplot for Queen contiguity
    print(
      ggplot(states) +
        geom_sf(aes(fill = !!sym(cluster_name_queen))) +
        scale_fill_manual(values = c("Significant" = "red", "Not Significant" = "grey")) +
        labs(title = paste("Local Spatial Autocorrelation -", var, "(Queen)"), fill = "Cluster")
    )
  } else {
    cat("Skipping", var, "due to missing or invalid data.\n")
  }
}

## 
## --- Global Moran's I for tourism_gdp ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 10.52, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      6.589445e-02     -1.901141e-03      4.152932e-05 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 2.8037, p-value = 0.002526
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.580368e-02     -1.901141e-03      3.987752e-05

## 
## --- Global Moran's I for crime_rate ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 27.072, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.717773e-01     -1.901141e-03      4.115763e-05 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 28.098, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.747368e-01     -1.901141e-03      3.952118e-05

## 
## --- Global Moran's I for college_education ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 25.741, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.644027e-01     -1.901141e-03      4.173859e-05 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 27.213, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.703747e-01     -1.901141e-03      4.007815e-05

## 
## --- Global Moran's I for unemployment ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 11.857, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      7.460807e-02     -1.901141e-03      4.163943e-05 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 11.411, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      7.025300e-02     -1.901141e-03      3.998308e-05

## 
## --- Global Moran's I for employment ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 14.873, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0940311498     -0.0019011407      0.0000416063 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 14.299, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      8.847970e-02     -1.901141e-03      3.995132e-05

## 
## --- Global Moran's I for business_activity ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 49.039, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      3.124443e-01     -1.901141e-03      4.108953e-05 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 49.286, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.3076820338     -0.0019011407      0.0000394559

## 
## --- Global Moran's I for real_wage ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 23.309, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.485904e-01     -1.901141e-03      4.168397e-05 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 22.319, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.393045e-01     -1.901141e-03      4.002578e-05

## 
## --- Global Moran's I for pop_density ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 71.396, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      4.553038e-01     -1.901141e-03      4.100861e-05 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 72.35, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      4.521094e-01     -1.901141e-03      3.937831e-05

## 
## --- Global Moran's I for good_governance ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 4.0102, p-value = 3.033e-05
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.608462e-02     -1.901141e-03      2.011524e-05 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 3.6914, p-value = 0.0001115
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      1.433605e-02     -1.901141e-03      1.934785e-05

## 
## --- Global Moran's I for ratio_public_investment ---
## 
## --- Global Moran's I with Rook contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_rook_subset    
## 
## Moran I statistic standard deviate = 9.4053, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      0.0582749306     -0.0019011407      0.0000409357 
## 
## 
## --- Global Moran's I with Queen contiguity ---
## 
##  Moran I test under randomisation
## 
## data:  data_for_moran[[var]]  
## weights: lw_queen_subset    
## 
## Moran I statistic standard deviate = 9.3924, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      5.698597e-02     -1.901141e-03      3.930841e-05

Tras analizar las variables de la base de datos, se determinó que todas las variables son significativamente estadísticas, sin embargo la mayoría cuenta con un Índice de Global Moran muy cercano a 0, indicando mínima autocorrelación espacial de los datos siendo la variable de employment, ratio_public_investment, business_activity la de mayor indice en autocorrelación espacial.

En cuestión de la autocorrelación local representada por los mapas de significancia, vemos como variables como tourism_gdp y employment consideran tanto la región norte, centro y sur significativa dejando los estados más al oeste y este fuera de estos clusters.

Por otro lado, pop_density, business_activity, y crime_rate consideran a casi el 80% de los estados de la republica mexicana como significativos. Cabe aclarar que este patrón tambien se puede revisar contrariamente donde good_governance y ratio_public_investmente consideran menos del 30% de los estados como significativos; estando muy segregados para formar varios clusters pequeños.

Conclusiones y Hallazgos

1. Presencia de autocorrelación espacial: Encontramos que hay autocorrelación espacial significativa en todas las variables analizadas (p-valor < 0.05), lo que indica que los valores en un estado si llegan a estar influenciados por los estados vecinos tanto en la matriz de conectividad rook como queen.

2. Valor del Índice de Global Moran’s I en general: Aunque los índices fueron significativos, habían variables muy bajas como turismo_gdp cuya correlación espacial era débil mientras que variables como employment y ratio_public_investment en matriz rook y business_activity en matriz queen mostraron una autocorrelación espacial más fuerte, aunque sin superar el valor de 0.5 de Global Moran’s

3. Turismo (tourism_gdp): La Ciudad de México destaca por tener consistentemente el mayor PIB turístico, mientras que estados como Jalisco, Nuevo León, Guanajuato y Quintana Roo han mostrado un crecimiento importante a lo largo de los años. Esto se demuestra en la significacia de los clusters, pues el centro y norte del pais fueron identificados como regiones significativas para la el análisis de autocorrelación espacial.

4. Educación avanzada (college_education): En cuanto a la educación, la proporción de educación superior están en el norte y centro del país (CDMX, NL, QRO). En contraste, estados del sur (Chiapas, Oaxaca, Guerrero) formanron los clusters de bajo nivel educativo por lo que valdría la pena realizar un análisis a profundidad de las posibles causas de este nivel. En cuestión de autocorrelación local, tanto la parte noroeste, centro y suroeste fueron significativas.

5. Salario real (real_wage): El centro y noroeste del país concentran los salarios más altos, destacando CDMX y Campeche. Hay una autocorrelación moderada entre estados vecinos con niveles similares de salario real posiblemente por términos de inversión extranjera o cantidad de empleos.

6. Gobernanza (good_governance): Si bien Nuevo León destaca con un valor extremadamente alto, esto no es mas que un outlier. A pesar de la significancia estadística, la autocorrelación es baja, sugiriendo que la gobernanza es más un fenómeno local y no es influenciado por los demas estados.

7. Densidad poblacional (pop_density): Esta variable mostró la más alta autocorrelación espacial. Los estados más densamente poblados (CDMX, EDOMEX) están agrupados, y lo mismo ocurre con los menos densos (BCS, Durango), podría ser interesante si está relacionada con el tamaño del terreno .

8. Inversión pública (ratio_public_investment): Sorpresivamente, Oaxaca destaca como el estado con mayor inversión pública respecto al PIB, formando un cluster significativo con vecinos del sur, aunque la autocorrelación no fue de las más altas.

9. Clusters positivos para el turismo: Los estados que han incrementado su PIB turístico están cerca de otros con desempeño similar (CDMX-EDOMEX-Guanajuato). Esto sugiere que políticas de turismo pueden tener efectos de derrame positivo regional y que probablemente la gente que visita ciertos estados, esta interesada en conocer los estados colindantes.

10. Zonas prioritarias para intervención: Estados con clusters de alto crimen, bajo empleo y baja educación (Chiapas, Oaxaca, Guerrero, Zacatecas) deberían ser priorizados para intervenciones integrales en seguridad, empleo y educación y poder mejorar estas estadísticas en el futuro.