Presentación

Este documento analiza el patrón espacial de los atropellos en la Región Metropolitana de Chile (datos de CONASET, 2018–2022) en tres niveles complementarios:

Nivel	Pregunta que responde	Técnica
(A) Descriptivo	¿Dónde se concentran los eventos?	Densidad de kernel (KDE)
(B) Inferencial	¿El agrupamiento es real o azar?	Función K / L de Ripley
(C) Particional	¿Qué grupos discretos existen?	DBSCAN / HDBSCAN

Advertencia de interpretación (importante). Trabajamos con puntos crudos, sin normalizar por exposición (flujo peatonal o vehicular). Por lo tanto, una “zona caliente” indica dónde ocurren más atropellos en términos absolutos, no necesariamente dónde son más probables por persona o por vehículo. Un hot-spot puede reflejar simplemente mayor actividad urbana.

Cómo leer este documento. Cada técnica abre con una idea intuitiva y una guía de lectura del gráfico, y cierra con un ejercicio propuesto. El menú lateral (TOC) funciona como índice; el botón Code permite mostrar u ocultar el código de cada bloque.

Glosario rápido

CRS: sistema de referencia de coordenadas. Geográfico = grados (lon/lat); proyectado (UTM) = metros, apto para medir distancias y áreas.
sf: formato moderno de datos espaciales vectoriales en R.
ppp: point pattern de spatstat (puntos + ventana de estudio).
KDE: Kernel Density Estimation; superficie continua de densidad.
bandwidth (sigma): ancho de banda; controla cuánto se suaviza el KDE.
CSR: Complete Spatial Randomness; hipótesis de puntos al azar.
eps: radio de búsqueda de DBSCAN (en metros, aquí).
minPts: número mínimo de vecinos para considerar una zona “densa”.
ruido: en DBSCAN, puntos que no pertenecen a ningún clúster (= 0).
convex hull: polígono convexo mínimo que encierra un conjunto de puntos.

1 Configuración y paquetes

Ejecuta el siguiente bloque una sola vez antes de compilar (Knit) el documento, para asegurarte de tener todos los paquetes instalados. Está como eval=FALSE para que la compilación no intente instalarlos cada vez.

paquetes <- c("here", "readxl", "sf", "tidyverse", "mapview", "ggspatial",
              "prettymapr", "rosm",   # requeridos por annotation_map_tile()
              "spatstat", "stars", "FNN", "dbscan", "raster")

faltantes <- paquetes[!paquetes %in% rownames(installed.packages())]
if (length(faltantes) > 0) install.packages(faltantes, dependencies = TRUE)

Cargamos raster antes que tidyverse para que dplyr::select() y dplyr::filter() no queden enmascaradas.

library(sf)
library(mapview)
library(ggspatial)
library(spatstat)
library(stars)
library(FNN)
library(raster)
library(tidyverse)   # se carga al final: dplyr gana los conflictos de nombres

Parámetros centralizados (evita “números mágicos” repartidos por el código):

CRS_UTM     <- 32719   # UTM 19S, métrico, adecuado para distancias en la RM
SEED        <- 2025    # semilla para reproducir las simulaciones
N_SIM       <- 39      # simulaciones para las envolventes (más = más lento)
N_SUB_ENV   <- 2000    # submuestra para la envolvente (rapidez)
EPS_DBSCAN  <- 150     # radio (m) para DBSCAN; se justifica con el gráfico kNN
EPS_DBSCAN2 <- 500     # radio (m) del segundo modelo
MINPTS_2    <- 30      # mínimo de puntos del segundo modelo
MINPTS_HDB  <- 40      # mínimo de puntos para HDBSCAN

set.seed(SEED)

2 Lectura de datos

Construimos la ruta con here::here() (buena práctica: trabajar dentro de un Proyecto de RStudio, sin setwd() ni rutas absolutas). Los datos pueden descargarse en https://mapas-conaset.opendata.arcgis.com/.

# Buscamos el .RDS en ubicaciones habituales. Añade/ajusta tu ruta si hace falta.
candidatas <- c(
  here::here("datos", "consolidado_atropellos.RDS"),  # <proyecto>/datos/
  here::here("consolidado_atropellos.RDS"),           # raíz del proyecto
  "consolidado_atropellos.RDS",                        # junto al .Rmd
  file.path(getwd(), "consolidado_atropellos.RDS")
)
ruta_datos <- candidatas[file.exists(candidatas)][1]

