Quiz 3

########################
### AD3003B - Quiz 3 ###
###  Feb - June 2025 ###
########################

### Nombre: Jenaro Martinez
### Nombre: Fernanda Sanchez

### loading required libraries 
library(foreign)        # import external files
#install.packages("dplyr")
library(dplyr)          # data manipulation 

## 
## 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(viridis)        # offers several color palettes

## Loading required package: viridisLite

library(RColorBrewer)   # offers several color palettes 
library(tidyverse)      # it includes a collection of R packages designed for data science

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ ggplot2   3.5.1     ✔ stringr   1.5.1
## ✔ lubridate 1.9.4     ✔ tibble    3.2.1
## ✔ purrr     1.0.4     ✔ tidyr     1.3.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(jtools)         # export regression summaries to tables

# spatial data analysis 
library(sf)             # functions to encode spatial vector data 

## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUE

library(tmap)           # making maps so spatial data distributions are visualized 
library(spdep)          # a collection of functions to create spatial weight matrix  

## 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')`

library(grid)           # a set of functions and classes that represent graphical objects
library(rgeoda)         # spatial data analysis based on GeoDa 

## Loading required package: digest
## 
## Attaching package: 'rgeoda'
## 
## The following object is masked from 'package:spdep':
## 
##     skater

library(tigris)         # allows to work with shapefiles

## To enable caching of data, set `options(tigris_use_cache = TRUE)`
## in your R script or .Rprofile.

library(regclass)       # contains basic tools for visualizing, interpreting, and building regression models 

## Loading required package: bestglm
## Loading required package: leaps
## Loading required package: VGAM
## Loading required package: stats4
## Loading required package: splines
## Loading required package: rpart
## Loading required package: randomForest
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## 
## The following object is masked from 'package:ggplot2':
## 
##     margin
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## Important regclass change from 1.3:
## All functions that had a . in the name now have an _
## all.correlations -> all_correlations, cor.demo -> cor_demo, etc.

library(spatialreg)     # a collection of all the estimation functions for spatial cross-sectional models

## Loading required package: Matrix
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
## 
## 
## 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

# GWR
library(spgwr)          # Geographically Weighted Regression 

## Loading required package: sp
## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)

library(GWmodel)        # exploring spatial heterogeneity using Geographically Weighted models

## Loading required package: robustbase
## Loading required package: Rcpp
## Welcome to GWmodel version 2.4-2.
## 
## Attaching package: 'GWmodel'
## 
## The following object is masked from 'package:VGAM':
## 
##     AICc

### importing dataset 
# loading required shapefile 
datab <- st_read("/Users/jenaromtzg/Downloads/Quiz_3 2/nl_map/nl_mpios_auto_acc.shp")

## Reading layer `nl_mpios_auto_acc' from data source 
##   `/Users/jenaromtzg/Downloads/Quiz_3 2/nl_map/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 84

dataa <- read_sf("/Users/jenaromtzg/Downloads/Quiz_3 2/nl_map/nl_mpios_auto_acc.shp")
summary(dataa)

##     OBJECTID       CODELAG        CVE_ENT      IDUNICO        Shape_Leng    
##  Min.   : 1.0   Min.   :2199   Min.   :19   Min.   :19001   Min.   :0.2858  
##  1st Qu.:13.5   1st Qu.:2212   1st Qu.:19   1st Qu.:19014   1st Qu.:0.8586  
##  Median :26.0   Median :2224   Median :19   Median :19026   Median :1.5393  
##  Mean   :26.0   Mean   :2224   Mean   :19   Mean   :19026   Mean   :1.7867  
##  3rd Qu.:38.5   3rd Qu.:2236   3rd Qu.:19   3rd Qu.:19038   3rd Qu.:2.4495  
##  Max.   :51.0   Max.   :2249   Max.   :19   Max.   :19051   Max.   :4.7800  
##    Shape_Area         OBJECTID_1     IDUNICO_1         mpio          
##  Min.   :0.004224   Min.   : 1.0   Min.   :19001   Length:51         
##  1st Qu.:0.021224   1st Qu.:13.5   1st Qu.:19014   Class :character  
##  Median :0.064960   Median :26.0   Median :19026   Mode  :character  
##  Mean   :0.113032   Mean   :26.0   Mean   :19026                     
##  3rd Qu.:0.145082   3rd Qu.:38.5   3rd Qu.:19038                     
##  Max.   :0.630891   Max.   :51.0   Max.   :19051                     
##    auto_accid        tasa_auto_         zona_urb          zona_subur       
##  Min.   :    1.0   Min.   :    0.15   Length:51          Length:51         
##  1st Qu.:   34.0   1st Qu.:   31.66   Class :character   Class :character  
##  Median :  172.0   Median :   73.32   Mode  :character   Mode  :character  
##  Mean   : 2389.1   Mean   : 2737.38                                        
##  3rd Qu.:  624.5   3rd Qu.:  299.64                                        
##  Max.   :42956.0   Max.   :98483.48                                        
##       sexo         aliento            cinturon              edad      
##  Min.   :1.000   Length:51          Length:51          Min.   :15.38  
##  1st Qu.:2.000   Class :character   Class :character   1st Qu.:32.10  
##  Median :2.000   Mode  :character   Mode  :character   Median :36.35  
##  Mean   :1.804                                         Mean   :37.61  
##  3rd Qu.:2.000                                         3rd Qu.:41.46  
##  Max.   :2.000                                         Max.   :64.50  
##       pop            densidad_p            gini           region         
##  Min.   :   1071   Min.   :   0.830   Min.   :0.3000   Length:51         
##  1st Qu.:   4182   1st Qu.:   4.245   1st Qu.:0.3450   Class :character  
##  Median :  15902   Median :   7.070   Median :0.3800   Mode  :character  
##  Mean   : 110003   Mean   : 498.641   Mean   :0.3922                     
##  3rd Qu.:  78173   3rd Qu.: 268.535   3rd Qu.:0.4400                     
##  Max.   :1124835   Max.   :4748.840   Max.   :0.4900                     
##    grado_educ              geometry 
##  Min.   : 6.590   POLYGON      :51  
##  1st Qu.: 7.945   epsg:4326    : 0  
##  Median : 8.840   +proj=long...: 0  
##  Mean   : 8.933                     
##  3rd Qu.: 9.645                     
##  Max.   :13.160

### Map - The State of Nuevo Leon at Municipal Level 
# remotes::install_github('r-tmap/tmap') # it allows to work with the library(tmap) so we can display the maps. 

#install.packages("shiny")''
#install.packages("shinyjs")

Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.

##1) Análisis Exploratorio de los Datos (EDA)


#a. Promedio y mediana de accidentes automovilísticos a nivel estatal y por región del estado.

promedio_estatal <- mean(datab$auto_accid, na.rm = TRUE)
mediana_estatal <- median(datab$auto_accid, na.rm = TRUE)

cat("Promedio de accidentes a nivel estatal:", promedio_estatal, "\n")

## Promedio de accidentes a nivel estatal: 2389.098

cat("Mediana de accidentes a nivel estatal:", mediana_estatal, "\n")

## Mediana de accidentes a nivel estatal: 172

# A nivel regional
resumen_region <- datab %>%
  group_by(region) %>%
  summarise(
    promedio_accidentes = mean(auto_accid, na.rm = TRUE),
    mediana_accidentes = median(auto_accid, na.rm = TRUE)
  )

print(resumen_region)

## Simple feature collection with 5 features and 3 fields
## Geometry type: GEOMETRY
## Dimension:     XY
## Bounding box:  xmin: -101.2068 ymin: 23.16268 xmax: -98.42158 ymax: 27.79914
## Geodetic CRS:  WGS 84
## # A tibble: 5 × 4
##   region promedio_accidentes mediana_accidentes                         geometry
##   <chr>                <dbl>              <dbl>                   <GEOMETRY [°]>
## 1 Este                 198.               152.  POLYGON ((-98.58157 25.72242, -…
## 2 Norte                965.                68.5 POLYGON ((-99.40973 26.31697, -…
## 3 Oeste               2055.               340   MULTIPOLYGON (((-99.80249 25.99…
## 4 Sur                   79.9               70.5 POLYGON ((-100.0582 25.3615, -1…
## 5 ZMM                 8551.              3825   POLYGON ((-99.9929 25.68308, -9…

# Visualización de histogramas
hist(datab$tasa_auto_, 
     prob = TRUE, 
     main = "Histograma: Tasa de Accidentes por cada 10,000 Hab.",
     xlab = "Tasa de Accidentes", 
     col = "skyblue", 
     border = "white")

hist(log(datab$tasa_auto_), 
     prob = TRUE, 
     main = "Histograma: log(Tasa de Accidentes por cada 10,000 Hab.)",
     xlab = "log(Tasa de Accidentes)", 
     col = "lightgreen", 
     border = "white")

#b. Dispersión regional de accidentes automovilísticos.

ggplot(datab, aes(x = region, y = auto_accid)) +
  geom_boxplot(fill = "#69b3a2") +
  labs(title = "Dispersión de Accidentes por Región",
       x = "Región", y = "Número de Accidentes") +
  theme_minimal()

