O pacote COVIDIBGE

O pacote COVIDIBGE foi desenvolvido para facilitar o download, importação e análise dos dados amostrais da Pesquisa Nacional por Amostra de Domicílios - COVID19 (Estatísticas Experimentais) realizada pelo Instituto Brasileiro de Geografia e Estatística - IBGE.

A PNAD COVID19 possui somente um tipo de microdados:

Mensal, que contém a parte básica investigada pela pesquisa, contendo as informações sobre mercado de trabalho e síndrome gripal.

Maiores informações sobre a pesquisa e os temas investigados podem ser encontradas no site oficial do IBGE.

Através do objeto criado com este pacote, é possível utilizar o pacote survey para realizar análises considerando o efeito do esquema de seleção utilizado no plano amostral complexo da pesquisa e calcular corretamente as medidas de erro das estimativas, considerando o estimador de pós-estratificação utilizado na pesquisa.

Instalação do pacote

O Pacote está disponível no repositório CRAN do R, onde pode ser acessada sua documentação. Tais informações também podem ser acessadas pelo repositório no GitHub do mantenedor do pacote através do endereço: https://github.com/Gabriel-Assuncao/COVIDIBGE.
É importante sempre conferir se a versão instalada em seu computador é a mais recente disponível no CRAN.
A instalação do pacote pode ser feita pelo seguinte comando:

install.packages("COVIDIBGE")

Antes de utilizar o pacote é necessário carregá-lo no R através do comando:

library(COVIDIBGE)

Importação de microdados

Para a importação dos microdados da PNAD COVID19, há duas opções. A primeira é uma opção online, que exige que o computador esteja conectado à Internet. A segunda opção não necessita de conexão a Internet, e é utilizada para a leitura de microdados que estejam no disco rígido.

Importação online

A função get_covid permite o download, leitura, rotulação, inserção das variáveis para deflacionamento e criação do objeto do plano amostral da pesquisa.
Esta função pode ser usada para os microdados mensais disponíveis.

help("get_covid")

Microdados mensais

A importação dos microdados mensais através da função get_covid é muito simples. Basta indicar o ano e o mês dos microdados desejados nos argumentos da função. Abaixo um exemplo de como importar os microdados do mês de maio de 2020:

dadosPNADCOVID19 <- get_covid(year=2020, month=5)

Apenas com esse comando, os microdados da pesquisa já estão disponíveis no objeto dadosPNADCOVID19 para análise.

Além dos argumentos year e month, que indicam, respectivamente, o ano e o mês dos microdados a serem baixados, a função get_covid ainda possui outros cinco argumentos que podem ser modificados no download para este tipo de microdados:

vars: Este argumento recebe um vetor de caracteres com o nome das variáveis a serem baixadas. Caso nenhuma variável seja passada, todas as variáveis disponíveis na pesquisa são baixadas. É útil caso deseje trabalhar com poucas variáveis, pois assim o objeto ocupará um espaço menor na memória do computador;
labels: Um argumento lógico que indica se os níveis das variáveis categóricas devam ser rotuladas de acordo com o dicionário da pesquisa. O default é rotulá-los;
deflator: Um argumento lógico que indica se as variáveis para deflacionamento serão incluídas nos microdados. O default é incluí-las;
design: Um argumento lógico que indica se a função deve retornar um objeto do plano amostral para análise com o pacote survey. Caso design=FALSE, a função retorna apenas um data-frame com os microdados. É altamente recomendado que mantenha essa opção como TRUE, caso contrário suas análises poderão ser feitas de forma incorreta;
reload: Um argumento lógico que indica se a função deve baixar novamente os arquivos no diretório de referência mesmo que estes arquivos já tenham sido baixados;
curlopts: Uma lista com indicação das opções para utilização nas funções do pacote RCurl;
savedir: O endereço onde devem ser salvos os arquivos baixados. Padrão é utilizar uma pasta temporária.

Por exemplo, um usuário que deseje trabalhar apenas com as variáveis de Realização de trabalho ou algum bico na semana de referência (C001) e Situação de afastamento temporário de algum trabalho (C002) do mês de maio de 2020, pode utilizar o argumento vars para selecionar apenas estas variáveis:

dadosPNADCOVID19 <- get_covid(year=2020, month=5, vars=c("C001","C002"))

Através do código abaixo, podemos ver que o objeto retornado pela função é da classe survey.design ou svyrep.design, sendo o tipo de objeto utilizado para análises de dados amostrais complexos através do pacote survey.

dadosPNADCOVID19

## Stratified 1 - level Cluster Sampling design (with replacement)
## With (14906) clusters.
## survey::postStratify(design = data_prior, strata = ~posest, population = popc.types)

class(dadosPNADCOVID19)

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

Caso o usuário não queira trabalhar com este objeto, ele pode escolher a opção design=FALSE para baixar os microdados brutos:

dadosPNADCOVID19_brutos <- get_covid(year=2020, month=5, vars=c("C001","C002"), design=FALSE)

O Objeto resultante é um data-frame com as variáveis selecionadas e as variáveis relacionadas ao processo de amostragem:

dadosPNADCOVID19_brutos

