Econometria com Microdados da PNAD contínua de 2018 com R::PNADcIBGE, de BRAGA (2017)

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. 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.

1 Introdução

Aplicação baseada em Braga (2017), para análise com microdados da Pesquisa Nacionalde Amostragem Domiciliar Contínua (PNAD Contínua), realizada pelo IBGE. Esta aplicação utiliza os pacotes PNADcIBGE (BRAGA, 2018) e survey (LUMLEY, 2019). Portanto, é interessante que o usuário tenha os pacotes previamente instalados.

O pacote PNADcIBGE permite a importação online dos dados, o que geralmente demorará dependendo da quantidade de variáveis e períodos solicitados (que podem ser para os microdados trimestrais e anuais). A função get_pnadc permite o download, leitura, rotulação e criação do objeto do plano amostral da pesquisa.

library(PNADcIBGE)
help("get_pnadc")

2 Microdados trimestrais

A importação dos microdados através da função get_pnadc é feita indicando-se o ano e o trimestre dos dados desejados nos argumentos da função. Por exemplo, para o 3º trimestre de 2017:

• para todas as variáveis

dadosPNADc = get_pnadc(year = 2017, quarter = 3)

ou

• para algumas variáveis apenas (de nomes ‘VD4001’ e ‘VD4002’, por exemplo)

dadosPNADc <- get_pnadc(year = 2017, quarter = 3, 
              vars=c("VD4001","VD4002"))
dadosPNADc
class(dadosPNADc)

3 Microdados do Ano de 2018

Utilizarei nesse exercício os dados anuais da 5ª visita. Para baixar os dados anuais, é preciso colocar o número da entrevista (interview), entre 1 e 5 visitas, conforme expresso em https://www.ibge.gov.br/estatisticas/sociais/habitacao/17270-pnad-continua.html?=&t=o-que-e. Questionam-se temas especificos em trimestres distintos (entrevistas de 1 a 5), de modo que nem todas as variáveis estão em todas as entrevistas, conforme o tema.

# https://rpubs.com/BragaDouglas/335574 vou tentar pegar os dados anuais se não
# estiver instalado, fazer a partir do CRAN install.packages('PNADcIBGE')


######## Temas e tópicos suplementares pesquisados em trimestres específicos do
######## ano: Educação (2o trimestre); e Acesso à televisão e à Internet e
######## posse de telefone móvel celular para uso pessoal (4o trimestre).
######## Temas e tópicos pesquisados ao longo do ano em determinada visita:
######## Habitação (1a visita); Características gerais dos moradores (1a
######## visita); Informações adicionais da força de trabalho (1a visita);
######## Outras formas de trabalho (afazeres domésticos, cuidados de pessoas,
######## produção para o próprio consumo e trabalho voluntário) (5a visita);
######## Trabalho de crianças e adolescentes (5a visita); e Rendimentos de
######## outras fontes (1a e 5a visitas).  UF Unidade da Federação V2007\t
######## \tSexo V2009\t \tIdade do morador na data de referência V2010\t \tCor
######## ou raça V3007\t já concluiu algum outro curso de graduação?
######## VD3004\t\tNível de instrução mais elevado alcançado (pessoas de 5 anos
######## ou mais de idade) padronizado para o Ensino fundamental sistema de 9
######## anos VD4001\t\tCondição em relação à força de trabalho na semana de
######## referência para pessoas de 14 anos ou mais de idade VD4002\t\tCondição
######## de ocupação na semana de referência para pessoas de 14 anos ou mais de
######## idade VD4020\t\tRendimento 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)
######## VD4035\t\tHoras efetivamente trabalhadas na semana de referência em
######## todos os trabalhos para pessoas de 14 anos ou mais de idade

4 Download dos dados selecionados

library(PNADcIBGE)
# variaveis selecionadas
variaveis_selecionadas <- c("UF", "V2007", "V2009", "V2010", "V3007", "VD3004", "VD4001",
    "VD4002", "VD4020", "VD4035")
dadosPNADc <- get_pnadc(year = 2018, interview = 5, vars = variaveis_selecionadas)
# nos dados trimestrais, posso usar labels, nos anuais não!  se indicar labels
# = FALSE, os dados virão, por exemplo, com 1 e 2 para homem e mulher se
# indicar labels = TRUE, virão Homem e Mulher invés de números
class(dadosPNADc)

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

4.1 Primeiros resultados

library(convey)
library(survey)
library(PNADcIBGE)
totalrenda <- svytotal(~VD4020, dadosPNADc, na.rm = T)
options(scipen = 999)  # desativando notacao cientifica
totalrenda

              total         SE
