Ejercicio de Práctica # 5 - Parte 2

 

Introducción

Durante esta actividad, se trabajará con la base de datos precios_2023, la cual contiene información sobre el precio del combustible en Colombia durante 2023. Esta base de datos está disponible en la página oficial del Ministerio de Minas y Energía. Los datos incluyen registros sobre el precio del combustible en diversas estaciones de servicio del país, junto con la marca, el tipo de producto, la fecha de registro, el departamento, el municipio y el precio del combustible.

df <- read_excel("C:/Users/jcami/Downloads/precios_2023.xlsx")
df %<>% clean_names

El análisis de las dimensiones muestra que la base de datos está compuesta por 267,943 observaciones y siete variables:

df %>% glimpse

## Rows: 267,943
## Columns: 7
## $ bandera          <chr> "TERPEL", "TERPEL", "TERPEL", "TERPEL", "TERPEL", "TE…
## $ nombre_comercial <chr> "ESTACIONDESERVICIOFLUVIALBALSAJUNIN", "ESTACIONDESER…
## $ producto         <chr> "GASOLINAMOTOR", "GASOLINAMOTOR", "DIESEL", "DIESEL",…
## $ fecha_registro   <chr> "01-Apr-2023", "01-Apr-2023", "01-Apr-2023", "01-Apr-…
## $ departamento     <chr> "AMAZONAS", "AMAZONAS", "AMAZONAS", "AMAZONAS", "AMAZ…
## $ municipio        <chr> "LETICIA", "LETICIA", "LETICIA", "LETICIA", "LETICIA"…
## $ valor_precio     <dbl> 12968, 13000, 11800, 11980, 12605, 11800, 12968, 1180…

La tabla a continuación resume el tipo de cada variable:

Nombre <- c("bandera","nombre_comercial","producto","fecha_registro","departamento","municipio","valor_precio")
Tipo_original <- c("Cadena","Cadena","Cadena","Cadena","Cadena","Cadena","Numerico")
Tipo_modificado <- c("Factor","Factor","Factor","Fecha","Factor","Factor","Numerico")
as.data.frame(cbind(Nombre,Tipo_original,Tipo_modificado))

##             Nombre Tipo_original Tipo_modificado
## 1          bandera        Cadena          Factor
## 2 nombre_comercial        Cadena          Factor
## 3         producto        Cadena          Factor
## 4   fecha_registro        Cadena           Fecha
## 5     departamento        Cadena          Factor
## 6        municipio        Cadena          Factor
## 7     valor_precio      Numerico        Numerico

A continuación, se presentan las primeras 6 observaciones de la base de datos:

head(df) 

## # A tibble: 6 × 7
##   bandera nombre_comercial        producto fecha_registro departamento municipio
##   <chr>   <chr>                   <chr>    <chr>          <chr>        <chr>    
## 1 TERPEL  ESTACIONDESERVICIOFLUV… GASOLIN… 01-Apr-2023    AMAZONAS     LETICIA  
## 2 TERPEL  ESTACIONDESERVICIOBALS… GASOLIN… 01-Apr-2023    AMAZONAS     LETICIA  
## 3 TERPEL  ESTACIONDESERVICIOFLUV… DIESEL   01-Apr-2023    AMAZONAS     LETICIA  
## 4 TERPEL  ESTACIONDESERVICIOBALS… DIESEL   01-Apr-2023    AMAZONAS     LETICIA  
## 5 TERPEL  ESTACIONDESERVICIODIST… GASOLIN… 01-Apr-2023    AMAZONAS     LETICIA  
## 6 TERPEL  ESTACIONDESERVICIODIST… DIESEL   01-Apr-2023    AMAZONAS     LETICIA  
## # ℹ 1 more variable: valor_precio <dbl>

Ahora procedemos a ajustar el tipo de datos, a cambiar a minuscula el texto, y a reducir nombres para hacerlos más manejables en analisis posteriores:

selected_columns <- c("bandera","nombre_comercial","producto","departamento","municipio")
df[selected_columns] <- lapply(df[selected_columns], tolower)
df$bandera <- str_replace(df$bandera , "zapatayvelasquez", "zapata")
df[selected_columns] <- lapply(df[selected_columns], as.factor)
df$fecha_registro <- dmy(df$fecha_registro)

Finalmente concluimos esta primera sección, con una descripción general de la base de datos:

summary(df)

##          bandera                                 nombre_comercial 
##  terpel      :94823   edsautomotrizsanjosemina           :   546  
##  primax      :36112   estaciondeserviciojltsanalberto    :   459  
##  biomax      :32225   estaciondeserviciosanfrancisco     :   441  
##  texaco      :27092   estaciondeservicioautomotrizfullleo:   361  
##  petromil    :17249   estaciondeserviciootun             :   353  
##  coomulpinort: 8918   combustibleyhotellasmulas          :   344  
##  (Other)     :51524   (Other)                            :265439  
##           producto      fecha_registro                 departamento   
##  diesel       :108205   Min.   :2023-01-01   narino          : 31054  
##  extra        : 32400   1st Qu.:2023-04-01   antioquia       : 25009  
##  gasolinamotor:127338   Median :2023-06-07   nortedesantander: 21752  
##                         Mean   :2023-06-20   valledelcauca   : 18754  
##                         3rd Qu.:2023-09-22   cundinamarca    : 16631  
##                         Max.   :2023-12-31   bogotad.c.      : 16031  
##                                              (Other)         :138712  
##              municipio       valor_precio     
##  bogota,d.c.      : 16031   Min.   :       0  
##  cali             :  7388   1st Qu.:    9350  
##  sanandresdetumaco:  6540   Median :   10880  
##  sanjosedecucuta  :  6183   Mean   :   12023  
##  medellin         :  5429   3rd Qu.:   13849  
##  barranquilla     :  4366   Max.   :14750147  
##  (Other)          :222006

