Actividades prácticas en R

Author

Nicoll Fontalvo

Actividad #3

Análisis de los datos

Veamos cuales son las variables presentes en el DataFrame master.

master <- read_csv("master.csv", col_types = cols(sex = col_factor(levels = c("male", 
    "female")), `country-year` = col_skip()))
names(master)
 [1] "country"            "year"               "sex"               
 [4] "age"                "suicides_no"        "population"        
 [7] "suicides/100k pop"  "HDI for year"       "gdp_for_year ($)"  
[10] "gdp_per_capita ($)" "generation"

Antes de analizar cada una de las variables por separado, renombremos algunas variables que tienen nombres algo extraños.

colnames(master)[5] <- "suicidesno"
colnames(master)[8] <- "IDH"
colnames(master)[9] <- "PIB"
colnames(master)[10] <- "PIB_PC"
colnames(master)[6] <- "population"
colnames(master)[7] <- "suicides100k"

names(master)
 [1] "country"      "year"         "sex"          "age"          "suicidesno"  
 [6] "population"   "suicides100k" "IDH"          "PIB"          "PIB_PC"      
[11] "generation"  
knitr:: kable(head(master),"html", caption = "DataFrame Master") %>%
  row_spec(0, bold = TRUE, color = "black", background = "#C1CDCD", extra_css = "text-align: center") %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(2:5, width = "10em", extra_css = "padding-left: 15px; padding-right: 10px;")
DataFrame Master
country year sex age suicidesno population suicides100k IDH PIB PIB_PC generation
Albania 1987 male 15-24 years 21 312900 6.71 NA 2156624900 796 Generation X
Albania 1987 male 35-54 years 16 308000 5.19 NA 2156624900 796 Silent
Albania 1987 female 15-24 years 14 289700 4.83 NA 2156624900 796 Generation X
Albania 1987 male 75+ years 1 21800 4.59 NA 2156624900 796 G.I. Generation
Albania 1987 male 25-34 years 9 274300 3.28 NA 2156624900 796 Boomers
Albania 1987 female 75+ years 1 35600 2.81 NA 2156624900 796 G.I. Generation

Ahora si pasemos a analizar cada una de las variables de DataFrame.

Country

summary(master$country)
   Length     Class      Mode 
    27820 character character 
length(unique(master$country))
[1] 101

Vemos que tenemos registros de un total de 101 paises. Indaguemos que pasa con la variable Year

Year

unique(master$year)
 [1] 1987 1988 1989 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003
[16] 2004 2005 2006 2007 2008 2009 2010 1985 1986 1990 1991 2012 2013 2014 2015
[31] 2011 2016
length(unique(master$year))
[1] 32

Por parte de esta variable contamos con un total de 32 años de registros, empezando desde 1987 hasta 2016.

Sex

summary(master$sex)
  male female 
 13910  13910

Se evidencia que tenemos la misma cantidad de observaciones tanto de hombres como de mujeres.

Número de suicidios

Como primera variable cuantitativa tenemos es número de suicidios, veremos a continuación que los datos tienen una media de 242.6 (902.04) con una dispersión bastante grande, lo que nos lleva a pensar en la presencia de datos atípicos.

kable(descr(master$suicidesno))
suicidesno
Mean 242.574407
Std.Dev 902.047917
Min 0.000000
Q1 3.000000
Median 25.000000
Q3 131.000000
Max 22338.000000
MAD 37.065000
IQR 128.000000
CV 3.718644
Skewness 10.351794
SE.Skewness 0.014685
Kurtosis 157.128868
N.Valid 27820.000000
Pct.Valid 100.000000 
ggp1 <- ggplot(master,aes(x=master$suicidesno)) +
  geom_histogram(bins = 30 ,fill="#548B54") +
  xlab("Suicidios por año") + 
  ylab("")   

ggp2 <- ggplot(master, aes(x=master$suicidesno)) +
  geom_boxplot(fill="#43B047") +
  ggtitle("Numero de suicidios ") +
  xlab("Número de suicidios") + ylab("") 

grid.arrange(ggp1, ggp2, ncol = 2)

Como era de esperarse, hay una gran cantidad de datos atípicos, detectemos cuantos hay utilizando las puntuaciones de Z. Veamos cuantos datos atipicos existen en la variable Suicidios por año, usando la función detect_outliers_zscore.

detect_outliers_zscore <- function(data) {
  thres <- 3
  mean_value <- mean(data, na.rm = TRUE)
  std_value <- sd(data, na.rm = TRUE)
  z_scores <- (data - mean_value) / std_value
  outliers <- data[abs(z_scores) > thres]
  return(outliers)
}
val_atipicos<- detect_outliers_zscore(master$suicidesno)

cat("Los valores atípicos de la variable [Números de suicidios] fueron en total:", length(val_atipicos), "\n")
Los valores atípicos de la variable [Números de suicidios] fueron en total: 409

PIB y PIB percapital

knitr::kable(descr(master$PIB))
PIB
Mean 4.455810e+11
Std.Dev 1.453610e+12
Min 4.691962e+07
Q1 8.985353e+09
Median 4.811469e+10
Q3 2.602024e+11
Max 1.812071e+13
MAD 6.929537e+10
IQR 2.512171e+11
CV 3.262280e+00
Skewness 7.232975e+00
SE.Skewness 1.468500e-02
Kurtosis 6.421703e+01
N.Valid 2.782000e+04
Pct.Valid 1.000000e+02 
knitr::kable(descr(master$PIB_PC))
PIB_PC
Mean 1.686646e+04
Std.Dev 1.888758e+04
Min 2.510000e+02
Q1 3.447000e+03
Median 9.372000e+03
Q3 2.487400e+04
Max 1.263520e+05
MAD 1.087339e+04
IQR 2.142700e+04
CV 1.119830e+00
Skewness 1.963258e+00
SE.Skewness 1.468500e-02
Kurtosis 4.936084e+00
N.Valid 2.782000e+04
Pct.Valid 1.000000e+02 
gg1 <- ggplot(master,aes(x=PIB))+ geom_boxplot(fill="#548B54")
gg2 <- ggplot(master,aes(x=PIB_PC))+ geom_boxplot(fill="#548B54")

grid.arrange(gg1, gg2, ncol = 2)

Podemos evidenciar que estas dos variables tienen una gran presecia de datos atípicos, por lo cual la dispersión es bastante grande, aunque vemos que la desviación típica en el PIV_PC es menor.

Veamos cuantos valores atípicos estan presentes en estas variables:

val_atipicosPIB<- detect_outliers_zscore(master$PIB)
val_atipicosPIB_PC<- detect_outliers_zscore(master$PIB_PC)

cat("Los valores atípicos de la variable [PIB] fueron en total:", length(val_atipicosPIB), "\n")
Los valores atípicos de la variable [PIB] fueron en total: 492 
cat("Los valores atípicos de la variable [PIB_PC] fueron en total:", length(val_atipicosPIB_PC), "\n")
Los valores atípicos de la variable [PIB_PC] fueron en total: 492

IDH

summary(master$IDH)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  0.483   0.713   0.779   0.777   0.855   0.944   19456

Vemos que esta es la variable que presenta datos faltantes (19456), veamos su proporción.

aggr(master, numbers = TRUE, col = c("#548B54", "#43CD80"), 
     cex.axis = 0.7, gap = 3, ylab = c("Proportion of missingness", "Missingness pattern"))

Con el gráfico anterior no solo vemos que la proporción de los datos faltantes del IDH si no que tambien confirmamos que es la única variable que cuenta con datos faltantes. Veamos la distribución que sigue la variable.

ggplot(master,aes(x=IDH))+ geom_histogram(fill="#548B54")
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Warning: Removed 19456 rows containing non-finite outside the scale range
(`stat_bin()`).

