rm(list = ls())
# paquetes que necesitan ser instalados ----
# install.packages("remotes")
# install.packages("rstan", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
# remotes::install_github("timriffe/DemoTools")

# librerias necesarias ----
library(tidyverse)
library(ggpubr)
library(DemoTools)
library(stringr)
# directorio de trabajo (debe ser modificado) ----
setwd("~/Dropbox/C\ O\ N\ S\ U\ T\ O\ R\ I\ A\ S/UNFPA-INE_Paraguay/bases")

# carga de base de datos de censos y proyecciones ----
load("censos_proy.RData")

setwd("~/Dropbox/C\ O\ N\ S\ U\ T\ O\ R\ I\ A\ S/UNFPA-INE_Paraguay/bases/figuras_analisis_edades")

# no especificados por edad ----

names(table(base0$sex))

## total nacional ----
base0 %>% 
  filter(type == "cens", year == 2022, sex != "NE", age == "NE") %>% 
  group_by(age) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  select(pob)/base0 %>% 
  filter(type == "cens", year == 2022, sex != "NE", age != "NE") %>% 
  group_by(age) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  summarise(pob = sum(pob)) %>% 
  select(pob)

## total nacional por sexo ----
base0 %>% 
  filter(type == "cens", year == 2022, sex != "NE", age == "NE") %>% 
  group_by(age, sex) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  select(pob)/base0 %>% 
  filter(type == "cens", year == 2022, sex != "NE", age != "NE") %>% 
  group_by(sex) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  select(pob)

## total por departamento por sexo ----
prop.ne <- full_join(base0 %>% 
  filter(type == "cens", year == 2022, sex != "NE", age == "NE") %>% 
  group_by(age, sex, depto) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  mutate(id = paste(sex, depto, sep = "")),
base0 %>% 
  filter(type == "cens", year == 2022, sex != "NE", age != "NE") %>% 
  group_by(sex, depto) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  mutate(id = paste(sex, depto, sep = "")), by = "id") %>% 
  select(sex = sex.x, depto = depto.x, pobne = pob.x, pob_ne = pob.y) %>% 
  mutate(propne = 100*pobne/pob_ne)

write.csv(prop.ne, file = "no_esp_depto.csv")

# no especificados por sexo ----

names(table(base0$sex))

## total nacional ----
base0 %>% 
  filter(type == "cens", year == 2022, sex == "NE", age != "NE") %>% 
  group_by(age) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  summarise(pob = sum(pob)) %>% 
  select(pob)/base0 %>% 
  filter(type == "cens", year == 2022, sex != "NE", age != "NE") %>% 
  group_by(age) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  summarise(pob = sum(pob)) %>% 
  select(pob)

## total nacional por sexo ----
base0 %>% 
  filter(type == "cens", year == 2022, sex == "NE", age != "NE") %>% 
  type_convert() %>% 
  group_by(age) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  select(pob)/base0 %>% 
  filter(type == "cens", year == 2022, sex != "NE", age != "NE") %>% 
  type_convert() %>% 
  group_by(age) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  select(pob)

## total por departamento por sexo ----
prop.ne <- full_join(base0 %>% 
                       filter(type == "cens", year == 2022, sex == "NE", age != "NE") %>% 
                       group_by(depto) %>% 
                       summarise(pob = sum(pob)) %>% 
                       ungroup() %>% 
                       mutate(id = paste(depto, sep = "")),
                     base0 %>% 
                       filter(type == "cens", year == 2022, sex != "NE", age != "NE") %>% 
                       group_by(depto) %>% 
                       summarise(pob = sum(pob)) %>% 
                       ungroup() %>% 
                       mutate(id = paste(depto, sep = "")), by = "id") %>% 
  select(depto = depto.x, pobne = pob.x, pob_ne = pob.y) %>% 
  mutate(propne = 100*pobne/pob_ne)

write.csv(prop.ne, file = "no_esp_depto.csv")



