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.

Além dos argumentos 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. 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 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.  138. 512631 6.34e6
##  2 2017  4         Rondônia 11000… 1110011 01    06     121.  138. 512631 8.65e6
##  3 2017  4         Rondônia 11000… 1110011 03    06     121.  148. 512631 7.52e6
##  4 2017  4         Rondônia 11000… 1110011 03    06     121.  148. 512631 8.52e6
##  5 2017  4         Rondônia 11000… 1110011 03    06     121.  148. 512631 7.43e6
##  6 2017  4         Rondônia 11000… 1110011 04    06     121.  177. 512631 8.65e6
##  7 2017  4         Rondônia 11000… 1110011 04    06     121.  177. 512631 8.53e6
##  8 2017  4         Rondônia 11000… 1110011 05    06     121.  132. 512631 8.04e6
##  9 2017  4         Rondônia 11000… 1110011 05    06     121.  132. 512631 6.34e6
## 10 2017  4         Rondônia 11000… 1110011 05    06     121.  132. 512631 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.  138. 512631 6.34e6
##  2 2017  4         11    110000016 1110011 01    06     121.  138. 512631 8.65e6
##  3 2017  4         11    110000016 1110011 03    06     121.  148. 512631 7.52e6
##  4 2017  4         11    110000016 1110011 03    06     121.  148. 512631 8.52e6
##  5 2017  4         11    110000016 1110011 03    06     121.  148. 512631 7.43e6
##  6 2017  4         11    110000016 1110011 04    06     121.  177. 512631 8.65e6
##  7 2017  4         11    110000016 1110011 04    06     121.  177. 512631 8.53e6
##  8 2017  4         11    110000016 1110011 05    06     121.  132. 512631 8.04e6
##  9 2017  4         11    110000016 1110011 05    06     121.  132. 512631 6.34e6
## 10 2017  4         11    110000016 1110011 05    06     121.  132. 512631 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.

Neste exemplo, mostraremos como importar os microdados anuais da 1ª visita de 2017, que compreende os temas:
  • 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.

Neste exemplo, mostraremos como importar os microdados anuais do 2º trimestre de 2017, que compreende os temas:
  • 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ção get_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 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_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 é 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 mensal efetiva:

totalrenda <- svytotal(x=~VD4020, design=dadosPNADc, na.rm=TRUE)
totalrenda
##               total         SE
## VD4020 196383378821 2904924488

Alé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.01479211
confint(object=totalrenda)
##               2.5 %       97.5 %
## VD4020 190689831446 202076926195
confint(object=totalrenda, level=0.99)
##               0.5 %       99.5 %
## VD4020 188900789200 203865968442

Variá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  101168151 0.0008
## V2007Mulher 105660589 0.0008

E 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    101168151      0.0008
## V2007Mulher   105660589      0.0008
## V2010Branca    88962855 405609.7265
## V2010Preta     18042962 200447.8939
## V2010Amarela    1177942  63965.7389
## V2010Parda     98070916 390227.4130
## V2010Indígena    553711  30014.0190
## V2010Ignorado     20354   5623.0329

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(V2007,V2010), design=dadosPNADc, na.rm=TRUE)
ftable(x=totalsexoEraca)
##                                                     A            B
## interaction(V2007, V2010)Homem.Branca    42785223.917   219669.092
## interaction(V2007, V2010)Mulher.Branca   46177630.546   226174.723
## interaction(V2007, V2010)Homem.Preta      9137903.722   108004.561
## interaction(V2007, V2010)Mulher.Preta     8905058.613   117212.913
## interaction(V2007, V2010)Homem.Amarela     532742.237    33412.608
## interaction(V2007, V2010)Mulher.Amarela    645199.753    37185.801
## interaction(V2007, V2010)Homem.Parda     48449054.043   212701.133
## interaction(V2007, V2010)Mulher.Parda    49621862.115   214874.798
## interaction(V2007, V2010)Homem.Indígena    252469.831    16168.029
## interaction(V2007, V2010)Mulher.Indígena   301240.820    18624.701
## interaction(V2007, V2010)Homem.Ignorado     10757.250     2905.930
## interaction(V2007, V2010)Mulher.Ignorado     9597.154     3581.625

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 efetiva pode ser facilmente reescrito para médias:

mediarenda <- svymean(x=~VD4020, design=dadosPNADc, na.rm=TRUE)
mediarenda
##          mean     SE
## VD4020 2182.6 31.166

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

cv(object=mediarenda)
##     VD4020 
## 0.01427965
confint(object=mediarenda)
##           2.5 %   97.5 %
## VD4020 2121.499 2243.669

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=~V2007, design=dadosPNADc, na.rm=TRUE)
propsexo
##                mean SE
## V2007Homem  0.48914  0
## V2007Mulher 0.51086  0

Da 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.489139715 0.0000
## V2007Mulher   0.510860285 0.0000
## V2010Branca   0.430128107 0.0020
## V2010Preta    0.087236244 0.0010
## V2010Amarela  0.005695253 0.0003
## V2010Parda    0.474164839 0.0019
## V2010Indígena 0.002677146 0.0001
## V2010Ignorado 0.000098412 0.0000

E 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.20686304968 0.00106208205
## interaction(V2007, V2010)Mulher.Branca   0.22326505758 0.00109353624
## interaction(V2007, V2010)Homem.Preta     0.04418101528 0.00052219320
## interaction(V2007, V2010)Mulher.Preta    0.04305522827 0.00056671482
## interaction(V2007, V2010)Homem.Amarela   0.00257576504 0.00016154722
## interaction(V2007, V2010)Mulher.Amarela  0.00311948791 0.00017979030
## interaction(V2007, V2010)Homem.Parda     0.23424720396 0.00102839254
## interaction(V2007, V2010)Mulher.Parda    0.23991763483 0.00103890203
## interaction(V2007, V2010)Homem.Indígena  0.00122067093 0.00007817110
## interaction(V2007, V2010)Mulher.Indígena 0.00145647466 0.00009004890
## interaction(V2007, V2010)Homem.Ignorado  0.00005201042 0.00001404993
## interaction(V2007, V2010)Mulher.Ignorado 0.00004640145 0.00001731686

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 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.1189622
## SEs=
##             [,1]
## [1,] 0.001161197

Cá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.009761057
confint(object=txdesocup)
##                                                                              2.5 %
## VD4002 == "Pessoas desocupadas"/VD4001 == "Pessoas na força de trabalho" 0.1166863
##                                                                             97.5 %
## VD4002 == "Pessoas desocupadas"/VD4001 == "Pessoas na força de trabalho" 0.1212381

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 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.97.5
## 0.5     1300   1300    1400 25.35553
## 
## attr(,"hasci")
## [1] TRUE
## attr(,"class")
## [1] "newsvyquantile"
SE(object=medianarenda)
##   VD4020 
## 25.35553
cv(object=medianarenda)
##     VD4020 
## 0.01950425

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=~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 268854596417 4007968380

E 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 2988 43.014

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.

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 1860.4 22.572

Ou 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.08557168
## SEs=
##             [,1]
## [1,] 0.001012151

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=~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.105920 0.0015
## VD3004Fundamental incompleto ou equivalente    0.423641 0.0029
## VD3004Fundamental completo ou equivalente      0.088819 0.0016
## VD3004Médio incompleto ou equivalente          0.043352 0.0011
## VD3004Médio completo ou equivalente            0.240204 0.0025
## VD3004Superior incompleto ou equivalente       0.022557 0.0010
## VD3004Superior completo                        0.075507 0.0016

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.

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

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=~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.4980244   0.5019756 0.0026381409
## Fundamental incompleto ou equivalente     0.5124954   0.4875046 0.0009569814
## Fundamental completo ou equivalente       0.5043056   0.4956944 0.0029622978
## Médio incompleto ou equivalente           0.5043189   0.4956811 0.0033950945
## Médio completo ou equivalente             0.4709257   0.5290743 0.0014756623
## Superior incompleto ou equivalente        0.4705091   0.5294909 0.0046038732
## Superior completo                         0.4177301   0.5822699 0.0024185658
##                                                   se2
## Sem instrução e menos de 1 ano de estudo 0.0026381409
## Fundamental incompleto ou equivalente    0.0009569814
## Fundamental completo ou equivalente      0.0029622978
## Médio incompleto ou equivalente          0.0033950945
## Médio completo ou equivalente            0.0014756623
## Superior incompleto ou equivalente       0.0046038732
## Superior completo                        0.0024185658

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=~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.08269097
## Mulher Mulher                                     0.07925675
##        VD3004Fundamental incompleto ou equivalente
## Homem                                    0.3863812
## Mulher                                   0.3495031
##        VD3004Fundamental completo ou equivalente
## Homem                                 0.08567548
## Mulher                                0.08007979
##        VD3004Médio incompleto ou equivalente
## Homem                             0.07245218
## Mulher                            0.06771656
##        VD3004Médio completo ou equivalente
## Homem                            0.2324130
## Mulher                           0.2482966
##        VD3004Superior incompleto ou equivalente VD3004Superior completo
## Homem                                0.04281874              0.09756842
## Mulher                               0.04582164              0.12932549
##                 se1         se2          se3          se4         se5
## Homem  0.0007134664 0.001595710 0.0007758974 0.0007343163 0.001340821
## Mulher 0.0006855442 0.001348694 0.0007215427 0.0006784040 0.001188696
##                 se6         se7
## Homem  0.0006334596 0.001565445
## Mulher 0.0006847164 0.001542099