## # A tibble: 349,306 × 18
##    Ano   V1013 UF      V1008 V1012 Estrato UPA    V1030 V1031 V1032 posest  A001
##    <chr> <chr> <fct>   <chr> <chr> <chr>   <chr>  <int> <dbl> <dbl> <chr>  <int>
##  1 2020  5     Rondôn… 1     1     1110011 1100… 124609  181.  179. 1125       1
##  2 2020  5     Rondôn… 1     1     1110011 1100… 154290  181.  203. 1123       2
##  3 2020  5     Rondôn… 1     1     1110011 1100… 139963  181.  203. 1122       3
##  4 2020  5     Rondôn… 1     1     1110011 1100… 140194  181.  235. 1111       4
##  5 2020  5     Rondôn… 10    4     1110011 1100… 122285  181.  222. 1115       1
##  6 2020  5     Rondôn… 10    4     1110011 1100… 150235  181.  214. 1124       2
##  7 2020  5     Rondôn… 10    4     1110011 1100… 154290  181.  203. 1123       3
##  8 2020  5     Rondôn… 10    4     1110011 1100… 146160  181.  251. 1112       4
##  9 2020  5     Rondôn… 11    3     1110011 1100…  54453  181.  162. 1127       1
## 10 2020  5     Rondôn… 13    1     1110011 1100… 122285  181.  222. 1115       1
## # ℹ 349,296 more rows
## # ℹ 6 more variables: C001 <fct>, C002 <fct>, ID_DOMICILIO <chr>,
## #   Habitual <dbl>, Efetivo <dbl>, CO3 <dbl>

Também é possível baixar os dados sem incluir os rótulos dos níveis, através do argumento labels.

dadosPNADCOVID19_brutos_sem <- get_covid(year=2020, month=5, vars=c("C001","C002"), labels=FALSE, design=FALSE)

Perceba que agora os níveis das categorias são representados por números:

dadosPNADCOVID19_brutos_sem

## # A tibble: 349,306 × 18
##    Ano   V1013 UF    V1008 V1012 Estrato UPA      V1030 V1031 V1032 posest  A001
##    <chr> <chr> <chr> <chr> <chr> <chr>   <chr>    <int> <dbl> <dbl> <chr>  <int>
##  1 2020  5     11    1     1     1110011 110000… 124609  181.  179. 1125       1
##  2 2020  5     11    1     1     1110011 110000… 154290  181.  203. 1123       2
##  3 2020  5     11    1     1     1110011 110000… 139963  181.  203. 1122       3
##  4 2020  5     11    1     1     1110011 110000… 140194  181.  235. 1111       4
##  5 2020  5     11    10    4     1110011 110000… 122285  181.  222. 1115       1
##  6 2020  5     11    10    4     1110011 110000… 150235  181.  214. 1124       2
##  7 2020  5     11    10    4     1110011 110000… 154290  181.  203. 1123       3
##  8 2020  5     11    10    4     1110011 110000… 146160  181.  251. 1112       4
##  9 2020  5     11    11    3     1110011 110000…  54453  181.  162. 1127       1
## 10 2020  5     11    13    1     1110011 110000… 122285  181.  222. 1115       1
## # ℹ 349,296 more rows
## # ℹ 6 more variables: C001 <chr>, C002 <chr>, ID_DOMICILIO <chr>,
## #   Habitual <dbl>, Efetivo <dbl>, CO3 <dbl>

Com relação ao processo de inserção das variáveis para deflacionamento através do argumento deflator, existe um documento disponível no site oficial do IBGE para auxílio no processo de deflacionamento das variáveis de rendimento.

Importação offline

Caso não seja possível utilizar a importação online através da função get_covid, é possível criar o mesmo objeto para análise com arquivos que estejam em disco. Para isto, são utilizadas quatro funções:

read_covid: Para a leitura do arquivo .csv dos microdados;
covid_labeller: Opcional. Coloca os rótulos dos níveis nas variáveis categóricas;
covid_deflator: Opcional. Insere as variáveis para deflacionamento nos microdados;
covid_design: Cria o objeto do plano amostral para a análise com o pacote survey.

Para utilizar estas funções, primeiramente é necessário ter os microdados e sua documentação no disco e extraídos dos arquivos compactados, se for o caso.
Estes arquivos podem ser baixados diretamente do FTP do IBGE.

Para a função read_covid é utilizado o arquivo de CSV contendo os microdados.

Exemplo de leitura para o mês de maio de 2020:

dadosPNADCOVID19 <- read_covid(microdata="PNAD_COVID_052020.csv")

Assim como na importação online, pode ser utilizado o argumento vars para definir quais variáveis serão lidas.

A função covid_labeller utiliza o objeto criado pela função read_covid e o arquivo de dicionário das variáveis presente na documentação da pesquisa.

dadosPNADCOVID19 <- covid_labeller(data_covid=dadosPNADCOVID19, dictionary.file="Dicionario_PNAD_COVID_052020.xls")

A utilização de covid_labeller é opcional e não interfere nos resultados obtidos.

A função covid_deflator utiliza o objeto criado pela função read_covid e o arquivo de deflatores presente na documentação da pesquisa.

dadosPNADCOVID19 <- covid_deflator(data_covid=dadosPNADCOVID19, deflator.file="deflator_PNAD_COVID.xls")

A utilização de covid_deflator é opcional e não interfere nos resultados obtidos, desde que não seja o objetivo a análise das variáveis de rendimento em valor real.

A criação do objeto do plano amostral é feita diretamente pela função covid_design aplicada no objeto que contém os microdados:

dadosPNADCOVID19 <- covid_design(data_covid=dadosPNADCOVID19)

