Solución Tarea 1

Code

rm(list = ls())
library(pacman)
p_load(tidyverse, data.table, readr, tidylog, readxl,
       skimr, Hmisc, lubridate, foreign, ggtext,
       gridExtra, modelsummary, sf, stringi
       )

# tidyverse for obious reasons
# data.table
# readr for reading files easily
# tidylog for captions about data wrangling commands
# readxl for reading excel files
# skimr for skim_without_charts function
# Hmisc for describe function 
# lubridate for date processing 
# foreign for dbf databases reading 
# ggtext for text format for Gini plot 
# modelsummary for pretty tables
# sf for shapefiles processing
# stringi for accent string processing 

Problema 1

a) Muertes diarias 1998-2021

Calcule el número de muertes diarias de 1994 a 2021 y en una gráfica muestre su evolución. Ponga en la gráfica eventos importantes (e.g., cuando hay cambios de presidente, la epidemia del N1H1, la epidemia del COVID-19, y otros que considere importantes).

Leemos las bases de mortalidad en formato dbf. A partir de 1998 es cuando hay día de ocurrencia en la muerte.

Code

#lista de archivos
filenames <- list.files("tarea1/p1/raw_data/bases_defunciones", pattern="*.dbf", full.names = TRUE)

#leemos dbfs
ldf <- lapply(filenames, read.dbf)

#creamos string con nombres
#df.names <- paste("df",1998:2021,sep="") 

#df vacío
df <- data.frame()

#loop para asignar nombres y pegar bases en una sola base
for (i in 1:24){
   aux <- ldf[[i]] #creamos df en blanco
   #nos quedamos con las variables que necesitamos
   aux <- aux %>% select(DIA_OCURR, MES_OCURR, ANIO_OCUR, LISTA_MEX) 
   # assign(df.names[i], df) #asignamos nombre
   df <- rbind(df, aux)
}

rm(aux, ldf)

Cargamos catálogo de tipos de muerte para identificar homicidios. De 1994 a 1997 es la base LISBASMIN con número de catálogo 55. De 1998 en adelante es la LISTAMEX con número de catálogo 55 también.

Code

lisbasmin <- read.dbf("/Users/emiliano/Documents/Mauricio/Tareas/tarea1/p1/raw_data/defunciones_generales_base_datos_1990_1994/defunciones_base_datos_1990_dbf/LISBASMIN.dbf")

listamex <- read.dbf("/Users/emiliano/Documents/Mauricio/Tareas/tarea1/p1/raw_data/defunciones_generales_base_datos_2000_2004/defunciones_base_datos_2004_dbf/LISTAMEX.dbf")

Revisamos los rangos de anio, mes y día.

Code

describe(df$DIA_OCURR)

df$DIA_OCURR 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
14549180        0       32    0.999    15.82    10.38        2        4 
     .25      .50      .75      .90      .95 
       8       16       23       28       30 

lowest :  1  2  3  4  5, highest: 28 29 30 31 99

Code

describe(df$MES_OCURR)

df$MES_OCURR 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
14549180        0       13    0.993    6.552     4.23        1        1 
     .25      .50      .75      .90      .95 
       3        7       10       11       12 
                                                                          
Value         1.00    1.98    2.96    3.94    4.92    5.90    6.88    7.86
Frequency  1459295 1199125 1198909 1113072 1162845 1114126 1178047 1199700
Proportion   0.100   0.082   0.082   0.077   0.080   0.077   0.081   0.082
                                                  
Value         8.84    9.82   10.80   11.78   99.00
Frequency  1146128 1180762 1203437 1382261   11473
Proportion   0.079   0.081   0.083   0.095   0.001

For the frequency table, variable is rounded to the nearest 0.98

Code

describe(df$ANIO_OCUR)

df$ANIO_OCUR 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
14549180        0      123    0.998     2013    12.05     1999     2001 
     .25      .50      .75      .90      .95 
    2005     2012     2018     2020     2021 

lowest : 1900 1901 1902 1903 1904, highest: 2018 2019 2020 2021 9999

Descubrimos que los missing están como ‘99’. Para los que no tienen día asignémosles un día del 1 al 30 de forma aleatoria. Asimismo, para los que no tienen mes y año.

Code

#asignamos día aleatorio
df <- df %>% mutate(DIA_OCURR=ifelse(DIA_OCURR==99, 
                                     sample(1:30,nrow(filter(df, DIA_OCURR==99)), 
                                            replace=TRUE), DIA_OCURR))