També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 1802.134  62.19958
## Acre                               Acre 1623.126  86.92702
## Amazonas                       Amazonas 1910.835 182.58939
## Roraima                         Roraima 2099.078 136.35270
## Pará                               Pará 1421.336  47.42050
## Amapá                             Amapá 1973.229 153.13325
## Tocantins                     Tocantins 1799.273  81.39674
## Maranhão                       Maranhão 1264.903  69.47827
## Piauí                             Piauí 1336.593  67.00924
## Ceará                             Ceará 1445.663  61.05556
## Rio Grande do Norte Rio Grande do Norte 1497.004  98.14735
## Paraíba                         Paraíba 1607.444  98.03307
## Pernambuco                   Pernambuco 1623.816  75.33920
## Alagoas                         Alagoas 1347.245  48.88096
## Sergipe                         Sergipe 1503.512  82.42209
## Bahia                             Bahia 1564.988 127.87174
## Minas Gerais               Minas Gerais 1932.488  40.58339
## Espírito Santo           Espírito Santo 2051.034  63.49199
## Rio de Janeiro           Rio de Janeiro 2336.392  53.21622
## São Paulo                     São Paulo 2823.983 112.90861
## Paraná                           Paraná 2379.388  62.65090
## Santa Catarina           Santa Catarina 2437.859  38.22309
## Rio Grande do Sul     Rio Grande do Sul 2376.719  56.88476
## Mato Grosso do Sul   Mato Grosso do Sul 2145.462  65.52613
## Mato Grosso                 Mato Grosso 2186.392  52.71229
## Goiás                             Goiás 2140.620  59.71019
## Distrito Federal       Distrito Federal 4048.095 235.21853

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

confint(object=mediaRendaUF)
##                        2.5 %   97.5 %
## Rondônia            1680.225 1924.043
## Acre                1452.752 1793.500
## Amazonas            1552.966 2268.703
## Roraima             1831.832 2366.325
## Pará                1328.393 1514.278
## Amapá               1673.094 2273.365
## Tocantins           1639.739 1958.808
## Maranhão            1128.729 1401.078
## Piauí               1205.258 1467.929
## Ceará               1325.997 1565.330
## Rio Grande do Norte 1304.639 1689.369
## Paraíba             1415.303 1799.585
## Pernambuco          1476.154 1771.478
## Alagoas             1251.440 1443.049
## Sergipe             1341.967 1665.056
## Bahia               1314.364 1815.612
## Minas Gerais        1852.946 2012.030
## Espírito Santo      1926.592 2175.476
## Rio de Janeiro      2232.091 2440.694
## São Paulo           2602.686 3045.280
## Paraná              2256.594 2502.181
## Santa Catarina      2362.943 2512.775
## Rio Grande do Sul   2265.227 2488.211
## Mato Grosso do Sul  2017.034 2273.891
## Mato Grosso         2083.078 2289.706
## Goiás               2023.590 2257.650
## Distrito Federal    3587.075 4509.115