O objeto resultante é o mesmo retornado pela função get_covid e pode ser utilizado para análises com o pacote survey.

Análise com pacote survey

Devido ao plano amostral da PNAD Contínua, é necessário que sejam utilizadas ferramentas específicas para a análise de dados amostrais complexos. O pacote survey é um pacote criado especificamente para análise e modelagem de dados provenientes de pesquisas com estes tipos de planos amostrais. Maiores detalhes sobre o pacote podem ser encontrados no site do autor.

Para os exemplos, utilizaremos microdados da PNAD COVID19 do mês de maio de 2020, com algumas variáveis selecionadas de acordo com a tabela abaixo:

Variável	Descrição
UF	Unidade da Federação
A001A	Condição no domicílio
A002	Idade do morador na data de referência
A003	Sexo
A004	Cor ou raça
A005	Escolaridade
C001	Realização de trabalho ou algum bico na semana de referência para pessoas de 14 anos ou mais de idade
C002	Situação de afastamento temporário de algum trabalho para pessoas de 14 anos ou mais de idade
C003	Motivo do afastamento temporário do trabalho para pessoas de 14 anos ou mais de idade
C008	Número de horas normalmente trabalhadas para pessoas de 14 anos ou mais de idade
C01012	Rendimento normalmente recebido em dinheiro para pessoas de 14 anos ou mais de idade
F001	Tipo do domicílio

Importando os dados:

variaveis_selecionadas <- c("UF","A001A","A002","A003","A004","A005","C001","C002","C003","C008","C01012","F001")
dadosPNADCOVID19 <- get_covid(year=2020, month=5, vars=variaveis_selecionadas)

O objeto dadosPNADCOVID19 vai ser utilizado com o pacote survey para a análise dos microdados.

Para iniciar a análise dos dados é necessário carregar o pacote survey.

library(survey)

Estimando totais

A função do pacote para a estimação de totais populacionais é a svytotal. Sua sintaxe precisa de três parâmetros principais:

O nome da variável que se deseja calcular o total, precedido por um ~;
O nome do objeto do plano amostral;
A opção na.rm=TRUE, que remove as observações onde a variável é não-aplicável.

Variáveis numéricas

Ela pode ser utilizada para a estimação do total de uma variável numérica, como a renda normalmente recebida em dinheiro:

totalrenda <- svytotal(x=~C01012, design=dadosPNADCOVID19, na.rm=TRUE)
totalrenda

##               total         SE
## C01012 192161334027 2298127060

Além da estimativa do total da renda normalmente recebida em dinheiro, o comando também retorna o erro padrão (SE) dessa estimativa.
A partir desse resultado, também é possível calcular coeficientes de variação e intervalos de confiança para a estimativa:

cv(object=totalrenda)

##            C01012
## C01012 0.01195936

confint(object=totalrenda)

##               2.5 %       97.5 %
## C01012 187657087758 196665580296

confint(object=totalrenda, level=0.99)

##               0.5 %       99.5 %
## C01012 186241751004 198080917050

Variáveis categóricas

Também é possível estimar totais populacionais de categorias, utilizando variáveis categóricas, como o sexo:

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

##                total SE
## A003Homem  103090474  0
## A003Mulher 107778927  0

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

