Microdados com R: exemplo da PNADCovid de Assuncao (2021)

Abstract

This is an undergrad student level instruction for class use.

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Microdados com R: exemplo da PNADCovid de Assuncao (2021). Campo Grande-MS,Brasil: RStudio/Rpubs, 2021. Disponível em http://rpubs.com/amrofi/survey_pnadcovid.

1 Introdução

Este é um script baseado em Assunção (2021). Ele utiliza o pacote COVIDIBGE (ASSUNÇÃO; HIDALGO, 2021) entre outros.

Como exposto nos vídeos mais teóricos, a PNADC apresenta uma amostragem complexa, e portanto o usuário não pode simplesmente usar os microdados sem adotar o devido plano amostral e seus pesos. O pacote COVIDIBGE facilitará o trabalho e O srvyr apresenta algumas vantagens como o de adequar os objetos de dados ao ambiente tidyverse. O usuário deve também olhar Figueiredo (2021) onde coloco mais informações sobre o

Será interessante o usuário já ter instalados os pacotes srvyr, tidyverse, COVIDIBGE, convey e survey.

Como as rotinas levam um certo tempo para baixar da web, rodar e gerar o rmd, fiz em script e salvei os datasets em rds, depois rodei as rotinas e gerei os objetos de saídas para colocar no rmd. Recomendo o leitor que rode cada script colocado, de modo a compreender o procedimento. Para comparar os resultados, olhar o post original de Assunção (2021) em <https://rpubs.com/gabriel-assuncao-ibge/covid>.

2 INSTALAR O PACOTE COVIDIBGE

install.packages("COVIDIBGE")

3 CARREGAR O PACOTE

# PACOTE -----
library(COVIDIBGE)
help("get_covid")

4 DADOS

O pacote puxa os dados, vincula o dicionário e os deflatores e incorpora design amostral.

dadosPNADCOVID19 <- get_covid(year = 2020, month = 5)
class(dadosPNADCOVID19)
# vou salvar para não precisar chamar novamente
saveRDS(dadosPNADCOVID19, file = "dados_2020_05.rds")

Uma vez que os dados estejam salvos em rds, poderei trabalhar a partir deste, de modo offline, que é mais rápido.

# para chamar do rds, fazer:
dadosPNADCOVID_2020_05 <- readRDS("D:/Documentos/disciplinas/econometria microdados/laboratorio microdados R/pnad_covid_gabrielassuncao/dados_2020_05.rds")

É possível chamar apenas algumas variáveis de modo a reduzir o tamanho do arquivo. Ele virá como objeto survey.design, já com deflação e pós-estratificação, o que facilita o trabalho. Da mesma forma, fiz no script separado e salvei rds. O leitor pode carregar o dataset da mesma forma que o chunk anterior, com readRDS, especificando o nome do arquivo “dados_2020_05_srv.rds”. Coloquei o sufixo srv para indicar que é um objeto survey.design.

# se quiser somente algumas variaveis ----
dadosPNADCOVID19 <- get_covid(year = 2020, month = 5, vars = c("C001", "C002"))
# o comando abaixo traz as caracteristicas do plano amostral
dadosPNADCOVID19
# vendo a classe do objeto: 'survey.design2' 'survey.design'
class(dadosPNADCOVID19)  # 'survey.design2' 'survey.design'
# vou salvar para não precisar chamar novamente
saveRDS(dadosPNADCOVID19, file = "dados_2020_05_srv.rds")

5 OPÇÃO DATA.FRAME E LABELS

Neste caso, o leitor pode chamar variáveis específicas (aqui estão “C001”,“C002”), e especificar design = FALSE para retornar um objeto dataframe (classes “tbl_df”,“tbl”,“data.frame”), e depois a opção de não colocar textos das variáveis qualitativas (labels=FALSE), e desta forma, no lugar dos textos teremos os números indicativos de cada categoria. Por exemplo, onde estava UF = Rondônia, aparecerá UF = 11.

# com desing false: retorna tibble design
dadosPNADCOVID19_brutos <- get_covid(year = 2020, month = 5, vars = c("C001", "C002"),
    design = FALSE)
class(dadosPNADCOVID19_brutos)  # 'tbl_df'     'tbl'        'data.frame'
dadosPNADCOVID19_brutos
saveRDS(dadosPNADCOVID19_brutos, file = "dados_2020_05_df.rds")

# sem labels
dadosPNADCOVID19_brutos_sem <- get_covid(year = 2020, month = 5, vars = c("C001",
    "C002"), labels = FALSE, design = FALSE)
# os níveis das categorias são representados por números: UF Rondônia = 11
dadosPNADCOVID19_brutos_sem
saveRDS(dadosPNADCOVID19_brutos_sem, file = "dados_2020_05_df_sem.rds")

dadosPNADCOVID19_brutos

dadosPNADCOVID19_brutos_sem

