Skip to content

Instantly share code, notes, and snippets.

@arthurwelle

arthurwelle/PNADCT_NovosPesos.Rmd

Last active December 6, 2021 01:06
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save arthurwelle/f036300efa70f3329b6d45ff7adfc3f7 to your computer and use it in GitHub Desktop.
Save arthurwelle/f036300efa70f3329b6d45ff7adfc3f7 to your computer and use it in GitHub Desktop.
Download ZIP
Raw
PNADCT_NovosPesos.Rmd
---
title: "Análise PNADCT"
output:
html_document:
toc: no
---
# Análise PNADC Trimestral
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)
library(PNADcIBGE)
library(data.table)
library(survey)
gc()
```
Quantos homens e mulheres no Brasil no 3T2021?
Sidra diz: https://sidra.ibge.gov.br/tabela/5917#resultado (site atualizado em 30/11/2021)
homens 104020 x 10e3
mulheres 108788 x 10e3
## Leitura base 3T2021 a partir dos microdados TXT disponibilizados pelo IBGE.
```{r}
diretorio_dados_txt <- "C:/Users/arthur/Desktop/PNADCT/txt/"
diretorio_input_sas <- "C:/Users/arthur/Desktop/PNADCT/dicionarios/"
# leitura
d <- PNADcIBGE::read_pnadc(
microdata = paste0(diretorio_dados_txt, "PNADC_032021.txt"),
input_txt = paste0(diretorio_input_sas, "input_PNADC_trimestral.txt"))
# transforma em data.table
d <- data.table::as.data.table(d)
```
## 1. Data.table
```{r}
# quantidade por sexo
d[,.(t = sum(V1028)), by = .(V2007)]
```
2 108787836
1 104020393
Soma simples de pesos usando data.table (poderia ser feito com qualquer outro método, ColSums, tidyverse etc).
O resultado bate com o Sidra.
## 2. Método survey antigo V1027 (ERRADO!!)
Usando peso sem pós-estratificação (V1027) em conjunção com a função postStratify.
Assim era feito antes e assim retornava os resultados corretos.
```{r}
# preliminary survey design
pre_stratified <-
survey::svydesign(
ids = ~ UPA ,
strata = ~ Estrato ,
weights = ~ V1027 ,
data = d ,
nest = TRUE
)
# post-stratification targets
df_pos <- data.frame( posest = unique( d$posest ) , Freq = unique( d$V1029 ) ) # V1029 projeção da população
# final survey design object
p1 <- survey::postStratify( pre_stratified , ~ posest , df_pos )
result <- survey::svytotal(x = ~as.factor(V2007),
design = p1,
na.rm=TRUE,
level=0.95)
result
```
total SE
as.factor(V2007)1 101205560 180824
as.factor(V2007)2 111602669 180824
Porém agora os resultados não mais batem.
A correção dos pesos ocorreu somente na variável V1028. Ela deveria ter sido feita também na variável V1027? Não a V1027 é antes de pos-estratificação e antes de calibração. Então ela é um instrumento basico para ser usado quando quisermos aplicar outros métodos.
## 3. Se trocarmos V1027 por V1028 (alterando o método antigo).
Pra mim não deveria, mas o resultado bate.
```{r}
# preliminary survey design
pre_stratified <-
survey::svydesign(
ids = ~ UPA ,
strata = ~ Estrato ,
weights = ~ V1028 , # aqui era pra ser V1027, mas só funcionou pra mim com o peso ERRADO com pos-estratificação
data = d ,
nest = TRUE
)
# post-stratification targets
df_pos <- data.frame( posest = unique( d$posest ) , Freq = unique( d$V1029 ) ) # V1029 projeção da população
# final survey design object
p1 <- survey::postStratify( pre_stratified , ~ posest , df_pos )
result <- survey::svytotal(x = ~as.factor(V2007),
design = p1,
na.rm=TRUE,
level=0.95)
result
```
total SE
as.factor(V2007)1 104020393 184545
as.factor(V2007)2 108787836 184545
É o mesmo que usar direto o V1028 com srvyr::as_survey_design
```{r}
#install.packages("srvyr")
p2 <- srvyr::as_survey_design(
ids = UPA,
strata = Estrato,
weights = V1028,
.data = d ,
nest = TRUE
)
result <- survey::svytotal(x = ~as.factor(V2007),
design = p2,
na.rm=TRUE,
level=0.95)
result
```
total SE
as.factor(V2007)1 104020393 595444
as.factor(V2007)2 108787836 597805
Resultado correto, com alterações no SE.
## 4. Usando os pesos de replicação
Usando exemplo de código de Guilherme Jacob (veja https://guilhermejacob.github.io/2021/12/pnadc-raking-bootstrap/)
```{r}
# coleta matriz de pesos
bootweights.cols <- grep( "V1028[0-9]{3}" , colnames( d ) , value = TRUE )
bootweights.mat <- d[ , ..bootweights.cols ]
bootweights.mat <- as.matrix( bootweights.mat )
d[ , bootweights.cols ] <- NULL
# declarando o plano amostral
pnadc.boot <- svrepdesign( data = d , type = "bootstrap" , weights = ~V1028 , repweights = bootweights.mat )
pnadc.boot_result <- survey::svytotal(x = ~as.factor(V2007),
design = pnadc.boot,
na.rm=TRUE,
level=0.95)
pnadc.boot_result
```
total SE
as.factor(V2007)1 104020393 0.0987
as.factor(V2007)2 108787836 0.0797
Resultado é o mesmo, correto, porém SE se torna muito baixa (porém não zero) pois a calibração foi feita com sexo mas se trata de boostrap.
## 5. Método dos Conglomerados Primários
Novamente, usando exemplo de código de Guilherme Jacob (veja https://guilhermejacob.github.io/2021/12/pnadc-raking-bootstrap/)
```{r}
# ajusta formatos
d[ , c( "posest" , "posest_sxi" ) ] <- lapply( d[ , c( "posest" , "posest_sxi" ) ] , factor )
# declarando o plano amostral
pnadc.ucm <- svydesign( ids = ~UPA+V1008 , strata = ~Estrato , data = d , weights = ~V1028 )
# coleta tabela com os totais das marginais na população
pop.posest <- d[ !duplicated( d$posest ) , c( "posest" , "V1029" ) ]
pop.posest_sxi <- d[ !duplicated( d$posest_sxi ) , c( "posest_sxi" , "V1033" ) ]
colnames( pop.posest )[2] <-"Freq"
colnames( pop.posest_sxi )[2] <-"Freq"
# aplica calibração via raking com limites e pesos fixos no domicílio
pnadc.ucm <- calibrate( pnadc.ucm ,
formula = list( ~posest , ~posest_sxi ), population = list( pop.posest, pop.posest_sxi ) ,
bounds = c(.2,5) , bounds.const = FALSE ,
calfun = "raking" , aggregate.stage = 2 ,
sparse = TRUE )
pnadc.ucm_result <- survey::svytotal(x = ~as.factor(V2007),
design = pnadc.ucm,
na.rm=TRUE,
level=0.95)
pnadc.ucm_result
```
total SE
as.factor(V2007)1 104020393 0
as.factor(V2007)2 108787836 0
Resultado é o mesmo, correto, porém SE vai à zero o que é esperado visto que a calibração foi feita justamente com as variáveis de sexo.
Adicionalmente é muito mais demorado para retornar o resultado.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment