Introducción

AE-1.2. Modelos de Regresión Espacial (actividad en equipo). BD utilizada: regional_price_dataset.csv

  • Ángel Santiago Díaz Castillo
  • Angie Michelle Zerón Monge
  • Eli Gabriel Hernánez Medina
  • Patricio Javier Pérez Fajardo

Importación de los datos

Revisamos la base de datos regional_price_dataset.csv, que contiene 1,491 registros diarios entre enero de 2018 y mayo de 2019. Las variables incluidas son: fecha, zona, precio y cantidad. La información está distribuida en tres regiones: Centro, Norte y Sur.

Los precios registran un promedio de $11.55, con valores que van desde $5.10 hasta $18.50. En cuanto a las cantidades vendidas, el promedio diario es de 195 unidades, con una dispersión que alcanza hasta 405 unidades. Las observaciones están equilibradas entre las tres zonas, lo que permite comparaciones directas.

En general, observamos cierta estabilidad en los precios y mayor variación en las cantidades. Para el análisis espacial, será necesario asignar coordenadas representativas a cada zona. Más adelante, sería ideal trabajar con datos más desagregados, como a nivel municipal, para obtener mayor precisión.

Transformación de los datos

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(lubridate)
## 
## Attaching package: 'lubridate'
## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union
# Agregamos por mes y zona
data$fecha <- as.Date(data$fecha)
data <- data %>%
  mutate(mes = floor_date(fecha, "month")) %>%
  group_by(zona, mes) %>%
  summarise(precio_promedio = mean(precio),
            cantidad_total = sum(cantidad), .groups = "drop")

Visualización de los datos y EDA

library(leaflet)

zonas_coord <- data.frame(
  zona = c("Centro", "Norte", "Sur"),
  lat = c(19.4326, 25.6866, 16.7539),
  lng = c(-99.1332, -100.3161, -93.1131)
)

latest_month <- max(data$mes)
latest_data <- data %>% filter(mes == latest_month)

map_data <- merge(latest_data, zonas_coord, by = "zona")

leaflet(map_data) %>%
  addTiles() %>%
  addCircleMarkers(~lng, ~lat,
                   radius = ~sqrt(cantidad_total) / 2,
                   color = ~ifelse(precio_promedio > 14, "red", "blue"),
                   label = ~paste(zona, "<br>Precio:", round(precio_promedio,2), "<br>Cantidad:", cantidad_total),
                   fillOpacity = 0.7)

En el mapa interactivo se muestran las tres zonas geográficas consideradas en el análisis: Norte, Centro y Sur, representadas por círculos proporcionales a la cantidad total vendida en el mes más reciente disponible. Como se observa, las regiones Centro y Sur presentan volúmenes de venta comparables, mientras que la zona Norte muestra un tamaño ligeramente menor. Esto sugiere que, en ese periodo, la demanda se concentró principalmente en el centro y sur del país. Además, los colores indican que el precio promedio se mantiene relativamente estable entre zonas, con ligeras variaciones. Esta visualización permite identificar rápidamente patrones regionales de consumo y precios.

Autocorrelación espacial y clústers

# Cargar librerías necesarias
library(dplyr)
library(sf)
## Linking to GEOS 3.11.0, GDAL 3.5.3, PROJ 9.1.0; sf_use_s2() is TRUE
library(spdep)
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
# Asignamos coordenadas representativas de cada zona
zonas_coord <- data.frame(
  zona = c("Centro", "Norte", "Sur"),
  lat = c(19.4326, 25.6866, 16.7539),
  lng = c(-99.1332, -100.3161, -93.1131)
)

library(readr)

data <- read_csv("/Users/gabrielmedina/Downloads/regional_price_dataset.csv")
## Rows: 1491 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (1): zona
## dbl  (2): precio, cantidad
## date (1): fecha
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Calculamos el precio promedio general por zona
zona_promedios <- data %>%
  group_by(zona) %>%
  summarise(precio_promedio = mean(precio),
            cantidad_total = sum(cantidad),
            .groups = "drop")

