Microdados com R: exemplo da PNADC com srvyr

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 PNADC com srvyr. Campo Grande-MS,Brasil: RStudio/Rpubs, 2021. Disponível em http://rpubs.com/amrofi/microdados_pnadc_srvyr.

1 Introdução

Este é um script baseado em Barone (2020) e https://github.com/leobarone/cebrap_lab_cetic_programacao_r. Ele utiliza o pacote PNADcIBGE (BRAGA e ASSUNÇÃO, 2021) e o pacote srvyr de Ellis e Schneider (2021). Recomendo que olhe os vídeos no Youtube sobre microdados (video mais teorico em https://youtu.be/k3rLK8V_38s; sobre a NHANES https://youtu.be/UFIFnWnjYeM; e https://youtu.be/hhNA8MR__aA para a POF 2017-2018; ) e especialmente sobre a PNADC: https://youtu.be/qC22jz0FA0U.

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 PNADcIBGE facilitará o trabalho e apresenta algumas vantagens como o de adequar os objetos de dados ao ambiente tidyverse. O tutorial 07 de Barone (2021) auxiliará neste exercício.

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

2 Carregando os dados com design

Chamarei os dados da PNADC para o segundo trimestre de 2021 (IBGE, 2021). Farei uma contagem do tempo para o leitor ter uma ideia de quanto tempo demora para rodar o chunk. O código está preparado para cache = TRUE de modo a reduzir o tempo quando os chunks não sofrerem novas alterações. O leitor atento verá que ao executar a rotina, a função puxa as informações dos microdados, o dicionário e deflatores todos automaticamente, colocando o objeto já em formato de survey design. É importante entender que as rotinas com microdados e demandantes de conexão estável com internet podem gerar erros e “travamentos” indesejados, associados ainda ao ambiente de software e hardware do momento. A opção savedir pode ser curinga para salvar os dados baixados a fim de ganhar tempo no futuro.

# pacotes iniciais
library(PNADcIBGE)
start.time1 <- Sys.time()
# variaveis selecionadas
variaveis_selecionadas <- c("UF", "V2007", "V2009", "V2010", "V3007", "VD3004", "VD4001",
    "VD4002", "VD4020", "VD4035")
pnad_srv <- get_pnadc(year = 2021, quarter = 2, vars = variaveis_selecionadas, savedir = "C:\\Users\\amrof\\dados_pnadc")
end.time1 <- Sys.time()
time.taken1 <- end.time1 - start.time1
time.taken1

Time difference of 30.68214 secs

# Time difference of 34.59961 secs (primeira tentativa)

É possível verificar que o objeto é da classe de survey design:

class(pnad_srv)

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

Conforme instrução de Barone (2020), os dados carregados com a função get_pnadc já foram, por default, transformados com a dependência do pacote survey de Lumley (2019), com o plano amostral incorporado via funções svydesign e postStratify (ver FIGUEIREDO, 2019) e, desta forma, não vêm como data frame nem table data frame.

Como Barone (2020) indica, >[…] o objetivo destas funções é transformar data frames em objetos que contenham peso, estrato, identifiquem a unidade de ponderação e quaisquer outros aspectos do desenho da amostra necessários para a produção de estatística.

3 Carregando os dados sem design

Se o pesquisador quiser dados tipo data frame utilizando a função get_pnadc, deve alterar a rotina fazendo uso da opção design = FALSE. Recomendamos bastante cuidado pois o pesquisador deverá então utilizar o peso dentro da especificação que será feita no srvyr.

# comentario: por alguma razao está travando quando Knit mas executa no console
# executei no console e salvei o objeto pnad_df e carrego no próximo chunk
start.time2 <- Sys.time()
#
library(PNADcIBGE)
pnad_df <- get_pnadc(year = 2021, quarter = 2, design = FALSE, vars = variaveis_selecionadas)
end.time2 <- Sys.time()
time.taken2 <- end.time2 - start.time2
time.taken2

Time difference of 27.25782 secs

# Time difference of 23.58122 secs (primeira tentativa)

# load('pnad_df.Rdata')

É possível ver que agora o objeto pnad_df é um data frame.

class(pnad_df)

[1] "tbl_df"     "tbl"        "data.frame"

Nesta altura, é interessante ver como as variáveis são carregadas. A função glimpse do tidyverse permitirá visualizar que as variáveis qualitativas já são importadas como fatores (tipo <fct>).

library(tidyverse)
glimpse(pnad_df)