Para este caso usaremos la técnica te imputación mas simple que es reemplazar los datos faltantes por la mediana, esto es:

mediana_IDH <- median(master$IDH, na.rm = TRUE)
master$IDH[is.na(master$IDH)] <- mediana_IDH
summary(master$IDH)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.4830  0.7790  0.7790  0.7783  0.7790  0.9440

Filtrando el DataFrame

Supongamos que nos interesa estudiar los suicidios que fueron registrados en los 101 paises en el año 2010.

info_2010<-master %>% filter(year =="2010")

knitr:: kable(head(info_2010),"html", caption = "DataFrame Master filtrado para el año 2010") %>%
  row_spec(0, bold = TRUE, color = "black", background = "#C1CDCD", extra_css = "text-align: center") %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(2:5, width = "10em", extra_css = "padding-left: 15px; padding-right: 10px;")
DataFrame Master filtrado para el año 2010
country year sex age suicidesno population suicides100k IDH PIB PIB_PC generation
Albania 2010 male 55-74 years 20 241852 8.27 0.722 11926953259 4359 Silent
Albania 2010 male 35-54 years 20 371611 5.38 0.722 11926953259 4359 Generation X
Albania 2010 male 25-34 years 9 179720 5.01 0.722 11926953259 4359 Generation X
Albania 2010 male 75+ years 2 50767 3.94 0.722 11926953259 4359 Silent
Albania 2010 male 15-24 years 10 279508 3.58 0.722 11926953259 4359 Millenials
Albania 2010 female 25-34 years 6 183579 3.27 0.722 11926953259 4359 Generation X 
dim(info_2010)
[1] 1056   11

Veamos como se comportan las variables en este nuevo DataFrame filtrado por el año 2010.

gg1_10 <- ggplot(info_2010,aes(x=PIB))+ geom_boxplot(fill="#548B54")
gg2_10 <- ggplot(info_2010,aes(x=PIB_PC))+ geom_boxplot(fill="#548B54")

grid.arrange(gg1_10, gg2_10, ncol = 2)

Vemos que en este nuevo DataFrame las variables que anteriormente tenian una gran cantidad de datos atípicos presentan una proporción mas bajas de los mismos. Vemos mas de cerca la distribución de la variable PIB

ggplot(info_2010,aes(x=info_2010$PIB)) +
  geom_histogram(bins = 10 ,fill="#548B54") +
  xlab("Suicidios por año") + 
  ylab("")   

summary(info_2010$sex)
  male female 
   528    528

Existe la misma proporción de hombres y mujeres.

Veamos que cambios ha sufrido los suicidios reportados en el año 2010.

ggp1 <- ggplot(info_2010,aes(x=info_2010$suicidesno)) +
  geom_histogram(bins = 10 ,fill="#548B54") +
  xlab("Suicidios por año") + 
  ylab("")   

ggp2 <- ggplot(info_2010, aes(x=info_2010$`suicides100k`)) +
  geom_boxplot(fill="#43B047")  +
  xlab("suicides/100k") + ylab("") 

grid.arrange(ggp1, ggp2, ncol = 2)

Analisemos la variable IDH.

summary(info_2010$IDH)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.6110  0.7342  0.7855  0.7932  0.8708  0.9400 
ggplot(info_2010, aes(x=info_2010$IDH)) +
  geom_boxplot(fill="#43B047")  +
  xlab("IDH") + ylab("")

Análisis complementarios:

Veamos cuales son los siete países que cuentan con mayores registros de suicidios en el 2010.

total_suicides <- info_2010 %>%
  group_by(country) %>%
  summarise(total_suicides = sum(suicidesno)) %>%
  arrange(desc(total_suicides)) %>%
  slice_head(n = 7)

ggplot(data = total_suicides, aes(x = reorder(country, -total_suicides), y = total_suicides, fill = info_2010$sex)) +
  geom_bar(stat = "identity", fill = "#43CD80") +
  labs(title = "Top 7 países con mayores registros de suicidios en 2010", x = "País", y = "Total de suicidios") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

En la siguiente gráfica se muestra los registros de los países con mayores suicidios en 2010 agrupados por el sexo.

top_countries <- info_2010 %>%
  group_by(country) %>%
  summarise(total_suicides = sum(suicidesno), .groups = 'drop') %>%
  arrange(desc(total_suicides)) %>%
  slice_head(n = 7) %>%
  pull(country)  

total_suicides <- info_2010 %>%
  filter(country %in% top_countries) %>%
  group_by(country, sex) %>%
  summarise(total_suicides = sum(suicidesno), .groups = 'drop') %>%
  arrange(desc(total_suicides))


ggplot(data = total_suicides, aes(x = reorder(country, -total_suicides), y = total_suicides, fill = sex)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = c("#2E8B57", "#66CDAA")) +
  labs(title = "Top 7 países con mayores registros de suicidios en 2010", x = "País", y = "Total de suicidios") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Veamos la cantidad de suicidios registrados por cada generación agrupados por el sexo.

suicides_generation <- info_2010 %>% 
  group_by(generation,sex) %>%
  summarise(total_s = sum(suicidesno), .groups = 'drop')

ggplot(data = suicides_generation, aes(x = generation, y = total_s, fill=sex)) +
  geom_bar(stat = "identity", ) +
  scale_fill_manual(values = c("#2E8B57", "#66CDAA"))+
  labs(title = "Proporción de suicidios registrados por cada generación", x = "Generación", y = "Número de suicidios")

kable(suicides_generation)
generation sex total_s
Generation X male 99850
Generation X female 25831
Generation Z male 1085
Generation Z female 675
Millenials male 21520
Millenials female 6164
Silent male 62025
Silent female 21552

Actividad #4

Considere la base de datos master y realice lo siguiente (utilice los operadores pipe de continuidad y compuesto):

Filtro de los datos de Colombia y EEUU:

Genere dos bases de datos, llamadas master_col y master_eu.

master_col <- master %>% filter(master$country == "Colombia")
knitr:: kable(head(master_col))
country year sex age suicidesno population suicides100k IDH PIB PIB_PC generation
Colombia 1985 male 75+ years 21 123400 17.02 0.573 34894411352 1393 G.I. Generation
Colombia 1985 male 55-74 years 113 1015200 11.13 0.573 34894411352 1393 G.I. Generation
Colombia 1985 male 25-34 years 193 2323700 8.31 0.573 34894411352 1393 Boomers
Colombia 1985 male 15-24 years 256 3190200 8.02 0.573 34894411352 1393 Generation X
Colombia 1985 male 35-54 years 188 2451100 7.67 0.573 34894411352 1393 Silent
Colombia 1985 female 15-24 years 117 3140700 3.73 0.573 34894411352 1393 Generation X 
master_eu <- master %>% filter(master$country == "United States")
knitr:: kable(head(master_eu))
country year sex age suicidesno population suicides100k IDH PIB PIB_PC generation
United States 1985 male 75+ years 2177 4064000 53.57 0.841 4.346734e+12 19693 G.I. Generation
United States 1985 male 55-74 years 5302 17971000 29.50 0.841 4.346734e+12 19693 G.I. Generation
United States 1985 male 25-34 years 5134 20986000 24.46 0.841 4.346734e+12 19693 Boomers
United States 1985 male 35-54 years 6053 26589000 22.77 0.841 4.346734e+12 19693 Silent
United States 1985 male 15-24 years 4267 19962000 21.38 0.841 4.346734e+12 19693 Generation X
United States 1985 female 35-54 years 2105 27763000 7.58 0.841 4.346734e+12 19693 Silent

Análisis de la evolución de los suicidios por cada 100.000 habitantes, del PIB per cápita y del IDH a lo largo de los años.