totalsexoraca <- svytotal(x=~A003+A004, design=dadosPNADCOVID19, 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

Também é possível estimar o total resultante do cruzamento de duas ou mais variáveis categóricas, com a função interaction:

totalsexoEraca <- svytotal(x=~interaction(A003,A004), design=dadosPNADCOVID19, 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

A função ftable possibilita uma saída mais organizada da tabela no caso de cruzamento de variáveis. Para todos esses objetos também podem ser utilizadas as funções cv e confint.

Estimando médias

A média de uma variável numérica é estimada através da função svymean, que possui uma sintaxe idêntica à svytotal. O exemplo do total da renda normalmente recebida em dinheiro pode ser facilmente reescrito para médias:

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

##          mean     SE
## C01012 2325.2 25.882

E podemos calcular coeficientes de variação e intervalos de confiança da mesma forma:

cv(object=mediarenda)

##            C01012
## C01012 0.01113082

confint(object=mediarenda)

##           2.5 %   97.5 %
## C01012 2274.522 2375.977

Estimando proporções

Utilizando variáveis categóricas, é possível estimar a frequência relativa de cada categoria. Isso pode ser feito também através da função svymean, com uma sintaxe análoga à utilizada para estimar os totais das categorias.

Para estimar a proporção de cada sexo na população:

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

##               mean SE
## A003Homem  0.48888  0
## A003Mulher 0.51112  0

Da mesma forma, podemos estimar a proporção de mais de uma variável:

propsexoraca <- svymean(x=~A003+A004, design=dadosPNADCOVID19, 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

E estimando a proporção de um cruzamento de duas ou mais variáveis com a função interaction:

propsexoEraca <- svymean(x=~interaction(A003,A004), design=dadosPNADCOVID19, 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

Neste caso, as frequências relativas calculadas são em relação ao total, não a uma marginal.

Estimando razões

Além das proporções, também é possível estimar razões entre duas variáveis. Um exemplo de razão é a taxa de afastamento do trabalho por motivo de quarentena, isolamento, distanciamento social ou férias coletivas: ela é a razão entre o total de pessoas afastadas do trabalho por motivo de quarentena, isolamento, distanciamento social ou férias coletivas pelo total de pessoas temporariamente afastadas.

A função para estimar razões é a svyratio. Sua sintaxe utiliza quatro argumentos: a variável cujo total estará no numerador, a variável cujo total estará no denominador, o objeto do plano amostral e a opção de remover valores não aplicáveis.

Um exemplo para a estimativa dessa taxa é:

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

## Warning in `[.survey.design2`(design, !nas, ): 7 strata have only one PSU in
## this subset.

txafastquarentena

## Ratio estimator: svyratio.survey.design2(numerator = ~(C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"), 
##     denominator = ~(C002 == "Sim"), design = dadosPNADCOVID19, 
##     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

Cálculos de coeficiente de variação e intervalos de confiança para essa taxa:

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

Estimando medianas e quantis

Medianas e quantis de variáveis numéricas são estimados através da função svyquantile. Além dos argumentos utilizados para estimar a média, é necessário definir os quantis a serem calculados no argumento quantiles.

Para calcular a mediana da renda normalmente recebida em dinheiro, basta utilizar quantiles=0.5:

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

## $C01012
##       0.5
## [1,] 1500
## 
## attr(,"hasci")
## [1] FALSE
## attr(,"class")
## [1] "newsvyquantile"

Para obtenção do intervalo de confiança pela função svyquantile, é necessário utilizar o parâmetro ci=TRUE e a partir dele, utilizar as funções SE e cv para estimar o erro padrão e coeficiente de variação, respectivamente.

medianarenda <- svyquantile(x=~C01012, design=dadosPNADCOVID19, 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"

SE(object=medianarenda)

##   C01012 
## 24.74323

cv(object=medianarenda)

##     C01012 
## 0.01649548

Além disso, é possível calcular vários quantis simultaneamente, colocando um vetor no argumento quantiles.
No código abaixo, estimamos, além da mediana, o primeiro e nono decis e primeiro e terceiro quartis:

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

## $C01012
##      0.1 0.25  0.5 0.75  0.9
## [1,] 600 1045 1500 2500 4500
## 
## attr(,"hasci")
## [1] FALSE
## attr(,"class")
## [1] "newsvyquantile"

Estimação para variáveis de rendimento deflacionadas

Para estimar totais e médias de variáveis de rendimento deflacionadas, ou seja, ao invés de obter as estimativas sob o valor nominal destas variáveis, passar a obter as estimativas sob o valor real. Primeiramente, é necessário realizar a criação da variável de rendimento com os valores deflacionados nos microdados.

Para criação desta nova variável nos microdados utilizamos a função transform.

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

Após a criação da variável de rendimento deflacionada, para obter a estimativa do total desta variável basta utilizar a função svytotal conforme visto anteriormente:

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

##                    total         SE
## C01012_real 198491400251 2372672593

E para obtenção da média desta nova variável, basta utilizar a função svymean:

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

##               mean     SE
## C01012_real 2401.8 26.721

O processo de deflacionamento é essencial no processo de comparação de uma série histórica dos microdados quando estiver sendo analisada uma variável de rendimento.

Estimação para um domínio

Muitas vezes queremos fazer uma análise para um domínio específico da população. Com o pacote survey, esse domínio pode ser selecionado utilizando a função subset no objeto do plano amostral, aplicando a condição que define aquele domínio. É necessário que a seleção desse domínio nos dados seja feita somente após a criação do objeto que define o plano amostral, caso contrário poderá ser obtido resultados incorretos.

Para a seleção dos domínios podem ser utilizadas condicionais com igualdade, como por exemplo para estimar a média da renda normalmente recebida em dinheiro de mulheres:

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

## Warning in `[.survey.design2`(design, nas == 0, ): 4 strata have only one PSU in
## this subset.

mediarendaM

##        mean     SE
## C01012 2074 29.328

Ou condicionais com desigualdade, como por exemplo estimar a taxa de afastamento do trabalho por motivo de quarentena, isolamento, distanciamento social ou férias coletivas para pessoas com idade igual ou superior a 25 anos:

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

## Warning in `[.survey.design2`(design, !nas, ): 8 strata have only one PSU in
## this subset.

txafastquarentena25

## Ratio estimator: svyratio.survey.design2(numerator = ~(C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"), 
##     denominator = ~(C002 == "Sim"), design = subset(dadosPNADCOVID19, 
##         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

Além disso, é possível utilizar múltiplas condições com os operados lógicos & (“e”) e | (“ou”), para estimar a frequência relativa de cada nível de instrução para homens pardos com mais de 30 anos, por exemplo.

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

## Warning in `[.survey.design2`(x, r, ): 2 strata have only one PSU in this
## subset.

## Warning in `[.survey.design2`(design, nas == 0, ): 2 strata have only one PSU in
## this subset.

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

Além disso, caso deseje realizar diversas análises para um mesmo domínio, é possível criar um objeto que contenha apenas os dados daquele domínio a partir de um subset do objeto original do plano amostral.

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

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

Outro ponto importante de ser mencionado diante de análises por domínio, é sobre a obtenção de indicadores em nível de domicílio, que para analisar corretamente é necessário realizar o recorte na variável de condição no domicílio, conforme exemplo apresentado abaixo da proporção do tipo de domicílio:

proptipodom <- svymean(x=~F001, design=subset(dadosPNADCOVID19, A001A=="Pessoa responsável pelo domicílio"), na.rm=TRUE)
proptipodom

##                                  mean     SE
## F001Próprio - já pago       0.6245074 0.0023
## F001Próprio - ainda pagando 0.0784109 0.0018
## F001Alugado                 0.1779809 0.0017
## F001Cedido por empregador   0.0125954 0.0005
## F001Cedido por familiar     0.0929706 0.0013
## F001Cedido de outra forma   0.0095728 0.0004
## F001Outra condição          0.0039620 0.0003

Estimação para vários domínios

Em outros casos há interesse em estimar quantidades de interesse para diversos domínios mutuamente exclusivos, a fim de possibilitar comparações. Para isto, podemos utilizar a função svyby.

Para utilizá-la são necessários os seguintes argumentos:

A variável da qual se deseja calcular a quantidade;
A variável que define os domínios;
O objeto do plano amostral;
A função utilizada para calcular a quantidade de interesse (svytotal, svymean, svyratio, svyquantile, …).

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(formula=~A003, by=~A005, design=dadosPNADCOVID19, FUN=svymean, na.rm=TRUE)

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 1 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, nas == 0, ): 1 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 1 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, nas == 0, ): 1 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 5 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, nas == 0, ): 5 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 3 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, nas == 0, ): 3 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 19 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, nas == 0, ): 19 strata have only one PSU
## in this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 11 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, nas == 0, ): 11 strata have only one PSU
## in this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 73 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, nas == 0, ): 73 strata have only one PSU
## in this subset.