# Unimos coordenadas a los datos agregados
map_data <- merge(zona_promedios, zonas_coord, by = "zona")

# Convertimos a objeto espacial
zona_sf <- st_as_sf(map_data, coords = c("lng", "lat"), crs = 4326)

# Creamos vecinos por distancia entre puntos (hasta 1000 km)
coords <- st_coordinates(zona_sf)
nb <- dnearneigh(coords, 0, 1000)
lw <- nb2listw(nb, style = "W", zero.policy = TRUE)

# Variabl
# Usamos la cantidad total como variable dependiente para el análisis espacial
spatial_values <- map_data$cantidad_total

Debido a la limitada cantidad de zonas geográficas disponibles en la base de datos (solo tres: Centro, Norte y Sur), no fue posible aplicar de manera válida los indicadores de autocorrelación espacial global (Moran’s I) ni local (LISA). Ambos requieren una cantidad mínima de unidades espaciales para garantizar la robustez estadística de los resultados. Sin embargo, esta limitación no impide el desarrollo de modelos de regresión espacial ajustados a los datos disponibles.

Modelos de regresion especial

SAR

# Creamos dummies manualmente
map_data$zona_norte <- ifelse(map_data$zona == "Norte", 1, 0)
map_data$zona_sur <- ifelse(map_data$zona == "Sur", 1, 0)

# Modelo SAR con variables dummy
library(spatialreg)
## Loading required package: Matrix
## 
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
## 
##     get.ClusterOption, get.coresOption, get.mcOption,
##     get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption
modelo_sar <- lagsarlm(cantidad_total ~ zona_norte + zona_sur, 
                       data = as.data.frame(map_data), 
                       listw = lw, 
                       zero.policy = TRUE)
## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value
## Warning in lagsarlm(cantidad_total ~ zona_norte + zona_sur, data = as.data.frame(map_data), : inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16 
##   system is exactly singular - using numerical Hessian.
summary(modelo_sar)
## 
## Call:
## lagsarlm(formula = cantidad_total ~ zona_norte + zona_sur, data = as.data.frame(map_data), 
##     listw = lw, zero.policy = TRUE)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##      0      0      0      0      0 
## 
## Type: lag 
## Coefficients: (numerical Hessian approximate standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)    16408        NaN     NaN      NaN
## zona_norte    -11814        NaN     NaN      NaN
## zona_sur      -37412        NaN     NaN      NaN
## 
## Rho: 1, LR test value: NaN, p-value: NA
## Approximate (numerical Hessian) standard error: NaN
##     z-value: NaN, p-value: NA
## Wald statistic: NaN, p-value: NA
## 
## Log likelihood: NaN for lag model
## ML residual variance (sigma squared): 0, (sigma: 0)
## Number of observations: 3 
## Number of parameters estimated: 5 
## AIC: NaN, (AIC for lm: -Inf)

Se intentó estimar un modelo SAR (Spatial Autoregressive Regression) utilizando variables dummy por zona como explicativas de la cantidad total vendida. Sin embargo, debido a la muy baja cantidad de observaciones (solo tres zonas), el modelo presentó problemas numéricos graves: no fue posible calcular errores estándar, los coeficientes estimados no son significativos y el sistema resultó ser singular, es decir, matemáticamente indeterminado. Esto provocó que se reemplazaran valores no definidos (NA/Inf) por valores máximos, sin lograr una solución estable. En resumen, el modelo no es confiable ni interpretable en estas condiciones, lo cual es una limitación natural al trabajar con un conjunto de datos tan reducido y homogéneo.

SEM

library(spatialreg)

# de que sea un data.frame plano (no sf)
map_df <- as.data.frame(map_data)

# Modelo SEM
modelo_sem <- errorsarlm(cantidad_total ~ zona_norte + zona_sur, 
                         data = map_df, 
                         listw = lw, 
                         zero.policy = TRUE)
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
## Warning in log(s2): NaNs produced
## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value