suicides_100k_col <- master_col %>% 
  group_by(year) %>%
  summarise(total100k = mean(`suicides100k`), .groups = 'drop')

suicides_100k_eu <- master_eu %>% 
  group_by(year) %>%
  summarise(total100k = mean(`suicides100k`), .groups = 'drop')

g1<-ggplot(suicides_100k_col, aes(x=year, y =total100k)) +
  geom_line(stat = "identity", color="#4A708B",linewidth = 1.2)+
  labs(title = "Evolucion de suicidios/100K en Colombia", x = "Años", y = "Suicidios/100k")
g2<- ggplot(suicides_100k_eu, aes(x=year, y =total100k)) +
  geom_line(stat = "identity", color="#4A708B",linewidth = 1.2)+
  labs(title = "Evolucion de suicidios/100K en Estados Unidos", x = "Años", y = "Suicidios/100k")

  grid.arrange(g1, g2, ncol = 2)

PIB_col <- master_col %>% 
  group_by(year) %>%
  summarise(pib = mean(PIB_PC, na.rm = TRUE), .groups = 'drop')

ggplot(PIB_col, aes(x=year, y =pib)) +
  geom_line(color = "#4A708B", linewidth = 1.2, linetype = "solid") + 
  geom_point(color = "#4A708B", size = 3) +
  labs(title = "Evolucion del PIB per capital a traves de los años en Colombia", x = "Años", y = "PIB per capital")

PIB_eu <- master_eu %>% 
  group_by(year) %>%
  summarise(pib = mean(PIB_PC, na.rm = TRUE), .groups = 'drop')

ggplot(PIB_eu, aes(x=year, y =pib)) +
  geom_line(color = "#4A708B", linewidth = 1.2, linetype = "solid") + 
  geom_point(color = "#4A708B", size = 3) +
  labs(title = "Evolucion del PIB per capital a traves de los años en Estados Unidos", x = "Años", y = "PIB per capital")

idh_col <- master_col %>% 
  group_by(year) %>%
  summarise(idh = mean(IDH), .groups = 'drop')

idh_eu <- master_eu %>% 
  group_by(year) %>%
  summarise(idh = mean(IDH), .groups = 'drop')

g1<-ggplot(idh_col, aes(x=year, y =idh)) +
  geom_line(stat = "identity", color="#4A708B",size = 1.2)+
  labs(title = "Evolucion de suicidios/100K en Colombia", x = "Años", y = "Suicidios/100k")
Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
g2<- ggplot(idh_eu, aes(x=year, y =idh)) +
  geom_line(stat = "identity", color="#4A708B",size = 1.2)+
  labs(title = "Evolucion de suicidios/100K en Estados Unidos", x = "Años", y = "Suicidios/100k")

  grid.arrange(g1, g2, ncol = 2)

Análisis de la evolución de los suicidios por cada 100.000 habitantes, del PIB per cápita y del IDH por género.

suicides_100k_col <- master_col %>% 
  group_by(year,sex) %>%
  summarise(total100k = mean(`suicides100k`), .groups = 'drop')

ggplot(suicides_100k_col, aes(x=year, y =total100k, color=sex)) +
  geom_line(stat = "identity")+
  geom_line(linewidth = 1)+
  scale_color_manual(values = c("#4A708B", "#CD2990"))+
  labs(title = "Evolucion de suicidios/100K a traves de los años en Colombia por Sexo", x = "Años", y = "Suicidios/100k")

suicides_100k_eu <- master_eu %>% 
  group_by(year,sex) %>%
  summarise(total100k = mean(`suicides100k`), .groups = 'drop')

ggplot(suicides_100k_eu, aes(x=year, y =total100k, color=sex)) +
  geom_line(stat = "identity")+
  geom_line(linewidth = 1)+
  scale_color_manual(values = c("#4A708B", "#CD2990"))+
  
  labs(title = "Evolucion de suicidios/100K a traves de los años en Estados Unidos por Sexo", x = "Años", y = "Suicidios/100k")

PIB_col <- master_col %>% 
  group_by(year, sex) %>%
  summarise(pib = mean(PIB_PC, na.rm = TRUE), .groups = 'drop')

ggplot(PIB_col, aes(x = year, y = pib, color = sex, linetype = sex)) +
  geom_line(linewidth = 1) +  # Usar linewidth para el grosor de la línea
  geom_point(size = 3) +
  labs(title = "Evolución del PIB per cápita a través de los años en Colombia por sexo", 
       x = "Años", y = "PIB per cápita") +
  scale_color_manual(values = c("#CD2990", "#27408B"))

PIB_eu <- master_eu %>% 
  group_by(year, sex) %>%
  summarise(pib = mean(PIB_PC, na.rm = TRUE), .groups = 'drop')

ggplot(PIB_eu, aes(x = year, y = pib, color = sex, linetype = sex)) +
  geom_line(linewidth = 1) +  # Usar linewidth para el grosor de la línea
  geom_point(size = 3) +
  labs(title = "Evolución del PIB per cápita a través de los años en Estados Unidos por sexo", 
       x = "Años", y = "PIB per cápita") +
  scale_color_manual(values = c("#CD2990", "#27408B"))

idh_col <- master_col %>% 
  group_by(year,sex) %>%
  summarise(idh = mean(IDH), .groups = 'drop')

idh_eu <- master_eu %>% 
  group_by(year,sex) %>%
  summarise(idh = mean(IDH), .groups = 'drop')

g1<-ggplot(idh_col, aes(x=year, y =idh,color=sex)) +
  geom_line(stat = "identity",size = 1.2)+
  labs(title = "Evolucion de IDH en Colombia por sexo", x = "Años", y = "IDH")+
  scale_color_manual(values = c("#CD2990", "#27408B"))

g2<- ggplot(idh_eu, aes(x=year, y =idh,color=sex)) +
  geom_line(stat = "identity",size = 1.2)+
  labs(title = "Evolucion deIDH en Estados Unidos por sexo", x = "Años", y = "IDH")+
  scale_color_manual(values = c("#CD2990", "#27408B"))

  grid.arrange(g1, g2, ncol = 2)

Análisis de la evolución de los suicidios por cada 100.000 habitantes, del PIB per cápita y del IDH por grupo de edad.

suicides_100k_col <- master_col %>% 
  group_by(year,age) %>%
  summarise(total100k = mean(`suicides100k`), .groups = 'drop')

ggplot(suicides_100k_col, aes(x = year, y = total100k, color = age, group = age)) +
  geom_line(linewidth = 1) +
  scale_color_manual(values = c("#4A708B", "#6A5ACD", "#5CACEE", "#CD2990", "#4169E1", "#27408B")) +  # Colores para las líneas
  labs(title = "Evolución de suicidios/100k a través de los años en Colombia edad", x = "Años", y = "Suicidios/100k") 

suicides_100k_eu <- master_eu %>% 
  group_by(year,age) %>%
  summarise(total100k = mean(`suicides100k`), .groups = 'drop')

ggplot(suicides_100k_eu, aes(x = year, y = total100k, color = age, group = age)) +
  geom_line(size = 1) +
  scale_color_manual(values = c("#4A708B", "#6A5ACD", "#5CACEE", "#CD2990", "#4169E1", "#27408B")) + 
  labs(title = "Evolución de suicidios/100k a través de los años en Estados unidos por edad", 
       x = "Años", y = "Suicidios/100k")

PIB_col <- master_col %>% 
  group_by(year, age) %>%
  summarise(pib = mean(PIB_PC, na.rm = TRUE), .groups = 'drop')

