## Exemplos básicos de uso do pacote PNADcIBGE ---------------------------------
## Autor: Prof. Wagner Hugo Bonat · Ômega Data Academy -------------------------
## Curso para o IPARDES --------------------------------------------------------
rm(list = ls())
gc(reset = TRUE)

## Vamos dar uma olhada na documentação
# https://cran.r-project.org/web/packages/PNADcIBGE/refman/PNADcIBGE.html#get_pnadc

## Certifique-se de estar trabalhando no local CORRETO
setwd("~/GitProjects/Curso_R_IPARDES/Scripts")

## Instalação dos pacotes
# install.packages("survey")
# install.packages("PNADcIBGE")

## Carregando os pacotes adicionais
library(survey)
library(PNADcIBGE)
library(tidyverse)
library(scales)

## Carregando funções auxiliares
source('functions.R')

## Download dos dados e organização dos dados
## Crie um pasta onde os dados vão ser baixados
diretorio <- "/home/wagner/GitProjects/Curso_R_IPARDES/Scripts/Dados_PNAD"

## Passo 1: Baixar os dados e o delineamento amostral
# Vai fazer o download e gravar na pasta diretorio (apenas na primeira vez)
# Vou usar o ano de 2025 o quarto trimestre.

dadosPNADc <- get_pnadc(year = 2025, quarter = 4, 
                        reload = FALSE,
                        savedir = diretorio)
## Use o argumento reload = FALSE se já tiver feito o download dos dados


## Diretório ftp do IBGE
# https://ftp.ibge.gov.br/Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_continua/Trimestral/Microdados/
# data_pnadc <- PNADcIBGE::read_pnadc(microdata = microdatafile, input_txt = inputfile, vars = vars)
# Exemplo básico de delineamento de plano amostral
# data_prior <- survey::svydesign(ids = ~UPA + ID_DOMICILIO, strata = ~Estrato, 
# data = data_pnadc, weights = ~V1035, nest = TRUE)

## Classe do objeto
class(dadosPNADc)

## Variáveis que foram extraídas
View(dadosPNADc$variables)
names(dadosPNADc$variables)
table(dadosPNADc$variables$UF)


## Mostrar o diretório e quais informações foram baixadas

################################################################################
## Estimando totais ############################################################
################################################################################

## Podemos selecionar subpopulações conforme seja de interesse.
## Vou focar apenas na UF == PR
PR <- subset(dadosPNADc, UF == 'Paraná') 
View(PR$variables)

## Exemplo para variável numérica
## Consultar o dicionário é fundamental!
## VD4020 - Rendimento mensal efetivo de todos os 
# trabalhadores para pessoas de 14 anos ou mais

total_renda <- svytotal(x= ~VD4020, design = PR, na.rm = TRUE)
total_renda
cv(object = total_renda)
confint(object = total_renda, level = 0.99)

## Limpando a área
gc(reset = TRUE)

## Exemplo variável categórica
## V2007 - Sexo
total_sexo <- svytotal(x = ~V2007, design = PR, na.rm = TRUE)
total_sexo
confint(object = total_sexo, level = 0.99)

## Categórica com vários níveis
total_raca <- svytotal(x = ~V2010, design = PR, na.rm = TRUE)
total_raca
confint(object = total_raca)

## Tabelas de dupla entrada
## Sexo vs Raça
sexo_raca <- svytotal(x = ~interaction(V2007, V2010), 
                      design = PR, na.rm = TRUE)
sexo_raca
organizar_interaction(as.data.frame(sexo_raca), nomes_variaveis = c('sexo', 'raça'))

## Sexo vs Urbana e Rural
sexo_urb <- svytotal(x = ~interaction(V2007, V1022), 
                      design = PR, na.rm = TRUE)
sexo_urb
organizar_interaction(as.data.frame(sexo_Urb), nomes_variaveis = c('sexo', 'urbana'))

## Tabelas de tripla entrada. Cuidado que vai ficar caro calcular!
sexo_raca_urb <- svytotal(x = ~interaction(V2007, V2010, V1022), 
                         design = dadosPNADc, na.rm = TRUE)