Rows: 356,239
Columns: 21
$ Ano       <chr> "2021", "2021", "2021", "2021", "2021", "2021", "2021", "202~
$ Trimestre <chr> "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", ~
$ UF        <fct> Rondônia, Rondônia, Rondônia, Rondônia, Rondônia, Rondônia, ~
$ UPA       <chr> "110000016", "110000016", "110000016", "110000016", "1100000~
$ Estrato   <chr> "1110011", "1110011", "1110011", "1110011", "1110011", "1110~
$ V1008     <chr> "02", "02", "03", "03", "03", "03", "04", "04", "04", "04", ~
$ V1027     <dbl> 77.44494, 77.44494, 77.44494, 77.44494, 77.44494, 77.44494, ~
$ V1029     <dbl> 546949, 546949, 546949, 546949, 546949, 546949, 546949, 5469~
$ posest    <chr> "111", "111", "111", "111", "111", "111", "111", "111", "111~
$ V2003     <chr> "01", "02", "01", "02", "03", "04", "01", "02", "03", "04", ~
$ V2007     <fct> Mulher, Homem, Mulher, Homem, Homem, Homem, Homem, Mulher, H~
$ V2009     <dbl> 37, 31, 45, 20, 8, 66, 51, 45, 26, 23, 55, 18, 14, 71, 24, 4~
$ V2010     <fct> Parda, Parda, Branca, Branca, Branca, Parda, Parda, Parda, P~
$ V3007     <fct> NA, NA, NA, NA, NA, NA, NA, NA, Não, NA, NA, NA, NA, NA, NA,~
$ VD3004    <fct> Sem instrução e menos de 1 ano de estudo, Fundamental incomp~
$ VD4001    <fct> Pessoas fora da força de trabalho, Pessoas na força de traba~
$ VD4002    <fct> NA, Pessoas desocupadas, Pessoas desocupadas, Pessoas desocu~
$ VD4020    <dbl> NA, NA, NA, NA, NA, NA, 1600, NA, 1045, NA, 1500, NA, NA, NA~
$ VD4035    <dbl> NA, NA, NA, NA, NA, NA, 40, 36, 48, 36, 48, NA, NA, NA, 48, ~
$ Habitual  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ~
$ Efetivo   <dbl> 1.005071, 1.005071, 1.005071, 1.005071, 1.005071, 1.005071, ~

4 Aplicando o design ao data frame

O pacote PNADcIBGE aplicar o plano amostral ao data frame com a função pnadc_design. Neste caso, ele já considerará os nomes das variáveis para o estrato, UPA etc de modo a aplicar o design. Recomendo olhar também: https://cran.r-project.org/web/packages/srvyr/vignettes/srvyr-database.html

library(srvyr)
library(PNADcIBGE)
pnad_srv_2 <- pnadc_design(pnad_df)
class(pnad_srv_2)

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

Como não usaremos mais o objeto pand_srv_2, removerei da memória.

rm(pnad_srv_2)

5 O pacote srvyr na PNAD Contínua

É possível utilizar funções do dplyr na PNADC já carregada do objeto de survey design. Por exemplo, filtraremos para a variável sexo, V2007, para apenas “Mulher”:

library(srvyr)
library(tidyverse)
pnad_srv_fem <- pnad_df %>%
    filter(V2007 == "Mulher")
glimpse(pnad_srv_fem)

Rows: 184,659
Columns: 21
$ Ano       <chr> "2021", "2021", "2021", "2021", "2021", "2021", "2021", "202~
$ Trimestre <chr> "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", "2", ~
$ UF        <fct> Rondônia, Rondônia, Rondônia, Rondônia, Rondônia, Rondônia, ~
$ UPA       <chr> "110000016", "110000016", "110000016", "110000016", "1100000~
$ Estrato   <chr> "1110011", "1110011", "1110011", "1110011", "1110011", "1110~
$ V1008     <chr> "02", "03", "04", "04", "05", "06", "07", "07", "07", "08", ~
$ V1027     <dbl> 77.44494, 77.44494, 77.44494, 77.44494, 77.44494, 77.44494, ~
$ V1029     <dbl> 546949, 546949, 546949, 546949, 546949, 546949, 546949, 5469~
$ posest    <chr> "111", "111", "111", "111", "111", "111", "111", "111", "111~
$ V2003     <chr> "01", "01", "02", "04", "03", "01", "01", "03", "04", "02", ~
$ V2007     <fct> Mulher, Mulher, Mulher, Mulher, Mulher, Mulher, Mulher, Mulh~
$ V2009     <dbl> 37, 45, 45, 23, 14, 71, 40, 8, 22, 56, 46, 61, 47, 12, 57, 6~
$ V2010     <fct> Parda, Branca, Parda, Parda, Parda, Parda, Parda, Parda, Par~
$ V3007     <fct> NA, NA, NA, NA, NA, NA, NA, NA, Não, NA, NA, NA, NA, NA, NA,~
$ VD3004    <fct> Sem instrução e menos de 1 ano de estudo, Médio completo ou ~
$ VD4001    <fct> Pessoas fora da força de trabalho, Pessoas na força de traba~
$ VD4002    <fct> NA, Pessoas desocupadas, Pessoas ocupadas, Pessoas ocupadas,~
$ VD4020    <dbl> NA, NA, NA, NA, NA, NA, 1045, NA, 1045, NA, NA, NA, 1100, NA~
$ VD4035    <dbl> NA, NA, 36, 36, NA, NA, 36, NA, 36, NA, NA, NA, 40, NA, 30, ~
$ Habitual  <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, ~
$ Efetivo   <dbl> 1.005071, 1.005071, 1.005071, 1.005071, 1.005071, 1.005071, ~

