4

O governo anunciou ano passado o resultado da PNAD 2013, nele aparece que há 85,3% das casa tinham água encanada.

Baixei os microdados da PNAD e com o R tentei reproduzir a análise segue abaixo o código:

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 que​m tem​​ agua encanada (0.853)
 # mas retorna 0.83964
 svymean(~V4624, design = sample.pnad, na.rm=TRUE)

Com isso obtenho o resultado de 84%. 1,3% a menos que o governo.

Essa é a primeira vez que uso o pacote survey, então queria saber o que estou fazendo de errado.

6
  • Não sei se é o caso, mas o dicionário da PNAD tem um problema: duas variáveis compartilham caracteres do banco de dados (a UF e Cidade). Confere se no teu arquivo não tem isso. Se tiver, apaga a linha do UF (tem que ser essa) do dicionário e tenta de novo!
    – Rcoster
    26/09/2015 às 14:20
  • Pequena correção: É UF e Número de Controle (V0102), e não UF e Município
    – Rcoster
    26/09/2015 às 14:25
  • 1
    Rodei o teu código e não é esse o problema... talvez nem haja problema: no próprio IBGE diz que o percentual é 83,96% (goo.gl/qAzi1X)
    – Rcoster
    26/09/2015 às 14:49
  • Muito obrigado @Rcoster. Você tem alguma ideia porque esse número diverge do dados do governo?
    – celacanto
    26/09/2015 às 15:13
  • 1
    Sinalizei a pergunta para ser fechada pois não há problema no código, o resultado bate com a referência oficial. O erro está nas outras fontes, não há dúvida sobre programação aqui.
    – Molx
    28/09/2015 às 23:05

0

Pesquise outras perguntas com a tag .