if (is.na(ruta_datos)) {
  stop("No encontré 'consolidado_atropellos.RDS'. Rutas probadas:\n",
       paste(" -", candidatas, collapse = "\n"),
       "\n\nColoca el archivo en una de esas carpetas (lo más limpio: ",
       "una carpeta 'datos/' en la raíz del proyecto) o define 'ruta_datos' ",
       "a mano. Carpeta de trabajo actual (getwd): ", getwd())
}
message("Leyendo datos desde: ", ruta_datos)

siniestros_rm <- readRDS(ruta_datos)

class(siniestros_rm)

## [1] "sf"         "data.frame"

st_crs(siniestros_rm)

## Coordinate Reference System:
##   User input: WGS 84 / Pseudo-Mercator 
##   wkt:
## PROJCRS["WGS 84 / Pseudo-Mercator",
##     BASEGEOGCRS["WGS 84",
##         ENSEMBLE["World Geodetic System 1984 ensemble",
##             MEMBER["World Geodetic System 1984 (Transit)"],
##             MEMBER["World Geodetic System 1984 (G730)"],
##             MEMBER["World Geodetic System 1984 (G873)"],
##             MEMBER["World Geodetic System 1984 (G1150)"],
##             MEMBER["World Geodetic System 1984 (G1674)"],
##             MEMBER["World Geodetic System 1984 (G1762)"],
##             MEMBER["World Geodetic System 1984 (G2139)"],
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ENSEMBLEACCURACY[2.0]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4326]],
##     CONVERSION["Popular Visualisation Pseudo-Mercator",
##         METHOD["Popular Visualisation Pseudo Mercator",
##             ID["EPSG",1024]],
##         PARAMETER["Latitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",0,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Web mapping and visualisation."],
##         AREA["World between 85.06°S and 85.06°N."],
##         BBOX[-85.06,-180,85.06,180]],
##     ID["EPSG",3857]]

head(siniestros_rm)

3 Exploración visual inicial

Primero observamos los datos en bruto sobre un mapa base de OpenStreetMap.

ggplot(data = siniestros_rm) +
  annotation_map_tile(type = "osm", zoom = 12) +
  geom_sf(inherit.aes = FALSE, color = "darkblue", size = 0.5, alpha = 0.6) +
  labs(
    title    = "Atropellos en la Región Metropolitana (2018–2022)",
    subtitle = "Datos crudos de localización (sin normalizar por exposición)",
    caption  = "Fuente: CONASET"
  ) +
  theme_minimal()

Versión interactiva (se incrusta como mapa navegable en el HTML):

mapview(siniestros_rm, cex = 1, col.regions = "red", layer.name = "Atropellos")

4 Limpieza de datos

na.omit() sobre un objeto sf borra filas con NA en cualquier columna, lo que puede eliminar más de lo deseado. Por eso primero descartamos las geometrías vacías y luego, si hace falta, quitamos NA.

siniestros_rm_clean <- siniestros_rm[!st_is_empty(siniestros_rm), ]
siniestros_rm_clean <- na.omit(siniestros_rm_clean)

cat("Registros originales:", nrow(siniestros_rm),
    "| tras limpieza:", nrow(siniestros_rm_clean), "\n")

## Registros originales: 10683 | tras limpieza: 3710

5 Proyección a UTM y ventana de observación

Proyectamos a UTM para medir densidades y distancias en metros (no en grados). La ventana es la caja envolvente de los puntos: el área de estudio para spatstat. Calculamos las coordenadas una sola vez y las reutilizamos.

siniestros_utm <- st_transform(siniestros_rm_clean, CRS_UTM)

win <- siniestros_utm |>
  st_bbox() |>
  st_as_sfc()

coords <- st_coordinates(siniestros_utm)
x <- coords[, 1]
y <- coords[, 2]

6 (A) Densidad de kernel (KDE)

¿Qué pregunta responde? ¿Dónde se concentran los atropellos? Convierte una nube de puntos discreta en una superficie continua de intensidad.