sexo_raca_urb
organizar_interaction(as.data.frame(sexo_raca_urb), nomes_variaveis = c('sexo', 'raça','urbana'))

## Interação entre numérica e categórica tem que usar outra estrutura

################################################################################
## Estimando médias ############################################################
################################################################################
gc(reset = TRUE)
# Renda média
# Cuidado com o deflacionador! Consultar a documentação. 
renda_media <- svymean(x = ~VD4020, design = PR, na.rm = TRUE)
renda_media
confint(renda_media)

# Lembre-se que proporção é uma média
# Use de memória RAM é intenso!!!
prop_sexo <- svymean(x = ~V2007, design = PR, na.rm = TRUE)
prop_sexo

## Use o comando abaixo para tentar manter a memória RAM o mais limpa possível
gc(reset = TRUE)

## Tabela de dupla entrada em proporção
## TOMAR MUITO CUIDADO AQUI QUE O DENOMINADOR É O TOTAL
prop_sexo_raca <- svymean(x = ~interaction(V2007,V2010), 
                          design = PR, na.rm = TRUE)
tab <- organizar_interaction(as.data.frame(prop_sexo_raca), 
                             nomes_variaveis = c('Sexo', 'Raça'))
tab <- tab |>
  mutate(
    percentual = percent(mean)
  )

tab |> group_by(Sexo) |>
  summarize('Total' = sum(mean))

tab |> group_by(Raça) |>
  summarize('Total' = sum(mean))

## Depois vamos ver como calcular condicionais!

################################################################################
## Tabelas de frequência absoluta e relativa ###################################
################################################################################
svytable(~V2007, design = PR) ## Absoluto
svytable(~V2007, design = PR) |> prop.table() ## Relativa

svytable(~V4002, design = PR) ## Absoluto
svytable(~V4002, design = PR) |> prop.table() ## Absoluto

svytable(~V2007 + V4002, design = PR) ## Absoluto
svytable(~V2007 + V4002, design = PR) |> prop.table() ## Relativa

################################################################################
## Estimando razões ou taxas ###################################################
################################################################################
gc(reset = TRUE)
tx_desocupado <- svyratio(numerator= ~(VD4002=='Pessoas desocupadas'),
                          denominator = ~(VD4001=='Pessoas na força de trabalho'),
                          design = PR, na.rm = TRUE)
tx_desocupado
confint(tx_desocupado)

################################################################################
## Estimando medianas e quantis ################################################
################################################################################
gc(reset = TRUE)
mediana_renda <- svyquantile(x = ~VD4020, design = PR,
                             quantiles = 0.5, ci=TRUE, na.rm = TRUE)
mediana_renda

quantis_renda <- svyquantile(x = ~VD4020, design = PR,
                             quantiles = c(0.25, 0.5, 0.75), 
                             ci=TRUE, na.rm = TRUE)
quantis_renda

quantis_renda_mulher <- svyquantile(x = ~VD4020, 
                                    design = subset(PR, V2007=="Mulher"),
                                    quantiles = c(0.25, 0.5, 0.75), 
                                    ci=TRUE, na.rm = TRUE)
quantis_renda_mulher


quantis_renda_homem <- svyquantile(x = ~VD4020, 
                                    design = subset(dadosPNADc, V2007=="Homem"),
                                    quantiles = c(0.25, 0.5, 0.75), 
                                    ci=TRUE, na.rm = TRUE)
quantis_renda_homem


################################################################################
## Estimação com filtros adicionais ############################################
################################################################################

## Renda médio por UF
renda_sexo_raca <- svyby(formula = ~VD4020, by = ~interaction(V2007, V2010), 
                         design = PR, 
                         FUN = svymean, 
                         na.rm = TRUE)
renda_sexo_raca |> arrange(desc(VD4020))


## Podemos obter os quantis
## Renda médio por UF
renda_sexo_raca_med <- svyby(formula = ~VD4020, by = ~interaction(V2007, V2010), 
                         design = PR, 
                         FUN = svyquantile,
                         quantiles = c(0.25, 0.5, 0.75),
                         na.rm = TRUE)