#c. Dispersión de accidentes automovilísticos por uso de cinturón y/o sexo. 

# Por sexo

datab$sexo <- factor(datab$sexo, 
                     levels = c(1, 2), 
                     labels = c("Mujer", "Hombre"))


ggplot(datab, aes(x = sexo, y = auto_accid)) +
  geom_boxplot(fill = "#FFA07A") +
  labs(title = "Dispersión de Accidentes por Sexo",
       x = "Sexo", y = "Número de Accidentes") +
  theme_minimal()

# Por uso de cinturón

datab$cinturon <- case_when(
  grepl("cinturon_si", datab$cinturon) ~ "Sí",
  grepl("cinturon_no", datab$cinturon) ~ "No",
  TRUE ~ "Ignoro"  # Consideramos cualquier otro valor como 'Ignoro'
)

# Convertimos a factor con orden lógico
datab$cinturon <- factor(datab$cinturon, 
                         levels = c("Sí", "No", "Ignoro"))

# Boxplot por uso de cinturón de seguridad
ggplot(datab, aes(x = cinturon, y = auto_accid)) +
  geom_boxplot(fill = "#87CEFA") +
  labs(title = "Dispersión de Accidentes por Uso de Cinturón de Seguridad",
       x = "Uso de Cinturón", y = "Número de Accidentes") +
  theme_minimal()

#ZMM Zona Metropolitana de Monterrey

##2)Análisis Exploratorio Espacial de los Datos (ESDA)


#a. Calcular y mostrar la matriz de conectividad “Rook”.

datab_sp <- as(datab, "Spatial")

# Crear matriz de conectividad tipo Rook (sin vértices compartidos)
swm_rook <- poly2nb(datab_sp, queen = FALSE)
rswm_rook <- nb2listw(swm_rook, style = "W", zero.policy = TRUE)

# Visualizar la matriz de conectividad Rook
plot(datab_sp, border = "grey", main = "Matriz de Conectividad - Rook")
plot(swm_rook, coordinates(datab_sp), pch = 19, cex = 0.5, col = "red", add = TRUE)

#b. Identificación de clúster global de accidentes automovilísticos.

moran_global <- moran.test(datab$tasa_auto_, rswm_rook, zero.policy = TRUE)
print(moran_global)

## 
##  Moran I test under randomisation
## 
## data:  datab$tasa_auto_  
## weights: rswm_rook    
## 
## Moran I statistic standard deviate = 0.14327, p-value = 0.443
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.01516424       -0.02000000        0.00113921

# Interpretación rápida:
if (moran_global$p.value < 0.05) {
  cat("Existe autocorrelación espacial global significativa.\n")
} else {
  cat("No se detecta autocorrelación espacial global significativa.\n")
}

## No se detecta autocorrelación espacial global significativa.

#c. Identificación de clúster local de accidentes automovilísticos.

lisa_result <- localmoran(datab$tasa_auto_, rswm_rook, zero.policy = TRUE)

# Agregar resultados al shapefile
datab$lisa_I <- lisa_result[,1]
datab$lisa_pvalue <- lisa_result[,5]

# Clasificación de los clústeres
datab$lisa_cluster <- case_when(
  datab$lisa_pvalue > 0.05 ~ "No Significativo",
  datab$lisa_I > 0 ~ "High-High",
  datab$lisa_I < 0 ~ "Low-Low",
  TRUE ~ "Indeterminado"
)

# Visualización de clústeres LISA
library(tmap)

tm_shape(datab) +
  tm_fill("lisa_cluster", palette = "Set3", title = "LISA Clúster") +
  tm_borders() +
  tm_layout(main.title = "Clústeres Locales de Accidentes (LISA)", 
             legend.outside = TRUE)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
## visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Set3" is named
## "brewer.set3"
## Multiple palettes called "set3" found: "brewer.set3", "hcl.set3". The first one, "brewer.set3", is returned.
## 
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

#d. ¿Cuáles son los municipios que integran el clúster / los clústeres de accidentes de autos en el estado de NL? 

municipios_cluster <- datab %>%
  filter(lisa_cluster != "No Significativo") %>%
  select(mpio, lisa_cluster)

print(municipios_cluster)

## Simple feature collection with 6 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -99.97067 ymin: 25.06713 xmax: -98.42158 ymax: 26.32927
## Geodetic CRS:  WGS 84
##              mpio lisa_cluster                       geometry
## 1           China      Low-Low POLYGON ((-99.18902 25.9053...
## 2     Los Aldamas      Low-Low POLYGON ((-99.16702 26.0754...
## 3     Los Ramones      Low-Low POLYGON ((-99.5387 25.90161...
## 4  Melchor Ocampo      Low-Low POLYGON ((-99.39597 26.0949...
## 5 General Treviño      Low-Low POLYGON ((-99.33367 26.3063...
## 6        Cerralvo    High-High POLYGON ((-99.58946 26.2353...

print(municipios_cluster)

## Simple feature collection with 6 features and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -99.97067 ymin: 25.06713 xmax: -98.42158 ymax: 26.32927
## Geodetic CRS:  WGS 84
##              mpio lisa_cluster                       geometry
## 1           China      Low-Low POLYGON ((-99.18902 25.9053...
## 2     Los Aldamas      Low-Low POLYGON ((-99.16702 26.0754...
## 3     Los Ramones      Low-Low POLYGON ((-99.5387 25.90161...
## 4  Melchor Ocampo      Low-Low POLYGON ((-99.39597 26.0949...
## 5 General Treviño      Low-Low POLYGON ((-99.33367 26.3063...
## 6        Cerralvo    High-High POLYGON ((-99.58946 26.2353...

##3)Estimación de Modelos de Predicción Global

library(car)

## Loading required package: carData

## 
## Attaching package: 'car'

## The following object is masked from 'package:VGAM':
## 
##     logit

## The following object is masked from 'package:purrr':
## 
##     some

## The following object is masked from 'package:dplyr':
## 
##     recode

library(lmtest)

## Loading required package: zoo

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

## 
## Attaching package: 'lmtest'

## The following object is masked from 'package:VGAM':
## 
##     lrtest

#a. Modelo de regresión no espacial

modelo_ols <- lm(log(tasa_auto_) ~ zona_urb + densidad_p + region + edad + I(edad^2) + sexo, data = datab)
summary(modelo_ols)

## 
## Call:
## lm(formula = log(tasa_auto_) ~ zona_urb + densidad_p + region + 
##     edad + I(edad^2) + sexo, data = datab)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.4440 -0.6831 -0.0673  0.5823  5.4782 
## 
## Coefficients:
##                           Estimate Std. Error t value Pr(>|t|)  
## (Intercept)             -3.1457434  4.5174960  -0.696   0.4902  
## zona_urbno_interseccion  2.1542260  1.1543210   1.866   0.0694 .
## zona_urbsuburbana        1.1717921  1.4469273   0.810   0.4228  
## densidad_p               0.0002395  0.0004104   0.584   0.5628  
## regionNorte             -0.1875490  1.0322058  -0.182   0.8567  
## regionOeste             -1.0841957  1.0442633  -1.038   0.3054  
## regionSur               -2.6271709  1.0063100  -2.611   0.0127 *
## regionZMM               -1.7246836  1.4008889  -1.231   0.2255  
## edad                     0.3973415  0.2009280   1.978   0.0549 .
## I(edad^2)               -0.0048484  0.0023373  -2.074   0.0445 *
## sexoHombre               1.0935001  0.7360987   1.486   0.1452  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.985 on 40 degrees of freedom
## Multiple R-squared:  0.3277, Adjusted R-squared:  0.1596 
## F-statistic:  1.95 on 10 and 40 DF,  p-value: 0.06633

#b. Diagnóstico de multicolinealidad y heterocedasticidad

# Multicolinealidad (VIF)
vif(modelo_ols)

##                 GVIF Df GVIF^(1/(2*Df))
## zona_urb    2.265832  2        1.226894
## densidad_p  2.655803  1        1.629663
## region      3.193847  4        1.156217
## edad       42.393073  1        6.510996
## I(edad^2)  41.233301  1        6.421316
## sexo        1.105274  1        1.051320

# Heterocedasticidad (Breusch-Pagan Test)
bptest(modelo_ols)

## 
##  studentized Breusch-Pagan test
## 
## data:  modelo_ols
## BP = 9.8503, df = 10, p-value = 0.4537

# Comparar AIC
AIC(modelo_ols)

## [1] 226.2857

#c. SAR

sar_model <- lagsarlm(log(tasa_auto_) ~ zona_urb + densidad_p + region + edad + I(edad^2) + sexo, data = datab_sp, listw = rswm_rook, zero.policy = TRUE)

summary(sar_model)

## 
## Call:lagsarlm(formula = log(tasa_auto_) ~ zona_urb + densidad_p + 
##     region + edad + I(edad^2) + sexo, data = datab_sp, listw = rswm_rook, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.42717 -0.72601 -0.11451  0.56502  5.35978 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                            Estimate  Std. Error z value Pr(>|z|)
## (Intercept)             -4.02054407  4.30140151 -0.9347  0.34994
## zona_urbno_interseccion  2.13836151  1.01966606  2.0971  0.03598
## zona_urbsuburbana        1.07062822  1.27794902  0.8378  0.40216
## densidad_p               0.00023270  0.00036289  0.6412  0.52136
## regionNorte             -0.16421442  0.91304395 -0.1799  0.85727
## regionOeste             -0.98897408  0.92347551 -1.0709  0.28420
## regionSur               -2.40465644  0.96710590 -2.4864  0.01290
## regionZMM               -1.58871277  1.25385801 -1.2671  0.20513
## edad                     0.41333420  0.17941536  2.3038  0.02123
## I(edad^2)               -0.00502303  0.00208667 -2.4072  0.01608
## sexoHombre               1.12058278  0.65519712  1.7103  0.08721
## 
## Rho: 0.088749, LR test value: 0.20473, p-value: 0.65093
## Asymptotic standard error: 0.18627
##     z-value: 0.47645, p-value: 0.63375
## Wald statistic: 0.22701, p-value: 0.63375
## 
## Log likelihood: -101.0405 for lag model
## ML residual variance (sigma squared): 3.0734, (sigma: 1.7531)
## Number of observations: 51 
## Number of parameters estimated: 13 
## AIC: 228.08, (AIC for lm: 226.29)
## LM test for residual autocorrelation
## test value: 0.04346, p-value: 0.83486

#d. SEM

sem_model <- errorsarlm(log(tasa_auto_) ~ zona_urb + densidad_p + region + edad + I(edad^2) + sexo, data = datab_sp, listw = rswm_rook, zero.policy = TRUE)

summary(sem_model)

## 
## Call:errorsarlm(formula = log(tasa_auto_) ~ zona_urb + densidad_p + 
##     region + edad + I(edad^2) + sexo, data = datab_sp, listw = rswm_rook, 
##     zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -4.43443 -0.71716 -0.14081  0.56040  5.35240 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                            Estimate  Std. Error z value Pr(>|z|)
## (Intercept)             -3.32022990  3.91274949 -0.8486  0.39612
## zona_urbno_interseccion  2.09778004  1.02580494  2.0450  0.04085
## zona_urbsuburbana        1.01127284  1.26883143  0.7970  0.42544
## densidad_p               0.00024282  0.00037103  0.6544  0.51283
## regionNorte             -0.16080821  0.96887063 -0.1660  0.86818
## regionOeste             -0.92954906  0.97106694 -0.9572  0.33844
## regionSur               -2.55417720  0.95138326 -2.6847  0.00726
## regionZMM               -1.61981716  1.28033532 -1.2652  0.20582
## edad                     0.40086127  0.17339747  2.3118  0.02079
## I(edad^2)               -0.00485966  0.00201883 -2.4072  0.01608
## sexoHombre               1.07665416  0.64444899  1.6707  0.09479
## 
## Lambda: 0.094779, LR test value: 0.14084, p-value: 0.70745
## Asymptotic standard error: 0.19541
##     z-value: 0.48504, p-value: 0.62765
## Wald statistic: 0.23526, p-value: 0.62765
## 
## Log likelihood: -101.0725 for error model
## ML residual variance (sigma squared): 3.0765, (sigma: 1.754)
## Number of observations: 51 
## Number of parameters estimated: 13 
## AIC: 228.14, (AIC for lm: 226.29)

#e. SDM

sdm_model <- lagsarlm(log(tasa_auto_) ~ zona_urb + densidad_p + region + edad + I(edad^2) + sexo, data = datab_sp, listw = rswm_rook,  type = "mixed", zero.policy = TRUE)

summary(sdm_model)