VD4020 200041565112 3011018309

cv(totalrenda)

           VD4020
VD4020 0.01505196

confint(totalrenda)  #intervalo de confiança de 95% (padrão)

              2.5 %       97.5 %
VD4020 194140077669 205943052554

confint(totalrenda, level = 0.99)  #intervalo de confiança de 99%

              0.5 %       99.5 %
VD4020 192285695918 207797434305

4.2 Contagens com svytotal

# contagem por sexo 1\tHomem e 2\tMulher
totalsexo <- svytotal(~V2007, dadosPNADc, na.rm = T)
totalsexo

                total     SE
V2007Homem  100260773 162783
V2007Mulher 107592520 162783

# contagem por sexo e raca sem interacao 1\tBranca 2\tPreta 3\tAmarela 4\tParda
# 5\tIndígena 9\tIgnorado

totalsexoraca <- svytotal(~V2007 + V2010, dadosPNADc, na.rm = T)
totalsexoraca

                    total       SE
V2007Homem    100260772.8 162782.6
V2007Mulher   107592520.2 162782.6
V2010Branca    89242075.1 468319.2
V2010Preta     17973280.5 233563.1
V2010Amarela    1387923.0  72921.0
V2010Parda     98602844.5 435722.1
V2010Indígena    639573.5  38665.4
V2010Ignorado      7596.5   3860.1

# contagem por sexo, raca e com interacao
totalsexoEraca <- svytotal(~interaction(V2007, V2010), dadosPNADc, na.rm = T)
ftable(totalsexoEraca)  # coluna A com valor e B com SE

                                                     A             B
interaction(V2007, V2010)Homem.Branca    42259053.9401   250303.2851
interaction(V2007, V2010)Mulher.Branca   46983021.1773   276661.1094
interaction(V2007, V2010)Homem.Preta      8919156.5371   122427.2642
interaction(V2007, V2010)Mulher.Preta     9054123.9518   135056.7166
interaction(V2007, V2010)Homem.Amarela     646442.6129    36998.8905
interaction(V2007, V2010)Mulher.Amarela    741480.3431    44929.1217
interaction(V2007, V2010)Homem.Parda     48130172.4300   240676.6372
interaction(V2007, V2010)Mulher.Parda    50472672.0982   253767.2288
interaction(V2007, V2010)Homem.Indígena    300173.0925    20839.5180
interaction(V2007, V2010)Mulher.Indígena   339400.3645    22466.5828
interaction(V2007, V2010)Homem.Ignorado      5774.1799     3436.2441
interaction(V2007, V2010)Mulher.Ignorado     1822.2727      899.6151

4.3 Médias com svymean

# média da renda
mediarenda <- svymean(~VD4020, dadosPNADc, na.rm = T)
mediarenda

         mean    SE
VD4020 2253.7 32.13

cv(mediarenda)

           VD4020
VD4020 0.01425679

confint(mediarenda)

          2.5 %   97.5 %
VD4020 2190.678 2316.625

# proporção em variável categórica sexo
propsexo <- svymean(~V2007, dadosPNADc, na.rm = T)
propsexo

               mean     SE
V2007Homem  0.48236 0.0008
V2007Mulher 0.51764 0.0008

# proporção em variável categórica sexo e raca sem interacao
propsexoraca <- svymean(~V2007 + V2010, dadosPNADc, na.rm = T)
propsexoraca

                     mean     SE
V2007Homem    0.482363168 0.0008
V2007Mulher   0.517636832 0.0008
V2010Branca   0.429351269 0.0023
V2010Preta    0.086470992 0.0011
V2010Amarela  0.006677416 0.0004
V2010Parda    0.474386733 0.0021
V2010Indígena 0.003077043 0.0002
V2010Ignorado 0.000036547 0.0000

# proporção em variável categórica sexo e raca com cruzamento (interaction)
propsexoEraca <- svymean(~interaction(V2007, V2010), dadosPNADc, na.rm = T)
ftable(propsexoEraca)

                                                      A              B
interaction(V2007, V2010)Homem.Branca    0.203311928958 0.001204230549
interaction(V2007, V2010)Mulher.Branca   0.226039340052 0.001331040300
interaction(V2007, V2010)Homem.Preta     0.042910826229 0.000589008057
interaction(V2007, V2010)Mulher.Preta    0.043560166024 0.000649769434
interaction(V2007, V2010)Homem.Amarela   0.003110090793 0.000178004832
interaction(V2007, V2010)Mulher.Amarela  0.003567325455 0.000216157854
interaction(V2007, V2010)Homem.Parda     0.231558383008 0.001157915921
interaction(V2007, V2010)Mulher.Parda    0.242828349600 0.001220895879
interaction(V2007, V2010)Homem.Indígena  0.001444158465 0.000100260706
interaction(V2007, V2010)Mulher.Indígena 0.001632884231 0.000108088655
interaction(V2007, V2010)Homem.Ignorado  0.000027780074 0.000016532065
interaction(V2007, V2010)Mulher.Ignorado 0.000008767110 0.000004328125

4.4 Teste entre médias de renda para grupos de sexo

# ttest entre médias de VD4020 (renda) para os grupos de V2007 (sexo)
tt <- svyttest(VD4020 ~ V2007, dadosPNADc)
tt

    Design-based t-test

data:  VD4020 ~ V2007
t = -22.479, df = 11486, p-value < 0.00000000000000022
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -616.3543 -517.4853
sample estimates:
difference in mean 
         -566.9198 

confint(tt, level = 0.9)

[1] -608.4055 -525.4341
attr(,"conf.level")
[1] 0.9

# média com critério/domínio (# 1\tHomem e 2\tMulher)
mediarendaM <- svymean(~VD4020, subset(dadosPNADc, V2007 == "2"), na.rm = T)
mediarendaM

       mean  SE
VD4020  NaN NaN

mediarendaH <- svymean(~VD4020, subset(dadosPNADc, V2007 == "1"), na.rm = T)
mediarendaH

       mean  SE
VD4020  NaN NaN

# Outro exemplo para testar diferenças entre médias de horas trabalhadas, entre
# concluintes e não concluintes de graduação:
svyttest(as.numeric(VD4035) ~ V3007, dadosPNADc)

    Design-based t-test

data:  as.numeric(VD4035) ~ V3007
t = -1.5093, df = 4896, p-value = 0.1313
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -2.0600330  0.2678703
sample estimates:
difference in mean 
        -0.8960814 

4.5 Taxa de desocupação com svyratio

# taxa de desocupação VD4002 1\tPessoas ocupadas 2\tPessoas desocupadas ' ' Não
# aplicável VD4001 1\tPessoas na força de trabalho 2\tPessoas fora da força de
# trabalho ' ' Não aplicável

txdesocup <- svyratio(~VD4002 == "2", ~VD4001 == "1", dadosPNADc, na.rm = T)
txdesocup

Ratio estimator: svyratio.survey.design2(~VD4002 == "2", ~VD4001 == "1", dadosPNADc, 
    na.rm = T)
Ratios=
              VD4001 == "1"
VD4002 == "2"           NaN
SEs=
              VD4001 == "1"
VD4002 == "2"           NaN

cv(txdesocup)

              VD4001 == "1"
VD4002 == "2"           NaN

confint(txdesocup)

                            2.5 % 97.5 %
VD4002 == "2"/VD4001 == "1"   NaN    NaN

# taxa de desocupação com desigualdade (V2009 idade >= 25 anos)
txdesocup25 <- svyratio(~VD4002 == "2", ~VD4001 == "1", subset(dadosPNADc, V2009 >=
    25), na.rm = T)
txdesocup25

Ratio estimator: svyratio.survey.design2(~VD4002 == "2", ~VD4001 == "1", subset(dadosPNADc, 
    V2009 >= 25), na.rm = T)
Ratios=
              VD4001 == "1"
VD4002 == "2"           NaN
SEs=
              VD4001 == "1"
VD4002 == "2"           NaN

# subset mulher (V2007 = 2)
dadosPNADc_mulheres <- subset(dadosPNADc, V2007 == "2")
dadosPNADc_mulheres

Stratified 1 - level Cluster Sampling design (with replacement)
With (12068) clusters.
subset(dadosPNADc, V2007 == "2")

# 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, dadosPNADc, svymean, na.rm = T)
freqSexoInstr

4.6 Renda de cada estado (Unidade da Federação, UF)

# média da renda por estado
mediaRendaUF <- svyby(~VD4020, ~UF, dadosPNADc, svymean, na.rm = T)
mediaRendaUF

confint(mediaRendaUF)

                       2.5 %   97.5 %
