Análise de microdados da PNAD Contínua
Com os pacotes PNADcIBGE e
survey
Gabriel Assuncao1
01 de abril de 2026
O pacote PNADcIBGE
O pacote PNADcIBGE foi desenvolvido para facilitar o download, importação e análise dos dados amostrais da Pesquisa Nacional por Amostra de Domicílios Contínua realizada pelo Instituto Brasileiro de Geografia e Estatística - IBGE.
A PNAD Contínua possui três tipos de microdados:- Trimestral, que contém a parte básica investigada pela pesquisa, contendo variáveis conjunturais de mercado de trabalho referentes a um trimestre civil;
- Anual acumulados em determinada visita, que contém temas e tópicos suplementares pesquisados ao longo do ano em determinada visita;
- Anual concentrados em determinado trimestre, que contém temas e tópicos suplementares pesquisados em trimestres específicos do ano.
Além disso, a partir de 2021 foram incluídos módulos aplicados
somente para um morador selecionado na PNAD Contínua, para temas e
tópicos suplementares específicos.
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/PNADcIBGE.
É
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("PNADcIBGE")Antes de utilizar o pacote é necessário carregá-lo no R através do comando:
library(PNADcIBGE)Importação de microdados
Para a importação dos microdados da PNAD Contínua, 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_pnadc 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 microdados trimestrais e anuais (acumulados em determinada visita
ou concentrados em determinado trimestre).
help("get_pnadc")Microdados trimestrais
A importação dos microdados trimestrais através da função
get_pnadc é muito simples. Basta indicar o ano e o
trimestre dos microdados desejados nos argumentos da função. Abaixo um
exemplo de como importar os microdados do 4º trimestre de 2017:
dadosPNADc <- get_pnadc(year=2017, quarter=4)Apenas com esse comando, os microdados da pesquisa já estão
disponíveis no objeto dadosPNADc para análise.
year e quarter, que
indicam, respectivamente, o ano e o trimestre dos microdados a serem
baixados, a função get_pnadc 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. Odefaulté rotulá-los; -
deflator: Um argumento lógico que indica se as variáveis para deflacionamento serão incluídas nos microdados. Odefaulté 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. Casodesign=FALSE, a função retorna apenas um data-frame com os microdados. É altamente recomendado que mantenha essa opção comoTRUE, 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 Condição em relação à força de trabalho (VD4001) e
Condição de ocupação (VD4002) do 4º trimestre de 2017,
pode utilizar o argumento vars para selecionar apenas estas
variáveis:
dadosPNADc <- get_pnadc(year=2017, quarter=4, vars=c("VD4001","VD4002"))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.
dadosPNADc## Call: svrepdesign.default(data = data_pnadc, weight = ~V1028, type = "bootstrap",
## repweights = "V1028[0-9]+", mse = TRUE, replicates = length(sprintf("V1028%03d",
## seq(1:200))), df = length(sprintf("V1028%03d", seq(1:200))))
## Survey bootstrap with 200 replicates and MSE variances.class(dadosPNADc)## [1] "svyrep.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:
dadosPNADc_brutos <- get_pnadc(year=2017, quarter=4, vars=c("VD4001","VD4002"), design=FALSE)O Objeto resultante é um data-frame com as variáveis selecionadas e as variáveis relacionadas ao processo de amostragem:
dadosPNADc_brutos## # A tibble: 561,288 × 219
## Ano Trimestre UF UPA Estrato V1008 V1014 V1027 V1028 V1029 V1033
## <chr> <chr> <fct> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2017 4 Rondônia 11000… 1110011 01 06 121. 129. 488201 6.36e6
## 2 2017 4 Rondônia 11000… 1110011 01 06 121. 129. 488201 8.39e6
## 3 2017 4 Rondônia 11000… 1110011 03 06 121. 140. 488201 7.43e6
## 4 2017 4 Rondônia 11000… 1110011 03 06 121. 140. 488201 8.34e6
## 5 2017 4 Rondônia 11000… 1110011 03 06 121. 140. 488201 7.43e6
## 6 2017 4 Rondônia 11000… 1110011 04 06 121. 161. 488201 8.39e6
## 7 2017 4 Rondônia 11000… 1110011 04 06 121. 161. 488201 8.38e6
## 8 2017 4 Rondônia 11000… 1110011 05 06 121. 128. 488201 8.06e6
## 9 2017 4 Rondônia 11000… 1110011 05 06 121. 128. 488201 6.36e6
## 10 2017 4 Rondônia 11000… 1110011 05 06 121. 128. 488201 1.00e7
## # ℹ 561,278 more rows
## # ℹ 208 more variables: posest <chr>, posest_sxi <chr>, V2003 <chr>,
## # VD4001 <fct>, VD4002 <fct>, V1028001 <dbl>, V1028002 <dbl>, V1028003 <dbl>,
## # V1028004 <dbl>, V1028005 <dbl>, V1028006 <dbl>, V1028007 <dbl>,
## # V1028008 <dbl>, V1028009 <dbl>, V1028010 <dbl>, V1028011 <dbl>,
## # V1028012 <dbl>, V1028013 <dbl>, V1028014 <dbl>, V1028015 <dbl>,
## # V1028016 <dbl>, V1028017 <dbl>, V1028018 <dbl>, V1028019 <dbl>, …Também é possível baixar os dados sem incluir os rótulos dos níveis,
através do argumento labels.
dadosPNADc_brutos_sem <- get_pnadc(year=2017, quarter=4, vars=c("VD4001","VD4002"), labels=FALSE, design=FALSE)Perceba que agora os níveis das categorias são representados por números:
dadosPNADc_brutos_sem## # A tibble: 561,288 × 219
## Ano Trimestre UF UPA Estrato V1008 V1014 V1027 V1028 V1029 V1033
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 2017 4 11 110000016 1110011 01 06 121. 129. 488201 6.36e6
## 2 2017 4 11 110000016 1110011 01 06 121. 129. 488201 8.39e6
## 3 2017 4 11 110000016 1110011 03 06 121. 140. 488201 7.43e6
## 4 2017 4 11 110000016 1110011 03 06 121. 140. 488201 8.34e6
## 5 2017 4 11 110000016 1110011 03 06 121. 140. 488201 7.43e6
## 6 2017 4 11 110000016 1110011 04 06 121. 161. 488201 8.39e6
## 7 2017 4 11 110000016 1110011 04 06 121. 161. 488201 8.38e6
## 8 2017 4 11 110000016 1110011 05 06 121. 128. 488201 8.06e6
## 9 2017 4 11 110000016 1110011 05 06 121. 128. 488201 6.36e6
## 10 2017 4 11 110000016 1110011 05 06 121. 128. 488201 1.00e7
## # ℹ 561,278 more rows
## # ℹ 208 more variables: posest <chr>, posest_sxi <chr>, V2003 <chr>,
## # VD4001 <chr>, VD4002 <chr>, V1028001 <dbl>, V1028002 <dbl>, V1028003 <dbl>,
## # V1028004 <dbl>, V1028005 <dbl>, V1028006 <dbl>, V1028007 <dbl>,
## # V1028008 <dbl>, V1028009 <dbl>, V1028010 <dbl>, V1028011 <dbl>,
## # V1028012 <dbl>, V1028013 <dbl>, V1028014 <dbl>, V1028015 <dbl>,
## # V1028016 <dbl>, V1028017 <dbl>, V1028018 <dbl>, V1028019 <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.
Microdados anuais acumulados em determinada visita
Além dos microdados trimestrais, a função get_pnadc
também permite a importação dos microdados anuais.
A relação dos
temas e tópicos pesquisados ao longo do ano em determinada visita, pode
ser acessada na página
da pesquisa.
- Características Adicionais do Mercado de Trabalho;
- Características Gerais dos Moradores;
- Características Gerais dos Domicílios;
- Rendimentos de Outras Fontes.
Para a importação dos microdados anuais de uma entrevista basta
informar o número da entrevista desejada no argumento
interview:
dadosPNADc_anual_visita <- get_pnadc(year=2017, interview=1)
dadosPNADc_anual_visita## Call: svrepdesign.default(data = data_pnadc, weight = ~V1032, type = "bootstrap",
## repweights = "V1032[0-9]+", mse = TRUE, replicates = length(sprintf("V1032%03d",
## seq(1:200))), df = length(sprintf("V1032%03d", seq(1:200))))
## Survey bootstrap with 200 replicates and MSE variances.O objeto retornado é da classe survey.design ou
svyrep.design, e os argumentos vars,
labels, deflator, design,
reload, curlopts e savedir podem
ser utilizados da mesma forma dos microdados trimestrais.
Sendo ainda importante ressaltar que para o caso de inserção das
variáveis de deflacionamento, para este tipo de microdados, existe um
argumento adicional defyear que é essencial para realização
desta etapa, que se encontra mais detalhado no documento disponível no
site oficial do IBGE para auxílio no processo de deflacionamento
das variáveis de rendimento.
Microdados anuais concentrados em determinado trimestre
Conforme mencionado, além dos microdados trimestrais, a função
get_pnadc também permite a importação dos microdados
anuais.
A relação dos temas e tópicos suplementares pesquisados em
trimestres específicos do ano, pode ser acessada na página
da pesquisa.
- Educação.
Para a importação dos microdados anuais de um trimestre basta
informar o número do trimestre desejado no argumento
topic:
dadosPNADc_anual_trimestre <- get_pnadc(year=2017, topic=2)
dadosPNADc_anual_trimestre## Call: svrepdesign.default(data = data_pnadc, weight = ~V1028, type = "bootstrap",
## repweights = "V1028[0-9]+", mse = TRUE, replicates = length(sprintf("V1028%03d",
## seq(1:200))), df = length(sprintf("V1028%03d", seq(1:200))))
## Survey bootstrap with 200 replicates and MSE variances.O objeto retornado é da classe survey.design ou
svyrep.design, e os argumentos vars,
labels, deflator, design,
reload, curlopts e savedir podem
ser utilizados da mesma forma dos microdados trimestrais.
Sendo ainda importante ressaltar que para o caso de inserção das
variáveis de deflacionamento, para este tipo de microdados, existem dois
argumentos adicionais defyear e defperiod que
são essenciais para realização desta etapa, que se encontra mais
detalhado no documento disponível no site oficial do IBGE para auxílio
no processo de deflacionamento
das variáveis de rendimento.
Microdados de temas e tópicos suplementares de morador selecionado
Além dos tipos de microdados apresentados anteriormente, a função
get_pnadc também permite a importação dos microdados de
temas e tópicos suplementares de morador selecionado.
Neste exemplo, mostraremos como importar os microdados do módulo de Sensação de Segurança aplicado para morador selecionado da PNAD Contínua nos microdados anuais concentrados no 4º trimestre de 2021.
Para a importação dos microdados de temas e tópicos suplementares de
morador selecionado basta utilizar o valor TRUE no
argumento selected:
dadosPNADc_anual_trimestre_morador_selecionado <- get_pnadc(year=2021, topic=4, selected=TRUE)
dadosPNADc_anual_trimestre_morador_selecionado## Call: svrepdesign.default(data = data_pnadc, weight = ~V1036, type = "bootstrap",
## repweights = "V1036[0-9]+", mse = TRUE, replicates = length(sprintf("V1036%03d",
## seq(1:200))), df = length(sprintf("V1036%03d", seq(1:200))))
## Survey bootstrap with 200 replicates and MSE variances.O objeto retornado é da classe survey.design ou
svyrep.design, e os argumentos vars,
labels, deflator, design,
reload, curlopts e savedir podem
ser utilizados da mesma forma dos microdados trimestrais.
Já sobre o processo de inserção das variáveis de deflacionamento, devem ser seguidas as mesmas orientações que o tipo de microdados em que o tema e tópico suplementar de morador selecionado está inserido. Então, para o exemplo fornecido deve se levar em consideração o mesmo processo dos microdados anuais concentrados em determinado trimestre.
Importação offline
Caso não seja possível utilizar a importação online através da funçãoget_pnadc, é possível criar o mesmo objeto para
análise com arquivos que estejam em disco. Para isto, são utilizadas
quatro funções:
-
read_pnadc: Para a leitura do arquivo .txt dos microdados; -
pnadc_labeller: Opcional. Coloca os rótulos dos níveis nas variáveis categóricas; -
pnadc_deflator: Opcional. Insere as variáveis para deflacionamento nos microdados; -
pnadc_design: Cria o objeto do plano amostral para a análise com o pacotesurvey.
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_pnadc são utilizados os arquivos de
texto contendo os microdados e o input para SAS.
Exemplo de leitura para o 4º trimestre de 2017:
dadosPNADc <- read_pnadc(microdata="PNADC_042017.txt", input_txt="input_PNADC_trimestral.txt")Assim como na importação online, pode ser utilizado o
argumento vars para definir quais variáveis serão
lidas.
A função pnadc_labeller utiliza o objeto criado pela
função read_pnadc e o arquivo de dicionário das variáveis
presente na documentação da pesquisa.
dadosPNADc <- pnadc_labeller(data_pnadc=dadosPNADc, dictionary.file="dicionario_PNADC_microdados_trimestral.xls")A utilização de pnadc_labeller é opcional e não
interfere nos resultados obtidos.
A função pnadc_deflator utiliza o objeto criado pela
função read_pnadc e o arquivo de deflatores presente na
documentação da pesquisa.
dadosPNADc <- pnadc_deflator(data_pnadc=dadosPNADc, deflator.file="deflator_PNADC_trimestral.xls")A utilização de pnadc_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 pnadc_design aplicada no objeto que
contém os microdados:
dadosPNADc <- pnadc_design(data_pnadc=dadosPNADc)O objeto resultante é o mesmo retornado pela função
get_pnadc e pode ser utilizado para análises com o pacote
survey.
A importação offline é feita da mesma forma tanto para microdados trimestrais, quanto para microdados anuais (acumulados em determinada visita ou concentrados em determinado trimestre), e também para os microdados de temas e tópicos suplementares de morador selecionado (sendo que para este caso específico é necessário realizar o ajuste das bases com o recorte do morador selecionado e removendo os fatores de expansão que não correspondem à parte de morador selecionado).
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 Contínua do 4º trimestre de 2017, com algumas variáveis selecionadas de acordo com a tabela abaixo:
| Variável | Descrição |
|---|---|
| UF | Unidade da Federação |
| V2001 | Número de pessoas no domicílio |
| V2005 | Condição no domicílio |
| V2007 | Sexo |
| V2009 | Idade do morador na data de referência |
| V2010 | Cor ou raça |
| V3007 | Já concluiu algum outro curso de graduação? |
| VD3004 | Nível de instrução mais elevado alcançado (pessoas de 5 anos ou mais de idade) |
| VD4001 | Condição em relação à força de trabalho na semana de referência para pessoas de 14 anos ou mais de idade |
| VD4002 | Condição de ocupação na semana de referência para pessoas de 14 anos ou mais de idade |
| 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) |
| VD4035 | Horas efetivamente trabalhadas na semana de referência em todos os trabalhos para pessoas de 14 anos ou mais de idade |
Importando os dados:
variaveis_selecionadas <- c("UF","V2001","V2005","V2007","V2009","V2010","V3007","VD3004","VD4001","VD4002","VD4020","VD4035")
dadosPNADc <- get_pnadc(year=2017, quarter=4, vars=variaveis_selecionadas)O objeto dadosPNADc 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 é asvytotal. 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 mensal efetiva:
totalrenda <- svytotal(x=~VD4020, design=dadosPNADc, na.rm=TRUE)
totalrenda## total SE
## VD4020 193667566247 2850406391Além da estimativa do total da renda efetiva, 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)## VD4020
## 0.01471804confint(object=totalrenda)## 2.5 % 97.5 %
## VD4020 188080872378 199254260115confint(object=totalrenda, level=0.99)## 0.5 % 99.5 %
## VD4020 186325405937 201009726557Variáveis categóricas
Também é possível estimar totais populacionais de categorias, utilizando variáveis categóricas, como o sexo:
totalsexo <- svytotal(x=~V2007, design=dadosPNADc, na.rm=TRUE)
totalsexo## total SE
## V2007Homem 100158028 0.0015
## V2007Mulher 104841068 0.0016E estimar o total de mais de uma variável categórica no mesmo código,
separando com o operador +:
totalsexoraca <- svytotal(x=~V2007+V2010, design=dadosPNADc, na.rm=TRUE)
totalsexoraca## total SE
## V2007Homem 100158028 0.0015
## V2007Mulher 104841068 0.0016
## V2010Branca 88337400 402169.4110
## V2010Preta 17809091 197222.1726
## V2010Amarela 1165276 63090.6709
## V2010Parda 97118150 386763.9396
## V2010Indígena 548942 29839.4417
## V2010Ignorado 20238 5584.5975També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(V2007,V2010), design=dadosPNADc, na.rm=TRUE)
ftable(x=totalsexoEraca)## A B
## interaction(V2007, V2010)Homem.Branca 42451217.472 217179.138
## interaction(V2007, V2010)Mulher.Branca 45886182.314 224683.585
## interaction(V2007, V2010)Homem.Preta 9004498.809 106021.159
## interaction(V2007, V2010)Mulher.Preta 8804591.952 115580.825
## interaction(V2007, V2010)Homem.Amarela 526629.500 32865.509
## interaction(V2007, V2010)Mulher.Amarela 638646.616 36842.319
## interaction(V2007, V2010)Homem.Parda 47914827.987 210288.245
## interaction(V2007, V2010)Mulher.Parda 49203321.635 213324.056
## interaction(V2007, V2010)Homem.Indígena 250264.613 16092.220
## interaction(V2007, V2010)Mulher.Indígena 298677.441 18490.685
## interaction(V2007, V2010)Homem.Ignorado 10589.619 2857.812
## interaction(V2007, V2010)Mulher.Ignorado 9648.042 3563.405A 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 efetiva pode ser
facilmente reescrito para médias:
mediarenda <- svymean(x=~VD4020, design=dadosPNADc, na.rm=TRUE)
mediarenda## mean SE
## VD4020 2180.2 30.954E podemos calcular coeficientes de variação e intervalos de confiança da mesma forma:
cv(object=mediarenda)## VD4020
## 0.0141979confint(object=mediarenda)## 2.5 % 97.5 %
## VD4020 2119.493 2240.83Estimando 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=~V2007, design=dadosPNADc, na.rm=TRUE)
propsexo## mean SE
## V2007Homem 0.48858 0
## V2007Mulher 0.51142 0Da mesma forma, podemos estimar a proporção de mais de uma variável:
propsexoraca <- svymean(x=~V2007+V2010, design=dadosPNADc, na.rm=TRUE)
propsexoraca## mean SE
## V2007Homem 0.488577901 0.0000
## V2007Mulher 0.511422099 0.0000
## V2010Branca 0.430916046 0.0020
## V2010Preta 0.086873997 0.0010
## V2010Amarela 0.005684299 0.0003
## V2010Parda 0.473749160 0.0019
## V2010Indígena 0.002677778 0.0001
## V2010Ignorado 0.000098721 0.0000E estimando a proporção de um cruzamento de duas ou mais variáveis
com a função interaction:
propsexoEraca <- svymean(x=~interaction(V2007,V2010), design=dadosPNADc, na.rm=TRUE)
ftable(x=propsexoEraca)## A B
## interaction(V2007, V2010)Homem.Branca 0.20708002279 0.00105941510
## interaction(V2007, V2010)Mulher.Branca 0.22383602274 0.00109602232
## interaction(V2007, V2010)Homem.Preta 0.04392457813 0.00051717867
## interaction(V2007, V2010)Mulher.Preta 0.04294941843 0.00056381139
## interaction(V2007, V2010)Homem.Amarela 0.00256893572 0.00016032026
## interaction(V2007, V2010)Mulher.Amarela 0.00311536308 0.00017971942
## interaction(V2007, V2010)Homem.Parda 0.23373189893 0.00102580084
## interaction(V2007, V2010)Mulher.Parda 0.24001726151 0.00104060974
## interaction(V2007, V2010)Homem.Indígena 0.00122080837 0.00007849898
## interaction(V2007, V2010)Mulher.Indígena 0.00145696955 0.00009019886
## interaction(V2007, V2010)Homem.Ignorado 0.00005165691 0.00001394061
## interaction(V2007, V2010)Mulher.Ignorado 0.00004706383 0.00001738254Neste 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 desocupação divulgada pelo IBGE: ela é a razão entre o total de pessoas desocupadas pelo total de pessoas na força de trabalho.
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 é:
txdesocup <- svyratio(numerator=~(VD4002=="Pessoas desocupadas"), denominator=~(VD4001=="Pessoas na força de trabalho"), design=dadosPNADc, na.rm=TRUE)
txdesocup## Ratio estimator: svyratio.svyrep.design(numerator = ~(VD4002 == "Pessoas desocupadas"),
## denominator = ~(VD4001 == "Pessoas na força de trabalho"),
## design = dadosPNADc, na.rm = TRUE)
## Ratios=
## VD4001 == "Pessoas na força de trabalho"
## VD4002 == "Pessoas desocupadas" 0.1188146
## SEs=
## [,1]
## [1,] 0.001161042Cálculos de coeficiente de variação e intervalos de confiança para essa taxa:
cv(object=txdesocup)## VD4001 == "Pessoas na força de trabalho"
## VD4002 == "Pessoas desocupadas" 0.009771884confint(object=txdesocup)## 2.5 %
## VD4002 == "Pessoas desocupadas"/VD4001 == "Pessoas na força de trabalho" 0.116539
## 97.5 %
## VD4002 == "Pessoas desocupadas"/VD4001 == "Pessoas na força de trabalho" 0.1210902Estimando 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 efetiva, basta utilizar
quantiles=0.5:
medianarenda <- svyquantile(x=~VD4020, design=dadosPNADc, quantiles=0.5, ci=FALSE, na.rm=TRUE)
medianarenda## $VD4020
## 0.5
## [1,] 1300
##
## 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=~VD4020, design=dadosPNADc, quantiles=0.5, ci=TRUE, na.rm=TRUE)
medianarenda## $VD4020
## quantile ci.2.5 ci.97.5 se
## 0.5 1300 1300 1400 25.35553
##
## attr(,"hasci")
## [1] TRUE
## attr(,"class")
## [1] "newsvyquantile"SE(object=medianarenda)## VD4020
## 25.35553cv(object=medianarenda)## VD4020
## 0.01950425Alé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=~VD4020, design=dadosPNADc, quantiles=c(0.1,0.25,0.5,0.75,0.9), ci=FALSE, na.rm=TRUE)
quantisrenda## $VD4020
## 0.1 0.25 0.5 0.75 0.9
## [1,] 400 937 1300 2200 4100
##
## 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.
dadosPNADc$variables <- transform(dadosPNADc$variables, VD4020_real=VD4020*Efetivo)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=~VD4020_real, design=dadosPNADc, na.rm=TRUE)
totalrenda_real## total SE
## VD4020_real 292829615132 4368919904E para obtenção da média desta nova variável, basta utilizar a função
svymean:
mediarenda_real <- svymean(x=~VD4020_real, design=dadosPNADc, na.rm=TRUE)
mediarenda_real## mean SE
## VD4020_real 3296.5 47.467O 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.
Além disso, também está disponível um script de utilização no software R para estimação do Rendimento Domiciliar Per Capita (RDPC) real, ou seja, com a realização do processo de deflacionamento, a partir dos microdados anuais acumulados em determinada visita, através do endereço: https://github.com/Gabriel-Assuncao/PNADcIBGE-RDPC/blob/main/PNADcIBGE-RDPC.R.
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 efetiva de mulheres:
mediarendaM <- svymean(x=~VD4020, design=subset(dadosPNADc, V2007=="Mulher"), na.rm=TRUE)
mediarendaM## mean SE
## VD4020 1857.3 22.38Ou condicionais com desigualdade, como por exemplo estimar a taxa de desocupação para pessoas com idade igual ou superior a 25 anos:
txdesocup25 <- svyratio(numerator=~(VD4002=="Pessoas desocupadas"), denominator=~(VD4001=="Pessoas na força de trabalho"), design=subset(dadosPNADc, V2009>=25), na.rm=TRUE)
txdesocup25## Ratio estimator: svyratio.svyrep.design(numerator = ~(VD4002 == "Pessoas desocupadas"),
## denominator = ~(VD4001 == "Pessoas na força de trabalho"),
## design = subset(dadosPNADc, V2009 >= 25), na.rm = TRUE)
## Ratios=
## VD4001 == "Pessoas na força de trabalho"
## VD4002 == "Pessoas desocupadas" 0.08537702
## SEs=
## [,1]
## [1,] 0.00100833Alé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=~VD3004, design=subset(dadosPNADc, V2007=="Homem" & V2010=="Parda" & V2009>30), na.rm=TRUE)
nivelinstrHP30## mean SE
## VD3004Sem instrução e menos de 1 ano de estudo 0.106253 0.0015
## VD3004Fundamental incompleto ou equivalente 0.424568 0.0029
## VD3004Fundamental completo ou equivalente 0.088863 0.0016
## VD3004Médio incompleto ou equivalente 0.043282 0.0011
## VD3004Médio completo ou equivalente 0.239637 0.0025
## VD3004Superior incompleto ou equivalente 0.022390 0.0010
## VD3004Superior completo 0.075007 0.0016Alé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.
dadosPNADc_mulheres <- subset(dadosPNADc, V2007=="Mulher")
dadosPNADc_mulheres## Call: subset(dadosPNADc, V2007 == "Mulher")
## Survey bootstrap with 200 replicates and MSE variances.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 média de pessoas no domicílio:
mediapesdom <- svymean(x=~V2001, design=subset(dadosPNADc, V2005=="Pessoa responsável pelo domicílio"), na.rm=TRUE)
mediapesdom## mean SE
## V2001 3.0247 0.0045Estimaçã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.
- 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=~V2007, by=~VD3004, design=dadosPNADc, FUN=svymean, na.rm=TRUE)
freqSexoInstr## VD3004
## Sem instrução e menos de 1 ano de estudo Sem instrução e menos de 1 ano de estudo
## Fundamental incompleto ou equivalente Fundamental incompleto ou equivalente
## Fundamental completo ou equivalente Fundamental completo ou equivalente
## Médio incompleto ou equivalente Médio incompleto ou equivalente
## Médio completo ou equivalente Médio completo ou equivalente
## Superior incompleto ou equivalente Superior incompleto ou equivalente
## Superior completo Superior completo
## V2007Homem V2007Mulher se1
## Sem instrução e menos de 1 ano de estudo 0.4988095 0.5011905 0.0026488762
## Fundamental incompleto ou equivalente 0.5123485 0.4876515 0.0009531932
## Fundamental completo ou equivalente 0.5031311 0.4968689 0.0029573754
## Médio incompleto ou equivalente 0.5022278 0.4977722 0.0033824838
## Médio completo ou equivalente 0.4696166 0.5303834 0.0014753103
## Superior incompleto ou equivalente 0.4680143 0.5319857 0.0046047459
## Superior completo 0.4167383 0.5832617 0.0024194022
## se2
## Sem instrução e menos de 1 ano de estudo 0.0026488762
## Fundamental incompleto ou equivalente 0.0009531932
## Fundamental completo ou equivalente 0.0029573754
## Médio incompleto ou equivalente 0.0033824838
## Médio completo ou equivalente 0.0014753103
## Superior incompleto ou equivalente 0.0046047459
## Superior completo 0.0024194022Caso 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=~VD3004, by=~V2007, design=dadosPNADc, FUN=svymean, na.rm=TRUE)
freqInstrSexo## V2007 VD3004Sem instrução e menos de 1 ano de estudo
## Homem Homem 0.08328064
## Mulher Mulher 0.07936966
## VD3004Fundamental incompleto ou equivalente
## Homem 0.3889687
## Mulher 0.3511568
## VD3004Fundamental completo ou equivalente
## Homem 0.08569784
## Mulher 0.08027363
## VD3004Médio incompleto ou equivalente
## Homem 0.07214586
## Mulher 0.06782404
## VD3004Médio completo ou equivalente
## Homem 0.2308892
## Mulher 0.2473389
## VD3004Superior incompleto ou equivalente VD3004Superior completo
## Homem 0.04215686 0.09686084
## Mulher 0.04545183 0.12858507
## se1 se2 se3 se4 se5
## Homem 0.0007195373 0.001594971 0.0007749615 0.0007254636 0.001327625
## Mulher 0.0006874389 0.001348767 0.0007220973 0.0006780805 0.001183179
## se6 se7
## Homem 0.0006213938 0.001559208
## Mulher 0.0006797198 0.001532940Também pode ser utilizado para estimar a renda média efetiva por Unidade da Federação:
mediaRendaUF <- svyby(formula=~VD4020, by=~UF, design=dadosPNADc, FUN=svymean, na.rm=TRUE)
mediaRendaUF## UF VD4020 se
## Rondônia Rondônia 1799.668 62.35688
## Acre Acre 1615.063 85.78414
## Amazonas Amazonas 1916.221 185.61486
## Roraima Roraima 2114.125 139.34732
## Pará Pará 1415.618 46.81313
## Amapá Amapá 1983.087 155.04907
## Tocantins Tocantins 1804.530 82.18981
## Maranhão Maranhão 1267.407 69.68086
## Piauí Piauí 1338.874 67.14455
## Ceará Ceará 1438.505 60.19346
## Rio Grande do Norte Rio Grande do Norte 1492.551 95.66656
## Paraíba Paraíba 1619.595 99.76616
## Pernambuco Pernambuco 1621.235 74.81305
## Alagoas Alagoas 1351.000 49.05916
## Sergipe Sergipe 1493.832 80.81880
## Bahia Bahia 1536.066 117.68865
## Minas Gerais Minas Gerais 1930.869 40.55018
## Espírito Santo Espírito Santo 2037.891 63.22958
## Rio de Janeiro Rio de Janeiro 2339.637 53.47363
## São Paulo São Paulo 2821.676 113.08951
## Paraná Paraná 2368.360 62.08609
## Santa Catarina Santa Catarina 2440.376 38.47701
## Rio Grande do Sul Rio Grande do Sul 2369.333 56.15529
## Mato Grosso do Sul Mato Grosso do Sul 2146.987 65.90625
## Mato Grosso Mato Grosso 2189.062 52.87659
## Goiás Goiás 2134.678 59.28424
## Distrito Federal Distrito Federal 4060.471 237.89469É possível também calcular o intervalo de confiança para cada uma dessas estimativas:
confint(object=mediaRendaUF)## 2.5 % 97.5 %
## Rondônia 1677.451 1921.885
## Acre 1446.929 1783.197
## Amazonas 1552.423 2280.019
## Roraima 1841.009 2387.241
## Pará 1323.865 1507.370
## Amapá 1679.197 2286.978
## Tocantins 1643.441 1965.619
## Maranhão 1130.835 1403.979
## Piauí 1207.274 1470.475
## Ceará 1320.528 1556.482
## Rio Grande do Norte 1305.048 1680.054
## Paraíba 1424.057 1815.133
## Pernambuco 1474.604 1767.865
## Alagoas 1254.845 1447.154
## Sergipe 1335.430 1652.234
## Bahia 1305.401 1766.732
## Minas Gerais 1851.392 2010.346
## Espírito Santo 1913.963 2161.819
## Rio de Janeiro 2234.830 2444.443
## São Paulo 2600.025 3043.327
## Paraná 2246.673 2490.046
## Santa Catarina 2364.962 2515.789
## Rio Grande do Sul 2259.270 2479.395
## Mato Grosso do Sul 2017.813 2276.161
## Mato Grosso 2085.426 2292.698
## Goiás 2018.483 2250.872
## Distrito Federal 3594.206 4526.736É 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 desocupação por sexo e
cor ou raça utilizando as funções svyratio e
svyby e o respectivo cv:
txdesocupSexoRaca <- svyby(formula=~(VD4002=="Pessoas desocupadas"), denominator=~(VD4001=="Pessoas na força de trabalho"), by=~interaction(V2007,V2010), design=dadosPNADc, FUN=svyratio, vartype="cv", na.rm=TRUE)
txdesocupSexoRaca## interaction(V2007, V2010)
## 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
## VD4002 == "Pessoas desocupadas"/VD4001 == "Pessoas na força de trabalho"
## Homem.Branca 0.08447461
## Mulher.Branca 0.10873399
## Homem.Preta 0.12857753
## Mulher.Preta 0.16865089
## Homem.Amarela 0.05860442
## Mulher.Amarela 0.09522533
## Homem.Parda 0.11831124
## Mulher.Parda 0.16112811
## Homem.Indígena 0.08123925
## Mulher.Indígena 0.11519625
## Homem.Ignorado 0.00000000
## Mulher.Ignorado 0.10581381
## cv
## Homem.Branca 0.02020736
## Mulher.Branca 0.01981370
## Homem.Preta 0.02993692
## Mulher.Preta 0.03245521
## Homem.Amarela 0.22882259
## Mulher.Amarela 0.15705325
## Homem.Parda 0.01611358
## Mulher.Parda 0.01623548
## Homem.Indígena 0.20344305
## Mulher.Indígena 0.16448319
## Homem.Ignorado NaN
## Mulher.Ignorado 0.97762693Grá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 trabalhadas para cada um dos
casos:
svyhist(formula=~as.numeric(VD4035), design=dadosPNADc, main="Histograma", xlab="Número de Horas Trabalhadas")
svyhist(formula=~as.numeric(VD4035), design=dadosPNADc, freq=TRUE, main="Histograma", xlab="Número de Horas 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=VD4035~1, design=dadosPNADc, main="Boxplot do Número de Horas Trabalhadas")
svyboxplot(formula=VD4035~V2007, design=dadosPNADc, main="Boxplot do Número de Horas Trabalhadas por Sexo")
svyboxplot(formula=VD4035~V2007, design=dadosPNADc, all.outliers=TRUE, main="Boxplot do Número de Horas Trabalhadas por Sexo")
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=VD4020~VD4035, design=dadosPNADc, style="bubble", xlab="Número de Horas Efetivamente Trabalhadas", ylab="Rendimento Efetivo")
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=VD4020~VD4035, design=dadosPNADc, style="transparent", xlab="Número de Horas Efetivamente Trabalhadas", ylab="Rendimento Efetivo")
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 efetivos entre sexos:
svyttest(formula=VD4020~V2007, design=dadosPNADc)##
## Design-based t-test
##
## data: VD4020 ~ V2007
## t = -16.174, df = 198, p-value < 0.00000000000000022
## alternative hypothesis: true difference in mean is not equal to 0
## 95 percent confidence interval:
## -627.7154 -491.2814
## sample estimates:
## difference in mean
## -559.4984Outro exemplo para testar diferenças entre médias de horas trabalhadas, entre concluintes e não concluintes de graduação:
svyttest(formula=as.numeric(VD4035)~V3007, design=dadosPNADc)##
## Design-based t-test
##
## data: as.numeric(VD4035) ~ V3007
## t = -1.6556, df = 198, p-value = 0.09939
## alternative hypothesis: true difference in mean is not equal to 0
## 95 percent confidence interval:
## -2.0488732 0.1787247
## sample estimates:
## difference in mean
## -0.9350743Modelos 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
efetivo 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=VD4020~VD3004+V2010+V2009, design=dadosPNADc)
summary(object=modeloLin)##
## Call:
## svyglm(formula = VD4020 ~ VD3004 + V2010 + V2009, design = dadosPNADc)
##
## Survey design:
## svrepdesign.default(data = data_pnadc, weight = ~V1028, type = "bootstrap",
## repweights = "V1028[0-9]+", mse = TRUE, replicates = length(sprintf("V1028%03d",
## seq(1:200))), df = length(sprintf("V1028%03d", seq(1:200))))
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -673.879 78.707 -8.562
## VD3004Fundamental incompleto ou equivalente 491.895 29.271 16.805
## VD3004Fundamental completo ou equivalente 886.616 33.502 26.465
## VD3004Médio incompleto ou equivalente 1087.402 43.949 24.742
## VD3004Médio completo ou equivalente 1321.380 36.091 36.612
## VD3004Superior incompleto ou equivalente 1889.396 78.590 24.041
## VD3004Superior completo 4382.838 122.735 35.710
## V2010Preta -583.732 32.774 -17.811
## V2010Amarela 751.191 308.476 2.435
## V2010Parda -531.222 31.361 -16.939
## V2010Indígena -583.175 100.238 -5.818
## V2010Ignorado -931.322 301.047 -3.094
## V2009 38.398 1.793 21.410
## Pr(>|t|)
## (Intercept) 0.00000000000000402 ***
## 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.01582 *
## V2010Parda < 0.0000000000000002 ***
## V2010Indígena 0.00000002542545115 ***
## V2010Ignorado 0.00228 **
## V2009 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 13338647)
##
## Number of Fisher Scoring iterations: 2Além disso, é possível computar intervalos de confiança para os parâmetros:
confint(object=modeloLin)## 2.5 % 97.5 %
## (Intercept) -829.14612 -518.61224
## VD3004Fundamental incompleto ou equivalente 434.15146 549.63911
## VD3004Fundamental completo ou equivalente 820.52509 952.70596
## VD3004Médio incompleto ou equivalente 1000.70274 1174.10120
## VD3004Médio completo ou equivalente 1250.18155 1392.57915
## VD3004Superior incompleto ou equivalente 1734.35846 2044.43372
## VD3004Superior completo 4140.71544 4624.96032
## V2010Preta -648.38559 -519.07886
## V2010Amarela 142.65104 1359.73094
## V2010Parda -593.08959 -469.35516
## V2010Indígena -780.91796 -385.43288
## V2010Ignorado -1525.20630 -337.43695
## V2009 34.85989 41.93585Regressã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” ou
family=“quasibinomial”.
No exemplo abaixo modelamos a não-conclusão de graduação pelo sexo, cor ou raça e idade:
modeloLog <- svyglm(formula=V3007~V2007+V2010+V2009, design=dadosPNADc, family="quasibinomial")
summary(object=modeloLog)##
## Call:
## svyglm(formula = V3007 ~ V2007 + V2010 + V2009, design = dadosPNADc,
## family = "quasibinomial")
##
## Survey design:
## svrepdesign.default(data = data_pnadc, weight = ~V1028, type = "bootstrap",
## repweights = "V1028[0-9]+", mse = TRUE, replicates = length(sprintf("V1028%03d",
## seq(1:200))), df = length(sprintf("V1028%03d", seq(1:200))))
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.254639 0.115971 45.310 < 0.0000000000000002 ***
## V2007Mulher -0.021892 0.078058 -0.280 0.77943
## V2010Preta 0.463357 0.147931 3.132 0.00201 **
## V2010Amarela -0.675711 0.474186 -1.425 0.15578
## V2010Parda 0.355725 0.088203 4.033 0.0000794 ***
## V2010Indígena -0.095953 0.573004 -0.167 0.86719
## V2010Ignorado 9.852572 3.426597 2.875 0.00449 **
## V2009 -0.103320 0.003512 -29.423 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for quasibinomial family taken to be 0.9144395)
##
## Number of Fisher Scoring iterations: 11Sendo possível também computar intervalos de confiança para os parâmetros:
confint(object=modeloLog)## 2.5 % 97.5 %
## (Intercept) 5.0258983 5.48338063
## V2007Mulher -0.1758538 0.13206975
## V2010Preta 0.1715791 0.75513482
## V2010Amarela -1.6109948 0.25957217
## V2010Parda 0.1817534 0.52969606
## V2010Indígena -1.2261440 1.03423756
## V2010Ignorado 3.0939647 16.61117878
## V2009 -0.1102463 -0.09639394Aná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)
dadosPNADc <- convey_prep(design=dadosPNADc)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 efetiva:
giniBR <- svygini(formula=~VD4020, design=dadosPNADc, na.rm=TRUE)
giniBR## gini SE
## VD4020 0.51644 0.0053Assim como nas funções do pacote survey, é possível, por
exemplo, estimar o coeficiente de variação desta estimativa:
cv(object=giniBR)## VD4020
## 0.01025818Também é possível utilizar a função svyby para estimar o
índice de Gini da renda efetiva por Unidade da Federação:
giniUF <- svyby(formula=~VD4020, by=~UF, design=dadosPNADc, FUN=svygini, na.rm=TRUE)
giniUF## UF VD4020 se.VD4020
## Rondônia Rondônia 0.4387176 0.013307043
## Acre Acre 0.5351462 0.016819203
## Amazonas Amazonas 0.5808917 0.032323156
## Roraima Roraima 0.5266756 0.011719801
## Pará Pará 0.5265942 0.009844522
## Amapá Amapá 0.4968721 0.016590571
## Tocantins Tocantins 0.4799456 0.014632895
## Maranhão Maranhão 0.5438715 0.019909493
## Piauí Piauí 0.5532688 0.015755322
## Ceará Ceará 0.5458276 0.012258728
## Rio Grande do Norte Rio Grande do Norte 0.5134428 0.021239750
## Paraíba Paraíba 0.5655298 0.018328293
## Pernambuco Pernambuco 0.5128341 0.014738716
## Alagoas Alagoas 0.4409763 0.013890622
## Sergipe Sergipe 0.5454992 0.015362465
## Bahia Bahia 0.5810676 0.026843697
## Minas Gerais Minas Gerais 0.4941199 0.007122953
## Espírito Santo Espírito Santo 0.4774252 0.010215114
## Rio de Janeiro Rio de Janeiro 0.4659357 0.007544267
## São Paulo São Paulo 0.5129375 0.014295864
## Paraná Paraná 0.4681082 0.009527394
## Santa Catarina Santa Catarina 0.4064057 0.005987333
## Rio Grande do Sul Rio Grande do Sul 0.4800229 0.007047832
## Mato Grosso do Sul Mato Grosso do Sul 0.4688161 0.010049016
## Mato Grosso Mato Grosso 0.4351670 0.008644748
## Goiás Goiás 0.4673377 0.009561833
## Distrito Federal Distrito Federal 0.5701457 0.015118407E estimar intervalos de confiança para cada UF:
confint(object=giniUF)## 2.5 % 97.5 %
## Rondônia 0.4126363 0.4647989
## Acre 0.5021812 0.5681112
## Amazonas 0.5175395 0.6442439
## Roraima 0.5037053 0.5496460
## Pará 0.5072993 0.5458891
## Amapá 0.4643551 0.5293890
## Tocantins 0.4512656 0.5086255
## Maranhão 0.5048496 0.5828934
## Piauí 0.5223889 0.5841487
## Ceará 0.5218010 0.5698543
## Rio Grande do Norte 0.4718136 0.5550719
## Paraíba 0.5296070 0.6014526
## Pernambuco 0.4839467 0.5417214
## Alagoas 0.4137512 0.4682014
## Sergipe 0.5153893 0.5756091
## Bahia 0.5284549 0.6336802
## Minas Gerais 0.4801592 0.5080807
## Espírito Santo 0.4574040 0.4974465
## Rio de Janeiro 0.4511492 0.4807222
## São Paulo 0.4849181 0.5409568
## Paraná 0.4494348 0.4867815
## Santa Catarina 0.3946708 0.4181407
## Rio Grande do Sul 0.4662094 0.4938364
## Mato Grosso do Sul 0.4491204 0.4885118
## Mato Grosso 0.4182236 0.4521104
## Goiás 0.4485968 0.4860785
## Distrito Federal 0.5405142 0.5997772Curva 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 efetiva:
curvaLorenz <- svylorenz(formula=~VD4020, design=dadosPNADc, quantiles=seq(from=0, to=1, by=0.05), na.rm=TRUE)
Instituto Brasileiro de Geografia e Estatística - pacotesipd@ibge.gov.br↩︎