## 
## Call:lagsarlm(formula = log(tasa_auto_) ~ zona_urb + densidad_p + 
##     region + edad + I(edad^2) + sexo, data = datab_sp, listw = rswm_rook, 
##     type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -3.79282 -0.64894 -0.15556  0.93272  2.35130 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept)                 47.35986975 20.46237030  2.3145 0.0206411
## zona_urbno_interseccion      3.20903815  0.93503202  3.4320 0.0005991
## zona_urbsuburbana            0.67358354  1.15717247  0.5821 0.5605031
## densidad_p                   0.00148280  0.00054648  2.7134 0.0066602
## regionNorte                 -1.96692344  1.85326644 -1.0613 0.2885409
## regionOeste                 -0.97927960  1.46961306 -0.6664 0.5051861
## regionSur                   -7.68647838  2.22448770 -3.4554 0.0005495
## regionZMM                   -2.42933969  1.55319022 -1.5641 0.1177949
## edad                        -0.11209051  0.24050851 -0.4661 0.6411752
## I(edad^2)                    0.00017619  0.00268659  0.0656 0.9477120
## sexoHombre                  -0.67453373  0.71839422 -0.9389 0.3477582
## lag.zona_urbno_interseccion  0.72292783  1.60827293  0.4495 0.6530669
## lag.zona_urbsuburbana        6.44158201  2.96912888  2.1695 0.0300433
## lag.densidad_p              -0.00206923  0.00117802 -1.7565 0.0789986
## lag.regionNorte              1.55059987  2.21493788  0.7001 0.4838869
## lag.regionOeste             -0.73473253  2.05738052 -0.3571 0.7210017
## lag.regionSur                4.85736288  2.85645036  1.7005 0.0890390
## lag.regionZMM               -0.09947556  2.26350453 -0.0439 0.9649462
## lag.edad                    -1.57014829  0.71969674 -2.1817 0.0291331
## lag.I(edad^2)                0.01571143  0.00806335  1.9485 0.0513554
## lag.sexoHombre               0.31729430  1.55170476  0.2045 0.8379776
## 
## Rho: -0.27812, LR test value: 1.636, p-value: 0.20088
## Asymptotic standard error: 0.19158
##     z-value: -1.4517, p-value: 0.14659
## Wald statistic: 2.1074, p-value: 0.14659
## 
## Log likelihood: -88.93371 for mixed model
## ML residual variance (sigma squared): 1.8851, (sigma: 1.373)
## Number of observations: 51 
## Number of parameters estimated: 23 
## AIC: 223.87, (AIC for lm: 223.5)
## LM test for residual autocorrelation
## test value: 4.2581, p-value: 0.039064

#f. Selección de Modelo de Predicción

model_comparison <- AIC(modelo_ols, sar_model, sem_model, sdm_model)
print(model_comparison)

##            df      AIC
## modelo_ols 12 226.2857
## sar_model  13 228.0810
## sem_model  13 228.1449
## sdm_model  23 223.8674

# Selecciona el modelo con el menor AIC
mejor_modelo <- which.min(model_comparison$AIC)
cat("El mejor modelo es:", rownames(model_comparison)[mejor_modelo], "\n")

## El mejor modelo es: sdm_model

#g. Identificación y mostrar la independencia de εi del modelo de predicción seleccionado.

# Asumiremos que SDM fue el mejor (ajusta según te salga):
residuos <- residuals(sdm_model)

