PEC1-raullendines-MD
Raul Lendines Ramos
2025-10-18
Ejercicio 1
- Comprensión del negocio - ¿Qué necesita el negocio?
Imaginemos un juego de datos de una empresa que se dedica a gestionar comedores escolares. Dicha empresa tiene unas tarifas fijas impuestas por la Generalitat de Catalunya para las escuelas públicas. El negocio necesita maximizar su margen de beneficio, controlando los gastos sin empeorar la calidad del servicio y minimizar la merma de comida. Para conseguir estos objetivos, debemos predecir el número exacto de comensales diarios por escuela con una tareas de Predicción/Regresión y una tarea de Clustering para identificar patrones de gasto ineficiente.
- Comprensión de los datos - ¿Qué datos tenemos/necesitamos? ¿Están limpios?
Para los suministros sabemos el costo diario y semanal de los alimentos por categoría, el proveedor, precio de compra o cantidad de merma por escuela.
Para la demanda sabemos el número de comensales diarios en cada escuela ya que los alumnos pagan al hacer la matricula en el colegio, sabemos el tipo de menú que cada alumno comerá y deberíamos predecir el número de ausencias/bajas en base a años anteriores.
Como datos de operación, sabemos el número de empleados y las horas trabajadas por escuela. Para la limpieza, debemos verificar inconsistencias de costo por kilo, valores nulos en mermas y errores en la cuenta de comensales.
En esta fase, deberíamos identificar picos de gasto que no correspondan con el número de comensales para encontrar posibles ineficiencias de compra.
- Preparación de los datos - ¿Cómo organizamos los datos para el modelado?
En esta fase debemos limpiar los datos, añadiendo datos faltantes, eliminando redundancias, arreglando errores, etc. Tendremos que crear variables de tiempo para “dia de la semana”, “proximo a festivo”, “menú específico del dia”, tambien tendremos que crear una variable que contenga el numero total de comensales servidos cada dia para saber la demanda de este dia de la semana dias, semanas u meses anteriores.
Una vez limpiados los datos y añadidias las nuevas variables se dividirian los datos en conjuntos para poder entrenar al modelo y probar su validez como modelo de predicción.
- Modelado - ¿Qué técnicas de modelado debemos aplicar?
En esta fase queremos construir modelos que resuelvan el problema de la demanda de alimentos.
Para el objetivo principal de predicción/regresión deberiamos aplicar un modelo predictivo ya que buscamos predecir un valor numerico del numero real de comensales. Para el objetivo secundario, la tarea de encontrar similitudes y agrupar escuelas para generar patrones, el modelo de agregación agruparía escuelas con patrones similares de costo, merma o ineficiencia.
Por último, para identificar tendencias que clasifiquen escuelas por nivel de merma podemos utilizar arboles de decisión.
- Evaluación - ¿Qué modelo cumple mejor con los objetivos del negocio?
Una vez creados los modelos, podemos comenzar a evaluarlos con los objetivos del negocio y determinar su calidad y validez. Tambien tendremos que comparar el ahorro potencial del modelo de predicción contra el costo de sobreestimar o subestimar suministros para evitar merma.
- Implementación - ¿Cómo acceden los interesados a los resultados?
El modelo debe integrarse en los procesos de la empresa y requeriría de una transformación a un entorno de producción al que se pueda acceder, como por ejemplo un Dashboard al que logisitica pueda acceder. Debemos tener en cuenta que los modelos envejecen y hay que monitorizarlo para que no sea menos preciso.
Ejercicio 2
Descripción del origen del conjunto de datos
El conjunto de datos “TABLA_ACCIDENTES_22.xlsx” pertenece a la Dirección General de Tráfico, un organismo autónomo del Gobierno de España dependiente del Ministerio del Interior. El fichero “microdatos de accidentes con víctimas 2022” contiene información de los accidentes con víctimas en vias urbanas e interurbanas en el estado en 2022. Cada línea del fichero corresponde a un accidente con víctimas. El objetivo analítico es el de predecir la gravedad de los accidentes en las vias urbanas e interurbanas de Cataluña para determinar los factores de mayor riesgo.
Análisis exploratorio
Queremos hacer una aproximación al conjunto de datos escogido y responder a las preguntas más básicas: ¿Cuánto registros tiene? ¿Cuántas variables? ¿De qué tipología son? ¿Cómo se distribuyen los valores de las variables? ¿Hay problemas con los datos, por ejemplo, campos vacíos? ¿Puedo intuir ya el valor analítico de los datos? ¿Qué primeras conclusiones puedo extraer?
Lo primero que haremos será cargar el CSV en nuestro entorno R.
library(readxl)
accidentes <- read_excel("~/Downloads/TABLA_ACCIDENTES_22.xlsx")Exploración del conjunto de datos
Verificamos si la estructura del juego de datos principal está correcta.
structure = str(accidentes)## tibble [97,916 × 73] (S3: tbl_df/tbl/data.frame)
## $ ID_ACCIDENTE : num [1:97916] 1 2 3 4 5 6 7 8 9 10 ...
## $ ANYO : num [1:97916] 2022 2022 2022 2022 2022 ...
## $ MES : num [1:97916] 1 1 1 1 1 1 1 1 1 1 ...
## $ DIA_SEMANA : num [1:97916] 6 5 6 7 7 5 4 5 6 3 ...
## $ HORA : num [1:97916] 3 16 21 12 12 0 18 7 18 17 ...
## $ COD_PROVINCIA : num [1:97916] 1 1 1 1 1 1 1 1 1 1 ...
## $ COD_MUNICIPIO : chr [1:97916] "00000" "01059" "00000" "00000" ...
## $ ISLA : logi [1:97916] NA NA NA NA NA NA ...
## $ ZONA : num [1:97916] 1 1 1 1 1 1 1 1 1 1 ...
## $ ZONA_AGRUPADA : num [1:97916] 1 1 1 1 1 1 1 1 1 1 ...
## $ CARRETERA : chr [1:97916] "A-1" "A-1" "N-104" "BI-636" ...
## $ KM : num [1:97916] 326 355 362 14 21 24 26 30 9 50 ...
## $ SENTIDO_1F : num [1:97916] 2 1 1 2 3 1 2 1 1 2 ...
## $ TITULARIDAD_VIA : num [1:97916] 3 3 5 3 5 3 2 2 5 3 ...
## $ TIPO_VIA : num [1:97916] 3 8 8 3 6 14 2 2 5 6 ...
## $ TIPO_ACCIDENTE : num [1:97916] 12 12 4 19 12 12 8 15 19 16 ...
## $ TOTAL_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOTAL_HG24H : num [1:97916] 1 0 0 0 0 0 0 0 0 0 ...
## $ TOTAL_HL24H : num [1:97916] 1 1 3 1 3 1 1 1 2 1 ...
## $ TOTAL_VICTIMAS_24H : num [1:97916] 2 1 3 1 3 1 1 1 2 1 ...
## $ TOTAL_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOTAL_HG30DF : num [1:97916] 1 0 0 0 0 0 0 0 0 0 ...
## $ TOTAL_HL30DF : num [1:97916] 1 1 3 1 3 1 1 1 2 1 ...
## $ TOTAL_VICTIMAS_30DF : num [1:97916] 2 1 3 1 3 1 1 1 2 1 ...
## $ TOTAL_VEHICULOS : num [1:97916] 1 1 2 1 2 1 2 1 1 1 ...
## $ TOT_PEAT_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_BICI_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_CICLO_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_MOTO_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_TUR_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_FURG_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_CAM_MENOS3500_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_CAM_MAS3500_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_BUS_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_OTRO_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_SINESPECIF_MU24H : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_PEAT_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_BICI_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_CICLO_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_MOTO_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_TUR_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_FURG_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_CAM_MENOS3500_MU30DF: num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_CAM_MAS3500_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_BUS_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_VMP_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_OTRO_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ TOT_SINESPECIF_MU30DF : num [1:97916] 0 0 0 0 0 0 0 0 0 0 ...
## $ NUDO : num [1:97916] 2 2 1 2 2 2 2 2 1 2 ...
## $ NUDO_INFO : num [1:97916] NA NA 2 NA NA NA NA NA 999 NA ...
## $ CARRETERA_CRUCE : chr [1:97916] NA NA NA NA ...
## $ PRIORI_NORMA : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_AGENTE : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_SEMAFORO : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_VERT_STOP : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_VERT_CEDA : num [1:97916] 999 999 1 999 999 999 999 999 0 999 ...
## $ PRIORI_HORIZ_STOP : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_HORIZ_CEDA : num [1:97916] 999 999 1 999 999 999 999 999 0 999 ...
## $ PRIORI_MARCAS : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_PEA_NO_ELEV : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_PEA_ELEV : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_MARCA_CICLOS : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_CIRCUNSTANCIAL : num [1:97916] 999 999 0 999 999 999 999 999 0 999 ...
## $ PRIORI_OTRA : num [1:97916] 999 999 0 999 999 999 999 999 1 999 ...
## $ CONDICION_NIVEL_CIRCULA : num [1:97916] 1 1 1 1 1 1 1 1 1 1 ...
## $ CONDICION_FIRME : num [1:97916] 1 3 1 3 1 1 1 1 1 1 ...
## $ CONDICION_ILUMINACION : num [1:97916] 6 1 6 1 1 6 6 1 5 1 ...
## $ CONDICION_METEO : num [1:97916] 1 3 2 3 1 1 1 1 1 1 ...
## $ CONDICION_NIEBLA : num [1:97916] NA NA NA NA NA NA NA NA NA NA ...
## $ CONDICION_VIENTO : num [1:97916] NA NA NA NA NA NA NA NA NA NA ...
## $ VISIB_RESTRINGIDA_POR : num [1:97916] 1 1 1 1 1 1 1 1 1 1 ...
## $ ACERA : num [1:97916] 998 998 998 998 998 998 998 998 998 998 ...
## $ TRAZADO_PLANTA : num [1:97916] 3 2 2 1 3 2 1 3 1 3 ...Una vez revisado, organizamos los datos que vamos a querer utilizar para nuestro estudio.
- COD_PROVINCIA - Columna que usaremos para filtrar.
- GRAV_ACC - Utilizaremos las variables TOTAL_MU24H, TOTAL_HG24H y TOTAL_HL24H para generar la variable.
- Predictoras de interés:
- TIPO_VIA
- CONDICION_LUMINICA
- TIPO_ACCIDENTE
- DIA_SEMANA
Preprocesamiento y gestión de características
Limpieza
En este paso revisaremos si tenemos datos vacios o nulos.
print('Valores nulos')## [1] "Valores nulos"colSums(is.na(accidentes))## ID_ACCIDENTE ANYO MES
## 0 0 0
## DIA_SEMANA HORA COD_PROVINCIA
## 0 0 0
## COD_MUNICIPIO ISLA ZONA
## 0 92709 0
## ZONA_AGRUPADA CARRETERA KM
## 0 0 46972
## SENTIDO_1F TITULARIDAD_VIA TIPO_VIA
## 0 0 0
## TIPO_ACCIDENTE TOTAL_MU24H TOTAL_HG24H
## 0 0 0
## TOTAL_HL24H TOTAL_VICTIMAS_24H TOTAL_MU30DF
## 0 0 0
## TOTAL_HG30DF TOTAL_HL30DF TOTAL_VICTIMAS_30DF
## 0 0 0
## TOTAL_VEHICULOS TOT_PEAT_MU24H TOT_BICI_MU24H
## 0 0 0
## TOT_CICLO_MU24H TOT_MOTO_MU24H TOT_TUR_MU24H
## 0 0 0
## TOT_FURG_MU24H TOT_CAM_MENOS3500_MU24H TOT_CAM_MAS3500_MU24H
## 0 0 0
## TOT_BUS_MU24H TOT_OTRO_MU24H TOT_SINESPECIF_MU24H
## 0 0 0
## TOT_PEAT_MU30DF TOT_BICI_MU30DF TOT_CICLO_MU30DF
## 0 0 0
## TOT_MOTO_MU30DF TOT_TUR_MU30DF TOT_FURG_MU30DF
## 0 0 0
## TOT_CAM_MENOS3500_MU30DF TOT_CAM_MAS3500_MU30DF TOT_BUS_MU30DF
## 0 0 0
## TOT_VMP_MU30DF TOT_OTRO_MU30DF TOT_SINESPECIF_MU30DF
## 0 0 0
## NUDO NUDO_INFO CARRETERA_CRUCE
## 0 60207 95248
## PRIORI_NORMA PRIORI_AGENTE PRIORI_SEMAFORO
## 0 0 0
## PRIORI_VERT_STOP PRIORI_VERT_CEDA PRIORI_HORIZ_STOP
## 0 0 0
## PRIORI_HORIZ_CEDA PRIORI_MARCAS PRIORI_PEA_NO_ELEV
## 0 0 0
## PRIORI_PEA_ELEV PRIORI_MARCA_CICLOS PRIORI_CIRCUNSTANCIAL
## 0 0 0
## PRIORI_OTRA CONDICION_NIVEL_CIRCULA CONDICION_FIRME
## 0 0 0
## CONDICION_ILUMINACION CONDICION_METEO CONDICION_NIEBLA
## 0 0 90136
## CONDICION_VIENTO VISIB_RESTRINGIDA_POR ACERA
## 97653 0 0
## TRAZADO_PLANTA
## 0print('Valores en blanco')## [1] "Valores en blanco"colSums(accidentes == "")## ID_ACCIDENTE ANYO MES
## 0 0 0
## DIA_SEMANA HORA COD_PROVINCIA
## 0 0 0
## COD_MUNICIPIO ISLA ZONA
## 0 NA 0
## ZONA_AGRUPADA CARRETERA KM
## 0 0 NA
## SENTIDO_1F TITULARIDAD_VIA TIPO_VIA
## 0 0 0
## TIPO_ACCIDENTE TOTAL_MU24H TOTAL_HG24H
## 0 0 0
## TOTAL_HL24H TOTAL_VICTIMAS_24H TOTAL_MU30DF
## 0 0 0
## TOTAL_HG30DF TOTAL_HL30DF TOTAL_VICTIMAS_30DF
## 0 0 0
## TOTAL_VEHICULOS TOT_PEAT_MU24H TOT_BICI_MU24H
## 0 0 0
## TOT_CICLO_MU24H TOT_MOTO_MU24H TOT_TUR_MU24H
## 0 0 0
## TOT_FURG_MU24H TOT_CAM_MENOS3500_MU24H TOT_CAM_MAS3500_MU24H
## 0 0 0
## TOT_BUS_MU24H TOT_OTRO_MU24H TOT_SINESPECIF_MU24H
## 0 0 0
## TOT_PEAT_MU30DF TOT_BICI_MU30DF TOT_CICLO_MU30DF
## 0 0 0
## TOT_MOTO_MU30DF TOT_TUR_MU30DF TOT_FURG_MU30DF
## 0 0 0
## TOT_CAM_MENOS3500_MU30DF TOT_CAM_MAS3500_MU30DF TOT_BUS_MU30DF
## 0 0 0
## TOT_VMP_MU30DF TOT_OTRO_MU30DF TOT_SINESPECIF_MU30DF
## 0 0 0
## NUDO NUDO_INFO CARRETERA_CRUCE
## 0 NA NA
## PRIORI_NORMA PRIORI_AGENTE PRIORI_SEMAFORO
## 0 0 0
## PRIORI_VERT_STOP PRIORI_VERT_CEDA PRIORI_HORIZ_STOP
## 0 0 0
## PRIORI_HORIZ_CEDA PRIORI_MARCAS PRIORI_PEA_NO_ELEV
## 0 0 0
## PRIORI_PEA_ELEV PRIORI_MARCA_CICLOS PRIORI_CIRCUNSTANCIAL
## 0 0 0
## PRIORI_OTRA CONDICION_NIVEL_CIRCULA CONDICION_FIRME
## 0 0 0
## CONDICION_ILUMINACION CONDICION_METEO CONDICION_NIEBLA
## 0 0 NA
## CONDICION_VIENTO VISIB_RESTRINGIDA_POR ACERA
## NA 0 0
## TRAZADO_PLANTA
## 0Vemos que ISLA tiene 92709 valores en NA ya que no hay isla en esos lugares y se podria modificar y poner a 0, KM tiene 46972 valores en NA, NUDO_INFO tiene 60207 valores en NA, CARRETERA_CRUCE tiene 95248 valores en NA y CONDICION_NIEBLA tiene 90136 en NA. En cambio, no hay ningun valor en blanco. Como ninguno de estos valores va a ser utilizado por el momento podemos seguir adelante.
Lo primero que haremos será crear la variable objetivo y la factorización de variables clave.
if (!require('dplyr')) install.packages('dplyr'); library('dplyr')## Loading required package: dplyr##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union# Barcelona (8), Girona (17), Lleida (25), Tarragona (43).
codigos_catalunya <- c(8, 17, 25, 43)
# Filtrar el conjunto de datos solo para Cataluña
accidentData_CAT <- accidentes %>%
filter(COD_PROVINCIA %in% codigos_catalunya)
# Creamos la variable GRAV_ACC
# Grave/Mortal = 1: TOTAL_MU24H > 0 o TOTAL_HG24H > 0
# Leve = 1: Resto de casos
accidentData_CAT <- accidentData_CAT %>%
mutate(
GRAV_ACC = ifelse(TOTAL_MU24H > 0 | TOTAL_HG24H > 0, 1, 0),
GRAV_ACC = factor(GRAV_ACC, levels = c(0,1), labels = c("Leve", "Grave/Mortal")),
DIA_SEMANA = factor(DIA_SEMANA),
TIPO_VIA = factor(TIPO_VIA),
CONDICION_ILUMINACION = factor(CONDICION_ILUMINACION)
)Ahora queremos confirmar que, realmente hay un gran desbalance entre gravedad Leve y Grave/Mortal y para ello vamos a crear histogramas y describir los valores para ver los datos en general de estos atributos para hacer una primera aproximación a los datos:
if (!require('ggplot2')) install.packages('ggplot2'); library('ggplot2')## Loading required package: ggplot2# Visualización de Leve vs Grave/Mortal
gravedad_summary <- accidentData_CAT %>%
group_by(GRAV_ACC) %>%
summarise( Conteo = n(), Porcentaje = Conteo / nrow(accidentData_CAT) * 100
)
print(gravedad_summary)## # A tibble: 2 × 3
## GRAV_ACC Conteo Porcentaje
## <fct> <int> <dbl>
## 1 Leve 21476 93.1
## 2 Grave/Mortal 1583 6.86# Gráfico de barras de GRAV_ACC
ggplot(gravedad_summary, aes(x = GRAV_ACC, y = Porcentaje, fill = GRAV_ACC)) +
geom_bar(stat = "identity") +
theme_minimal() +
labs(title = "Distribución de Gravedad de Accidentes en Cataluña (2022)", x = "Gravedad", y = "Porcentaje de Accidentes") + scale_fill_manual(values = c("Leve" = "steelblue", "Grave/Mortal" = "firebrick"))
Observaciones:
Como se puede observar se confirma un desbalance de clases en el que Grave/Mortal es minoritario, representando un pequeño porcentaje total. Esto confirma que se debe usar la técnica SMOTE en la fase de modelado.
Vamos a ver la correlacion de GRAV_ACC con factores de riesgo como el dia de la semana o la condición lumínica.
# 1. Porcentaje de accidentes Graves/Mortales por dia de la semana
gravedad_dia <- accidentData_CAT %>%
group_by(DIA_SEMANA) %>%
summarise(
Total = n(),
Graves_Mortales = sum(GRAV_ACC == "Grave/Mortal")
) %>%
mutate(Pct_Grave_Mortal = 100 * Graves_Mortales / Total)
# Gráfico de la proporción de accidentes graves por día
ggplot(gravedad_dia, aes(x = DIA_SEMANA, y = Pct_Grave_Mortal, fill = DIA_SEMANA)) +
geom_bar(stat = "identity") +
theme_minimal() +
labs(title = "Proporción de Accidentes Graves/Mortales por Día (Cataluña)",
x = "Día de la Semana (Código)",
y = "% Grave/Mortal sobre el total del día") +
theme(legend.position = "none")
Vemos que las dos columnas más importantes son las del fin de semana ya que se va a una mayor velocidad debido a que puede haber un menor trafico en carretera e induce a ir a mayor velocidad.
# 2. Porcentaje de accidentes Graves/Mortales por CONDICIÓN DE ILUMINACIÓN
gravedad_ilum <- accidentData_CAT %>%
group_by(CONDICION_ILUMINACION) %>%
summarise(
Total = n(),
Graves_Mortales = sum(GRAV_ACC == "Grave/Mortal")
) %>%
mutate(Pct_Grave_Mortal = 100 * Graves_Mortales / Total)
gravedad_ilum_plot <- gravedad_ilum %>%
mutate(
CONDICION_ILUMINACION = factor(CONDICION_ILUMINACION)
) %>%
filter(CONDICION_ILUMINACION != 999)
ggplot(gravedad_ilum_plot,
aes(x = reorder(CONDICION_ILUMINACION, -Pct_Grave_Mortal),
y = Pct_Grave_Mortal,
fill = CONDICION_ILUMINACION)) +
geom_bar(stat = "identity") +
theme_minimal() +
labs(title = "Riesgo Relativo de Gravedad por Condición Lumínica en Cataluña",
subtitle = "Porcentaje de accidentes Grave/Mortal sobre el total de cada condición\nCódigos DGT: 1=Día, 2=Crepusuclo sin luz artificial, 3=Crepúsculo con luz artificial, \n4=Sin luz natural y con luz artificial encendida, 5=Sin luz natural y con luz artificial apagada, \n6=Sin luz natural ni artificial",
x = "Código de Condición de Iluminación",
y = "% de Accidentes Grave/Mortal") +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5),
legend.position = "none") +
scale_fill_brewer(palette = "Reds")
Vemos que la mayor cantidad de accidentes Graves/Mortales son de noche y al atardecer debido a la reducción de la visibilidad o condicionantes externos como el alcohol y drogas.
# Heatmap de accidentes por día y hora (similar al ejemplo)
# Filtramos HORA para evitar códigos de error (ej. 99)
accidentes_heatmap_CAT <- accidentData_CAT %>%
filter(HORA >= 0 & HORA <= 23) %>%
group_by(DIA_SEMANA, HORA) %>%
summarise(num_accidentes = n(), .groups = "drop")
# El día de la semana (1-7) se ordena automáticamente por ser factor/numérico.
ggplot(accidentes_heatmap_CAT, aes(x = HORA, y = DIA_SEMANA, fill = num_accidentes)) +
geom_tile(color = "white") +
scale_fill_gradient(low = "lightyellow", high = "red") +
theme_minimal() +
labs(title = "Heatmap de Frecuencia de Accidentes en Cataluña por Día y Hora",
x = "Hora del día",
y = "Día de la semana (código)",
fill = "Nº de accidentes")
Este gráfico revela los momentos de mayor frecuencia de accidentes. Se observa alta intensidad durante las horas punta de las mañanas y tardes de los días laborables, asi como lo que parece a la salida del colegio (15-16) de los dias entre semana y en las madrugadas de los fines de semana.
Construcción de conjunto de datos final
Codificación
Seguidamente vayamos a asignar un 1 por accidentes que se producen de madrugada (01h a 07h) y un 0 para el resto de franja horaria, la denominaremos madrugada. Después la utilizaremos para ver cómo se distribuyen los accidentes en las dos franjas horarias.
accidentData_CAT_AUX=subset(accidentData_CAT, accidentData_CAT$HORA <= 24)
accidentData_CAT$madrugada <- NA
accidentData_CAT$madrugada[accidentData_CAT_AUX$HORA >=1 & accidentData_CAT_AUX$HORA <= 7] <- 1
accidentData_CAT$madrugada[accidentData_CAT_AUX$HORA ==0 | accidentData_CAT_AUX$HORA >7 ] <- 0
counts <- table(accidentData_CAT$madrugada)
barplot(prop.table(counts),col=c("green","red"),legend.texto=c("Resto del día","Madrugada"),ylim=c(0,1), main="Distribución de accidentes la madrugada y resto del día",xlab="0 Resto del día 1 Madrugada",ylab="Porcentaje" )## Warning in plot.window(xlim, ylim, log = log, ...): "legend.texto" is not a
## graphical parameter## Warning in axis(if (horiz) 2 else 1, at = at.l, labels = names.arg, lty =
## axis.lty, : "legend.texto" is not a graphical parameter## Warning in title(main = main, sub = sub, xlab = xlab, ylab = ylab, ...):
## "legend.texto" is not a graphical parameter## Warning in axis(if (horiz) 1 else 2, cex.axis = cex.axis, ...): "legend.texto"
## is not a graphical parameter
Discretización
Ahora añadiremos un campo nuevo a los datos. Este campo contendrá el valor de la hora del accidente discretizada.
summary(accidentData_CAT_AUX[,"HORA"])## HORA
## Min. : 0.0
## 1st Qu.:10.0
## Median :14.0
## Mean :13.7
## 3rd Qu.:18.0
## Max. :23.0Discretizamos con intervalos.
accidentData_CAT_AUX["segmento_horario"] <- cut(accidentData_CAT_AUX$HORA, breaks = c(0,4,11,14,18,22), labels = c("Madrugada", "Mañana", "Mediodía", "Anochecer","Noche"))Observamos los datos discretizados y construimos un gráfico para analizar cómo se agrupan los accidentes.
head(accidentData_CAT_AUX$segmento_horario)## [1] Anochecer Noche Noche Mañana Mediodía Madrugada
## Levels: Madrugada Mañana Mediodía Anochecer Nocheplot(accidentData_CAT_AUX$segmento_horario,main="Número de accidentes por segmento horario",xlab="Segmento horario", ylab="Cantidad",col = "ivory")
Por último en la discretización, aplicamos la discretización basada en el algoritmo k-means a la variable TOTAL_VEHICULOS para crear grupos de complejidad vehicular.
if (!require('arules')) install.packages('arules'); library('arules')## Loading required package: arules## Loading required package: Matrix##
## Attaching package: 'arules'## The following object is masked from 'package:dplyr':
##
## recode## The following objects are masked from 'package:base':
##
## abbreviate, writeset.seed(2)
# Generar la tabla de frecuencias de los clusters
print("Frecuencia de Accidentes por Cluster de Vehículos Implicados:")## [1] "Frecuencia de Accidentes por Cluster de Vehículos Implicados:"table(discretize(accidentData_CAT_AUX$TOTAL_VEHICULOS, "cluster" ))##
## [1,2.48) [2.48,6.25) [6.25,27]
## 21236 1804 19hist(accidentData_CAT_AUX$TOTAL_VEHICULOS,
main="Frecuencia de Accidentes por Vehículos Implicados con k-means",
xlab="Vehículos implicados",
ylab="Cantidad",
col = "ivory")
abline(v = discretize(accidentData_CAT_AUX$TOTAL_VEHICULOS, method = "cluster", onlycuts = TRUE),
col = "red")
# Asignar la nueva variable al dataset auxiliar
accidentData_CAT_AUX$VEHICULOS_KM <- discretize(accidentData_CAT_AUX$TOTAL_VEHICULOS, "cluster" )Como se puede observar tanto en la tabla como en el grafico, la mayoria de accidentes implican de 1 a 3 coches.
Normalización
Normalizamos las variables numericas de victimas y vehiculos para que contribuyan por igual en el analisis de PCA.
nor <- function(x) { (x -min(x, na.rm = TRUE))/(max(x, na.rm = TRUE)-min(x, na.rm = TRUE))}
# columnas numéricas a normalizar
n_norm <- c("TOTAL_MU24H", "TOTAL_HG24H", "TOTAL_HL24H",
"TOTAL_VICTIMAS_24H", "TOTAL_VEHICULOS")
# creamos un nuevo dataset normalizado solo con las columnas clave
accidentData_CAT_nor <- accidentData_CAT_AUX %>%
select(all_of(n_norm)) %>%
as.data.frame(lapply(., nor))hist(accidentData_CAT_nor$TOTAL_MU24H,xlab="Muertos 24H", col="ivory",ylab="Cantidad", main="Distribución Normalizada de Muertos 24H")
hist(accidentData_CAT_nor$TOTAL_VICTIMAS_24H,xlab="Muertos",ylab="Cantidad",col="ivory", main="Distribución Normalizada de Víctimas 24H")
Proceso de PCA
El objetivo del PCA es reducir la dimensionalidad de las variables de impacto de los accidentes (Víctimas, Heridos, Muertes, Vehículos). Como estas variables suelen aumentar juntas, son redundantes. El PCA nos permite crear un “Índice de Severidad” (PC1) que encapsula casi toda esa información en un solo predictor, mejorando el modelado.
pca.acc <- prcomp(accidentData_CAT_nor)
summary(pca.acc)## Importance of components:
## PC1 PC2 PC3 PC4 PC5
## Standard deviation 1.1240 0.6900 0.35761 0.12068 3.845e-13
## Proportion of Variance 0.6713 0.2530 0.06795 0.00774 0.000e+00
## Cumulative Proportion 0.6713 0.9243 0.99226 1.00000 1.000e+00Hemos aplicado el PCA a los datos normalizados y el resumen nos muestra que PC1 explica el 67% de la varianza lo cual es un muy buen indicio.
Utilizando el metodo de Kaiser para decidir cuales de las variables obtenidas serán escogidas, mantendremos todas aquellas cuya varianza sea superior a 1.
var_acc <- pca.acc$sdev^2
var_acc## [1] 1.263443e+00 4.760798e-01 1.278883e-01 1.456263e-02 1.478720e-25Con los resultados obtenidos es muy complicado decidir cuáles son los componentes principales componentes a escoger. Por lo tanto, el siguiente paso es escalar los datos y volver a calcular la varianza para ver qué datos selecciona.
Aunque los datos ya estaban normalizados, para que el criterio de Kaiser sea valido para que PCA funcione óptimamente, los datos deben estar escalados y, por lo tanto, debemos escalar los datos.
# Escalamos los datos
acc_scale <- scale(accidentData_CAT_nor)
# Calculamos las componentes principales
pca.acc_scale <- prcomp(acc_scale)
# Mostramos la varianza de dichas variables:
var_acc_scale <- pca.acc_scale$sdev^2
head(var_acc_scale)## [1] 2.039316e+00 1.162058e+00 9.302256e-01 8.684004e-01 1.967194e-26Aplicando el método de Kaiser solo retendremos los componentes principales 1 y 2, ya que ambos tienen una varianza superior a 1. Estos dos nuevos indices son suficientes para representar la complejidad de las cinco variables de impacto.
Continuamos con el análisis de los componentes principales. Después de aplicar el método Káiser se han seleccionado los 2 componentes principales.
if (!require('factoextra')) install.packages('factoextra'); library('factoextra')## Loading required package: factoextra## Welcome! Want to learn more? See two factoextra-related books at https://goo.gl/ve3WBavar <- get_pca_var(pca.acc_scale)
var## Principal Component Analysis Results for variables
## ===================================================
## Name Description
## 1 "$coord" "Coordinates for the variables"
## 2 "$cor" "Correlations between variables and dimensions"
## 3 "$cos2" "Cos2 for the variables"
## 4 "$contrib" "contributions of the variables"Retuvimos PC1 y PC2 basandonos en el criterio de Kaiser. Vamos a analizar su contenido para darles significado.
head(var$coord[,1:2])## Dim.1 Dim.2
## TOTAL_MU24H 0.01669318 0.659605338
## TOTAL_HG24H 0.04322462 0.808284903
## TOTAL_HL24H 0.95214018 -0.214389187
## TOTAL_VICTIMAS_24H 0.96160676 0.166361691
## TOTAL_VEHICULOS 0.45377331 -0.003955029Esto que acabamos de hacer, muestra las coordenadas de cada variable en el nuevo sistema. Esto nos dice si una variable está fuertemente relacionada con un componente.
head(var$contrib[,1:2])## Dim.1 Dim.2
## TOTAL_MU24H 0.01366450 37.440394799
## TOTAL_HG24H 0.09161738 56.221319055
## TOTAL_HL24H 44.45466218 3.955285617
## TOTAL_VICTIMAS_24H 45.34303099 2.381654446
## TOTAL_VEHICULOS 10.09702494 0.001346082En cambio contrib muestra la contribución de cada variable en porcentaje a la formación de PC1 y PC2. En nuestro caso, PC1 es el componente más importante y representa principalmente la Escala del Accidente medida por la siniestralidad leve.
En cambio, PC2 representa el Grado de Severidad del Accidente definido por HG24H y MU24H y en el contexto de la DGT, es el indice de Daño/Costo Humano.
PCA ha logrado separar dos conceptos: - PC1 mide la frecuencia y magnitud de las colisiones. - PC2 mide la calidad o letalidad del accidente.
Calidad de representación
La calidad de representación de las variables en el mapa de factores se denomina cos2.
head(var$cos2[,1:2])## Dim.1 Dim.2
## TOTAL_MU24H 0.0002786623 4.350792e-01
## TOTAL_HG24H 0.0018683677 6.533245e-01
## TOTAL_HL24H 0.9065709193 4.596272e-02
## TOTAL_VICTIMAS_24H 0.9246875642 2.767621e-02
## TOTAL_VEHICULOS 0.2059102180 1.564226e-05Un valor alto significa que PC1 y PC2 son casi perfectos sustitutos de esa variable. Por lo que, PC1 es casi un sustituto perfecto para las variables de magnitud (TOTAL_VICTIMAS_24H y TOTAL_HL24H) ya que están cerca del 1 y pueden ser eliminadas sin perdida de información relevante. PC2 en cambio, no tiene unas buenas métricas y aunque sea la mejor dimension para capturar las variables de gravedad, perderiamos mucha información.
Aqui podemos verlo visualmente:
if(!require("corrplot")) install.packages("corrplot"); library("corrplot")## Loading required package: corrplot## corrplot 0.95 loadedcorrplot(var$cos2[,1:2], is.corr=FALSE,
tl.col="black", tl.srt=45,
main="Calidad de Representación (cos²) en PC1 y PC2")
fviz_cos2(pca.acc_scale, choice = "var", axes = 1:2,
title = "Varianza Explicada por Variables (PC1 + PC2)")
El proceso de PCA ha transformado cinco variables en dos nuevos predictores PC1 (67% de la varianza) y PC2 (23% de la varianza).
Contribución
Las contribuciones de las variables en la contabilización de la variabilidad de un determinado componente principal se expresan en porcentaje.
Las variables que están correlacionadas con PC1 y PC2 son las más importantes para explicar la variabilidad en el conjunto de datos.
print("Contribución de las variables a los PC1 y PC2:")## [1] "Contribución de las variables a los PC1 y PC2:"head(var$contrib[, 1:2]) # Mostramos solo los dos componentes retenidos## Dim.1 Dim.2
## TOTAL_MU24H 0.01366450 37.440394799
## TOTAL_HG24H 0.09161738 56.221319055
## TOTAL_HL24H 44.45466218 3.955285617
## TOTAL_VICTIMAS_24H 45.34303099 2.381654446
## TOTAL_VEHICULOS 10.09702494 0.001346082Cuanto más grande sea el valor de la contribucion, mas contribucion habra en el componente.
if(!require("corrplot")) install.packages("corrplot"); library("corrplot")
corrplot(var$contrib[, 1:2], is.corr=FALSE,
tl.col="black", tl.srt=45,
main="Contribución de Variables a PC1 y PC2")
if (!require('factoextra')) install.packages('factoextra'); library('factoextra')
fviz_pca_var(pca.acc_scale,
col.var = "contrib",
gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
repel = TRUE,
title = "Contribución de Variables a PC1 y PC2")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## ℹ The deprecated feature was likely used in the ggpubr package.
## Please report the issue at <https://github.com/kassambara/ggpubr/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the factoextra package.
## Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
El análisis confirma que laPC1es dominada por las contribuciones de TOTAL_VICTIMAS_24H y TOTAL_HL24H, por lo que este componente representa la Magnitud o Escala del accidente. Por otro lado, la PC2 es dominada por TOTAL_HG24H y TOTAL_MU24H, representando la verdadera Letalidad/Severidad del siniestro. Esto justifica retener ambos componentes.
Interpretación de los resultados
El análisis exploratorio y el preprocesamiento de los microdatos de accidentes con víctimas de la DGT de 2022, enfocados en Cataluña y en el objetivo de Clasificación de la Gravedad, han sentado las bases para el modelado.
Desbalance de clases en la variable objetivo, siendo la clase Grave/Motal minoritaria, requiriendo de tecnicas de remuestreo como SMOTE.
El análisis de riesgo relativo ha revelado que la probabilidad de que un accidente sea grave o mortal es significativamente mayor en condiciones de oscuridad y durante los fines de semana. El heatmap ha mostrado que los ficos de frecuencia absoluta ocurren en horas punta.
Se crearon variables categoricas de tiempo para mejorar la segmentación horaria.
Las cinco variables de impacto originales demostraron una alta correlación, lo que justifica la reducción de dimensionalidad para evitar multicolinealidad.
El PCA ha separado los conceptos de PC1 representadno la magnitud y PC2 representando la letalidad.
Ejercicio 3
Las funciones explore_all y explote_tbl son
funciones de analisis exploratorio de datos automatizado.
if(!require('explore')) install.packages('explore'); library('explore')## Loading required package: explorepath = 'accident.CSV'
accidentData <- read.csv(path, row.names=NULL)
# Inspecciona las 81 variables del dataset y genera un resumen tabular y grafico de todas las variables
explore_all(accidentData)
# Inspecciona la variable y su relacion con el dataset, en este caso con STATE
explore_tbl(accidentData, "STATE")
Estas dos funciones son útiles en la fase de comprension de datos ya que permiten realizar el diagnostico inicial de la calidad y distribución de los datos de forma automatizada.