In tomasbarcellos/pof: Leitura da Pesquisa de Orçamentos Familiares (POF)

knitr::opts_chunk$set(
  collapse = TRUE, warning = FALSE, message = FALSE,
  comment = "#>"
)

library(pof)
library(tidyverse)

setwd("..")
# v9001 é código do rendimento
produtos <- readxl::read_excel("dados/2018/Cadastro de Produtos.xls", 
                               col_names = c("quadro", "v9001",
                                             "descricao_produto"), 
                               skip = 1) 

# Rendimento - Nível rendimento
rendimentos_trabalho <- pof::ler_rend_trab(2018) %>% 
  janitor::clean_names() %>% 
  mutate(
    cod_uc = paste0(cod_upa, num_dom, num_uc),
    cod_pessoa = paste0(cod_uc, cod_informante),
    cod_rendimento = paste0(cod_pessoa, quadro, sub_quadro, seq),
  )

# Rendimento - Nível rendimento
rendimentos_outros <- pof::ler_rend_outros(2018) %>% 
  janitor::clean_names() %>% 
  mutate(
    cod_uc = paste0(cod_upa, num_dom, num_uc),
    cod_pessoa = paste0(cod_uc, cod_informante),
    cod_rendimento = paste0(cod_pessoa, quadro, seq),
  ) 

rendas_uc_trab <- rendimentos_trabalho %>% 
  filter(v8500 != 9999999) %>% 
  group_by(cod_uc, v9001 = as.character(v9001)) %>% 
  summarise(rendimento = sum(v8500 * v9011, na.rm = TRUE)/12, 
            peso_final = first(peso_final),
            renda_total = first(renda_total)) 

rendas_uc_outros <- rendimentos_outros %>% 
  filter(v8500 != 9999999) %>% 
  group_by(cod_uc, v9001 = as.character(v9001)) %>% 
  summarise(rendimento = sum(v8500 * v9011, na.rm = TRUE)/12, 
            peso_final = first(peso_final),
            renda_total = first(renda_total)) 

rendas_uc <- rendas_uc_outros %>% 
  bind_rows(rendas_uc_trab) %>%
  arrange(cod_uc, desc(rendimento)) %>% 
  left_join(produtos) %>% 
  mutate(partic = round(rendimento * 100 / renda_total)) %>% 
  ungroup()

# Daqui para baixo nessa chunck é tudo furada

rendas_uc %>% 
  group_by(cod_uc) %>% 
  summarise(rendimento = max(rendimento),
            rendimento2 = cut(rendimento, c(-Inf, 1000, 2000, 5000, 10000, Inf), c("Mis", "Massa", "B", "M", "A"))) %>% 
  count(rendimento2) %>% 
  mutate(prop = n / sum(n))

rendimentos <- rendimentos_trabalho %>% 
  bind_rows(rendimentos_outros, .id = "tipo") %>% 
  filter(v8500 != 9999999) %>% 
  mutate(tipo = c("trabalho", "outros")[as.integer(tipo)]) %>% 
  group_by(cod_uc, tipo) %>% 
  summarise(rendimento = sum(v8500, na.rm = TRUE), 
            peso_final = first(peso_final),
            renda_total = first(renda_total)) 

unidades <- rendimentos %>% 
  select(-renda_total) %>% 
  ungroup() %>% 
  pivot_wider(names_from = tipo, values_from = rendimento, 
              values_fill = list(rendimento = 0)) %>% 
  mutate(total = outros + trabalho,
         p_outros = outros/total,
         p_trabalho = trabalho/total)


unidades %>% 
  arrange(desc(p_trabalho)) %>% 
  group_by(classe = ifelse(p_trabalho <= 0.1, "baixa", "alta")) %>%
  summarise_if(is.numeric, mean)

grupos <- kmeans(unidades$p_trabalho, c(0.1, 0.9))
grupos$cluster %>% table()
grupos$centers

grupos2 <- unidades %>% 
  select(p_trabalho, trabalho) %>% 
  kmeans(2)

grupos2$cluster %>% table()
grupos2$centers
grupos2$totss

unidades %>% 
  mutate(grupo1 = as.factor(grupos$cluster),
         grupo2 = as.factor(grupos2$cluster)) %>% 
  filter(total < 100000) %>% 
  ggplot(aes(p_trabalho * 100, total, col = grupo1)) + 
  geom_point(alpha = 0.3, size = 0.3) + 
  geom_vline(xintercept = c(80, 90, 95), lty = 1:3) +
  theme_classic()