É 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.08475640
## Mulher.Branca                                                                 0.10895098
## Homem.Preta                                                                   0.12863903
## Mulher.Preta                                                                  0.16877014
## Homem.Amarela                                                                 0.05842629
## Mulher.Amarela                                                                0.09516543
## Homem.Parda                                                                   0.11841542
## Mulher.Parda                                                                  0.16108434
## Homem.Indígena                                                                0.08126455
## Mulher.Indígena                                                               0.11657355
## Homem.Ignorado                                                                0.00000000
## Mulher.Ignorado                                                               0.10761555
##                         cv
## Homem.Branca    0.02026444
## Mulher.Branca   0.01993706
## Homem.Preta     0.02990095
## Mulher.Preta    0.03245517
## Homem.Amarela   0.22660061
## Mulher.Amarela  0.15784742
## Homem.Parda     0.01602689
## Mulher.Parda    0.01617196
## Homem.Indígena  0.20303502
## Mulher.Indígena 0.16543479
## Homem.Ignorado         NaN
## Mulher.Ignorado 0.97330024

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 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.005, df = 198, p-value < 0.00000000000000022
## alternative hypothesis: true difference in mean is not equal to 0
## 95 percent confidence interval:
##  -626.0983 -488.7370
## sample estimates:
## difference in mean 
##          -557.4177

Outro 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.6373, df = 198, p-value = 0.1032
## alternative hypothesis: true difference in mean is not equal to 0
## 95 percent confidence interval:
##  -2.0393618  0.1891236
## sample estimates:
## difference in mean 
##         -0.9251191

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 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)                                  -678.62      79.22  -8.566
## VD3004Fundamental incompleto ou equivalente   491.69      29.66  16.577
## VD3004Fundamental completo ou equivalente     887.00      33.87  26.189
## VD3004Médio incompleto ou equivalente        1089.30      44.28  24.597
## VD3004Médio completo ou equivalente          1321.88      36.42  36.298
## VD3004Superior incompleto ou equivalente     1898.70      84.23  22.541
## VD3004Superior completo                      4387.09     122.88  35.703
## V2010Preta                                   -583.89      32.80 -17.799
## V2010Amarela                                  758.64     309.53   2.451
## V2010Parda                                   -529.65      31.75 -16.684
## V2010Indígena                                -581.63      99.96  -5.819
## V2010Ignorado                                -933.58     294.71  -3.168
## V2009                                          38.55       1.80  21.420
##                                                         Pr(>|t|)    
## (Intercept)                                  0.00000000000000391 ***
## 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.01517 *  
## V2010Parda                                  < 0.0000000000000002 ***
## V2010Indígena                                0.00000002531158936 ***
## V2010Ignorado                                            0.00179 ** 
## V2009                                       < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 3061666267046)
## 
## 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)                                  -834.90403 -522.34267
## VD3004Fundamental incompleto ou equivalente   433.17851  550.20704
## VD3004Fundamental completo ou equivalente     820.18543  953.81745
## VD3004Médio incompleto ou equivalente        1001.93326 1176.65835
## VD3004Médio completo ou equivalente          1250.04065 1393.72613
## VD3004Superior incompleto ou equivalente     1732.52700 2064.86746
## VD3004Superior completo                      4144.68707 4629.49792
## V2010Preta                                   -648.60781 -519.17937
## V2010Amarela                                  148.01140 1369.26180
## V2010Parda                                   -592.27935 -467.02867
## V2010Indígena                                -778.82455 -384.44452
## V2010Ignorado                               -1514.96806 -352.18801
## V2009                                          34.99628   42.09649

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 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="binomial")
## Warning in eval(family$initialize): non-integer #successes in a binomial glm!

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

## Warning in eval(family$initialize): non-integer #successes in a binomial glm!
summary(object=modeloLog)
## 
## Call:
## svyglm(formula = V3007 ~ V2007 + V2010 + V2009, design = dadosPNADc, 
##     family = "binomial")
## 
## 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.252620   0.116202  45.203 < 0.0000000000000002 ***
## V2007Mulher   -0.019232   0.078109  -0.246              0.80577    
## V2010Preta     0.458235   0.149045   3.074              0.00242 ** 
## V2010Amarela  -0.695098   0.481922  -1.442              0.15083    
## V2010Parda     0.352718   0.088074   4.005            0.0000887 ***
## V2010Indígena -0.086270   0.573350  -0.150              0.88056    
## V2010Ignorado  9.866664   3.442063   2.866              0.00461 ** 
## V2009         -0.103383   0.003497 -29.563 < 0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 18164.51)
## 
## Number of Fisher Scoring iterations: 11

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

