Trabajo Final - Métodos de Regresión.
Yeison Buitrago
2025-12-12
Trabajo 3 - Métodos de Regresión.
library(sf)## Linking to GEOS 3.13.0, GDAL 3.8.5, PROJ 9.5.1; sf_use_s2() is TRUElibrary(readxl)
library(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, unionlibrary(ggplot2)Introducción
Vamos a utilizar el Anuario de Siniestralidad Vial de Bogotá 2023, el documento es un informe estadístico oficial que consolida el análisis de los siniestros viales en Bogotá, utilizando registros del sistema SIGAT e informes policiales (IPAT) para caracterizar la mortalidad y morbilidad en el tránsito. Los datos abarcan series históricas (2015-2024) y desgloses detallados del último año, clasificando a las víctimas (fallecidos y lesionados) por actor vial (con énfasis en el aumento de motociclistas), género, rango etario, localidad y variables temporales como hora y día de la semana. Además, el reporte incluye tasas de mortalidad por cada 100.000 habitantes, matrices de interacción de riesgos (incluyendo siniestros con el SITP) y estudios de caso específicos que aplican un enfoque diferencial y de género para evaluar el impacto social de la accidentalidad.
SIGAT_ANUARIO_2023 <- read_excel("/Users/usermac/Downloads/SIGAT_ANUARIO_2023.xlsx")## Warning: Expecting logical in AA2651 / R2651C27: got 'SI'## Warning: Expecting logical in AA3219 / R3219C27: got 'SI'## Warning: Expecting logical in AA9357 / R9357C27: got 'SI'## Warning: Expecting logical in AA13730 / R13730C27: got 'SI'head(SIGAT_ANUARIO_2023)## # A tibble: 6 × 38
## Codigo_Accidente Formulario Longitud Latitud Direccion Fecha_Acc
## <dbl> <chr> <dbl> <dbl> <chr> <dttm>
## 1 10591851 A001571139 -74.2 4.60 AV AVENIDA C… 2023-05-26 00:00:00
## 2 10591854 A001571298 -74.2 4.57 KR 44 D - DG… 2023-05-31 00:00:00
## 3 10592149 A001570489 -74.1 4.62 CL 26 - KR 2… 2023-06-08 00:00:00
## 4 10592358 A001572190 -74.1 4.59 CL 6 - KR 5 … 2023-06-14 00:00:00
## 5 10592361 A001570128 -74.0 4.76 CL 187 A - K… 2023-06-14 00:00:00
## 6 10592370 A001571474 -74.1 4.74 CL 139 - KR … 2023-06-13 00:00:00
## # ℹ 32 more variables: AA_Acc <dbl>, MM_Acc <chr>, DD_Mes_Acc <dbl>,
## # Dia_Semana_Acc <chr>, Hora_Acc <dbl>, Min_Acc <dbl>, Localidad <chr>,
## # Clase_Acc <chr>, Elemento_Choque <chr>, Tipo_Objeto_Fijo <chr>,
## # Gravedad_Indicador_Tradicional <chr>, Gravedad_indicador_30d <chr>,
## # Con_Bicicleta <chr>, Con_Carga <chr>, Con_Embriaguez <chr>,
## # Con_Huecos <chr>, Con_Menores <chr>, Con_Moto <chr>, Con_Peaton <chr>,
## # Con_Persona_Mayor <chr>, Con_Rutas <lgl>, Con_Tpi <chr>, Con_Tpp <chr>, …Variables a usar.
- Con_Embriaguez: Señala si en el siniestro hubo alguna hipótesis relacionada con embriaguez
- Con_Moto: Señala si en el siniestro hubo una motocicleta involucrada
- Con_Huecos: Señala si en el siniestro hubo alguna hipótesis relacionada con huecos
- Con_Menores: Señala si en el siniestro hubo una persona menor de edad involucrada
- Con_Peaton: Señala si en el siniestro hubo un peatón involucrado
- Gravedad_indicador_30d: Gravedad del siniestro (Con Muertos, Con Heridos ó Solo daños) tomando en cuenta los fallecimientos que ocurren dentro de los 30 días posteriores al siniestro
- Clase_Acc: Clase de accidente (Choque, atropello, Volcamiento, caída de ocupante, incendio, otro, autolesión)
- Localidad : Nombre de la localidad del distrito capital donde ocurrió el siniestro vial
df_to_use <- SIGAT_ANUARIO_2023 %>%
select(Con_Embriaguez, Con_Moto,Con_Huecos, Con_Menores, Con_Peaton, Gravedad_indicador_30d,Clase_Acc, Localidad, )
head(df_to_use)## # A tibble: 6 × 8
## Con_Embriaguez Con_Moto Con_Huecos Con_Menores Con_Peaton
## <chr> <chr> <chr> <chr> <chr>
## 1 <NA> SI <NA> <NA> <NA>
## 2 <NA> SI <NA> <NA> SI
## 3 <NA> SI <NA> <NA> <NA>
## 4 <NA> <NA> <NA> SI <NA>
## 5 <NA> SI <NA> <NA> <NA>
## 6 <NA> SI <NA> <NA> SI
## # ℹ 3 more variables: Gravedad_indicador_30d <chr>, Clase_Acc <chr>,
## # Localidad <chr>Vamos a definir dos conjuntos de variables, unas categoricas, que serán a quellas categoricas que más de dos valores.
categoricas <- c('Gravedad_indicador_30d','Clase_Acc', 'Localidad')Binarias, que serán aquellas que solo toman dos variables categoricas. Revismamos carateristicas generales sobre las variables.
binarias <- c('Con_Embriaguez','Con_Moto', 'Con_Huecos', 'Con_Menores', 'Con_Peaton')
summary(df_to_use[, binarias])## Con_Embriaguez Con_Moto Con_Huecos Con_Menores
## Length:14106 Length:14106 Length:14106 Length:14106
## Class :character Class :character Class :character Class :character
## Mode :character Mode :character Mode :character Mode :character
## Con_Peaton
## Length:14106
## Class :character
## Mode :characterTransformamos las variables a 1 si es un SI, es decir si el siniestro ocurrio Con_Embriaguez, aparecera un si.
conteo <- sapply(df_to_use[, binarias], function(x) {
sum(x == "Si" | x == "SI" | x == 1, na.rm = TRUE)
})
print(conteo)## Con_Embriaguez Con_Moto Con_Huecos Con_Menores Con_Peaton
## 554 8417 171 1325 2891convert <- function(columna) {
as.integer(columna %in% "SI")
}binarias <- c('Con_Embriaguez','Con_Moto', 'Con_Huecos', 'Con_Menores', 'Con_Peaton')for (i in binarias) {
df_to_use[[i]] <- as.factor(convert(df_to_use[[i]]))
}conteo <- sapply(df_to_use[, binarias], function(x) {
sum(x == "Si" | x == "SI" | x == 1, na.rm = TRUE)
})
print(conteo)## Con_Embriaguez Con_Moto Con_Huecos Con_Menores Con_Peaton
## 554 8417 171 1325 2891EDA
colores <- c("lightgrey", "red")
for (i in binarias) {
tabla <- table(df_to_use$Con_Embriaguez, df_to_use[[i]])
barplot(tabla,
beside = TRUE,
main = paste("Embriaguez según:", i),
xlab = i,
ylab = "Cantidad de Accidentes",
col = colores,
legend = rownames(tabla),
args.legend = list(title = "Con Embriaguez", x = "topright"))
}
for (i in binarias) {
print(i)
plot(df_to_use$Con_Embriaguez ~df_to_use[[i]], main = paste("Con Embriaguez vs", i), ylab = i)
}## [1] "Con_Embriaguez"
## [1] "Con_Moto"
## [1] "Con_Huecos"
## [1] "Con_Menores"
## [1] "Con_Peaton"
data_embriaguez <- df_to_use %>%
filter(Con_Embriaguez == 1)
ggplot(data_embriaguez, aes(x = Localidad, fill = Clase_Acc)) +
geom_bar() +
facet_wrap(~ Gravedad_indicador_30d, scales = "free_x") +
coord_flip() +
labs(title = "Total de Accidentes con Embriaguez",
subtitle = "Distribución por Localidad, Clase y Gravedad",
x = "Localidad",
y = "Cantidad de casos positivos",
fill = "Clase de Accidente") +
theme_bw()
### Modelo. Ajustamos un modelo logit, ya que queremos saber la
proabilidad de que el accidente dados unas variables sea con
embriaguez.
formula_logit <- Con_Embriaguez ~ Con_Moto + Con_Huecos + Con_Menores + Con_Peaton + Gravedad_indicador_30d + Clase_Acc + Localidad
modelo_logit <- glm(formula_logit,
data = df_to_use,
family = binomial(link = "logit"))
summary(modelo_logit)##
## Call:
## glm(formula = formula_logit, family = binomial(link = "logit"),
## data = df_to_use)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.10364 0.68505 -7.450 9.34e-14 ***
## Con_Moto1 -0.02813 0.11179 -0.252 0.8013
## Con_Huecos1 -1.27031 0.74672 -1.701 0.0889 .
## Con_Menores1 -0.17492 0.25451 -0.687 0.4919
## Con_Peaton1 -0.26383 0.57694 -0.457 0.6475
## Gravedad_indicador_30dCon Muertos 1.20223 0.23507 5.114 3.15e-07 ***
## Gravedad_indicador_30dSolo Daños 3.29987 0.11906 27.717 < 2e-16 ***
## Clase_AccCaida de ocupante -1.31945 0.92441 -1.427 0.1535
## Clase_AccChoque 0.68870 0.58610 1.175 0.2400
## Clase_AccIncendio -15.34376 2477.89393 -0.006 0.9951
## Clase_AccOtro -13.66810 297.56415 -0.046 0.9634
## Clase_AccVolcamiento 0.80665 0.64666 1.247 0.2122
## LocalidadBARRIOS UNIDOS 0.37158 0.41451 0.896 0.3700
## LocalidadBOSA 0.44419 0.38155 1.164 0.2444
## LocalidadCANDELARIA 0.05520 0.84666 0.065 0.9480
## LocalidadCHAPINERO 0.36696 0.41881 0.876 0.3809
## LocalidadCIUDAD BOLIVAR 0.68821 0.39083 1.761 0.0783 .
## LocalidadENGATIVA 0.50599 0.36985 1.368 0.1713
## LocalidadFONTIBON 0.26578 0.38470 0.691 0.4896
## LocalidadKENNEDY 0.50947 0.35585 1.432 0.1522
## LocalidadLOS MARTIRES 0.58170 0.41288 1.409 0.1589
## LocalidadPUENTE ARANDA 0.47991 0.36945 1.299 0.1940
## LocalidadRAFAEL URIBE URIBE 0.35456 0.41136 0.862 0.3887
## LocalidadSAN CRISTOBAL 0.18440 0.40709 0.453 0.6506
## LocalidadSANTA FE 0.14743 0.42816 0.344 0.7306
## LocalidadSUBA 0.74696 0.36647 2.038 0.0415 *
## LocalidadSUMAPAZ -12.62562 2774.65838 -0.005 0.9964
## LocalidadTEUSAQUILLO -0.28598 0.43890 -0.652 0.5147
## LocalidadTUNJUELITO 0.27899 0.42457 0.657 0.5111
## LocalidadUSAQUEN 0.23499 0.39327 0.598 0.5502
## LocalidadUSME 0.37994 0.43644 0.871 0.3840
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 4672.8 on 14105 degrees of freedom
## Residual deviance: 3349.4 on 14075 degrees of freedom
## AIC: 3411.4
##
## Number of Fisher Scoring iterations: 16best_model <- step(modelo_logit, direction = "both", trace = 1)## Start: AIC=3411.45
## Con_Embriaguez ~ Con_Moto + Con_Huecos + Con_Menores + Con_Peaton +
## Gravedad_indicador_30d + Clase_Acc + Localidad
##
## Df Deviance AIC
## - Localidad 19 3369.8 3393.8
## - Con_Moto 1 3349.5 3409.5
## - Con_Peaton 1 3349.7 3409.7
## - Con_Menores 1 3349.9 3409.9
## <none> 3349.4 3411.4
## - Con_Huecos 1 3353.5 3413.5
## - Clase_Acc 5 3379.0 3431.0
## - Gravedad_indicador_30d 2 4235.6 4293.6
##
## Step: AIC=3393.84
## Con_Embriaguez ~ Con_Moto + Con_Huecos + Con_Menores + Con_Peaton +
## Gravedad_indicador_30d + Clase_Acc
##
## Df Deviance AIC
## - Con_Moto 1 3370.0 3392.0
## - Con_Peaton 1 3370.1 3392.1
## - Con_Menores 1 3370.2 3392.2
## <none> 3369.8 3393.8
## - Con_Huecos 1 3374.2 3396.2
## + Localidad 19 3349.4 3411.4
## - Clase_Acc 5 3399.3 3413.3
## - Gravedad_indicador_30d 2 4252.7 4272.7
##
## Step: AIC=3391.99
## Con_Embriaguez ~ Con_Huecos + Con_Menores + Con_Peaton + Gravedad_indicador_30d +
## Clase_Acc
##
## Df Deviance AIC
## - Con_Peaton 1 3370.2 3390.2
## - Con_Menores 1 3370.3 3390.3
## <none> 3370.0 3392.0
## + Con_Moto 1 3369.8 3393.8
## - Con_Huecos 1 3374.3 3394.3
## + Localidad 19 3349.5 3409.5
## - Clase_Acc 5 3399.7 3411.7
## - Gravedad_indicador_30d 2 4443.6 4461.6
##
## Step: AIC=3390.23
## Con_Embriaguez ~ Con_Huecos + Con_Menores + Gravedad_indicador_30d +
## Clase_Acc
##
## Df Deviance AIC
## - Con_Menores 1 3370.6 3388.6
## <none> 3370.2 3390.2
## + Con_Peaton 1 3370.0 3392.0
## + Con_Moto 1 3370.1 3392.1
## - Con_Huecos 1 3374.5 3392.5
## + Localidad 19 3349.7 3407.7
## - Clase_Acc 5 3416.5 3426.5
## - Gravedad_indicador_30d 2 4449.8 4465.8
##
## Step: AIC=3388.6
## Con_Embriaguez ~ Con_Huecos + Gravedad_indicador_30d + Clase_Acc
##
## Df Deviance AIC
## <none> 3370.6 3388.6
## + Con_Menores 1 3370.2 3390.2
## + Con_Peaton 1 3370.3 3390.3
## + Con_Moto 1 3370.5 3390.5
## - Con_Huecos 1 3374.9 3390.9
## + Localidad 19 3350.2 3406.2
## - Clase_Acc 5 3417.7 3425.7
## - Gravedad_indicador_30d 2 4469.7 4483.7summary(best_model)##
## Call:
## glm(formula = Con_Embriaguez ~ Con_Huecos + Gravedad_indicador_30d +
## Clase_Acc, family = binomial(link = "logit"), data = df_to_use)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.9684 0.2157 -23.033 < 2e-16 ***
## Con_Huecos1 -1.3009 0.7466 -1.742 0.08144 .
## Gravedad_indicador_30dCon Muertos 1.2143 0.2345 5.179 2.23e-07 ***
## Gravedad_indicador_30dSolo Daños 3.3037 0.1027 32.174 < 2e-16 ***
## Clase_AccCaida de ocupante -1.0543 0.7403 -1.424 0.15442
## Clase_AccChoque 0.9352 0.2251 4.154 3.27e-05 ***
## Clase_AccIncendio -14.9987 2474.2348 -0.006 0.99516
## Clase_AccOtro -13.3937 299.7653 -0.045 0.96436
## Clase_AccVolcamiento 1.0483 0.3522 2.977 0.00291 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 4672.8 on 14105 degrees of freedom
## Residual deviance: 3370.6 on 14097 degrees of freedom
## AIC: 3388.6
##
## Number of Fisher Scoring iterations: 16model_optimo<- glm(formula = Con_Embriaguez ~ Con_Huecos + Gravedad_indicador_30d +
Clase_Acc, family = binomial(link = "logit"), data = df_to_use)
summary(model_optimo)##
## Call:
## glm(formula = Con_Embriaguez ~ Con_Huecos + Gravedad_indicador_30d +
## Clase_Acc, family = binomial(link = "logit"), data = df_to_use)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.9684 0.2157 -23.033 < 2e-16 ***
## Con_Huecos1 -1.3009 0.7466 -1.742 0.08144 .
## Gravedad_indicador_30dCon Muertos 1.2143 0.2345 5.179 2.23e-07 ***
## Gravedad_indicador_30dSolo Daños 3.3037 0.1027 32.174 < 2e-16 ***
## Clase_AccCaida de ocupante -1.0543 0.7403 -1.424 0.15442
## Clase_AccChoque 0.9352 0.2251 4.154 3.27e-05 ***
## Clase_AccIncendio -14.9987 2474.2348 -0.006 0.99516
## Clase_AccOtro -13.3937 299.7653 -0.045 0.96436
## Clase_AccVolcamiento 1.0483 0.3522 2.977 0.00291 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 4672.8 on 14105 degrees of freedom
## Residual deviance: 3370.6 on 14097 degrees of freedom
## AIC: 3388.6
##
## Number of Fisher Scoring iterations: 16diagnostico_logistico <- function(y_real, y_pred_prob, umbral = 0.5) {
y_pred_clase <- ifelse(y_pred_prob > umbral, 1, 0)
tabla_confusion <- table(Real = y_real, Predicho = y_pred_clase)
VP <- tabla_confusion[2, 2]
VN <- tabla_confusion[1, 1]
FP <- tabla_confusion[1, 2]
FN <- tabla_confusion[2, 1]
exactitud <- (VP + VN) / (VP + VN + FP + FN)
sensibilidad <- VP / (VP + FN)
especificidad <- VN / (VN + FP)
precision <- VP / (VP + FP)
f1_score <- 2 * (precision * sensibilidad) / (precision + sensibilidad)
library(pROC)
curva_roc <- roc(y_real, y_pred_prob)
auc_valor <- auc(curva_roc)
resultados <- list(
Matriz_Confusion = tabla_confusion,
Metricas = data.frame(
Exactitud = round(exactitud, 4),
Sensibilidad = round(sensibilidad, 4),
Especificidad = round(especificidad, 4),
Precision = round(precision, 4),
F1_Score = round(f1_score, 4),
AUC = round(auc_valor, 4)
)
)
return(resultados)
}set.seed(1234567)
n_prueba <- round(0.7*nrow(df_to_use))
id_muestra <- sample(1:nrow(df_to_use), n_prueba, replace = F)
datos_entren <- df_to_use[id_muestra, ]
datos_prueba <- df_to_use[-id_muestra, ]
m.final.logit <- glm(formula = Con_Embriaguez ~ Con_Huecos + Gravedad_indicador_30d +
Clase_Acc, family = binomial(link = "logit"), data = datos_entren)
y_logit.pred <- predict(m.final.logit, newdata = datos_prueba,type = "response")
y_observado <- datos_prueba$Con_EmbriaguezDiagnostico del modelo
diagnostico_logistico(y_observado, y_logit.pred, umbral=0.3)## Type 'citation("pROC")' for a citation.##
## Attaching package: 'pROC'## The following objects are masked from 'package:stats':
##
## cov, smooth, var## Setting levels: control = 0, case = 1## Setting direction: controls < cases## $Matriz_Confusion
## Predicho
## Real 0 1
## 0 3829 237
## 1 59 107
##
## $Metricas
## Exactitud Sensibilidad Especificidad Precision F1_Score AUC
## 1 0.9301 0.6446 0.9417 0.311 0.4196 0.8308Con_Huecos <- as.factor(c(0,0,1))
Gravedad_indicador_30d <- c('Con Muertos', 'Con Muertos', 'Solo Daños')
Clase_Acc <- c("Choque", "Choque", "Volcamiento")Predicciones.
nuevos_datos <- data.frame(
Con_Huecos = Con_Huecos,
Gravedad_indicador_30d = Gravedad_indicador_30d,
Clase_Acc = Clase_Acc
)
nuevos_datos$Probabilidad <- predict(m.final.logit,
newdata = nuevos_datos,
type = "response")
nuevos_datos$Prediccion_Clase <- ifelse(nuevos_datos$Probabilidad > 0.3,
"Si (Embriaguez)",
"No")
print(nuevos_datos)## Con_Huecos Gravedad_indicador_30d Clase_Acc Probabilidad Prediccion_Clase
## 1 0 Con Muertos Choque 0.0597447 No
## 2 0 Con Muertos Choque 0.0597447 No
## 3 1 Solo Daños Volcamiento 0.1069208 No