students <- listed_df %>%
  unnest(cols = '.') %>% # unlist and merge data.frames
  janitor::clean_names() %>%
  select(co_ies, nu_ano_censo, no_regiao, no_uf, starts_with('qt_in')) %>%
  right_join(regions_states, by = c('co_ies', 'no_regiao', 'no_uf')) %>%
  group_by(nu_ano_censo, no_regiao, no_uf) %>%
  summarise_all(sum, na.rm = T) %>%
  ungroup()

head(students) 

students %>%
  mutate(data = zoo::as.yearmon(nu_ano_censo)) %>%
  group_by(data) %>%
  summarise_at(c('qt_ing_fem', 'qt_ing_masc'), sum) %>%
  gather('var', 'value', 2:3) %>%
  ggplot(aes(x = data, 
             y = value,
             linetype = factor(var, 
                        labels = c('Female', 'Male')),
             fill = factor(var, 
                        labels = c('Female', 'Male')))) +
  geom_bar(stat = 'identity', position = 'dodge') +
  geom_smooth(method = 'loess', col = 'black') +
  ggthemes::scale_fill_ptol() +
  ggthemes::theme_clean() +
  guides(fill = guide_legend('Sex'),
         linetype = guide_legend('Sex')) +
  xlab('') +
  scale_y_continuous(name = 'Number of students admitted', 
                     labels = scales::label_number(suffix = "K", scale = 1e-3, big.mark = '.'),
                     breaks = seq(0,1500000, by = 250000)) +
  ggtitle('Number of Students by sex')

students %>%
  mutate(data = zoo::as.yearmon(nu_ano_censo)) %>%
  group_by(data, no_regiao) %>%
  summarise_at(c('qt_ing_fem', 'qt_ing_masc'), sum) %>%
  gather('var', 'value', 3:4) %>%
  ggplot(aes(x = data, 
             y = value,
             linetype = factor(var, 
                        labels = c('Female', 'Male')),
             fill = factor(var, 
                        labels = c('Female', 'Male')))) +
  geom_bar(stat = 'identity', position = 'dodge') +
  geom_smooth(method = 'loess', col = 'black') +
  facet_wrap(~ no_regiao, scales = 'free_y') +
  ggthemes::scale_fill_ptol() +
  ggthemes::theme_clean() +
  guides(fill = guide_legend('Sex'),
         linetype = guide_legend('Sex')) +
  xlab('') +
  scale_y_continuous(name = 'Number of students admitted', 
                     labels = scales::label_number(suffix = " K", scale = 1e-3)) +
  ggtitle('Number of Students by sex and Region')

students %>%
  mutate(data = zoo::as.yearmon(nu_ano_censo)) %>%
  group_by(data, no_uf) %>%
  summarise_at(c('qt_ing_fem', 'qt_ing_masc'), sum) %>%
  gather('var', 'value', 3:4) %>%
  ggplot(aes(x = data, 
             y = value,
             linetype = factor(var, 
                        labels = c('Female', 'Male')),
             fill = factor(var, 
                        labels = c('Female', 'Male')))) +
  geom_bar(stat = 'identity', position = 'dodge') +
  geom_smooth(method = 'loess', col = 'black') +
  facet_wrap(~ no_uf, scales = 'free_y') +
  ggthemes::scale_fill_ptol() +
  ggthemes::theme_clean() +
  guides(fill = guide_legend('Sex'),
         linetype = guide_legend('Sex')) +
  xlab('') +
  scale_y_continuous(name = 'Number of students admitted', 
                     labels = scales::label_number(suffix = " K", scale = 1e-3)) +
  ggtitle('Number of Students by sex and States')

