Quiz 3
Cesar Romeo Vega Rodriguez
Rodrigo Arroyo Oropeza
2025-05-12
Importar Librerías y Base de Datos
Librerías
### loading required libraries
library(foreign) # import external files
library(dplyr) # data manipulation
library(viridis) # offers several color palettes
library(RColorBrewer) # offers several color palettes
library(tidyverse) # it includes a collection of R packages designed for data science
library(jtools) # export regression summaries to tables
# spatial data analysis
library(sf) # functions to encode spatial vector data
library(tmap) # making maps so spatial data distributions are visualized
library(spdep) # a collection of functions to create spatial weight matrix
library(grid) # a set of functions and classes that represent graphical objects
library(rgeoda) # spatial data analysis based on GeoDa
library(tigris) # allows to work with shapefiles
library(regclass) # contains basic tools for visualizing, interpreting, and building regression models
library(spatialreg) # a collection of all the estimation functions for spatial cross-sectional models
library(tmap)
# GWR
library(spgwr) # Geographically Weighted Regression
library(GWmodel) # exploring spatial heterogeneity using Geographically Weighted modelsBases de Datos
map <- st_read("nl_mpios_auto_acc.shp")## Reading layer `nl_mpios_auto_acc' from data source
## `C:\Users\rodio\OneDrive\Escritorio\8semestre\Bloque Final\Quiz 3\nl_mpios_auto_acc.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 51 features and 22 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -101.2068 ymin: 23.16268 xmax: -98.42158 ymax: 27.79914
## Geodetic CRS: WGS 84data <- read_sf("nl_mpios_auto_acc.shp")1) Análisis Exploratorio de los Datos (EDA)
#remotes::install_github('r-tmap/tmap')
#tmaptools::palette_explorer() # palette of colors
#tm_shape(data) +
# tm_polygons("densidad_p", palette = "Purples", style="quantile", title="NL Population Density across Mpios") +
#tmap_mode("view")
#tmap_last()
dataa <- data |>
mutate(
## variables categóricas (niveles ordenados explícitamente)
mpio = factor(mpio), # 51 niveles – no cambia pero queda como factor
zona_urb = factor(zona_urb,
levels = c("interseccion", "no_interseccion", "suburbana")),
zona_subur = factor(zona_subur,
levels = c("camino_rural", "carretera_estatal","otro_camino", "urbana")),
aliento = factor(aliento,
levels = c("negativo", "positivo", "se_ignora")),
cinturon = factor(cinturon,
levels = c("cinturon_no", "cinturon_si", "se_ignora")),
region = factor(region,
levels = c("Norte", "Este", "Oeste", "Sur", "ZMM")),
sexo_lbl = factor(sexo, levels = c(1, 2),
labels = c("Hombres", "Mujeres")),
# etiqueta legible para el cinturón
cinturon_lbl = factor(cinturon,
levels = c("cinturon_si", "cinturon_no", "se_ignora"),
labels = c("Sí usa", "No usa", "Se ignora")),
## variable continua transformada (sólo la usaremos en histogramas y ESDA)
log_tasa = log(tasa_auto_), # sin ceros; min = 0.15
log_acc = log(auto_accid)
)a. Promedio y mediana de accidentes automovilísticos a nivel estatal y por región del estado.
hist(dataa$tasa_auto_, prob=TRUE, main = "Histogram Ratio de Accidentes de Coches per 10,000 Pop", xlab = "Car Accidents")
hist(dataa$log_tasa, prob=TRUE, main = "Histogram Ratio de Accidentes de Coches per 10,000 Pop", xlab = "log Car Accidents")
hist(dataa$auto_accid, prob=TRUE, main = "Cantidad de Accidentes", xlab = "Car Accidents")
hist(dataa$log_acc, prob=TRUE, main = "Cantidad de Accidentes", xlab = "log Car Accidents")
Como se puede apreciar, en los municipios de nuevo León existe una gran disparidad en la cantidad de accidentes automovilisticos teniendo regiones con una muy alta concentración de accidentes y otras regiones con muy pocos accidentes. Por lo que se deberá utilizar la transformación logaritmica para obtener un resultado con mayor interpretación y con una distribución más normal.
eda_tbl <- dataa |>
st_drop_geometry() |>
group_by(region) |>
summarise(
media_region_tasa = mean(tasa_auto_, na.rm = TRUE),
mediana_region_tasa = median(tasa_auto_, na.rm = TRUE),
media_region_cantidad = mean(auto_accid, na.rm = TRUE),
mediana_region_cantidad = median(auto_accid, na.rm = TRUE),
n_mpios = n()
) |>
ungroup()
media_estatal <- mean(dataa$tasa_auto_, na.rm = TRUE)
mediana_estatal <- median(dataa$tasa_auto_, na.rm = TRUE)
print(eda_tbl)## # A tibble: 5 × 6
## region media_region_tasa mediana_region_tasa media_region_cantidad
## <fct> <dbl> <dbl> <dbl>
## 1 Norte 2872. 167. 965.
## 2 Este 413. 325. 198.
## 3 Oeste 9095. 66.7 2055.
## 4 Sur 78.1 28.1 79.9
## 5 ZMM 168. 84.6 8551.
## # ℹ 2 more variables: mediana_region_cantidad <dbl>, n_mpios <int>cat("\nMedia estatal:", round(media_estatal, 2),
"\nMediana estatal:", round(mediana_estatal, 2), "\n")##
## Media estatal: 2737.38
## Mediana estatal: 73.32Como se puede apreciar la tasa promedio de accidentes por región sucede en las regiones oeste y norte, sin embargo esto se ve influenciado por la densidad poblacional, ya que como se aprecia en la cantidad promedio de accidentes de la región Metropolitana, este supera todas las regiones.
b. Dispersión regional de accidentes automovilísticos.
ggplot(dataa, aes(region, log_tasa, fill = region)) +
geom_boxplot(alpha = 0.85, outlier.colour = "red") +
scale_y_continuous("log(Tasa de accidentes)") +
scale_fill_brewer(palette = "Set3", guide = "none") +
ggtitle("Dispersión regional de accidentes (escala log)") +
theme_minimal(base_size = 12)
Tal como se vió en la tabla anterior, existe una diferencia notable en la dispersión en la variable de tasas de accidentes. Por lo que un analisis geoespacial podría ser útil para verificar si hay diferencias notables entre regiones o los mismos municipios.
c. Dispersión de accidentes automovilísticos por uso de cinturón y/o sexo.
ggplot(dataa, aes(cinturon_lbl, log_tasa, fill = sexo_lbl)) +
geom_boxplot(alpha = 0.80,
position = position_dodge(width = 0.8),
outlier.colour = "red") +
scale_fill_viridis_d("Sexo") +
scale_y_continuous("log(Tasa de accidentes)") +
xlab("Uso de cinturón") +
ggtitle("Accidentes por uso de cinturón y sexo (log)") +
theme_minimal(base_size = 12)
Una de las variables más notorias es la variable sexo y el uso del cinturon, pues si revisamos su despersión en el boxplot anterior, se puede identificar que la tasa de accidentes es mayor en mujeres, sin embargo en los hombres hay un comportamiento más aleatorio entre si se usó o no el cinturón
2) Análisis Exploratorio Espacial de los Datos (ESDA)
a. Calcular y mostrar la matriz de conectividad “Rook”.
swm_rook <- poly2nb(map, queen = FALSE)
rswm_rook <- nb2listw(swm_rook, style = "W", zero.policy = TRUE)
map_a <- as(map, "Spatial")
map_a_centroid <- coordinates(map_a)
plot(map_a,border="blue",axes=FALSE,las=1, main="SWM Rook - Nuevo Leon Mpios")
plot(map_a,col="grey",border=grey(0.9),axes=T,add=T)
plot(rswm_rook,coords=map_a_centroid,pch=19,cex=0.1,col="red",add=T) 
b. Identificación de clúster global de accidentes automovilísticos.
moran_global <- moran.test(dataa$log_tasa, rswm_rook, zero.policy = TRUE)
moran_global # se imprime en el Rmd##
## Moran I test under randomisation
##
## data: dataa$log_tasa
## weights: rswm_rook
##
## Moran I statistic standard deviate = 2.029, p-value = 0.02123
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.162597771 -0.020000000 0.008099056c. Identificación de clúster local de accidentes automovilísticos.
localI <- localmoran(dataa$log_tasa, rswm_rook, zero.policy = TRUE)
# identificar automáticamente la columna p‑valor
p_col <- grep("^Pr", colnames(localI), value = TRUE)[1]
dataa$Ii <- localI[, "Ii"]
dataa$P_Ii <- localI[, p_col]
mean_log <- mean(dataa$log_tasa, na.rm = TRUE)
lag_log <- lag.listw(rswm_rook, dataa$log_tasa)
dataa$quadrant <- with(dataa,
ifelse(log_tasa >= mean_log & lag_log >= mean_log & P_Ii <= 0.05, "Alto‑Alto",
ifelse(log_tasa < mean_log & lag_log < mean_log & P_Ii <= 0.05, "Bajo‑Bajo",
ifelse(log_tasa >= mean_log & lag_log < mean_log & P_Ii <= 0.05, "Alto‑Bajo",
ifelse(log_tasa < mean_log & lag_log >= mean_log & P_Ii <= 0.05, "Bajo‑Alto",
"No significativo")))))
# mapa LISA inline
lisa_cols <- c("Alto‑Alto"="#BB4444","Bajo‑Bajo"="#3E8FC5",
"Alto‑Bajo"="#E9842C","Bajo‑Alto"="#66C2A5",
"No significativo"="grey90")
tm_shape(dataa) +
tm_fill("quadrant", palette = lisa_cols, title = "Clúster LISA (log)") +
tm_borders(col = "white") +
tm_layout(frame = FALSE, legend.outside = TRUE)
# lista de municipios en clúster Alto‑Alto
mpios_HH <- dataa %>% filter(quadrant == "Alto‑Alto") %>% pull(mpio)
mpios_HH # se imprime en el Rmd## [1] China Los Aldamas Paras Los Herreras
## [5] General Treviño
## 51 Levels: Abasolo Agualeguas Allende Anahuac Apodaca Aramberri ... Villaldamad. ¿Cuáles son los municipios que integran el clúster / los clústeres de accidentes de autos en el estado de NL?
# Tabla con todos los municipios por cuadrante
cluster_mpios <- dataa |>
st_drop_geometry() |>
group_by(quadrant) |>
summarise(
n_mpios = n(),
municipios = paste(sort(mpio), collapse = ", ")
) |>
ungroup()
cluster_mpios # se imprime completo en el R Markdown## # A tibble: 5 × 3
## quadrant n_mpios municipios
## <chr> <int> <chr>
## 1 Alto‑Alto 5 China, General Treviño, Los Aldamas, Los Herreras, P…
## 2 Alto‑Bajo 2 General Teran, Rayones
## 3 Bajo‑Alto 2 Los Ramones, Melchor Ocampo
## 4 Bajo‑Bajo 1 Allende
## 5 No significativo 41 Abasolo, Agualeguas, Anahuac, Apodaca, Aramberri, Bu…# Si solo quieres los Alto‑Alto, por ejemplo:
mpios_HH <- cluster_mpios |> dplyr::filter(quadrant == "Alto‑Alto") |> pull(municipios)
cat("\nMunicipios clúster Alto‑Alto:\n", mpios_HH, "\n")##
## Municipios clúster Alto‑Alto:
## China, General Treviño, Los Aldamas, Los Herreras, ParasComo se puede observar, China, General Treviño, Los Aldamas, Los Herreras y Paras son los municipios con mas álto valor de autocorrelación espacial positivo cuentan. Por otro lado, General Terán, Allende, Los Ramones y Melchor Ocampo son otro cluster, un poco más separado con un valor de autocorrelación esoacial negativa.
3) Estimación de Modelos de Predicción Global
a. Modelo de regresión no espacial
# Modelo de regresión no espacial
model_NS <- lm(log(tasa_auto_) ~ gini + edad + log(pop) + grado_educ + cinturon + aliento + sexo, data = dataa)
summary(model_NS)##
## Call:
## lm(formula = log(tasa_auto_) ~ gini + edad + log(pop) + grado_educ +
## cinturon + aliento + sexo, data = dataa)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0750 -0.9253 -0.0174 1.1204 3.2702
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.95603 5.11587 -1.360 0.1814
## gini 13.41848 6.58086 2.039 0.0479 *
## edad 0.04889 0.03064 1.595 0.1183
## log(pop) -0.38670 0.21632 -1.788 0.0812 .
## grado_educ 0.68459 0.34758 1.970 0.0557 .
## cinturoncinturon_si 1.22477 0.66756 1.835 0.0738 .
## cinturonse_ignora -0.59691 0.72790 -0.820 0.4169
## alientopositivo 1.11685 0.67000 1.667 0.1031
## alientose_ignora -1.08824 1.26799 -0.858 0.3958
## sexo 0.72277 0.66045 1.094 0.2802
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.756 on 41 degrees of freedom
## Multiple R-squared: 0.4609, Adjusted R-squared: 0.3426
## F-statistic: 3.895 on 9 and 41 DF, p-value: 0.001229# Valores Reales
valores_reales <- dataa$tasa_auto_
# Predicciones del modelo
predicciones_NS <- exp(predict(model_NS))
# Calcular el RMSE y AIC
mse_NS <- mean((valores_reales - predicciones_NS)^2)
RMSE_NS <- sqrt(mse_NS)
AIC_NS <- AIC(model_NS)b. Diagnóstico de multicolinealidad y heterocedasticidad
# Multicolinealidad
VIF(model_NS)## GVIF Df GVIF^(1/(2*Df))
## gini 2.074122 1 1.440181
## edad 1.260473 1 1.122708
## log(pop) 2.943965 1 1.715799
## grado_educ 3.778599 1 1.943862
## cinturon 1.735424 2 1.147761
## aliento 1.775338 2 1.154304
## sexo 1.137351 1 1.066467Con los resultados del Variance Inflation Factor, ninguna de las variables independientes cuenta con un muy alto grado de autocorrelación; por lo que no existe multicolinealidad en los datos.
# Heterocedasticidad
lmtest::bptest(model_NS)##
## studentized Breusch-Pagan test
##
## data: model_NS
## BP = 8.4088, df = 9, p-value = 0.4935Conforme al resultado de la prueba, se obtuvo un valor mayor al 0.05 por lo que no se puede decir que no existe heterocedasticidad en los datos seleccionados en el modelo.
c. SAR
# SAR - Spatial Autoregressive Model
model_SAR <- lagsarlm(log(tasa_auto_) ~ gini + edad + log(pop) + grado_educ + cinturon + aliento + sexo, data = dataa, listw = rswm_rook)
summary(model_SAR)##
## Call:lagsarlm(formula = log(tasa_auto_) ~ gini + edad + log(pop) +
## grado_educ + cinturon + aliento + sexo, data = dataa, listw = rswm_rook)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.16407 -0.94335 -0.21126 1.06519 3.31322
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -7.852421 4.603028 -1.7059 0.08802
## gini 12.855179 5.875876 2.1878 0.02868
## edad 0.052923 0.027137 1.9502 0.05115
## log(pop) -0.322374 0.195326 -1.6504 0.09885
## grado_educ 0.624198 0.307863 2.0275 0.04261
## cinturoncinturon_si 1.228917 0.590680 2.0805 0.03748
## cinturonse_ignora -0.625627 0.645280 -0.9695 0.33227
## alientopositivo 1.145120 0.598772 1.9124 0.05582
## alientose_ignora -0.927115 1.135066 -0.8168 0.41405
## sexo 0.774951 0.584937 1.3248 0.18522
##
## Rho: 0.16052, LR test value: 1.0596, p-value: 0.30331
## Asymptotic standard error: 0.16254
## z-value: 0.98755, p-value: 0.32337
## Wald statistic: 0.97525, p-value: 0.32337
##
## Log likelihood: -94.98281 for lag model
## ML residual variance (sigma squared): 2.4139, (sigma: 1.5537)
## Number of observations: 51
## Number of parameters estimated: 12
## AIC: 213.97, (AIC for lm: 213.03)
## LM test for residual autocorrelation
## test value: 0.058093, p-value: 0.80954# Predicciones del modelo
predicciones_SAR <- exp(predict(model_SAR))
# Calcular el RMSE y AIC
mse_SAR <- mean((valores_reales - predicciones_SAR)^2)
RMSE_SAR <- sqrt(mse_SAR)
AIC_SAR <- AIC(model_SAR)d. SEM
model_SEM <- errorsarlm(log(tasa_auto_) ~ gini + edad + log(pop) + grado_educ + cinturon + aliento + sexo, data = dataa, listw = rswm_rook)
summary(model_SEM)##
## Call:errorsarlm(formula = log(tasa_auto_) ~ gini + edad + log(pop) +
## grado_educ + cinturon + aliento + sexo, data = dataa, listw = rswm_rook)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.98958 -0.93837 -0.23086 0.99236 3.26822
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.839982 4.702474 -1.4545 0.14579
## gini 13.663502 6.178335 2.2115 0.02700
## edad 0.053252 0.026488 2.0104 0.04439
## log(pop) -0.331270 0.206253 -1.6061 0.10825
## grado_educ 0.576526 0.322356 1.7885 0.07370
## cinturoncinturon_si 1.255976 0.584082 2.1503 0.03153
## cinturonse_ignora -0.569830 0.650700 -0.8757 0.38118
## alientopositivo 0.940593 0.584813 1.6084 0.10776
## alientose_ignora -1.088079 1.100181 -0.9890 0.32266
## sexo 0.785695 0.578762 1.3575 0.17461
##
## Lambda: 0.20081, LR test value: 0.85901, p-value: 0.35402
## Asymptotic standard error: 0.1861
## z-value: 1.079, p-value: 0.28058
## Wald statistic: 1.1643, p-value: 0.28058
##
## Log likelihood: -95.08309 for error model
## ML residual variance (sigma squared): 2.4155, (sigma: 1.5542)
## Number of observations: 51
## Number of parameters estimated: 12
## AIC: 214.17, (AIC for lm: 213.03)# Predicciones del modelo
predicciones_SEM <- exp(predict(model_SEM))
# Calcular el RMSE y AIC
mse_SEM <- mean((valores_reales - predicciones_SEM)^2)
RMSE_SEM <- sqrt(mse_SEM)
AIC_SEM <- AIC(model_SEM)e. SDM
model_DURB <- lagsarlm(log(tasa_auto_) ~ gini + edad + log(pop) + grado_educ + cinturon + aliento + sexo, data = dataa, listw = rswm_rook, type="mixed")
summary(model_DURB)##
## Call:lagsarlm(formula = log(tasa_auto_) ~ gini + edad + log(pop) +
## grado_educ + cinturon + aliento + sexo, data = dataa, listw = rswm_rook,
## type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.67482 -0.76156 -0.12165 1.01951 2.44957
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -22.201686 11.146273 -1.9918 0.04639
## gini 16.648601 7.080650 2.3513 0.01871
## edad 0.059199 0.028595 2.0703 0.03843
## log(pop) -0.379936 0.238139 -1.5954 0.11061
## grado_educ 0.343210 0.403150 0.8513 0.39459
## cinturoncinturon_si 0.278806 0.678404 0.4110 0.68109
## cinturonse_ignora -1.134573 0.683844 -1.6591 0.09709
## alientopositivo 1.043866 0.610573 1.7096 0.08733
## alientose_ignora -0.773761 1.139986 -0.6787 0.49730
## sexo 0.327764 0.572940 0.5721 0.56727
## lag.gini 22.851090 13.645754 1.6746 0.09401
## lag.edad 0.061754 0.085633 0.7211 0.47082
## lag.log(pop) -0.644146 0.419995 -1.5337 0.12510
## lag.grado_educ 2.014882 0.818958 2.4603 0.01388
## lag.cinturoncinturon_si -1.676897 2.440245 -0.6872 0.49197
## lag.cinturonse_ignora -1.857984 1.791394 -1.0372 0.29966
## lag.alientopositivo 2.682559 1.611342 1.6648 0.09595
## lag.alientose_ignora -1.765826 2.335457 -0.7561 0.44959
## lag.sexo -2.371249 1.712656 -1.3845 0.16619
##
## Rho: -0.27439, LR test value: 1.4339, p-value: 0.23112
## Asymptotic standard error: 0.20193
## z-value: -1.3589, p-value: 0.17419
## Wald statistic: 1.8465, p-value: 0.17419
##
## Log likelihood: -87.5661 for mixed model
## ML residual variance (sigma squared): 1.7874, (sigma: 1.3369)
## Number of observations: 51
## Number of parameters estimated: 21
## AIC: 217.13, (AIC for lm: 216.57)
## LM test for residual autocorrelation
## test value: 12.862, p-value: 0.0003353# Predicciones del modelo
predicciones_DURB <- exp(predict(model_DURB))
# Calcular el RMSE y AIC
mse_DURB <- mean((valores_reales - predicciones_DURB)^2)
RMSE_DURB <- sqrt(mse_DURB)
AIC_DURB <- AIC(model_DURB)f. Selección de Modelo de Predicción
table1 <- data.frame(Variable = c("No Espacial", "SAR", "SEM", "DURB"),
RMSE = c(RMSE_NS, RMSE_SAR, RMSE_SEM, RMSE_DURB),
AIC = c(AIC_NS,AIC_SAR,AIC_SEM,AIC_DURB))
table1## Variable RMSE AIC
## 1 No Espacial 13686.78 213.0252
## 2 SAR 13578.19 213.9656
## 3 SEM 13540.66 214.1662
## 4 DURB 12985.52 217.1322Conforme a la tabla comparativa anterior el mejor modelo de predicción fue el DURB al contar con un menor RMSE y un AIC menor que los demás modelos. Como se puede observar, la mayoría de los modelos cuenta con un RMSE similar, sin embargo destaca el DURBIN por su menor RMSE, apesar de tener un mayor AIC.
Las variables más significativas para este modelo resultó ser el gini, la edad y el resago del grado de educación con coeficientes positivos.
g. Identificación y mostrar la independencia de εi del modelo de predicción seleccionado.
#Generar variable que contenga los residuales de los modelos seleccionados
data2 <- dataa
data2$NS <- model_NS$residuals
data2$SAR <- model_SAR$residuals
data2$SEM <- model_SEM$residuals
data2$SDM <- model_DURB$residualsmoran.test(data2$SDM, rswm_rook, zero.policy = TRUE)##
## Moran I test under randomisation
##
## data: data2$SDM
## weights: rswm_rook
##
## Moran I statistic standard deviate = -0.91303, p-value = 0.8194
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.104461380 -0.020000000 0.008557436# Generar mapas de distribución de residuales
map1 <- tm_shape(data2) +
tm_polygons(
fill = "NS",
fill.scale = tm_scale_intervals(style = "quantile", n = 8, values = "brewer.yl_or_rd"),
fill.legend = tm_legend(title = "Residuals")
) +
tm_title("Non-Spatial - Residuals", position = c("right", "top"))
map2 <- tm_shape(data2) +
tm_polygons(
fill = "SAR",
fill.scale = tm_scale_intervals(style = "quantile", n = 8, values = "brewer.yl_or_rd"),
fill.legend = tm_legend(title = "Residuals")
) +
tm_title("SAR - Residuals", position = c("right", "top"))
map3 <- tm_shape(data2) +
tm_polygons(
fill = "SEM",
fill.scale = tm_scale_intervals(style = "quantile", n = 8, values = "brewer.yl_or_rd"),
fill.legend = tm_legend(title = "Residuals")
) +
tm_title("SEM - Residuals", position = c("right", "top"))
map4 <- tm_shape(data2) +
tm_polygons(
fill = "SDM",
fill.scale = tm_scale_intervals(style = "quantile", n = 8, values = "brewer.yl_or_rd"),
fill.legend = tm_legend(title = "Residuals")
) +
tm_title("DURB - Residuals", position = c("right", "top"))
tmap_arrange(map1, map2, map3, map4, ncol = 2)
Como se puede apreciar en todos los modelos, existe una independiencia en los residuales ya que no hay zonas donde exista una alta concentración de valores similares. Si bien sí existe poca autocorrelación espacial negativa entre los residuales este no es significativo.
h. Visualizar los valores estimados a nivel municipal de la variable dependiente del modelo de predicción seleccionado.
#Mapeo de la variable Dependiente
dataa$pred <- predicciones_DURB
tm_shape(dataa) +
tm_polygons(
fill = "pred",
fill.scale = tm_scale_intervals(style = "quantile", n = 8, values = "brewer.yl_or_rd"),
fill.legend = tm_legend(title = "MXN")
) +
tm_title("Predicciones de Accidentes Automovilisticos", position = c("right", "top"))
pred <- data.frame( Municipio = dataa$mpio,
Predicción = dataa$pred,
Original = dataa$auto_accid)
pred## Municipio Predicción.fit Predicción.trend Predicción.signal
## 1 Guadalupe 241.793753 6.785913 -1.2978282
## 2 Juarez 35.291101 4.756476 -1.1928450
## 3 Pesqueria 24.438155 4.240524 -1.0443779
## 4 Santiago 7.242709 2.869478 -0.8894822
## 5 Garcia 23.060772 4.533273 -1.3951398
## 6 Santa Catarina 521.605926 7.619614 -1.3627021
## 7 Monterrey 69.203870 5.438983 -1.2019263
## 8 San Pedro Garza Garcia 304.809484 7.057185 -1.3374978
## 9 El Carmen 89.328411 5.656109 -1.1637894
## 10 General Escobedo 291.067150 6.909590 -1.2360359
## 11 Salinas Victoria 51.907777 5.261941 -1.3124723
## 12 General Zuazua 55.256137 5.160925 -1.1489456
## 13 Apodaca 11.341705 3.630420 -1.2019337
## 14 Cienega de Flores 46.187948 5.068374 -1.2356547
## 15 Hidalgo 100.168243 5.876989 -1.2701377
## 16 Abasolo 354.670774 7.020435 -1.1492452
## 17 China 87.446482 6.318178 -1.8471511
## 18 General Bravo 113.945419 6.220788 -1.4850683
## 19 Los Aldamas 2034.201493 9.660367 -2.0425085
## 20 Doctor Coss 691.718020 8.131998 -1.5928196
## 21 Villaldama 140.927735 5.996189 -1.0479419
## 22 Sabinas Hidalgo 25.293551 4.720246 -1.4896960
## 23 San Nicolas de los Garza 71.099978 5.568673 -1.3045860
## 24 Bustamante 85.384569 5.738916 -1.2917511
## 25 Linares 20.246618 3.918356 -0.9103679
## 26 Vallecillo 2719.712113 8.788426 -0.8801446
## 27 Rayones 213.285633 5.978482 -0.6158495
## 28 Los Ramones 208.694604 7.024789 -1.6839167
## 29 Paras 280.357796 7.730336 -2.0942696
## 30 Montemorelos 1.267474 1.322024 -1.0849976
## 31 Mina 265.871951 6.925841 -1.3428260
## 32 Mier y Noriega 4.047833 2.409547 -1.0113651
## 33 Melchor Ocampo 87.696557 7.023056 -2.5491734
## 34 Marin 70.747894 5.706864 -1.4477417
## 35 Lampazos de Naranjo 136.467470 6.355450 -1.4393634
## 36 Iturbide 10.676658 3.474111 -1.1060511
## 37 Hualahuises 17.493782 3.675344 -0.8134981
## 38 Higueras 183.967301 6.586719 -1.3719614
## 39 Los Herreras 8502.189481 10.717837 -1.6697580
## 40 General Zaragoza 269.003390 6.566023 -0.9712990
## 41 General Treviño 370.440586 7.944515 -2.0298223
## 42 General Teran 247.400104 6.235181 -0.7241742
## 43 Galeana 43.205456 4.928347 -1.1623800
## 44 Doctor Gonzalez 515.649719 7.639907 -1.3944789
## 45 Doctor Arroyo 49.821201 4.892440 -0.9839991
## 46 Cerralvo 944.849254 8.387642 -1.5366169
## 47 Cadereyta Jimenez 12.471696 3.429393 -0.9059313
## 48 Aramberri 480.180241 7.234058 -1.0598964
## 49 Anahuac 27.074366 4.853020 -1.5544325
## 50 Allende 14.244878 3.092517 -0.4361196
## 51 Agualeguas 219.094984 7.075099 -1.6855933
## Original
## 1 7384
## 2 1167
## 3 555
## 4 125
## 5 2205
## 6 1445
## 7 42956
## 8 7828
## 9 212
## 10 4472
## 11 659
## 12 403
## 13 3178
## 14 340
## 15 106
## 16 10
## 17 126
## 18 179
## 19 179
## 20 88
## 21 179
## 22 10
## 23 14563
## 24 91
## 25 172
## 26 4675
## 27 176
## 28 31
## 29 13
## 30 1
## 31 590
## 32 11
## 33 8
## 34 19
## 35 2
## 36 13
## 37 12
## 38 27
## 39 19677
## 40 88
## 41 206
## 42 583
## 43 133
## 44 37
## 45 151
## 46 6228
## 47 310
## 48 53
## 49 98
## 50 24
## 51 46Como se puede apreciar, las predicciones muestran una alta concentración de accidentes cerca de la zona Metropolitana y Norte, mientras que en la región Sur y Oeste no existe mucha concentración de ellos.
4) Estimación de Modelos de Predicción Local
a. GWR
bw1 <- bw.gwr(log(tasa_auto_) ~ gini + edad + log(pop) + grado_educ + cinturon + aliento + sexo,
approach = "AIC", adaptive = T, data=data) ## Adaptive bandwidth (number of nearest neighbours): 39 AICc value: 247.5901
## Error in eval(expr, envir) : inv(): matrix is singular
## Adaptive bandwidth (number of nearest neighbours): 32 AICc value: Inf
## Adaptive bandwidth (number of nearest neighbours): 44 AICc value: 234.9491
## Adaptive bandwidth (number of nearest neighbours): 46 AICc value: 231.4571
## Adaptive bandwidth (number of nearest neighbours): 49 AICc value: 226.844
## Adaptive bandwidth (number of nearest neighbours): 49 AICc value: 226.844model_gwr <- gwr.basic(log(tasa_auto_) ~ gini + edad + log(pop) + grado_educ + cinturon + aliento + sexo,
adaptive = T, data = data, bw = bw1)
model_gwr## ***********************************************************************
## * Package GWmodel *
## ***********************************************************************
## Program starts at: 2025-05-12 10:44:51.798706
## Call:
## gwr.basic(formula = log(tasa_auto_) ~ gini + edad + log(pop) +
## grado_educ + cinturon + aliento + sexo, data = data, bw = bw1,
## adaptive = T)
##
## Dependent (y) variable: tasa_auto_
## Independent variables: gini edad pop grado_educ cinturon aliento sexo
## Number of data points: 51
## ***********************************************************************
## * Results of Global Regression *
## ***********************************************************************
##
## Call:
## lm(formula = formula, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0750 -0.9253 -0.0174 1.1204 3.2702
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.95603 5.11587 -1.360 0.1814
## gini 13.41848 6.58086 2.039 0.0479 *
## edad 0.04889 0.03064 1.595 0.1183
## log(pop) -0.38670 0.21632 -1.788 0.0812 .
## grado_educ 0.68459 0.34758 1.970 0.0557 .
## cinturoncinturon_si 1.22477 0.66756 1.835 0.0738 .
## cinturonse_ignora -0.59691 0.72790 -0.820 0.4169
## alientopositivo 1.11685 0.67000 1.667 0.1031
## alientose_ignora -1.08824 1.26799 -0.858 0.3958
## sexo 0.72277 0.66045 1.094 0.2802
##
## ---Significance stars
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Residual standard error: 1.756 on 41 degrees of freedom
## Multiple R-squared: 0.4609
## Adjusted R-squared: 0.3426
## F-statistic: 3.895 on 9 and 41 DF, p-value: 0.001229
## ***Extra Diagnostic information
## Residual sum of squares: 126.4113
## Sigma(hat): 1.606182
## AIC: 213.0252
## AICc: 219.7944
## BIC: 226.5254
## ***********************************************************************
## * Results of Geographically Weighted Regression *
## ***********************************************************************
##
## *********************Model calibration information*********************
## Kernel function: bisquare
## Adaptive bandwidth: 49 (number of nearest neighbours)
## Regression points: the same locations as observations are used.
## Distance metric: Euclidean distance metric is used.
##
## ****************Summary of GWR coefficient estimates:******************
## Min. 1st Qu. Median 3rd Qu. Max.
## Intercept -1.0897e+01 -1.0044e+01 -9.0528e+00 -7.9402e+00 -3.4586
## gini 9.4842e+00 1.5475e+01 1.7160e+01 1.8317e+01 20.0551
## edad -5.3734e-04 5.4904e-02 6.6064e-02 6.8929e-02 0.0778
## log(pop) -4.5054e-01 -2.5649e-01 -2.3277e-01 -2.2359e-01 -0.2110
## grado_educ 5.4334e-01 5.6625e-01 5.7728e-01 6.1741e-01 0.7933
## cinturoncinturon_si 5.2726e-01 1.3058e+00 1.4440e+00 1.5178e+00 1.6677
## cinturonse_ignora -6.8257e-01 -5.1685e-01 -4.2145e-01 -2.1443e-01 0.1108
## alientopositivo 5.9780e-01 8.2391e-01 1.0271e+00 1.2102e+00 1.7917
## alientose_ignora -2.4305e+00 -1.9115e+00 -1.7438e+00 -1.5416e+00 -1.2247
## sexo 3.6618e-01 4.7065e-01 5.4081e-01 6.9996e-01 1.0052
## ************************Diagnostic information*************************
## Number of data points: 51
## Effective number of parameters (2trace(S) - trace(S'S)): 18.65232
## Effective degrees of freedom (n-2trace(S) + trace(S'S)): 32.34768
## AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 226.844
## AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 191.2964
## BIC (GWR book, Fotheringham, et al. 2002,GWR p. 61, eq. 2.34): 186.4029
## Residual sum of squares: 93.36351
## R-square value: 0.6018355
## Adjusted R-square value: 0.3649218
##
## ***********************************************************************
## Program stops at: 2025-05-12 10:44:52.04059b. Visualizar variable dependiente del modelo de predicción local.
gwr_sf = st_as_sf(model_gwr$SDF)#Mapeo de la variable Dependiente
gwr_sf$y_predicted <- exp(gwr_sf$yhat)
tm_shape(gwr_sf) +
tm_polygons(
fill = "y_predicted",
fill.scale = tm_scale_intervals(style = "quantile", n = 8, values = "brewer.yl_or_rd"),
fill.legend = tm_legend(title = "MXN")
) +
tm_title("Tasa de Accidentes por Municipio", position = c("right", "top"))
Como se esperaba para el analisis, la zona Norte y Metropolitana cuenta con la mayor tasa de accidentes, dejando sus alrededores con un tasa menor. Algo a destacar es que únicamente hay un cluster notable que abarca 7 estados con una tasa de accidentes mayor al valor 378, los demás municipios cuentan con tasas distintas como para formar clusters.
c. Visualizar coeficiente R2 del modelo de predicción local.
#Mapeo de la variable Dependiente
gwr_sf$R2 <- gwr_sf$Local_R2
tm_shape(gwr_sf) +
tm_polygons(
fill = "R2",
fill.scale = tm_scale_intervals(style = "quantile", n = 8, values = "brewer.yl_or_rd"),
fill.legend = tm_legend(title = "MXN")
) +
tm_title("R2 por Municipio", position = c("right", "top"))
Por la parte del R2 vemos que la región Este es la que mayor R2 tiene, dejando los municipios del Oeste y los extremos del Norte y Sur con menor valor de R2.
5) Hallazgos identificados.
Hallazgos Análisis Exploratorio de los Datos (EDA)
- La zona Norte es la que tiene la mayor mediana y tiene un rango más amplio por el tamaño. Mientras que la zona Sur es la menor en ambos.
- En cinturon las mujeres representan una dispersion mas alta mientras que la de los hombres es menor
- Esta la existencia de muchos outliers en los datos por lo que se podria recomendar quitarlos o maniobrarlos
Hallazgos Análisis Exploratorio Espacial de los Datos (ESDA)
- La matriz Rook no deja ningun municipio aislado por lo que es posible que haya autocorrelacion espacial.
- El global Moran I = 0.163 (con un p = 0.02) por lo que se confirma el agrupamiento
- Los cinco municipios mas importantes fueron China, G. Treviño, Los Aldamas, Los Herreras y Parás
- Allende es el unico municipio Bajo-Bajo
- La mayoria de los estados quedo como “No significativo” por lo que no siguen ningun patron
Hallazgos Estimación de Modelos de Predicción Global
- El modelo DURBIN fue el mejor modelo estimado donde la variable gini, edad y el resago de la variable grad_educacion fueron las más significativas con un impacto positivo en la tasa de accidentes. Por lo que a mayor edad y valor de gini se aumenta la tasa de accidentes.
- El valor de Rho resultó negativo con un valor de -0.27, por lo que lo que sucede en los municipios cercanos, disminuye la tasa de accidentes en el municipio seleccionado.
- Una diferencia entre el modelo espacial y no espacial estimado es que la variable cinturon en el modelo no espacial cuando su valor es SI es significativa y positiva, sin embargo cuando esta entra al modelo espacial el resago de esta variable termina siendo negativa; implicando que valores más altos en municipios vecinos llegan a disminuir el impacto de la variable dependiente.
Hallazgos Estimación de Modelos de Predicción Local
- Se mejoro el ajute del modelo al pasar a GWR, el AIC se reduce a 213 y se eleva la R2 global
- Los valores más altos de la tasa_auto_ terminan ocurriendo en la ZMM
- Existe más explicacion en el sureste (R2 0.71) que en el norte (R2 0.68) por lo que hay mas precision en una region que otra
- La mayor tasa de accidentes se concentran en la región Noreste y en la Zona Metropolitana.