library(gridExtra)
library(dplyr)
library(readr)
library(ggplot2)
library(knitr)
library(leaflet)
library(stringr)
library(kableExtra)
library(patchwork)
library(DT)
library(ggmap)
library(EnvStats)
library(httr)
library(jsonlite)
library(sf)
library(naniar)
library(VIM)
library(tidyr)
library(summarytools)
library(lubridate)
library(mice)
library(skimr)Quinta entrega
Instrucciones
Realice un análisis descriptivo completo y actualizado al julio 31 2024, incluyendo graficos, tablas y georreferenciación de la base de datos Accidentalidad_en_Barranquilla.csv. Deben entregar un script .R con los códigos usados.
Contextualizar tanto la base de datos como las variables describiendo en qué consiste cada una de ellas. La base de datos está en el sitio web: https://www.datos.gov.co/Transporte/Accidentalidad-en-Barranquilla/yb9r-2dsi
Analizar las características de la base de datos. Estas pueden incluir: número de filas, número de columnas, nombres de las variables, tipos de variables, entre otros.
Analizar cada una de las variables según su tipo: numéricas y categóricas.
Filtrar la base de datos para entender mejor su estructura. Aplique filtros en al menos cinco oportunidades.
Utilice la función table para explorar la base de datos.
Identifique los valores NA (Not Available) en la base de datos.
Analice la presencia de posibles valores atípicos.
Aplique imputación de datos usando las siguientes opciones para method: method=“pmm” method=“norm.predict” method=“norm.nob” method=“norm”
Identifique datos atípicos para cada variable en el dataset usando las técnicas estudiadas en clase. Además, realice imputación de los datos atípicos con base en lo desarrollado en el ítem anterior. Imputación/Capping/Predicción Test de Rosner, Dixon, Grubbs, Hampel, Percentiles, Boxplots, Histogramas, Descriptivos.
Realice una exploración y un análisis descriptivo completo (incluyendo tablas de resumen y gráficos) de la base de datos disponible en este enlace. Incluya gráficos que presenten los datos sobre un mapa de Colombia, para visualizar la distribución geográfica de los precios de los combustibles en el país (2018-2024).
Librerías
Desarrollo de la primera parte
Sobre los datos
| FECHA_ACCIDENTE | HORA_ACCIDENTE | GRAVEDAD_ACCIDENTE | CLASE_ACCIDENTE | SITIO_EXACTO_ACCIDENTE | CANT_HERIDOS_EN._SITIO_ACCIDENTE | CANT_MUERTOS_EN._SITIO_ACCIDENTE | CANTIDAD_ACCIDENTES | AÑO_ACCIDENTE | MES_ACCIDENTE | DIA_ACCIDENTE |
|---|---|---|---|---|---|---|---|---|---|---|
| 2018-01-01T00:00:00.000 | 01:30:00:am | Con heridos | Atropello | CL 87 9H 24 | 1 | NA | 1 | 2018 | January | Mon |
| 2018-01-01T00:00:00.000 | 02:00:00:pm | Solo daños | Choque | CL 110 CR 46 | NA | NA | 1 | 2018 | January | Mon |
| 2018-01-01T00:00:00.000 | 04:00:00:am | Solo daños | Choque | AV CIRCUNVALAR CR 9G | NA | NA | 1 | 2018 | January | Mon |
| 2018-01-01T00:00:00.000 | 04:30:00:am | Solo daños | Choque | CLLE 72 CRA 29 | NA | NA | 1 | 2018 | January | Mon |
| 2018-01-01T00:00:00.000 | 05:20:00:pm | Solo daños | Choque | VIA 40 CALLE 75 | NA | NA | 1 | 2018 | January | Mon |
| 2018-01-01T00:00:00.000 | 06:00:00:pm | Con heridos | Choque | CR 8 CL 41 | 3 | NA | 1 | 2018 | January | Mon |
La base de datos proporcionada contiene 25,600 registros con información detallada sobre los accidentes de tránsito ocurridos en Barranquilla, según los informes policiales de accidentes de tránsito (IPAT). Para cada incidente, se reportan datos como la fecha, hora, gravedad, clase del accidente, sitio exacto, y cantidades de heridos y muertos en el lugar del accidente. La información es proporcionada por la Alcaldía Distrital de Barranquilla.
La base de datos cuenta con 11 variables tanto numéricas como categóricas. Las variables que contienen datos de texto incluyen la fecha completa del incidente, la hora exacta en la que ocurrió, la gravedad del accidente, tipo del incidente, y el lugar exacto donde sucedió. Por otro lado, las variables numéricas abarcan el número de heridos y fallecidos en el sitio del accidente, así como el total de accidentes registrados en ese momento. Además, componentes temporales como el año, mes y día de la semana también se almacenan como texto.
Enfocándonos un poco en las variables numéricas, veamos ahora la distribución de el número de heridos en sitio del accidente. En promedio se registraron, 1.47 heridos por accidente, también se destaca la presencia de 15626 datos faltantes.
summary(be5$CANT_HERIDOS_EN._SITIO_ACCIDENTE) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1.000 1.000 1.000 1.472 2.000 42.000 15626 ggplot(be5, aes(y = CANT_HERIDOS_EN._SITIO_ACCIDENTE)) +
geom_boxplot() +
labs(title = "Heridos en el Sitio del Accidente",
y = "Cantidad de Heridos",
x = "") +
theme_minimal()
Al graficar esta variable, observamos la presencia de muchos datos atípicos, con valores extremos hacia arriba. El registro con más heridos presentados fue con 42 heridos. Ahora, para la variable cantidad de muertos en accidente.
summary(be5$CANT_MUERTOS_EN._SITIO_ACCIDENTE) Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
1.000 1.000 1.000 1.036 1.000 2.000 25358 ggplot(be5, aes(y = CANT_MUERTOS_EN._SITIO_ACCIDENTE)) +
geom_boxplot() +
labs(title = "Muertos en el Sitio del Accidente",
y = "Cantidad de muertos",
x = "") +
theme_minimal()
Para aquellos registros donde se presentaron decesos, se tiene un promedio de una persona fallecida por accidente, sin embargo, existen outliers donde la cantidad de muertos asciende a dos personas. De igual manera que con la variable anterior, se destaca la presencia de 25358 datos faltantes.Y por último, para la variable relacionada al total de accidentes reportados en ese momento
summary(be5$CANTIDAD_ACCIDENTES) Min. 1st Qu. Median Mean 3rd Qu. Max.
1 1 1 1 1 2 ggplot(be5, aes(y = CANTIDAD_ACCIDENTES)) +
geom_boxplot() +
labs(title = "Accidentes en el mismo momento",
y = "Cantidad de accidentes",
x = "") +
theme_minimal()
De manera similar a la variable anterior, se tiene una media de 1 (sd: 0.013971) y mediana de 1 (RIQ: 0). También se destaca la presencia de outliers que indican 2 accidentes simultaneos, en el mismo sitio y a la misma hora. A diferencia de las otras dos variables numéricas, no se encuentran datos faltantes en esta variable.
Para las variables categóricas, comenzando con la fecha del accidente:
kable(head(sort(table(be5$FECHA_ACCIDENTE), decreasing = TRUE))) %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| Var1 | Freq |
|---|---|
| 2018-06-08T00:00:00.000 | 34 |
| 2018-11-16T00:00:00.000 | 33 |
| 2019-05-22T00:00:00.000 | 31 |
| 2019-10-05T00:00:00.000 | 30 |
| 2018-10-19T00:00:00.000 | 29 |
| 2019-11-20T00:00:00.000 | 29 |
Se tiene información de 2357 días, siendo el 8 de junio del 2018 el día con más accidentes (34 accidentes), seguido del 16 de noviembre del 2018 con 33 y el 22 de mayo del 2019 con 31. Realizando el mismo análisis con la hora de accidente.
kable(head(sort(table(be5$HORA_ACCIDENTE), decreasing = TRUE))) %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| Var1 | Freq |
|---|---|
| 03:00:00:pm | 408 |
| 04:00:00:pm | 408 |
| 12:30:00:pm | 387 |
| 05:00:00:pm | 384 |
| 08:00:00:am | 384 |
| 02:00:00:pm | 371 |
Se nota que las horas con más accidentes registrados son las 3pm (408), 4pm (408) y 12:30pm (387).
kable(table(be5$GRAVEDAD_ACCIDENTE))%>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| Var1 | Freq |
|---|---|
| Con heridos | 9901 |
| Con muertos | 252 |
| Solo daños | 15457 |
Afortunadamente, la mayoría de los accidentes fueron solo daños materiales, seguido de heridos y en menor medida con muertos (0.98%).
kable(table(be5$CLASE_ACCIDENTE))%>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| Var1 | Freq |
|---|---|
| Atropello | 1344 |
| Caida Ocupante | 194 |
| Choque | 23819 |
| Incendio | 13 |
| Otro | 123 |
| Volcamiento | 117 |
De los 25610 accidentes reportados, 23819 (93%) son debido a choques, 1344 (5,24%) por atropellos, el 1,76% restante está compuesto entre incendios, volcamientos, caída del ocupante, entre otros. Para la variable año, encontramos que 2018 es el año con más accidentes registrados y 2023 el que menos accidentes presentó, excluyendo 2024 dado que no se ha finalizado el año.
kable(table(be5$AÑO))%>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| Var1 | Freq |
|---|---|
| 2018 | 5898 |
| 2019 | 5645 |
| 2020 | 3281 |
| 2021 | 4700 |
| 2022 | 3683 |
| 2023 | 1662 |
| 2024 | 741 |
En cuanto a la variable relacionada al mes del accidente, se puede observar que febrero es el mes con más accidentes reportados en la ciudad de Barranquilla, seguido de marzo y enero. Agosto, con 1918 registros, es el mes con menos accidentes.
kable(table(be5$MES_ACCIDENTE))%>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| Var1 | Freq |
|---|---|
| April | 2010 |
| August | 1918 |
| December | 2189 |
| February | 2477 |
| January | 2349 |
| July | 1932 |
| June | 2103 |
| March | 2446 |
| May | 2121 |
| November | 1995 |
| October | 2090 |
| September | 1980 |
Y por último, el día con más accidentes reportados es el día de martes, seguido por viernes y miércoles. En cambio, los domingos, son los días con menos accidentes reportados.
kable(table(be5$DIA_ACCIDENTE))%>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| Var1 | Freq |
|---|---|
| Fri | 3920 |
| Mon | 3774 |
| Sat | 3735 |
| Sun | 2577 |
| Thu | 3756 |
| Tue | 4009 |
| Wed | 3839 |
Filtros
Veamos ahora aquellos accidentes realizados los fines de semana por choques según el mes.
be5 %>%
filter(DIA_ACCIDENTE %in% c('Sat', 'Sun')) %>%
filter(CLASE_ACCIDENTE == 'Choque') %>%
pull(MES_ACCIDENTE) %>%
table() %>%
kable() %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| . | Freq |
|---|---|
| April | 443 |
| August | 446 |
| December | 492 |
| February | 557 |
| January | 511 |
| July | 436 |
| June | 498 |
| March | 578 |
| May | 478 |
| November | 423 |
| October | 495 |
| September | 465 |
Siendo Marzo (578) el mes con más accidentes por choques en los fines de semana, seguido por febrero (557) y enero (511).Noviembre fue el mes que presentó menos accidentes los fines de semana (423).
be5 %>%
filter(!DIA_ACCIDENTE %in% c('Sat', 'Sun')) %>%
filter(CLASE_ACCIDENTE == 'Choque') %>%
pull(MES_ACCIDENTE) %>%
table() %>%
kable() %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| . | Freq |
|---|---|
| April | 1428 |
| August | 1354 |
| December | 1536 |
| February | 1753 |
| January | 1670 |
| July | 1395 |
| June | 1454 |
| March | 1690 |
| May | 1490 |
| November | 1426 |
| October | 1435 |
| September | 1366 |
Para el resto de días de la semana, febrero es el mes con más choques (1753), seguido de marzo (1690) y enero (1670). Agosto fue el mes con menos accidentes entre semana (1354).
be5 %>%
filter(MES_ACCIDENTE %in% c('January', 'February', 'March', 'April', 'May', 'June')) %>%
filter(GRAVEDAD_ACCIDENTE == 'Con heridos') %>%
pull(AÑO_ACCIDENTE) %>%
table()%>%
kable() %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| . | Freq |
|---|---|
| 2018 | 865 |
| 2019 | 795 |
| 2020 | 533 |
| 2021 | 659 |
| 2022 | 861 |
| 2023 | 833 |
| 2024 | 722 |
be5 %>%
filter(!MES_ACCIDENTE %in% c('January', 'February', 'March', 'April', 'May', 'June')) %>%
filter(GRAVEDAD_ACCIDENTE == 'Con heridos') %>%
pull(AÑO_ACCIDENTE) %>%
table()%>%
kable() %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| . | Freq |
|---|---|
| 2018 | 866 |
| 2019 | 795 |
| 2020 | 524 |
| 2021 | 782 |
| 2022 | 892 |
| 2023 | 774 |
Durante la primera mitad del año, 2020 fue el año con menos accidentes con heridos, mientras 2018 fue el año con más accidentes graves. Para el segundo semestre del año, el año con más accidentes graves fue 2022, nuevamente 2020 siendo el año con menos de este tipo de accidentes.
be5 %>%
filter(DIA_ACCIDENTE %in% c('Sat', 'Sun')) %>%
filter(GRAVEDAD_ACCIDENTE == 'Con heridos') %>%
pull(CANT_HERIDOS_EN._SITIO_ACCIDENTE) %>%
table()%>%
kable() %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| . | Freq |
|---|---|
| 1 | 1798 |
| 2 | 694 |
| 3 | 154 |
| 4 | 49 |
| 5 | 23 |
| 6 | 9 |
| 7 | 1 |
| 9 | 1 |
| 11 | 1 |
be5 %>%
filter(!DIA_ACCIDENTE %in% c('Sat', 'Sun')) %>%
filter(GRAVEDAD_ACCIDENTE == 'Con heridos') %>%
pull(CANT_HERIDOS_EN._SITIO_ACCIDENTE) %>%
table()%>%
kable() %>%
kable_styling(full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed"))| . | Freq |
|---|---|
| 1 | 5061 |
| 2 | 1588 |
| 3 | 309 |
| 4 | 103 |
| 5 | 34 |
| 6 | 20 |
| 7 | 13 |
| 8 | 8 |
| 9 | 5 |
| 10 | 7 |
| 11 | 5 |
| 12 | 5 |
| 13 | 2 |
| 14 | 2 |
| 15 | 1 |
| 16 | 1 |
| 18 | 1 |
| 19 | 1 |
| 20 | 1 |
| 21 | 1 |
| 22 | 1 |
| 23 | 1 |
| 42 | 1 |
Se puede observar como el rango de heridos en el sitio del accidente aumentó en los días de semana, manteniéndose entre 1 y 11 los sábados y domingo y entre 1 y 42 los días de semana. Aún así, son muchísimo más común que solo hayan entre 1 y 3 heridos.
Valores NA
summary(be5, na.rm = FALSE) FECHA_ACCIDENTE HORA_ACCIDENTE GRAVEDAD_ACCIDENTE CLASE_ACCIDENTE
Length:25610 Length:25610 Length:25610 Length:25610
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
SITIO_EXACTO_ACCIDENTE CANT_HERIDOS_EN._SITIO_ACCIDENTE
Length:25610 Min. : 1.000
Class :character 1st Qu.: 1.000
Mode :character Median : 1.000
Mean : 1.472
3rd Qu.: 2.000
Max. :42.000
NA's :15626
CANT_MUERTOS_EN._SITIO_ACCIDENTE CANTIDAD_ACCIDENTES AÑO_ACCIDENTE
Min. :1.000 Min. :1 Min. :2018
1st Qu.:1.000 1st Qu.:1 1st Qu.:2019
Median :1.000 Median :1 Median :2020
Mean :1.036 Mean :1 Mean :2020
3rd Qu.:1.000 3rd Qu.:1 3rd Qu.:2021
Max. :2.000 Max. :2 Max. :2024
NA's :25358
MES_ACCIDENTE DIA_ACCIDENTE
Length:25610 Length:25610
Class :character Class :character
Mode :character Mode :character
Se observa la presencia de NA’s solo en las columnas que brindan información sobre la cantidad de heridos y muertos en accidente, veamos ahora si el número de valores faltantes coincide con aquellos registros donde la gravedad del accidente sugiere la ausencia de personas heridas y muertas.
be5 %>%
# Heridos
filter(GRAVEDAD_ACCIDENTE == 'Con heridos') %>%
summary(na.rm = FALSE) FECHA_ACCIDENTE HORA_ACCIDENTE GRAVEDAD_ACCIDENTE CLASE_ACCIDENTE
Length:9901 Length:9901 Length:9901 Length:9901
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
SITIO_EXACTO_ACCIDENTE CANT_HERIDOS_EN._SITIO_ACCIDENTE
Length:9901 Min. : 1.000
Class :character 1st Qu.: 1.000
Mode :character Median : 1.000
Mean : 1.469
3rd Qu.: 2.000
Max. :42.000
CANT_MUERTOS_EN._SITIO_ACCIDENTE CANTIDAD_ACCIDENTES AÑO_ACCIDENTE
Min. : NA Min. :1 Min. :2018
1st Qu.: NA 1st Qu.:1 1st Qu.:2019
Median : NA Median :1 Median :2021
Mean :NaN Mean :1 Mean :2021
3rd Qu.: NA 3rd Qu.:1 3rd Qu.:2022
Max. : NA Max. :1 Max. :2024
NA's :9901
MES_ACCIDENTE DIA_ACCIDENTE
Length:9901 Length:9901
Class :character Class :character
Mode :character Mode :character
Tomando los accidentes con heridos, se puede observar que el número de NA’s en el número de muertos coincide con el número de registros de este dataframe filtrado, es decir, no hay problema con la presencia de estos datos dado que representan el que no hayan habido muertos en estos accidentes
be5 %>%
# Muertos
filter(GRAVEDAD_ACCIDENTE == 'Con muertos') %>%
summary(na.rm = FALSE) FECHA_ACCIDENTE HORA_ACCIDENTE GRAVEDAD_ACCIDENTE CLASE_ACCIDENTE
Length:252 Length:252 Length:252 Length:252
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
SITIO_EXACTO_ACCIDENTE CANT_HERIDOS_EN._SITIO_ACCIDENTE
Length:252 Min. : 1.000
Class :character 1st Qu.: 1.000
Mode :character Median : 1.000
Mean : 1.831
3rd Qu.: 2.000
Max. :12.000
NA's :169
CANT_MUERTOS_EN._SITIO_ACCIDENTE CANTIDAD_ACCIDENTES AÑO_ACCIDENTE
Min. :1.000 Min. :1 Min. :2018
1st Qu.:1.000 1st Qu.:1 1st Qu.:2019
Median :1.000 Median :1 Median :2021
Mean :1.036 Mean :1 Mean :2021
3rd Qu.:1.000 3rd Qu.:1 3rd Qu.:2023
Max. :2.000 Max. :1 Max. :2024
MES_ACCIDENTE DIA_ACCIDENTE
Length:252 Length:252
Class :character Class :character
Mode :character Mode :character
be5 %>%
# Muertos con heridos
filter(GRAVEDAD_ACCIDENTE == 'Con muertos') %>%
filter(!is.na(CANT_HERIDOS_EN._SITIO_ACCIDENTE)) %>%
summary(na.rm = FALSE) FECHA_ACCIDENTE HORA_ACCIDENTE GRAVEDAD_ACCIDENTE CLASE_ACCIDENTE
Length:83 Length:83 Length:83 Length:83
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
SITIO_EXACTO_ACCIDENTE CANT_HERIDOS_EN._SITIO_ACCIDENTE
Length:83 Min. : 1.000
Class :character 1st Qu.: 1.000
Mode :character Median : 1.000
Mean : 1.831
3rd Qu.: 2.000
Max. :12.000
CANT_MUERTOS_EN._SITIO_ACCIDENTE CANTIDAD_ACCIDENTES AÑO_ACCIDENTE
Min. :1.00 Min. :1 Min. :2018
1st Qu.:1.00 1st Qu.:1 1st Qu.:2019
Median :1.00 Median :1 Median :2021
Mean :1.06 Mean :1 Mean :2021
3rd Qu.:1.00 3rd Qu.:1 3rd Qu.:2022
Max. :2.00 Max. :1 Max. :2024
MES_ACCIDENTE DIA_ACCIDENTE
Length:83 Length:83
Class :character Class :character
Mode :character Mode :character
En los 7 años de registros del dataframe se presentaron 252 accidentes con al menos un deceso, donde 169 no presentaron heridos. Es decir, en los 83 accidentes restantes se presentó por lo menos una persona herida.
Para entender de mejor manera lo que está pasando con los datos faltantes, se hizo uso de la función aggr() de VIMla cual muestran la proporción y el patrón de estos datos.
aggr(be5,
numbers = TRUE,
col = c("navyblue", "red"),
cex.axis = 0.5,
gap = 3,
ylab = c("Proportion of missingness", "Missingness pattern"),
labels = c(
"Fecha",
"Hora",
"gravedad",
"clase",
"sitio",
"n_heridos",
"n_muertos",
"n_acc",
"Año",
"Mes",
"Día"
)
)
prop.table(table(be5$GRAVEDAD_ACCIDENTE))*100
Con heridos Con muertos Solo daños
38.6606794 0.9839906 60.3553299 Esto arroja un 60,36% registros con la ausencia de datos en en el número de heridos y número de muertos al tiempo, este porcentaje coincide con el 60,36% de accidentes relacionados a solo daño. De igual manera, corresponde también el 38,66% faltante de numero de muertos al 38,66% de aquellos accidentes con heridos. El porcentaje correspondiente al 0.9% de accidentes con muertes se reparte entre una ausencia de heridos de 0.66% y registros completos de 0.32%.
Esto quiere decir que los datos faltantes no requieren una imputación dado que no están relacionados a un error o datos no tomados, si no que hacen referencia a un valor que representa la ausencia de heridos o muertos como lo podría ser un 0.
be5s <- be5 %>%
mutate_all(~replace(., is.na(.), 0))
aggr(be5s,
numbers = TRUE,
col = c("navyblue", "red"),
cex.axis = 0.5,
gap = 3,
ylab = c("Proportion of missingness", "Missingness pattern"),
labels = c(
"Fecha",
"Hora",
"gravedad",
"clase",
"sitio",
"n_heridos",
"n_muertos",
"n_acc",
"Año",
"Mes",
"Día"
)
)
par(mfrow = c(1, 2), mar = c(5, 4, 4, 2) + 0.1)
hist(be5$CANT_HERIDOS_EN._SITIO_ACCIDENTE,
main = "BE5 (con NAs)",
xlab = "Cantidad de Heridos",
ylab = "Frecuencia",
col = "#b0ed9a",
border = "black",
xlim = c(0, 45),
breaks = 10)
hist(be5s$CANT_HERIDOS_EN._SITIO_ACCIDENTE,
main = "BE5S (con 0)",
xlab = "Cantidad de Heridos",
ylab = "Frecuencia",
col = "#b0ed9a",
border = "black",
xlim = c(0, 45),
breaks = 10)
par(mfrow = c(1, 1))La distribución no parece cambiar mucho, excepto por la gran concentración ahora de los datos en la parte izquierda de la gráfica de BE5S.
par(mfrow = c(1, 2), mar = c(5, 4, 4, 2) + 0.1)
hist(be5$CANT_MUERTOS_EN._SITIO_ACCIDENTE,
main = "BE5 (con NAs)",
xlab = "Cantidad de muertos",
ylab = "Frecuencia",
col = "#b0ed9a",
border = "black",
xlim = c(0, 3),
breaks =3)
hist(be5s$CANT_MUERTOS_EN._SITIO_ACCIDENTE,
main = "BE5S (con 0)",
xlab = "Cantidad de muertos",
ylab = "Frecuencia",
col = "#b0ed9a",
border = "black",
xlim = c(0, 3),
breaks = 3)
par(mfrow = c(1, 1))No se puede decir lo mismo de la cantidad de muertos, debido a que su distribución si cambia notoriamente. Esto era de esperarse debido al porcentaje de missingness que presenta esta variable (99,1%).
Debido al gran procentaje de datos faltantes también se podría proceder a eliminar las columnas involucradas, sin embargo para proceder con el análisis exploratorio se trabajará con la base de datos con los NAs. ### Gráficas
ggplot(be5, aes(x = AÑO_ACCIDENTE, fill = GRAVEDAD_ACCIDENTE)) +
geom_bar(position = "dodge") +
scale_x_continuous(breaks = seq(min(be5$AÑO_ACCIDENTE), max(be5$AÑO_ACCIDENTE), by = 1)) +
labs(title = "Relación entre Año del Accidente y Gravedad del Accidente",
x = "Año del Accidente",
y = "Número de Accidentes",
fill = "Gravedad del Accidente") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Se puede observar como 2018 fue un año con una recurrencia de accidentes bastante alta, por suerte, predominan aquellos accidentes donde solo hubo daños materiales en la mayoría de los años. En lo transcurrido del 2024, pasa algo muy curioso, y es que la categoría dominante de accidentes fueron aquellos con heridos, seguido de con muertos. Además, se notó una disminución bastante notable de accidentes en el 2020 en comparación al 2019 y en el 2021 en comparación al 2023. Los accidentes con muertes tienen recurrencias similares en todos los años.
be5 %>%
group_by(AÑO_ACCIDENTE, CLASE_ACCIDENTE) %>%
summarise(n = n(), .groups = 'drop') %>%
ggplot(aes(x = AÑO_ACCIDENTE, y = n, color = CLASE_ACCIDENTE, group = CLASE_ACCIDENTE)) +
geom_line(size = 0.5) +
scale_x_continuous(breaks = seq(min(be5$AÑO_ACCIDENTE), max(be5$AÑO_ACCIDENTE), by = 1)) +
labs(title = "Relación entre Año del Accidente y Tipo del Accidente",
x = "Año del Accidente",
y = "Número de Accidentes",
color = "Tipo de accidente") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Tal como se observó al momento de analizar cada variable individualmente, la proporción de accidentes por choques es muchísimo mayor en comparación a las otras razones. Se puede notar una disminución también los casos de accidentes por choque, mientras las otras categorías se mantienen constantes.
ggplot(be5, aes(x = CLASE_ACCIDENTE, y = CANT_HERIDOS_EN._SITIO_ACCIDENTE)) +
geom_boxplot() +
labs(title = "Distribución de la Cantidad de Heridos en el Sitio por Clase de Accidente",
x = "Clase de Accidente",
y = "Cantidad de Heridos en el Sitio del Accidente") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Aquí vemos por primera vez los atípicos, como se puede ver son muchísimos. Ya se había observado que la variable de la cantidad de heridos tenía un rango bastante amplio, desde, ahora por la imputación realizada, 0 hasta 42.
out <- boxplot.stats(be5$CANT_HERIDOS_EN._SITIO_ACCIDENTE)$out
out_ind <- which(be5$CANT_HERIDOS_EN._SITIO_ACCIDENTE %in% c(out))
datatable(be5s[out_ind, c(3,4,6)])Al tomar los datos atípicos, se observa la existencia de 303 datos atípicos en la variable CANT_HERIDOS_EN._SITIO_ACCIDENTE, tomando como atípicos los registros con más de 4 heridos en el accidente. Al realizar la misma observación en CANT_MUERTOS_EN._SITIO_ACCIDENTE se obtiene lo siguiente:
out <- boxplot.stats(be5$CANT_MUERTOS_EN._SITIO_ACCIDENTE)$out
out_ind <- which(be5$CANT_MUERTOS_EN._SITIO_ACCIDENTE %in% c(out))
datatable(be5s[out_ind, c(3,4,7)])Aquí se puede observar que tomó como atípicos todos aquellos accidentes donde hubo más de un fallecido. Es curioso como la imputación realizada en be5s, al momento de hallarle los valores atípicos, se obtienen muchísimos más debido a la concentración que está generando en el cero de estas variables. Incluso, toma como atípicos en el caso de CANT_MUERTOS_EN._SITIO_ACCIDENTE a todos aquellos valores difentes de cero, es decir, los que originalmente se encontraban en el dataframe.
Desarrollo de la segunda parte
La base de datos analizada proviene del Sistema de Información de la Cadena de Distribución de Combustibles (SICOM), un sistema gestionado por el Ministerio de Minas y Energía. Este sistema integra a los diferentes agentes en la cadena de distribución de combustibles a nivel nacional, facilitando la gestión y sistematización de datos relacionados con la comercialización, distribución, transporte y almacenamiento de combustibles.
Importar dataframes:
t32019 <- read.csv("C:/Users/Valen/Downloads/precios_t32019.csv")
t32018 <- read.csv("C:/Users/Valen/Downloads/precios_t32018.csv")
t32021 <- read.csv("C:/Users/Valen/Downloads/precios_t32021.csv")
t32020 <- read.csv("C:/Users/Valen/Downloads/precios_t32020.csv")
t32022 <- read.csv("C:/Users/Valen/Downloads/precios_t32022.csv")
t32023 <- read.csv("C:/Users/Valen/Downloads/precios_t32023.csv")
t42019 <- read.csv("C:/Users/Valen/Downloads/precios_t42019.csv")
t42018 <- read.csv("C:/Users/Valen/Downloads/precios_t42018.csv")
t42021 <- read.csv("C:/Users/Valen/Downloads/precios_t42021.csv")
t42020 <- read.csv("C:/Users/Valen/Downloads/precios_t42020.csv")
t42022 <- read.csv("C:/Users/Valen/Downloads/precios_t42022.csv")
t42023 <- read.csv("C:/Users/Valen/Downloads/precios_t42023.csv")
t12019 <- read.csv("C:/Users/Valen/Downloads/precios_t12019.csv")
t12018 <- read.csv("C:/Users/Valen/Downloads/precios_t12018.csv")
t12021 <- read.csv("C:/Users/Valen/Downloads/precios_t12021.csv")
t12020 <- read.csv("C:/Users/Valen/Downloads/precios_t12020.csv")
t12022 <- read.csv("C:/Users/Valen/Downloads/precios _t12022.csv")
t12023 <- read.csv("C:/Users/Valen/Downloads/precios_t12023.csv")
t12024 <- read.csv("C:/Users/Valen/Downloads/precios_t12024.csv")
t22019 <- read.csv("C:/Users/Valen/Downloads/precios_t22019.csv")
t22018 <- read.csv("C:/Users/Valen/Downloads/precios_t22018.csv")
t22021 <- read.csv("C:/Users/Valen/Downloads/precios_t22021.csv")
t22020 <- read.csv("C:/Users/Valen/Downloads/precios_t22020.csv")
t22022 <- read.csv("C:/Users/Valen/Downloads/precios_t22022.csv")
t22023 <- read.csv("C:/Users/Valen/Downloads/precios_t22023.csv")
t22024 <- read.csv("C:/Users/Valen/Downloads/precios_t22024.csv")Juntemos ahora los dataframes por años.
d2018 <- bind_rows(
t32018,
t42018,
t12018,
t22018
)
d2019 <- bind_rows(
t32019,
t42019,
t12019,
t22019
)
d2020 <- bind_rows(
t32020,
t42020,
t12020,
t22020
)
d2021 <- bind_rows(
t32021,
t42021,
t12021,
t22021
)
d2022 <- bind_rows(
t32022,
t42022,
t12022,
t22022
)
d2023 <- bind_rows(
t32023,
t42023,
t12023,
t22023
)
d2024 <- bind_rows(
t12024,
t22024
)
d2018 <- d2018 %>% mutate(AÑO = 2018)
d2019 <- d2019 %>% mutate(AÑO = 2019)
d2020 <- d2020 %>% mutate(AÑO = 2020)
d2021 <- d2021 %>% mutate(AÑO = 2021)
d2022 <- d2022 %>% mutate(AÑO = 2022)
d2023 <- d2023 %>% mutate(AÑO = 2023)
d2024 <- d2024 %>% mutate(AÑO = 2024)
df <- bind_rows(d2018, d2019, d2020, d2021, d2022, d2023, d2024)
head(df) BANDERA NOMBRE.COMERCIAL PRODUCTO
1 TERPEL ESTACION DE SERVICIO DISTRIBUIDORA EL PORVENIR DIESEL
2 TERPEL ESTACION DE SERVICIO DISTRIBUIDORA EL PORVENIR GASOLINA MOTOR
3 TERPEL SERVICENTRO TERPEL ABEJORRAL GASOLINA MOTOR
4 TERPEL SERVICENTRO TERPEL ABEJORRAL DIESEL
5 TERPEL ESTACION DE SERVICIO TERPEL LA PALOMA GASOLINA MOTOR
6 PRIMAX ESTACION DE SERVICIO PRIMAX AUTOPISTA GASOLINA MOTOR
FECHA.REGISTRO DEPARTAMENTO MUNICIPIO VALOR.PRECIO AÑO
1 01-Jul-2018 AMAZONAS LETICIA 8550 2018
2 01-Jul-2018 AMAZONAS LETICIA 11180 2018
3 01-Jul-2018 ANTIOQUIA ABEJORRAL 10150 2018
4 01-Jul-2018 ANTIOQUIA ABEJORRAL 9450 2018
5 01-Jul-2018 ANTIOQUIA BELLO 8910 2018
6 01-Jul-2018 ANTIOQUIA BELLO 9190 2018Sobre los datos
Debido a la gran cantidad de datos, se limitó a trabajar con el dataframe del 2023.
La base de datos incluye un total de 267,943 registros, cada uno representando registros diarios relacionados a la venta combustibles. Las variables incluyen información sobre el producto, su precio, y su ubicación y fecha de registro.
names(d2023)[1] "BANDERA" "NOMBRE.COMERCIAL" "PRODUCTO" "FECHA.REGISTRO"
[5] "DEPARTAMENTO" "MUNICIPIO" "VALOR.PRECIO" "AÑO" Esta consta de 9 variables, 8 categóricas y 1 numérica. La variable numérica corresponde a VALOR.PRECIO, una variable que contiene el precio del galón de gasolina el día de ese registro. `
summary(d2023$VALOR.PRECIO, na.rm = FALSE) Min. 1st Qu. Median Mean 3rd Qu. Max.
0 9350 10880 12023 13849 14750147 ggplot(d2023, aes(y = VALOR.PRECIO)) +
geom_boxplot() +
labs(title = "Precio del galón",
y = "Precio",
x = "") +
theme_minimal()
Se destaca la presencia de muchos atípicos, con valor extremadamente altos como 14750147, el cual no tiene sentido para el precio de un galón de gasolina. Más adelante, se verán más a detalle estos valores.
En cuanto a las variables categóricas restantes, tenemos `BANDERA que hace referencia a la compañía que tiene la sede especificada en el registro. Esta está liderada por Terpel (35,39%), Primax (13,47%), Biomax (12,02%), Texaco (10,11%) y Petromil (6,43%).
kable((sort(prop.table(table(d2023$BANDERA))*100, decreasing = TRUE))) %>%
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover"),
position = "center",
font_size = 12,
fixed_thead = TRUE
)| Var1 | Freq |
|---|---|
| TERPEL | 35.3892432 |
| PRIMAX | 13.4774933 |
| BIOMAX | 12.0268117 |
| TEXACO | 10.1111057 |
| PETROMIL | 6.4375632 |
| COOMULPINORT | 3.3283198 |
| AYATAWACOOP | 3.0069828 |
| ZEUSS | 2.7946242 |
| PETROBRAS | 1.8179240 |
| DISCOWACOOP | 1.4536674 |
| ECOS | 1.3484211 |
| ESSO | 1.3069944 |
| PETRDECOL | 1.3043819 |
| PETRODECOL | 1.2454141 |
| PUMA | 1.1416607 |
| DISCOM | 1.0886644 |
| OCTANO | 0.8289076 |
| PLUS MAS | 0.6740240 |
| P Y B | 0.6284919 |
| BRIO | 0.2683407 |
| ZAPATA Y VELASQUEZ | 0.1522712 |
| SAVE | 0.1063659 |
| PROXXON | 0.0623267 |
También se puede encontrar a la variable *NOMBRE.COMERCIAL* que indica el nombre de la estación de gasolina del registro. Esta cuenta con información sobre 5950 establecimientos en 911 municipios, guardados en la variable *MUNICIPIOS* de 32 departamentos, además de Bogotá D.C. como un departamento más representados en la variable `*DEPARTAMENTO*
length(unique(d2023$NOMBRE.COMERCIAL))[1] 5950length(unique(d2023$MUNICIPIO))[1] 911length(unique(d2023$DEPARTAMENTO))[1] 33kable((sort(prop.table(table(d2023$DEPARTAMENTO))*100, decreasing = TRUE))) %>%
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover"),
position = "center",
font_size = 12,
fixed_thead = TRUE
)| Var1 | Freq |
|---|---|
| NARIÑO | 11.5897784 |
| ANTIOQUIA | 9.3337016 |
| NORTE DE SANTANDER | 8.1181445 |
| VALLE DEL CAUCA | 6.9992498 |
| CUNDINAMARCA | 6.2069171 |
| BOGOTA D.C. | 5.9829889 |
| CESAR | 5.8605748 |
| LA GUAJIRA | 5.0891421 |
| SANTANDER | 3.8243208 |
| ATLANTICO | 3.2077718 |
| TOLIMA | 2.9909346 |
| BOYACA | 2.8994973 |
| CORDOBA | 2.7767100 |
| BOLIVAR | 2.6225727 |
| PUTUMAYO | 2.4669426 |
| HUILA | 2.3717731 |
| META | 2.3325857 |
| CAUCA | 2.0963414 |
| RISARALDA | 2.0284165 |
| MAGDALENA | 1.9925880 |
| CALDAS | 1.6858063 |
| SUCRE | 1.6264653 |
| CASANARE | 1.0240984 |
| CHOCO | 0.9852842 |
| QUINDIO | 0.9535610 |
| CAQUETA | 0.8919808 |
| ARAUCA | 0.8628701 |
| GUAVIARE | 0.3892619 |
| VICHADA | 0.3441030 |
| AMAZONAS | 0.1914586 |
| ARCHIPIELAGO DE SAN ANDRES, SANTA CATALINA Y PROVIDENCIA | 0.0996481 |
| GUAINIA | 0.0974088 |
| VAUPES | 0.0571017 |
Los departamentos con mayor número de registros son Nariño (11,58%), Antioquia con un 9,3% seguido de Norte de Santander (8,11%). Departamentos como Vaupés, Guianía y Archipielago de San Andrés y Providencia representan menos del 0.1% cada uno.
kable((sort(table(d2023$MES), decreasing = TRUE))) %>%
kable_styling(
full_width = FALSE,
bootstrap_options = c("striped", "hover"),
position = "center",
font_size = 12,
fixed_thead = TRUE
)| x |
|---|
| NA |
| --: |
El número de registros por mes no se mantuvo constante, siendo junio el mes con más registros (26887) y julio (18453) el mes con menos registros, más adelante veremos mejor el comportamiento de los datos según el mes.
d2023 <- d2023 %>%
mutate(MES = month(dmy(FECHA.REGISTRO), label = TRUE))
resumen_precios <- d2023 %>%
group_by(BANDERA, MES) %>%
summarise(PROMEDIO_PRECIO = mean(VALOR.PRECIO, na.rm = TRUE), .groups = 'drop')
ggplot(resumen_precios, aes(x = MES, y = PROMEDIO_PRECIO, color = BANDERA, group = BANDERA)) +
geom_line() +
geom_point() +
labs(title = "Precio Promedio por Mes y Bandera",
x = "Mes",
y = "Precio Promedio",
color = "Bandera") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Según un artículo publicado en septiembre del 2022 en portafolio.co, en Colombia, más del 90 % de la torta en la distribución de combustibles en el territorio nacional sigue concentrada en cinco grandes compañías: Terpel (42,6%), Primax (20,2%), Biomax (8,7%), Chevron (Texaco 10,7%) y Petromil (4,4%). Enfoquémonos en las principales 5 empresas.
d2023 <- d2023 %>%
mutate(BANDERA = str_squish(BANDERA))
resumen_precios <- d2023 %>%
filter(BANDERA %in% c("TERPEL", "PRIMAX", "BIOMAX", "TEXACO", "PETROMIL")) %>%
group_by(BANDERA, MES) %>%
summarise(PROMEDIO_PRECIO = mean(VALOR.PRECIO, na.rm = TRUE), .groups = 'drop')
ggplot(resumen_precios, aes(x = MES, y = PROMEDIO_PRECIO, color = BANDERA, group = BANDERA)) +
geom_line() +
geom_point() +
labs(title = "Precio Promedio por Mes y Bandera",
x = "Mes",
y = "Precio Promedio",
color = "Bandera") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Se puede notar como en el transcurso del año, se mantiene un precio promedio similar en todas las compañías gasolineras, sin embargo en septiembre hubo un disparo en el precio de BIOMAX en comparación a las demás empresas. También se puede notar que en general PRIMAX maneja un precio más alto, seguido de TEXACO, luego TERPEL y por últim BIOMAX.
d2023_mensual <- d2023 %>%
group_by(MES) %>%
summarise(PRECIO_PROMEDIO = mean(VALOR.PRECIO, na.rm = TRUE))
ggplot(d2023, aes(x = MES, y = VALOR.PRECIO)) +
geom_boxplot() +
labs(title = "Distribución del Precio por Mes",
x = "Mes",
y = "Precio") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
Se pueden observar muchísimos datos atípicos, pero el que más llama la atención es el que se encuentra en septiembre, el que se vio cuando se analizaron las variables individualmente como un valor muy exagerado. Veamos entonces todos los posibles valores atípicos que tiene nuestra base de datos.
out <- boxplot.stats(d2023$VALOR.PRECIO)$out
out_ind <- which(d2023$VALOR.PRECIO %in% c(out))
datatable(d2023[out_ind, ])Se puede observar una gran cantidad de datos atípicos (2907) tanto por arriba como por debajo, siendo 14750147 el precio que sobresalía tanto en la distribución, correspondiente a EDS DOÑA DANIELA en Bolívar. Lo más seguro es que la existencia de estos valores muy lejanos sean debido a errores al momento de digitar por tanto se realizará imputación de datos,
lb <- quantile(d2023$VALOR.PRECIO, 0.025)
up <- quantile(d2023$VALOR.PRECIO, 0.975)outlier_ind <- which(d2023$VALOR.PRECIO < lb | d2023$VALOR.PRECIO > up)
d2023$VALOR.PRECIO[outlier_ind] <- median(d2023$VALOR.PRECIO, na.rm = TRUE)
summary(d2023) BANDERA NOMBRE.COMERCIAL PRODUCTO FECHA.REGISTRO
Length:267943 Length:267943 Length:267943 Length:267943
Class :character Class :character Class :character Class :character
Mode :character Mode :character Mode :character Mode :character
DEPARTAMENTO MUNICIPIO VALOR.PRECIO AÑO
Length:267943 Length:267943 Min. : 7550 Min. :2023
Class :character Class :character 1st Qu.: 9410 1st Qu.:2023
Mode :character Mode :character Median :10880 Median :2023
Mean :11806 Mean :2023
3rd Qu.:13589 3rd Qu.:2023
Max. :20120 Max. :2023
MES
jun : 26887
nov : 25872
may : 25590
ene : 22369
feb : 22206
mar : 22125
(Other):122894 ggplot(d2023, aes(x = VALOR.PRECIO)) +
geom_histogram(binwidth = 1000, fill = "blue", color = "black") +
labs(title = "Distribución del Precio del Galón",
x = "Precio",
y = "Frecuencia") +
theme_minimal()
ggplot(d2023, aes(y = VALOR.PRECIO)) +
geom_boxplot() +
labs(title = "Precio del galón",
y = "Precio",
x = "") +
theme_minimal()
d2023 <- d2023 %>%
mutate(BANDERA = str_squish(BANDERA))
resumen_precios <- d2023 %>%
filter(BANDERA %in% c("TERPEL", "PRIMAX", "BIOMAX", "TEXACO", "PETROMIL")) %>%
group_by(BANDERA, MES) %>%
summarise(PROMEDIO_PRECIO = mean(VALOR.PRECIO, na.rm = TRUE), .groups = 'drop')
ggplot(resumen_precios, aes(x = MES, y = PROMEDIO_PRECIO, color = BANDERA, group = BANDERA)) +
geom_line() +
geom_point() +
labs(title = "Precio Promedio por Mes y Bandera",
x = "Mes",
y = "Precio Promedio",
color = "Bandera") +
theme_minimal() +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
En comparación a la gráfica realizada anteriormente con la base de datos sin imputar, podemos notar la gran diferencia con el pico que se apreciaba en el mes de septiembre puesto que ahí estaba el outlier más alejado por mucho.
Veamos ahora la distribución del precio promedio a nivel departamental.Por motivos estéticos se realizó con ggplot para evitar afectaciones en el formato del quarto. Se dejan los códigos correspondientes al mapa con leaf
Creación del mapa
api_key <- 'AIzaSyATpl4evKr2bgAOUhDrGz81NW7LeOfD9jY'
register_google(key = api_key)
mapa_base <- get_googlemap("Colombia",
maptype = "hybrid",
scale = 2,
zoom = 7)ℹ <https://maps.googleapis.com/maps/api/staticmap?center=Colombia&zoom=7&size=640x640&scale=2&maptype=hybrid&key=xxx>ℹ <https://maps.googleapis.com/maps/api/geocode/json?address=Colombia&key=xxx>mapa_col <- st_read("C:/Users/Valen/Downloads/COLOMBIA.shp", quiet = FALSE)Reading layer `COLOMBIA' from data source `C:\Users\Valen\Downloads\COLOMBIA.shp' using driver `ESRI Shapefile'
Simple feature collection with 33 features and 11 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: -81.73575 ymin: -4.227907 xmax: -66.84735 ymax: 13.39453
Geodetic CRS: WGS 84mapa_col <- st_make_valid(mapa_col) departamentos_shapefile <- unique(mapa_col$DPTO_CNMBR)
departamentos_dataframe <- unique(d2023$DEPARTAMENTO)
departamentos_no_en_dataframe <- setdiff(departamentos_shapefile, departamentos_dataframe)
departamentos_no_en_shapefile <- setdiff(departamentos_dataframe, departamentos_shapefile)mapa_col$DPTO_CNMBR <- str_replace_all(mapa_col$DPTO_CNMBR, "\\?", "Ñ")
mapa_col$DPTO_CNMBR <- str_replace_all(mapa_col$DPTO_CNMBR, "
ARCHIPIELAGO DE SAN ANDRES, SANTA CATALINA Y PROVIDENCIA", "ARCHIPIELAGO DE SAN ANDRES")
mapa_col <- mapa_col %>%
rename(DEPARTAMENTO = DPTO_CNMBR)precios_promedio <- d2023 %>%
group_by(DEPARTAMENTO) %>%
summarise(promedio = mean(VALOR.PRECIO, na.rm = TRUE))
datos_unidos <- merge(mapa_col, precios_promedio, by = "DEPARTAMENTO", all = TRUE, sort = FALSE)
datos_unidos_sf <- st_as_sf(datos_unidos)
datos_unidos_sf <- st_simplify(datos_unidos_sf, dTolerance = 0.001)ggplot(data = datos_unidos_sf) +
geom_sf(aes(fill = promedio), color = "white", size = 0.2) +
scale_fill_distiller(palette = "Spectral", name = "Precio Promedio", direction = 1) +
geom_text(aes(geometry = geometry, label = round(promedio, 2)), size = 2, color = "black", stat = "sf_coordinates") +
labs(title = "Precio promedio del galón por departamento") +
theme_minimal() +
theme(
plot.title = element_text(hjust = 0.5),
axis.title = element_blank(),
axis.text = element_blank(),
axis.ticks = element_blank()
)
# Error in `to_ring.default()`: ! Don't know how to get polygon data from object of class XY,GEOMETRYCOLLECTION,sfg
leaflet(datos_unidos_sf) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
fillColor = ~colorNumeric(palette = "Spectral", domain = datos_unidos_sf$promedio)(promedio),
weight = 1,
opacity = 1,
color = 'black',
fillOpacity = 0.7,
label = ~paste0(DEPARTAMENTO, ": ", round(promedio, 2)),
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "15px",
direction = "auto"
)
) %>%
addLegend(
pal = colorNumeric(palette = "Spectral", domain = datos_unidos_sf$promedio),
values = ~promedio,
opacity = 0.7,
title = "Precio Promedio",
position = "bottomright"
)