ggplot(PIB_col, aes(x = year, y = pib, color = age, linetype = age)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  labs(title = "Evolución del PIB per cápita a través de los años en Colombia por edad", 
       x = "Años", y = "PIB per cápita") +
  scale_color_manual(values = c("#4A708B", "#6A5ACD", "#5CACEE", "#CD2990", "#4169E1", "#27408B"))

PIB_eu <- master_eu %>% 
  group_by(year, age) %>%
  summarise(pib = mean(PIB_PC, na.rm = TRUE), .groups = 'drop')

ggplot(PIB_eu, aes(x = year, y = pib, color = age, linetype = age)) +
  geom_line(linewidth = 1.2) +  # Usar linewidth para el grosor de la línea
  geom_point(size = 3) +
  labs(title = "Evolución del PIB per cápita a través de los años en Colombia, por edad", 
       x = "Años", y = "PIB per cápita") +
  scale_color_manual(values = c("#4A708B", "#6A5ACD", "#5CACEE", "#CD2990", "#4169E1", "#27408B"))

idh_col <- master_col %>% 
  group_by(year,age) %>%
  summarise(idh = mean(IDH), .groups = 'drop')

idh_eu <- master_eu %>% 
  group_by(year,age) %>%
  summarise(idh = mean(IDH), .groups = 'drop')

g1<-ggplot(idh_col, aes(x=year, y =idh,color=age)) +
  geom_line(stat = "identity",size = 1.2)+
  geom_point(size = 3)+
  labs(title = "Evolucion de IDH en Colombia por edad", x = "Años", y = "IDH")+
  scale_color_manual(values = c("#4A708B", "#6A5ACD", "#5CACEE", "#CD2990", "#4169E1", "#27408B"))

g2<- ggplot(idh_eu, aes(x=year, y =idh,color=age)) +
  geom_line(stat = "identity",size = 1.2)+
  geom_point(size = 3)+
  labs(title = "Evolucion de IDH en Estados Unidos por edad", x = "Años", y = "IDH")+
  scale_color_manual(values = c("#4A708B", "#6A5ACD", "#5CACEE", "#CD2990", "#4169E1", "#27408B"))

  grid.arrange(g1, g2, ncol = 2)

Actividad #5

Parte 1

accidentes_bq <- read_csv("Accidentalidad_en_Barranquilla_20240830.csv", show_col_types = FALSE)

Iniciemos analizando la estructura de el DataFrame accidentes_bq

dim(accidentes_bq)
[1] 25610    11
glimpse(accidentes_bq)
Rows: 25,610
Columns: 11
$ FECHA_ACCIDENTE                    <dttm> 2018-01-01, 2018-01-01, 2018-01-01…
$ HORA_ACCIDENTE                     <chr> "01:30:00:am", "02:00:00:pm", "04:0…
$ GRAVEDAD_ACCIDENTE                 <chr> "Con heridos", "Solo daños", "Solo …
$ CLASE_ACCIDENTE                    <chr> "Atropello", "Choque", "Choque", "C…
$ SITIO_EXACTO_ACCIDENTE             <chr> "CL 87 9H 24", "CL 110 CR 46", "AV …
$ `CANT_HERIDOS_EN _SITIO_ACCIDENTE` <dbl> 1, NA, NA, NA, NA, 3, 1, NA, NA, NA…
$ `CANT_MUERTOS_EN _SITIO_ACCIDENTE` <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA,…
$ CANTIDAD_ACCIDENTES                <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,…
$ AÑO_ACCIDENTE                      <dbl> 2018, 2018, 2018, 2018, 2018, 2018,…
$ MES_ACCIDENTE                      <chr> "January", "January", "January", "J…
$ DIA_ACCIDENTE                      <chr> "Mon", "Mon", "Mon", "Mon", "Mon", …

Veamos cada una de las variables para asegurarnos que no haya un tipo de problema en los datos.

Fechas de registros:

cat("fecha mínima en los registros: \n")
fecha mínima en los registros: 
print(min(accidentes_bq$FECHA_ACCIDENTE))
[1] "2018-01-01 UTC"
cat("fecha máxima en los registros: \n")
fecha máxima en los registros: 
print(max(accidentes_bq$FECHA_ACCIDENTE))
[1] "2024-06-30 UTC"

Nuestro data frame accidentes_bq cuenta con registros desde enero del 2018 hasta junio de 2024. Haremos los análisis pertinentes solo hasta el 2023, puesto que solo tenemos informacion preeliminar del 2024.

accidentes_bq <- accidentes_bq %>%
  filter(AÑO_ACCIDENTE != 2024)

Gravedad del accidente y clase de accidente:

Veamos la proporción de cada tipo de accidente y su gravedad que fueron registrados en el df.

accidentes_bq$GRAVEDAD_ACCIDENTE<- as.factor(accidentes_bq$GRAVEDAD_ACCIDENTE)
summary(accidentes_bq$GRAVEDAD_ACCIDENTE) 
Con heridos Con muertos  Solo daños 
       9179         235       15455 
accidentes_bq$CLASE_ACCIDENTE<- as.factor(accidentes_bq$CLASE_ACCIDENTE)
summary(accidentes_bq$CLASE_ACCIDENTE) 
     Atropello Caida Ocupante         Choque       Incendio           Otro 
          1236            181          23214             13            120 
   Volcamiento 
           105 
ggplot(accidentes_bq, aes(x=GRAVEDAD_ACCIDENTE))+ 
  geom_bar(fill="#7AC5CD")+
  labs(x = "Gravedad de Accidente", y = "Frecuencia") +
  theme_minimal()

count_acc <- accidentes_bq %>%
  count(CLASE_ACCIDENTE) %>%
  arrange(desc(n))  # Ordenar de mayor a menor

# Crear el diagrama de barras descendentes
ggplot(count_acc, aes(x = reorder(CLASE_ACCIDENTE, -n), y = n)) +
  geom_bar(stat = "identity", fill = "#7AC5CD") +
  coord_flip() +  
  labs(x = "Clase de Accidente", y = "Frecuencia") +
  theme_minimal()

Cantidad de accidentes por año:

accidentes_bq$AÑO_ACCIDENTE <- as.factor(accidentes_bq$AÑO_ACCIDENTE)

acc_años<- accidentes_bq %>%
      group_by(AÑO_ACCIDENTE)%>%
      summarise(total_acc =sum(CANTIDAD_ACCIDENTES))

ggplot(acc_años, aes(reorder(AÑO_ACCIDENTE, -total_acc), y = total_acc)) +
  geom_bar(stat = "identity", fill = "#7AC5CD") +
  labs(x = "Año de Accidente", y = "Total de Accidentes", title = "Total de Accidentes por Año") +
  theme_minimal()

El gráfico anterior es muy sutil, ya que nos muestra que al parecer la cantidad de accidentes ha disminuido desde el 2018 hasta el 2023 de manera proporcional en la ciudad de Barranquilla.

acc_años1<- accidentes_bq %>%
      group_by(AÑO_ACCIDENTE,GRAVEDAD_ACCIDENTE)%>%
      summarise(total_acc =sum(CANTIDAD_ACCIDENTES), .groups = 'drop')

ggplot(data = acc_años1, aes(x = reorder(AÑO_ACCIDENTE, -total_acc), y = total_acc, fill = GRAVEDAD_ACCIDENTE)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = c("#53868B","#5F9EA0","#98F5FF")) +
  labs(title = "Top 7 países con mayores registros de suicidios en 2010", x = "País", y = "Total de suicidios") +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))+
  theme_minimal()

En el gráfico anterior se puede notar que la mayoría de los accidentes que se registraron entre 2018 y 2023 solo han presentados daños, mientras que el número de accidentes que solo han tenido heridos se ha mantenido constante durante estos años.

accidentes_bq$MES_ACCIDENTE <- as.factor(accidentes_bq$MES_ACCIDENTE)