Idea intuitiva. Imagina que sobre cada punto colocas una pequeña “montañita” de probabilidad (el kernel). Donde hay muchos puntos cercanos las montañitas se suman y forman picos; donde hay pocos, el terreno queda plano. El parámetro clave es el ancho de banda (sigma): define cuán anchas son esas montañitas. Sigma pequeño da un mapa detallado pero ruidoso; sigma grande da un mapa suave que muestra tendencias generales. Es un método no paramétrico: no asume una forma matemática previa para la distribución.

Cómo leer el gráfico. Colores más intensos = mayor densidad de eventos.

box <- as.owin(win)
pp1 <- ppp(x = x, y = y, window = box)
plot(pp1, main = "Patrón de puntos (atropellos)")

kernel <- density(pp1, sigma = bw.diggle)
plot(kernel, main = "Densidad de Kernel — ancho de banda automático")

6.1 Sensibilidad al ancho de banda

El resultado del KDE depende del ancho de banda. Lo comparamos lado a lado para que la elección no parezca neutral: a la izquierda, el sigma automático (bw.diggle); a la derecha, el doble.

sigma_auto   <- bw.diggle(pp1)
sigma_grande <- 2 * sigma_auto

kernel_fino  <- density(pp1, sigma = sigma_auto)
kernel_suave <- density(pp1, sigma = sigma_grande)

op <- par(mfrow = c(1, 2))
plot(kernel_fino,  main = paste0("KDE fino (sigma=", round(sigma_auto), " m)"))
plot(kernel_suave, main = paste0("KDE suave (sigma=", round(sigma_grande), " m)"))

par(op)

El KDE también puede verse como un ráster en un mapa interactivo:

kernel_stars <- kernel |>
  st_as_stars() |>
  st_set_crs(CRS_UTM)

mapview(kernel_stars, layer.name = "Densidad (KDE)")

Ejercicio 1 (densidad de kernel).

Cambia sigma_grande a 0.5 * sigma_auto y vuelve a compilar. ¿Aparecen más o menos hot-spots? ¿Algunos parecen ruido?

Prueba sigma = 4 * sigma_auto. ¿Qué información pierdes al suavizar tanto?

Si tuvieras que recomendar un solo mapa a una autoridad de tránsito, ¿qué ancho de banda elegirías y por qué? Recuerda la advertencia de exposición de la presentación.

7 (B) Función K / L de Ripley con envolventes

¿Qué pregunta responde? El agrupamiento que vemos, ¿es real o podría deberse al azar? La K de Ripley compara la distribución observada con la esperada bajo aleatoriedad completa (CSR). Usamos la función L (la K linealizada, más fácil de leer) y envolventes de simulación para tener un criterio estadístico.

Cómo leer el gráfico.

La banda sombreada = rango esperado bajo azar (intervalo de simulación).
Curva observada por encima de la banda → agrupamiento significativo.
Curva dentro de la banda → compatible con el azar.
Curva por debajo → dispersión / regularidad.

k_ripley <- Kest(pp1)
plot(k_ripley, main = "Función K de Ripley (observada vs. teórica)")

Nota de rendimiento. Las envolventes simulan N_SIM patrones; con decenas de miles de puntos esto es muy lento. Por eso usamos una submuestra aleatoria (N_SUB_ENV) solo para esta prueba. La submuestra es reproducible gracias a set.seed(). El chunk usa cache=TRUE para no recalcular en cada compilación.

n_sub  <- min(N_SUB_ENV, pp1$n)
idx    <- sample.int(pp1$n, n_sub)
pp_sub <- pp1[idx]

env_L <- envelope(pp_sub, fun = Lest, nsim = N_SIM, verbose = FALSE)
plot(env_L, main = paste0("Envolvente de la función L (n=", n_sub, ", ",
                          N_SIM, " sim.)"))

Ejercicio 2 (función L de Ripley).

Sube N_SIM a 99 y N_SUB_ENV a 4000, y vuelve a compilar. ¿Se mantiene la conclusión de agrupamiento? ¿La banda se angosta?

¿Por qué usamos una submuestra en vez de todos los puntos? ¿Qué riesgo tiene una submuestra demasiado pequeña?