renda_sexo_raca_med |> arrange(desc(VD4020.0.5))

## Comparando urbana e rural
renda_sexo_raca_urb <- svyby(formula = ~VD4020, 
                             by = ~interaction(V2007, V2010,  V1022), 
                             design = PR, 
                             FUN = svymean, 
                             na.rm = TRUE)
renda_sexo_raca_urb

renda_sexo_raca_urb_quantile <- svyby(formula = ~VD4020, 
                                      by = ~interaction(V2007, V2010,  V1022), 
                                      design = PR, 
                                      FUN = svyquantile,
                                      quantiles = c(0.25, 0.5, 0.75),
                                      na.rm = TRUE)
renda_sexo_raca_urb_quantile |> arrange(desc(VD4020.0.5))

################################################################################
## Gráficos básicos ############################################################
################################################################################

## Variáveis quantitativas (histogramas)
svyhist(~VD4020, design = PR) ## Renda geral

## Histograma suavizado (kernel)
plot(svysmooth(~VD4020, design = PR, bandwidth = 500))

## Histograma por nível de categórica
par(mfrow=c(1,2))
hist_mulher <- svyhist(~VD4020, design = subset(PR, V2007 == "Mulher")) ## Renda Mulher
hist_homem <- svyhist(~VD4020, design = subset(PR, V2007 == "Homem")) ## Renda Mulher

## Boxplot
svyboxplot(VD4020 ~ V2007, design = PR, all.outliers = TRUE)
svyboxplot(VD4020 ~ V2010, design = PR, all.outliers = TRUE)
svyboxplot(VD4020 ~ V1022, design = PR, all.outliers = TRUE)

## Distribuição acumulada
cdf_renda <- svycdf(~VD4020, design = PR)
plot(cdf_renda)

cdf_renda_mulher <- svycdf(~VD4020, design = subset(PR, V2007 == 'Mulher'))
cdf_renda_homem <- svycdf(~VD4020, design = subset(PR, V2007 == 'Homem'))
str(cdf_renda_homem)

homem <- cdf_renda_homem$VD4020(0:10000)
mulher <- cdf_renda_mulher$VD4020(0:10000)

plot(homem ~ c(0:10000), type = 'l', ylab = 'Acumulado')
lines(c(0:10000), mulher, col = 'red')

## Gráficos para relação entre duas númericas
install.packages('hexbin')
library(hexbin)
svycoplot(VD4020 ~ V2009|V2007, design = PR)
svycoplot(VD4020 ~ V2009|V2010, design = PR)
svycoplot(VD4020 ~ V2009|V1022 + V2007, design = PR)


################################################################################
## Modelos Lineares Generalizados via IPW ######################################
################################################################################

svytable(~V2007, design = subset(PR, V2009 > 25)) ## Sexo
svytable(~V2010, design = subset(PR, V2009 > 25)) ## Raça


## Fatores que influenciam a renda
## Apenas efeitos principais
modelo <- svyglm(VD4020 ~ V2009 + V2007 + V2010, 
                 design = subset(PR, V2009 > 25) )
summary(modelo)
plot(modelo)


## Interação
modelo_2 <- svyglm(VD4020 ~ V2009 + V2007*V2010, 
                 design = subset(PR, V2009 > 25) )
summary(modelo_2)
plot(modelo_2)
AIC(modelo, modelo_2)

## Interpretação é um pouco diferente da usual devido aos pesos amostrais.

# eff.p: Número efetivo de parâmetros do modelo.
# AIC: Menor AIC → melhor equilíbrio entre ajuste e complexidade.
# deltabar: É uma medida associada ao ajuste de Rao-Scott utilizado na 
# construção do AIC para amostras complexas.
# Em termos práticos: não costuma ser utilizada diretamente para escolher modelos;
# serve para calcular a penalização efetiva do AIC sob o desenho amostral;
# valores menores indicam menor efeito do desenho na correção aplicada.

## FIM -------------------------------------------------------------------------