ggplot(unidades, aes(100 * p_outros)) + 
  geom_density() + 
  geom_vline(xintercept = c(80, 90, 95), lty = 1:3) +
  labs(title = "Percentual da renda de outros", x = "")

ggplot(unidades %>% filter(total < 100000), 
       aes(100 * p_trabalho, total)) + 
  geom_point(alpha = 0.5, shape = ".") + 
  geom_vline(xintercept = c(80, 90, 95), lty = 1:3) +
  labs(title = "Percentual da renda do trabalho x renda total", 
       x = "", y = "renda total") + 
  theme_classic()

ggplot(unidades %>% filter(total < 100000), 
       aes(100 * p_trabalho, trabalho)) + 
  geom_point(alpha = 0.5, shape = ".") + 
  geom_vline(xintercept = c(80, 90, 95), lty = 1:3) +
  labs(title = "Percentual da renda do trabalho x renda trabalho",
       x = "", y = "renda do trabalho") + 
  theme_classic()

ggplot(unidades %>% filter(total < 100000), 
       aes(100 * p_trabalho, outros)) + 
  geom_point(alpha = 0.5, shape = ".") + 
  geom_vline(xintercept = c(80, 90, 95), lty = 1:3) +
  labs(title = "Percentual da renda do trabalho x renda outros",
       x = "", y = "renda do outros") + 
  theme_classic()

categs_outros <- rendas_uc_outros %>% 
  ungroup() %>% 
  # count(v9001, sort = T) %>% 
  # select(v9001) %>%
  group_by(v9001) %>% 
  summarise(total = sum(rendimento)) %>%
  arrange(desc(total)) %>% 
  mutate(tipo = "outros") %>% 
  left_join(produtos)

categs_trab <- rendas_uc_trab %>% 
  ungroup() %>% 
  # count(v9001, sort = T) %>% 
  # select(v9001) %>%
  group_by(v9001) %>% 
  summarise(total = sum(rendimento)) %>%
  arrange(desc(total)) %>% 
  mutate(tipo = "trabalho") %>% 
  left_join(produtos)


dicionario <- read_csv2("codigos.csv", na = "-") %>% 
  mutate(v9001 = as.character(v9001))

rendas_cv_mv <- rendas_uc_outros %>% 
  bind_rows(rendas_uc_trab) %>% 
  # left_join(dicionario, by = "v9001")
  left_join(dicionario %>% mutate(tipo = c(tipo[1:37], "cv", tipo[39:41])))

rendas_cv_mv %>% 
  filter(str_detect(descricao_produto, "VENDA DE|POUP|EMPRESTIMO", negate = TRUE)) %>% 
  group_by(tipo) %>% 
  summarise(qtd = n(),
            peso = sum(peso_final),
            massa = sum(rendimento),
            media = weighted.mean(rendimento, peso_final)) %>% 
  mutate(massa_peso = massa * (peso  / sum(peso)),
         media_peso = massa / peso)


unidade_cv_mv <- rendas_cv_mv %>% 
  mutate(tipo = replace_na(tipo, "x")) %>% 
  group_by(cod_uc, tipo, peso_final) %>% 
  summarise(rendimento = sum(rendimento)) %>% 
  pivot_wider(names_from = tipo, values_from = rendimento, 
              values_fill = list(rendimento = 0)) %>% 
  mutate(total = cv + mv + x,
         tx_cv = cv / total) %>% 
  ungroup()

grupos3 <- unidade_cv_mv %>% 
  select(total, tx_cv) %>% 
  mutate_all(scale) %>% 
  kmeans(2)


unidade_cv_mv %>% 
  ggplot(aes(tx_cv)) + 
  geom_density() + 
  geom_vline(xintercept = c(0.1, 0.15, 0.2), lty = 2:4)

unidade_cv_mv %>%
  mutate(grupo = grupos3$cluster) %>% 
  ggplot(aes(cv, total)) +
  geom_point(size = 0.3, alpha = 0.1, aes(col = factor(grupo))) + 
  scale_y_log10() +
  scale_x_log10() +
  theme_classic()

unidade_cv_mv %>%
  mutate(grupo = grupos3$cluster,
         grupo2 = tx_cv > 0.15) %>% 
  ggplot(aes(factor(grupo2), mv, fill = factor(grupo2))) +
  geom_boxplot() + 
  scale_y_log10()

rendas_uc_outros %>% 
  summarise(rendimento = max(rendimento),
            rendimento2 = cut(rendimento,
                              c(0, 1000, 2000, 5000, 10000, Inf), 
                              c("Mis", "Massa", "B", "M", "A")),
            peso = sum(peso_final)) %>% 
  group_by(rendimento2) %>% 
  summarise(n = n(), peso = sum(peso)) %>% 
  mutate(prop = n / sum(n), prop2 = peso/sum(peso))

# cad_coletiva_uc <- cad_coletiva %>% 
#   group_by(cod_uc) %>% 
#   select(starts_with("v"), renda_total) %>% 
#   summarise_all(sum, na.rm = TRUE)

usando tradutor

# Tudo errado.
# trabalho em progresso
tradutor_2018 <- ler_tradutor_rendimento(2018) %>% 
  mutate(descricao_3 = ifelse(is.na(descricao_3), descricao_2, descricao_3)) %>% 
  select(codigo, nivel_3, tipo_renda = descricao_3)

# Tentativa de reproduzir dados UBGE
peso_uc <- rendas_uc %>% 
  select(cod_uc, peso_final) %>% 
  unique()

rendas_uc %>% 
  mutate(codigo = str_sub(v9001, 1, 5)) %>% 
  filter(str_sub(cod_uc, 1, 1) == "2") %>% 
  left_join(tradutor_2018, "codigo") %>% 
  select(cod_uc, rendimento, nivel_3, tipo_renda) %>% 
  group_by(cod_uc, nivel_3, tipo_renda) %>% 
  summarise(rendimento = sum(rendimento, na.rm = TRUE)) %>% 
  left_join(peso_uc, "cod_uc") %>% 
  ungroup() %>% 
  # group_by(nivel_3, tipo_renda) %>% 
  summarise(media = collapse::fmean(rendimento, w = peso_final),
            media2 = mean(rendimento),
            cv = collapse::fsd(rendimento, w = peso_final)/media,
            cv2 = sd(rendimento)/mean(rendimento),
            soma = collapse::fsum(rendimento * peso_final)/1e9
            ) 






rendimento_prof <- rendimentos_trabalho %>% 
  bind_rows(rendimentos_outros) %>% 
  mutate(codigo = str_sub(v9001, 1, 5)) %>% 
  left_join(tradutor_2018, by = "codigo") %>% 
  group_by(cod_uc, cod_pessoa, tipo_renda, v5302, v5303, v5304) %>% 
  summarise(rendimento = sum(v8500, na.rm = TRUE),
            peso_final = first(peso_final)) %>% 
  ungroup()

rend_categs <- rendimento_prof %>% 
  mutate(forma = case_when(
    v5302 == 1 ~ "cv",
    v5302 == 2 ~ "mv",
    v5302 == 3 ~ "cv",
    v5302 == 4 & v5303 == 1 ~ "mv", # setor estatutário
    # supoe não ter informalidade no setor publico
    v5302 == 4 & v5304 == 1 ~ "cv", # trabalhador de estatal
    # Tem casos (2491) de servidor sem carteria assinada ou estatuto. Pq?
    v5302 == 5 ~ "mv", # empregador
    v5302 == 6 ~ "cv/mv", # conta própria # não encontrei setor trabalhado
    v5302 == 7 ~ "cv", # trab não remunerado ?
    TRUE ~ NA_character_,
  ))

rend_categs %>% 
  group_by(forma) %>% 
  summarise(rendimento = sum(rendimento))

todos_rendimentos <- rendimentos_trabalho %>% 
  bind_rows(rendimentos_outros)

todos_rendimentos %>% 
  mutate(codigo = str_sub(v9001, 1, 5)) %>% 
  left_join(tradutor_2018, by = "codigo") %>% 
  filter(v8500 != 9999999) %>%
  group_by(cod_uc, tipo_renda) %>%
  summarise(v8500 = sum(v8500, na.rm = TRUE), 
            peso_final = first(peso_final)) %>% 
  # ungroup() %>% 
  group_by(tipo_renda) %>%
  summarise(rendimento = weighted.mean(v8500, peso_final, na.rm = TRUE))
  # summarise(rendimento = sum(v8500, na.rm = TRUE)) %>% 
  # mutate(prop = rendimento * 100 / sum(rendimento))