The government announced last year the result of the 2013 PNAD, it appears that there is #
I downloaded the PNAD microdata and with the R attempt to reproduce the analysis below the code:
setwd( "/saneamento_basico/")
# pacotes
# -------------------------------------------------------------------------------
library(xlsx)
library(descr)
library(survey)
# paths
# -------------------------------------------------------------------------------
dicionarioPath <- "PNAD/Dicionаrios e input/Dicionаrio de variаveis de domicбlios - PNAD 2013.xls"
pnadPath <- "PNAD/Dados/DOM2013.txt"
pnadCsvOut <- "PNAD/DOM2013.txt"
# traduz arquivo
# -------------------------------------------------------------------------------
# dicionario
dicionario <- read.xlsx(dicionarioPath, sheetIndex = 1,
stringsAsFactors = FALSE,
startRow = 4)[,1:3]
dicionario <- dicionario[complete.cases(dicionario),]
colnames(dicionario) <- c("start", "size", "var")
dicionario$end <- as.numeric(dicionario$start) + dicionario$size - 1
# traducao
fwf2csv(fwffile = pnadPath,
csvfile = pnadCsvOut,
names = dicionario$var,
begin = dicionario$start,
end = dicionario$end)
# le e formata pnad
# -------------------------------------------------------------------------------
pnad <- read.csv(pnadCsvOut, sep = '\t')
# "Forma de abastecimento de água"
coluna_formaDeAbastecimentoDeAgua <- "V4624"
colunas_pesoDomicilio <- "V4611" # "Peso do domicílio"
strat <- "V4602" # "Estrato"
id <- "V4618" # "PSU - Unidade primária de amostragem"
# retira quem tem NA nos pesos
pnad <- pnad[!is.na(pnad[,colunas_pesoDomicilio]),]
# retira quem tem NA em forma de abastecimento
pnad <- pnad[!is.na(pnad[,coluna_formaDeAbastecimentoDeAgua]),]
# tranforma quem tem agua encanada em 1 e os outros em 0
temAguaEncanada <- as.numeric(pnad[,coluna_formaDeAbastecimentoDeAgua] == 1)
pnad[,coluna_formaDeAbastecimentoDeAgua] <- temAguaEncanada
sample.pnad <- svydesign(ids = ~V4618,
strata = ~V4602,
weights = ~V4611,
data = pnad,
nest = TRUE)
# para evitar erro "has only one PSU at stage 1"
options(survey.lonely.psu = "adjust")
# deveria retornar a porcentagem de quem tem agua encanada (0.853)
# mas retorna 0.83964
svymean(~V4624, design = sample.pnad, na.rm=TRUE)
With this I get the result of 84%. 1.3% less than the government.
This is the first time I use the survey package, so I wanted to know what I'm doing wrong.