6 IMPORTAÇÃO OFFLINE

O leitor deve primeiro baixar os arquivos de https://ftp.ibge.gov.br/Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_PNAD_COVID19/ e colocar na mesma pasta deste rmarkdown. O .csv dos Microdados estará dentro de Dados e em formato zip: https://ftp.ibge.gov.br/Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_PNAD_COVID19/Microdados/Dados/. O arquivo .xls do dicionário estará em https://ftp.ibge.gov.br/Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_PNAD_COVID19/Microdados/Documentacao/. Descompactar antes e colocar no mesmo local deste rmd. Cuidado com o nome do arquivo, pois o IBGE costuma atualizar e alterar os nomes constando a identificação da atualização. A rotina covid_labeller já vai associar os dados com o dicionário e depois com o deflator, compondo um objeto tibble.

dadosPNADCOVID19 <- read_covid(microdata = "PNAD_COVID_052020.csv")
#
dadosPNADCOVID19 <- covid_labeller(data_covid = dadosPNADCOVID19, dictionary.file = "Dicionario_PNAD_COVID_052020_20210726.xls")
#
class(dadosPNADCOVID19)  # [1] 'tbl_df'     'tbl'        'data.frame'

Agora vai associar com deflator que foi baixado de https://ftp.ibge.gov.br/Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_PNAD_COVID19/Microdados/Documentacao/ e extraído o .zip para o mesmo local do rmd. Cuidado com o nome do arquivo!

dadosPNADCOVID19 <- covid_deflator(data_covid = dadosPNADCOVID19, deflator.file = "Deflator_PNAD_COVID_2020_11.xls")
class(dadosPNADCOVID19)  # 'tbl_df'     'tbl'        'data.frame'

Agora passa para survey.design onde incorpora plano amostral. Agora estamos prontos e o comando covid_design já fez a pós-estratificação. O leitor poderá verificar o design solicitando o objeto dadosPNADCOVID19. Salvarei em .rds para uso futuro, caso deseje.

dadosPNADCOVID19 <- covid_design(data_covid = dadosPNADCOVID19)
class(dadosPNADCOVID19)  # 'survey.design2' 'survey.design' 
#
dadosPNADCOVID19
# Stratified 1 - level Cluster Sampling design (with replacement) With (14906)
# clusters.  survey::postStratify(design = data_prior, strata = ~posest,
# population = popc.types) posso salvar para depois
saveRDS(dadosPNADCOVID19, file = "dados_2020_05_final.rds")

Para construir o dataset selecionando algumas variáveis dentro do dataset completo, posso especificar vars=variaveis_selecionadas. A classe do objeto é idêntica ao procedimento com todas as variáveis.

variaveis_selecionadas <- c("UF", "A002", "A003", "A004", "A005", "C001", "C002",
    "C003", "C008", "C01012")
dadosPNADCOVID19.var <- get_covid(year = 2020, month = 5, vars = variaveis_selecionadas)
class(dadosPNADCOVID19.var)
saveRDS(dadosPNADCOVID19.var, file = "dados_2020_05_var.rds")

7 ALTERNATIVA 2 - srvyr

Uma alternativa 2 é separar as variáveis do objeto completo dadosPNADCOVID19 após design, ou seja, já deflacionado e com pós-estratificação. A funcao get_covid faz tudo isso mas precisa de web estável. Rodar com o dataset completo deixará mais pesado! Para separação das variáveis específicas a partir do dadosPNADCOVID19 completo, preciso das variáveis:

names(dadosPNADCOVID19.var[["variables"]])
# [1] 'Ano' 'V1013' 'UF' 'V1008' [5] 'V1012' 'Estrato' 'UPA' 'V1030' [9]
# 'V1031' 'V1032' 'posest' 'A001' [13] 'A002' 'A003' 'A004' 'A005' [17] 'C001'
# 'C002' 'C003' 'C008' [21] 'C01012' 'ID_DOMICILIO' 'Habitual' 'Efetivo' [25]
# 'CO3'