#asignamos mes aleatorio
df <- df %>% mutate(MES_OCURR=ifelse(MES_OCURR==99, 
                                     sample(1:12,nrow(filter(df, MES_OCURR==99)), 
                                            replace=TRUE), MES_OCURR))

df <- df %>% mutate(ANIO_OCUR=ifelse(ANIO_OCUR<1998 | ANIO_OCUR>2021, 
                                     sample(1998:2021,nrow(filter(df, ANIO_OCUR<1998 | ANIO_OCUR>2021)), 
                                            replace=TRUE), ANIO_OCUR))

Graficamos muertes diarias. Para ello necesitamos colapsar la base a nivel día

Code

#base a nivel día 
df_dia <- df %>% group_by(DIA_OCURR, MES_OCURR, ANIO_OCUR) %>%
  mutate(muertes_diarias=n()) %>%
  ungroup() %>%
  dplyr::distinct(DIA_OCURR, MES_OCURR, ANIO_OCUR, .keep_all=TRUE) %>%
  arrange(ANIO_OCUR, DIA_OCURR, MES_OCURR)

#convertimos a character dia, mes y anio
df_dia <- df_dia %>% mutate(DIA_OCURR=as.character(DIA_OCURR),
                            MES_OCURR=as.character(MES_OCURR),
                            ANIO_OCUR=as.character(ANIO_OCUR))

#pegamos dia, mes y anio para hacer id de fecha
df_dia <- df_dia %>% mutate(id_fecha=paste(DIA_OCURR,MES_OCURR,ANIO_OCUR,sep = "-"))

df_dia <- df_dia %>% mutate(id_fecha=dmy(id_fecha))

#observamos fechas que tienen NA de su id de fecha
df_dia_na <- df_dia %>% filter(is.na(id_fecha))

#tiramos fechas que no existen (por eso tienen NA)
df_dia <- df_dia %>% filter(!is.na(id_fecha))

Code

ggplot(data=df_dia, aes(x=as.factor(id_fecha), 
                                  y=muertes_diarias,
                        group=1)) +
  geom_line(alpha=0.6, color="blue") +
  ggtitle('Muertes diarias 1998-2021 en México') +
  scale_x_discrete(breaks=c('2000-12-01','2006-12-01', "2009-04-23",
                            "2012-12-01", "2018-12-01", "2020-03-20")) +
  labs(y='Muertes diarias', x='Día') +
  geom_vline(xintercept = '2000-12-01', colour="red", linetype="longdash") +
  annotate(x='2000-12-01',y=+Inf,label="Pdte. V. Fox",
           vjust=15, hjust=-0.1, geom="text", size=3)+
  geom_vline(xintercept = '2006-12-01', colour="red", linetype="longdash") +
  annotate(x='2006-12-01',y=+Inf,label="Pdte. F. Calderón",
           vjust=20, hjust=1, geom="text", size=3)+
  geom_vline(xintercept = "2009-04-23", colour="red", linetype="longdash") +
  annotate(x="2009-04-23",y=+Inf,label="H1N1",
           vjust=15, hjust=-0.01, geom="text", size=3)+
  geom_vline(xintercept = "2012-12-01", colour="red", linetype="longdash") +
  annotate(x="2012-12-01",y=+Inf,label="Pdte. Peña Nieto",
           vjust=15, hjust=-0.01, geom="text", size=3)+
  geom_vline(xintercept = "2018-12-01", colour="red", linetype="longdash") +
  annotate(x="2018-12-01",y=+Inf,label="Pdte. AMLO",
           vjust=20, hjust=1, geom="text", size=3)+
  geom_vline(xintercept = "2020-03-20", colour="red", linetype="longdash") +
  annotate(x="2020-03-20",y=+Inf,label="COVID 19",
           vjust=10, hjust=1, geom="text", size=3)+
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45)) #center

b) Muertes en los 365 días del año

Graficar el número de muertos en los 365 días del calendario.

Code

#obtenemos el día del año que corresponde a la fecha 
df_dia <- df_dia %>%
  mutate(dia_anio=yday(id_fecha))

Code

ggplot(data=df_dia, aes(x=as.factor(dia_anio), 
                                  y=muertes_diarias,
                        group=ANIO_OCUR)) +
  geom_line(alpha=0.9, aes(color=ANIO_OCUR)) +
  ggtitle('Muertes diarias por año') +
  labs(y='Muertes diarias', x='Día del año') +
  scale_x_discrete(breaks=c(1, seq(0,365,25))) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) #center

c) Exceso de muertes por COVID 19