El KDE ya “mostraba” agrupamiento. ¿Qué añade la prueba de Ripley que el mapa de calor por sí solo no permite afirmar?

7.1 Alternativa descriptiva con ggplot: `stat_density2d`

Misma idea que el KDE, pero dentro de la gramática de ggplot2. Usamos after_stat(level) (sintaxis moderna; ..level.. está obsoleto).

ggplot() +
  geom_sf(data = siniestros_utm, size = 0.001, alpha = 0.3) +
  stat_density2d(
    aes(x = x, y = y,
        fill  = after_stat(level),
        alpha = after_stat(level)),
    geom = "polygon"
  ) +
  scale_fill_gradient(name = "Nivel de densidad",
                      low = "lightblue", high = "darkblue",
                      guide = guide_colorbar()) +
  scale_alpha_continuous(name = "Nivel de densidad", guide = guide_colorbar()) +
  labs(title = "Densidad de atropellos (stat_density2d)") +
  theme_minimal()

8 (C) Clustering por densidad: DBSCAN

¿Qué pregunta responde? ¿Existen grupos discretos de eventos y cuáles puntos son “ruido”? DBSCAN define clústeres por densidad local y no exige fijar el número de grupos de antemano.

Idea intuitiva. DBSCAN recorre los puntos preguntándose: “¿tengo al menos minPts vecinos dentro de un radio eps?”. Si la respuesta es sí, ese punto es un núcleo y arrastra a sus vecinos al mismo clúster, que crece como una mancha de aceite mientras siga habiendo densidad suficiente. Los puntos aislados se etiquetan como ruido (cluster = 0). Frente a k-means, sus ventajas son: no hay que decir cuántos grupos queremos, detecta formas irregulares e identifica el ruido. Su desventaja: depende de eps y minPts.

8.1 Elección de los parámetros

Para minPts seguimos la regla de Ester et al. (1996): minPts = log(n).

min_pts <- floor(log(nrow(coords)))
min_pts

## [1] 8

Para eps usamos el gráfico de distancia al k-ésimo vecino más cercano. El “codo” de la curva marca la transición entre zonas densas y ruido. La línea roja muestra el eps elegido, para que se vea de dónde sale.

plot_kNN_codo <- function(coords, k, eps_elegido) {
  knn   <- FNN::get.knn(coords, k = k)
  dists <- apply(knn$nn.dist, 1, max)
  plot(sort(dists), type = "l",
       main = paste("Distancia al", k, "vecino más cercano"),
       ylab = "Distancia (m)", xlab = "Punto (ordenado)")
  abline(h = eps_elegido, col = "red", lty = 2, lwd = 2)
  legend("topleft", legend = paste("eps elegido =", eps_elegido, "m"),
         col = "red", lty = 2, lwd = 2, bty = "n")
}

plot_kNN_codo(coords, k = min_pts, eps_elegido = EPS_DBSCAN)

8.2 Modelo y resultados

modelo <- dbscan::dbscan(coords, eps = EPS_DBSCAN, minPts = min_pts)
modelo

## DBSCAN clustering for 3710 objects.
## Parameters: eps = 150, minPts = 8
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 59 cluster(s) and 3106 noise points.
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
## 3106   12    8   15   20    9   16    8    9   10    8   13    9    8    8   11 
##   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31 
##    8   20   12   15   11   10   12   12    8   19    8   10    8    8    8    8 
##   32   33   34   35   36   37   38   39   40   41   42   43   44   45   46   47 
##    8    8   12    8    9   11   11    8    8   15    9    8   14   10   10    9 
##   48   49   50   51   52   53   54   55   56   57   58   59 
##   12    8    8    8   11    8    8    8    8    9    8    9 
## 
## Available fields: cluster, eps, minPts, metric, borderPoints

siniestros_utm$cluster <- as.factor(modelo$cluster)  # 0 = ruido

ggplot(siniestros_utm) +
  geom_sf(aes(color = cluster), size = 1, show.legend = FALSE) +
  scale_color_viridis_d(option = "plasma", end = 0.85) +
  theme_minimal() +
  labs(title = "Clústeres detectados por DBSCAN",
       subtitle = paste0("eps = ", EPS_DBSCAN, " m, minPts = ", min_pts))