Usarei o pacote srvyr para operar com pipe (ver https://rpubs.com/amrofi/microdados_pnadc_srvyr). Fazendo o as_survey, o objeto que era survey.design agora ficará também tbl_svy. Da mesma forma que anteriormente, salvarei o dataset completo como pnad_srvyr.rds. Ao transformar em tbl_svy, especifico o conjunto de variáveis desejados no objeto de lista variaveis e uso a função select para separar as variáveis de interesse. Salvei o dataset pnad_srvyr_select.rds.

library(srvyr)
#
pnad_srvyr <- srvyr::as_survey(dadosPNADCOVID19)
saveRDS(pnad_srvyr, file = "pnad_srvyr.rds")
class(pnad_srvyr)
# 'tbl_svy' 'survey.design2' 'survey.design' posso selecionar as colunas
# simplesmente fazendo
variaveis <- c("Ano", "V1013", "UF", "V1008", "V1012", "Estrato", "UPA", "V1030",
    "V1031", "V1032", "posest", "A001", "A002", "A003", "A004", "A005", "C001", "C002",
    "C003", "C008", "C01012", "ID_DOMICILIO", "Habitual", "Efetivo", "CO3")
# usar select() in dplyr
pnad_srvyr_select <- pnad_srvyr %>%
    select(variaveis)
class(pnad_srvyr_select)
saveRDS(pnad_srvyr_select, file = "pnad_srvyr_select.rds")

Verifique que estou com mais variáveis em variaveis do que em variaveis_selecionadas, pois a função get_covid contêm variáveis que se referem ao design e não foram incluídas em variaveis_selecionadas.

names(pnad_srvyr_select[["variables"]])
# [1] 'Ano' 'V1013' 'UF' 'V1008' 'V1012' [6] 'Estrato' 'UPA' 'V1030' 'V1031'
# 'V1032' [11] 'posest' 'A001' 'A002' 'A003' 'A004' [16] 'A005' 'C001' 'C002'
# 'C003' 'C008' [21] 'C01012' 'ID_DOMICILIO' 'Habitual' 'Efetivo' 'CO3'
# comparar com o names de pnad_srvyr
names(pnad_srvyr[["variables"]])  # retornam 117 variaveis

7.1 Comparação entre objetos

dadosPNADCOVID19.var <- readRDS("D:/Documentos/disciplinas/econometria microdados/laboratorio microdados R/pnad_covid_gabrielassuncao/dados_2020_05_var.rds")
# Manipulando as variáveis -----
pnad_srvyr_select <- readRDS("D:/Documentos/disciplinas/econometria microdados/laboratorio microdados R/pnad_covid_gabrielassuncao/pnad_srvyr_select.rds")
library(survey)
# compararei as duas abordagens: survey.design x 'tbl_svy'
totalrenda1 <- svytotal(x = ~C01012, design = dadosPNADCOVID19.var, na.rm = TRUE)
totalrenda2 <- svytotal(x = ~C01012, design = pnad_srvyr_select, na.rm = TRUE)
totalrenda1

            total         SE
C01012 1.9216e+11 2298127060

totalrenda2

            total         SE
C01012 1.9216e+11 2298127060

cv(object = totalrenda1)

           C01012
C01012 0.01195936

cv(object = totalrenda2)

           C01012
C01012 0.01195936

confint(object = totalrenda1)

              2.5 %       97.5 %
C01012 187657087758 196665580296

confint(object = totalrenda2)

              2.5 %       97.5 %
C01012 187657087758 196665580296

confint(object = totalrenda1, level = 0.99)

              0.5 %       99.5 %
C01012 186241751004 198080917050

confint(object = totalrenda2, level = 0.99)

              0.5 %       99.5 %
C01012 186241751004 198080917050

# resultados identicos, ok

8 CORINGA

O uso deste coringa auxiliará para a alteração de dados para outro mês, por exemplo. O uso das opções do chunck colocando cache=TRUE em chunks mais lentos, auxilia a otimizar o tempo de procedimento.

# usarei um coringa x para depois compor uma rotina iterativa entre arquivos
x <- dadosPNADCOVID19.var

9 VARIÁVEIS CATEGÓRICAS

Estimação dos totais populacionais de categorias, utilizando variáveis categóricas, como o sexo:

totalsexo <- svytotal(x = ~A003, design = x, na.rm = TRUE)
totalsexo

               total SE
A003Homem  103090474  0
A003Mulher 107778927  0

9.1 Total de mais de uma variável categórica no mesmo código, separando com o operador +:

totalsexoraca <- svytotal(x = ~A003 + A004, design = x, na.rm = TRUE)
totalsexoraca

                 total     SE
A003Homem    103090474      0
A003Mulher   107778927      0
A004Branca    92966520 457236
A004Preta     18485175 227872
A004Amarela    1667234  78960
A004Parda     96989351 423539
A004Indígena    682710  42740
A004Ignorado     78411  17457

9.2 Total resultante do cruzamento de duas ou mais variáveis categóricas, com a função interaction:

totalsexoEraca <- svytotal(x = ~interaction(A003, A004), design = x, na.rm = TRUE)
ftable(x = totalsexoEraca)

                                                  A            B
interaction(A003, A004)Homem.Branca    44608026.970   250154.323
interaction(A003, A004)Mulher.Branca   48358492.853   250845.009
interaction(A003, A004)Homem.Preta      9344035.193   128048.173
interaction(A003, A004)Mulher.Preta     9141139.996   131489.011
interaction(A003, A004)Homem.Amarela     812101.042    45751.346
interaction(A003, A004)Mulher.Amarela    855133.370    42373.693
interaction(A003, A004)Homem.Parda     47948839.576   235865.787
interaction(A003, A004)Mulher.Parda    49040510.903   232020.589
interaction(A003, A004)Homem.Indígena    336751.697    22757.907
interaction(A003, A004)Mulher.Indígena   345958.158    24812.989
interaction(A003, A004)Homem.Ignorado     40719.523     9571.067
interaction(A003, A004)Mulher.Ignorado    37691.720     9879.737

10 ESTIMANDO MÉDIAS

mediarenda <- svymean(x = ~C01012, design = x, na.rm = TRUE)
mediarenda

         mean     SE
C01012 2325.2 25.882

10.1 Cálculo dos coeficientes de variação e intervalos de confiança

cv(object = mediarenda)

           C01012
C01012 0.01113082

confint(object = mediarenda)

          2.5 %   97.5 %
C01012 2274.522 2375.977

11 ESTIMAÇÃO DE PROPORÇÕES

propsexo <- svymean(x = ~A003, design = x, na.rm = TRUE)
propsexo

              mean SE
A003Homem  0.48888  0
A003Mulher 0.51112  0

11.1 Proporção de mais de uma variável:

propsexoraca <- svymean(x = ~A003 + A004, design = x, na.rm = TRUE)
propsexoraca

                   mean     SE
A003Homem    0.48888304 0.0000
A003Mulher   0.51111696 0.0000
A004Branca   0.44087250 0.0022
A004Preta    0.08766172 0.0011
A004Amarela  0.00790648 0.0004
A004Parda    0.45994986 0.0020
A004Indígena 0.00323760 0.0002
A004Ignorado 0.00037185 0.0001

11.2 Cruzamento de duas ou mais variáveis com a função interaction:

propsexoEraca <- svymean(x = ~interaction(A003, A004), design = x, na.rm = TRUE)
ftable(x = propsexoEraca)

                                                  A            B
interaction(A003, A004)Homem.Branca    0.2115433854 0.0011862998
interaction(A003, A004)Mulher.Branca   0.2293291138 0.0011895752
interaction(A003, A004)Homem.Preta     0.0443119540 0.0006072392
interaction(A003, A004)Mulher.Preta    0.0433497698 0.0006235566
interaction(A003, A004)Homem.Amarela   0.0038512038 0.0002169653
interaction(A003, A004)Mulher.Amarela  0.0040552748 0.0002009476
interaction(A003, A004)Homem.Parda     0.2273864266 0.0011185397
interaction(A003, A004)Mulher.Parda    0.2325634287 0.0011003047
interaction(A003, A004)Homem.Indígena  0.0015969681 0.0001079242
interaction(A003, A004)Mulher.Indígena 0.0016406276 0.0001176699
interaction(A003, A004)Homem.Ignorado  0.0001931030 0.0000453886
interaction(A003, A004)Mulher.Ignorado 0.0001787444 0.0000468524

12 ESTIMAÇÃO DE RAZÕES (RATIOS) OU TAXAS

txafastquarentena <- svyratio(numerator = ~(C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"),
    denominator = ~(C002 == "Sim"), design = x, na.rm = TRUE)
txafastquarentena

Ratio estimator: svyratio.survey.design2(numerator = ~(C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"), 
    denominator = ~(C002 == "Sim"), design = x, na.rm = TRUE)
Ratios=
                                                                                      C002 == "Sim"
C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"     0.8196935
SEs=
                                                                                      C002 == "Sim"
C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"   0.003123121

12.1 Intervalos de confiança e coeficiente de variação:

cv(object = txafastquarentena)

                                                                                      C002 == "Sim"
C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"   0.003810108

confint(object = txafastquarentena)

                                                                                                        2.5 %
C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"/C002 == "Sim" 0.8135723
                                                                                                       97.5 %
C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"/C002 == "Sim" 0.8258147

13 ESTIMAÇÃO DE MEDIANAS E QUANTIS

13.1 Mediana

medianarenda <- svyquantile(x = ~C01012, design = x, quantiles = 0.5, na.rm = TRUE)
medianarenda

$C01012
    quantile ci.2.5 ci.97.5       se
0.5     1500   1500    1597 24.74323

attr(,"hasci")
[1] TRUE
attr(,"class")
[1] "newsvyquantile"

13.2 Intervalo de confiança

medianarenda <- svyquantile(x = ~C01012, design = x, quantiles = 0.5, ci = TRUE,
    na.rm = TRUE)
medianarenda

$C01012
    quantile ci.2.5 ci.97.5       se
0.5     1500   1500    1597 24.74323

attr(,"hasci")
[1] TRUE
attr(,"class")
[1] "newsvyquantile"

13.3 Desvio-padrão e coeficiente de variação

SE(object = medianarenda)

  C01012 
24.74323 

cv(object = medianarenda)

    C01012 
0.01649548 

13.4 Múltiplos quantis

quantisrenda <- svyquantile(x = ~C01012, design = x, quantiles = c(0.1, 0.25, 0.5,
    0.75, 0.9), na.rm = TRUE)
quantisrenda

$C01012
     quantile ci.2.5 ci.97.5         se
0.1       600    600     650  12.754240
0.25     1045   1045    1050   1.275424
0.5      1500   1500    1597  24.743225
0.75     2500   2500    2700  51.016959
0.9      4500   4500    5000 127.542399

attr(,"hasci")
[1] TRUE
attr(,"class")
[1] "newsvyquantile"

14 ESTIMAÇÃO PARA VARIÁVEIS DE RENDIMENTO DEFLACIONADAS

14.1 Criar variável rendimento deflacionada

x$variables <- transform(x$variables, C01012_real = C01012 * Habitual)

14.2 Estimativa do total

totalrenda_real <- svytotal(x = ~C01012_real, design = x, na.rm = TRUE)
totalrenda_real

                 total         SE
C01012_real 1.9849e+11 2372672593

14.3 Média

mediarenda_real <- svymean(x = ~C01012_real, design = x, na.rm = TRUE)
mediarenda_real

              mean     SE
C01012_real 2401.8 26.721

15 ESTIMAÇÃO PARA UM DOMÍNIO

mediarendaM <- svymean(x = ~C01012, design = subset(x, A003 == "Mulher"), na.rm = TRUE)
mediarendaM

       mean     SE
C01012 2074 29.328

15.1 Condicionais com desigualdade

txafastquarentena25 <- svyratio(numerator = ~(C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"),
    denominator = ~(C002 == "Sim"), design = subset(x, A002 >= 25), na.rm = TRUE)
txafastquarentena25

Ratio estimator: svyratio.survey.design2(numerator = ~(C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"), 
    denominator = ~(C002 == "Sim"), design = subset(x, A002 >= 
        25), na.rm = TRUE)
Ratios=
                                                                                      C002 == "Sim"
C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"     0.8131341
SEs=
                                                                                      C002 == "Sim"
C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"    0.00330684

15.2 Utilizando múltiplas condições com os operadores lógicos & (“e”) e | (“ou”)

Esse caso serve para obter frequências relativas de cada nível (aqui utilizamos para o nível de instrução) para sexo = “Homem” e raça = “Parda”, com Idade (A002) “>30” anos.

nivelinstrHP30 <- svymean(x = ~A005, design = subset(x, A003 == "Homem" & A004 ==
    "Parda" & A002 > 30), na.rm = TRUE)
nivelinstrHP30

                                             mean     SE
A005Sem instrução                        0.070465 0.0014
A005Fundamental incompleto               0.344785 0.0032
A005Fundamental completa                 0.096410 0.0019
A005Médio incompleto                     0.072562 0.0017
A005Médio completo                       0.280868 0.0031
A005Superior incompleto                  0.032775 0.0012
A005Superior completo                    0.082917 0.0020
A005Pós-graduação, mestrado ou doutorado 0.019218 0.0009

15.3 Múltiplas análises em um mesmo domínio

Neste caso, faz-se uma subamostra usando a função subset.

dadosPNADCOVID19_mulheres <- subset(x, A003 == "Mulher")
dadosPNADCOVID19_mulheres

Stratified 1 - level Cluster Sampling design (with replacement)
With (14906) clusters.
subset(x, A003 == "Mulher")

16 ESTIMAÇÕES PARA VÁRIOS DOMÍNIOS

16.1 Estimação de quantidades em vários domínios mutuamente exclusivos

Homem x Mulher em cada nível de instrução:

freqSexoInstr <- svyby(formula = ~A003, by = ~A005, design = x, FUN = svymean, na.rm = TRUE)
freqSexoInstr

16.2 Frequência relativa de cada nível de instrução para cada sexo

freqInstrSexo <- svyby(formula = ~A005, by = ~A003, design = x, FUN = svymean, na.rm = TRUE)
freqInstrSexo

16.3 Renda média normalmente recebida em dinheiro (C01012) por Unidade da Federação (UF)

mediaRendaUF <- svyby(formula = ~C01012, by = ~UF, design = x, FUN = svymean, na.rm = TRUE)
mediaRendaUF

# farei tambem pelo objeto srvyr
library(srvyr)
mediaRendaUF2 <- pnad_srvyr_select %>%
    group_by(UF) %>%
    summarise(media = survey_mean(C01012, na.rm = T, vartype = "ci"))
mediaRendaUF2

16.4 Intervalo de confiança de (C01012) x (UF)

confint(object = mediaRendaUF)

                       2.5 %   97.5 %
Rondônia            1764.308 1974.814
Acre                1719.254 2017.290
Amazonas            1589.647 1926.892
Roraima             1793.859 2296.123
Pará                1622.586 1905.745
Amapá               1498.115 1876.096
Tocantins           1751.628 2043.755
Maranhão            1350.497 1501.071
Piauí               1514.803 1843.733
Ceará               1537.942 1755.738
Rio Grande do Norte 1654.700 2120.327
Paraíba             1623.194 1976.442
Pernambuco          1643.727 1979.794
Alagoas             1424.344 1685.165
Sergipe             1547.314 1930.720
Bahia               1427.324 1707.994
Minas Gerais        1917.897 2090.484
Espírito Santo      1987.506 2271.582
Rio de Janeiro      2632.150 2901.583
São Paulo           2735.098 3095.183
Paraná              2397.803 2659.572
Santa Catarina      2393.512 2551.635
Rio Grande do Sul   2379.729 2609.127
Mato Grosso do Sul  2124.166 2402.514
Mato Grosso         2122.413 2464.428
Goiás               1982.341 2223.232
Distrito Federal    3570.541 4451.600

library(ggplot2)
ggplot(as.data.frame(mediaRendaUF2), aes(x = UF, y = media)) + geom_point() + geom_errorbar(aes(ymin = media_low,
    ymax = media_upp), width = 0.2) + theme(axis.text.x = element_text(angle = 90)) +
    labs(title = "Rendimento Médio Domiciliar Estadual", caption = "Fonte: IBGE-PNADC, elaboração própria.")

16.5 Cruzamentos de variáveis categóricas com a função interaction

txafastquarentenaSexoRaca <- svyby(formula = ~(C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"),
    denominator = ~(C002 == "Sim"), by = ~interaction(A003, A004), design = x, FUN = svyratio,
    vartype = "cv", na.rm = TRUE)
txafastquarentenaSexoRaca

17 GRÁFICOS PARA DADOS AMOSTRAIS via survey

17.1 Histograma

svyhist(formula = ~as.numeric(C008), design = x, main = "Histograma com densidade",
    xlab = "Número de Horas Normalmente Trabalhadas")

svyhist(formula = ~as.numeric(C008), design = x, freq = TRUE, main = "Histograma com frequências absolutas",
    xlab = "Número de Horas Normalmente Trabalhadas")

17.2 Boxplot

17.2.1 Sem usar grupos

svyboxplot(formula = C008 ~ 1, design = x, main = "Boxplot do Número de Horas Normalmente Trabalhadas")

17.2.2 Com grupos definidos por variável

svyboxplot(formula = C008 ~ A003, design = x, main = "Boxplot do Número de Horas Normalmente Trabalhadas por Sexo")

17.3 Boxplot com todos os outliers

A opção anterior plotava apenas os outliers extremos. Agora têm-se todos os outliers.

svyboxplot(formula = C008 ~ A003, design = x, all.outliers = TRUE, main = "Boxplot do Número de Horas Normalmente Trabalhadas por Sexo")

17.4 Gráficos de dispersão

17.4.1 Estilo bolha

svyplot(formula = C01012 ~ C008, design = x, style = "bubble", xlab = "Número de Horas Normalmente Trabalhadas",
    ylab = "Rendimento Normalmente Recebido em Dinheiro")

17.4.2 Estilo transparente

svyplot(formula = C01012 ~ C008, design = x, style = "transparent", xlab = "Número de Horas Normalmente Trabalhadas",
    ylab = "Rendimento Normalmente Recebido em Dinheiro")

18 MODELAGEM COM PACOTE survey

18.1 Testes de Hipóteses

18.1.1 Teste-t para médias de rendimentos entre sexos

svyttest(formula = C01012 ~ A003, design = x)

    Design-based t-test

data:  C01012 ~ A003
t = -19.732, df = 14078, p-value < 2.2e-16
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -481.7773 -394.7078
sample estimates:
difference in mean 
         -438.2425 

18.1.2 Teste-t para médias de horas trabalhadas (C008) entre quem trabalhou ou não (C001)

svyttest(formula = as.numeric(C008) ~ C001, design = x)

    Design-based t-test

data:  as.numeric(C008) ~ C001
t = -21.681, df = 14099, p-value < 2.2e-16
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -2.520826 -2.102821
sample estimates:
difference in mean 
         -2.311824 

18.2 Modelos lineares

18.2.1 Regressão linear

modeloLin <- svyglm(formula = C01012 ~ A005 + A004 + A002, design = x)
summary(object = modeloLin)

Call:
svyglm(formula = C01012 ~ A005 + A004 + A002, design = x)

Survey design:
survey::postStratify(design = data_prior, strata = ~posest, population = popc.types)

Coefficients:
                                          Estimate Std. Error t value Pr(>|t|)
(Intercept)                               -482.700     70.679  -6.830 8.87e-12
A005Fundamental incompleto                 484.605     37.630  12.878  < 2e-16
A005Fundamental completa                   823.335     43.204  19.057  < 2e-16
A005Médio incompleto                       997.392     44.743  22.291  < 2e-16
A005Médio completo                        1219.472     41.736  29.219  < 2e-16
A005Superior incompleto                   1738.514     59.486  29.225  < 2e-16
A005Superior completo                     3240.245     65.525  49.451  < 2e-16
A005Pós-graduação, mestrado ou doutorado  5845.701    228.253  25.611  < 2e-16
A004Preta                                 -549.097     38.149 -14.394  < 2e-16
A004Amarela                                127.499    181.763   0.701  0.48303
A004Parda                                 -510.387     29.768 -17.146  < 2e-16
A004Indígena                              -630.133    110.201  -5.718 1.10e-08
A004Ignorado                             -1235.340    438.179  -2.819  0.00482
A002                                        34.720      1.039  33.404  < 2e-16
                                            
(Intercept)                              ***
A005Fundamental incompleto               ***
A005Fundamental completa                 ***
A005Médio incompleto                     ***
A005Médio completo                       ***
A005Superior incompleto                  ***
A005Superior completo                    ***
A005Pós-graduação, mestrado ou doutorado ***
A004Preta                                ***
A004Amarela                                 
A004Parda                                ***
A004Indígena                             ***
A004Ignorado                             ** 
A002                                     ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 3468448)

Number of Fisher Scoring iterations: 2

confint(object = modeloLin)

                                               2.5 %     97.5 %
(Intercept)                               -621.23912 -344.16034
A005Fundamental incompleto                 410.84619  558.36401
A005Fundamental completa                   738.64937  908.02020
A005Médio incompleto                       909.68925 1085.09473
A005Médio completo                        1137.66482 1301.27944
A005Superior incompleto                   1621.91339 1855.11524
A005Superior completo                     3111.80835 3368.68191
A005Pós-graduação, mestrado ou doutorado  5398.29368 6293.10753
A004Preta                                 -623.87349 -474.32110
A004Amarela                               -228.78123  483.77846
A004Parda                                 -568.73569 -452.03814
A004Indígena                              -846.14206 -414.12314
A004Ignorado                             -2094.22933 -376.45146
A002                                        32.68255   36.75731

18.2.2 Regressão Logística

modeloLog <- svyglm(formula = C002 ~ A003 + A004 + A002, design = x, family = "binomial")

summary(object = modeloLog)

Call:
svyglm(formula = C002 ~ A003 + A004 + A002, design = x, family = "binomial")

Survey design:
survey::postStratify(design = data_prior, strata = ~posest, population = popc.types)

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.006666   0.022279  45.185  < 2e-16 ***
A003Mulher    0.297449   0.015164  19.615  < 2e-16 ***
A004Preta    -0.112571   0.031502  -3.573 0.000354 ***
A004Amarela  -0.188918   0.098943  -1.909 0.056235 .  
A004Parda     0.010171   0.018533   0.549 0.583155    
A004Indígena -0.281620   0.136807  -2.059 0.039558 *  
A004Ignorado  2.161287   0.590344   3.661 0.000252 ***
A002          0.006435   0.000298  21.595  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 0.4859236)

Number of Fisher Scoring iterations: 5

confint(object = modeloLog)

                    2.5 %       97.5 %
(Intercept)   0.962996939  1.050335254
A003Mulher    0.267725574  0.327172653
A004Preta    -0.174319282 -0.050822352
A004Amarela  -0.382858513  0.005022564
A004Parda    -0.026156686  0.046498635
A004Indígena -0.549778851 -0.013461248
A004Ignorado  1.004137360  3.318437522
A002          0.005850794  0.007018923

19 ANÁLISE DE CONCENTRAÇÃO DE RENDA COM O PACOTE convey

O pacote convey requer uso da função convey_prep para adaptar o objeto do plano amostral do survey ao objeto que o pacote convey requer para as estimações (JACOB et al, 2021):

# para fazer o gini, precisa aplicar o convey_prep para usar o convey
library(convey)
dadosPNADCOVID19 <- convey_prep(design = x)
x <- dadosPNADCOVID19  # vou reatribuir
class(x)  # 'convey.design'  'survey.design2' 'survey.design' 

[1] "convey.design"  "survey.design2" "survey.design" 

19.1 Índice de Gini

giniHab <- svygini(formula = ~C01012, design = x, na.rm = TRUE)
giniHab

          gini     SE
C01012 0.47399 0.0043

cv(object = giniHab)

            C01012
C01012 0.009041141

19.2 Índice de Gini da renda normalmente recebida em dinheiro por Unidade da Federação:

giniUF <- svyby(formula = ~C01012, by = ~UF, design = x, FUN = svygini, na.rm = TRUE)
giniUF

confint(object = giniUF)

                        2.5 %    97.5 %
Rondônia            0.3852342 0.4296240
Acre                0.4084939 0.4695756
Amazonas            0.4470538 0.5061581
Roraima             0.4161803 0.5022714
Pará                0.4538041 0.5164472
Amapá               0.4282381 0.4931463
Tocantins           0.4189697 0.4732731
Maranhão            0.4232838 0.4634399
Piauí               0.4523957 0.5177020
Ceará               0.4438874 0.4898357
Rio Grande do Norte 0.4619424 0.5356331
Paraíba             0.4468704 0.5113244
Pernambuco          0.4529214 0.5223176
Alagoas             0.4106529 0.4745413
Sergipe             0.4606127 0.5300124
Bahia               0.4489249 0.5181641
Minas Gerais        0.4155318 0.4501937
Espírito Santo      0.4314812 0.4766895
Rio de Janeiro      0.4609274 0.4933300
São Paulo           0.4530767 0.4977770
Paraná              0.4201416 0.4555347
Santa Catarina      0.3700423 0.3957532
Rio Grande do Sul   0.4253530 0.4609469
Mato Grosso do Sul  0.3940728 0.4404393
Mato Grosso         0.3776868 0.4358869
Goiás               0.3962929 0.4409802
Distrito Federal    0.5073187 0.5414538

19.3 Curva de Lorenz

curvaLorenz <- svylorenz(formula = ~C01012, design = x, quantiles = seq(0, 1, 0.05),
    na.rm = TRUE)

20 REFLEXÕES

O arquivo tbl_svy fica mais pesado que o survey.design puro obtido pela get_covid. Uma vantagem do tbl_svy é que o formato tibble facilita gráficos em ggplot2 e manipulações do ambiente tidyverse. A desvantagem é que o arquivo é mais pesado que o formato survey.design puro. Vou continuar os procedimentos com o objeto dadosPNADCOVID19.var.

REFERÊNCIAS

ASSUNÇÃO, Gabriel. Análise de microdados da PNAD COVID19 Com os Pacotes COVIDIBGE e survey. RStudio/Rpubs, 2021. Disponível em: https://rpubs.com/gabriel-assuncao-ibge/covid.

ASSUNÇÃO, Gabriel; HIDALGO, Luna. COVIDIBGE: Downloading, Reading and Analysing PNAD COVID19 Microdata. 2021. R package version 0.1.7. Disponível em: https://CRAN.R-project.org/package=COVIDIBGE.

BARONE, Leonardo Sangali. PNAD Contínua no R com srvyr. RStudio/Rpubs, 2020. Disponível em: https://rpubs.com/leobarone/pnadc_srvyr.

BARONE, Leonardo Sangali. Cebrap.lab/CETIC - Introdução à Programação em R. Github, 2021. Disponível em: https://github.com/leobarone/cebrap_lab_cetic_programacao_r.

BRAGA, Douglas; ASSUNÇÃO, Gabriel. PNADcIBGE: Downloading, Reading and Analysing PNADC Microdata. R package version 0.6.5. 2021. Disponível em: https://CRAN.R-project.org/package=PNADcIBGE.

ELLIS, Greg Freedman; SCHNEIDER, Ben. srvyr: ‘dplyr’-Like Syntax for Summary Statistics of Survey Data. R package version 1.1.0. 2021. Disponível em: https://CRAN.R-project.org/package=srvyr.

FIGUEIREDO, Adriano Marcos Rodrigues. Econometria com Microdados da PNAD contínua de 2018 com R::PNADcIBGE, de BRAGA (2017). Campo Grande-MS,Brasil: RStudio/Rpubs, 2019. Disponível em http://rpubs.com/amrofi/pnadcibge_survey.

IBGE, PNAD Contínua - Pesquisa Nacional por Amostra de Domicílios Contínua: Microdados. Rio de Janeiro: IBGE, 2021. https://www.ibge.gov.br/estatisticas/downloads-estatisticas.html?caminho=Trabalho_e_Rendimento/Pesquisa_Nacional_por_Amostra_de_Domicilios_continua/Trimestral/Microdados/2021.

JACOB, Guilherme; PESSOA, Djalma; DAMICO, Anthony. Poverty and Inequality with Complex Survey Data. 2021.

LUMLEY, T. Survey: analysis of complex survey samples. R package version 3.35-1. 2019.

LUMLEY, T. Analysis of complex survey samples. Journal of Statistical Software, 9(1): 1-19. 2004.

PESSOA, Djalma; DAMICO, Anthony; JACOB, Guilherme. convey: Income Concentration Analysis with Complex Survey Samples. R package version 0.2.3. 2021.

NOTA:

Coloquei o tempo de execução abaixo, mas usei cache = TRUE em vários chunks de modo que o tempo total parece pequeno mas em vários chunks a rotina demora para executar, conforme seu micro.

end.time <- Sys.time()
time.taken <- end.time - start.time
time.taken

Time difference of 7.077375 secs

O primeiro procedimento levou (Time difference of) 7.640449 minutos, ou seja, considerando que usou os arquivos .rds previamente salvos.