students %>%
  mutate(data = zoo::as.yearmon(nu_ano_censo)) %>%
  group_by(data) %>%
  summarise_at(c('qt_ing_branca', 'qt_ing_preta', 'qt_ing_indigena',
                 'qt_ing_parda', 'qt_ing_amarela', 'qt_ing_cornd'), sum) %>%
  gather('var', 'value', 2:7) %>%
  mutate(var = case_when(var == 'qt_ing_branca' ~ 'White',
                         var == 'qt_ing_preta' ~ 'Black',
                         var == 'qt_ing_indigena' ~ 'Indigenous',
                         var == 'qt_ing_parda' ~ 'Brown', 
                         var == 'qt_ing_amarela' ~ 'Asian', 
                         var == 'qt_ing_cornd' ~ 'Not\nDefined')) %>%
  ggplot(aes(x = data, 
             y = value,
             linetype = var,
             fill = var)) +
  geom_bar(stat = 'identity', position = 'dodge') +
  geom_smooth(method = 'loess', col = 'black') +
  ggthemes::scale_fill_ptol() +
  ggthemes::theme_clean() +
  facet_wrap( ~ var, scales = 'free_y') +
  guides(fill = F,
         linetype = F) +
  xlab('') +
  scale_y_continuous(name = 'Number of students admitted', 
                     labels = scales::label_number(suffix = "K", scale = 1e-3, big.mark = '.')) +
  ggtitle('Number of Students by ethnic group')

students %>%
  mutate(data = zoo::as.yearmon(nu_ano_censo)) %>%
  group_by(data, no_regiao) %>%
  summarise_at(c('qt_ing_branca', 'qt_ing_preta', 'qt_ing_indigena',
                 'qt_ing_parda', 'qt_ing_amarela', 'qt_ing_cornd'), sum) %>%
  gather('var', 'value', 3:8) %>%
  mutate(var = case_when(var == 'qt_ing_branca' ~ 'White',
                         var == 'qt_ing_preta' ~ 'Black',
                         var == 'qt_ing_indigena' ~ 'Indigenous',
                         var == 'qt_ing_parda' ~ 'Brown', 
                         var == 'qt_ing_amarela' ~ 'Asian', 
                         var == 'qt_ing_cornd' ~ 'Not\nDefined')) %>%
  ggplot(aes(x = data, 
             y = value,
             linetype = var,
             fill = var)) +
  geom_bar(stat = 'identity', position = 'dodge') +
  geom_smooth(method = 'loess', col = 'black') +
  ggthemes::scale_fill_ptol() +
  ggthemes::theme_clean() +
  facet_grid(var ~ no_regiao, scales = 'free_y') +
  guides(fill = F,
         linetype = F) +
  xlab('') +
  scale_y_continuous(name = 'Number of students admitted', 
                     labels = scales::label_number(suffix = "K", scale = 1e-3, big.mark = '.')) +
  ggtitle('Number of Students by ethnic group and Region')

students %>%
  mutate(data = zoo::as.yearmon(nu_ano_censo)) %>%
  group_by(data) %>%
  summarise_at(c("qt_ing_0_17", "qt_ing_18_24", "qt_ing_25_29", "qt_ing_30_34",
                 "qt_ing_35_39", "qt_ing_40_49", "qt_ing_50_59", "qt_ing_60_mais"), sum) %>%
  gather('var', 'value', 2:9) %>%
  mutate(var = case_when(var == 'qt_ing_0_17' ~ '<18',
                         var == 'qt_ing_18_24' ~ '18~24',
                         var == 'qt_ing_25_29' ~ '25~29',
                         var == 'qt_ing_30_34' ~ '30~34', 
                         var == 'qt_ing_35_39' ~ '35~39', 
                         var == 'qt_ing_40_49' ~ '40~49',
                         var == 'qt_ing_50_59' ~ '50~59',
                         var == 'qt_ing_60_mais' ~ '>60'),
         var = fct_relevel(as.factor(var), '<18', '18~24', '25~29',
                           '30~34', '35~39', '40~49', 
                           '50~59', '>60')) %>%
  ggplot(aes(x = data, 
             y = value,
             linetype = var,
             fill = var)) +
  geom_bar(stat = 'identity', position = 'dodge') +
  geom_smooth(method = 'loess', col = 'black') +
  ggthemes::scale_fill_ptol() +
  ggthemes::theme_clean() +
  facet_wrap( ~ var, scales = 'free_y') +
  guides(fill = F,
         linetype = F) +
  xlab('') +
  scale_y_continuous(name = 'Number of students admitted', 
                     labels = scales::label_number(suffix = "K", scale = 1e-3, big.mark = '.')) +
  ggtitle('Number of Students by age')

