Quiz 3

library(foreign)        # import external files
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) 
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

library(shinyjs)

## 
## Attaching package: 'shinyjs'
## 
## 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

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 
df <- read.csv("/Users/hectordelagarzatrevino/Desktop/Quiz_3/nl_mpios_auto_acc.csv")
datab <- st_read("/Users/hectordelagarzatrevino/Desktop/Quiz_3/nl_map/nl_mpios_auto_acc.shp")

## Reading layer `nl_mpios_auto_acc' from data source 
##   `/Users/hectordelagarzatrevino/Desktop/Quiz_3/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/hectordelagarzatrevino/Desktop/Quiz_3/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. 
tm_shape(dataa) + 
  tm_polygons("densidad_p", palette = "Purples", style="quantile", title="NL Population Density across Mpios") +
tmap_mode("view") 

## 

## ── 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>)'
## ℹ tmap mode set to "view".
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Purples" is named
## "brewer.purples"
## Multiple palettes called "purples" found: "brewer.purples", "matplotlib.purples". The first one, "brewer.purples", is returned.

tmap_last()

## ℹ tmap mode set to "view".
## 
## 
## ── 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>)'
## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Purples" is named
## "brewer.purples"
## Multiple palettes called "purples" found: "brewer.purples", "matplotlib.purples". The first one, "brewer.purples", is returned.

1) Exploratory Data Analysis (EDA)

# a. Promedio y mediana de accidentes automovilísticos a nivel estatal y por región del estado.
mean_acc <- mean(dataa$auto_accid, na.rm = TRUE)
median_acc <- median(dataa$auto_accid, na.rm = TRUE)
mean_acc

## [1] 2389.098

median_acc

## [1] 172

# b. Dispersión regional de accidentes automovilísticos.
tm_shape(dataa) +
  tm_polygons("auto_accid", 
              palette = "Reds", 
              style = "quantile", 
              title = "Accidentes por municipio") +
  tm_borders() +
  tm_layout(title = "Dispersión Regional de Accidentes en 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 "Reds" is named
## "brewer.reds"
## Multiple palettes called "reds" found: "brewer.reds", "matplotlib.reds". The first one, "brewer.reds", is returned.

# c. Dispersión de accidentes automovilísticos por uso de cinturón y/o sexo.
boxplot1 <- boxplot(dataa$auto_accid ~ dataa$cinturon, main="Accidentes por uso de cinturón", ylab="Número de accidentes", xlab="Usa cinturón (0=No, 1=Sí)")

boxplot2 <- boxplot(dataa$auto_accid ~ dataa$sexo, main="Accidentes por sexo", ylab="Número de accidentes", xlab="Sexo (1=Mujer, 2=Hombre)")

2) Exploratory Spatial Data Analysis (ESDA)

# a. Calcular y mostrar la matriz de conectividad "Rook".
coords <- st_centroid(st_geometry(dataa))
nb_rook <- poly2nb(dataa, queen = FALSE)
lw_rook <- nb2listw(nb_rook, style = "W")
summary(lw_rook)

## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 51 
## Number of nonzero links: 244 
## Percentage nonzero weights: 9.381007 
## Average number of links: 4.784314 
## Link number distribution:
## 
##  1  2  3  4  5  6  7  8  9 12 
##  2  5  6 10 13  7  3  2  2  1 
## 2 least connected regions:
## 32 37 with 1 link
## 1 most connected region:
## 11 with 12 links
## 
## Weights style: W 
## Weights constants summary:
##    n   nn S0       S1       S2
## W 51 2601 51 24.56053 215.5498

# b. Identificación de clúster global de accidentes automovilísticos.
global_moran <- moran.test(dataa$auto_accid, lw_rook)
global_moran

## 
##  Moran I test under randomisation
## 
## data:  dataa$auto_accid  
## weights: lw_rook    
## 
## Moran I statistic standard deviate = 3.0394, p-value = 0.001185
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.181482128      -0.020000000       0.004394321

# c. Identificación de clúster local de accidentes automovilísticos.
local_moran <- localmoran(dataa$auto_accid, lw_rook)
dataa$local_moran_I <- local_moran[,1]
dataa$p_value <- local_moran[,5]

# d. ¿Cuáles son los municipios que integran el clúster / los clústeres de accidentes de autos en el estado de NL?
dataa$significant_cluster <- ifelse(dataa$p_value < 0.05, "Significativo", "No significativo")
tm_shape(dataa) +
  tm_polygons("significant_cluster", palette = c("gray", "red"), title = "Clúster Local (p < 0.05)") +
  tm_borders() +
  tm_layout(title = "Clústeres de Accidentes en NL")