Calculamos exceso de muertes. Primero obtenemos el promedio histórico de muertes pre pandemia.

Code

#filtramos años que no son de pandemia en méxico
df_dia_preCOVID <- df_dia %>% filter(!(ANIO_OCUR) %in% c('2020', '2021'))

#creamos promedio histórico de muertes por día del año
df_dia_preCOVID <- df_dia_preCOVID %>%
  group_by(dia_anio) %>%
  summarise(mean_historico_muertes=mean(muertes_diarias))

#creamos base con muertes diarias de años de covid
df_dia_COVID <- df_dia %>% filter(ANIO_OCUR %in% c('2020', '2021'))

#creamos promedio covid de muertes por día del año
df_dia_COVID <- df_dia_COVID %>%
  group_by(dia_anio) %>%
  summarise(mean_covid_muertes=mean(muertes_diarias))

#unimos bases para graficar 
df_exceso_muertes <- inner_join(df_dia_COVID, df_dia_preCOVID)

#creamos variable de exceso de muertes 
df_exceso_muertes <- df_exceso_muertes %>%
  mutate(exceso_muertes=mean_covid_muertes - mean_historico_muertes)

Graficamos exceso de muertes por Covid

Code

ggplot(data=df_exceso_muertes, aes(x=as.factor(dia_anio), 
                                  y=exceso_muertes)) +
  geom_line(alpha=0.9, aes(group=1), color='blue') +
  ggtitle('Exceso de muertes por COVID 19 ') +
  labs(y='Exceso de muertes', x='Día del año', subtitle = 'Los años 2020 y 2021 son considerados como años de COVID ') +
  scale_x_discrete(breaks=c(1, seq(0,365,25))) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) #center

d) Exceso de homicidios

Para los homicidios repitamos lo mismo solo que ahora solo filtramos las muertes con código 55 en LISTAMEX

Code

#base a nivel día 
df_dia_homicidios <- df %>%
  filter(LISTA_MEX=='55') %>%
  group_by(DIA_OCURR, MES_OCURR, ANIO_OCUR) %>%
  mutate(homicidios_diarios=n()) %>%
  ungroup() %>%
  dplyr::distinct(DIA_OCURR, MES_OCURR, ANIO_OCUR, .keep_all=TRUE) %>%
  arrange(ANIO_OCUR, DIA_OCURR, MES_OCURR)

#convertimos a character dia, mes y anio
df_dia_homicidios <- df_dia_homicidios %>% 
  mutate(DIA_OCURR=as.character(DIA_OCURR),
                            MES_OCURR=as.character(MES_OCURR),
                            ANIO_OCUR=as.character(ANIO_OCUR))

#pegamos dia, mes y anio para hacer id de fecha
df_dia_homicidios <- df_dia_homicidios %>% 
  mutate(id_fecha=paste(DIA_OCURR,MES_OCURR,ANIO_OCUR,sep = "-"))

df_dia_homicidios <- df_dia_homicidios %>% mutate(id_fecha=dmy(id_fecha))

#tiramos fechas que no existen (por eso tienen NA)
df_dia_homicidios <- df_dia_homicidios %>% filter(!is.na(id_fecha))

#obtenemos el día del año que corresponde a la fecha 
df_dia_homicidios <- df_dia_homicidios %>%
  mutate(dia_anio=yday(id_fecha))

Graficamos evolución homicidios diarios

Code

ggplot(data=df_dia_homicidios, aes(x=as.factor(id_fecha), 
                                  y=homicidios_diarios,
                        group=1)) +
  geom_line(alpha=0.5, color="red") +
  geom_smooth(color='red') +
  ggtitle('Homicidios diarios 1998-2021 en México') +
  scale_x_discrete(breaks=c('2000-12-01','2006-12-01', "2009-04-23",
                            "2012-12-01", "2018-12-01", "2020-03-20")) +
  labs(y='Homicidios diarios', x='Día') +
  geom_vline(xintercept = '2000-12-01', colour="blue", linetype="longdash") +
  annotate(x='2000-12-01',y=+Inf,label="Pdte. V. Fox",
           vjust=15, hjust=-0.1, geom="text", size=3)+
  geom_vline(xintercept = '2006-12-01', colour="blue", linetype="longdash") +
  annotate(x='2006-12-01',y=+Inf,label="Pdte. F. Calderón",
           vjust=20, hjust=1, geom="text", size=3)+
  geom_vline(xintercept = "2009-04-23", colour="blue", linetype="longdash") +
  annotate(x="2009-04-23",y=+Inf,label="H1N1",
           vjust=15, hjust=-0.01, geom="text", size=3)+
  geom_vline(xintercept = "2012-12-01", colour="blue", linetype="longdash") +
  annotate(x="2012-12-01",y=+Inf,label="Pdte. Peña Nieto",
           vjust=15, hjust=-0.01, geom="text", size=3)+
  geom_vline(xintercept = "2018-12-01", colour="blue", linetype="longdash") +
  annotate(x="2018-12-01",y=+Inf,label="Pdte. AMLO",
           vjust=20, hjust=1, geom="text", size=3)+
  geom_vline(xintercept = "2020-03-20", colour="blue", linetype="longdash") +
  annotate(x="2020-03-20",y=+Inf,label="COVID 19",
           vjust=10, hjust=1, geom="text", size=3)+
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45)) #center