## Warning in optimize(sar.error.f, interval = interval, maximum = TRUE, tol =
## con$tol.opt, : NA/Inf replaced by maximum positive value
summary(modelo_sem)
## 
## Call:errorsarlm(formula = cantidad_total ~ zona_norte + zona_sur, 
##     data = map_df, listw = lw, zero.policy = TRUE)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##      0      0      0      0      0 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   107936        NaN     NaN      NaN
## zona_norte     -7876        NaN     NaN      NaN
## zona_sur      -24941        NaN     NaN      NaN
## 
## Lambda: -0.60318, LR test value: -Inf, p-value: < 2.22e-16
## Asymptotic standard error: NaN
##     z-value: NaN, p-value: NA
## Wald statistic: NaN, p-value: NA
## 
## Log likelihood: 13.78073 for error model
## ML residual variance (sigma squared): 0, (sigma: 0)
## Number of observations: 3 
## Number of parameters estimated: 5 
## AIC: -17.561, (AIC for lm: -Inf)

El modelo de error espacial (SEM) fue estimado utilizando variables dummy por zona como regresores de la cantidad total vendida. Aunque el modelo logró ajustarse matemáticamente, los resultados reflejan las limitaciones de contar únicamente con tres observaciones. Los coeficientes estimados indican que, en comparación con la zona Centro (categoría base), la zona Norte registra en promedio 7,876 unidades menos y la zona Sur 24,941 unidades menos. Sin embargo, los errores estándar, los valores z y los p-values aparecen como NaN, lo que impide confirmar la significancia estadística de estos resultados.

La estimación del parámetro Lambda, que mide la correlación espacial en los errores, arrojó un valor de -0.60, lo cual sugiere una posible correlación negativa entre los errores de zonas vecinas. A pesar de esto, el estadístico de razón de verosimilitud (LR test) no fue interpretable debido a problemas numéricos (valor -Inf), y la varianza residual estimada fue cero, lo que implica que el modelo ajusta exactamente los valores observados, aunque probablemente por sobreajuste.

En resumen, aunque el modelo SEM pudo calcularse, los resultados no son estadísticamente válidos debido al número extremadamente reducido de observaciones. Esta limitación resalta la necesidad de contar con una mayor desagregación geográfica para que los modelos espaciales generen resultados sólidos.

SDM

# Modelo SDM (Durbin)
modelo_sdm <- lagsarlm(cantidad_total ~ zona_norte + zona_sur,
                       data = map_df,
                       listw = lw,
                       type = "mixed",
                       zero.policy = TRUE)
## Warning in lagsarlm(cantidad_total ~ zona_norte + zona_sur, data = map_df, :
## Aliased variables found: lag.zona_norte lag.zona_sur
## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value

## Warning in optimize(sar.lag.mixed.f, interval = interval, maximum = TRUE, :
## NA/Inf replaced by maximum positive value
## Warning in lagsarlm(cantidad_total ~ zona_norte + zona_sur, data = map_df, : inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16 
##   system is exactly singular - using numerical Hessian.
summary(modelo_sdm)
## 
## Call:lagsarlm(formula = cantidad_total ~ zona_norte + zona_sur, data = map_df, 
##     listw = lw, type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##      0      0      0      0      0 
## 
## Type: mixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##     (2 not defined because of singularities)
##                Estimate Std. Error z value Pr(>|z|)
## (Intercept)       16408        NaN     NaN      NaN
## zona_norte       -11814        NaN     NaN      NaN
## zona_sur         -37412        NaN     NaN      NaN
## lag.zona_norte       NA         NA      NA       NA
## lag.zona_sur         NA         NA      NA       NA
## 
## Rho: 1, LR test value: NaN, p-value: NA
## Approximate (numerical Hessian) standard error: NaN
##     z-value: NaN, p-value: NA
## Wald statistic: NaN, p-value: NA
## 
## Log likelihood: NaN for mixed model
## ML residual variance (sigma squared): 0, (sigma: 0)
## Number of observations: 3 
## Number of parameters estimated: 5 
## AIC: NaN, (AIC for lm: -Inf)

El modelo de regresión espacial Durbin (SDM), que incorpora tanto efectos directos como indirectos (lag de las variables explicativas en zonas vecinas), fue estimado con las variables dummy por zona como predictores. Sin embargo, debido al reducido número de observaciones (solo tres zonas), el modelo presentó múltiples problemas técnicos.