Estudio de las varibales

Bandera

La variable Bandera indica la marca del combustible. El análisis de esta variable iniciara construyendo una tabla de frecuencia que nos ayudara a determinar el número de observaciones en cada condición y la posible existencia de datos faltantes.

table(df$bandera)

## 
##  ayatawacoop       biomax         brio coomulpinort       discom  discowacoop 
##         8057        32225          719         8918         2917         3895 
##         ecos         esso       octano    petrdecol    petrobras   petrodecol 
##         3613         3502         2221         3495         4871         3337 
##     petromil      plusmas       primax      proxxon         puma          pyb 
##        17249         1806        36112          167         3059         1684 
##         save       terpel       texaco       zapata        zeuss 
##          285        94823        27092          408         7488

Debido a que la tabla de frencuencia muestra que no hay datos faltantes en esta variable, procederemos a realizar una tabla con un conteo y el porcentaje de observaciones según su clase:

Tabla_bandera <- df %>% 
  group_by(bandera) %>%
  count() %>% 
  ungroup() %>% 
  mutate(Proporcion = `n` / sum(`n`)) %>% 
  arrange(Proporcion) %>%
  mutate(Porcentaje = scales::percent(Proporcion))

Tabla_bandera

## # A tibble: 23 × 4
##    bandera        n Proporcion Porcentaje
##    <fct>      <int>      <dbl> <chr>     
##  1 proxxon      167   0.000623 0.0623%   
##  2 save         285   0.00106  0.1064%   
##  3 zapata       408   0.00152  0.1523%   
##  4 brio         719   0.00268  0.2683%   
##  5 pyb         1684   0.00628  0.6285%   
##  6 plusmas     1806   0.00674  0.6740%   
##  7 octano      2221   0.00829  0.8289%   
##  8 discom      2917   0.0109   1.0887%   
##  9 puma        3059   0.0114   1.1417%   
## 10 petrodecol  3337   0.0125   1.2454%   
## # ℹ 13 more rows

Culiminamos el analisis de variable Bandera visualizando el número de observaciones por bandera en un diagrama de barras.

ggplot(Tabla_bandera, aes(x=reorder(
bandera,-n),y=n))+
  geom_bar(stat="identity", position=position_dodge(),fill="darkolivegreen")+
  scale_y_continuous(labels = scales::label_number(scale = 1e-3, suffix = "K"))+
  theme(plot.title = element_text(hjust = 0.5, size = 12, vjust = 2),
    axis.text.x = element_text(angle = 90, hjust = 0.8, vjust = 0.5, ))+
  ggtitle("Observaciones por bandera de combustibles")+
  xlab("Bandera")+
  ylab("Observaciones")

Producto

La variable Producto indica el tipo de combustible. Tal como en la variable anterior, se construira una tabla de frecuencia para determinar el número de observaciones en cada condición y la posible existencia de datos faltantes.

table(df$producto)

## 
##        diesel         extra gasolinamotor 
##        108205         32400        127338

Debido a que la tabla de frencuencia muestra que no hay datos faltantes en esta variable, procederemos a realizar una tabla con un conteo y el porcentaje de observaciones según su clase:

Tabla_producto <- df %>% 
  group_by(producto) %>%
  count() %>% 
  ungroup() %>% 
  mutate(Proporcion = `n` / sum(`n`)) %>% 
  arrange(Proporcion) %>%
  mutate(Porcentaje = scales::percent(Proporcion))

Tabla_producto

## # A tibble: 3 × 4
##   producto           n Proporcion Porcentaje
##   <fct>          <int>      <dbl> <chr>     
## 1 extra          32400      0.121 12.1%     
## 2 diesel        108205      0.404 40.4%     
## 3 gasolinamotor 127338      0.475 47.5%

Culiminamos el analisis de variable Producto visualizando el número de observaciones por bandera en un diagrama de barras.

ggplot(Tabla_producto, aes(x=reorder(
producto,-n),y=n))+
geom_bar(width = 0.7, stat = "identity", position = position_dodge(0.7), fill = "cyan4") +  coord_cartesian(ylim = c(0, 180000)) + geom_text(aes(label = Porcentaje), vjust = -0.5, color = "black", hjust = 0.5,position = position_dodge(0.7),angle = 0,size = 4.5) +   scale_y_continuous(labels = scales::label_number(scale = 1e-3, suffix = "K"))+ theme(plot.title = element_text(hjust = 0.5, size = 12, vjust = 2), axis.text.x = element_text(angle = 0))+ ggtitle("Observaciones por tipo de producto")+ xlab("Producto")+ ylab("Observaciones")

Fecha

Tal como su nombre de indica, la variable fecha_registro ofrece información sobre la fecha en donde se recolecto la información. Debido a la naturaleza de esta varibale no es recomendable darle el mismo tratamiento aplicado en las secciones anteriores.En este caso, comenzaremos verificando el número de datos faltantes.

which(is.na(df$fecha_registro))

## integer(0)

Dado que no se han encontrado datos faltantes, se procedera visualizar el número de observaciones a lo largo del tiempo:

ggplot(df, aes(x = fecha_registro)) +
  geom_line(stat = "count", color = "blue") +
  labs(title = "Frecuencia de registros por fecha",
       x = "Fecha",
       y = "Numero de registros") +
  theme(
    plot.title = element_text(hjust = 0.5, vjust = 1, margin = margin(b = 20)),
    axis.title.x = element_text(margin = margin(t = 15)),
    axis.title.y = element_text(margin = margin(r = 15))
  )

De manera similar, podemos extraer los meses de la variable fecha_registro y estudiar el número de observaciones en cada uno de ellos. Esto se muestra acontinuación:

df$mes <- month(df$fecha_registro)
df$mes <- factor(df$mes,levels = 1:12,labels = c("Enero", "Febrero", "Marzo", "Abril", "Mayo", "Junio", "Julio", "Agosto", "Septiembre", "Octubre", "Noviembre", "Diciembre"))
table(df$mes)

## 
##      Enero    Febrero      Marzo      Abril       Mayo      Junio      Julio 
##      22369      22206      22125      21635      25590      26887      18453 
##     Agosto Septiembre    Octubre  Noviembre  Diciembre 
##      21959      21927      19487      25872      19433

Con estos datos podemos seguir el procedimiento aplicado en las primeras tres varibables y generar una tabla de frecuencia con datos descriptivos:

Tabla_mes <- df %>% 
  group_by(mes) %>%
  count() %>% 
  ungroup() %>% 
  mutate(Proporcion = `n` / sum(`n`)) %>% 
  arrange(Proporcion) %>%
  mutate(Porcentaje = scales::percent(Proporcion))

Tabla_mes

## # A tibble: 12 × 4
##    mes            n Proporcion Porcentaje
##    <fct>      <int>      <dbl> <chr>     
##  1 Julio      18453     0.0689 6.887%    
##  2 Diciembre  19433     0.0725 7.253%    
##  3 Octubre    19487     0.0727 7.273%    
##  4 Abril      21635     0.0807 8.074%    
##  5 Septiembre 21927     0.0818 8.183%    
##  6 Agosto     21959     0.0820 8.195%    
##  7 Marzo      22125     0.0826 8.257%    
##  8 Febrero    22206     0.0829 8.288%    
##  9 Enero      22369     0.0835 8.348%    
## 10 Mayo       25590     0.0955 9.551%    
## 11 Noviembre  25872     0.0966 9.656%    
## 12 Junio      26887     0.100  10.035%

ggplot(Tabla_mes, aes(x=mes,y=n))+
  geom_bar(stat="identity", position=position_dodge(),fill="cyan4")+
  coord_cartesian(ylim = c(0, 30000))+
  scale_y_continuous(labels = scales::label_number(scale = 1e-3, suffix = "K"))+
  theme(plot.title = element_text(hjust = 0.5, size = 12, vjust = 2),
    axis.text.x = element_text(angle = 90, hjust = 0.7, vjust = 0.5))+
  ggtitle("Observaciones por mes")+
  xlab("Mes")+
  ylab("Observaciones")

Municipio y Departamento

Las variables departamento y municipio contienen información sobre el lugar donde se registro la información. Debido a que estas variables son categoricas primer se examinaran sus tablas de frecuencia y se revisar la posibles presencias de datos faltantes.

table(df$departamento)

## 
##         amazonas        antioquia           arauca        atlantico 
##              513            25009             2312             8595 
##       bogotad.c.          bolivar           boyaca           caldas 
##            16031             7027             7769             4517 
##          caqueta         casanare            cauca            cesar 
##             2390             2744             5617            15703 
##            choco          cordoba     cundinamarca          guainia 
##             2640             7440            16631              261 
##         guaviare            huila        laguajira        magdalena 
##             1043             6355            13636             5339 
##             meta           narino nortedesantander         putumayo 
##             6250            31054            21752             6610 
##          quindio        risaralda san andres islas        santander 
##             2555             5435              267            10247 
##            sucre           tolima    valledelcauca           vaupes 
##             4358             8014            18754              153 
##          vichada 
##              922

head(table(df$municipio),30)

## 
##      abejorral         abrego       abriaqui        acacias        acevedo 
##             50           1061             33            459            178 
##           achi         agrado      aguachica        aguadas     aguadedios 
##             77             60           3660            134             24 
##        aguazul agustincodazzi           aipe          alban        albania 
##            280            188             99            151            667 
##         alcala         aldana     alejandria      algarrobo      algeciras 
##             90             97             28             35            131 
##       almaguer      alpujarra       altamira       alvarado          amaga 
##             79             46            208            116            124 
##         amalfi       ambalema       anapoima         ancuya      andalucia 
##            111             38             59             99            141

Debido a la gran número de categorias, la inspección de la tabla de frecuencia no es recomendable. Por lo cual, se solo nos hemos limitado a mostrar los 30 primeros municipios y procederemos a revisar los datos faltantes por otro metodo:

which(is.na(df$municipio))

## integer(0)

Tal que no se encontraon datos atipicos procederemos a presentar dos graficos de barras que resumen el numero de observaciones por condición:

Tabla_departamento <- df %>% 
  group_by(departamento) %>%
  count() %>% 
  ungroup() %>% 
  mutate(Proporcion = `n` / sum(`n`)) %>% 
  arrange(Proporcion) %>%
  mutate(Porcentaje = scales::percent(Proporcion))

Tabla_departamento

## # A tibble: 33 × 4
##    departamento         n Proporcion Porcentaje
##    <fct>            <int>      <dbl> <chr>     
##  1 vaupes             153   0.000571 0.0571%   
##  2 guainia            261   0.000974 0.0974%   
##  3 san andres islas   267   0.000996 0.0996%   
##  4 amazonas           513   0.00191  0.1915%   
##  5 vichada            922   0.00344  0.3441%   
##  6 guaviare          1043   0.00389  0.3893%   
##  7 arauca            2312   0.00863  0.8629%   
##  8 caqueta           2390   0.00892  0.8920%   
##  9 quindio           2555   0.00954  0.9536%   
## 10 choco             2640   0.00985  0.9853%   
## # ℹ 23 more rows