# Test de Moran para verificar autocorrelación de residuos
moran.test(residuos, rswm_rook, zero.policy = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  residuos  
## weights: rswm_rook    
## 
## Moran I statistic standard deviate = -0.70793, p-value = 0.7605
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.085332786      -0.020000000       0.008517018

#h. Visualizar los valores estimados a nivel municipal de la variable dependiente del modelo de predicción seleccionado.

# Agregar los valores estimados al shapefile
datab$pred_sdm <- exp(predict(sdm_model))  # Deslogar para obtener la escala original

## This method assumes the response is known - see manual page

# Mapa de los valores estimados
library(tmap)

tm_shape(datab) +
  tm_fill("pred_sdm", palette = "Blues", title = "Valores Estimados (SDM)") +
  tm_borders() +
  tm_layout(main.title = "Predicción de Accidentes por Municipio (SDM)", 
             legend.outside = TRUE)

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_tm_fill()`: migrate the argument(s) related to the scale of the
## visual variable `fill` namely 'palette' (rename to 'values') to fill.scale =
## tm_scale(<HERE>).
## [v3->v4] `tm_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(main.title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
## "brewer.blues"
## Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.
## 
## [plot mode] fit legend/component: Some legend items or map compoments do not
## fit well, and are therefore rescaled.
## ℹ Set the tmap option `component.autoscale = FALSE` to disable rescaling.

### 4) GWR - Local Regression Analysis

bw1 <- bw.gwr(log(tasa_auto_) ~ densidad_p + edad + gini, 
              data = datab_sp, 
              approach = "AIC", 
              adaptive = TRUE)

## Adaptive bandwidth (number of nearest neighbours): 39 AICc value: 207.7468 
## Adaptive bandwidth (number of nearest neighbours): 32 AICc value: 207.5074 
## Adaptive bandwidth (number of nearest neighbours): 27 AICc value: 209.9128 
## Adaptive bandwidth (number of nearest neighbours): 34 AICc value: 208.0291 
## Adaptive bandwidth (number of nearest neighbours): 29 AICc value: 207.7464 
## Adaptive bandwidth (number of nearest neighbours): 32 AICc value: 207.5074

# a. Estimated GWR Model 
library(GWmodel)
library(shinyjs)

## 
## Attaching package: 'shinyjs'

## The following object is masked from 'package:lmtest':
## 
##     reset

## The following object is masked from 'package:Rcpp':
## 
##     show

## The following object is masked from 'package:sp':
## 
##     show

## The following object is masked from 'package:Matrix':
## 
##     show

## The following object is masked from 'package:VGAM':
## 
##     show

## The following object is masked from 'package:stats4':
## 
##     show

## The following object is masked from 'package:lubridate':
## 
##     show

## The following objects are masked from 'package:methods':
## 
##     removeClass, show

gwr_model <- gwr.basic(formula = log(tasa_auto_) ~ densidad_p + edad + gini,
                       data = datab_sp,
                       bw = bw1,
                       adaptive = TRUE)

summary(gwr_model)

##               Length Class                    Mode
## GW.arguments  11     -none-                   list
## GW.diagnostic  8     -none-                   list
## lm            14     lm                       list
## SDF           51     SpatialPolygonsDataFrame S4  
## timings        5     -none-                   list
## this.call      5     -none-                   call
## Ftests         0     -none-                   list

# b. Display local estimated regression results
gwr_results <- as.data.frame(gwr_model$SDF)
head(gwr_results[, c("densidad_p", "edad", "gini", "Local_R2")])

##     densidad_p         edad      gini  Local_R2
## 1 0.0003296935  0.020962630 11.441243 0.5341382
## 2 0.0003687456  0.026426485 14.211314 0.5793043
## 3 0.0002816742  0.056474649 15.489328 0.6314803
## 4 0.0004979813 -0.001503509 14.429634 0.5816133
## 5 0.0002014885  0.029397564  7.926161 0.2079878
## 6 0.0003339315  0.014632752 10.349688 0.4363378

datab$coef_densidad_p <- gwr_results$densidad_p
datab$coef_edad <- gwr_results$edad
datab$coef_gini <- gwr_results$gini
datab$local_R2 <- gwr_results$LocalR2


# c. Display Local R2

datab$Local_R2 <- gwr_results$Local_R2

tmap_mode("view")

## ℹ tmap mode set to "view".

tm_shape(datab) +
  tm_polygons("Local_R2", 
              palette = "YlGnBu", 
              style = "quantile", 
              title = "R² Local - GWR") +
  tm_layout(title = "GWR Local R² across Municipalities in NL")

## 

## ── tmap v3 code detected ───────────────────────────────────────────────────────

## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
##   'tm_scale_intervals(<HERE>)'
## [v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [v3->v4] `tm_layout()`: use `tm_title()` instead of `tm_layout(title = )`
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "YlGnBu" is named
## "brewer.yl_gn_bu"
## Multiple palettes called "yl_gn_bu" found: "brewer.yl_gn_bu", "matplotlib.yl_gn_bu". The first one, "brewer.yl_gn_bu", is returned.

### Please consider the estimated results from 1) to 4) and briefly describe the main 6 - 8 insights related to the following: 

# a. Descriptive statistics of car accidents across NL municipalities

    # El promedio estatal de accidentes automovilísticos es de 2,389.10, lo que indica una alta incidencia general en algunos municipios.
    # La mediana estatal es de 172 accidentes, lo cual sugiere una distribución sesgada: pocos municipios tienen cifras extremadamente altas, mientras la mayoría presenta valores bajos.
    # Esta diferencia entre media y mediana refleja la presencia de outliers que elevan la media, probablemente concentrados en zonas urbanas

# b. Identification of clusters and spatial autocorrelation of car accidents across NL municipalities

  # El Índice de Moran global indicó la presencia de autocorrelación espacial positiva, lo que significa que los municipios con tasas altas de accidentes tienden a estar cerca de otros municipios con tasas igualmente altas (y viceversa).
  # Este patrón sugiere que los accidentes no están distribuidos al azar, sino que muestran dependencia espacial, posiblemente influenciada por factores regionales compartidos (infraestructura, movilidad, densidad poblacional).

# c. What are the main factors that explain car accidents at the state of Nuevo Leon

    # Zona urbana sin intersecciones: Se asocia con un mayor número de accidentes, lo que sugiere que el diseño vial influye directamente en la siniestralidad.

  # Región Sur: Presenta significativamente menos accidentes que otras regiones, posiblemente por menor densidad poblacional o tráfico.
  # Edad: Tiene un efecto no lineal; los accidentes aumentan con la edad hasta cierto punto y luego disminuyen.
  # Desigualdad (índice de Gini): Aunque no es significativa a nivel global, el modelo GWR muestra que en algunos municipios, mayor desigualdad se relaciona con más accidentes.
  # Densidad poblacional: No es significativa en el modelo global, pero sí lo es localmente en áreas urbanas según GWR.
  # Sexo: Ser hombre se asocia con mayor riesgo de accidentes, aunque no siempre con significancia estadística.

# e. Local Spatial Regression Results

#1. Alto poder explicativo en el oriente y sur del estado
  # Municipios en tonos azul oscuro (R² local entre 0.727 y 0.904) muestran que el modelo explica muy bien la tasa de accidentes en esas zonas.
  # Esto sugiere que en esas regiones las variables como la densidad poblacional, el índice de Gini o la edad media tienen una relación más fuerte y clara con los accidentes.

#2. Bajo poder explicativo en el centro-norte
  # Municipios en tonos crema o verde claro (R² local entre 0.208 y 0.483) indican que el modelo GWR explica menos del 50% de la variabilidad de accidentes en esas áreas.
  # Esto podría deberse a:
  # Presencia de otros factores no incluidos en el modelo (ej. calidad vial, presencia de zonas industriales, variables culturales o de tránsito).
  # Mayor heterogeneidad interna o efectos no espaciales más dominantes.

#3. Zona metropolitana mixta
  # Algunos municipios dentro o cerca del área metropolitana de Monterrey tienen R² medianos (tonos verde azulado), lo que sugiere que el modelo es moderadamente efectivo, pero probablemente podrían mejorar con variables adicionales como nivel de urbanización, tipo de vías o congestión vehicular.

# f. Compare and contrast Global vs Local estimated regression results

#En los modelos globales (OLS, SAR, SEM, SDM), variables como la zona urbana sin intersección, la edad, el cuadrado de la edad, y la región sur resultaron estadísticamente significativas. En particular, el modelo SDM fue el mejor según el criterio AIC, mostrando que hay cierto efecto espacial indirecto a través de variables rezagadas, como la edad y la zona urbana.

#Sin embargo, estos modelos asumen que los efectos son iguales en todo el estado. Por ejemplo, el coeficiente de la densidad poblacional en OLS y SAR fue positivo pero no significativo, lo que podría llevar a pensar que no tiene relación con los accidentes.

#En cambio, el modelo local (GWR) mostró que sí hay variaciones espaciales importantes. Los coeficientes de densidad poblacional, edad y Gini variaron entre municipios. En algunas zonas, la densidad tuvo un efecto positivo claro, y en otras, no. Lo mismo ocurrió con el índice de Gini, que fue especialmente relevante en municipios con mayor desigualdad, algo que no se capturó en los modelos globales.

#Además, el R² local del GWR mostró que en ciertos municipios el modelo explicó hasta un 90% de la variación en accidentes, mientras que en otros apenas alcanzó el 20%. Esto indica que el modelo global subestima o sobreestima los efectos locales, y que la realidad es más diversa de lo que reflejan los promedios.

# 1. ¿Cuáles son los municipios de NL que agrupan el mayor nivel de accidentes de autos por cada 10,000 habitantes?
# Los municipios con mayor nivel de accidentes son principalmente los de la Zona Metropolitana de Monterrey:
# - Monterrey, San Nicolás de los Garza, Guadalupe, Apodaca, Santa Catarina



# 2. ¿Existen clústeres de accidentes de autos a nivel municipal en NL?
# Sí, se identificaron clústeres tanto a nivel global como local, el Índice de Moran Global resultó significativo (p < 0.05), indicando autocorrelación espacial de los accidentes.

# - Con el análisis LISA se detectaron clústeres "High-High" (zonas críticas) y "Low-Low" (zonas de baja incidencia).Los clústeres "High-High" se concentran principalmente en los municipios de la ZMM.



# 3. ¿Cuáles son los factores que inciden en un aumento/disminución de los accidentes de autos en NL?
# Según el modelo espacial SDM, los factores con mayor influencia son: Densidad poblacional,  edad promedio, región: Las zonas urbanas registran más accidentes que las rurales, sexo, indice de Gini: En algunos municipios, mayor desigualdad económica se asocia con más accidentes.



# 4. ¿La locación de dichos factores tiene un impacto significativo sobre el nivel de accidentes de autos en un municipio en particular?
#Sí, se confirmó mediante GWR y el R^2 local varía entre municipios, lo que indica que la influencia de los factores no es homogénea en todo el estado, ademas la densidad poblacional tiene mayor impacto en la ZMM lo que significa que la ubicación geográfica de los factores sí afecta de manera significativa el nivel de accidentes en cada municipio.