Dos de las variables (lag.zona_norte y lag.zona_sur) fueron eliminadas automáticamente del modelo por colinealidad, es decir, porque eran combinaciones lineales de otras variables ya presentes. Además, los errores estándar, los valores z y p-values no pudieron calcularse (NaN), lo que impide evaluar la significancia estadística de los coeficientes. El parámetro rho, que mide la dependencia espacial, fue estimado como 1, lo cual en un contexto real indicaría dependencia perfecta, pero en este caso refleja un sobreajuste extremo.

En conclusión, aunque el modelo SDM fue ejecutado, sus resultados no pueden ser interpretados ni utilizados de forma confiable. Esta situación refuerza la necesidad de contar con un mayor número de unidades espaciales para que los modelos espaciales produzcan resultados estadísticamente válidos.

Selección

Dado que los tres modelos espaciales (SAR, SEM y SDM) fueron probados con base en las variables disponibles, el modelo SEM (Spatial Error Model) fue el único que logró completarse con un proceso de estimación relativamente estable, sin errores fatales. Aunque los resultados no son estadísticamente significativos por el número muy reducido de observaciones (3 zonas), el modelo SEM resultó el más funcional técnicamente.

Este modelo asume que la correlación espacial no afecta directamente la variable dependiente, sino que está presente en los errores no observados del modelo. En este caso, se observó un valor de lambda = -0.60, lo que sugiere una posible correlación espacial negativa en los residuos. A pesar de la falta de significancia estadística en los coeficientes, el modelo se considera el mejor entre los tres por ser el único que no presentó singularidades críticas ni exclusión de variables.

Reeepresentar en un mapa los valores estimados del modelo SEM para visualizar cómo se distribuye la cantidad estimada

library(leaflet)

# Agregar los valores ajustados al data frame
map_data$cantidad_estimado <- fitted(modelo_sem)
## This method assumes the response is known - see manual page
# Mapa interactivo
leaflet(map_data) %>%
  addTiles() %>%
  addCircleMarkers(
    lng = ~lng,
    lat = ~lat,
    radius = ~sqrt(cantidad_estimado) / 150,
    color = "darkgreen",
    label = ~paste(zona, "<br>Estimado:", round(cantidad_estimado)),
    fillOpacity = 0.7
  )

En el mapa se presentan los valores estimados de cantidad vendida por zona según el modelo SEM. Aunque las estimaciones coinciden exactamente con los valores reales debido al sobreajuste (varianza residual cero), este ejercicio permite visualizar cómo se distribuirían espacialmente los resultados del modelo. La visualización es útil como herramienta de comunicación, especialmente para respaldar decisiones de localización y estrategia comercial.

Propuestas

A partir de las limitaciones encontradas, se proponen las siguientes acciones para mejorar futuros análisis:

  • Mayor desagregación geográfica: Es fundamental contar con datos a nivel municipal o incluso por localidad. Esto incrementa el número de observaciones y permite que los modelos espaciales identifiquen patrones y clústers de manera robusta.

  • Incorporar variables socioeconómicas o logísticas: Agregar información como población, ingreso promedio, acceso a infraestructura, competencia regional, o ubicación de centros logísticos permitiría construir modelos más explicativos y realistas, con mayor poder predictivo.

Propuesta de localización de negocio

Con base en la demanda observada y la estimada en el modelo, se propone enfocar la estrategia de expansión en la zona Centro, que muestra la mayor cantidad vendida. Esta zona también se encuentra geográficamente bien conectada con otras regiones, lo que facilita la distribución.

  • Propuesta: Establecer un punto de distribución o centro logístico en el centro del país, preferentemente en la zona de Querétaro o Toluca, donde se combinan altos volúmenes de consumo, conectividad vial y cercanía con múltiples mercados.

Esta decisión no solo responde al comportamiento de la demanda, sino también a criterios logísticos y de escalabilidad, lo cual representa una oportunidad para optimizar operaciones y reducir costos de distribución.

Referencias