## 

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

## [v3->v4] `tm_tm_polygons()`: 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_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 = )`
## Multiple palettes called "gray" found: "matplotlib.gray", "tableau.gray". The first one, "matplotlib.gray", is returned.

3) Predictive Modeling using Linear Regression

dataa$region <- as.factor(dataa$region)
dataa$sexo <- as.factor(dataa$sexo)
dataa$zona_urb <- as.factor(dataa$zona_urb)

# a. Modelo de regresión no espacial
ols_model <- lm(auto_accid ~ densidad_p + sexo + cinturon, data = dataa)
summary(ols_model)

## 
## Call:
## lm(formula = auto_accid ~ densidad_p + sexo + cinturon, data = dataa)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
##  -7602  -2165   -135    274  32853 
## 
## Coefficients:
##                       Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          -254.5684  2102.4910  -0.121 0.904155    
## densidad_p              3.0632     0.7735   3.960 0.000258 ***
## sexo2                2598.9795  2233.5647   1.164 0.250586    
## cinturoncinturon_si  -379.8886  2214.6037  -0.172 0.864553    
## cinturonse_ignora   -2312.2521  2102.9440  -1.100 0.277258    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6075 on 46 degrees of freedom
## Multiple R-squared:  0.2838, Adjusted R-squared:  0.2215 
## F-statistic: 4.557 on 4 and 46 DF,  p-value: 0.003487

# c) SAR
sar_model <- lagsarlm(auto_accid ~ densidad_p + sexo + cinturon, data = dataa, listw = lw_rook, zero.policy = TRUE)
summary(sar_model)

## 
## Call:lagsarlm(formula = auto_accid ~ densidad_p + sexo + cinturon, 
##     data = dataa, listw = lw_rook, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -7976.01 -2157.89  -255.72   336.26 32619.97 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                       Estimate Std. Error z value  Pr(>|z|)
## (Intercept)          -150.2085  2009.2996 -0.0748 0.9404083
## densidad_p              3.2315     0.8395  3.8494 0.0001184
## sexo2                2571.6544  2118.2389  1.2141 0.2247275
## cinturoncinturon_si  -345.5266  2100.0339 -0.1645 0.8693109
## cinturonse_ignora   -2236.5753  1997.6769 -1.1196 0.2628893
## 
## Rho: -0.061488, LR test value: 0.11723, p-value: 0.73206
## Asymptotic standard error: 0.196
##     z-value: -0.31372, p-value: 0.75373
## Wald statistic: 0.09842, p-value: 0.75373
## 
## Log likelihood: -513.9853 for lag model
## ML residual variance (sigma squared): 33185000, (sigma: 5760.7)
## Number of observations: 51 
## Number of parameters estimated: 7 
## AIC: 1042, (AIC for lm: 1040.1)
## LM test for residual autocorrelation
## test value: 0.0074248, p-value: 0.93133

# d) SEM
sem_model <- errorsarlm(auto_accid ~ densidad_p + sexo + cinturon, data = dataa, listw = lw_rook, zero.policy = TRUE)
summary(sem_model)

## 
## Call:errorsarlm(formula = auto_accid ~ densidad_p + sexo + cinturon, 
##     data = dataa, listw = lw_rook, zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -7841.33 -2207.77  -245.16   312.71 32681.78 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                        Estimate  Std. Error z value Pr(>|z|)
## (Intercept)          -213.64131  1987.52386 -0.1075   0.9144
## densidad_p              3.09534     0.70881  4.3669 1.26e-05
## sexo2                2563.04856  2128.53841  1.2041   0.2285
## cinturoncinturon_si  -339.30642  2109.61358 -0.1608   0.8722
## cinturonse_ignora   -2287.10359  2007.38938 -1.1393   0.2546
## 
## Lambda: -0.054168, LR test value: 0.078806, p-value: 0.77892
## Asymptotic standard error: 0.20462
##     z-value: -0.26472, p-value: 0.79122
## Wald statistic: 0.070079, p-value: 0.79122
## 
## Log likelihood: -514.0045 for error model
## ML residual variance (sigma squared): 33216000, (sigma: 5763.3)
## Number of observations: 51 
## Number of parameters estimated: 7 
## AIC: 1042, (AIC for lm: 1040.1)