Separamos los clústeres del ruido:

siniestros_clustered <- siniestros_utm |>
  dplyr::filter(cluster != 0)

ggplot(siniestros_clustered) +
  geom_sf(aes(color = cluster), size = 1, show.legend = FALSE) +
  scale_color_viridis_d(option = "plasma", end = 0.85) +
  theme_minimal() +
  labs(title = "Clústeres DBSCAN (sin puntos de ruido)")

siniestros_ruido <- siniestros_utm |>
  dplyr::filter(cluster == 0)

ggplot(siniestros_ruido) +
  geom_sf(color = "grey40", size = 0.5) +
  theme_minimal() +
  labs(title = "DBSCAN: puntos clasificados como ruido")

8.3 Envolvente convexa de cada clúster

El convex hull es el polígono convexo más pequeño que encierra todos los puntos de un clúster. Delimita su extensión espacial y permite comparar el “territorio” (área) entre grupos. Como estamos en UTM (metros), el área sale en m². Pero a la escala de estos clústeres (eps = 150 m), expresarla en km² daría valores diminutos (0,0…); por eso la reportamos en hectáreas (1 ha = 10.000 m²) y añadimos el radio equivalente en metros para tener intuición de tamaño, comparable con el propio eps.

Antes de calcular, un breve diagnóstico: ¿cuántos clústeres reales hay y cómo se llama la columna de geometría? (Si todo quedó como ruido, no habrá envolventes; si la columna no se llama geometry, conviene no nombrarla.)

table(siniestros_utm$cluster)                 # 0 = ruido

## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12   13   14   15 
## 3106   12    8   15   20    9   16    8    9   10    8   13    9    8    8   11 
##   16   17   18   19   20   21   22   23   24   25   26   27   28   29   30   31 
##    8   20   12   15   11   10   12   12    8   19    8   10    8    8    8    8 
##   32   33   34   35   36   37   38   39   40   41   42   43   44   45   46   47 
##    8    8   12    8    9   11   11    8    8   15    9    8   14   10   10    9 
##   48   49   50   51   52   53   54   55   56   57   58   59 
##   12    8    8    8   11    8    8    8    8    9    8    9

nrow(siniestros_clustered)                    # puntos en algún clúster

## [1] 604

attr(siniestros_clustered, "sf_column")       # nombre de la geometría activa

## [1] "geometry"

Calculamos el hull dejando que summarise() una la geometría activa sin nombrarla: así funciona sin importar cómo se llame la columna. Luego nos quedamos solo con los resultados que son polígonos (un clúster casi colineal puede dar una línea o un punto, con área 0).

hulls <- siniestros_clustered |>
  dplyr::group_by(cluster) |>
  dplyr::summarise(.groups = "drop") |>        # une la geometría activa de cada clúster
  st_convex_hull()

hulls <- hulls[st_geometry_type(hulls) %in% c("POLYGON", "MULTIPOLYGON"), ]

st_geometry_type(hulls)                        # verificación: deberían ser POLYGON

##  [1] POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON
## [10] POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON
## [19] POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON
## [28] POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON
## [37] POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON
## [46] POLYGON POLYGON POLYGON POLYGON POLYGON POLYGON
## 18 Levels: GEOMETRY POINT LINESTRING POLYGON MULTIPOINT ... TRIANGLE

nrow(hulls)                                    # nº de polígonos a dibujar

## [1] 51

# Área en HECTÁREAS (1 ha = 10.000 m²): unidad adecuada a la escala urbana.
# A eps = 150 m, los clústeres son de barrio/intersección; en km² casi todo
# saldría como 0,0... con muchos ceros, difícil de comparar.
area_m2 <- as.numeric(st_area(hulls))
hulls$area_ha       <- area_m2 / 1e4
# Radio equivalente (m): radio de un círculo de igual área. Da intuición de
# tamaño y es comparable con el eps = 150 m de DBSCAN.
hulls$radio_equiv_m <- sqrt(area_m2 / pi)

hulls |>
  sf::st_drop_geometry() |>
  dplyr::select(cluster, area_ha, radio_equiv_m) |>
  dplyr::arrange(dplyr::desc(area_ha))