students %>%
  mutate(data = zoo::as.yearmon(nu_ano_censo)) %>%
  group_by(data, no_regiao) %>%
  summarise_at(c("qt_ing_0_17", "qt_ing_18_24", "qt_ing_25_29", "qt_ing_30_34",
                 "qt_ing_35_39", "qt_ing_40_49", "qt_ing_50_59", "qt_ing_60_mais"), sum) %>%
  gather('var', 'value', 3:10) %>%
  mutate(var = case_when(var == 'qt_ing_0_17' ~ '<18',
                         var == 'qt_ing_18_24' ~ '18~24',
                         var == 'qt_ing_25_29' ~ '25~29',
                         var == 'qt_ing_30_34' ~ '30~34', 
                         var == 'qt_ing_35_39' ~ '35~39', 
                         var == 'qt_ing_40_49' ~ '40~49',
                         var == 'qt_ing_50_59' ~ '50~59',
                         var == 'qt_ing_60_mais' ~ '>60'),
         var = fct_relevel(as.factor(var), '<18', '18~24', '25~29',
                           '30~34', '35~39', '40~49', 
                           '50~59', '>60')) %>%
  ggplot(aes(x = data, 
             y = value,
             linetype = var,
             fill = var)) +
  geom_bar(stat = 'identity', position = 'dodge') +
  geom_smooth(method = 'loess', col = 'black') +
  ggthemes::scale_fill_ptol() +
  ggthemes::theme_clean() +
  facet_grid(var ~ no_regiao, scales = 'free_y') +
  guides(fill = F,
         linetype = F) +
  xlab('') +
  scale_y_continuous(name = 'Number of students admitted', 
                     labels = scales::label_number(suffix = "K", scale = 1e-3, big.mark = '.')) +
  ggtitle('Number of Students by age')

stats_df <- students %>%
  mutate(data = zoo::as.yearmon(nu_ano_censo)) %>%
  group_by(data, no_regiao, no_uf) %>%
  summarise_if(is.numeric, sum) %>%
  ungroup() %>%
  select(data, no_regiao, no_uf, qt_ing_fem:qt_ing_noturno, qt_ing_0_17:qt_ing_cornd, qt_ing_financ:qt_ing_apoio_social) %>%
  gather('var', 'value', 4:42) %>%
  group_by(data, no_regiao, no_uf, var) %>%
  summarise_if(is.numeric, sum) %>%
  ungroup()

stats_df %>%
  mutate(value = ifelse(value == 0, 1, value)) %>%
  lm(log(value) ~ as.factor(var) + data + as.factor(no_regiao) + as.factor(no_uf), data = .) %>%
  sjPlot::tab_model(show.ci = F, show.est = T, show.se = T, show.r2 = T, show.stat = T,
                    pred.labels = labels_table_model,
                    dv.labels = 'Relationship between number of students admitted')

cluster_regiao_sex <- stats_df %>%
  pivot_wider(names_from = var, values_from = value) %>%
  group_by(no_regiao) %>%
  summarise_if(is.numeric, sum) %>%
  select(no_regiao, Male = qt_ing_masc, Female = qt_ing_fem)

names_cluster_regiao_sex <- cluster_regiao_sex$no_regiao

cluster_regiao_sex <- scale(cluster_regiao_sex[,2:3])