Graficamos exceso de muertes por sexenio de Felipe Calderón: pre 2006-2012 y post 2006-2012.

Code

df_dia_homicidios <- df_dia_homicidios %>% mutate(ANIO_OCUR=as.numeric(ANIO_OCUR))

#filtramos años que no son pre FECAL (FELIPE CALDERON)
df_dia_preFECAL <- df_dia_homicidios %>% filter(ANIO_OCUR < 2006)

#creamos promedio histórico de homicidios por día del año
df_dia_preFECAL <- df_dia_preFECAL %>%
  group_by(dia_anio) %>%
  summarise(mean_preFECAL_homicidios=mean(homicidios_diarios))

#filtramos años que no son post FECAL (FELIPE CALDERON)
df_dia_postFECAL <- df_dia_homicidios %>% filter(ANIO_OCUR > 2012)

#creamos promedio histórico de homicidios por día del año
df_dia_postFECAL <- df_dia_postFECAL %>%
  group_by(dia_anio) %>%
  summarise(mean_postFECAL_homicidios=mean(homicidios_diarios))

#creamos base con homicidios diarios de años de FECAL
df_dia_FECAL <- df_dia_homicidios %>% filter(ANIO_OCUR >=2006 & ANIO_OCUR<=2012)

#creamos promedio FECAL de homicidios por día del año
df_dia_FECAL <- df_dia_FECAL %>%
  group_by(dia_anio) %>%
  summarise(mean_FECAL_homicidios=mean(homicidios_diarios))

#unimos bases para graficar 
df_exceso_homicidios <- inner_join(df_dia_FECAL, df_dia_preFECAL)
df_exceso_homicidios <- inner_join(df_exceso_homicidios, df_dia_postFECAL)

#creamos variable de exceso de muertes 
df_exceso_homicidios <- df_exceso_homicidios %>%
  mutate(exceso_muertes_pre=mean_FECAL_homicidios - mean_preFECAL_homicidios,
         exceso_muertes_post=mean_FECAL_homicidios - mean_postFECAL_homicidios)

Graficamos exceso de homicidios pre sexenio Calderón

Code

ggplot(data=df_exceso_homicidios, aes(x=as.factor(dia_anio), 
                                  y=exceso_muertes_pre,
                                  group=1)) +
  geom_line(alpha=0.7, aes(group=1), color='red') +
  geom_smooth(color='red') +
  ggtitle('Exceso de homicidios por sexenio de Calderón 1998-2005') +
  labs(y='Exceso de homicidios', x='Día del año') +
  scale_x_discrete(breaks=c(1, seq(0,365,25))) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) #center

Graficamos exceso de homicidios post sexenio Calderón

Code

ggplot(data=df_exceso_homicidios, aes(x=as.factor(dia_anio), 
                                  y=exceso_muertes_post,
                                  group=1)) +
  geom_line(alpha=0.7, aes(group=1), color='red') +
  geom_smooth(color='red') +
  ggtitle('Exceso de homicidios por sexenio de Calderón 2007-2021') +
  labs(y='Exceso de homicidios', x='Día del año') +
  scale_x_discrete(breaks=c(1, seq(0,365,25))) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5)) #center

Code

rm(list=ls())

Problema 2

Calculamos número de crímenes diarioas

Code

crimenes_CDMX <- read_csv("tarea1/p2/raw_data/victimasFGJ_2023_07.csv")

#observamos las fechas de los crímenes
describe(crimenes_CDMX$anio_hecho)

crimenes_CDMX$anio_hecho 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
 1125335      386       69    0.959     2021     1.71     2019     2019 
     .25      .50      .75      .90      .95 
    2019     2021     2022     2023     2023 