# Zoom: recortamos a la caja envolvente de los clústeres + un margen, para que
# los polígonos no se vean diminutos frente al extent completo de la RM.
margen_m <- 1000
bb <- st_bbox(siniestros_clustered)

ggplot() +
  geom_sf(data = siniestros_clustered, aes(color = cluster),
          size = 0.6, alpha = 0.5, show.legend = FALSE) +
  geom_sf(data = hulls, aes(fill = cluster),
          alpha = 0.30, color = "black", linewidth = 0.5, show.legend = FALSE) +
  scale_fill_viridis_d(option = "plasma", end = 0.85) +
  scale_color_viridis_d(option = "plasma", end = 0.85) +
  coord_sf(xlim = c(bb["xmin"] - margen_m, bb["xmax"] + margen_m),
           ylim = c(bb["ymin"] - margen_m, bb["ymax"] + margen_m),
           expand = FALSE) +
  theme_minimal() +
  labs(title = "Envolventes convexas de los clústeres DBSCAN",
       subtitle = "Cada polígono delimita la extensión espacial de un clúster")

Si los clústeres están repartidos por toda la región, el recorte anterior ayuda poco. Para inspeccionar un clúster en detalle, esta función hace zoom a uno solo (cambia el id por el clúster que quieras ver):

zoom_cluster <- function(id, margen = 300) {
  pts  <- dplyr::filter(siniestros_clustered, cluster == id)
  hull <- dplyr::filter(hulls, cluster == id)
  bb   <- st_bbox(hull)
  ggplot() +
    geom_sf(data = pts,  color = "#3b0f70", size = 1.2) +
    geom_sf(data = hull, fill = "#fca50a", alpha = 0.30,
            color = "black", linewidth = 0.7) +
    coord_sf(xlim = c(bb["xmin"] - margen, bb["xmax"] + margen),
             ylim = c(bb["ymin"] - margen, bb["ymax"] + margen),
             expand = FALSE) +
    theme_minimal() +
    labs(title = paste("Detalle del clúster", id))
}

zoom_cluster(1)   # cambia el número por el clúster que quieras examinar

Mapa interactivo con puntos y envolventes:

mapview(siniestros_clustered, zcol = "cluster", cex = 1,
        legend = FALSE, layer.name = "Clústeres DBSCAN") +
  mapview(hulls, zcol = "cluster", alpha.regions = 0.2,
          legend = FALSE, layer.name = "Envolventes convexas")

Nota. Un clúster con menos de 3 puntos no colineales no forma polígono (su hull es una línea o un punto, con área 0); con minPts esto es poco frecuente.

Ejercicio 3 (DBSCAN y envolventes).

Cambia EPS_DBSCAN a 100 y luego a 300, manteniendo minPts. Anota en una tabla el número de clústeres y de puntos de ruido en cada caso.

Mira el gráfico kNN: ¿el eps elegido cae cerca del “codo” de la curva? Justifica tu elección.

¿Cuál es el clúster de mayor área? ¿Coincide con el hot-spot más intenso del KDE? Comenta las diferencias.

8.4 DBSCAN con otros parámetros

Radio mayor (500 m) y más puntos mínimos (30): clústeres más grandes y conservadores. Comparar con el modelo anterior muestra la sensibilidad a los parámetros.

modelo_1 <- dbscan::dbscan(coords, eps = EPS_DBSCAN2, minPts = MINPTS_2)
modelo_1

## DBSCAN clustering for 3710 objects.
## Parameters: eps = 500, minPts = 30
## Using euclidean distances and borderpoints = TRUE
## The clustering contains 12 cluster(s) and 3149 noise points.
## 
##    0    1    2    3    4    5    6    7    8    9   10   11   12 
## 3149   47   36   48   95   47   34   49   39   55   33   44   34 
## 
## Available fields: cluster, eps, minPts, metric, borderPoints

siniestros_utm$cluster_1 <- as.factor(modelo_1$cluster)

ggplot(siniestros_utm) +
  geom_sf(aes(color = cluster_1), size = 1, show.legend = FALSE) +
  scale_color_viridis_d(option = "plasma", end = 0.85) +
  theme_minimal() +
  labs(title = "Clústeres DBSCAN — parámetros más conservadores",
       subtitle = paste0("eps = ", EPS_DBSCAN2, " m, minPts = ", MINPTS_2))