Rondônia            1701.948 2033.838
Acre                1529.899 1857.053
Amazonas            1507.315 1901.412
Roraima             1852.872 2611.356
Pará                1327.691 1560.463
Amapá               1796.217 2485.678
Tocantins           1588.482 1990.057
Maranhão            1178.403 1372.702
Piauí               1218.331 1546.512
Ceará               1330.282 1667.728
Rio Grande do Norte 1358.361 1735.170
Paraíba             1491.559 1937.565
Pernambuco          1498.189 1927.964
Alagoas             1366.425 1617.456
Sergipe             1333.608 1758.711
Bahia               1290.937 1868.100
Minas Gerais        1856.328 2040.821
Espírito Santo      1907.016 2183.382
Rio de Janeiro      2295.842 2540.547
São Paulo           2735.184 3183.695
Paraná              2279.642 2532.827
Santa Catarina      2429.163 2608.902
Rio Grande do Sul   2368.897 2611.567
Mato Grosso do Sul  2077.927 2415.445
Mato Grosso         2173.148 2511.155
Goiás               1967.921 2200.717
Distrito Federal    3638.691 4729.577

# taxa de desocupação com sexo (1 e 2) e raça (1 a 5 e 9) e interacao
txdesocupSexoRaca <- svyby(~VD4002 == "2", ~interaction(V2007, V2010), dadosPNADc,
    svyratio, denominator = ~VD4001 == "1", na.rm = T, vartype = "cv")
txdesocupSexoRaca

4.7 Regressão para determinar a renda

# o sistema escolar mudou da VD3001 (8 anos) para a VD3004 (9 anos) modelo
# linear
modeloLin <- svyglm(VD4020 ~ VD3004 + V2010 + V2009, dadosPNADc)
(summary(modeloLin))

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

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

Coefficients:
                                            Estimate Std. Error t value
(Intercept)                                 -711.169     70.630 -10.069
VD3004Fundamental incompleto ou equivalente  514.885     28.516  18.056
VD3004Fundamental completo ou equivalente    924.614     32.956  28.056
VD3004Médio incompleto ou equivalente       1147.284     41.066  27.937
VD3004Médio completo ou equivalente         1372.317     36.185  37.925
VD3004Superior incompleto ou equivalente    1860.788     55.356  33.615
VD3004Superior completo                     4491.540    113.957  39.414
V2010Preta                                  -631.888     36.667 -17.233
V2010Amarela                                 504.112    274.165   1.839
V2010Parda                                  -585.071     29.536 -19.808
V2010Indígena                               -413.749    172.428  -2.400
V2010Ignorado                               3291.189    940.931   3.498
V2009                                         38.708      1.365  28.347
                                                        Pr(>|t|)    
(Intercept)                                 < 0.0000000000000002 ***
VD3004Fundamental incompleto ou equivalente < 0.0000000000000002 ***
VD3004Fundamental completo ou equivalente   < 0.0000000000000002 ***
VD3004Médio incompleto ou equivalente       < 0.0000000000000002 ***
VD3004Médio completo ou equivalente         < 0.0000000000000002 ***
VD3004Superior incompleto ou equivalente    < 0.0000000000000002 ***
VD3004Superior completo                     < 0.0000000000000002 ***
V2010Preta                                  < 0.0000000000000002 ***
V2010Amarela                                            0.065983 .  
V2010Parda                                  < 0.0000000000000002 ***
V2010Indígena                                           0.016432 *  
V2010Ignorado                                           0.000471 ***
V2009                                       < 0.0000000000000002 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 4431486)

Number of Fisher Scoring iterations: 2

confint(modeloLin)

                                                 2.5 %     97.5 %
(Intercept)                                 -849.61605 -572.72136
VD3004Fundamental incompleto ou equivalente  458.98788  570.78138
VD3004Fundamental completo ou equivalente    860.01480  989.21401
VD3004Médio incompleto ou equivalente       1066.78668 1227.78054
VD3004Médio completo ou equivalente         1301.38858 1443.24560
VD3004Superior incompleto ou equivalente    1752.27991 1969.29567
VD3004Superior completo                     4268.16562 4714.91442
V2010Preta                                  -703.76036 -560.01503
V2010Amarela                                 -33.29889 1041.52233
V2010Parda                                  -642.96771 -527.17485
V2010Indígena                               -751.73753  -75.75953
V2010Ignorado                               1446.80369 5135.57410
V2009                                         36.03148   41.38467

library(stargazer)
stargazer(modeloLin, title = "Título: Resultado da Regressão com Survey Data", align = TRUE,
    type = "text", style = "all")

Título: Resultado da Regressão com Survey Data
===============================================================================
                                                    Dependent variable:        
                                            -----------------------------------
                                                          VD4020               
-------------------------------------------------------------------------------
VD3004Fundamental incompleto ou equivalente             514.885***             
                                                         (28.516)              
                                                        t = 18.056             
                                                         p = 0.000             