lowest :  222 1917 1950 1952 1955, highest: 2019 2020 2021 2022 2023

Code

#para no complicarnos tiramos observaciones que no corresponden al intervalo de 2019-2021
crimenes_CDMX <- crimenes_CDMX %>% filter(anio_hecho>=2019 & anio_hecho<=2021)

#base a nivel día 
df_dia <- crimenes_CDMX %>% group_by(fecha_hecho) %>%
  mutate(crimenes_diarios=n()) %>%
  dplyr::distinct(fecha_hecho, .keep_all = TRUE)

a) Crímenes diarios

Graficamos

Code

ggplot(data=df_dia, aes(x=as.factor(fecha_hecho), 
                                  y=crimenes_diarios,
                        group=1)) +
  geom_line(alpha=0.6, color="blue") +
  geom_smooth(color="blue")+
  ggtitle('Crímenes diarios 2019-2021 en México') +
  scale_x_discrete(breaks=c("2019-12-24", "2020-03-20", "2020-12-24")) +
  labs(y='Crímenes diarios', x='Día') +
  geom_vline(xintercept = "2020-03-20", colour="red", linetype="longdash") +
  annotate(x="2020-03-20",y=+Inf,label="COVID 19",
           vjust=10, hjust=-0.1, geom="text", size=3)+
  geom_vline(xintercept = "2019-12-25", colour="red", linetype="longdash") +
  annotate(x="2019-12-24",y=+Inf,label="Dic 24 2019",
           vjust=10, hjust=1, geom="text", size=3)+
  geom_vline(xintercept = "2020-12-25", colour="red", linetype="longdash") +
  annotate(x="2020-12-24",y=+Inf,label="Dic 24 2020",
           vjust=10, hjust=-0.1, geom="text", size=3)+
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45)) #center

b) Crímenes más comunes

Obtengamos los crímenes más comunes

Code

top_5 <- crimenes_CDMX %>% 
  group_by(categoria) %>%
  mutate(num_delitos=n()) %>%
  ungroup() %>% arrange(-num_delitos) %>%
  dplyr::distinct(categoria, .keep_all = TRUE) %>%
  top_n(5)

print(top_5$categoria)

[1] "DELITO DE BAJO IMPACTO"                              
[2] "ROBO A TRANSEUNTE EN VÍA PÚBLICA CON Y SIN VIOLENCIA"
[3] "ROBO DE VEHÍCULO CON Y SIN VIOLENCIA"                
[4] "ROBO A NEGOCIO CON VIOLENCIA"                        
[5] "HECHO NO DELICTIVO"                                  

Code

#creamos base de top 5 dleitos más comunes 
crimenes_CDMX_top5 <- crimenes_CDMX %>%
  filter(categoria %in% top_5$categoria) %>%
  filter(!is.na(sexo))

#base a nivel día TOP 5 
df_dia_top5 <- crimenes_CDMX_top5 %>% group_by(fecha_hecho) %>%
  mutate(crimenes_diarios=n()) %>%
  dplyr::distinct(fecha_hecho, .keep_all = TRUE) 

Graficamos top 5 delitos

Code

ggplot(data=df_dia_top5, aes(x=as.factor(fecha_hecho), 
                                  y=crimenes_diarios,
                        group=1)) +
  geom_line(alpha=0.6, color="blue") +
  geom_smooth(color="blue")+
  facet_wrap(~categoria)+
  ggtitle('Crímenes diarios 2019-2021 en México') +
  scale_x_discrete(breaks=c("2019-01-01", "2020-01-01", "2021-01-01")) +
  labs(y='Crímenes diarios', x='Día') +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5),
        axis.text.x = element_text(angle = 45)) #center

c) Distribuciones de tipos de crímenes por género

Code

ggplot(df_dia_top5, aes(x=crimenes_diarios, color=sexo, fill=sexo)) +
  facet_wrap(~categoria) +
  geom_density(size=1, alpha=0.2) +
  xlab('Número crímenes') + 
  ylab('Densidad') + 
  guides(alpha='none') + 
  theme_minimal()

Problema 3

a) Visualización de cómo se hace coeficiente Gini

Solución de visualización de coeficiente Gini tomada de Ulises Quevedo y Sebastián Dulong.

Code

rm(list = ls())

Code

#Jalamos los datos
ingresos <- read_csv("https://datos.cdmx.gob.mx/dataset/52f97506-ef52-449a-b2c9-29b6197abaa6/resource/7f0d7073-6861-4d8c-9ca5-17ebe3ff2388/download/remuneraciones_da_qna_10_23.csv")