Como não usaremos mais o objeto pand_srv_fem, removerei da memória.

rm(pnad_srv_fem)

Transformaremos o objeto de survey com a função as_survey, do pacote srvyr:

class(pnad_srv)

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

pnad_srvyr <- srvyr::as_survey(pnad_srv)
class(pnad_srvyr)

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

Agora aparece a classe “tbl_svy” que nos auxiliará nos filtros. Podemos agora usar o filter com o objeto de design:

pnad_mulheres_srvyr <- pnad_srvyr %>%
    filter(V2007 == "Mulher")
glimpse(pnad_mulheres_srvyr)

List of 10
 $ cluster   :'data.frame': 356239 obs. of  1 variable:
  ..$ UPA: Factor w/ 14438 levels "1110011.110000016",..: 1 1 1 1 1 1 1 1 1 1 ...
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language ~UPA
  .. .. ..- attr(*, "variables")= language list(UPA)
  .. .. ..- attr(*, "factors")= int [1, 1] 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. ..- attr(*, "term.labels")= chr "UPA"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 0
  .. .. ..- attr(*, ".Environment")=<environment: 0x00000000161f1980> 
  .. .. ..- attr(*, "predvars")= language list(UPA)
  .. .. ..- attr(*, "dataClasses")= Named chr "character"
  .. .. .. ..- attr(*, "names")= chr "UPA"
 $ strata    :'data.frame': 356239 obs. of  1 variable:
  ..$ Estrato: chr [1:356239] "1110011" "1110011" "1110011" "1110011" ...
  ..- attr(*, "terms")=Classes 'terms', 'formula'  language ~Estrato
  .. .. ..- attr(*, "variables")= language list(Estrato)
  .. .. ..- attr(*, "factors")= int [1, 1] 1
  .. .. .. ..- attr(*, "dimnames")=List of 2
  .. .. ..- attr(*, "term.labels")= chr "Estrato"
  .. .. ..- attr(*, "order")= int 1
  .. .. ..- attr(*, "intercept")= int 1
  .. .. ..- attr(*, "response")= int 0
  .. .. ..- attr(*, ".Environment")=<environment: 0x00000000161f1980> 
  .. .. ..- attr(*, "predvars")= language list(Estrato)
  .. .. ..- attr(*, "dataClasses")= Named chr "character"
  .. .. .. ..- attr(*, "names")= chr "Estrato"
 $ has.strata: logi TRUE
 $ prob      : Named num [1:356239] 0.0088 Inf 0.0088 Inf Inf ...
  ..- attr(*, "names")= chr [1:356239] "1" "2" "3" "4" ...
 $ allprob   :'data.frame': 356239 obs. of  1 variable:
  ..$ V1027: num [1:356239] 0.0129 0.0129 0.0129 0.0129 0.0129 ...
 $ call      : 'quoteless_text' chr "Called via srvyr"
 $ variables : tibble [356,239 x 21] (S3: tbl_df/tbl/data.frame)
  ..$ Ano      : chr [1:356239] "2021" "2021" "2021" "2021" ...
  ..$ Trimestre: chr [1:356239] "2" "2" "2" "2" ...
  ..$ UF       : Factor w/ 27 levels "Rondônia","Acre",..: 1 1 1 1 1 1 1 1 1 1 ...
  ..$ UPA      : chr [1:356239] "110000016" "110000016" "110000016" "110000016" ...
  ..$ Estrato  : chr [1:356239] "1110011" "1110011" "1110011" "1110011" ...
  ..$ V1008    : chr [1:356239] "02" "02" "03" "03" ...
  ..$ V1027    : num [1:356239] 77.4 77.4 77.4 77.4 77.4 ...
  ..$ V1029    : num [1:356239] 546949 546949 546949 546949 546949 ...
  ..$ posest   : chr [1:356239] "111" "111" "111" "111" ...
  ..$ V2003    : chr [1:356239] "01" "02" "01" "02" ...
  ..$ V2007    : Factor w/ 2 levels "Homem","Mulher": 2 1 2 1 1 1 1 2 1 2 ...
  ..$ V2009    : num [1:356239] 37 31 45 20 8 66 51 45 26 23 ...
  ..$ V2010    : Factor w/ 6 levels "Branca","Preta",..: 4 4 1 1 1 4 4 4 4 4 ...
  ..$ V3007    : Factor w/ 2 levels "Sim","Não": NA NA NA NA NA NA NA NA 2 NA ...
  ..$ VD3004   : Factor w/ 7 levels "Sem instrução e menos de 1 ano de estudo",..: 1 2 5 2 2 2 5 5 6 5 ...
  ..$ VD4001   : Factor w/ 2 levels "Pessoas na força de trabalho",..: 2 1 1 1 NA 2 1 1 1 1 ...
  ..$ VD4002   : Factor w/ 2 levels "Pessoas ocupadas",..: NA 2 2 2 NA NA 1 1 1 1 ...
  ..$ VD4020   : num [1:356239] NA NA NA NA NA ...
  ..$ VD4035   : num [1:356239] NA NA NA NA NA NA 40 36 48 36 ...
  ..$ Habitual : num [1:356239] 1 1 1 1 1 1 1 1 1 1 ...
  ..$ Efetivo  : num [1:356239] 1.01 1.01 1.01 1.01 1.01 ...
 $ fpc       :List of 2
  ..$ popsize : NULL
  ..$ sampsize: int [1:356239, 1] 51 51 51 51 51 51 51 51 51 51 ...
  ..- attr(*, "class")= chr "survey_fpc"
 $ pps       : logi FALSE
 $ postStrata:List of 1
  ..$ : int [1:356239] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "oldweights")= Named num [1:356239] 77.4 77.4 77.4 77.4 77.4 ...
  .. .. ..- attr(*, "names")= chr [1:356239] "1" "2" "3" "4" ...
  .. ..- attr(*, "weights")= Named num [1:356239] 114 114 114 114 114 ...
  .. .. ..- attr(*, "names")= chr [1:356239] "1" "2" "3" "4" ...
 - attr(*, "class")= chr [1:3] "tbl_svy" "survey.design2" "survey.design"
 - attr(*, "survey_vars")= list()
  ..- attr(*, "class")= chr "survey_vars"