freqSexoInstr

##                                                                      A005
## Sem instrução                                               Sem instrução
## Fundamental incompleto                             Fundamental incompleto
## Fundamental completa                                 Fundamental completa
## Médio incompleto                                         Médio incompleto
## Médio completo                                             Médio completo
## Superior incompleto                                   Superior incompleto
## Superior completo                                       Superior completo
## Pós-graduação, mestrado ou doutorado Pós-graduação, mestrado ou doutorado
##                                      A003Homem A003Mulher se.A003Homem
## Sem instrução                        0.5035142  0.4964858  0.002554670
## Fundamental incompleto               0.5128943  0.4871057  0.001430836
## Fundamental completa                 0.4867680  0.5132320  0.003782277
## Médio incompleto                     0.5206277  0.4793723  0.003143817
## Médio completo                       0.4822497  0.5177503  0.001762822
## Superior incompleto                  0.4775020  0.5224980  0.004402266
## Superior completo                    0.4186103  0.5813897  0.002912055
## Pós-graduação, mestrado ou doutorado 0.4103251  0.5896749  0.005833542
##                                      se.A003Mulher
## Sem instrução                          0.002554670
## Fundamental incompleto                 0.001430836
## Fundamental completa                   0.003782277
## Médio incompleto                       0.003143817
## Médio completo                         0.001762822
## Superior incompleto                    0.004402266
## Superior completo                      0.002912055
## Pós-graduação, mestrado ou doutorado   0.005833542

Caso desejemos estimar o contrário, a frequência relativa de cada nível de instrução para cada sexo, basta inverter as variáveis correspondentes na sintaxe:

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

##          A003 A005Sem instrução A005Fundamental incompleto
## Homem   Homem         0.1097536                  0.3174499
## Mulher Mulher         0.1035139                  0.2883734
##        A005Fundamental completa A005Médio incompleto A005Médio completo
## Homem                0.06540027           0.10151699          0.2371308
## Mulher               0.06595625           0.08940651          0.2435124
##        A005Superior incompleto A005Superior completo
## Homem               0.05653695            0.09058213
## Mulher              0.05917339            0.12033296
##        A005Pós-graduação, mestrado ou doutorado se.A005Sem instrução
## Homem                                0.02162928         0.0008130148
## Mulher                               0.02973113         0.0007939639
##        se.A005Fundamental incompleto se.A005Fundamental completa
## Homem                    0.001658020                0.0008207212
## Mulher                   0.001548506                0.0008300213
##        se.A005Médio incompleto se.A005Médio completo se.A005Superior incompleto
## Homem             0.0009892119           0.001496970               0.0008482823
## Mulher            0.0008810969           0.001390631               0.0007893702
##        se.A005Superior completo se.A005Pós-graduação, mestrado ou doutorado
## Homem               0.001326250                                0.0006608044
## Mulher              0.001399994                                0.0006509593

Também pode ser utilizado para estimar a renda média normalmente recebida em dinheiro por Unidade da Federação:

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

##                                      UF   C01012        se
## Rondônia                       Rondônia 1869.561  53.70151
## Acre                               Acre 1868.272  76.03116
## Amazonas                       Amazonas 1758.269  86.03340
## Roraima                         Roraima 2044.991 128.13089
## Pará                               Pará 1764.165  72.23597
## Amapá                             Amapá 1687.106  96.42551
## Tocantins                     Tocantins 1897.691  74.52361
## Maranhão                       Maranhão 1425.784  38.41231
## Piauí                             Piauí 1679.268  83.91220
## Ceará                             Ceará 1646.840  55.56131
## Rio Grande do Norte Rio Grande do Norte 1887.514 118.78468
## Paraíba                         Paraíba 1799.818  90.11591
## Pernambuco                   Pernambuco 1811.761  85.73306
## Alagoas                         Alagoas 1554.754  66.53707
## Sergipe                         Sergipe 1739.017  97.80963
## Bahia                             Bahia 1567.659  71.60076
## Minas Gerais               Minas Gerais 2004.190  44.02823
## Espírito Santo           Espírito Santo 2129.544  72.46981
## Rio de Janeiro           Rio de Janeiro 2766.867  68.73415
## São Paulo                     São Paulo 2915.141  91.86012
## Paraná                           Paraná 2528.688  66.77902
## Santa Catarina           Santa Catarina 2472.574  40.33833
## Rio Grande do Sul     Rio Grande do Sul 2494.428  58.52093
## Mato Grosso do Sul   Mato Grosso do Sul 2263.340  71.00848
## Mato Grosso                 Mato Grosso 2293.420  87.25054
## Goiás                             Goiás 2102.786  61.45284
## Distrito Federal       Distrito Federal 4011.070 224.76406