#Nos quedamos con las variables relevantes
ingresos <- ingresos %>% select(n_puesto,n_cabeza_sector,sueldo_tabular_bruto)

#Coeficiente de gini
DescTools::Gini(ingresos$sueldo_tabular_bruto)

[1] 0.3076823

Code

#Gini plot
ggplot() + geom_abline(slope=1,intercept=0) + 
  coord_cartesian(ylim = c(0,1),xlim=c(0,1))+
  theme_classic()+
  labs(title=str_wrap("1 In a perfectly equal society, X% of people accumulate X% of total wealth",60),
    x="Percentage of people",y="Percentage of wealth",
    subtitle = "People are sorted based on their wealth") +
    geom_point(aes(x=0.25,y=0.25),color="#0072B2",size=3)+
    annotate("text",x=0.5,y=0.25,label=str_wrap("The 25%  has 25% of total wealth",30),size=3)+
    geom_point(aes(x=0.75,y=0.75),color="#0072B2",size=3)+
    annotate("text",x=0.8,y=0.5,label=str_wrap("The top 25% has 25% of total wealth",30),size=3)

Code

ggplot() + geom_abline(slope=1,intercept=0) + 
  coord_cartesian(ylim = c(0,1),xlim=c(0,1))+
  theme_classic()+
  labs(title=str_wrap("2 As societies become more unequal, the lower % amass a smaller percentage of wealth",60),
       x="Percentage of people",y="Percentage of wealth") +
  geom_point(aes(x=0.25,y=0.1),color="#0072B2",size=3)+
  annotate("text",x=0.45,y=0.15,label=str_wrap("The 25th percentile now has only 10% of total wealth",30),size=3)+
  geom_point(aes(x=0.75,y=0.60),color="#0072B2",size=3)+
  annotate("text",x=0.8,y=0.5,label=str_wrap("The top 25th percentile now has 40% of total wealth",30),size=3)

Code

x <- seq(0,1,0.0001)
lorenz.df <- data.frame(x=x,ymin=(0.2)*x + 0.8*x^2,ymax=x)

ggplot() + geom_abline(slope=1,intercept=0) + 
                    coord_cartesian(ylim = c(0,1),xlim=c(0,1))+
                    stat_function(fun=function(x)  (0.2)*x + 0.8*x^2,
                                  color="#0072B2",
                                  linewidth=1.5)+
                    theme_classic()+
                    geom_ribbon(data=lorenz.df,aes(x=x,ymin=ymin,ymax=ymax),
                                alpha=0.2,fill="navyblue")+
                    geom_ribbon(data=lorenz.df,aes(x=x,ymin=0,ymax=ymin),
                                alpha=0.2,fill="red")+
                    labs(title=str_wrap("3 Doing this for all points defines the  <span style='color:#0072B2;'>**Lorenz Curve**</span>",60),
                         subtitle="The futher the Lorenz Curve is from our equality line,\nless equal a society is",
                         x="Percentage of people",y="Percentage of wealth") +
                    scale_x_continuous(expand=c(0,0))+scale_y_continuous(expand=c(0,0))+
                    theme(plot.title = element_markdown(lineheight = 1.1))

Code

ggplot() + geom_abline(slope=1,intercept=0) + 
  coord_cartesian(ylim = c(0,1),xlim=c(0,1))+
  stat_function(fun=function(x)  (0.2)*x + 0.8*x^2,
                color="#0072B2",
                linewidth=1.5)+
  theme_classic()+
  geom_ribbon(data=lorenz.df,aes(x=x,ymin=ymin,ymax=ymax),
              alpha=0.2,fill="navyblue")+
  annotate("text",x=0.5,y=0.4,color="navyblue",label="A",size=7)+
  geom_ribbon(data=lorenz.df,aes(x=x,ymin=0,ymax=ymin),
              alpha=0.2,fill="red")+
  annotate("text",x=0.75,y=0.25,color="red",label="B",size=7)+
  annotate("text",x=0.3,y=0.75,label="Gini = A/(A+B)",size=4)+
  labs(title=stringr::str_wrap("4 The Gini Coefficient is the ratio between 
                      <span style='color:#000080;'>Area A</span> 
                      <br>and the total area under the equality line 
                      (<span style='color:#000080;'>A</span>+<span style='color:#FF0000;'>B</span>)",20),
       subtitle="The greater the Gini Coefficient, less equal a society",
       x="Percentage of people",y="Percentage of wealth") +
  scale_x_continuous(expand=c(0,0))+scale_y_continuous(expand=c(0,0)) +
  theme(plot.title = element_markdown(lineheight = 1.1))

b) Histograma de ingresos