confint(object=modeloLog)
##                    2.5 %      97.5 %
## (Intercept)    5.0234235  5.48181558
## V2007Mulher   -0.1732937  0.13482959
## V2010Preta     0.1642597  0.75221032
## V2010Amarela  -1.6456381  0.25544281
## V2010Parda     0.1790017  0.52643376
## V2010Indígena -1.2171432  1.04460400
## V2010Ignorado  3.0775504 16.65577691
## V2009         -0.1102810 -0.09648594

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)
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.51662 0.0053

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

cv(object=giniBR)
##     VD4020 
## 0.01032032

També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.4383580 0.013274791
## Acre                               Acre 0.5337446 0.016942518
## Amazonas                       Amazonas 0.5797876 0.031786151
## Roraima                         Roraima 0.5265372 0.011727658
## Pará                               Pará 0.5263289 0.009915294
## Amapá                             Amapá 0.4958021 0.016589781
## Tocantins                     Tocantins 0.4790207 0.014615431
## Maranhão                       Maranhão 0.5450098 0.019778752
## Piauí                             Piauí 0.5530826 0.015760267
## Ceará                             Ceará 0.5456317 0.012357868
## Rio Grande do Norte Rio Grande do Norte 0.5145410 0.021915709
## Paraíba                         Paraíba 0.5639860 0.018218409
## Pernambuco                   Pernambuco 0.5127469 0.014798173
## Alagoas                         Alagoas 0.4414840 0.013905890
## Sergipe                         Sergipe 0.5466702 0.015427913
## Bahia                             Bahia 0.5838889 0.028499272
## Minas Gerais               Minas Gerais 0.4936558 0.007116192
## Espírito Santo           Espírito Santo 0.4787428 0.010079671
## Rio de Janeiro           Rio de Janeiro 0.4651384 0.007519689
## São Paulo                     São Paulo 0.5127561 0.014253134
## Paraná                           Paraná 0.4685151 0.009574032
## Santa Catarina           Santa Catarina 0.4056347 0.005938358
## Rio Grande do Sul     Rio Grande do Sul 0.4796909 0.007126207
## Mato Grosso do Sul   Mato Grosso do Sul 0.4681657 0.010003687
## Mato Grosso                 Mato Grosso 0.4344578 0.008621995
## Goiás                             Goiás 0.4671751 0.009601496
## Distrito Federal       Distrito Federal 0.5693400 0.014915770

E estimar intervalos de confiança para cada UF:

confint(object=giniUF)
##                         2.5 %    97.5 %
## Rondônia            0.4123399 0.4643761
## Acre                0.5005379 0.5669513
## Amazonas            0.5174879 0.6420873
## Roraima             0.5035514 0.5495230
## Pará                0.5068953 0.5457626
## Amapá               0.4632867 0.5283175
## Tocantins           0.4503750 0.5076664
## Maranhão            0.5062441 0.5837754
## Piauí               0.5221931 0.5839722
## Ceará               0.5214107 0.5698526
## Rio Grande do Norte 0.4715870 0.5574950
## Paraíba             0.5282785 0.5996934
## Pernambuco          0.4837430 0.5417508
## Alagoas             0.4142289 0.4687390
## Sergipe             0.5164320 0.5769084
## Bahia               0.5280314 0.6397465
## Minas Gerais        0.4797083 0.5076033
## Espírito Santo      0.4589870 0.4984986
## Rio de Janeiro      0.4504001 0.4798767
## São Paulo           0.4848205 0.5406918
## Paraná              0.4497503 0.4872798
## Santa Catarina      0.3939958 0.4172737
## Rio Grande do Sul   0.4657238 0.4936581
## Mato Grosso do Sul  0.4485588 0.4877725
## Mato Grosso         0.4175590 0.4513566
## Goiás               0.4483566 0.4859937
## Distrito Federal    0.5401056 0.5985744

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 efetiva:

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


  1. Instituto Brasileiro de Geografia e Estatística - ↩︎