acc_mes<- accidentes_bq %>%
      group_by(MES_ACCIDENTE)%>%
      summarise(total_acc =sum(CANTIDAD_ACCIDENTES))

ggplot(acc_mes, aes(reorder(MES_ACCIDENTE, -total_acc), y = total_acc)) +
  geom_bar(stat = "identity", fill = "#7AC5CD") +
  scale_y_continuous(limits = c(0, 2500)) +  # Ajusta el rango del eje y
  labs(x = "Mes de Accidente", y = "Total de Accidentes", title = "Total de Accidentes por Mes") +
  theme_minimal()

Según lo que se evidencia en el gráfico anterior podemos ver que la cantidad de los accidentes registrados en cada mes no es muy diferente, los accidentes registrados por mes estan estre 1900 y 2300.

accidentes_bq$DIA_ACCIDENTE <- as.factor(accidentes_bq$DIA_ACCIDENTE)

acc_dia<- accidentes_bq %>%
      group_by(DIA_ACCIDENTE,GRAVEDAD_ACCIDENTE)%>%
      summarise(total_acc =sum(CANTIDAD_ACCIDENTES), .groups = 'drop')
acc_dia
# A tibble: 21 × 3
   DIA_ACCIDENTE GRAVEDAD_ACCIDENTE total_acc
   <fct>         <fct>                  <dbl>
 1 Fri           Con heridos             1360
 2 Fri           Con muertos               32
 3 Fri           Solo daños              2419
 4 Mon           Con heridos             1410
 5 Mon           Con muertos               39
 6 Mon           Solo daños              2222
 7 Sat           Con heridos             1358
 8 Sat           Con muertos               34
 9 Sat           Solo daños              2228
10 Sun           Con heridos             1166
# ℹ 11 more rows
ggplot(acc_dia, aes(x=reorder(DIA_ACCIDENTE, -total_acc), y = total_acc,fill=GRAVEDAD_ACCIDENTE)) +
  geom_bar(stat = "identity", position = "dodge") +
  scale_fill_manual(values = c("#53868B","#5F9EA0","#98F5FF"))+ 
  labs(x = "Día de accidente", y = "Total de Accidentes", title = "Total de Accidentes por día") +
  theme_minimal()

El diagrama anterior nos perimite ver como es la gravedad de accidentes por día, viendo así que se ve la misma tendencia que se veía en las gráficas anteriores, donde la mayor cantidad de accidentes son registrados solo con daños, podemos notar que el sabado es la excepción, parece que el los sabados se registran una cantidad similar de accidentes con muertos y heridos.

Detección de datos faltantes y atípicos.

names(accidentes_bq) <- c("fecha", "hora", "gravedad", 
                          "clase", "sitio", "no_heridos", 
                          "no_muertos", "cantidad", 
                          "año", "mes", "dia")

Inicialmente veamos la proporcion de datos faltantes.

aggr(accidentes_bq, numbers = TRUE, col = c("#5F9EA0", "#7AC5CD"), 
     cex.axis = 0.7, gap = 3, ylab = c("Proportion of missingness", "Missingness pattern"))

missmap(accidentes_bq, col = c("#5F9EA0", "#7AC5CD"), legend = TRUE)
Warning: Unknown or uninitialised column: `arguments`.
Unknown or uninitialised column: `arguments`.
Warning: Unknown or uninitialised column: `imputations`.

summary(accidentes_bq$no_heridos)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
  1.000   1.000   1.000   1.474   2.000  42.000   15612 
summary(accidentes_bq$no_muertos)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
   1.00    1.00    1.00    1.03    1.00    2.00   24634

Con los análisis anteriores, vemos que los datos faltantes están en las variables no_heridos y no_muertos, donde la proporción en la variable no_muertos es del 99%, por tanto, no es posible hacer la imputación para tratarlos, por otro lado, la variable no_heridos es del 62%.

Si estudiamos bien la variable de no_heridos llegamos a la conclusión de que esta variable puede ser tomada como una variable factor, puesto que las categrías serían la cantidad de heridos en cada accidente, sería erroneo tratarla como una variable numérica continua.

accidentes_bq$no_heridos <- as.factor(accidentes_bq$no_heridos)
cate_heridos<- as.data.frame(table(accidentes_bq$no_heridos))
cate_heridos
   Var1 Freq
1     1 6421
2     2 2106
3     3  443
4     4  147
5     5   52
6     6   29
7     7   13
8     8    9
9     9    6
10   10    7
11   11    5
12   12    7
13   13    2
14   14    1
15   15    1
16   16    1
17   18    1
18   19    1
19   20    1
20   21    1
21   22    1
22   23    1
23   42    1

De esta forma si se puede hacer una mejor interpretación acerca de los registros de los números de heridos registrados en cada accidente.

Filtrando el DataFrame para el año 2023.

accidentes_bq2023 <-accidentes_bq %>%
  filter(año == 2023)

head(accidentes_bq2023)
# A tibble: 6 × 11
  fecha               hora   gravedad clase sitio no_heridos no_muertos cantidad
  <dttm>              <chr>  <fct>    <fct> <chr> <fct>           <dbl>    <dbl>
1 2023-01-01 00:00:00 07:50… Con her… Choq… CALL… 1                  NA        1
2 2023-01-01 00:00:00 11:10… Con her… Choq… CALL… 1                  NA        1
3 2023-01-02 00:00:00 02:15… Con her… Choq… CALL… 4                  NA        1
4 2023-01-02 00:00:00 07:00… Con her… Atro… CALL… 2                  NA        1
5 2023-01-02 00:00:00 09:00… Con her… Choq… CALL… 2                  NA        1
6 2023-01-02 00:00:00 10:15… Con her… Choq… CALL… 1                  NA        1
# ℹ 3 more variables: año <fct>, mes <fct>, dia <fct>
dim(accidentes_bq2023)
[1] 1662   11
accidentes_bq2023$no_heridos <-as.factor(accidentes_bq2023$no_heridos)

cate_heridos2023<- as.data.frame(table(accidentes_bq2023$no_heridos))
cate_heridos2023
   Var1 Freq
1     1 1059
2     2  417
3     3   88
4     4   33
5     5    7
6     6    6
7     7    3
8     8    1
9     9    0
10   10    1
11   11    1
12   12    2
13   13    0
14   14    0
15   15    1
16   16    0
17   18    1
18   19    0
19   20    0
20   21    1
21   22    0
22   23    1
23   42    0
num_na <- sum(is.na(accidentes_bq2023$no_heridos))

print(num_na)
[1] 40
num_na <- sum(is.na(accidentes_bq2023$no_muertos))

print(num_na)
[1] 1613

Vemos que al filtrar la base de datos accidentes_bq solo para el año 2023 accidentes_bq2023 la cantidad de datos faltantes en la variable no_heridos se redujo a solo 40, a diferencia de la variable no_muertos que mantiene una gran cantidad con 1613.

Parte 2

En esta parte del trabajo se hará un análisis de los precios del combustible en colombia en el año 2023.

combustible_1 <- read_csv("2023_I.csv",show_col_types = FALSE)
combustible_2 <- read_csv("2023_II.csv",show_col_types = FALSE)
combustible_3 <- read_csv("2023_III.csv",show_col_types = FALSE)
combustible_4 <- read_csv("2023_IV.csv",show_col_types = FALSE)
df_2023 <- rbind( combustible_1, combustible_2, combustible_3, combustible_4)

Luego de cargar nuestro dataset y unir toda la información en uno nuevo llamado df_2023 podemos empezar los análisis.