ggplot(Tabla_departamento, aes(x=reorder(
departamento,-n),y=n))+
  geom_bar(stat="identity", position=position_dodge(),fill="darkolivegreen")+
  scale_y_continuous(labels = scales::label_number(scale = 1e-3, suffix = "K"))+
  theme(plot.title = element_text(hjust = 0.5, size = 12, vjust = 2),
    axis.text.x = element_text(angle = 90, hjust = 0.8, vjust = 0.5, ))+
  ggtitle("Observaciones por departamento")+
  xlab("Departamento")+
  ylab("Observaciones")

Una vez más debido al gran numero de municipios considerados nos limitamos a graficas solo los 30 con más observaciones:

Tabla_municipio <- df %>% 
  group_by(municipio) %>%
  count() %>% 
  ungroup() %>% 
  mutate(Proporcion = `n` / sum(`n`)) %>% 
  arrange(Proporcion) %>%
  mutate(Porcentaje = scales::percent(Proporcion))

Tabla_municipio

## # A tibble: 910 × 4
##    municipio      n Proporcion Porcentaje
##    <fct>      <int>      <dbl> <chr>     
##  1 quipama        2 0.00000746 0.00075%  
##  2 mutiscua       4 0.0000149  0.00149%  
##  3 california     6 0.0000224  0.00224%  
##  4 guapota        8 0.0000299  0.00299%  
##  5 coper         10 0.0000373  0.00373%  
##  6 chima         12 0.0000448  0.00448%  
##  7 chinavita     12 0.0000448  0.00448%  
##  8 sanjose       16 0.0000597  0.00597%  
##  9 cacota        18 0.0000672  0.00672%  
## 10 ubaque        20 0.0000746  0.00746%  
## # ℹ 900 more rows

Tabla_municipio <- Tabla_municipio %>%
  arrange(desc(Proporcion))

ggplot(head(Tabla_municipio,30), aes(x=reorder(
municipio,-n),y=n))+
  geom_bar(stat="identity", position=position_dodge(),fill="cyan4")+
  scale_y_continuous(labels = scales::label_number(scale = 1e-3, suffix = "K"))+
  theme(plot.title = element_text(hjust = 0.5, size = 12, vjust = 2),
    axis.text.x = element_text(angle = 90, hjust = 0.8, vjust = 0.5, ))+
  ggtitle("Observaciones por municipio")+
  xlab("Municipio")+
  ylab("Observaciones")

Comercio

La variable nombre_comercial contiene información sobre las estaciones de servicios en donde se recolecto la información del valor de la gasolina. Esta varibale contiene un gran cantidad de niveles que no permiten un análisis directo de cada una de las categorias. Por lo cual, comenzaremos transformando la variable a factor y examinando si contiene datos faltantes:

df$nombre_comercial <- as.factor(df$nombre_comercial)
length(which(is.na(df$nombre_comercial)))

## [1] 0

Debido a que no se encontraron datos faltantes, se procedera a constrir una tabla de frecuencia con datos descriptivos:

Tabla_comercio <- df %>% 
  group_by(nombre_comercial) %>%
  count() %>% 
  ungroup() %>% 
  mutate(Proporcion = `n` / sum(`n`)) %>% 
  arrange(Proporcion) %>%
  mutate(Porcentaje = scales::percent(Proporcion))

Tabla_comercio

## # A tibble: 5,910 × 4
##    nombre_comercial                n Proporcion Porcentaje
##    <fct>                       <int>      <dbl> <chr>     
##  1 e.d.scombustibleslaperlita      1 0.00000373 0.00037%  
##  2 ancombustiblespacificrs         2 0.00000746 0.00075%  
##  3 conexioncuarentaycincosas       2 0.00000746 0.00075%  
##  4 e.d.s.playamar                  2 0.00000746 0.00075%  
##  5 e.d.saguaclara                  2 0.00000746 0.00075%  
##  6 edscombustiblesbecerras.a.s     2 0.00000746 0.00075%  
##  7 edseledencampestre              2 0.00000746 0.00075%  
##  8 edsnuevaesperanza               2 0.00000746 0.00075%  
##  9 estaciondegasolinaoronegro      2 0.00000746 0.00075%  
## 10 estaciondeservicio4m            2 0.00000746 0.00075%  
## # ℹ 5,900 more rows

Una vez más debido al gran numero de categoria nos limitamos a graficas solo los 10 con más observaciones:

Tabla_comercio <- Tabla_comercio %>%
  arrange(desc(Proporcion))

ggplot(head(Tabla_comercio, 10), aes(x = reorder(nombre_comercial, -n), y = n)) +
  geom_bar(stat = "identity", position = position_dodge(), fill = "cyan4") +
  coord_flip() +
  theme(
    plot.title = element_text(hjust = 0.5),
    axis.title.x = element_text(margin = margin(t = 15)),  
    axis.title.y = element_text(margin = margin(r = 15)),
    axis.text.y = element_text(angle = 0, hjust = 0.5),  
    plot.margin = margin(10, 10, 10, 20)
  ) +
  ggtitle("Observaciones por Comercio") +
  xlab("Comercio") +
  ylab("Observaciones")

Precio

Finalmente, se examinara la variable valor_precio la cual contiene el valor de la gasolina en cada una de las estaciones del país. Primeramente comenzaremos graficando un histograma y un diagramas de caja que nos ayudara estudiar la distribución de la variable y posible datos atipicos:

A <- ggplot(df, aes(x = valor_precio)) +
  geom_histogram(binwidth = 3600, fill = "cyan4", color = "black") +
  labs(title = "Histograma", x = "Valor de la gasolina", y = "Frecuencia") + 
  theme_minimal() +
    theme(
    plot.title = element_text(hjust = 0.5),  # Centrar el título
    axis.title.x = element_text(margin = margin(t = 15)),  
    axis.title.y = element_text(margin = margin(r = 15)),
    plot.margin = margin(10, 20, 10, 10) 
  )