É possível também calcular o intervalo de confiança para cada uma dessas estimativas:

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

É possível definir domínios que sejam cruzamentos de variáveis categóricas com a função interaction.
Na função svyby também é possível utilizar o argumento vartype=“cv” caso desejemos que no output apareça o coeficiente de variação ao invés do erro padrão da estimativa. Para estimar a taxa de afastamento do trabalho por motivo de quarentena, isolamento, distanciamento social ou férias coletivas por sexo e cor ou raça utilizando as funções svyratio e svyby e o respectivo cv:

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

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 5 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 63 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 4 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 46 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 18 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 143 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 23 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 108 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 156 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 58 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 163 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 88 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 4 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 40 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 2 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 35 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 154 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 45 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 164 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 51 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 35 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 3 strata have only one PSU in
## this subset.

## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 24 strata
## have only one PSU in this subset.

## Warning in `[.survey.design2`(design, !nas, ): 1 strata have only one PSU in
## this subset.

txafastquarentenaSexoRaca

##                 interaction(A003, A004)
## Homem.Branca               Homem.Branca
## Mulher.Branca             Mulher.Branca
## Homem.Preta                 Homem.Preta
## Mulher.Preta               Mulher.Preta
## Homem.Amarela             Homem.Amarela
## Mulher.Amarela           Mulher.Amarela
## Homem.Parda                 Homem.Parda
## Mulher.Parda               Mulher.Parda
## Homem.Indígena           Homem.Indígena
## Mulher.Indígena         Mulher.Indígena
## Homem.Ignorado           Homem.Ignorado
## Mulher.Ignorado         Mulher.Ignorado
##                 C003 == "Estava em quarentena, isolamento, distanciamento social ou férias coletivas"/C002 == "Sim"
## Homem.Branca                                                                                              0.7929560
## Mulher.Branca                                                                                             0.8163735
## Homem.Preta                                                                                               0.8264842
## Mulher.Preta                                                                                              0.8494372
## Homem.Amarela                                                                                             0.8573920
## Mulher.Amarela                                                                                            0.8489387
## Homem.Parda                                                                                               0.8060785
## Mulher.Parda                                                                                              0.8462040
## Homem.Indígena                                                                                            0.7526980
## Mulher.Indígena                                                                                           0.9010053
## Homem.Ignorado                                                                                            1.0000000
## Mulher.Ignorado                                                                                           1.0000000
##                          cv
## Homem.Branca    0.008962047
## Mulher.Branca   0.007705869
## Homem.Preta     0.015419888
## Mulher.Preta    0.012504395
## Homem.Amarela   0.042440989
## Mulher.Amarela  0.047796216
## Homem.Parda     0.007676511
## Mulher.Parda    0.006359161
## Homem.Indígena  0.108320049
## Mulher.Indígena 0.046924872
## Homem.Ignorado  0.000000000
## Mulher.Ignorado 0.000000000

Gráficos para dados amostrais

O pacote survey possui funções específicas para gerar gráficos que incorporem os pesos amostrais das observações. Para a criação dos gráficos, é utilizado o mesmo objeto de plano amostral utilizado para estimar as quantidades populacionais.

Entre os gráficos disponíveis no pacote estão o histograma, boxplot e gráficos de dispersão.

Histograma

A função svyhist permite a criação de histogramas para variáveis numéricas que consideram os pesos amostrais para computar as frequências ou densidades das barras.
A sintaxe básica é semelhante as de estimação apresentadas anteriormente. Além disto, possui os mesmos parâmetros que a função hist: breaks, xlab, main, […].
Para estimativas das frequências absolutas, deve-se utilizar freq=TRUE. Abaixo, um exemplo do histograma do número de horas normalmente trabalhadas para cada um dos casos:

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

svyhist(formula=~as.numeric(C008), design=dadosPNADCOVID19, freq=TRUE, main="Histograma", xlab="Número de Horas Normalmente Trabalhadas")

Boxplot

Para a construção de boxplots que considerem os pesos amostrais, a função é svyboxplot.

A sintaxe dela difere um pouco das demais. É necessário declarar a variável numérica, sucedida por um ~ e a variável de agrupamento do boxplot. Caso não deseje usar grupos, basta colocar o número 1. Além disto, há a opção de plotar apenas os outliers mais extremos ou todos, através do argumento all.outliers=TRUE. Exemplos:

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

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

## Warning in `[.survey.design2`(design, nas == 0, ): 3 strata have only one PSU in
## this subset.

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

## Warning in `[.survey.design2`(design, nas == 0, ): 3 strata have only one PSU in
## this subset.

Gráficos de dispersão