dim(df_2023)
[1] 270276      7
glimpse(df_2023)
Rows: 270,276
Columns: 7
$ BANDERA            <chr> "TERPEL", "TERPEL", "TERPEL", "TERPEL", "TERPEL", "…
$ `NOMBRE COMERCIAL` <chr> "ESTACION DE SERVICIO SERVICENTRO LA PEDRERA", "EST…
$ PRODUCTO           <chr> "DIESEL", "GASOLINA MOTOR", "GASOLINA MOTOR", "DIES…
$ `FECHA REGISTRO`   <chr> "01-Jan-2023", "01-Jan-2023", "01-Jan-2023", "01-Ja…
$ DEPARTAMENTO       <chr> "AMAZONAS", "AMAZONAS", "AMAZONAS", "AMAZONAS", "AM…
$ MUNICIPIO          <chr> "LA PEDRERA", "LA PEDRERA", "LETICIA", "LETICIA", "…
$ `VALOR PRECIO`     <dbl> 15000, 15500, 11380, 10840, 11380, 11380, 10671, 11…

Se tiene un total de 270.276 observaciones y 7 variables.

names(df_2023)
[1] "BANDERA"          "NOMBRE COMERCIAL" "PRODUCTO"         "FECHA REGISTRO"  
[5] "DEPARTAMENTO"     "MUNICIPIO"        "VALOR PRECIO"    
names(df_2023)<- c("bandera", "nombreCo", "producto","fecha","depto","munic","valor")
names(df_2023)
[1] "bandera"  "nombreCo" "producto" "fecha"    "depto"    "munic"    "valor"

Variables:

unique(df_2023$bandera)
 [1] "TERPEL"              "TEXACO"              "PRIMAX"             
 [4] "ZEUSS "              "BIOMAX"              "PETROMIL "          
 [7] "PETROBRAS"           "ESSO"                "BRIO "              
[10] "PUMA"                "PLUS MAS"            "ECOS"               
[13] "OCTANO"              "AYATAWACOOP"         "DISCOWACOOP"        
[16] "COOMULPINORT"        "DISCOM"              "P Y B"              
[19] "PETRODECOL"          "PETRDECOL"           "SAVE"               
[22] "PROXXON"             "ZAPATA Y VELASQUEZ "

En la variable BANDERA estan guaradas las empresas de distribución del combustible.

head(df_2023$nombreCo)
[1] "ESTACION DE SERVICIO SERVICENTRO LA PEDRERA"     
[2] "ESTACION DE SERVICIO SERVICENTRO LA PEDRERA"     
[3] "BALSA EL CONDOR"                                 
[4] "BALSA EL CONDOR"                                 
[5] "ESTACION DE SERVICIO DISTRIBUIDORA LOS COMUNEROS"
[6] "ESTACION DE SERVICIO DISTRIBUIDORA LOS COMUNEROS"

Se tienen registros del nombre de la estación de distribución en la variable NOMBRE COMERCIAL .

df_2023$PRODUCTO <- as.factor(df_2023$producto)
table(df_2023$producto)

        DIESEL          EXTRA GASOLINA MOTOR 
        112801          33357         124118

Se tienen regintros de los 32 departamentos de colombia.

head(unique(df_2023$depto))
[1] "AMAZONAS"    "ANTIOQUIA"   "ARAUCA"      "ATLANTICO"   "BOGOTA D.C."
[6] "BOLIVAR"

Valor de combustible por producto:

gasolina <- df_2023%>%
  filter(producto =="GASOLINA MOTOR")
head(gasolina)
# A tibble: 6 × 8
  bandera nombreCo                     producto fecha depto munic valor PRODUCTO
  <chr>   <chr>                        <chr>    <chr> <chr> <chr> <dbl> <fct>   
1 TERPEL  ESTACION DE SERVICIO SERVIC… GASOLIN… 01-J… AMAZ… LA P… 15500 GASOLIN…
2 TERPEL  BALSA EL CONDOR              GASOLIN… 01-J… AMAZ… LETI… 11380 GASOLIN…
3 TERPEL  ESTACION DE SERVICIO DISTRI… GASOLIN… 01-J… AMAZ… LETI… 11380 GASOLIN…
4 TERPEL  ESTACION DE SERVICIO DISTRI… GASOLIN… 01-J… AMAZ… LETI… 11380 GASOLIN…
5 TEXACO  EDS COMDECOM ABRIAQUI        GASOLIN… 01-J… ANTI… ABRI… 11870 GASOLIN…
6 TEXACO  ESTACIÓN DE SERVICIO Y MALL… GASOLIN… 01-J… ANTI… AMAGÁ 11200 GASOLIN…
diesel <- df_2023 %>%
  filter(producto =="DIESEL")
head(diesel)
# A tibble: 6 × 8
  bandera nombreCo                     producto fecha depto munic valor PRODUCTO
  <chr>   <chr>                        <chr>    <chr> <chr> <chr> <dbl> <fct>   
1 TERPEL  ESTACION DE SERVICIO SERVIC… DIESEL   01-J… AMAZ… LA P… 15000 DIESEL  
2 TERPEL  BALSA EL CONDOR              DIESEL   01-J… AMAZ… LETI… 10840 DIESEL  
3 TERPEL  ESTACION DE SERVICIO DISTRI… DIESEL   01-J… AMAZ… LETI… 10671 DIESEL  
4 TEXACO  EDS COMDECOM ABRIAQUI        DIESEL   01-J… ANTI… ABRI… 10910 DIESEL  
5 TEXACO  ESTACIÓN DE SERVICIO Y MALL… DIESEL   01-J… ANTI… AMAGÁ  9610 DIESEL  
6 TERPEL  ESTACION DE SERVICIO POPALI… DIESEL   01-J… ANTI… BARB…  8950 DIESEL  
extra <- df_2023 %>%
  filter(producto =="EXTRA")
head(extra)
# A tibble: 6 × 8
  bandera nombreCo                     producto fecha depto munic valor PRODUCTO
  <chr>   <chr>                        <chr>    <chr> <chr> <chr> <dbl> <fct>   
1 TERPEL  ESTACION DE SERVICIO POPALI… EXTRA    01-J… ANTI… BARB… 18760 EXTRA   
2 TERPEL  ESTACION DE SERVICIO PUERTA… EXTRA    01-J… ANTI… BARB… 19400 EXTRA   
3 PRIMAX  ESTACION DE SERVICIO PRIMAX… EXTRA    01-J… ANTI… BELLO 19580 EXTRA   
4 TERPEL  TERPEL FONTIDUENO            EXTRA    01-J… ANTI… BELLO 20520 EXTRA   
5 ZEUSS   ESTACION DE SERVICIO LA CRI… EXTRA    01-J… ANTI… CISN… 20000 EXTRA   
6 TERPEL  EDS LA COLINA S.A.S          EXTRA    01-J… ANTI… CIUD… 17000 EXTRA
Gasolina
summary(gasolina$valor)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1000    9721   11200   11512   13148   26000 
# Asegúrate de que los nombres de las columnas no tengan errores de sintaxis
ggplot(gasolina, aes(x =valor)) +
  geom_boxplot(fill = "#7AC5CD") +
  labs(x = "Producto", y = "Valor Precio") +
  theme_minimal()

Detectemos los valores atípicos usando la técnica de z-scores .

detect_outliers_zscore <- function(data) {
  thres <- 3
  mean_value <- mean(data, na.rm = TRUE)
  std_value <- sd(data, na.rm = TRUE)
  z_scores <- (data - mean_value) / std_value
  outlier_indices <- which(abs(z_scores) > thres)  # Obtiene los índices de los outliers
  return(outlier_indices)
}
out <- detect_outliers_zscore(gasolina$valor)
print(length(out))
[1] 194
cat("Indices\n")
Indices
print(head(out,10))
 [1]   916  7265  7267  9980 10030 10031 16254 17119 17120 19132