Como não usaremos mais o objeto pnad_mulheres_srvyr, removerei da memória. Deixarei apenas o objeto pnad_srvyr.

rm(pnad_mulheres_srvyr)
rm(pnad_srv)
rm(pnad_df)

Voltando para o objeto completo, pnad_srvyr, para nossas análises. Agora calcularei a média de VD4020 - Rendimento mensal efetivo de todos os trabalhos para pessoas de 14 anos ou mais de idade (apenas para pessoas que receberam em dinheiro, produtos ou mercadorias em qualquer trabalho) pelo survey e pelo srvyr. Na sequência farei a média de cada Nível de instrução (VD3004) mais elevado alcançado (pessoas de 5 anos ou mais de idade) padronizado para o Ensino fundamental sistema de 9 anos para o subset sexo (V2007 == “Homem”) e raça (V2010 == “Parda”) e Idade (V2009 > 30). O resultado é a proporção por grau de instrução para homens pardos com mais de 30 anos.

library(survey)
svymean(~VD4020, pnad_srvyr, na.rm = T, multicore = TRUE)

         mean     SE
VD4020 2496.5 27.948

pnad_srvyr %>%
    summarise(media_renda = survey_mean(VD4020, na.rm = T))

svymean(~VD3004, subset(pnad_srvyr, V2007 == "Homem" & V2010 == "Parda" & V2009 >
    30), na.rm = T, multicore = TRUE)

                                                   mean     SE
VD3004Sem instrução e menos de 1 ano de estudo 0.090754 0.0021
VD3004Fundamental incompleto ou equivalente    0.390448 0.0039
VD3004Fundamental completo ou equivalente      0.086470 0.0022
VD3004Médio incompleto ou equivalente          0.047356 0.0015
VD3004Médio completo ou equivalente            0.261310 0.0040
VD3004Superior incompleto ou equivalente       0.023322 0.0012
VD3004Superior completo                        0.100339 0.0026