# armamos nuestra base del censo 2022 sin sexo no especificado ni edades no especificadas ----
pob2022 <- rbind(base0 %>% 
                   filter(type == "cens", sex != "NE", age != "NE") %>% 
                   type_convert() %>% 
                   filter(year == 2022, age < 100) %>% 
                   select(depto, sex, age, pob),
                 base0 %>% 
                   filter(type == "cens", sex != "NE", age != "NE") %>% 
                   type_convert() %>% 
                   filter(year == 2022, age >= 100) %>% 
                   group_by(depto, sex) %>% 
                   summarise(pob = sum(pob)) %>% 
                   mutate(age = 100, .before = pob) %>% 
                   ungroup()) %>% 
  arrange(depto, sex, age)

## Poblacion total por edades simples (deptos) ---- 
ptot22 <- pob2022 %>% 
  filter(sex != "NE", age != "NE") %>% 
  type_convert() %>% 
  group_by(depto, sex, age) %>%
  summarise(pob = sum(pob)) %>% 
  ungroup()

## Poblacion total por edades quinquenales (deptos) ---- 

ptot22_5 <- ptot22 %>% 
  mutate(
    age_group = factor(case_when(
      age >= 0 & age <= 4 ~ "0-4",
      age >= 5 & age <= 9 ~ "5-9",
      age >= 10 & age <= 14 ~ "10-14",
      age >= 15 & age <= 19 ~ "15-19",
      age >= 20 & age <= 24 ~ "20-24",
      age >= 25 & age <= 29 ~ "25-29",
      age >= 30 & age <= 34 ~ "30-34",
      age >= 35 & age <= 39 ~ "35-39",
      age >= 40 & age <= 44 ~ "40-44",
      age >= 45 & age <= 49 ~ "45-49",
      age >= 50 & age <= 54 ~ "50-54",
      age >= 55 & age <= 59 ~ "55-59",
      age >= 60 & age <= 64 ~ "60-64",
      age >= 65 & age <= 69 ~ "65-69",
      age >= 70 & age <= 74 ~ "70-74",
      age >= 75 & age <= 79 ~ "75-79",
      age >= 80 & age <= 84 ~ "80-84",
      age >= 85 & age <= 89 ~ "85-89",
      age >= 90 & age <= 94 ~ "90-94",
      age >= 95 & age <= 99 ~ "95-99",
      age >= 100 ~ "100+",
    ), levels = c(paste(seq(0, 95, 5), seq(4, 99, 5), sep = "-"), "100+"))
  ) %>% 
  group_by(depto, sex, age_group) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  arrange(depto, sex, age_group)

## primeras graficas que muestren la estructura ----

dep <- names(table(ptot22$depto))
depa <- str_replace_all(dep, " ", "")