Hagamos la respectiva imputación de los datos faltantes con el método de pmm .

gasolina$valor[out] <- NaN
imp <- mice(gasolina, m=5, maxit=50, method ='pmm', seed=500, printFlag = FALSE)
Warning: Number of logged events: 7
imp_gasolina <- complete(imp)
head(imp_gasolina)
  bandera                                         nombreCo       producto
1  TERPEL      ESTACION DE SERVICIO SERVICENTRO LA PEDRERA GASOLINA MOTOR
2  TERPEL                                  BALSA EL CONDOR GASOLINA MOTOR
3  TERPEL ESTACION DE SERVICIO DISTRIBUIDORA LOS COMUNEROS GASOLINA MOTOR
4  TERPEL ESTACION DE SERVICIO DISTRIBUIDORA LOS COMUNEROS GASOLINA MOTOR
5  TEXACO                            EDS COMDECOM ABRIAQUI GASOLINA MOTOR
6  TEXACO   ESTACIÓN DE SERVICIO Y MALL SANTA LUCIA S.A.S. GASOLINA MOTOR
        fecha     depto      munic valor       PRODUCTO
1 01-Jan-2023  AMAZONAS LA PEDRERA 15500 GASOLINA MOTOR
2 01-Jan-2023  AMAZONAS    LETICIA 11380 GASOLINA MOTOR
3 01-Jan-2023  AMAZONAS    LETICIA 11380 GASOLINA MOTOR
4 01-Jan-2023  AMAZONAS    LETICIA 11380 GASOLINA MOTOR
5 01-Jan-2023 ANTIOQUIA   ABRIAQUÍ 11870 GASOLINA MOTOR
6 01-Jan-2023 ANTIOQUIA      AMAGÁ 11200 GASOLINA MOTOR
ggplot(imp_gasolina, aes(x =valor)) +
  geom_boxplot(fill = "#7AC5CD") +
  labs(x = "Producto", y = "Valor Precio") +
  theme_minimal()

Veamos otra grafica importante:

gas_bandera <- imp_gasolina %>%
  group_by(bandera) %>%
  summarise(valor = mean(valor, na.rm = TRUE)) %>%
  arrange(desc(valor)) %>%
  slice_head(n = 12)

ggplot(data = gas_bandera, aes(x = reorder(bandera, -valor), y = valor)) +
  geom_bar(stat = "identity", fill = "#7AC5CD") +
  labs(x = "Bandera", y = "Precio por Bandera") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Diesel
summary(diesel$valor)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      0    8750    9245    9145    9580   71590 
plot1 <-ggplot(diesel, aes(x =valor)) +
  geom_histogram(fill = "#7AC5CD") +
  labs(x = "Producto", y = "Valor Precio") +
  theme_minimal()
plot2 <- ggplot(diesel, aes(x =valor)) +
  geom_boxplot(fill = "#7AC5CD") +
  labs(x = "Producto", y = "Valor Precio") +
  theme_minimal()
grid.arrange(plot1, plot2, ncol = 2)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

La presencia de datos atípicos en el valor de diesel parece ser mayor, veamos la prueba:

out <- detect_outliers_zscore(diesel$valor)
print(length(out))
[1] 1058
cat("Indices\n")
Indices
print(head(out,10))
 [1]    1  941 1281 1820 2676 2750 2751 2813 2814 2817

Procedemos de nuevo a hacer la imputación debida de los datos atípicos.

diesel$valor[out] <- NaN
imp <- mice(diesel, m=5, maxit=50, method ='pmm', seed=500, printFlag = FALSE)
Warning: Number of logged events: 7
imp_diesel <- complete(imp)
head(imp_diesel)
  bandera                                         nombreCo producto       fecha
1  TERPEL      ESTACION DE SERVICIO SERVICENTRO LA PEDRERA   DIESEL 01-Jan-2023
2  TERPEL                                  BALSA EL CONDOR   DIESEL 01-Jan-2023
3  TERPEL ESTACION DE SERVICIO DISTRIBUIDORA LOS COMUNEROS   DIESEL 01-Jan-2023
4  TEXACO                            EDS COMDECOM ABRIAQUI   DIESEL 01-Jan-2023
5  TEXACO   ESTACIÓN DE SERVICIO Y MALL SANTA LUCIA S.A.S.   DIESEL 01-Jan-2023
6  TERPEL                    ESTACION DE SERVICIO POPALITO   DIESEL 01-Jan-2023
      depto      munic valor PRODUCTO
1  AMAZONAS LA PEDRERA  9999   DIESEL
2  AMAZONAS    LETICIA 10840   DIESEL
3  AMAZONAS    LETICIA 10671   DIESEL
4 ANTIOQUIA   ABRIAQUÍ 10910   DIESEL
5 ANTIOQUIA      AMAGÁ  9610   DIESEL
6 ANTIOQUIA    BARBOSA  8950   DIESEL
plot1 <-ggplot(imp_diesel, aes(x =valor)) +
  geom_histogram(fill = "#7AC5CD") +
  labs(x = "Producto", y = "Valor Precio") +
  theme_minimal()
plot2 <- ggplot(imp_diesel, aes(x =valor)) +
  geom_boxplot(fill = "#7AC5CD") +
  labs(x = "Producto", y = "Valor Precio") +
  theme_minimal()
grid.arrange(plot1, plot2, ncol = 2)
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Vemos que en este caso, a diferencia de la gasolina, parece ser que no se detectaron todos los datos, puesto que nuestra función para detectarlos solo está tomando los que están por fuera de 3 veces el RIQ. Pero aun así hubo una gran mejoría con la imputación de los datos más extremos.

Precio del diesel por bandera:

diesel_bandera <- imp_diesel %>%
  group_by(bandera) %>%
  summarise(valor = mean(valor, na.rm = TRUE)) %>%
  arrange(desc(valor)) %>%
  slice_head(n = 12)

ggplot(data = diesel_bandera, aes(x = reorder(bandera, -valor), y = valor)) +
  geom_bar(stat = "identity", fill = "#7AC5CD") +
  labs(x = "Bandera", y = "Precio por Bandera") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Vemos con este análisis que el precio del diesel por bandera no cambia mucho, estas 12 empresas con el diesel mas alto tienen precios cercanos a 9.300.

Extra
summary(extra$valor)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   1000   18890   19490   18819   19980  200000

Desde este punto ya podemos hacernos una idea de que ese valor máximo en el valor del combustible extra puede ser un error de digitación.

ggplot(extra, aes(x =valor)) +
  geom_boxplot(fill = "#7AC5CD") +
  labs(x = "Producto", y = "Valor Precio") +
  theme_minimal()

out <- detect_outliers_zscore(extra$valor)
print(length(out))
[1] 897
cat("Indices\n")
Indices
print(head(out,10))
 [1] 115 120 122 123 138 266 269 695 748 749
extra$valor[out] <- NaN
imp <- mice(extra, m=5, maxit=50, method ='pmm', seed=500, printFlag = FALSE)
Warning: Number of logged events: 7
imp_extra <- complete(imp)
head(imp_extra)
  bandera                                  nombreCo producto       fecha
1  TERPEL             ESTACION DE SERVICIO POPALITO    EXTRA 01-Jan-2023
2  TERPEL ESTACION DE SERVICIO PUERTAS DEL NORDESTE    EXTRA 01-Jan-2023
3  PRIMAX     ESTACION DE SERVICIO PRIMAX AUTOPISTA    EXTRA 01-Jan-2023
4  TERPEL                         TERPEL FONTIDUENO    EXTRA 01-Jan-2023
5  ZEUSS         ESTACION DE SERVICIO LA CRISTALINA    EXTRA 01-Jan-2023
6  TERPEL                       EDS LA COLINA S.A.S    EXTRA 01-Jan-2023
      depto          munic valor PRODUCTO