A mesma conta é feita com o pacote srvyr e o operador pipe:

pnad_srvyr %>%
    filter(V2007 == "Homem", V2010 == "Parda", V2009 > 30) %>%
    group_by(VD3004) %>%
    summarise(freq = survey_mean())

Como alertado por Barone (2020), uma vantagem importante do pacote srvyr comparativamente ao procedimento pelo survey é a classe do objeto resultante do srvyr ser data frame. Isto facilitará o uso de outros pacotes amigáveis como o dplyr e ggplot2.

Um terceiro exemplo comparativo entre os pacotes srvyr e survey, utilizando uma tabela de duas entradas: de sexo (V2007) e nível de instrução (VD3004), adaptando de Barone (2020) e Figueiredo (2019). Primeiro, computamos pelo survey e na sequência pelo srvyr:

start.time3 <- Sys.time()
#
require(survey)
# Se desejamos estimar a frequência relativa de homens e mulheres em cada nível
# de instrução, usamos o seguinte código:

freqSexoInstr <- svyby(~V2007, ~VD3004, pnad_srvyr, svymean, na.rm = T, multicore = TRUE)
freqSexoInstr

end.time3 <- Sys.time()
time.taken3 <- end.time3 - start.time3
time.taken3

Time difference of 22.1365 secs

start.time4 <- Sys.time()
# require(srvyr)
pnad_srvyr %>%
    group_by(V2007, VD3004) %>%
    summarise(freq = survey_mean())

end.time4 <- Sys.time()
time.taken4 <- end.time4 - start.time4
time.taken4  # Time difference of 48.83992 secs

Time difference of 50.57616 secs

5.1 Renda média de cada estado

start.time5 <- Sys.time()
# require(srvyr)
mediaRendaUF <- pnad_srvyr %>%
    group_by(UF) %>%
    summarise(media = survey_mean(VD4020, na.rm = T, vartype = "ci"))
mediaRendaUF

end.time5 <- Sys.time()
time.taken5 <- end.time5 - start.time5
time.taken5  # Time difference of 1.232956 mins

Time difference of 1.556456 mins

library(ggplot2)
ggplot(as.data.frame(mediaRendaUF), 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.")

ggsave("grafico1.png", plot = last_plot())

5.2 Regressão

require(srvyr)
pnad_srvyr %>%
    {
        survey::svyglm(VD4020 ~ VD3004 + V2010 + V2009, .)
    } %>%
    summary()

Call:
svyglm(formula = VD4020 ~ VD3004 + V2010 + V2009, design = .)

Survey design:
Called via srvyr

Coefficients:
                                            Estimate Std. Error t value
(Intercept)                                 -670.854     97.505  -6.880
VD3004Fundamental incompleto ou equivalente  458.202     52.265   8.767
VD3004Fundamental completo ou equivalente    802.783     53.365  15.043
VD3004Médio incompleto ou equivalente       1074.966     55.261  19.452
VD3004Médio completo ou equivalente         1256.136     55.258  22.732
VD3004Superior incompleto ou equivalente    1934.721     88.172  21.943
VD3004Superior completo                     4329.309     93.053  46.525
V2010Preta                                  -731.067     40.554 -18.027
V2010Amarela                                 314.649    231.683   1.358
V2010Parda                                  -650.844     36.260 -17.949
V2010Indígena                               -689.579    150.170  -4.592
V2010Ignorado                               4137.078   2576.035   1.606
V2009                                         41.415      1.767  23.445
                                            Pr(>|t|)    
(Intercept)                                 6.24e-12 ***
VD3004Fundamental incompleto ou equivalente  < 2e-16 ***
VD3004Fundamental completo ou equivalente    < 2e-16 ***
VD3004Médio incompleto ou equivalente        < 2e-16 ***
VD3004Médio completo ou equivalente          < 2e-16 ***
VD3004Superior incompleto ou equivalente     < 2e-16 ***
VD3004Superior completo                      < 2e-16 ***
V2010Preta                                   < 2e-16 ***
V2010Amarela                                   0.174    
V2010Parda                                   < 2e-16 ***
V2010Indígena                               4.43e-06 ***
V2010Ignorado                                  0.108    
V2009                                        < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 5850787)

Number of Fisher Scoring iterations: 2

Referências

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.

6 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 16.53238 secs