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)
# 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")

# 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 (nacional) ---- 
ptot22 <- pob2022 %>% 
  filter(sex != "NE", age != "NE") %>% 
  type_convert() %>% 
  group_by(sex, age) %>%
  summarise(pob = sum(pob)) %>% 
  ungroup()

## Poblacion total por edades quinquenales (nacional) ---- 
ptot22_5 <- rbind(ptot22 %>% 
                    filter(sex == "Hombres") %>%
                    select(pob) %>% 
                    as.matrix() %>% 
                    groupAges(N = 5, 0:100) %>% 
                    as_tibble() %>% 
                    mutate(sex = "Hombres",
                           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+")),
                           .before = "value"),
                  ptot22 %>% 
                    filter(sex == "Mujeres") %>%
                    select(pob) %>% 
                    as.matrix() %>% 
                    groupAges(N = 5, 0:100) %>% 
                    as_tibble() %>% 
                    mutate(sex = "Mujeres",
                           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+")),
                           .before = "value"))

## primeras graficas que muestren la estructura ----
ptot22.plot <- ptot22 %>% 
  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") +
  theme_light() +
  scale_y_continuous(n.breaks = 10) +
  scale_x_continuous(limits = c(0, 100), breaks = (seq(0, 100, 5)))

ggsave(filename = "px_cens2022.png", dpi = 320, width = 16.5, height = 11.5)

ptot22_5.plot <- ptot22_5 %>% 
  ggplot(aes(x = age, y = value/10000, group = sex, color = sex)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2.5) +
  labs(x = "Grupo de edad", y = "", color = "Sexo") +
  theme_light() +
  scale_y_continuous(n.breaks = 10)

ggsave(filename = "px5_cens2022.png", dpi = 320, width = 16.5, height = 11.5)

ptot22_nal.plot <- ggarrange(ptot22.plot, ptot22_5.plot, ncol = 2, nrow = 1, common.legend = T,
                             legend = "bottom")

ggsave(ptot22_nal.plot, file = "ptot22_nal.png", dpi = 320, )


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

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

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]) %>% 
  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 = "age_test.csv")

# pruebas de consistencia de edades (quinquenales) y sexo ----
## poblaciones por sexo (edades quinquenales) ----
ptot22_5F <- ptot22_5 %>% 
  filter(sex == "Mujeres") %>% 
  select(pob = value)

ptot22_5M <- ptot22_5 %>% 
  filter(sex == "Hombres") %>% 
  select(pob = value)

## 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") +
  theme_light()
ggsave(filename = "srx.png", 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") +
  theme_light()

ggsave(filename = "arx.png", 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 = "accuracy.csv")

# 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(sex, age) %>%
  summarise(pob = sum(pob)) %>% 
  ungroup()

## Poblacion total por edades quinquenales (nacional) ---- 
ptot02_5 <- rbind(ptot02 %>% 
                    filter(sex == "Hombres") %>%
                    select(pob) %>% 
                    as.matrix() %>% 
                    groupAges(N = 5, 0:100) %>% 
                    as_tibble() %>% 
                    mutate(sex = "Hombres",
                           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+")),
                           .before = "value"),
                  ptot02 %>% 
                    filter(sex == "Mujeres") %>%
                    select(pob) %>% 
                    as.matrix() %>% 
                    groupAges(N = 5, 0:100) %>% 
                    as_tibble() %>% 
                    mutate(sex = "Mujeres",
                           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+")),
                           .before = "value"))

## tasas de crecimiento ----
rt_02_22 <- 100*log(ptot22_5 %>% 
                      summarise(pob = sum(value))/ptot02_5 %>% 
                      summarise(pob = sum(value)))/20

rt_02_22M <- 100*log(ptot22_5 %>% 
                      filter(sex == "Hombres") %>% 
                      summarise(pob = sum(value))/ptot02_5 %>% 
                      filter(sex == "Hombres") %>% 
                      summarise(pob = sum(value)))/20
rt_02_22F <- 100*log(ptot22_5 %>% 
                       filter(sex == "Mujeres") %>% 
                       summarise(pob = sum(value))/ptot02_5 %>% 
                       filter(sex == "Mujeres") %>% 
                       summarise(pob = sum(value)))/20

write.csv(tibble(sexo = c("Total", "Hombres", "Mujeres"), 
                 tasa = c(rt_02_22$pob, rt_02_22F$pob, rt_02_22M$pob)), file = "tasas_nal.csv")

rxt_02_22 <- tibble(age = ptot22_5$age[1:21], 100*log((ptot22_5 %>% 
                      group_by(age) %>% 
                      summarise(pob = sum(value)) %>% 
                      ungroup() %>% 
                        select(pob))/(ptot02_5 %>% 
                      group_by(age) %>% 
                      summarise(pob = sum(value)) %>% 
                      ungroup() %>% 
                        select(pob)))/20) %>% 
  rename(rxt = pob)

## tasas de crecimiento por edades y sexo ----
rxt_02_22M <- tibble(age = ptot22_5$age[1:21], 100*log((ptot22_5 %>% 
                                                          filter(sex == "Hombres") %>% 
                                                         group_by(age) %>% 
                                                         summarise(pob = sum(value)) %>% 
                                                         ungroup() %>% 
                                                         select(pob))/(ptot02_5 %>% 
                                                                         filter(sex == "Hombres") %>% 
                                                                         group_by(age) %>% 
                                                                         summarise(pob = sum(value)) %>% 
                                                                         ungroup() %>% 
                                                                         select(pob)))/20) %>% 
  rename(rxt = pob)

rxt_02_22F <- tibble(age = ptot22_5$age[1:21], 100*log((ptot22_5 %>% 
                                                          filter(sex == "Mujeres") %>% 
                                                          group_by(age) %>% 
                                                          summarise(pob = sum(value)) %>% 
                                                          ungroup() %>% 
                                                          select(pob))/(ptot02_5 %>% 
                                                                          filter(sex == "Mujeres") %>% 
                                                                          group_by(age) %>% 
                                                                          summarise(pob = sum(value)) %>% 
                                                                          ungroup() %>% 
                                                                          select(pob)))/20) %>% 
  rename(rxt = pob)

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

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") +
  scale_y_continuous(n.breaks = 10) +
  scale_fill_manual(values = c("darkblue", "purple")) +
  theme_light()
  
ggsave(filename = "rxt_nal.png", dpi = 320)

## razones de sobreviviencia de cohortes ----

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

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

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

csrx.dat <- tibble(age = factor(levels(ptot22_5$age)[-c(1:4)],
                                levels = levels(ptot22_5$age)[-c(1:4)]),
  Hombres = csrxM$pob, Mujeres = csrxF$pob) %>% 
  gather(key = "sex", value = csrx, -age)

csrx.dat %>% 
  ggplot() +
  geom_line(aes(x = age, y = 100*csrx, group = sex, color = sex), linewidth = 1) +
  geom_hline(yintercept = 100) +
  scale_y_continuous(limits = c(0, 100), name = "Razones de Sobreviviencia de Cohortes (%)", n.breaks = 10) +
  labs(x = "Edad", color = "Sexo") +
  theme_light()
ggsave(filename = "csrx_nal.png", dpi = 320)