B <- ggplot(df, aes(x = "", y = valor_precio)) +
  geom_boxplot(fill = "cyan4", color = "black") +
  labs(title = "Boxplot", x = "", y = "Valor de la gasolina") +
  theme_minimal() +
    theme(
    plot.title = element_text(hjust = 0.5),  # Centrar el título
    axis.title.x = element_text(margin = margin(t = 15)),  
    axis.title.y = element_text(margin = margin(r = 15)),
    plot.margin = margin(10, 10, 10, 20)
  )

grid.arrange(A, B, ncol = 2)

Los graficos revelan la existencia de valores atipicos en la variable valor_precio . Asi mismo, se encuentran valores no aceptables en la varibale (p.e valor de gasolina 0). Por lo cual, antes de proseguir utilizaremos el filtro de Hampel y el método de los percentiles para estudiar más a profundidad estas observaciones.

# Calcular los límites usando el filtro de Hampel
lower_bound <- median(df$valor_precio, na.rm = TRUE) - 3 * mad(df$valor_precio, constant = 1, na.rm = TRUE)
upper_bound <- median(df$valor_precio, na.rm = TRUE) + 3 * mad(df$valor_precio, constant = 1, na.rm = TRUE)

# Identificar los índices de los outliers
outlier_ind <- which(df$valor_precio < lower_bound | df$valor_precio > upper_bound)

# Crear una columna que marque los outliers
df$Hampel <- ifelse(1:nrow(df) %in% outlier_ind, 1, 0)

# Ver el número de outliers
length(outlier_ind)

## [1] 30771

lower_bound <- quantile(df$valor_precio, 0.025)
upper_bound <- quantile(df$valor_precio, 0.975)
outlier_ind <- which(df$valor_precio < lower_bound | df$valor_precio > upper_bound)

# Crear una columna que marque los outliers
df$Perc <- ifelse(1:nrow(df) %in% outlier_ind, 1, 0)

# Ver el número de outliers
length(outlier_ind)

## [1] 13386

table(df$Hampel, df$Perc)

##    
##          0      1
##   0 230611   6561
##   1  23946   6825

Ahora se examinaran el número de datos atipicos en función a la diferentes categorias presentes en la base de datos:

df %>%
  group_by(bandera) %>%
  summarise(
    total = n(),
    Casos_Hampel = sum(Hampel),
    Casos_perc = sum(Perc),
    Procentaje_Hampel = round((sum(Hampel)/30771) * 100,3 ),
    Procentaje_perc = round((sum(Perc)/13386) * 100,3 ))

## # A tibble: 23 × 6
##    bandera      total Casos_Hampel Casos_perc Procentaje_Hampel Procentaje_perc
##    <fct>        <int>        <dbl>      <dbl>             <dbl>           <dbl>
##  1 ayatawacoop   8057           44        524             0.143           3.92 
##  2 biomax       32225         2052        809             6.67            6.04 
##  3 brio           719           83         34             0.27            0.254
##  4 coomulpinort  8918          120       1958             0.39           14.6  
##  5 discom        2917          106          8             0.344           0.06 
##  6 discowacoop   3895            0       1724             0              12.9  
##  7 ecos          3613          283          4             0.92            0.03 
##  8 esso          3502          622        294             2.02            2.20 
##  9 octano        2221          251          0             0.816           0    
## 10 petrdecol     3495            2         13             0.006           0.097
## # ℹ 13 more rows

df %>%
  group_by(producto) %>%
  summarise(
    total = n(),
    Casos_Hampel = sum(Hampel),
    Casos_perc = sum(Perc),
    Procentaje_Hampel = round((sum(Hampel)/30771) * 100,3 ),
    Procentaje_perc = round((sum(Perc)/13386) * 100,3 ))

## # A tibble: 3 × 6
##   producto       total Casos_Hampel Casos_perc Procentaje_Hampel Procentaje_perc
##   <fct>          <int>        <dbl>      <dbl>             <dbl>           <dbl>
## 1 diesel        108205          154       6519              0.5            48.7 
## 2 extra          32400        29714       6618             96.6            49.4 
## 3 gasolinamotor 127338          903        249              2.94            1.86

Notamos que el procentaje de datos atipicos depende del tipo de producto, espacialmente al utlizar el filtro de Hampel. Debido a esto procederemos a construir un diagrama de cajas por grupos del valor de la gasolina:

ggplot(df[df$valor_precio < 300000
,], aes(x = "", y = valor_precio, group = producto)) +
  geom_boxplot(fill = "cyan4", color = "black") +
  labs(title = "Boxplot", x = "", y = "Valor de la gasolina") +
  theme_minimal() +
    theme(
    plot.title = element_text(hjust = 0.5),  # Centrar el título
    axis.title.x = element_text(margin = margin(t = 15)),  
    axis.title.y = element_text(margin = margin(r = 15)),
    plot.margin = margin(10, 10, 10, 20)
  )

Los datos muestran que aun cuando se considera el tipo de producto, los valores mayores a 50.000 se desvian considerablemente de las distribuciones. Para confirmar estas observaciones se aplicara la prueba de Rosner:

test <- rosnerTest(df$valor_precio, k = 100)
test <- test$all.stats
outlier_ind <- test[test$Outlier == TRUE,]$Obs.Num
df$Rosner <- ifelse(1:nrow(df) %in% outlier_ind, 1, 0)

Los resultados del test de Rosener indica que los valores mayores a 30.000 como datos atipicos, pero no se marcaron los valores en cero. En vista de esto, se procedera a reemplazar los valores atipicos y aquellos los menores a 1000 como faltantes.

df$valor_precio <- ifelse(df$valor_precio < 1000 | df$Rosner == TRUE, NA, df$valor_precio)