rownames(cluster_regiao_sex) <- names_cluster_regiao_sex

clusGap(cluster_regiao_sex, FUN = kmeans, K.max = 4, B = 1000)

fviz_cluster(kmeans(cluster_regiao_sex, 4, nstart = 25), cluster_regiao_sex[,1:2], ellipse.type = "norm", 
             repel = T, ggtheme = ggthemes::theme_clean(), palette = c('#4477AA', '#CC6677', '#117733', '#AA4499'),
             main = 'Grouping by Regions and Sex')

cluster_uf_sex <- stats_df %>%
  pivot_wider(names_from = var, values_from = value) %>%
  group_by(no_uf) %>%
  summarise_if(is.numeric, sum) %>%
  select(no_uf, Male = qt_ing_masc, Female = qt_ing_fem) %>%
  ungroup()

names_cluster_uf_sex <- cluster_uf_sex$no_uf

cluster_uf_sex <- scale(cluster_uf_sex[,2:3])

rownames(cluster_uf_sex) <- names_cluster_uf_sex

clusGap(cluster_uf_sex, FUN = kmeans, K.max = 10, B = 1000)

fviz_cluster(kmeans(cluster_uf_sex, 7, nstart = 25), cluster_uf_sex, ellipse.type = "norm", 
             repel = T, ggtheme = ggthemes::theme_clean(), palette = c('#4477AA', '#CC6677', '#117733', '#DDCC77', '#AA4499', '#661100', '#44AA99'),
             main = 'Grouping by States and Sex')

cluster_regiao_ethnic <- stats_df %>%
  pivot_wider(names_from = var, values_from = value) %>%
  group_by(no_regiao) %>%
  summarise_if(is.numeric, sum) %>%
  select(no_regiao, 
         White = qt_ing_branca, 
         Black = qt_ing_preta,
         Brown = qt_ing_parda,
         Asian = qt_ing_amarela,
         Not_declared = qt_ing_cornd) %>%
  ungroup()

names_cluster_regiao_ethnic <- cluster_regiao_ethnic$no_regiao

cluster_regiao_ethnic <- scale(cluster_regiao_ethnic[,2:6])

rownames(cluster_regiao_ethnic) <- names_cluster_regiao_ethnic

clusGap(cluster_regiao_ethnic, FUN = kmeans, K.max = 4, B = 1000)

fviz_cluster(kmeans(cluster_regiao_ethnic, 4, nstart = 25), cluster_regiao_ethnic, ellipse.type = "norm", 
             repel = T, ggtheme = ggthemes::theme_clean(), palette = c('#4477AA', '#CC6677', '#117733', '#AA4499'),
             main = 'Grouping by Regions and Ethnic group')

cluster_uf_ethnic <- stats_df %>%
  pivot_wider(names_from = var, values_from = value) %>%
  group_by(no_uf) %>%
  summarise_if(is.numeric, sum) %>%
  select(no_uf, 
         White = qt_ing_branca, 
         Black = qt_ing_preta,
         Brown = qt_ing_parda,
         Asian = qt_ing_amarela,
         Not_declared = qt_ing_cornd) %>%
  ungroup()

names_cluster_uf_ethnic <- cluster_uf_ethnic$no_uf

cluster_uf_ethnic <- scale(cluster_uf_ethnic[,2:6])

rownames(cluster_uf_ethnic) <- names_cluster_uf_ethnic

clusGap(cluster_uf_ethnic, FUN = kmeans, K.max = 10, B = 1000)

fviz_cluster(kmeans(cluster_uf_ethnic, 9, nstart = 25), cluster_uf_ethnic, ellipse.type = "norm", 
             repel = T, ggtheme = ggthemes::theme_clean(), palette = c('#4477AA', '#CC6677', '#117733', '#DDCC77', '#AA4499', '#661100', '#44AA99', '#332288', '#AA4466'),
             main = 'Grouping by States and Ethnic group')