# e) SDM
sdm_model <- lagsarlm(auto_accid ~ densidad_p + sexo + cinturon, data = dataa, listw = lw_rook, type = "mixed", zero.policy = TRUE)
summary(sdm_model)

## 
## Call:lagsarlm(formula = auto_accid ~ densidad_p + sexo + cinturon, 
##     data = dataa, listw = lw_rook, type = "mixed", zero.policy = TRUE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -7643.41 -2019.90  -864.54   429.32 32765.00 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                            Estimate  Std. Error z value Pr(>|z|)
## (Intercept)              3028.12068  5887.70309  0.5143   0.6070
## densidad_p                  3.91255     1.94164  2.0151   0.0439
## sexo2                    2230.81514  2166.27049  1.0298   0.3031
## cinturoncinturon_si      -395.85180  2143.46405 -0.1847   0.8535
## cinturonse_ignora       -2546.88654  2080.99614 -1.2239   0.2210
## lag.densidad_p             -0.95396     2.90230 -0.3287   0.7424
## lag.sexo2               -2774.72939  5142.92709 -0.5395   0.5895
## lag.cinturoncinturon_si   685.24989  5826.16654  0.1176   0.9064
## lag.cinturonse_ignora   -1799.11525  5377.73199 -0.3345   0.7380
## 
## Rho: -0.060908, LR test value: 0.096781, p-value: 0.75573
## Asymptotic standard error: 0.20428
##     z-value: -0.29817, p-value: 0.76558
## Wald statistic: 0.088904, p-value: 0.76558
## 
## Log likelihood: -513.6441 for mixed model
## ML residual variance (sigma squared): 32745000, (sigma: 5722.3)
## Number of observations: 51 
## Number of parameters estimated: 11 
## AIC: 1049.3, (AIC for lm: 1047.4)
## LM test for residual autocorrelation
## test value: 0.62496, p-value: 0.42921

# f) Model Selection
AIC(ols_model, sar_model, sem_model, sdm_model)

##           df      AIC
## ols_model  6 1040.088
## sar_model  7 1041.971
## sem_model  7 1042.009
## sdm_model 11 1049.288

# g) Independencia de errores
lm.morantest(ols_model, lw_rook, zero.policy = TRUE)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = auto_accid ~ densidad_p + sexo + cinturon, data =
## dataa)
## weights: lw_rook
## 
## Moran I statistic standard deviate = 0.0026497, p-value = 0.4989
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     -0.028576666     -0.028821928      0.008567955

# h) Visualización de valores estimados
dataa$predicted_sar <- sar_model$fitted.values  

tm_shape(dataa) +
  tm_polygons("predicted_sar", palette = "Greens", title = "Predicción (SAR)") +
  tm_borders() +
  tm_layout(title = "Predicción del Modelo SAR por Municipio")

## 

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

## [v3->v4] `tm_tm_polygons()`: 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_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 "Greens" is named
## "brewer.greens"
## Multiple palettes called "greens" found: "brewer.greens", "matplotlib.greens". The first one, "brewer.greens", is returned.

4) Estimación de Modelos de Predicción Local

# a. Estimated GWR Model
# Convertir a objeto Spatial
dataa_sp <- as_Spatial(dataa)
dataa_sp

## class       : SpatialPolygonsDataFrame 
## features    : 51 
## extent      : -101.2068, -98.42158, 23.16268, 27.79914  (xmin, xmax, ymin, ymax)
## crs         : +proj=longlat +datum=WGS84 +no_defs 
## variables   : 26
## names       : OBJECTID, CODELAG, CVE_ENT, IDUNICO,     Shape_Leng,       Shape_Area, OBJECTID_1, IDUNICO_1,       mpio, auto_accid, tasa_auto_,     zona_urb,   zona_subur, sexo,   aliento, ... 
## min values  :        1,    2199,      19,   19001, 0.285842845914, 0.00422425985785,          1,     19001,    Abasolo,          1,       0.15, interseccion, camino_rural,    1,  negativo, ... 
## max values  :       51,    2249,      19,   19051,  4.77998092377,   0.630890761882,         51,     19051, Villaldama,      42956,   98483.48,    suburbana,       urbana,    2, se_ignora, ...

# Extraer coordenadas
coords <- coordinates(dataa_sp)

# Estimar el ancho de banda
gwr_bandwidth <- bw.gwr(formula = auto_accid ~ densidad_p + sexo + cinturon,
                        data = dataa_sp,
                        approach = "AICc",
                        kernel = "bisquare",
                        adaptive = TRUE) 