for(i in 1:18){
assign(paste(depa[i], ".plot", sep = ""),
  ptot22 %>% 
    filter(depto == dep[i]) %>% 
  ggplot(aes(x = age, y = pob/10000, group = sex, color = sex)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  labs(x = "Edad", y = "Población (miles)", color = "Sexo", title = dep[i]) +
  theme_light() +
  scale_y_continuous(n.breaks = 10) +
  scale_x_continuous(limits = c(0, 100), breaks = (seq(0, 100, 5))))

  assign(paste(depa[i], "5", ".plot", sep = ""),
         ptot22_5 %>% 
           filter(depto == dep[i]) %>% 
  ggplot(aes(x = age_group, y = pob/10000, group = sex, color = sex)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  labs(x = "Grupo de edad", y = "", color = "Sexo", title = " ") +
  theme_light() +
  scale_y_continuous(n.breaks = 10))

  assign(paste(depa[i], "fin", ".plot", sep = ""),
    ggarrange(get(paste(depa[i], ".plot", sep = "")),
              get(paste(depa[i], "5", ".plot", sep = "")),
              ncol = 2, nrow = 1, common.legend = T, legend = "bottom")
  )
  
  ggsave(get(paste(depa[i], "fin", ".plot", sep = "")),
    filename = paste(depa[i], ".png", sep = ""), dpi = 320, width = 16.5, height = 11.5)
}

# pruebas de age-heaping por sexo----

sx <- names(table(ptot22$sex))

for (j in 1:18) {
age.test <- tibble(test = c("Whipple", "Spoorenberg", "Noumbussi_0", "Noumbussi_5", "Myers", "Bachi", "Coale_Li_0", "Coale_Li_2", "Coale_Li_5", "Sawtooth", "Roughness"), hombres = 0, mujeres = 0)

for (i in 1:2) {
  ptot22F <- ptot22 %>% 
    filter(sex == sx[i], depto == dep[j]) %>% 
    select(pob)
  
  ## preferencia digital edades centrales ----
  
  ### whipple varia de 1 (no concentracion alrededor de digitos 0 y 5) a 5 (total concentracion en digitos 0 y 5) ###
  age.test[1, i + 1] <- check_heaping_whipple(ptot22F$pob, Age = 0:100)
  ### spoorenberg es whipple modificado ###
  age.test[2, i + 1] <- check_heaping_spoorenberg(ptot22F$pob, Age = 0:100)  
  ### noumbissi es como whipple pero concentrado en un solo digito ###
  age.test[3, i + 1] <- check_heaping_noumbissi(ptot22F$pob, Age = 0:100, digit = 0)
  age.test[4, i + 1] <- check_heaping_noumbissi(ptot22F$pob, Age = 0:100, digit = 5)
  
  ### myers y bachi miden la proporcion de poblacion que tiene preferencia por un digito ###
  age.test[5, i + 1] <- check_heaping_myers(ptot22F$pob, Age = 0:100)  
  age.test[6, i + 1] <- check_heaping_bachi(ptot22F$pob, Age = 0:100)  
  
  ## preferencia digital edades avanzadas ----
  ###  coale y li miden preferencias digitales a edades avanzadas, su rango va de 1 (baja preferencia digital) a 5 (alta preferencia digital) ###
  age.test[7, i + 1] <- check_heaping_coale_li(ptot22F$pob, Age = 0:100, ageMin = 60,  ageMax = 100, digit = 0)  
  age.test[8, i + 1] <- check_heaping_coale_li(ptot22F$pob, Age = 0:100, ageMin = 60,  ageMax = 100, digit = 2)  
  age.test[9, i + 1] <- check_heaping_coale_li(ptot22F$pob, Age = 0:100, ageMin = 60,  ageMax = 100, digit = 5)  
  
  ## concentracion en edades quinquenales ----
  ### debe ser mayor a 0 ###
  age.test[10, i + 1] <- check_heaping_sawtooth(ptot22F$pob, Age = 0:100)
  ### debe ser mayor a 0.1 ###
  age.test[11, i + 1] <- check_heaping_roughness(ptot22F$pob, Age = 0:100)
  ### si sawtooth > 0 y roughness > 0.1 es posible hacer suavizamiento sin suavizar edades quinquenales ###
}

write.csv(age.test, file = paste(depa[j], "_agetest.csv", sep = ""))

# pruebas de consistencia de edades (quinquenales) y sexo ----
## poblaciones por sexo (edades quinquenales) ----
ptot22_5F <- ptot22_5 %>% 
  filter(sex == "Mujeres", depto == dep[j]) %>% 
  select(pob)

ptot22_5M <- ptot22_5 %>% 
  filter(sex == "Hombres", depto == dep[j]) %>% 
  select(pob)

## razon de sexo por edades ----
srx <- tibble(age = factor(c(paste(seq(0, 95, 5), seq(4, 99, 5), sep = "-"), "100+"), 
                           levels = c(paste(seq(0, 95, 5), seq(4, 99, 5), sep = "-"), "100+")), 
              ptot22_5M/ptot22_5F) %>% 
  select(age, srx = pob)

### grafica de razon de sexo por edades ----
srx %>% 
  ggplot(aes(x = age, y = srx, group = F)) +
  geom_line(linewidth = 1, color = "blue") +
  geom_point(size = 2.5, color = "blue") +
  geom_hline(yintercept = 1) +
  labs(x = "Edad", y = "Razón de sexo por edades", title = dep[j]) +
  theme_light()
ggsave(filename = paste(depa[j], "srx.png", sep = ""), 
       dpi = 320, width = 16.5, height = 11.5)

## razon por edades para cada sexo ----
arx <- data.frame()
for (i in 1:19) {
  arx[i, 1] <- 3*ptot22_5M$pob[i + 1]/(ptot22_5M$pob[i] + ptot22_5M$pob[i + 1] + ptot22_5M$pob[i + 2])
  arx[i, 2] <- 3*ptot22_5F$pob[i + 1]/(ptot22_5F$pob[i] + ptot22_5F$pob[i + 1] + ptot22_5F$pob[i + 2])
}

### grafica razon por edades para cada sexo ----
arx <- tibble(age = factor(c(paste(seq(5, 90, 5), seq(9, 95, 5), sep = "-"), "100+"), 
                           levels = c(paste(seq(5, 90, 5), seq(9, 95, 5), sep = "-"), "100+")), 
              arx) %>% 
  rename(Hombres = V1, Mujeres = V2) %>% 
  gather(key = "sex", value = "arx", -age)

arx %>% 
  ggplot(aes(x = age, y = arx, group = sex, color = sex)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  geom_hline(yintercept = 1) +
  labs(x = "Edad", y = "Razón de edades", color = "Sexo", title = dep[j]) +
  theme_light()

ggsave(filename = paste(depa[j], "arx.png", sep = ""), 
       dpi = 320, width = 16.5, height = 11.5)

## medidas resumen de las razones por edad, razones por sexo y razones por edad y sexo ----
accur <- tibble(accur_ind = c("ARSM", "ARSF", "SRS", "ASA"), value = 0)
### las razones por sexo y edad deben ser cercanas a 0 ###
accur[1, 2] <- ageRatioScore(ptot22_5M$pob, Age = seq(0, 100, 5), method = "UN")
accur[2, 2] <- ageRatioScore(ptot22_5F$pob, Age = seq(0, 100, 5), method = "UN")

accur[3, 2] <- sexRatioScore(Males = ptot22_5M$pob, Females = ptot22_5F$pob, Age = seq(0, 100, 5))

### <= 20 preciso, (20, 40] impreciso, >= 40 muy impreciso ###
accur[4, 2] <- ageSexAccuracy(Males = ptot22_5M$pob, Females = ptot22_5F$pob, Age = seq(0, 100, 5), method = "UN")

write.csv(accur, file = paste(depa[j], "accuracy.csv", sep = ""))
}
# Analisis de razones de cohorte (dos censos) ----

# armamos nuestra base del censo 2002 sin sexo no especificado ni edades no especificadas ----
pob2002 <- rbind(base0 %>% 
                   filter(type == "cens", sex != "NE", age != "NE") %>% 
                   type_convert() %>% 
                   filter(year == 2002, age < 100) %>% 
                   select(depto, sex, age, pob),
                 base0 %>% 
                   filter(type == "cens", sex != "NE", age != "NE") %>% 
                   type_convert() %>% 
                   filter(year == 2002, age >= 100) %>% 
                   group_by(depto, sex) %>% 
                   summarise(pob = sum(pob)) %>% 
                   mutate(age = 100, .before = pob) %>% 
                   ungroup()) %>% 
  arrange(depto, sex, age)

## Poblacion total por edades simples (nacional) ---- 
ptot02 <- pob2002 %>% 
  filter(sex != "NE", age != "NE") %>% 
  type_convert() %>% 
  group_by(depto, sex, age) %>%
  summarise(pob = sum(pob)) %>% 
  ungroup()

## Poblacion total por edades quinquenales (deptos) ---- 
ptot02_5 <- ptot02 %>% 
  mutate(
    age_group = factor(case_when(
      age >= 0 & age <= 4 ~ "0-4",
      age >= 5 & age <= 9 ~ "5-9",
      age >= 10 & age <= 14 ~ "10-14",
      age >= 15 & age <= 19 ~ "15-19",
      age >= 20 & age <= 24 ~ "20-24",
      age >= 25 & age <= 29 ~ "25-29",
      age >= 30 & age <= 34 ~ "30-34",
      age >= 35 & age <= 39 ~ "35-39",
      age >= 40 & age <= 44 ~ "40-44",
      age >= 45 & age <= 49 ~ "45-49",
      age >= 50 & age <= 54 ~ "50-54",
      age >= 55 & age <= 59 ~ "55-59",
      age >= 60 & age <= 64 ~ "60-64",
      age >= 65 & age <= 69 ~ "65-69",
      age >= 70 & age <= 74 ~ "70-74",
      age >= 75 & age <= 79 ~ "75-79",
      age >= 80 & age <= 84 ~ "80-84",
      age >= 85 & age <= 89 ~ "85-89",
      age >= 90 & age <= 94 ~ "90-94",
      age >= 95 & age <= 99 ~ "95-99",
      age >= 100 ~ "100+",
    ), levels = c(paste(seq(0, 95, 5), seq(4, 99, 5), sep = "-"), "100+"))
  ) %>% 
  group_by(depto, sex, age_group) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  arrange(depto, sex, age_group)


## tasas de crecimiento ----
for(j in 1:18){
rt_02_22 <- 100*log(ptot22_5 %>%
                      filter(depto == dep[j]) %>% 
                      summarise(pob = sum(pob))/ptot02_5 %>% 
                      filter(depto == dep[j]) %>% 
                      summarise(pob = sum(pob)))/20

rt_02_22M <- 100*log(ptot22_5 %>% 
                       filter(sex == "Hombres", depto == dep[j]) %>% 
                       summarise(pob = sum(pob))/ptot02_5 %>% 
                       filter(sex == "Hombres", depto == dep[j]) %>% 
                       summarise(pob = sum(pob)))/20
rt_02_22F <- 100*log(ptot22_5 %>% 
                       filter(sex == "Mujeres", depto == dep[j]) %>% 
                       summarise(pob = sum(pob))/ptot02_5 %>% 
                       filter(sex == "Mujeres", depto == dep[j]) %>% 
                       summarise(pob = sum(pob)))/20

write.csv(tibble(sexo = c("Total", "Hombres", "Mujeres"), 
                 tasa = c(rt_02_22$pob, rt_02_22F$pob, rt_02_22M$pob)),
          file = paste(depa[j], "rt", ".csv", sep = ""))

rxt_02_22 <- tibble(ptot22_5 %>%
                      filter(depto == dep[j], sex == "Hombres") %>% 
                      select(age_group), 100*log((ptot22_5 %>% 
                                       filter(depto == dep[j]) %>% 
                                       group_by(age_group) %>% 
                                       summarise(pob = sum(pob)) %>% 
                                       ungroup() %>% 
                                       select(pob))/(ptot02_5 %>% 
                                                       filter(depto == dep[j]) %>% 
                                                       group_by(age_group) %>% 
                                                       summarise(pob = sum(pob)) %>% 
                                                       ungroup() %>% 
                                                       select(pob)))/20) %>% 
  rename(rxt = pob)

## tasas de crecimiento por edades y sexo ----
rxt_02_22M <- tibble(ptot22_5 %>%
                       filter(depto == dep[j], sex == "Hombres") %>% 
                       select(age_group), 100*log((ptot22_5 %>% 
                                                          filter(sex == "Hombres", depto == dep[j]) %>% 
                                                          group_by(age_group) %>% 
                                                          summarise(pob = sum(pob)) %>% 
                                                          ungroup() %>% 
                                                          select(pob))/(ptot02_5 %>% 
                                                                          filter(sex == "Hombres", depto == dep[j]) %>% 
                                                                          group_by(age_group) %>% 
                                                                          summarise(pob = sum(pob)) %>% 
                                                                          ungroup() %>% 
                                                                          select(pob)))/20) %>% 
  rename(rxt = pob)

rxt_02_22F <- tibble(ptot22_5 %>%
                       filter(depto == dep[j], sex == "Hombres") %>% 
                       select(age_group), 100*log((ptot22_5 %>% 
                                                          filter(sex == "Mujeres", depto == dep[j]) %>% 
                                                          group_by(age_group) %>% 
                                                          summarise(pob = sum(pob)) %>% 
                                                          ungroup() %>% 
                                                          select(pob))/(ptot02_5 %>% 
                                                                          filter(sex == "Mujeres", depto == dep[j]) %>% 
                                                                          group_by(age_group) %>% 
                                                                          summarise(pob = sum(pob)) %>% 
                                                                          ungroup() %>% 
                                                                          select(pob)))/20) %>% 
  rename(rxt = pob)

rxt.long <- left_join(rxt_02_22F, rxt_02_22M, by = "age_group") %>% 
  rename(age = age_group, Mujeres = rxt.x, Hombres = rxt.y) %>% 
  gather(key = "sex", value = rxt, -age) 

rxt.long$rxt[is.infinite(rxt.long$rxt)] <- NA

rxt.long %>% 
  ggplot(aes(x = age, y = rxt, group = sex, fill = sex)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(x = "Edad", y = "Tasa de crecimiento media anual (%)", fill = "Sexo", title = dep[j]) +
  scale_y_continuous(n.breaks = 10) +
  scale_fill_manual(values = c("darkblue", "purple")) +
  theme_light()
ggsave(filename = paste(depa[j], "rxt", ".png", sep = ""), dpi = 320, width = 16.5, height = 11.5)
## razones de sobreviviencia de cohortes ----

csrx <- ptot22_5 %>% 
  filter(depto == dep[j]) %>% 
  group_by(age_group) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  filter(!age_group %in% c("0-4", "5-9", "10-14", "15-19")) %>% 
  select(pob)/ptot02_5 %>% 
  filter(depto == dep[j]) %>% 
  group_by(age_group) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  filter(!age_group %in% c("85-89", "90-94", "95-99", "100+")) %>% 
  select(pob)

csrxM <- ptot22_5 %>% 
  filter(sex == "Hombres", depto == dep[j]) %>% 
  group_by(age_group) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  filter(!age_group %in% c("0-4", "5-9", "10-14", "15-19")) %>% 
  select(pob)/ptot02_5 %>% 
  filter(sex == "Hombres", depto == dep[j]) %>% 
  group_by(age_group) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  filter(!age_group %in% c("85-89", "90-94", "95-99", "100+")) %>% 
  select(pob)

csrxF <- ptot22_5 %>% 
  filter(sex == "Mujeres", depto == dep[j]) %>% 
  group_by(age_group) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  filter(!age_group %in% c("0-4", "5-9", "10-14", "15-19")) %>% 
  select(pob)/ptot02_5 %>% 
  filter(sex == "Mujeres", depto == dep[j]) %>% 
  group_by(age_group) %>% 
  summarise(pob = sum(pob)) %>% 
  ungroup() %>% 
  filter(!age_group %in% c("85-89", "90-94", "95-99", "100+")) %>% 
  select(pob)

csrx.dat <- tibble(ptot22_5 %>% 
                     filter(sex == "Mujeres", depto == dep[j], !age_group %in% c("0-4", "5-9", "10-14", "15-19")) %>% 
                     select(age_group),
                   Hombres = csrxM$pob, Mujeres = csrxF$pob) %>% 
  gather(key = "sex", value = csrx, -age_group)

csrx.dat %>% 
  ggplot() +
  geom_line(aes(x = age_group, y = 100*csrx, group = sex, color = sex), linewidth = 1) +
  geom_hline(yintercept = 100) +
  scale_y_continuous(name = "Razones de Sobreviviencia de Cohortes (%)", n.breaks = 10) +
  labs(x = "Edad", color = "Sexo", title = dep[j]) +
  theme_light()
ggsave(filename = paste(depa[j], "csrxt", ".png", sep = ""), dpi = 320, width = 16.5, height = 11.5)
}