Code

#Histograma
ggplot(ingresos)+
  geom_histogram(aes(x=sueldo_tabular_bruto),bins=50, color="darkblue", fill="lightblue")+
  scale_x_continuous(labels=scales::dollar_format())+
  labs(x="Sueldo bruto ($ MXN)",y="Frecuencia",
       title="Histograma de sueldos")+
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5))

c) Percentiles de sueldo

Code

#Calcule con qué porcentaje del ingreso se queda el: 0.1 % con mayores ingresos, el 1 % con mayores

# ingresos, el 5 % con mayores ingresos y el 10 % con mayores ingresos.
percentil <- c(0.001,0.01,0.05,0.1)
total_ingresos <- sum(ingresos$sueldo_tabular_bruto)

top_punto_1 <- ingresos %>% top_frac(sueldo_tabular_bruto,n=0.001) %>% summarise(valor=sum(sueldo_tabular_bruto)) %>% pull()
top_1 <- ingresos %>% top_frac(sueldo_tabular_bruto,n=0.01) %>% summarise(valor=sum(sueldo_tabular_bruto)) %>% pull()
top_5 <- ingresos %>% top_frac(sueldo_tabular_bruto,n=0.05) %>% summarise(valor=sum(sueldo_tabular_bruto)) %>% pull()
top_10 <- ingresos %>% top_frac(sueldo_tabular_bruto,n=0.1) %>% summarise(valor=sum(sueldo_tabular_bruto)) %>% pull()


##Tablita: 
distribucion_ingreso <- data.frame(
  top_porcentaje = percentil,
  porcentaje_ingresos = c(
    top_punto_1/total_ingresos,
    top_1/total_ingresos,
    top_5/total_ingresos,
    top_10/total_ingresos
  )*100
)

datasummary_df(distribucion_ingreso, 
               title = "Acumulación de ingresos por percentil", 
               fmt = 3,
               align = 'cc')

Acumulación de ingresos por percentil

top_porcentaje	porcentaje_ingresos
0.001	0.983
0.010	5.449
0.050	16.104
0.100	25.398

Bono

Crearemos una visualización del coeficiente de Gini de las alcaldías de la CDMX.

Code

#leemos shapefile de alcaldía 
cdmx_map <- st_read("tarea1/p3/raw_data/alcaldias/alcaldias.shp")

Reading layer `alcaldias' from data source 
  `/Users/emiliano/Documents/Mauricio/Tareas/tarea1/p3/raw_data/alcaldias/alcaldias.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 16 features and 4 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 461655.7 ymin: 2106166 xmax: 506274.7 ymax: 2166422
Projected CRS: WGS 84 / UTM zone 14N

Code

#calculamos Gini para cada alcaldía

ingresos_alcaldia <- ingresos %>% 
  filter(grepl("ALCALDIA",n_cabeza_sector)) %>%
  group_by(n_cabeza_sector) %>%
  mutate(gini=DescTools::Gini(sueldo_tabular_bruto)) %>%
  dplyr::distinct(n_cabeza_sector, .keep_all = TRUE) %>%
  mutate(n_cabeza_sector=toupper(gsub("ALCALDIA DE ","", n_cabeza_sector)),
         n_cabeza_sector=stri_trans_general(str = n_cabeza_sector, 
                                   id = "Latin-ASCII")) %>%
  select(n_cabeza_sector, gini)

cdmx_map <- cdmx_map %>%
  mutate(NOMBRE=toupper(stri_trans_general(str = NOMBRE, 
                                   id = "Latin-ASCII")))

#unimos base de shapefile y base que tiene índice gini
cdmx_map <- inner_join(cdmx_map,ingresos_alcaldia, by=c('NOMBRE'='n_cabeza_sector'))

Code

#mapa cdmx 
cdmx_map %>%
  # usamos el aesthetic fill para indicar la columna de casos
  ggplot(aes(fill = gini)) +
  geom_sf() +
  geom_sf_label(aes(label = NOMBRE), colour = "white")+
  labs(y=' ', x=' ') +
  ggtitle('Desigualdad de ingresos en la CDMX') +
  guides(fill = guide_legend(title = "Coef. Gini")) +
  theme_minimal() +
  theme(axis.text.x = element_blank(),
        axis.ticks.x = element_blank(),
        axis.text.y = element_blank(), 
        axis.ticks.y = element_blank())

Problema 4