## Adaptive bandwidth (number of nearest neighbours): 39 AICc value: 1059.448 
## Adaptive bandwidth (number of nearest neighbours): 32 AICc value: 1074.38 
## Adaptive bandwidth (number of nearest neighbours): 44 AICc value: 1053.555 
## Adaptive bandwidth (number of nearest neighbours): 46 AICc value: 1051.819 
## Adaptive bandwidth (number of nearest neighbours): 49 AICc value: 1048.274 
## Adaptive bandwidth (number of nearest neighbours): 49 AICc value: 1048.274

# Ajustar el modelo GWR
gwr_model <- gwr.basic(auto_accid ~ densidad_p + sexo + cinturon,
                       data = dataa_sp,
                       bw = gwr_bandwidth,
                       kernel = "bisquare",
                       adaptive = TRUE)

# Agregar resultados a tu shapefile original (sf)
dataa$gwr_pred <- gwr_model$SDF$yhat
dataa$gwr_r2 <- gwr_model$SDF$Local_R2

# 4b. Visualizar los valores estimados (Simplified tmap v4 syntax)
tm_shape(dataa) +
  tm_fill("gwr_pred",
          style = "quantile",
          palette = "Blues",
          title = "Predicción (GWR)") +
  tm_borders() +
  tm_title("Predicción GWR de Accidentes Automovilísticos") +
  tm_layout(legend.position = c("right", "bottom")) +
  tm_compass(type = "8star", position = c("right", "top")) +
  tm_scalebar(position = c("left", "bottom"))

## 

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

## [v3->v4] `tm_fill()`: 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_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
## [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.

# 4c. R2 local (Simplified tmap v4 syntax)
tm_shape(dataa) +
  tm_fill("gwr_r2",
          style = "quantile",
          palette = "Oranges",
          title = "R² Local (GWR)") +
  tm_borders() +
  tm_title("Coeficiente de Determinación Local (GWR)") +
  tm_layout(legend.position = c("right", "bottom")) +
  tm_compass(type = "8star", position = c("right", "top")) +
  tm_scalebar(position = c("left", "bottom"))

## 
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_fill()`: 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_fill()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Oranges" is named
## "brewer.oranges"Multiple palettes called "oranges" found: "brewer.oranges", "matplotlib.oranges". The first one, "brewer.oranges", is returned.

5) Con base en los resultados encontrados en 1) – 4) realizar una breve descripción de los principales 6 – 8 hallazgos identificados.

• Los accidentes automovilísticos no están distribuidos uniformemente en Nuevo León.
• Municipios urbanos (Monterrey, San Pedro, Guadalupe) presentan mayor incidencia que zonas rurales.
• La media estatal (r mean_acc) es significativamente influenciada por los valores atípicos en áreas metropolitanas.
• El test de Moran global (p = r global_moran$p.value) confirma que los accidentes tienden a agruparse. No es un patrón aleatorio: municipios con alta siniestralidad están cercanos entre sí.
• Hotspots (alta incidencia): Área Metropolitana de Monterrey (Monterrey, San Nicolás, Apodaca). Corredores industriales (Pesquería, Ciénega de Flores).
• Coldspots (baja incidencia): Municipios rurales del sur (Aramberri, Galeana).
• Variables significativas en SAR (mejor modelo según AIC): Densidad poblacional (+) es Más población = más accidentes.
• Sexo (Hombre) → Mayor riesgo en conductores masculinos.
• Uso de cinturón (–) → Reduce gravedad/frecuencia.
• Efecto espacial: El coeficiente ρ (r sar_model$rho) indica que los accidentes en un municipio afectan a sus vecinos.
• El impacto de los predictores varía por región:
• Densidad poblacional: Efecto más fuerte en zonas urbanas.
• Cinturón de seguridad: Mayor efecto preventivo en carreteras rurales.
• R² local: Oscila entre r min(dataa\(gwr_r2) y r max(dataa\)gwr_r2), mostrando que el modelo explica mejor ciertas áreas (ej. zonas urbanas).
• OLS/SAR subestiman efectos en regiones específicas.
• GWR nos dice que: Algunas variables cambian de signo (ej. densidad puede tener efecto negativo en áreas periféricas).
• La relación entre variables no es constante en el espacio.
• Enfoque territorial diferenciado: Ciudades, deben tener mejor Control de velocidad y fiscalización de cinturón y en Carreteras rurales se deben crear campañas sobre manejo nocturno y alcohol.
• Para el monitoreo de clústeres críticos: Usar mapas de predicción GWR para priorizar recursos.