9 Modelo jerárquico: HDBSCAN

¿Qué aporta? HDBSCAN no necesita eps, detecta clústeres de densidad variable, elimina el ruido automáticamente y extrae los grupos más estables de una jerarquía de densidades.

modelo_hdb <- dbscan::hdbscan(coords, minPts = MINPTS_HDB)
table(modelo_hdb$cluster)   # cluster 0 = ruido

## 
##    0    1    2    3    4    5    6    7    8    9   10   11 
## 1956   47  173   62   60   62  127   95   70   76  933   49

siniestros_utm$cluster_hdb <- as.factor(modelo_hdb$cluster)

ggplot(siniestros_utm) +
  geom_sf(aes(color = cluster_hdb == "0"), size = 1) +
  scale_color_manual(values = c("FALSE" = "#3b0f70", "TRUE" = "grey70"),
                     name = NULL, labels = c("En clúster", "Ruido")) +
  theme_minimal() +
  labs(title = "Clústeres detectados por HDBSCAN",
       subtitle = paste0("minPts = ", MINPTS_HDB))

Ejercicio 4 (HDBSCAN y comparación de métodos).

Compara el número de clústeres de DBSCAN, DBSCAN conservador y HDBSCAN. ¿Cuál detecta más grupos? ¿Cuál clasifica más ruido?

HDBSCAN no necesita eps. ¿Qué ventaja práctica tiene esto cuando los clústeres tienen densidades muy distintas (centro denso vs. periferia)?

Para un informe a una autoridad de tránsito, ¿qué combinación de métodos (KDE + Ripley + clustering) presentarías y en qué orden?

10 Síntesis para discusión

Descriptivo (KDE): muestra dónde hay más eventos, pero depende del ancho de banda y no dice si el patrón es distinto del azar.
Inferencial (Ripley): permite afirmar si hay agrupamiento significativo, pero no localiza los grupos.
Particional (DBSCAN/HDBSCAN): identifica y delimita grupos discretos y el ruido; es sensible a sus parámetros.

Los tres son complementarios: juntos responden “¿dónde?”, “¿es real?” y “¿qué grupos?”. Y, como trabajamos con datos crudos, todo se interpreta en términos de conteo absoluto, no de riesgo por exposición.

11 Guardar resultados y registro de sesión

saveRDS(siniestros_utm, here::here("datos", "siniestros_utm_clusters.RDS"))

sessionInfo()

## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Chile.utf8  LC_CTYPE=Spanish_Chile.utf8   
## [3] LC_MONETARY=Spanish_Chile.utf8 LC_NUMERIC=C                  
## [5] LC_TIME=Spanish_Chile.utf8    
## 
## time zone: America/Santiago
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] lubridate_1.9.5        forcats_1.0.1          stringr_1.6.0         
##  [4] dplyr_1.1.4            purrr_1.2.2            readr_2.2.0           
##  [7] tidyr_1.3.2            tibble_3.3.0           ggplot2_4.0.3         
## [10] tidyverse_2.0.0        raster_3.6-32          sp_2.2-0              
## [13] FNN_1.1.4.1            stars_0.7-2            abind_1.4-8           
## [16] spatstat_3.6-0         spatstat.linnet_3.5-1  spatstat.model_3.7-1  
## [19] rpart_4.1.27           spatstat.explore_3.8-1 nlme_3.1-169          
## [22] spatstat.random_3.5-0  spatstat.geom_3.8-1    spatstat.univar_3.2-0 
## [25] spatstat.data_3.1-9    ggspatial_1.1.10       mapview_2.11.4        
## [28] sf_1.1-1              
## 
## loaded via a namespace (and not attached):
##  [1] DBI_1.3.0               deldir_2.0-4            leafpop_0.1.0          
##  [4] rlang_1.2.0             magrittr_2.0.3          otel_0.2.0             
##  [7] e1071_1.7-16            compiler_4.5.1          mgcv_1.9-4             
## [10] png_0.1-8               systemfonts_1.3.2       vctrs_0.6.5            
## [13] pkgconfig_2.0.3         fastmap_1.2.0           labeling_0.4.3         
## [16] leafem_0.2.5            rmarkdown_2.31          tzdb_0.5.0             
## [19] xfun_0.57               satellite_1.0.6         cachem_1.1.0           
## [22] jsonlite_2.0.0          goftest_1.2-3           rosm_0.3.1             
## [25] uuid_1.2-2              spatstat.utils_3.2-3    terra_1.8-60           
## [28] parallel_4.5.1          R6_2.6.1                bslib_0.11.0           
## [31] stringi_1.8.7           RColorBrewer_1.1-3      jquerylib_0.1.4        
## [34] Rcpp_1.1.0              knitr_1.51              tensor_1.5.1           
## [37] base64enc_0.1-3         leaflet.providers_3.0.0 Matrix_1.7-5           
## [40] splines_4.5.1           timechange_0.4.0        tidyselect_1.2.1       
## [43] rstudioapi_0.18.0       yaml_2.3.12             codetools_0.2-20       
## [46] lattice_0.22-9          plyr_1.8.9              withr_3.0.2            
## [49] S7_0.2.0                evaluate_1.0.5          isoband_0.3.0          
## [52] units_0.8-7             proxy_0.4-27            polyclip_1.10-7        
## [55] pillar_1.11.1           KernSmooth_2.23-26      stats4_4.5.1           
## [58] dbscan_1.2.4            generics_0.1.4          rprojroot_2.1.1        
## [61] hms_1.1.4               scales_1.4.0            class_7.3-23           
## [64] glue_1.8.0              tools_4.5.1             grid_4.5.1             
## [67] crosstalk_1.2.2         cli_3.6.5               brew_1.0-10            
## [70] spatstat.sparse_3.2-0   textshaping_1.0.5       prettymapr_0.2.5       
## [73] viridisLite_0.4.3       svglite_2.2.2           gtable_0.3.6           
## [76] sass_0.4.10             digest_0.6.37           classInt_0.4-11        
## [79] htmlwidgets_1.6.4       farver_2.1.2            htmltools_0.5.8.1      
## [82] lifecycle_1.0.5         leaflet_2.2.3           here_1.0.2             
## [85] MASS_7.3-65

Análisis espacial de puntos: densidad, autocorrelación y clustering

R para el análisis espacial — Clase N° 4

Elaborado por Dr. Cristian Escobedo Catalán

2026-06-03

Presentación

Glosario rápido

1 Configuración y paquetes

2 Lectura de datos

3 Exploración visual inicial

4 Limpieza de datos

5 Proyección a UTM y ventana de observación

6 (A) Densidad de kernel (KDE)

6.1 Sensibilidad al ancho de banda

7 (B) Función K / L de Ripley con envolventes

7.1 Alternativa descriptiva con ggplot: `stat_density2d`

8 (C) Clustering por densidad: DBSCAN

8.1 Elección de los parámetros

8.2 Modelo y resultados

8.3 Envolvente convexa de cada clúster

8.4 DBSCAN con otros parámetros

9 Modelo jerárquico: HDBSCAN

10 Síntesis para discusión

11 Guardar resultados y registro de sesión

Análisis espacial de puntos: densidad, autocorrelación y clustering

R para el análisis espacial — Clase N° 4

Elaborado por Dr. Cristian Escobedo Catalán

2026-06-03

Presentación

Glosario rápido

1 Configuración y paquetes

2 Lectura de datos

3 Exploración visual inicial

4 Limpieza de datos

5 Proyección a UTM y ventana de observación

6 (A) Densidad de kernel (KDE)

6.1 Sensibilidad al ancho de banda

7 (B) Función K / L de Ripley con envolventes

7.1 Alternativa descriptiva con ggplot: stat_density2d

8 (C) Clustering por densidad: DBSCAN

8.1 Elección de los parámetros

8.2 Modelo y resultados

8.3 Envolvente convexa de cada clúster

8.4 DBSCAN con otros parámetros

9 Modelo jerárquico: HDBSCAN

10 Síntesis para discusión

11 Guardar resultados y registro de sesión

7.1 Alternativa descriptiva con ggplot: `stat_density2d`