PNADcIBGE e
surveyO 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: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.
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)
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.
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")
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. 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.
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.
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.
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.
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.
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.
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).
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)
svytotal. Sua sintaxe precisa de três parâmetros
principais:
~;
na.rm=TRUE, que remove as observações onde a
variável é não-aplicável.
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
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.
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
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.
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
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"
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.
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
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.
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
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.
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")
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")
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")
surveyNo pacote survey ainda há a possibilidade de realizar
testes de hipóteses e estimar parâmetros de modelos, considerando o
plano amostral complexo.
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
Além de testes de hipóteses, é possível estimar modelos lineares
generalizados para dados amostrais complexos através da função
svyglm.
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
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
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.
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
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↩︎