Tutorial SINAPE
Natalia R. S. Pires e Pedro L. N. Silva
27 de setembro de 2018
Carregando Dados da PNAD Contínua
Antes de mais nada, é necessária a instalação dos pacotes que serão utilizados ao longo da análise dos microdados da PNAD Contínua.
# Limpa objetos da memória
rm(list = ls())
# Define pasta onde estão ou serão armazenados arquivos nesta sessão
setwd("F:/SINAPE/Tutorial")
#setwd("F:/SkyDrive/RJ/Eventos/SINAPE2018/Tutorial")
# Opção para permitir estimar variância quando houver apenas 1 observação na UPA
options( survey.lonely.psu = "adjust" )
# Troca . por , nos outputs
options(OutDec=",")
# Instala e Carrega pacotes necessários
# Pacote para facilitar manuseio dos dados
#install.packages("tidyverse")
library(tidyverse)## -- Attaching packages ------------------------------------------------------- tidyverse 1.2.1 --## v ggplot2 2.2.1 v purrr 0.2.4
## v tibble 1.4.2 v dplyr 0.7.4
## v tidyr 0.8.1 v stringr 1.2.0
## v readr 1.1.1 v forcats 0.2.0## -- Conflicts ---------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag() masks stats::lag()# Pacote para analisar dados amostrais
#install.packages("survey")
library(survey)## Loading required package: grid## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following object is masked from 'package:tidyr':
##
## expand## Loading required package: survival##
## Attaching package: 'survey'## The following object is masked from 'package:graphics':
##
## dotchart# Pacote para carregar microdados da PNADc
#install.packages("PNADcIBGE")
library(PNADcIBGE)
# Pacote para analisar dados amostrais
#install.packages("srvyr")
library(srvyr)##
## Attaching package: 'srvyr'## The following object is masked from 'package:stats':
##
## filterLogo, precisamos carregar a base de microdados.
###########################################
# Carrega microdados no diretório temporário
# Usa pacote PNADcIBGE
pnadc_dat<-
get_pnadc(2018, quarter = 2, interview = NULL, #informa ano e trmestre de interesse
vars = NULL,labels = F, #informa quais e como ver as variáveis
design = F, #informa o objeto para análise (base de microdados ou objeto com desenho amostral)
savedir = tempdir() #informa diretório para salvar arquivo
)
# Adiciona coluna de 1's ao arquivo de microdados
pnadc_dat$one <- 1
# Conta o numero de pessoas da amostra
sum(pnadc_dat$one)## [1] 556186Com o arquivo de microdados disponível, sabendo que o mesmo é resultado de uma pesquisa amostral, o próximo passo é construir um objeto de dados que permita o cálculo de estimativas considerando a estrutura do plano amostral complexo.
###########################################
# Declara Plano Amostral
# Usa pacote survey
#Declara estrutura do plano amostral complexo
pnadc_plano <-
svydesign(
ids = ~ UPA , # Declara a unidade amostral mais granular
strata = ~ Estrato , # Declara a variável que contém os estratos
weights = ~ V1027 , # Declara variável com pesos
data = pnadc_dat , # Declara base de microdados
nest = TRUE # Declara que os estratos podem conter identificações identicas para UPA's distintas
)
# Sumariza o objeto contendo as informações sobre estrutura do plano amostral
#summary(pnadc_plano)O objeto pnadc_plano estaria pronto para uso na produção de inferências se não considerássemos o fato de que o IBGE recomenda a calibração dos pesos com base nas estimativas de população produzidas pela Fundação. Como é desejável que as estimativas de indicadores sejam compatíveis com as informações oficiais divulgadas pelo IBGE - e também porque esta informação adicional não deve ser descartada, construímos um novo objeto, agora com os pesos devidamente calibrados.
###########################################
# Tabela com frequencias populacionais (estimativas IBGE para calibração)
df_pos <- data.frame( posest = unique( pnadc_dat$posest ) , Freq = unique( pnadc_dat$V1029 ) )
df_pos## posest Freq
## 1 111 526181
## 2 114 1285826
## 3 124 439068
## 4 121 388828
## 5 131 2160249
## 6 134 1447876
## 7 132 360879
## 8 141 333362
## 9 144 137365
## 10 154 6149358
## 11 151 1457373
## 12 152 802705
## 13 161 482833
## 14 164 205318
## 15 162 116898
## 16 174 1259416
## 17 171 292734
## 18 214 5382262
## 19 211 1097895
## 20 212 355392
## 21 213 168532
## 22 224 2185434
## 23 221 852812
## 24 223 186921
## 25 231 2643393
## 26 234 5126108
## 27 232 1297819
## 28 244 1989283
## 29 241 891743
## 30 242 653908
## 31 254 2693228
## 32 251 819612
## 33 252 512008
## 34 264 5481277
## 35 262 2346593
## 36 261 1640608
## 37 274 2099475
## 38 272 248167
## 39 271 1035501
## 40 284 1345276
## 41 282 304410
## 42 281 657606
## 43 292 1071020
## 44 294 11342509
## 45 291 2967274
## 46 314 15757458
## 47 311 2532592
## 48 312 2813630
## 49 313 109270
## 50 322 1605719
## 51 321 375405
## 52 324 2068181
## 53 332 5764170
## 54 334 4487905
## 55 331 6538995
## 56 351 12166814
## 57 352 9354644
## 58 354 23870136
## 59 414 7828318
## 60 412 1626742
## 61 411 1920723
## 62 424 6038990
## 63 421 492693
## 64 422 535582
## 65 434 7089110
## 66 432 2760547
## 67 431 1488513
## 68 504 1788712
## 69 501 882919
## 70 514 2426624
## 71 511 593586
## 72 512 307475
## 73 524 3090982
## 74 522 1042623
## 75 523 1235582
## 76 521 1481393
## 77 531 3090833###########################################
# Calibra pesos
pnadc_calib <- postStratify( pnadc_plano , ~ posest , df_pos )
# Outras opções: funções rake() e calibrate()
# Obtém fatores de calibração dos pesos da amostra
pnadc_fatores = weights(pnadc_calib) / weights(pnadc_plano)
boxplot(pnadc_fatores, horizontal = TRUE, xlab="Fatores de calibração")
Finalmente, o objeto pnadc_calib pode ser utilizado para construir estimativas a cerca da população brasileira. Uma forma de validar os procedimentos adotados até aqui é comparar a estimativa populacional encontrada utilizando o objeto pnadc_calib com a estimativa divulgada pelo IBGE (https://sidra.ibge.gov.br/home/pnadct/brasil).
###########################################
# Valida estimativas populacionais
svytotal( ~ one , pnadc_calib ) #estimativa populacional## total SE
## one 208409201 0Com as estimativas devidamente validadas, pode-se avaliar então o Efeito do Plano Amostral (EPA) a partir de alguma estatística de interesse. Considerando que a PNADc é uma referência para estatísticas relacionadas ao mercado de trabalho, podemos verificar o EPA para o indicador “Rendimento Habitual do Trabalho Principal”.
###########################################
# Avalia Efeito do Plano Amostral
# Compara estimativas do RHTP com e sem plano (+ calculo do EPA de Kish)
round(svymean( ~ VD4016 , subset(pnadc_plano, V2009 >= 14 & VD4015 %in% 1) , na.rm = TRUE ),2)## mean SE
## VD4016 2148,6 27,517round(svymean( ~ VD4016 , subset(pnadc_calib, V2009 >= 14 & VD4015 %in% 1) , na.rm = TRUE, deff="replace"),2)## mean SE DEff
## VD4016 2127,620 25,918 13,831Finalizada a etapa de criar um objeto que leve em consideração o plano amostral utilizado na pesquisa, podemos salvá-lo em algum diretório de interesse para posterior uso.
###########################################
# Salva objeto final
saveRDS(pnadc_calib,"pnadc_calib_20182")
# Limpa objetos da memória
rm(pnadc_dat,pnadc_plano)Estimando estatísticas de interesse
A partir do objeto R contendo os microdados e detalhes da estrutura do plano amostral, pode-se construir estimativas a partir de diferentes estatísticas de interesse sobre as condições de vida dos brasileiros. Para tanto, fazemos uso dos pacotes tidyverse para construir novas variáveis, fazer filtros, etc. e do pacote survey para realizar estimativas considerando o plano amostral.
###################################################
# Se a sessão do R foi encerrada e se deseja carregar o objeto com os microdados
# e o plano amostral, é necessário rodar as funções abaixo
# Limpa objetos da memória
rm(list = ls())
# Define pasta onde estão ou serão armazenados arquivos nesta sessão
setwd("F:/SINAPE/Tutorial")
#setwd("F:/SkyDrive/RJ/Eventos/SINAPE2018/Tutorial")
# Troca . por , nos outputs
options(OutDec=",")
# Carrega pacotes necessários
library(tidyverse)
library(survey)
library(srvyr)
###################################################
# Carrega dados da PNADC
pnadc_calib <- readRDS(file="pnadc_calib_20182")
# Modifica objeto para permitir sintaxe tipo tidyverse
pnadc_calib <- as_survey_design(pnadc_calib)Após selecionar os indicadores de interesse para calcular as estimativas, é necessário realizar as devidas manipulações algébricas para preparar as variáveis para uso nas funções.
###################################################
# Prepara variáveis para cálculo de estimativas
pnadc_calib <- update(pnadc_calib,
idade5 = factor( 1 + findInterval( V2009 , seq( 5 , 60 , 5 ) ) ) ,
nivel_renda = factor( 1 + findInterval( VD4019 , seq( 5 , 60 , 5 ))),
sexo = as.numeric( V2007 == 1 ) ,
pia = as.numeric( V2009 >= 14 ) ,
analfabeto = 1*(V3001==2),
ocupado = ifelse( pia == 1 , as.numeric( VD4002 %in% 1 ) , NA ) ,
desocup30 = ifelse( pia == 1 , as.numeric( VD4002 %in% 2 ) , NA ) ,
pea_c = as.numeric( ocupado == 1 | desocup30 == 1 ),
# (rendimento habitual do trabalho principal)
VD4016n = ifelse( pia %in% 1 & VD4015 %in% 1 , VD4016 , NA ) ,
# (rendimento habitual de todos os trabalhos)
VD4019n = ifelse( pia %in% 1 & VD4015 %in% 1 , VD4019 , NA ) ,
# (rendimento efetivo do todos os trabalhos)
VD4020n = ifelse( pia %in% 1 & VD4015 %in% 1 , VD4020 , NA ) ,
#indicador de nível superior
VD3001n = 1*( VD3001 == 7)
) Feito isto, se pode prosseguir com o cálculo das estimativas de interesse. O pacote survey permite o cálculo de estimativas pontuais através de funções como a svymean, svytotal, svyquantile etc. Por default, as funções retornam além das estimativas, o Erro padrão associado ao cálculo.
###################################################
# Estimativas pontuais
###################################################
# Calcula estimativas do rendimento médio nominal de todos os trabalhos - Tabela 5429
rend_tot<-svymean( ~ VD4020n + VD4019n,
subset(pnadc_calib, VD4002 %in% 1) ,
na.rm = TRUE )
round(coef(rend_tot))## VD4020n VD4019n
## 2189 2199###################################################
# Calcula estimativas do rendimento efetivo médio nominal de todos os trabalhos por grupos
rend_grupo <- svyby( ~ VD4020n ,
~ sexo ,
subset(pnadc_calib, VD4002 %in% 1),
svymean , na.rm = TRUE )
round(rend_grupo)## sexo VD4020n se
## 0 0 1864 22
## 1 1 2436 34###################################################
# Calcula estimativas de quantis selecionados
rend_quantil <- svyquantile( ~ VD4020n ,
subset(pnadc_calib, VD4002 %in% 1),
c(0.25, 0.5, 0.75, 0.9, 0.95, 0.99) ,
na.rm = TRUE )
round(rend_quantil)## 0,25 0,5 0,75 0,9 0,95 0,99
## VD4020n 954 1300 2200 4000 6900 16000###################################################
# Estimativas de proporção
###################################################
# Calcula estimativa da Taxa de Analfabetismo
Taxa_analf <- svymean(~analfabeto,
subset(pnadc_calib, V2009 >= 15),
na.rm = TRUE)
round(100*coef(Taxa_analf),2)## analfabeto
## 6,8round(100*SE(Taxa_analf),2)## analfabeto
## analfabeto 0,06###################################################
# Calcula estimativa da Taxa de Desocupação (razão de dois totais populacionais) - Tabela 4099
Taxa_desocup <- svyratio( numerator = ~ desocup30 ,
denominator = ~ pea_c ,
pnadc_calib ,
na.rm = TRUE )
round(100*coef(Taxa_desocup),2)## desocup30/pea_c
## 12,44round(100*SE(Taxa_desocup),2)## desocup30/pea_c
## 0,12Também é possível calcular indicadores que relacionam mais de uma variável com o uso das funções svyratio, svytable e svychisq.
###################################################
# Tabelas de contingência
###################################################
# Relação de Contribuição para a previdência e Nível de ensino
tabela_cont <- svytable(~VD4012+VD3001n,
pnadc_calib)
tabela_cont## VD3001n
## VD4012 0 1
## 1 42708989 15369814
## 2 30634163 2524369###################################################
# Relação de Contribuição para a previdência e Nível de ensino
tabela_perc <- svytable(~VD4012+VD3001n,
pnadc_calib,
Ntotal=100)
round(tabela_perc,2) ## VD3001n
## VD4012 0 1
## 1 46,81 16,85
## 2 33,58 2,77###################################################
# Teste Qui-quadrado de Pearson com ajuste de Rao-Scott para avaliar homogeneidade
teste_qui <- svychisq(~VD4012+VD3001n,
pnadc_calib,
statistic = "Chisq")
teste_qui##
## Pearson's X^2: Rao & Scott adjustment
##
## data: NextMethod()
## X-squared = 29001, df = 1, p-value < 2,2e-16###################################################
# Distribuição condicional da Contribuição para a previdência dado Nível de ensino
tabela_cond <- svyby(~VD4012,
~VD3001n,
pnadc_calib,
svymean, na.rm=TRUE)
round(100*tabela_cond, 1)## VD3001n VD40121 VD40122 se.VD40121 se.VD40122
## 0 0 58,2 41,8 0,2 0,2
## 1 100 85,9 14,1 0,3 0,3Para salvar as modificações realizadas no objeto pnadc_calib, utilizamos novamente a função save.rds
###################################################
# Salva objeto final
saveRDS(pnadc_calib,"pnadc_calib_20182")
# Limpa objetos da memória
rm(list = ls())Modelagem MPV
A outra abordagem possível para análise dos microdados da PNADC é o uso de modelos de regressão linear para representar relações entre variáveis. Isto é possível com o uso da função svyglm.
Uma variável de interesse do mercado de trabalho é o RHTP (Rendimento Habitual do Trabalho Principal por Hora Trabalhada). É possível, com base em um conhecimento mínimo da sociedade brasileira, fazer algumas suposiçõs sobre a variabilidade desta variável em razão de diferenças de idade, sexo e posto de trabalho. A significancia destas variáveis em pesquisas amostrais podem ser testadas utilizando Modelos de Regressão por Máxima Pseudo-Verossimilhança. Para tanto, é necessário conhecer a distribuição da variável de interesse.
###################################################
# Carrega dados da PNADC
pnadc_calib <- readRDS(file="pnadc_calib_20182")
# Modifica objeto para permitir sintaxe tipo tidyverse
pnadc_calib <- as_survey_design(pnadc_calib)
# Transformação de variáveis selecionadas para uso na modelagem
pnadc_calib <- update(pnadc_calib,
Sexo = as.factor(V2007),
CorRaca = as.factor(V2010),
TipoDomic = as.factor(V1023),
TempoTrab = as.factor(V4040),
Escolaridade = as.factor(VD3001),
TipoEmpr = as.factor(VD4008),
Gr_AtvPrinc = as.factor(VD4010),
Gr_Ocup = as.factor(VD4011),
QtdeTrab = as.factor(V4009),
Idade=V2009,
RHTP = 12*VD4016/(52*VD4031)
)
###################################################
# Seleciona pessoas com idades entre 18 e 65 anos de idade
pnadc_calib_f<- subset(pnadc_calib, V2009 >= 18 & V2009 <= 65 )
pnadc_calib_f<-as_survey_design(pnadc_calib_f)
# Salva novo arquivo de interesse e libera memória
saveRDS(pnadc_calib_f,"pnadc_calib_f")
rm(pnadc_calib)
###################################################
# Elabora histograma da distribuição dos dados
svyhist(~RHTP, pnadc_calib_f)
# Aplica transformação log para simetrizar distribuição
pnadc_calib_f <- mutate(pnadc_calib_f, logRHTP = ifelse(RHTP==0, NA, log(RHTP)))
###################################################
# Elabora histograma da distribuição dos dados log-transformados
svyhist(~logRHTP,pnadc_calib_f)
# Estima parâmetros da distribuição para comparar distribuições observada e esperada
media_logRHTP <- svymean(~logRHTP, pnadc_calib_f, na.rm=TRUE)
media_logRHTP## mean SE
## logRHTP 2,1378 0,0053var_logRHTP <- svyvar(~logRHTP, pnadc_calib_f, na.rm=TRUE)
var_logRHTP## variance SE
## logRHTP 0,68267 0,0084###################################################
# Histograma da distribuição dos dados log-transformados
# com sobreposição da curva normal ajustada
svyhist(~logRHTP, subset(pnadc_calib_f,!is.na(logRHTP)),
main="Distribuições Observada e Esperada do RHTP")
curve(dnorm(x, mean=media_logRHTP, sd=sqrt(var_logRHTP)), col=2, lty=2, lwd=2, add=TRUE)
Dado o pressuposto da distribuição Normal dos dados, se pode seguir com o ajuste do modelo de regressão para relacionar o RHTP com as variáveis explicativas de interesse.
###################################################
# Ajusta modelo de regressão saturado para RHTP
modelo0_fit <- svyglm(logRHTP~
Sexo+
CorRaca+
Idade+
TipoDomic+
QtdeTrab+
Escolaridade+
TipoEmpr+
Gr_AtvPrinc+
Gr_Ocup+
TempoTrab,
pnadc_calib_f)
# Obtém resumo do modelo ajustado
summary(modelo0_fit)##
## Call:
## svyglm(formula = logRHTP ~ Sexo + CorRaca + Idade + TipoDomic +
## QtdeTrab + Escolaridade + TipoEmpr + Gr_AtvPrinc + Gr_Ocup +
## TempoTrab, pnadc_calib_f)
##
## Survey design:
## Called via srvyr
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1,5378352 0,0283034 54,334 < 2e-16 ***
## Sexo2 -0,1866689 0,0043820 -42,599 < 2e-16 ***
## CorRaca2 -0,1871554 0,0074399 -25,156 < 2e-16 ***
## CorRaca3 0,0074607 0,0298472 0,250 0,802621
## CorRaca4 -0,1800135 0,0045955 -39,171 < 2e-16 ***
## CorRaca5 -0,1141735 0,0306868 -3,721 0,000199 ***
## CorRaca9 -0,4956935 0,1650800 -3,003 0,002680 **
## Idade 0,0092683 0,0002089 44,372 < 2e-16 ***
## TipoDomic2 -0,0892909 0,0091639 -9,744 < 2e-16 ***
## TipoDomic3 -0,1130261 0,0237806 -4,753 2,02e-06 ***
## TipoDomic4 -0,1726534 0,0072169 -23,924 < 2e-16 ***
## QtdeTrab2 -0,3141366 0,0121411 -25,874 < 2e-16 ***
## QtdeTrab3 -0,1504797 0,0594234 -2,532 0,011341 *
## Escolaridade2 0,2854370 0,0129581 22,028 < 2e-16 ***
## Escolaridade3 0,4416197 0,0144271 30,610 < 2e-16 ***
## Escolaridade4 0,4685677 0,0148962 31,455 < 2e-16 ***
## Escolaridade5 0,5351541 0,0136449 39,220 < 2e-16 ***
## Escolaridade6 0,6700609 0,0158106 42,380 < 2e-16 ***
## Escolaridade7 0,9994329 0,0161532 61,872 < 2e-16 ***
## TipoEmpr2 0,1347386 0,0139686 9,646 < 2e-16 ***
## TipoEmpr3 0,1250814 0,0101939 12,270 < 2e-16 ***
## TipoEmpr4 0,3372769 0,0135222 24,942 < 2e-16 ***
## TipoEmpr5 -0,1542993 0,0065317 -23,623 < 2e-16 ***
## Gr_AtvPrinc02 0,1331981 0,0131476 10,131 < 2e-16 ***
## Gr_AtvPrinc03 0,1931731 0,0133866 14,430 < 2e-16 ***
## Gr_AtvPrinc04 0,0566441 0,0122904 4,609 4,08e-06 ***
## Gr_AtvPrinc05 0,2062809 0,0140270 14,706 < 2e-16 ***
## Gr_AtvPrinc06 0,0291373 0,0147448 1,976 0,048161 *
## Gr_AtvPrinc07 0,1820679 0,0134425 13,544 < 2e-16 ***
## Gr_AtvPrinc08 0,2490840 0,0175163 14,220 < 2e-16 ***
## Gr_AtvPrinc09 0,0409851 0,0154720 2,649 0,008082 **
## Gr_AtvPrinc10 0,1491716 0,0154564 9,651 < 2e-16 ***
## Gr_AtvPrinc12 0,2327309 0,0970710 2,398 0,016519 *
## Gr_Ocup02 0,0202694 0,0169015 1,199 0,230444
## Gr_Ocup03 -0,2195514 0,0170524 -12,875 < 2e-16 ***
## Gr_Ocup04 -0,4525589 0,0168992 -26,780 < 2e-16 ***
## Gr_Ocup05 -0,4744503 0,0158018 -30,025 < 2e-16 ***
## Gr_Ocup06 -0,6694843 0,0213864 -31,304 < 2e-16 ***
## Gr_Ocup07 -0,4256619 0,0178232 -23,882 < 2e-16 ***
## Gr_Ocup08 -0,4154008 0,0174655 -23,784 < 2e-16 ***
## Gr_Ocup09 -0,5694208 0,0166977 -34,102 < 2e-16 ***
## Gr_Ocup10 -0,0136013 0,0254519 -0,534 0,593078
## Gr_Ocup11 -0,1073419 0,2404695 -0,446 0,655326
## TempoTrab2 0,1469671 0,0125001 11,757 < 2e-16 ***
## TempoTrab3 0,1964845 0,0132021 14,883 < 2e-16 ***
## TempoTrab4 0,3197096 0,0123808 25,823 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0,2374139)
##
## Number of Fisher Scoring iterations: 2anova(modelo0_fit)## Anova table: (Rao-Scott LRT)
## svyglm(formula = logRHTP ~ Sexo, pnadc_calib_f)
## stats DEff df ddf p
## Sexo 391,00 1,5685 1,0000 14506 < 2,2e-16 ***
## CorRaca 14556,70 3,2177 5,0000 14501 < 2,2e-16 ***
## Idade 3279,36 2,7413 1,0000 14500 < 2,2e-16 ***
## TipoDomic 10261,66 6,9457 3,0000 14497 < 2,2e-16 ***
## QtdeTrab 110,82 3,8179 2,0000 14495 4,749e-06 ***
## Escolaridade 56758,65 1,8794 6,0000 14489 < 2,2e-16 ***
## TipoEmpr 7140,02 1,9380 4,0000 14485 < 2,2e-16 ***
## Gr_AtvPrinc 3468,07 1,6955 10,0000 14475 < 2,2e-16 ***
## Gr_Ocup 6832,02 1,6547 10,0000 14465 < 2,2e-16 ***
## TempoTrab 2226,45 1,2964 3,0000 14462 < 2,2e-16 ***
## ---
## Signif. codes: 0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1No caso do modelo ajustado não ter variáveis ou mesmo fatores das variáveis significativos, é necessária uma revisão dos parâmetros para estimação, seja através de filtros ou mesmo de combinação de fatores e/ou variáveis.
###################################################
### Ajuste de modelo revisado
# Constroi variáveis agrupadas
pnadc_calib_f <- update(pnadc_calib_f,
CorRacaG = as.factor(ifelse(CorRaca=="3"|CorRaca=="1" , 1, CorRaca)),
Gr_OcupG = as.factor(ifelse(Gr_Ocup=="01"| Gr_Ocup=="02" | Gr_Ocup=="10",1, Gr_Ocup))
)
###################################################
# Ajusta modelo para RHTP após recodificar variáveis
modelo2_fit <- svyglm(logRHTP~
Sexo+
CorRacaG+
Idade+
TipoDomic+
Escolaridade+
TipoEmpr+
Gr_AtvPrinc+
Gr_OcupG+
TempoTrab+
QtdeTrab,
pnadc_calib_f)
###################################################
# Obtém resumo do modelo ajustado
summary(modelo2_fit)##
## Call:
## svyglm(formula = logRHTP ~ Sexo + CorRacaG + Idade + TipoDomic +
## Escolaridade + TipoEmpr + Gr_AtvPrinc + Gr_OcupG + TempoTrab +
## QtdeTrab, pnadc_calib_f)
##
## Survey design:
## update(pnadc_calib_f, CorRacaG = as.factor(ifelse(CorRaca ==
## "3" | CorRaca == "1", 1, CorRaca)), Gr_OcupG = as.factor(ifelse(Gr_Ocup ==
## "01" | Gr_Ocup == "02" | Gr_Ocup == "10", 1, Gr_Ocup)))
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1,5472873 0,0246550 62,757 < 2e-16 ***
## Sexo2 -0,1862979 0,0043710 -42,622 < 2e-16 ***
## CorRacaG2 -0,1873176 0,0074427 -25,168 < 2e-16 ***
## CorRacaG4 -0,1801600 0,0046080 -39,097 < 2e-16 ***
## CorRacaG5 -0,1141809 0,0306593 -3,724 0,000197 ***
## CorRacaG6 -0,4987752 0,1653647 -3,016 0,002564 **
## Idade 0,0092695 0,0002097 44,196 < 2e-16 ***
## TipoDomic2 -0,0892145 0,0091655 -9,734 < 2e-16 ***
## TipoDomic3 -0,1132609 0,0237409 -4,771 1,85e-06 ***
## TipoDomic4 -0,1724679 0,0072041 -23,940 < 2e-16 ***
## Escolaridade2 0,2853020 0,0129569 22,019 < 2e-16 ***
## Escolaridade3 0,4414583 0,0144246 30,605 < 2e-16 ***
## Escolaridade4 0,4683862 0,0148955 31,445 < 2e-16 ***
## Escolaridade5 0,5348396 0,0136443 39,199 < 2e-16 ***
## Escolaridade6 0,6700850 0,0158053 42,396 < 2e-16 ***
## Escolaridade7 1,0025195 0,0160398 62,502 < 2e-16 ***
## TipoEmpr2 0,1349774 0,0139757 9,658 < 2e-16 ***
## TipoEmpr3 0,1250918 0,0102403 12,216 < 2e-16 ***
## TipoEmpr4 0,3349823 0,0131793 25,417 < 2e-16 ***
## TipoEmpr5 -0,1537186 0,0065603 -23,431 < 2e-16 ***
## Gr_AtvPrinc02 0,1331351 0,0131511 10,124 < 2e-16 ***
## Gr_AtvPrinc03 0,1933266 0,0133785 14,451 < 2e-16 ***
## Gr_AtvPrinc04 0,0565521 0,0122812 4,605 4,16e-06 ***
## Gr_AtvPrinc05 0,2059624 0,0140432 14,666 < 2e-16 ***
## Gr_AtvPrinc06 0,0288236 0,0147194 1,958 0,050226 .
## Gr_AtvPrinc07 0,1833588 0,0133784 13,706 < 2e-16 ***
## Gr_AtvPrinc08 0,2460176 0,0171170 14,373 < 2e-16 ***
## Gr_AtvPrinc09 0,0448273 0,0153155 2,927 0,003429 **
## Gr_AtvPrinc10 0,1498325 0,0154556 9,694 < 2e-16 ***
## Gr_AtvPrinc12 0,2325578 0,0970369 2,397 0,016561 *
## Gr_OcupG3 -0,2306107 0,0102744 -22,445 < 2e-16 ***
## Gr_OcupG4 -0,4633579 0,0097524 -47,512 < 2e-16 ***
## Gr_OcupG5 -0,4843354 0,0096991 -49,936 < 2e-16 ***
## Gr_OcupG6 -0,6790795 0,0173600 -39,117 < 2e-16 ***
## Gr_OcupG7 -0,4351443 0,0117774 -36,947 < 2e-16 ***
## Gr_OcupG8 -0,4245759 0,0115500 -36,760 < 2e-16 ***
## Gr_OcupG9 -0,5791124 0,0105310 -54,991 < 2e-16 ***
## Gr_OcupG11 -0,1174815 0,2399639 -0,490 0,624439
## TempoTrab2 0,1469333 0,0124989 11,756 < 2e-16 ***
## TempoTrab3 0,1963310 0,0131938 14,881 < 2e-16 ***
## TempoTrab4 0,3194163 0,0123765 25,808 < 2e-16 ***
## QtdeTrab2 -0,3132691 0,0121605 -25,761 < 2e-16 ***
## QtdeTrab3 -0,1482312 0,0596260 -2,486 0,012929 *
## ---
## Signif. codes: 0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0,2374241)
##
## Number of Fisher Scoring iterations: 2###################################################
#Obtém ANOVA
anova(modelo2_fit)## Anova table: (Rao-Scott LRT)
## svyglm(formula = logRHTP ~ Sexo, pnadc_calib_f)
## stats DEff df ddf p
## Sexo 391,0 1,5685 1,0000 14506 < 2,2e-16 ***
## CorRacaG 14463,2 2,6917 4,0000 14502 < 2,2e-16 ***
## Idade 3286,9 2,7523 1,0000 14501 < 2,2e-16 ***
## TipoDomic 10317,6 6,9804 3,0000 14498 < 2,2e-16 ***
## Escolaridade 56073,6 1,9016 6,0000 14492 < 2,2e-16 ***
## TipoEmpr 6952,6 1,9358 4,0000 14488 < 2,2e-16 ***
## Gr_AtvPrinc 3555,0 1,7110 10,0000 14478 < 2,2e-16 ***
## Gr_OcupG 6684,3 1,6027 8,0000 14470 < 2,2e-16 ***
## TempoTrab 2220,4 1,3028 3,0000 14467 < 2,2e-16 ***
## QtdeTrab 1074,4 2,3444 2,0000 14465 < 2,2e-16 ***
## ---
## Signif. codes: 0 '***' 0,001 '**' 0,01 '*' 0,05 '.' 0,1 ' ' 1Para avaliar a qualidade do ajuste, além da significância dos parâmetros, pode-se realizar a análise dos resíduos.
###################################################
# Obtém estatística AIC dos modelos ajustados
AIC(modelo2_fit, modelo0_fit)## eff.p AIC deltabar
## [1,] 68,27157 132187,9 1,625514
## [2,] 74,44049 132194,6 1,654233###################################################
# Analisa resíduos
plot(modelo2_fit$residuals, main="a.Dispersão dos resíduos") 
# Analisa homoscedasticidade dos resíduos
plot(modelo2_fit$fitted.values, modelo2_fit$residuals, main="b.Variância dos resíduos")
# Limpa objetos da memória
rm(list = ls())Modelos de regressão logística também são avaliáveis no pacote survey. Teoria e aplicações de Análise de Dados Amostrais Complexos podem ser avaliados na página https://djalmapessoa.github.io/adac/index.html.