Ahora volveremos a graficar el valor del combustible:

ggplot(df, aes(x = "", y = valor_precio, group = producto)) +
  geom_boxplot(fill = "cyan4", color = "black") +
  labs(title = "Boxplot", x = "", y = "Valor de la gasolina") +
  theme_minimal() +
    theme(
    plot.title = element_text(hjust = 0.5),  # Centrar el título
    axis.title.x = element_text(margin = margin(t = 15)),  
    axis.title.y = element_text(margin = margin(r = 15)),
    plot.margin = margin(10, 10, 10, 20)
  )

Finalmente, calcularemos estadisticas descriptivas de la varibale valor_precio:

df %>%
  get_summary_stats(valor_precio, type = "five_number")

## # A tibble: 1 × 7
##   variable          n   min   max    q1 median    q3
##   <fct>         <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl>
## 1 valor_precio 267843  1000 30000  9350  10880 13849

Imputación

Ahora procederemos a tratar con los datos faltantes. Primero calcularemos el procentaje de datos faltantes en la base de datos.

aggr(df, numbers = TRUE, col = c("navyblue", "red"), 
     cex.axis = 0.7, gap = 3, ylab = c("Proportion of missingness", "Missingness pattern"))

(length(which(is.na(df$valor_precio)))/nrow(df)*100)

## [1] 0.03732137

El analisis muestra que la unica varibale que contiene datos faltantes es el valor del combustible, siendo aproximadamente un 0.04% del total de observaciones. Para esto utilizaremos la libreria mice:

imp <- mice(df, m=5, maxit=50, method ='pmm', seed=500, printFlag = F)

Los resultados indican que no hay espacio disponible para realizar la imputación, presumiblemente debido al tamaño de la base de datos. En vista de esto, utilizaremos la mediana para imputar los datos de la variable. Este método se seleccionó por dos razones: 1) la inspección de los datos indica que no siguen una distribución normal, y 2) la baja proporción de datos faltantes sugiere que la pérdida de varianza asociada con la imputación simple puede ser aceptable. Sin embargo, se reconocen las ventajas de la imputación múltiple, y se sugiere realizar este ejercicio cuando se disponga de mayores recursos computacionales.

df_antes_imput <- df$valor_precio
df[which(is.na(df$valor_precio)),]$valor_precio <- median(df$valor_precio, na.rm = T)

Para comprobar la efectividad de la imputación, comprobaremos si existe algún cambio en la distribución de la varibale antes y despues de la imputción:

df_A <- as.data.frame(cbind(df_antes_imput,rep("antes",length(df_antes_imput))))
colnames(df_A) <- c("valor","imput")
df_B <- as.data.frame(cbind(df$valor_precio,rep("despues",length(df_antes_imput))))
colnames(df_B) <- c("valor","imput")

df_imput <- rbind(df_A,df_B)

levene_test(df_imput, as.numeric(valor) ~ as.factor(imput)) # Los datos son homocedasticos

## # A tibble: 1 × 4
##     df1    df2 statistic     p
##   <int>  <int>     <dbl> <dbl>
## 1     1 535784    0.0219 0.882

wilcox.test(df_antes_imput,df$valor_precio)

## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  df_antes_imput and df$valor_precio
## W = 3.5883e+10, p-value = 0.9999
## alternative hypothesis: true location shift is not equal to 0

Los prueba de wilcoxon revela que la distribución de la variable no cambio despues de la imputación. Por lo cual se concluye que el proceso se llevo a cabo adecuadamente. Concluiremos esta sección graficando la variable valor_precio una ultima vez:

ggplot(df, aes(x = "", y = valor_precio)) +
  geom_boxplot(fill = "cyan4", color = "black") +
  labs(title = "Boxplot", x = "", y = "Valor de la gasolina") +
  theme_minimal() +
    theme(
    plot.title = element_text(hjust = 0.5),  # Centrar el título
    axis.title.x = element_text(margin = margin(t = 15)),  
    axis.title.y = element_text(margin = margin(r = 15)),
    plot.margin = margin(10, 10, 10, 20)
  )

Análisis Adicional

Antes de proceder a visualizar los precios de la gasolina en el mapa de Colombia, realizaremos algunos análisis y visualizaciones adicionales para obtener una comprensión más profunda de los datos. Primero, examinaremos la distribución de los precios del combustible por tipo utilizando un diagrama de caja.

Distribución de Precios por Tipo de Combustible

El diagrama de caja muestra la distribución de precios para diferentes tipos de combustibles. Este gráfico es útil para identificar la variabilidad y los valores atípicos en los precios de cada tipo de combustible.

ggplot(df, aes(x = producto, y = valor_precio, fill = producto)) + 
  geom_boxplot() +
  labs(title = "Distribucion de Precios por Tipo de Combustible",
       x = "Producto",
       y = "Precio") +
  theme_minimal() +
  theme(legend.position = "none",
        axis.text.x = element_text(angle = 45, hjust = 1),
        plot.title = element_text(hjust = 0.5))

  scale_y_continuous(labels = scales::comma)

## <ScaleContinuousPosition>
##  Range:  
##  Limits:    0 --    1

Este gráfico sugiere que la gasolina extra tiende a tener un valor promedio más alto en comparación con la gasolina motor y el diésel. Esto se alinea con las tendencias generales del mercado de los últimos años.

Precio Promedio del Diésel por Departamento

A continuación, visualizaremos el precio promedio del diésel en cada departamento del país. Esta visualización permitirá identificar las regiones donde el precio del diésel es más alto.

resum <- df[df$producto == "diesel",] %>% group_by(departamento) %>%
  summarise(precio=mean(valor_precio))