VD3004Fundamental completo ou equivalente               924.614***             
                                                         (32.956)              
                                                        t = 28.056             
                                                         p = 0.000             
VD3004Médio incompleto ou equivalente                  1,147.284***            
                                                         (41.066)              
                                                        t = 27.937             
                                                         p = 0.000             
VD3004Médio completo ou equivalente                    1,372.317***            
                                                         (36.185)              
                                                        t = 37.925             
                                                         p = 0.000             
VD3004Superior incompleto ou equivalente               1,860.788***            
                                                         (55.356)              
                                                        t = 33.615             
                                                         p = 0.000             
VD3004Superior completo                                4,491.540***            
                                                         (113.957)             
                                                        t = 39.414             
                                                         p = 0.000             
V2010Preta                                              -631.888***            
                                                         (36.667)              
                                                        t = -17.233            
                                                         p = 0.000             
V2010Amarela                                             504.112*              
                                                         (274.165)             
                                                         t = 1.839             
                                                         p = 0.066             
V2010Parda                                              -585.071***            
                                                         (29.536)              
                                                        t = -19.808            
                                                         p = 0.000             
V2010Indígena                                           -413.749**             
                                                         (172.428)             
                                                        t = -2.400             
                                                         p = 0.017             
V2010Ignorado                                          3,291.189***            
                                                         (940.931)             
                                                         t = 3.498             
                                                        p = 0.0005             
V2009                                                    38.708***             
                                                          (1.365)              
                                                        t = 28.347             
                                                         p = 0.000             
Constant                                                -711.169***            
                                                         (70.630)              
                                                        t = -10.069            
                                                         p = 0.000             
-------------------------------------------------------------------------------
Observations                                              175,910              
Log Likelihood                                        -1,698,047.000           
Akaike Inf. Crit.                                      3,396,120.000           
Residual Deviance                           1,943,598,666,188.000 (df = 11475) 
Null Deviance                               2,418,101,798,275.000 (df = 175909)
===============================================================================
Note:                                               *p<0.1; **p<0.05; ***p<0.01

4.8 Gini

# para os calculos de gini preciso preparar pelo convey
library(convey)
dados_pnadcc <- convey_prep(dadosPNADc)
giniHab <- svygini(~VD4020, dados_pnadcc, na.rm = TRUE)
giniHab

          gini     SE
VD4020 0.52202 0.0048

cv(giniHab)

            VD4020
VD4020 0.009150046

giniUF <- svyby(~VD4020, by = ~UF, dados_pnadcc, svygini, na.rm = TRUE)
giniUF

confint(giniUF)

                        2.5 %    97.5 %
Rondônia            0.4195730 0.4951908
Acre                0.5020427 0.5569024
Amazonas            0.5119566 0.5762671
Roraima             0.5066617 0.5715221
Pará                0.5006856 0.5551459
Amapá               0.4586638 0.5284547
Tocantins           0.4584070 0.5302670
Maranhão            0.4960746 0.5475850
Piauí               0.5574596 0.6185949
Ceará               0.5265586 0.5970560
Rio Grande do Norte 0.5041850 0.5747581
Paraíba             0.5388949 0.6171288
Pernambuco          0.4971035 0.5742745
Alagoas             0.4259214 0.4883830
Sergipe             0.5126857 0.5955621
Bahia               0.5170227 0.6288732
Minas Gerais        0.4792008 0.5102320
Espírito Santo      0.4562336 0.5049896
Rio de Janeiro      0.4574140 0.4893069
São Paulo           0.4968729 0.5434477
Paraná              0.4548393 0.4916522
Santa Catarina      0.3997182 0.4246232
Rio Grande do Sul   0.4721840 0.5036270
Mato Grosso do Sul  0.4552676 0.4999415
Mato Grosso         0.4302313 0.4852954
Goiás               0.4306950 0.4644831
Distrito Federal    0.5445286 0.5956048

curvaLorenz <- svylorenz(~VD4020, dados_pnadcc, quantiles = seq(0, 1, 0.05), na.rm = TRUE)

Referências

BRAGA, Douglas. Análise de microdados da PNAD Contínua com os Pacotes PNADcIBGE e survey. Brasil: RStudio/Rpubs e IBGE, 2017. Disponível em https://rpubs.com/BragaDouglas/335574.05 dezembro 2017.

BRAGA, Douglas. PNADcIBGE: Downloading, Reading and Analysing PNADc Microdata. R package version 0.4.3. 2018. Disponível em https://CRAN.R-project.org/package=PNADcIBGE.

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.

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

Time difference of 6.21822 mins