Esta solución es tomada de la solución hecha por Ernesto Anaya y Guillermo Naranjo

Code

#para este ejercicio, se crearon 3 funciones:
# 1. Para juntar los resultados de las simulaciones en un df
# 2. Crear y guardar las visualizaciones
# 3. Correr las simulaciones 


create_df <- function(norm_means, chisq_means, exp_means, unif_means){
  #guardar los simulacioness como df
  
  
  means_list <- list('Normal' = as.numeric(norm_means),
                     'Xicuadrada' = chisq_means,
                     'Exponencial' = exp_means,
                     'Uniforme' = unif_means)
  
  df_means <- as.data.frame(do.call(cbind, means_list))
  
  #normalizar datos
  df_means <- as.data.frame(lapply(df_means, as.numeric))
  df_means <- as.data.frame(scale(df_means))
  
  
  return(df_means)
}


# HACER GRAFICAS

histogramas_sim <- function(df_means, n){
  #Crear grid de los histogramas de las simulaciones (normalizadas) y la verdadera distribucion
  
  x_lab <- "Valor"
  y_lab <- "Densidad"
  l_width <- 1
  
  plot_norm <- ggplot(df_means, aes(x = as.numeric(Normal))) +
    geom_histogram(aes(y = ..density..),bins=60, color = "black", fill = "white", center=0) + 
    labs(x = x_lab, y = y_lab, title = 'Normal') +
    stat_function(fun=dnorm, 
                  args=list(mean=mean(df_means$Normal), 
                            sd=sd(df_means$Normal)),
                  col = '#a8324c', linewidth=l_width
    ) +
    theme_minimal()
  
  plot_chisq <- ggplot(df_means, aes(x = as.numeric(Xicuadrada))) +
    geom_histogram(aes(y = ..density..),bins=60, color = "black", fill = "white", center=0) + 
    labs(x = x_lab, y = y_lab, title = 'Xi-Cuadrada') +
    stat_function(fun=dnorm, 
                  args=list(mean=mean(df_means$Xicuadrada), 
                            sd=sd(df_means$Xicuadrada)),
                  col = '#a8324c', linewidth=l_width
    ) +
    theme_minimal()
  
  plot_exp <- ggplot(df_means, aes(x = as.numeric(Exponencial))) +
    geom_histogram(aes(y = ..density..),bins=60, color = "black", fill = "white", center=0) + 
    labs(x = x_lab, y = y_lab, title = 'Exponencial') +
    stat_function(fun=dnorm, 
                  args=list(mean=mean(df_means$Exponencial), 
                            sd=sd(df_means$Exponencial)),
                  col = '#a8324c', linewidth=l_width
    ) +
    theme_minimal()
  
  plot_unif <- ggplot(df_means, aes(x = as.numeric(Uniforme))) +
    geom_histogram(aes(y = ..density..),bins=60, color = "black", fill = "white", center=0) + 
    labs(x = x_lab, y = y_lab, title = 'Uniforme') +
    stat_function(fun=dnorm, 
                  args=list(mean=mean(df_means$Uniforme), 
                            sd=sd(df_means$Uniforme)),
                  col = '#a8324c', linewidth=l_width
    ) +
    theme_minimal()
  
  plot_grid <- grid.arrange(plot_norm, plot_chisq, plot_exp, plot_unif,  
                            ncol=2, top=paste("Histogramas de simulaciones con n =", n)
  )
  
  
}


#simulaciones
run_simulations <- function(){
  #correr las simulaciones y guardar las visualizaciones de los resultados
  
  ns <- c(1,10,100,1000,4000)
  
  for(i in ns) {
    
    norm_means <- list()
    chisq_means <- list()
    exp_means <- list()
    unif_means <- list()
    
    for(j in 1:i) {
      
      df_simulation <- data.frame(x_norm = rnorm(10000, mean = -2, sd = 3),
                                  x_chisq = rchisq(10000, df=15),
                                  x_exp = rexp(10000, rate = 3),
                                  x_unif = runif(10000, min=-5, max=8)
      )
      norm_means <- append(norm_means, mean(df_simulation$x_norm))
      chisq_means <- append(chisq_means, mean(df_simulation$x_chisq))
      exp_means <- append(exp_means, mean(df_simulation$x_exp))
      unif_means <- append(unif_means, mean(df_simulation$x_unif))
      
    }
    
    df_means <- create_df(norm_means, chisq_means, exp_means, unif_means)
    hist <- histogramas_sim(df_means, i)
  }
  
  return(hist)
}

res <- run_simulations()