ggplot(resum, aes(x=reorder(
departamento,-precio),y=precio))+
  geom_bar(stat="identity", position=position_dodge(),fill="darkolivegreen")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ggtitle("Precio promedio gasolina diesel segun departamento")+
  xlab("Departamentos")+
  ylab("precio promedio")

Se observa que los departamentos de Vaupés, Amazonas y Guainía presentan los precios más altos para el diésel.

Precio Promedio de la Gasolina Extra por Departamento

Ahora procederemos a visualizar el precio promedio de la gasolina extra en los diferentes departamentos del país.

resum <- df[df$producto == "extra",] %>% group_by(departamento) %>%
  summarise(precio=mean(valor_precio))

ggplot(resum, aes(x=reorder(
departamento,-precio),y=precio))+
  coord_cartesian(ylim = c(0, 25000)) +
  geom_bar(stat="identity", position=position_dodge(),fill="cyan4")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ggtitle("Precio promedio gasolina extra segun departamento")+
  xlab("Departamentos")+
  ylab("precio promedio")

Precio Promedio de la Gasolina Motor por Departamento

Finalmente, visualizaremos el precio promedio de la gasolina motor en cada departamento.

resum <- df[df$producto == "gasolinamotor",] %>% group_by(departamento) %>%
  summarise(precio=mean(valor_precio))

ggplot(resum, aes(x=reorder(
departamento,-precio),y=precio))+
  geom_bar(stat="identity", position=position_dodge(),fill="red")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ggtitle("Precio promedio gasolina motor segun departamento")+
  xlab("Departamentos")+
  ylab("precio promedio")

Precio Promedio de la Gasolina por Marca (o Bandera)

Por último, examinaremos el precio promedio de cada tipo de gasolina según la marca o bandera. Esto nos permitirá comparar cómo varían los precios entre diferentes marcas.

resum <- df[df$producto == "gasolinamotor",] %>% group_by(bandera) %>%
  summarise(precio=mean(valor_precio))

ggplot(resum, aes(x=reorder(
bandera,-precio),y=precio))+
  coord_cartesian(ylim = c(0, 20000)) +
  geom_bar(stat="identity", position=position_dodge(),fill="red")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ggtitle("Precio promedio gasolina motor segun bandera")+
  xlab("Departamentos")+
  ylab("precio promedio")

resum <- df[df$producto == "extra",] %>% group_by(bandera) %>%
  summarise(precio=mean(valor_precio))

ggplot(resum, aes(x=reorder(
bandera,-precio),y=precio))+
  coord_cartesian(ylim = c(0, 20000)) +
  geom_bar(stat="identity", position=position_dodge(),fill="cyan4")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ggtitle("Precio promedio gasolina motor segun bandera")+
  xlab("Departamentos")+
  ylab("precio promedio")

resum <- df[df$producto == "diesel",] %>% group_by(bandera) %>%
  summarise(precio=mean(valor_precio))

ggplot(resum, aes(x=reorder(
bandera,-precio),y=precio))+
  coord_cartesian(ylim = c(0, 12000)) +
  geom_bar(stat="identity", position=position_dodge(),fill="darkolivegreen")+
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  ggtitle("Precio promedio gasolina motor segun bandera")+
  xlab("Departamentos")+
  ylab("precio promedio")

Mapa

En la ultima sección de este ejercicio se graficara el precio de la gasolina motor en cada departamente de colombia utilizando un mapa de Colombia.

mapa_col <- st_read("C:/Users/jcami/OneDrive/Escritorio/DATAVIZ/Departamentos_Junio_2024_shp/Departamentos_Junio_2024_shp/Departamento.shp", quiet = FALSE)

## Reading layer `Departamento' from data source 
##   `C:\Users\jcami\OneDrive\Escritorio\DATAVIZ\Departamentos_Junio_2024_shp\Departamentos_Junio_2024_shp\Departamento.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 33 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 4047822 ymin: 1090467 xmax: 5684465 ymax: 3053707
## Projected CRS: MAGNA-SIRGAS 2018 / Origen-Nacional

mapa_col <- st_make_valid(mapa_col)

Ahora se extraen los nombres de los departamentos desde la base de datos y el archivo shapefile. Los nombres se convierten a minúsculas y se eliminan caracteres especiales para facilitar la comparación. Esto asegura que los nombres de los departamentos en ambas fuentes sean comparables.

# Extraer los nombres de los departamentos desde la base de datos y el shapefile
dpto_1 <- tolower(df$departamento)  # Base de datos
dpto_2 <- tolower(mapa_col$DeNombre)# Archivo shapefile
dpto_2 <- iconv(dpto_2, from = "UTF-8", to = "ASCII//TRANSLIT")
dpto_1 <- iconv(dpto_1, from = "UTF-8", to = "ASCII//TRANSLIT")
dpto_2 <- as.character(dpto_2)
dpto_1 <- as.character(dpto_1)

Se verifica qué nombres de la base de datos no están presentes en el shapefile. Se identifican los nombres de los departamentos que no coinciden y se corrigen manualmente para que coincidan con los del shapefile. Esto incluye ajustes como renombrar “bogotad.c.” a “cundinamarca” y “valledelcauca” a “valle del cauca”.

# Verificar qué nombres de la base de datos no están en el shapefile
M1 <- which(is.na(match(dpto_1, dpto_2)))
na_df <- df[M1,]
na_df$departamento <- droplevels(na_df$departamento)
table(na_df$departamento)

## 
##       bogotad.c.        laguajira nortedesantander san andres islas 
##            16031            13636            21752              267 
##    valledelcauca 
##            18754

dpto_1 <- ifelse(dpto_1 == "bogotad.c.", "cundinamarca", dpto_1)
dpto_1 <- ifelse(dpto_1 == "nortedesantander", "norte de santander", dpto_1)
dpto_1 <- ifelse(dpto_1 == "valledelcauca", "valle del cauca", dpto_1)
dpto_1 <- ifelse(dpto_1 == "laguajira", "la guajira", dpto_1)
dpto_2 <- ifelse(dpto_2 == "san andres providencia y santa catalina", "san andres islas", dpto_2)