Os gráficos de dispersão usuais não são muito úteis para dados provindo de amostras complexas, pois não conseguem diferenciar o peso de cada ponto representado no gráfico. Para isso, a função svyplot tem algumas opções de construção de gráficos de dispersão que representam esses pesos através do argumento style.

Uma opção é a utilização de style=“bubble”, que o tamanho do ponto no gráfico representa o peso das observações naquele ponto.

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

Outra opção é a style=“transparent”, onde regiões com menor peso são mais transparentes e com maior peso, mais sólidas:

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

Modelagem com pacote `survey`

No pacote survey ainda há a possibilidade de realizar testes de hipóteses e estimar parâmetros de modelos, considerando o plano amostral complexo.

Testes de hipóteses

Os testes de hipóteses incluídos no pacote survey incluem teste-t para médias, teste qui-quadrado e teste de postos.

Abaixo um exemplo do teste-t para diferenças de médias de rendimentos normalmente recebidos entre sexos:

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

## 
##  Design-based t-test
## 
## data:  C01012 ~ A003
## t = -19.732, df = 14078, p-value < 0.00000000000000022
## 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

Outro exemplo para testar diferenças entre médias de horas normalmente trabalhadas, entre pessoas que realizaram e não realizaram trabalho ou algum bico na semana de referência:

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

## 
##  Design-based t-test
## 
## data:  as.numeric(C008) ~ C001
## t = -21.681, df = 14099, p-value < 0.00000000000000022
## 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

Modelos lineares

Além de testes de hipóteses, é possível estimar modelos lineares generalizados para dados amostrais complexos através da função svyglm.

Regressão linear

Para regressão linear simples ou múltipla, basta descrever a fórmula do modelo desejado na função svyglm.

No exemplo abaixo, utilizamos um modelo de regressão com rendimento normalmente recebido em dinheiro como variável dependente e escolaridade, cor ou raça e idade como variáveis explicativas. Para obter as estatísticas do modelo, pode ser utilizada a função summary, assim como feito para a regressão linear convencional.

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

## 
## Call:
## svyglm(formula = C01012 ~ A005 + A004 + A002, design = dadosPNADCOVID19)
## 
## Survey design:
## survey::postStratify(design = data_prior, strata = ~posest, population = popc.types)
## 
## Coefficients:
##                                           Estimate Std. Error t value
## (Intercept)                               -482.700     70.679  -6.830
## A005Fundamental incompleto                 484.605     37.630  12.878
## A005Fundamental completa                   823.335     43.204  19.057
## A005Médio incompleto                       997.392     44.743  22.291
## A005Médio completo                        1219.472     41.736  29.219
## A005Superior incompleto                   1738.514     59.486  29.225
## A005Superior completo                     3240.245     65.525  49.451
## A005Pós-graduação, mestrado ou doutorado  5845.701    228.253  25.611
## A004Preta                                 -549.097     38.149 -14.394
## A004Amarela                                127.499    181.763   0.701
## A004Parda                                 -510.387     29.768 -17.146
## A004Indígena                              -630.133    110.201  -5.718
## A004Ignorado                             -1235.340    438.179  -2.819
## A002                                        34.720      1.039  33.404
##                                                      Pr(>|t|)    
## (Intercept)                                  0.00000000000887 ***
## A005Fundamental incompleto               < 0.0000000000000002 ***
## A005Fundamental completa                 < 0.0000000000000002 ***
## A005Médio incompleto                     < 0.0000000000000002 ***
## A005Médio completo                       < 0.0000000000000002 ***
## A005Superior incompleto                  < 0.0000000000000002 ***
## A005Superior completo                    < 0.0000000000000002 ***
## A005Pós-graduação, mestrado ou doutorado < 0.0000000000000002 ***
## A004Preta                                < 0.0000000000000002 ***
## A004Amarela                                           0.48303    
## A004Parda                                < 0.0000000000000002 ***
## A004Indígena                                 0.00000001099712 ***
## A004Ignorado                                          0.00482 ** 
## A002                                     < 0.0000000000000002 ***
## ---
## 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

Além disso, é possível computar intervalos de confiança para os parâmetros:

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

Regressão logística

Outro modelo pertencente a classe de modelos generalizados é a Regressão Logística. Para utilizar este modelo, basta colocar o argumento family=“binomial”.

No exemplo abaixo modelamos o afastamento do trabalho pelo sexo, cor ou raça e idade:

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

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

summary(object=modeloLog)

## 
## Call:
## svyglm(formula = C002 ~ A003 + A004 + A002, design = dadosPNADCOVID19, 
##     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 < 0.0000000000000002 ***
## A003Mulher    0.297449   0.015164  19.615 < 0.0000000000000002 ***
## 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 < 0.0000000000000002 ***
## ---
## 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

Sendo possível também computar intervalos de confiança para os parâmetros:

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

Análise de concentração de renda com pacote convey

O pacote convey permite estimar diversas medidas de concentração de renda para dados provenientes de pesquisas com planos amostrais complexos.

Este pacote segue uma sintaxe bem próxima à sintaxe do pacote survey, sendo possível utilizar funções do pacote survey em objetos do pacote convey. Uma extensa documentação sobre este pacote e sobre o processo de estimação dos índices pode ser encontrada no GitHub dos autores: https://guilhermejacob.github.io/context/index.html.

Para a utilização do pacote convey, primeiramente é necessário utilizar a função convey_prep, que transforma o objeto do plano amostral do survey no objeto que o pacote convey utiliza para as estimações:

library(convey)
dadosPNADCOVID19 <- convey_prep(design=dadosPNADCOVID19)

Aqui, apresentaremos apenas exemplo de duas das diversas medidas disponíveis no pacote: O índice de Gini e a Curva de Lorenz.

Índice de Gini

Para estimarmos o índice de Gini, basta utilizar a função svygini com a variável de renda desejada. Neste exemplo, estimamos o índice de Gini do Brasil para a renda normalmente recebida em dinheiro:

giniBR <- svygini(formula=~C01012, design=dadosPNADCOVID19, na.rm=TRUE)
giniBR

##           gini     SE
## C01012 0.47399 0.0043

Assim como nas funções do pacote survey, é possível, por exemplo, estimar o coeficiente de variação desta estimativa:

cv(object=giniBR)

##             C01012
## C01012 0.009041141

Também é possível utilizar a função svyby para estimar o índice de Gini da renda normalmente recebida em dinheiro por Unidade da Federação:

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

##                                      UF    C01012          se
## Rondônia                       Rondônia 0.4074291 0.011324144
## Acre                               Acre 0.4390348 0.015582367
## Amazonas                       Amazonas 0.4766059 0.015077890
## Roraima                         Roraima 0.4592258 0.021962424
## Pará                               Pará 0.4851257 0.015980674
## Amapá                             Amapá 0.4606922 0.016558531
## Tocantins                     Tocantins 0.4461214 0.013853156
## Maranhão                       Maranhão 0.4433619 0.010244075
## Piauí                             Piauí 0.4850488 0.016660073
## Ceará                             Ceará 0.4668615 0.011721721
## Rio Grande do Norte Rio Grande do Norte 0.4987877 0.018799000
## Paraíba                         Paraíba 0.4790974 0.016442643
## Pernambuco                   Pernambuco 0.4876195 0.017703430
## Alagoas                         Alagoas 0.4425971 0.016298338
## Sergipe                         Sergipe 0.4953125 0.017704319
## Bahia                             Bahia 0.4835445 0.017663378
## Minas Gerais               Minas Gerais 0.4328628 0.008842483
## Espírito Santo           Espírito Santo 0.4540853 0.011532936
## Rio de Janeiro           Rio de Janeiro 0.4771287 0.008266112
## São Paulo                     São Paulo 0.4754268 0.011403350
## Paraná                           Paraná 0.4378381 0.009028997
## Santa Catarina           Santa Catarina 0.3828977 0.006559015
## Rio Grande do Sul     Rio Grande do Sul 0.4431499 0.009080239
## Mato Grosso do Sul   Mato Grosso do Sul 0.4172560 0.011828405
## Mato Grosso                 Mato Grosso 0.4067868 0.014847231
## Goiás                             Goiás 0.4186365 0.011400029
## Distrito Federal       Distrito Federal 0.5243862 0.008708074

E estimar intervalos de confiança para cada UF:

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

Curva de Lorenz

A Curva de Lorenz é um gráfico utilizado para relacionar a distribuição relativa de renda pelas pessoas. A área entre essa curva e a reta identidade, é uma das formas de definir o coeficiente de Gini. No pacote convey, é possível fazer o gráfico desta curva com a função svylorenz. Os quantis da população, para os quais a renda acumulada será plotada no gráfico, são definidos no argumento quantiles.

Exemplo para a Curva de Lorenz da renda normalmente recebida em dinheiro:

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

Instituto Brasileiro de Geografia e Estatística - pacotesipd@ibge.gov.br ↩︎

Análise de microdados da PNAD COVID19

Com os pacotes `COVIDIBGE` e `survey`

Gabriel Assuncao¹

01 de fevereiro de 2024

O pacote COVIDIBGE

Instalação do pacote

Importação de microdados

Importação online

Microdados mensais

Importação offline

Análise com pacote survey

Estimando totais

Variáveis numéricas

Variáveis categóricas

Estimando médias

Estimando proporções

Estimando razões

Estimando medianas e quantis

Estimação para variáveis de rendimento deflacionadas

Estimação para um domínio

Estimação para vários domínios

Gráficos para dados amostrais

Histograma

Boxplot

Gráficos de dispersão

Modelagem com pacote `survey`

Testes de hipóteses

Modelos lineares

Regressão linear

Regressão logística

Análise de concentração de renda com pacote convey

Índice de Gini

Curva de Lorenz

Análise de microdados da PNAD COVID19

Com os pacotes COVIDIBGE e survey

Gabriel Assuncao1

01 de fevereiro de 2024

O pacote COVIDIBGE

Instalação do pacote

Importação de microdados

Importação online

Microdados mensais

Importação offline

Análise com pacote survey

Estimando totais

Variáveis numéricas

Variáveis categóricas

Estimando médias

Estimando proporções

Estimando razões

Estimando medianas e quantis

Estimação para variáveis de rendimento deflacionadas

Estimação para um domínio

Estimação para vários domínios

Gráficos para dados amostrais

Histograma

Boxplot

Gráficos de dispersão

Modelagem com pacote survey

Testes de hipóteses

Modelos lineares

Regressão linear

Regressão logística

Análise de concentração de renda com pacote convey

Índice de Gini

Curva de Lorenz

Com os pacotes `COVIDIBGE` e `survey`

Gabriel Assuncao¹

Modelagem com pacote `survey`