1 ANTIOQUIA        BARBOSA 18760    EXTRA
2 ANTIOQUIA        BARBOSA 19400    EXTRA
3 ANTIOQUIA          BELLO 19580    EXTRA
4 ANTIOQUIA          BELLO 20520    EXTRA
5 ANTIOQUIA       CISNEROS 20000    EXTRA
6 ANTIOQUIA CIUDAD BOLÍVAR 17000    EXTRA
ggplot(imp_extra, aes(x =valor)) +
  geom_boxplot(fill = "#7AC5CD") +
  labs(x = "Producto", y = "Valor Precio") +
  theme_minimal()

En el caso del combustible extra pasó algo similar a con el diesel, con la imputación mejoró la distribución de la variable, aun cuando se siguen presentando datos atípicos, esto es razonable, puesto que el valor del extra varía en gran manera.

extra_bandera <- imp_extra %>%
  group_by(bandera) %>%
  summarise(valor = mean(valor, na.rm = TRUE)) %>%
  arrange(valor) %>%
  slice_head(n = 12)

# Crear el gráfico de barras
ggplot(data = extra_bandera, aes(x = reorder(bandera, valor), y = valor)) +
  geom_bar(stat = "identity", fill = "#7AC5CD") +
  labs(x = "Bandera", y = "Precio por Bandera") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

Puesto que el combustible “Extra” es el que tiene el precio promedio mas alto, con la grafica anterior vemos que la bandera “SAVE” es el que tiene menor precio promedio 12000, a diferencia de las otras empresas que sus precios ya rodean los 20.000 pesos colombianos.

Precio promedio de la gasolina por departamento.

mapa_col <- st_read("COLOMBIA/COLOMBIA.shp", quiet = FALSE)
Reading layer `COLOMBIA' from data source 
  `C:\Users\fonta\OneDrive\Documentos\CDD\6. Sexto semestre\Visualizacion de datos\Corte 1\Entrega2\Entrega 2\COLOMBIA\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 84
mapa_col <- st_make_valid(mapa_col)

Habiendo cargado las coordenadas del mapa de Colombia, rectificamos que los nombres de los departamentos en ambos datasets coinciden.

dpto_1 <- tolower(imp_gasolina$depto)        
dpto_2 <- tolower(mapa_col$DPTO_CNMBR)    


M1 <- which(is.na(match(dpto_1, dpto_2)))

dpto_2 <- str_replace_all(dpto_2, "\\?", "ñ")
dpto_2 <- str_replace_all(dpto_2, "archipielago de san andres", 
                          "archipielago de san andres, santa catalina y providencia")


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

Con el proceso anterior listo, pasamos a calcular el precio promedio por departamento y lo guardamos en una nueva columna de nuestro df datos_unidos2.

datos_unidos2 <- datos_unidos %>%
  group_by(codigo, DPTO_NANO_, DPTO_NAREA, DPTO_CSMBL, DPTO_NANO, PAIS_PAIS_, SHAPE_Leng, SHAPE_Area) %>%
  summarize(precio_prom_dpto = mean(valor)) %>%
  ungroup() 
`summarise()` has grouped output by 'codigo', 'DPTO_NANO_', 'DPTO_NAREA',
'DPTO_CSMBL', 'DPTO_NANO', 'PAIS_PAIS_', 'SHAPE_Leng'. You can override using
the `.groups` argument.

Hallamos el máximo para poder hacer los intervalos de los precios, para así poder graficar mejor en el mapa.

max(datos_unidos2$precio_prom_dpto)
[1] 13586.14
precio_fac <- cut(datos_unidos2$precio_prom_dpto, 
                  breaks = c(0, 8000, 9000, 10000, 11000, 13600),
                  labels = c("a. 0-8.000", "b. 8.000-9.000", "c. 9.000-10.000", "d. 10.000-11.000", "e. 12.000-13.600"))
colores <- c("a. 0-8.000" = "#E0FFFF",
             "b. 8.000-9.000" = "#FFBBFF",
             "c. 9.000-10.000" = "#FF6EB4",
             "d. 10.000-11.000" = "#FF4500",
             "e. 12.000-13.000" = "#8B1C62")
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.

En los pasos anteriores se definieron los intervalos de los precios, se eligió la paleta de colores y se guardaron los colores en el df dependiendo de en que intervalo se encontraba el precio promedio de cada departamento. Con esto podemos hacer la última unión de data frames para luego graficar.

datos_unidos_final <- merge(mapa_col, datos_unidos2, by = "codigo")
names(datos_unidos_final)
 [1] "codigo"           "OBJECTID"         "DPTO_CCDGO"       "DPTO_NANO_.x"    
 [5] "DPTO_CNMBR"       "DPTO_CACTO"       "DPTO_NAREA.x"     "DPTO_CSMBL.x"    
 [9] "DPTO_NANO.x"      "PAIS_PAIS_.x"     "SHAPE_Leng.x"     "SHAPE_Area.x"    
[13] "DPTO_NANO_.y"     "DPTO_NAREA.y"     "DPTO_CSMBL.y"     "DPTO_NANO.y"     
[17] "PAIS_PAIS_.y"     "SHAPE_Leng.y"     "SHAPE_Area.y"     "precio_prom_dpto"
[21] "color"            "geometry"        
leaflet(datos_unidos_final) %>% 
  addTiles() %>%
  addPolygons(popup = ~DPTO_CNMBR, color = "black", fillColor = ~color, weight = 1.1) %>%
  addLegend(position = "topright", 
            pal = colorFactor(palette = c("#E0FFFF", "#FFBBFF", "#FF6EB4", "#FF4500", "#8B1C62"), 
                              domain = levels(precio_fac)),
            values = levels(precio_fac),
            title = "Precio promedio de gasolina corriente en 2019")

Gracias al mapa se puede visualizar que en la mayoría de los departamentos de nuestro país el precio promedio de la gasolina está entre 12.000 y 13.600 en el año 2023, a excepción del departamento de Nariño que su precio promedio de gasolina está entre 10.000 y 11.000 junto con el departamento del Cesar que su precio promedio está entre 9.000 y 10.000.

# Asegúrate de tener instaladas y cargadas las librerías necesarias
library(leaflet)
library(mapview)
library(webshot)
webshot::install_phantomjs()

precio_fac <- cut(datos_unidos2$precio_prom_dpto, 
                  breaks = c(0, 8000, 9000, 10000, 11000, 13600),
                  labels = c("a. 0-8.000", "b. 8.000-9.000", "c. 9.000-10.000", "d. 10.000-11.000", "e. 12.000-13.600"))

# Crear el mapa con tus datos y variables ajustadas
# Crear el mapa con tus datos y variables ajustadas
map <- leaflet(datos_unidos_final) %>%
  addTiles() %>%
  addPolygons(
    popup = ~paste(
      "Departamento: ", DPTO_CNMBR, "<br>",
      "Rango de precio promedio de gasolina corriente: ", precio_fac
    ),
    color = "black",
    fillColor = ~color,  # Utiliza el color especificado en tu dataframe
    weight = 1.1
  ) %>%
  addLegend(
    position = "topright",
    pal = colorFactor(palette = c("#E0FFFF", "#FFBBFF", "#FF6EB4", "#FF4500", "#8B1C62"), 
                      domain = levels(precio_fac)),
    values = precio_fac,
    title = "Precio promedio de gasolina corriente en 2019"
  )

# Captura del mapa como imagen
mapshot(map, file = "mapa.png", selfcontained = FALSE, vwidth = 1200, vheight = 900)

knitr::include_graphics("mapa.png")