Verificamos que los nombres la base de datos y el archivos del mapa concuerden:

M2 <- which(is.na(match(dpto_1, dpto_2)))
na_df <- df[M2,]
na_df$departamento <- droplevels(na_df$departamento)
table(na_df$departamento)

## < table of extent 0 >

Se asignan códigos a los departamentos tanto en la base de datos como en el shapefile para facilitar la combinación de datos. Luego, se fusionan los datos de la base de datos con el shapefile utilizando la función merge(), asegurando que se conserven todos los datos del shapefile.

df$codigo <- dpto_1
mapa_col$codigo <- dpto_2
datos_unidos <- merge(st_drop_geometry(mapa_col), df, by = "codigo", all = TRUE, sort = FALSE)

Ahora se filtran los datos para obtener solo la información relevante sobre el producto “gasolinamotor”. Los datos se agrupan por departamento y se calcula el precio promedio de la gasolina motor en cada uno.

datos_unidos2 <- datos_unidos %>%
  filter(producto == "gasolinamotor") %>%
  group_by(codigo, DeCodigo, DeNombre, DeArea, DeNorma, SHAPE_Leng, SHAPE_Area) %>%
  summarize(precio_prom_dpto = mean(valor_precio)) %>%
  ungroup() 

## `summarise()` has grouped output by 'codigo', 'DeCodigo', 'DeNombre', 'DeArea',
## 'DeNorma', 'SHAPE_Leng'. You can override using the `.groups` argument.

Se genera un resumen estadístico y se visualiza la distribución de los precios utilizando gráficos de caja y histogramas.

summary(datos_unidos2$precio_prom_dpto)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   10789   12628   12923   13160   13240   18459

breaks <- c(0, 12000, 13000, 14000, 16000, 20000)

# Crear el histograma
ggplot(datos_unidos2, aes(x = precio_prom_dpto)) +
  geom_histogram(bins = length(breaks), fill = "cyan4", color = "black") +
  scale_x_continuous(breaks = breaks) +
  labs(title = "Histograma de precios de gasolina",
       x = "Valor de la gasolina",
       y = "Frecuencia") + 
  theme_minimal() +
  theme(
    plot.title = element_text(hjust = 0.5),  # Centrar el título
    axis.title.x = element_text(margin = margin(t = 15)),  
    axis.title.y = element_text(margin = margin(r = 15)),
    plot.margin = margin(10, 20, 10, 10) 
  )

Se categorizan los precios promedio de la gasolina en rangos y se asignan colores a cada rango. Esto facilita la visualización de las diferencias en los precios de gasolina entre departamentos en el mapa.

datos_unidos2$precio_prom_dpto <- as.numeric(datos_unidos2$precio_prom_dpto)
precio_fac <- cut(datos_unidos2$precio_prom_dpto, 
                  breaks = c(0, 12000,13000, 14000, 16000, 20000),
                  labels = c("a. < 12.000", "b. 12.000-13.000", "c. 13.000-14.000", "d. 14.000-16.000", "e. + 16.000"))

colores <- c("a. < 12.000" = "#66FFFF",
             "b. 10.000-12.000" = "green",
             "b. 12.000-14.000" = "yellow",
             "c. 14.000-16.000" = "orange",
             "d. +16.000" = "red")

Se asignan los colores a los datos de acuerdo con las categorías de precios. Se verifica que la longitud de los datos de precios coincida con el número de filas en el conjunto de datos antes de asignar los colores.

if (length(precio_fac) == nrow(datos_unidos2)) {
  datos_unidos2$color <- colores[precio_fac]
  cat("Colores asignados correctamente.\n")
} else {
  stop("Error: La longitud de precio_fac no coincide con el número de filas en datos_unidos2.")
}

## Colores asignados correctamente.

Se fusionan los datos unificados con el shapefile, y se transforma el sistema de referencia espacial a WGS 84 (EPSG:4326) para la visualización en el mapa.

datos_unidos_final <- merge(mapa_col, datos_unidos2, by = "codigo")
datos_unidos_final <- st_transform(datos_unidos_final, crs = 4326)
st_crs(datos_unidos_final)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["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]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

Finalmente, se utiliza la función leaflet para crear un mapa interactivo que muestra los departamentos con un color que representa el rango de precios de gasolina promedio. Se añade una leyenda al mapa para identificar las categorías de precios.

leaflet(datos_unidos_final) %>% 
  addTiles() %>%
  addPolygons(popup = ~DeNombre.x, color = "black", fillColor = ~color, weight = 1.1) %>%
  addLegend(position = "topright", 
            pal = colorFactor(palette = c("#66FFFF","green", "yellow", "orange", "red"), 
                              domain = levels(precio_fac)),
            values = levels(precio_fac),
            title = "Precio promedio de gasolina motor en 2023")

Terminaremos el ejercicio, presentando el valor medio de la gasolina motor en cada uno de los departamentos de Colombia:

mean <- datos_unidos2%>%
  group_by(codigo)%>%
  get_summary_stats(precio_prom_dpto, type = "mean_sd")

print(mean[c("codigo","mean")])

## # A tibble: 32 × 2
##    codigo      mean
##    <chr>      <dbl>
##  1 amazonas  15236.
##  2 antioquia 13093.
##  3 arauca    12779.
##  4 atlantico 12351.
##  5 bolivar   12595.
##  6 boyaca    13239.
##  7 caldas    13046.
##  8 caqueta   13208.
##  9 casanare  13100.
## 10 cauca     12884.
## # ℹ 22 more rows