Análise de microdados da PNS
Com os pacotes PNSIBGE e
survey
Gabriel Assuncao1
01 de fevereiro de 2024
O pacote PNSIBGE
O pacote PNSIBGE foi desenvolvido para facilitar o download, importação e análise dos dados amostrais da Pesquisa Nacional de Saúde realizada pelo Instituto Brasileiro de Geografia e Estatística - IBGE.
A PNS possui três tipos de microdados:- Questionário básico, que contém a parte básica investigada pela pesquisa, obtendo informações de todos os moradores dos domicílios selecionados;
- Questionário de morador selecionado, que contém a parte aplicada somente para um morador selecionado;
- Questionário de antropometria, que contém a parte com as medidas antropométricas, podendo ser aplicado juntamente com o questionário de morador selecionado ou utilizando uma subamostra específica.
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/PNSIBGE.
É
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("PNSIBGE")Antes de utilizar o pacote é necessário carregá-lo no R através do comando:
library(PNSIBGE)Importação de microdados
Para a importação dos microdados da PNS, 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_pns 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 de todos os níveis investigados pela pesquisa.
help("get_pns")Microdados do questionário básico
A importação dos microdados do questionário básico através da função
get_pns é muito simples. Basta indicar o ano dos microdados
desejados no respectivo argumento da função. Abaixo um exemplo de como
importar os microdados do questionário básico da PNS 2019:
dadosPNS <- get_pns(year=2019)Apenas com esse comando, os microdados da pesquisa já estão
disponíveis no objeto dadosPNS para análise.
year, que indica o ano dos microdados a
serem baixados, a função get_pns ainda possui outros cinco
argumentos que podem ser modificados no download para este tipo de
microdados:
-
vars: Este argumento recebe um vetor de caracteres com o nome das variáveis a serem baixadas. Caso nenhuma variável seja passada, todas as variáveis disponíveis na pesquisa são baixadas. É útil caso deseje trabalhar com poucas variáveis, pois assim o objeto ocupará um espaço menor na memória do computador; -
labels: Um argumento lógico que indica se os níveis das variáveis categóricas devam ser rotuladas de acordo com o dicionário da pesquisa. Odefaulté rotulá-los; -
deflator: Um argumento lógico que indica se as variáveis para deflacionamento serão incluídas nos microdados. Odefaulté incluí-las; -
design: Um argumento lógico que indica se a função deve retornar um objeto do plano amostral para análise com o pacote survey. Casodesign=FALSE, a função retorna apenas um data-frame com os microdados. É altamente recomendado que mantenha essa opção comoTRUE, caso contrário suas análises poderão ser feitas de forma incorreta; -
reload: Um argumento lógico que indica se a função deve baixar novamente os arquivos no diretório de referência mesmo que estes arquivos já tenham sido baixados; -
curlopts: Uma lista com indicação das opções para utilização nas funções do pacote RCurl; -
savedir: O endereço onde devem ser salvos os arquivos baixados. Padrão é utilizar uma pasta temporária.
Por exemplo, um usuário que deseje trabalhar apenas com as variáveis
de Diagnóstico de alguma doença crônica, física ou mental, ou
doença de longa duração (J007) e Informação sobre a
procura de um mesmo médico ou serviço de saúde para atendimento
(J009) da PNS 2019, pode utilizar o argumento vars para
selecionar apenas estas variáveis:
dadosPNS <- get_pns(year=2019, vars=c("J007","J009"))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.
dadosPNS## Stratified 1 - level Cluster Sampling design (with replacement)
## With (8028) clusters.
## survey::postStratify(design = data_prior, strata = ~V00283, population = popc.types)class(dadosPNS)## [1] "survey.design2" "survey.design"Caso o usuário não queira trabalhar com este objeto, ele pode
escolher a opção design=FALSE para baixar os microdados
brutos:
dadosPNS_brutos <- get_pns(year=2019, vars=c("J007","J009"), design=FALSE)O Objeto resultante é um data-frame com as variáveis selecionadas e as variáveis relacionadas ao processo de amostragem:
dadosPNS_brutos## # A tibble: 279,382 × 17
## V0020 V0001 V0024 UPA_PNS V0006_PNS V0015 C00301 J007 J009 M001 W001
## <chr> <fct> <chr> <chr> <chr> <fct> <chr> <fct> <fct> <fct> <fct>
## 1 2019 Rondônia 1110011 110000… 0001 Real… 01 Sim Sim Real… <NA>
## 2 2019 Rondônia 1110011 110000… 0001 Real… 02 Sim Sim <NA> <NA>
## 3 2019 Rondônia 1110011 110000… 0001 Real… 03 Não Sim <NA> <NA>
## 4 2019 Rondônia 1110011 110000… 0001 Real… 04 Não Sim <NA> <NA>
## 5 2019 Rondônia 1110011 110000… 0001 Real… 05 Não Sim <NA> <NA>
## 6 2019 Rondônia 1110011 110000… 0001 Real… 06 Não Sim <NA> <NA>
## 7 2019 Rondônia 1110011 110000… 0002 Real… 01 Não Não <NA> <NA>
## 8 2019 Rondônia 1110011 110000… 0002 Real… 02 Não Não <NA> <NA>
## 9 2019 Rondônia 1110011 110000… 0002 Real… 03 Não Não <NA> <NA>
## 10 2019 Rondônia 1110011 110000… 0002 Real… 04 Não Não Real… <NA>
## # ℹ 279,372 more rows
## # ℹ 6 more variables: V0028 <dbl>, V00281 <dbl>, V00282 <dbl>, V00283 <chr>,
## # ID_DOMICILIO <chr>, Deflator <dbl>Também é possível baixar os dados sem incluir os rótulos dos níveis,
através do argumento labels.
dadosPNS_brutos_sem <- get_pns(year=2019, vars=c("J007","J009"), labels=FALSE, design=FALSE)Perceba que agora os níveis das categorias são representados por números:
dadosPNS_brutos_sem## # A tibble: 279,382 × 17
## V0020 V0001 V0024 UPA_PNS V0006_PNS V0015 C00301 J007 J009 M001 W001
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 2019 11 1110011 110000016 0001 01 01 1 1 1 <NA>
## 2 2019 11 1110011 110000016 0001 01 02 1 1 <NA> <NA>
## 3 2019 11 1110011 110000016 0001 01 03 2 1 <NA> <NA>
## 4 2019 11 1110011 110000016 0001 01 04 2 1 <NA> <NA>
## 5 2019 11 1110011 110000016 0001 01 05 2 1 <NA> <NA>
## 6 2019 11 1110011 110000016 0001 01 06 2 1 <NA> <NA>
## 7 2019 11 1110011 110000016 0002 01 01 2 2 <NA> <NA>
## 8 2019 11 1110011 110000016 0002 01 02 2 2 <NA> <NA>
## 9 2019 11 1110011 110000016 0002 01 03 2 2 <NA> <NA>
## 10 2019 11 1110011 110000016 0002 01 04 2 2 1 <NA>
## # ℹ 279,372 more rows
## # ℹ 6 more variables: V0028 <dbl>, V00281 <dbl>, V00282 <dbl>, V00283 <chr>,
## # ID_DOMICILIO <chr>, Deflator <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 do questionário de morador selecionado
Além dos microdados do questionário básico, a função
get_pns também permite a importação dos microdados do
questionário de morador selecionado.
Neste exemplo, mostraremos como importar os microdados do questionário de morador selecionado da PNS 2019.
Para a importação dos microdados do questionário de morador
selecionado basta utilizar o valor TRUE no argumento
selected:
dadosPNS_morador_selecionado <- get_pns(year=2019, selected=TRUE)
dadosPNS_morador_selecionado## Stratified 1 - level Cluster Sampling design (with replacement)
## With (8027) clusters.
## survey::postStratify(design = data_prior, strata = ~V00293, population = popc.types)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 do questionário básico.
Microdados do questionário de antropometria
Conforme mencionado, além dos microdados do questionário básico, a
função get_pns também permite a importação dos microdados
do questionário de antropometria, quando tal módulo for uma subamostra
não diretamente relacionada ao questionário de morador selecionado.
Neste exemplo, mostraremos como importar os microdados do questionário de antropometria da PNS 2019.
Para a importação dos microdados do questionário de antropometria
basta utilizar o valor TRUE no argumento
anthropometry:
dadosPNS_antropometria <- get_pns(year=2019, anthropometry=TRUE)
dadosPNS_antropometria## Stratified 1 - level Cluster Sampling design (with replacement)
## With (1150) clusters.
## survey::postStratify(design = data_prior, strata = ~V00303, population = popc.types)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 do questionário básico.
Importação offline
Caso não seja possível utilizar a importação online através da funçãoget_pns, é possível criar o mesmo objeto para
análise com arquivos que estejam em disco. Para isto, são utilizadas
quatro funções:
-
read_pns: Para a leitura do arquivo .txt dos microdados; -
pns_labeller: Opcional. Coloca os rótulos dos níveis nas variáveis categóricas; -
pns_deflator: Opcional. Insere as variáveis para deflacionamento nos microdados; -
pns_design: Cria o objeto do plano amostral para a análise com o pacotesurvey.
Para utilizar estas funções, primeiramente é necessário ter os
microdados e sua documentação no disco e extraídos dos arquivos
compactados, se for o caso.
Estes arquivos podem ser baixados
diretamente do FTP do
IBGE.
Para a função read_pns são utilizados os arquivos de
texto contendo os microdados e o input para SAS.
Exemplo de leitura para a PNS 2019:
dadosPNS <- read_pns(microdata="PNS_2019.txt", input_txt="input_PNS_2019.txt")Assim como na importação online, pode ser utilizado o
argumento vars para definir quais variáveis serão
lidas.
A função pns_labeller utiliza o objeto criado pela
função read_pns e o arquivo de dicionário das variáveis
presente na documentação da pesquisa.
dadosPNS <- pns_labeller(data_pns=dadosPNS, dictionary.file="dicionario_PNS_microdados_2019.xls")A utilização de pns_labeller é opcional e não interfere
nos resultados obtidos.
A função pns_deflator utiliza o objeto criado pela
função read_pns e o arquivo de deflatores presente na
documentação da pesquisa.
dadosPNS <- pns_deflator(data_pns=dadosPNS, deflator.file="deflator_PNS.xls")A utilização de pns_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 pns_design aplicada no objeto que
contém os microdados:
dadosPNS <- pns_design(data_pns=dadosPNS)O objeto resultante é o mesmo retornado pela função
get_pns e pode ser utilizado para análises com o pacote
survey.
A importação offline é feita da mesma forma tanto para microdados do questionário básico, quanto para os microdados dos demais tipos de questionários abrangidos pela pesquisa (morador selecionado e antropometria). No entanto, para os demais questionários é preciso utilizar as variáveis identificadoras de seleção da subamostra e também remover as variáveis de peso relativas aos outros tipos de questionário.
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 PNS de 2019, com algumas variáveis selecionadas de acordo com a tabela abaixo:
| Variável | Descrição |
|---|---|
| V0001 | Unidade da Federação |
| C001 | Número de pessoas no domicílio |
| C004 | Condição no domicílio |
| C006 | Sexo |
| C008 | Idade do morador na data de referência |
| C009 | Cor ou raça |
| E01602 | Rendimento bruto mensal normalmente recebido para pessoas de 14 anos ou mais de idade |
| E017 | Número de horas normalmente trabalhadas para pessoas de 14 anos ou mais de idade |
| J002 | Interrupção na realização de quaisquer atividades habituais por motivo da própria saúde |
| J00402 | Motivo de saúde que impediu a realização das atividades habituais |
| J007 | Diagnóstico de alguma doença crônica, física ou mental, ou doença de longa duração |
| VDD004A | Nível de instrução mais elevado alcançado |
Importando os dados:
variaveis_selecionadas <- c("V0001","C001","C004","C006","C008","C009","E01602","E017","J002","J00402","J007","VDD004A")
dadosPNS <- get_pns(year=2019, vars=variaveis_selecionadas)O objeto dadosPNS vai ser utilizado com o pacote
survey para a análise dos microdados.
Para iniciar a análise dos dados é necessário carregar o pacote
survey.
library(survey)Estimando totais
A função do pacote para a estimação de totais populacionais é asvytotal. Sua sintaxe precisa de três parâmetros
principais:
-
O nome da variável que se deseja calcular o total, precedido por um
~; - O nome do objeto do plano amostral;
-
A opção
na.rm=TRUE, que remove as observações onde a variável é não-aplicável.
Variáveis numéricas
Ela pode ser utilizada para a estimação do total de uma variável numérica, como a renda bruta mensal normalmente recebida:
totalrenda <- svytotal(x=~E01602, design=dadosPNS, na.rm=TRUE)
totalrenda## total SE
## E01602 213227874692 3604489769Além da estimativa do total da renda bruta mensal normalmente
recebida, 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)## E01602
## E01602 0.0169044confint(object=totalrenda)## 2.5 % 97.5 %
## E01602 206163204561 220292544822confint(object=totalrenda, level=0.99)## 0.5 % 99.5 %
## E01602 203943324320 222512425064Variáveis categóricas
Também é possível estimar totais populacionais de categorias, utilizando variáveis categóricas, como o sexo:
totalsexo <- svytotal(x=~C006, design=dadosPNS, na.rm=TRUE)
totalsexo## total SE
## C006Homem 100158109 240473
## C006Mulher 109431498 240473E estimar o total de mais de uma variável categórica no mesmo código,
separando com o operador +:
totalsexoraca <- svytotal(x=~C006+C009, design=dadosPNS, na.rm=TRUE)
totalsexoraca## total SE
## C006Homem 100158109 240473.3
## C006Mulher 109431498 240473.3
## C009Branca 91037722 686604.4
## C009Preta 21786515 338021.3
## C009Amarela 1674691 144663.7
## C009Parda 94086142 600647.8
## C009Indígena 990265 74430.0
## C009Ignorado 14272 5707.5També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(C006,C009), design=dadosPNS, na.rm=TRUE)
ftable(x=totalsexoEraca)## A B
## interaction(C006, C009)Homem.Branca 42682904.814 378164.778
## interaction(C006, C009)Mulher.Branca 48354816.819 404879.054
## interaction(C006, C009)Homem.Preta 10691163.972 186451.899
## interaction(C006, C009)Mulher.Preta 11095350.705 200050.616
## interaction(C006, C009)Homem.Amarela 754080.767 73347.569
## interaction(C006, C009)Mulher.Amarela 920610.596 81402.906
## interaction(C006, C009)Homem.Parda 45496265.464 337575.630
## interaction(C006, C009)Mulher.Parda 48589876.693 351852.002
## interaction(C006, C009)Homem.Indígena 523966.395 52076.796
## interaction(C006, C009)Mulher.Indígena 466298.332 33901.415
## interaction(C006, C009)Homem.Ignorado 9727.597 3936.467
## interaction(C006, C009)Mulher.Ignorado 4544.846 1978.737A 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 bruta mensal
normalmente recebida pode ser facilmente reescrito para médias:
mediarenda <- svymean(x=~E01602, design=dadosPNS, na.rm=TRUE)
mediarenda## mean SE
## E01602 2224.8 35.465E podemos calcular coeficientes de variação e intervalos de confiança da mesma forma:
cv(object=mediarenda)## E01602
## E01602 0.01594085confint(object=mediarenda)## 2.5 % 97.5 %
## E01602 2155.251 2294.27Estimando 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=~C006, design=dadosPNS, na.rm=TRUE)
propsexo## mean SE
## C006Homem 0.47788 0.0011
## C006Mulher 0.52212 0.0011Da mesma forma, podemos estimar a proporção de mais de uma variável:
propsexoraca <- svymean(x=~C006+C009, design=dadosPNS, na.rm=TRUE)
propsexoraca## mean SE
## C006Homem 0.477877269 0.0011
## C006Mulher 0.522122731 0.0011
## C009Branca 0.434361813 0.0033
## C009Preta 0.103948449 0.0016
## C009Amarela 0.007990336 0.0007
## C009Parda 0.448906525 0.0029
## C009Indígena 0.004724780 0.0004
## C009Ignorado 0.000068097 0.0000E estimando a proporção de um cruzamento de duas ou mais variáveis
com a função interaction:
propsexoEraca <- svymean(x=~interaction(C006,C009), design=dadosPNS, na.rm=TRUE)
ftable(x=propsexoEraca)## A B
## interaction(C006, C009)Homem.Branca 0.203649911011 0.001804310734
## interaction(C006, C009)Mulher.Branca 0.230711901753 0.001931770665
## interaction(C006, C009)Homem.Preta 0.051009991023 0.000889604697
## interaction(C006, C009)Mulher.Preta 0.052938458463 0.000954487291
## interaction(C006, C009)Homem.Amarela 0.003597891986 0.000349958045
## interaction(C006, C009)Mulher.Amarela 0.004392443925 0.000388391902
## interaction(C006, C009)Homem.Parda 0.217073098781 0.001610650616
## interaction(C006, C009)Mulher.Parda 0.231833426231 0.001678766457
## interaction(C006, C009)Homem.Indígena 0.002499963631 0.000248470316
## interaction(C006, C009)Mulher.Indígena 0.002224816101 0.000161751415
## interaction(C006, C009)Homem.Ignorado 0.000046412592 0.000018781784
## interaction(C006, C009)Mulher.Ignorado 0.000021684503 0.000009441005Neste 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 interrupção das atividades habituais por problemas respiratórios: ela é a razão entre o total de pessoas que indicaram problemas respiratórios como motivo da interrupção das atividades habituais pelo total de pessoas que interromperam as atividades habituais por motivo de saúde.
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 é:
txprobresp <- svyratio(numerator=~(J00402=="Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"), denominator=~(J002=="Sim"), design=dadosPNS, na.rm=TRUE)## Warning in `[.survey.design2`(design, !nas, ): 4 strata have only one PSU in
## this subset.txprobresp## Ratio estimator: svyratio.survey.design2(numerator = ~(J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"),
## denominator = ~(J002 == "Sim"), design = dadosPNS, na.rm = TRUE)
## Ratios=
## J002 == "Sim"
## J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)" 0.2096767
## SEs=
## J002 == "Sim"
## J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)" 0.005052547Cálculos de coeficiente de variação e intervalos de confiança para essa taxa:
cv(object=txprobresp)## J002 == "Sim"
## J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)" 0.02409685confint(object=txprobresp)## 2.5 %
## J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"/J002 == "Sim" 0.1997739
## 97.5 %
## J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"/J002 == "Sim" 0.2195795Estimando 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 bruta mensal normalmente recebida,
basta utilizar quantiles=0.5:
medianarenda <- svyquantile(x=~E01602, design=dadosPNS, quantiles=0.5, ci=FALSE, na.rm=TRUE)
medianarenda## $E01602
## 0.5
## [1,] 1350
##
## 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=~E01602, design=dadosPNS, quantiles=0.5, ci=TRUE, na.rm=TRUE)
medianarenda## $E01602
## quantile ci.2.5 ci.97.5 se
## 0.5 1350 1300 1400 25.50653
##
## attr(,"hasci")
## [1] TRUE
## attr(,"class")
## [1] "newsvyquantile"SE(object=medianarenda)## E01602
## 25.50653cv(object=medianarenda)## E01602
## 0.01889373Alé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=~E01602, design=dadosPNS, quantiles=c(0.1,0.25,0.5,0.75,0.9), ci=FALSE, na.rm=TRUE)
quantisrenda## $E01602
## 0.1 0.25 0.5 0.75 0.9
## [1,] 400 998 1350 2200 4080
##
## 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.
dadosPNS$variables <- transform(dadosPNS$variables, E01602_real=E01602*Deflator)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=~E01602_real, design=dadosPNS, na.rm=TRUE)
totalrenda_real## total SE
## E01602_real 213227874692 3604489769E para obtenção da média desta nova variável, basta utilizar a função
svymean:
mediarenda_real <- svymean(x=~E01602_real, design=dadosPNS, na.rm=TRUE)
mediarenda_real## mean SE
## E01602_real 2224.8 35.465O processo de deflacionamento é essencial no processo de comparação de uma série histórica dos microdados quando estiver sendo analisada uma variável de rendimento.
Estimação para um domínio
Muitas vezes queremos fazer uma análise para um domínio específico da
população. Com o pacote survey, esse domínio pode ser
selecionado utilizando a função subset no objeto do plano
amostral, aplicando a condição que define aquele domínio. É
necessário que a seleção desse domínio nos dados seja
feita somente após a criação do objeto que define o
plano amostral, caso contrário poderá ser obtido resultados
incorretos.
Para a seleção dos domínios podem ser utilizadas condicionais com igualdade, como por exemplo para estimar a média da renda bruta mensal normalmente recebida de mulheres:
mediarendaM <- svymean(x=~E01602, design=subset(dadosPNS, C006=="Mulher"), na.rm=TRUE)
mediarendaM## mean SE
## E01602 1924.3 30.762Ou condicionais com desigualdade, como por exemplo estimar a taxa de interrupção das atividades habituais por problemas respiratórios para pessoas com idade igual ou superior a 25 anos:
txprobresp25 <- svyratio(numerator=~(J00402=="Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"), denominator=~(J002=="Sim"), design=subset(dadosPNS, C008>=25), na.rm=TRUE)## Warning in `[.survey.design2`(design, !nas, ): 10 strata have only one PSU in
## this subset.txprobresp25## Ratio estimator: svyratio.survey.design2(numerator = ~(J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"),
## denominator = ~(J002 == "Sim"), design = subset(dadosPNS,
## C008 >= 25), na.rm = TRUE)
## Ratios=
## J002 == "Sim"
## J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)" 0.130056
## SEs=
## J002 == "Sim"
## J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)" 0.004654685Alé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=~VDD004A, design=subset(dadosPNS, C006=="Homem" & C009=="Parda" & C008>30), na.rm=TRUE)## Warning in `[.survey.design2`(x, r, ): 3 strata have only one PSU in this
## subset.## Warning in `[.survey.design2`(design, nas == 0, ): 3 strata have only one PSU in
## this subset.nivelinstrHP30## mean SE
## VDD004ASem instrução 0.110647 0.0026
## VDD004AFundamental incompleto ou equivalente 0.394160 0.0045
## VDD004AFundamental completo ou equivalente 0.083547 0.0027
## VDD004AMédio incompleto ou equivalente 0.054271 0.0020
## VDD004AMédio completo ou equivalente 0.248419 0.0040
## VDD004ASuperior incompleto ou equivalente 0.022456 0.0013
## VDD004ASuperior completo 0.086501 0.0026Alé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.
dadosPNS_mulheres <- subset(dadosPNS, C006=="Mulher")
dadosPNS_mulheres## Stratified 1 - level Cluster Sampling design (with replacement)
## With (8028) clusters.
## subset(dadosPNS, C006 == "Mulher")Outro ponto importante de ser mencionado diante de análises por domínio, é sobre a obtenção de indicadores em nível de domicílio, que para analisar corretamente é necessário realizar o recorte na variável de condição no domicílio, conforme exemplo apresentado abaixo da média de pessoas no domicílio:
mediapesdom <- svymean(x=~C001, design=subset(dadosPNS, C004=="Pessoa responsável pelo domicílio"), na.rm=TRUE)
mediapesdom## mean SE
## C001 2.8599 0.0084Estimação para vários domínios
Em outros casos há interesse em estimar quantidades de interesse para
diversos domínios mutuamente exclusivos, a fim de possibilitar
comparações. Para isto, podemos utilizar a função
svyby.
- A variável da qual se deseja calcular a quantidade;
- A variável que define os domínios;
- O objeto do plano amostral;
-
A função utilizada para calcular a quantidade de interesse
(
svytotal,svymean,svyratio,svyquantile, …).
Se desejamos estimar a frequência relativa de homens e mulheres em cada nível de instrução, usamos o seguinte código:
freqSexoInstr <- svyby(formula=~C006, by=~VDD004A, design=dadosPNS, FUN=svymean, na.rm=TRUE)## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 5 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, nas == 0, ): 5 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 2 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, nas == 0, ): 2 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 5 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, nas == 0, ): 5 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 45 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, nas == 0, ): 45 strata have only one PSU
## in this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 27 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, nas == 0, ): 27 strata have only one PSU
## in this subset.freqSexoInstr## VDD004A
## Sem instrução Sem instrução
## 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
## C006Homem C006Mulher se.C006Homem
## Sem instrução 0.4892280 0.5107720 0.004866419
## Fundamental incompleto ou equivalente 0.5013440 0.4986560 0.002143222
## Fundamental completo ou equivalente 0.4919528 0.5080472 0.005334081
## Médio incompleto ou equivalente 0.5078484 0.4921516 0.005893286
## Médio completo ou equivalente 0.4636716 0.5363284 0.003156167
## Superior incompleto ou equivalente 0.4635797 0.5364203 0.007745619
## Superior completo 0.4036566 0.5963434 0.003870122
## se.C006Mulher
## Sem instrução 0.004866419
## Fundamental incompleto ou equivalente 0.002143222
## Fundamental completo ou equivalente 0.005334081
## Médio incompleto ou equivalente 0.005893286
## Médio completo ou equivalente 0.003156167
## Superior incompleto ou equivalente 0.007745619
## Superior completo 0.003870122Caso 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=~VDD004A, by=~C006, design=dadosPNS, FUN=svymean, na.rm=TRUE)
freqInstrSexo## C006 VDD004ASem instrução VDD004AFundamental incompleto ou equivalente
## Homem Homem 0.08530701 0.3712752
## Mulher Mulher 0.08118063 0.3365992
## VDD004AFundamental completo ou equivalente
## Homem 0.07962800
## Mulher 0.07495462
## VDD004AMédio incompleto ou equivalente
## Homem 0.07775848
## Mulher 0.06868542
## VDD004AMédio completo ou equivalente
## Homem 0.2349345
## Mulher 0.2476960
## VDD004ASuperior incompleto ou equivalente VDD004ASuperior completo
## Homem 0.04310541 0.1079915
## Mulher 0.04546368 0.1454205
## se.VDD004ASem instrução se.VDD004AFundamental incompleto ou equivalente
## Homem 0.001275986 0.002694816
## Mulher 0.001254615 0.002414976
## se.VDD004AFundamental completo ou equivalente
## Homem 0.001347127
## Mulher 0.001194810
## se.VDD004AMédio incompleto ou equivalente
## Homem 0.001321463
## Mulher 0.001161618
## se.VDD004AMédio completo ou equivalente
## Homem 0.002458023
## Mulher 0.002120820
## se.VDD004ASuperior incompleto ou equivalente se.VDD004ASuperior completo
## Homem 0.001153861 0.002688042
## Mulher 0.001009624 0.002680386Também pode ser utilizado para estimar a renda média bruta mensal normalmente recebida por Unidade da Federação:
mediaRendaUF <- svyby(formula=~E01602, by=~V0001, design=dadosPNS, FUN=svymean, na.rm=TRUE)
mediaRendaUF## V0001 E01602 se
## Rondônia Rondônia 1885.032 71.64988
## Acre Acre 1742.511 64.25349
## Amazonas Amazonas 1719.220 110.27816
## Roraima Roraima 1927.564 122.10457
## Pará Pará 1466.821 51.59407
## Amapá Amapá 1734.906 92.12348
## Tocantins Tocantins 1817.790 109.95027
## Maranhão Maranhão 1234.153 45.17598
## Piauí Piauí 1391.025 101.42714
## Ceará Ceará 1503.629 89.10509
## Rio Grande do Norte Rio Grande do Norte 1644.863 112.95368
## Paraíba Paraíba 1571.395 107.32236
## Pernambuco Pernambuco 1578.079 72.79566
## Alagoas Alagoas 1560.797 78.40446
## Sergipe Sergipe 1481.270 83.33896
## Bahia Bahia 1429.089 84.62655
## Minas Gerais Minas Gerais 1926.340 68.30364
## Espírito Santo Espírito Santo 2069.731 79.10199
## Rio de Janeiro Rio de Janeiro 2583.055 120.77197
## São Paulo São Paulo 2890.492 125.16684
## Paraná Paraná 2463.707 88.90834
## Santa Catarina Santa Catarina 2547.369 61.69700
## Rio Grande do Sul Rio Grande do Sul 2487.252 85.95990
## Mato Grosso do Sul Mato Grosso do Sul 2396.016 117.05969
## Mato Grosso Mato Grosso 2370.637 97.89929
## Goiás Goiás 2003.799 66.90671
## Distrito Federal Distrito Federal 3740.516 224.96974É possível também calcular o intervalo de confiança para cada uma dessas estimativas:
confint(object=mediaRendaUF)## 2.5 % 97.5 %
## Rondônia 1744.601 2025.463
## Acre 1616.576 1868.445
## Amazonas 1503.079 1935.362
## Roraima 1688.243 2166.884
## Pará 1365.699 1567.944
## Amapá 1554.348 1915.465
## Tocantins 1602.291 2033.288
## Maranhão 1145.610 1322.696
## Piauí 1192.232 1589.819
## Ceará 1328.986 1678.272
## Rio Grande do Norte 1423.478 1866.248
## Paraíba 1361.047 1781.743
## Pernambuco 1435.402 1720.756
## Alagoas 1407.127 1714.467
## Sergipe 1317.929 1644.611
## Bahia 1263.224 1594.954
## Minas Gerais 1792.468 2060.213
## Espírito Santo 1914.694 2224.768
## Rio de Janeiro 2346.347 2819.764
## São Paulo 2645.170 3135.815
## Paraná 2289.450 2637.964
## Santa Catarina 2426.445 2668.293
## Rio Grande do Sul 2318.774 2655.730
## Mato Grosso do Sul 2166.583 2625.449
## Mato Grosso 2178.758 2562.516
## Goiás 1872.664 2134.934
## Distrito Federal 3299.583 4181.449É 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 interrupção das
atividades habituais por problemas respiratórios por sexo e cor ou raça
utilizando as funções svyratio e svyby e o
respectivo cv:
txprobrespSexoRaca <- svyby(formula=~(J00402=="Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"), denominator=~(J002=="Sim"), by=~interaction(C006,C009), design=dadosPNS, FUN=svyratio, vartype="cv", na.rm=TRUE)## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 1 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 94 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 1 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 54 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 30 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 165 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 32 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 138 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 131 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 38 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 132 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 52 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 1 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 67 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 3 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 68 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 125 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 41 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 124 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 58 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 9 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 2 strata have only one PSU in
## this subset.## Warning in `[.survey.design2`(design, byfactor %in% byfactor[i], ): 6 strata
## have only one PSU in this subset.## Warning in `[.survey.design2`(design, !nas, ): 1 strata have only one PSU in
## this subset.txprobrespSexoRaca## interaction(C006, C009)
## 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
## J00402 == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"/J002 == "Sim"
## Homem.Branca 0.2340982
## Mulher.Branca 0.2036258
## Homem.Preta 0.2246148
## Mulher.Preta 0.1882300
## Homem.Amarela 0.1083416
## Mulher.Amarela 0.1537182
## Homem.Parda 0.2292033
## Mulher.Parda 0.1960717
## Homem.Indígena 0.1449440
## Mulher.Indígena 0.2353266
## Homem.Ignorado 0.0000000
## Mulher.Ignorado 0.0000000
## cv
## Homem.Branca 0.05430611
## Mulher.Branca 0.05006883
## Homem.Preta 0.09722969
## Mulher.Preta 0.08769010
## Homem.Amarela 0.51185867
## Mulher.Amarela 0.33776215
## Homem.Parda 0.04561212
## Mulher.Parda 0.03988026
## Homem.Indígena 0.36405478
## Mulher.Indígena 0.33047963
## Homem.Ignorado NaN
## Mulher.Ignorado NaNGráficos para dados amostrais
O pacote survey possui funções específicas para gerar
gráficos que incorporem os pesos amostrais das observações. Para a
criação dos gráficos, é utilizado o mesmo objeto de plano amostral
utilizado para estimar as quantidades populacionais.
Entre os gráficos disponíveis no pacote estão o histograma, boxplot e gráficos de dispersão.
Histograma
A função svyhist permite a criação de histogramas para
variáveis numéricas que consideram os pesos amostrais para computar as
frequências ou densidades das barras.
A sintaxe básica é semelhante
as de estimação apresentadas anteriormente. Além disto, possui os mesmos
parâmetros que a função hist: breaks,
xlab, main, […].
Para estimativas das
frequências absolutas, deve-se utilizar freq=TRUE. Abaixo,
um exemplo do histograma do número de horas normalmente trabalhadas para
cada um dos casos:
svyhist(formula=~as.numeric(E017), design=dadosPNS, main="Histograma", xlab="Número de Horas Normalmente Trabalhadas")
svyhist(formula=~as.numeric(E017), design=dadosPNS, freq=TRUE, main="Histograma", xlab="Número de Horas Normalmente Trabalhadas")
Boxplot
Para a construção de boxplots que considerem os pesos
amostrais, a função é svyboxplot.
A sintaxe dela difere um pouco das demais. É necessário declarar a
variável numérica, sucedida por um ~ e a variável de
agrupamento do boxplot. Caso não deseje usar grupos, basta
colocar o número 1. Além disto, há a opção de plotar apenas
os outliers mais extremos ou todos, através do argumento
all.outliers=TRUE. Exemplos:
svyboxplot(formula=E017~1, design=dadosPNS, main="Boxplot do Número de Horas Normalmente Trabalhadas")
svyboxplot(formula=E017~C006, design=dadosPNS, main="Boxplot do Número de Horas Normalmente Trabalhadas por Sexo")
svyboxplot(formula=E017~C006, design=dadosPNS, all.outliers=TRUE, main="Boxplot do Número de Horas Normalmente 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=E01602~E017, design=dadosPNS, style="bubble", xlab="Número de Horas Normalmente Trabalhadas", ylab="Rendimento Mensal Bruto Normalmente Recebido")
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=E01602~E017, design=dadosPNS, style="transparent", xlab="Número de Horas Normalmente Trabalhadas", ylab="Rendimento Mensal Bruto Normalmente Recebido")
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 brutos mensais normalmente recebidos entre sexos:
svyttest(formula=E01602~C006, design=dadosPNS)##
## Design-based t-test
##
## data: E01602 ~ C006
## t = -16.307, df = 7451, p-value < 0.00000000000000022
## alternative hypothesis: true difference in mean is not equal to 0
## 95 percent confidence interval:
## -603.5173 -473.9895
## sample estimates:
## difference in mean
## -538.7534Outro exemplo para testar diferenças entre médias de horas normalmente trabalhadas, entre pessoas com e sem diagnóstico de alguma doença crônica, física ou mental, ou doença de longa duração:
svyttest(formula=as.numeric(E017)~J007, design=dadosPNS)##
## Design-based t-test
##
## data: as.numeric(E017) ~ J007
## t = 6.7729, df = 7451, p-value = 0.00000000001358
## alternative hypothesis: true difference in mean is not equal to 0
## 95 percent confidence interval:
## 0.7717192 1.4003928
## sample estimates:
## difference in mean
## 1.086056Modelos 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
bruto mensal normalmente recebido 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=E01602~VDD004A+C009+C008, design=dadosPNS)
summary(object=modeloLin)##
## Call:
## svyglm(formula = E01602 ~ VDD004A + C009 + C008, design = dadosPNS)
##
## Survey design:
## survey::postStratify(design = data_prior, strata = ~V00283, population = popc.types)
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) -658.631 93.557 -7.040
## VDD004AFundamental incompleto ou equivalente 420.194 36.784 11.423
## VDD004AFundamental completo ou equivalente 856.624 44.061 19.442
## VDD004AMédio incompleto ou equivalente 1050.529 50.686 20.726
## VDD004AMédio completo ou equivalente 1311.015 55.812 23.490
## VDD004ASuperior incompleto ou equivalente 1790.372 70.420 25.424
## VDD004ASuperior completo 4077.126 106.676 38.220
## C009Preta -670.738 46.133 -14.539
## C009Amarela -117.818 214.938 -0.548
## C009Parda -608.849 41.372 -14.716
## C009Indígena -746.442 97.356 -7.667
## C009Ignorado 363.825 1324.524 0.275
## C008 39.837 1.872 21.285
## Pr(>|t|)
## (Intercept) 0.0000000000020961 ***
## VDD004AFundamental incompleto ou equivalente < 0.0000000000000002 ***
## VDD004AFundamental completo ou equivalente < 0.0000000000000002 ***
## VDD004AMédio incompleto ou equivalente < 0.0000000000000002 ***
## VDD004AMédio completo ou equivalente < 0.0000000000000002 ***
## VDD004ASuperior incompleto ou equivalente < 0.0000000000000002 ***
## VDD004ASuperior completo < 0.0000000000000002 ***
## C009Preta < 0.0000000000000002 ***
## C009Amarela 0.584
## C009Parda < 0.0000000000000002 ***
## C009Indígena 0.0000000000000198 ***
## C009Ignorado 0.784
## C008 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 5416107)
##
## Number of Fisher Scoring iterations: 2Além disso, é possível computar intervalos de confiança para os parâmetros:
confint(object=modeloLin)## 2.5 % 97.5 %
## (Intercept) -842.02949 -475.23153
## VDD004AFundamental incompleto ou equivalente 348.08821 492.30056
## VDD004AFundamental completo ou equivalente 770.25288 942.99555
## VDD004AMédio incompleto ou equivalente 951.17111 1149.88770
## VDD004AMédio completo ou equivalente 1201.60830 1420.42206
## VDD004ASuperior incompleto ou equivalente 1652.32816 1928.41597
## VDD004ASuperior completo 3868.01043 4286.24173
## C009Preta -761.17246 -580.30380
## C009Amarela -539.15746 303.52177
## C009Parda -689.95080 -527.74780
## C009Indígena -937.28735 -555.59706
## C009Ignorado -2232.61649 2960.26564
## C008 36.16794 43.50557Regressão logística
Outro modelo pertencente a classe de modelos generalizados é a
Regressão Logística. Para utilizar este modelo, basta colocar o
argumento family=“binomial”.
No exemplo abaixo modelamos o diagnóstico de alguma doença crônica, física ou mental, ou doença de longa duração pelo sexo, cor ou raça e idade:
modeloLog <- svyglm(formula=J007~C006+C009+C008, design=dadosPNS, family="binomial")## Warning in eval(family$initialize): non-integer #successes in a binomial glm!summary(object=modeloLog)##
## Call:
## svyglm(formula = J007 ~ C006 + C009 + C008, design = dadosPNS,
## family = "binomial")
##
## Survey design:
## survey::postStratify(design = data_prior, strata = ~V00283, population = popc.types)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.1782216 0.0290104 109.554 < 0.0000000000000002 ***
## C006Mulher -0.3486645 0.0158377 -22.015 < 0.0000000000000002 ***
## C009Preta -0.0013958 0.0277866 -0.050 0.960
## C009Amarela -0.2138631 0.2078184 -1.029 0.303
## C009Parda 0.0919759 0.0189782 4.846 0.00000128 ***
## C009Indígena 0.0042507 0.1099412 0.039 0.969
## C009Ignorado 0.0345987 0.6965349 0.050 0.960
## C008 -0.0559756 0.0004777 -117.176 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1.006224)
##
## Number of Fisher Scoring iterations: 4Sendo possível também computar intervalos de confiança para os parâmetros:
confint(object=modeloLog)## 2.5 % 97.5 %
## (Intercept) 3.12135294 3.23509027
## C006Mulher -0.37971095 -0.31761805
## C009Preta -0.05586541 0.05307385
## C009Amarela -0.62124593 0.19351971
## C009Parda 0.05477334 0.12917855
## C009Indígena -0.21126498 0.21976647
## C009Ignorado -1.33080649 1.40000388
## C008 -0.05691205 -0.05503918Aná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)
dadosPNS <- convey_prep(design=dadosPNS)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 bruta mensal
normalmente recebida:
giniBR <- svygini(formula=~E01602, design=dadosPNS, na.rm=TRUE)
giniBR## gini SE
## E01602 0.50802 0.0056Assim como nas funções do pacote survey, é possível, por
exemplo, estimar o coeficiente de variação desta estimativa:
cv(object=giniBR)## E01602
## E01602 0.01094457Também é possível utilizar a função svyby para estimar o
índice de Gini da renda bruta mensal normalmente recebida por Unidade da
Federação:
giniUF <- svyby(formula=~E01602, by=~V0001, design=dadosPNS, FUN=svygini, na.rm=TRUE)
giniUF## V0001 E01602 se
## Rondônia Rondônia 0.4359608 0.013213467
## Acre Acre 0.4776399 0.011391624
## Amazonas Amazonas 0.5124953 0.020344196
## Roraima Roraima 0.4999963 0.017543810
## Pará Pará 0.4835167 0.012171980
## Amapá Amapá 0.5038926 0.013597470
## Tocantins Tocantins 0.5081798 0.018013146
## Maranhão Maranhão 0.5024393 0.012978459
## Piauí Piauí 0.5567044 0.022037635
## Ceará Ceará 0.5324006 0.023822933
## Rio Grande do Norte Rio Grande do Norte 0.5260489 0.020666344
## Paraíba Paraíba 0.5409785 0.019179843
## Pernambuco Pernambuco 0.4934237 0.016054071
## Alagoas Alagoas 0.4748417 0.017280017
## Sergipe Sergipe 0.5341702 0.017287821
## Bahia Bahia 0.5330190 0.017696438
## Minas Gerais Minas Gerais 0.4540088 0.014635889
## Espírito Santo Espírito Santo 0.4738796 0.014263516
## Rio de Janeiro Rio de Janeiro 0.5154285 0.017859004
## São Paulo São Paulo 0.5101548 0.013463160
## Paraná Paraná 0.4535436 0.011302062
## Santa Catarina Santa Catarina 0.3941068 0.009037718
## Rio Grande do Sul Rio Grande do Sul 0.4683207 0.012028102
## Mato Grosso do Sul Mato Grosso do Sul 0.4790102 0.019008431
## Mato Grosso Mato Grosso 0.4503289 0.014585596
## Goiás Goiás 0.4188121 0.009881148
## Distrito Federal Distrito Federal 0.5396514 0.011241314E estimar intervalos de confiança para cada UF:
confint(object=giniUF)## 2.5 % 97.5 %
## Rondônia 0.4100629 0.4618587
## Acre 0.4553127 0.4999671
## Amazonas 0.4726214 0.5523692
## Roraima 0.4656111 0.5343816
## Pará 0.4596601 0.5073734
## Amapá 0.4772420 0.5305431
## Tocantins 0.4728747 0.5434849
## Maranhão 0.4770020 0.5278766
## Piauí 0.5135114 0.5998973
## Ceará 0.4857085 0.5790927
## Rio Grande do Norte 0.4855436 0.5665541
## Paraíba 0.5033867 0.5785703
## Pernambuco 0.4619583 0.5248891
## Alagoas 0.4409735 0.5087099
## Sergipe 0.5002867 0.5680537
## Bahia 0.4983346 0.5677034
## Minas Gerais 0.4253229 0.4826946
## Espírito Santo 0.4459236 0.5018356
## Rio de Janeiro 0.4804255 0.5504315
## São Paulo 0.4837675 0.5365421
## Paraná 0.4313919 0.4756952
## Santa Catarina 0.3763931 0.4118204
## Rio Grande do Sul 0.4447460 0.4918953
## Mato Grosso do Sul 0.4417544 0.5162660
## Mato Grosso 0.4217417 0.4789162
## Goiás 0.3994454 0.4381788
## Distrito Federal 0.5176188 0.5616839Curva 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 bruta mensal normalmente recebida:
curvaLorenz <- svylorenz(formula=~E01602, design=dadosPNS, quantiles=seq(from=0, to=1, by=0.05), na.rm=TRUE)
Instituto Brasileiro de Geografia e Estatística - pacotesipd@ibge.gov.br↩︎