Análise de Dados PNS 2019 com R

Importando bibliotecas

#Importando bibliotecas
library(dplyr)
library(PNSIBGE)
library(survey)
library(convey)
library(ggplot2)

Baixando e tratando dados/variáveis (Parte 1)

#Baixando e tratando dados/variáveis

Colunas <- c("V0001","C001","C004","C006","C008","C009","E01602","E017","J002","J00402","J007","VDD004A", "I00102")

Dicionário de variáveis: https://basedosdados.org/dataset/86bac6cc-575f-4289-a857-13f3f52c9a1d?table=f8b6030c-3fb1-4f64-81ef-3c5c3c888b6c

Base <- get_pns(year=2019, vars=Colunas)

Base
Stratified 1 - level Cluster Sampling design (with replacement)
With (8028) clusters.
survey::postStratify(design = data_prior, strata = ~V00283, population = popc.types)

Baixando e tratando dados/variáveis (Parte 2)

Base$variables <- transform(Base$variables, Idade=C008)
Base$variables <- transform(Base$variables, Sexo=C006)
Base$variables <- transform(Base$variables, HorasSemanaisTrabalho=E017)
Base$variables <- transform(Base$variables, Raça=C009)
Base$variables <- transform(Base$variables, Rendimento_bruto=E01602)
Base$variables <- transform(Base$variables, Doença=J007)
Base$variables <- transform(Base$variables, Motivo_Saude=J00402)
Base$variables <- transform(Base$variables, Interrupção_Atividade=J002)
Base$variables <- transform(Base$variables, Nível_Instrução=VDD004A)
Base$variables <- transform(Base$variables, TemPlano=I00102)

Baixando e tratando dados/variáveis (Parte 3)

limpar_coluna_binaria <- function(base, nome_coluna) {
  vars <- base$variables

 # Converte para texto (remove restrição de fator)
  vars[[nome_coluna]] <- as.character(vars[[nome_coluna]])

 # Substitui "Sim" -> 1 e "Não" -> 0
  vars[[nome_coluna]][vars[[nome_coluna]] == "Sim"] <- 1
  vars[[nome_coluna]][vars[[nome_coluna]] == "Não"] <- 0

 # Converte para numérico
  vars[[nome_coluna]] <- as.numeric(vars[[nome_coluna]])

 # Filtra apenas 0 e 1
  vars <- vars[vars[[nome_coluna]] %in% c(0, 1), ]

 # Atualiza o objeto survey.design
  base$variables <- vars
  base <- subset(base, vars[[nome_coluna]] %in% c(0, 1))

  return(base)
}

Base = limpar_coluna_binaria(Base, "Doença")
  • 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

Estimações

Estimações: Total Amostral

## Total amostral:
#Variáveis Numéricas:

total_renda = svytotal(x=~Rendimento_bruto, design=Base, na.rm=TRUE)
total_renda
                      total         SE
Rendimento_bruto 2.1323e+11 3604489769
confint(total_renda)
                        2.5 %       97.5 %
Rendimento_bruto 206163204561 220292544822
#Variáveis categóricas - Utilizando mais de uma coluna
svytotal(x=~Raça+Sexo, design=Base, na.rm=TRUE)
                 total       SE
RaçaBranca    91037722 686604.4
RaçaPreta     21786515 338021.3
RaçaAmarela    1674691 144663.7
RaçaParda     94086142 600647.8
RaçaIndígena    990265  74430.0
RaçaIgnorado     14272   5707.5
SexoHomem    100158109 240473.3
SexoMulher   109431498 240473.3
#Variáveis categóricas - Estimação por cruzamento
svytotal(x=~interaction(Sexo,Raça), design=Base, na.rm=TRUE)
                                            total       SE
interaction(Sexo, Raça)Homem.Branca    42682904.8 378164.8
interaction(Sexo, Raça)Mulher.Branca   48354816.8 404879.1
interaction(Sexo, Raça)Homem.Preta     10691164.0 186451.9
interaction(Sexo, Raça)Mulher.Preta    11095350.7 200050.6
interaction(Sexo, Raça)Homem.Amarela     754080.8  73347.6
interaction(Sexo, Raça)Mulher.Amarela    920610.6  81402.9
interaction(Sexo, Raça)Homem.Parda     45496265.5 337575.6
interaction(Sexo, Raça)Mulher.Parda    48589876.7 351852.0
interaction(Sexo, Raça)Homem.Indígena    523966.4  52076.8
interaction(Sexo, Raça)Mulher.Indígena   466298.3  33901.4
interaction(Sexo, Raça)Homem.Ignorado      9727.6   3936.5
interaction(Sexo, Raça)Mulher.Ignorado     4544.8   1978.7

Estimações: Médias, Proporções e Variância

##Médias e proporções

#svymean() 
svymean(x=~Rendimento_bruto, design=Base, na.rm=TRUE)
                   mean     SE
Rendimento_bruto 2224.8 35.465
svymean(x=~Sexo, design=Base, na.rm=TRUE) # frequência relativa (VC)
              mean     SE
SexoHomem  0.47788 0.0011
SexoMulher 0.52212 0.0011
##Variância

#svyvar() 
svyvar(x=~Rendimento_bruto, design=Base, na.rm=TRUE)
                 variance      SE
Rendimento_bruto 14048208 1720046

Estimações: Razão

##Razão

#svyratio()
razão <- svyratio(numerator=~(Motivo_Saude=="Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"),
denominator=~(Interrupção_Atividade=="Sim"), design=Base, na.rm=TRUE)
razão
Ratio estimator: svyratio.survey.design2(numerator = ~(Motivo_Saude == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"), 
    denominator = ~(Interrupção_Atividade == "Sim"), design = Base, 
    na.rm = TRUE)
Ratios=
                                                                                                      Interrupção_Atividade == "Sim"
Motivo_Saude == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"                      0.2096767
SEs=
                                                                                                      Interrupção_Atividade == "Sim"
Motivo_Saude == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"                    0.005052547
cv(object=razão)
                                                                                                      Interrupção_Atividade == "Sim"
Motivo_Saude == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"                     0.02409685
confint(object=razão, level = 0.95)
                                                                                                                                         2.5 %
Motivo_Saude == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"/Interrupção_Atividade == "Sim" 0.1997739
                                                                                                                                        97.5 %
Motivo_Saude == "Problemas respiratórios (Resfriado / gripe /sinusite/ asma / bronquite / pneumonia)"/Interrupção_Atividade == "Sim" 0.2195795

Estimações: Mediana e Quantis

##Mediana e Quantis

#svyquantile(c(.25,.5,.75)) # ci=TRUE/FALSE para mostrar ou não SE

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.

mediana = svyquantile(x=~Rendimento_bruto, design=Base, quantiles=0.5, ci=TRUE, na.rm=TRUE) ##MEDIANA
mediana
$Rendimento_bruto
    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=mediana)
Rendimento_bruto 
        25.50653 
cv(object=mediana)
Rendimento_bruto 
      0.01889373 
svyquantile(x=~Rendimento_bruto, design=Base, quantiles=c(.25,.5,.75), ci=FALSE, na.rm=TRUE) #Múltiplos quantis
$Rendimento_bruto
     0.25  0.5 0.75
[1,]  998 1350 2200

attr(,"hasci")
[1] FALSE
attr(,"class")
[1] "newsvyquantile"

Estimações: Por Domínio (Filtro)

# Estimação por dominio (filtro)
svytotal(x=~Rendimento_bruto, design=subset(Base, Sexo=="Mulher"), na.rm=TRUE)
                      total         SE
Rendimento_bruto 8.1577e+10 1440084928
svytotal(x=~Nível_Instrução, design=subset(Base, Sexo=="Homem" & Raça=="Parda" & Idade>30), na.rm=TRUE)
                                                       total     SE
Nível_InstruçãoSem instrução                         2573476  63547
Nível_InstruçãoFundamental incompleto ou equivalente 9167518 129978
Nível_InstruçãoFundamental completo ou equivalente   1943169  64405
Nível_InstruçãoMédio incompleto ou equivalente       1262245  48120
Nível_InstruçãoMédio completo ou equivalente         5777822 107139
Nível_InstruçãoSuperior incompleto ou equivalente     522290  30839
Nível_InstruçãoSuperior completo                     2011870  62089
svymean(x=~Nível_Instrução, design=subset(Base, Sexo=="Homem" & Raça=="Parda" & Idade>30), na.rm=TRUE)
                                                         mean     SE
Nível_InstruçãoSem instrução                         0.110647 0.0026
Nível_InstruçãoFundamental incompleto ou equivalente 0.394160 0.0045
Nível_InstruçãoFundamental completo ou equivalente   0.083547 0.0027
Nível_InstruçãoMédio incompleto ou equivalente       0.054271 0.0020
Nível_InstruçãoMédio completo ou equivalente         0.248419 0.0040
Nível_InstruçãoSuperior incompleto ou equivalente    0.022456 0.0013
Nível_InstruçãoSuperior completo                     0.086501 0.0026
##Coeficiente de variação
cv(svytotal(x=~Rendimento_bruto, design=subset(Base, Sexo=="Mulher"), na.rm=TRUE))
                 Rendimento_bruto
Rendimento_bruto       0.01765309

Estimações: Por Vários Domínios (“Group by”)

#Estimação para vários dominios ("Groupby/Sumarize")

svyby(formula=~Sexo, by=~Nível_Instrução, design=Base, FUN=svymean, na.rm=TRUE)
                                                            Nível_Instrução
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
                                      SexoHomem SexoMulher se.SexoHomem
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.SexoMulher
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.003870122

PLOTS

Plots: Histogramas

##HISTOGRAMAS


svyhist(formula=~as.numeric(HorasSemanaisTrabalho), design=Base, main="Histograma", xlab="Número de Horas Normalmente Trabalhadas")

#formula: variavel
#main,xlab:titulo, subtitulo

svyhist(formula=~as.numeric(HorasSemanaisTrabalho), design=Base, freq=TRUE, main="Histograma", xlab="Número de Horas Normalmente Trabalhadas")

Plots: Boxplots

##BOXBLOTS

svyboxplot(formula=HorasSemanaisTrabalho~1, design=Base, main="Boxplot do Número de Horas Normalmente Trabalhadas")

#~ ->variavel de agrupamento, caso queira sem, colocar 1
svyboxplot(formula=HorasSemanaisTrabalho~Sexo, design=Base, main="Boxplot do Número de Horas Normalmente Trabalhadas por Sexo")

#APENAS OUTLIERS MAIS EXTREMOS
svyboxplot(formula=HorasSemanaisTrabalho~Sexo, design=Base, all.outliers=TRUE, main="Boxplot do Número de Horas Normalmente Trabalhadas por Sexo")

#TODOS

Plots: Gráfico de Dispersão

##GRÁFICO DE DISPERSÃO

svyplot(formula=Rendimento_bruto~HorasSemanaisTrabalho, design=Base, style="transparent", xlab="Número de Horas Normalmente Trabalhadas", ylab="Rendimento Mensal Bruto Normalmente Recebido")
#style transparent

Modelagens

Modelagens: Testes de Hipóteses

##TESTES DE HIPOTESES

# Teste 1: Horas trabalhadas por sexo

pns_teste1 <- subset(Base, !is.na(HorasSemanaisTrabalho) & !is.na(Sexo)) #filtrando daados validos
teste1 <- svyttest(HorasSemanaisTrabalho ~ Sexo, design = pns_teste1)
teste1

    Design-based t-test

data:  HorasSemanaisTrabalho ~ Sexo
t = -39.962, df = 7451, p-value < 2.2e-16
alternative hypothesis: true difference in mean is not equal to 0
95 percent confidence interval:
 -5.686913 -5.155071
sample estimates:
difference in mean 
         -5.420992 
cat("\nInterpretacao: Diferenca media de", round(abs(teste1$estimate), 2), "horas\n")

Interpretacao: Diferenca media de 5.42 horas
cat("Conclusao: Ha diferenca significativa de horas trabalhadas entre sexos.\n\n")
Conclusao: Ha diferenca significativa de horas trabalhadas entre sexos.

Modelagens: Regressão Linear

#REGRESSÃO LINEAR

modeloLin <- svyglm(formula=Rendimento_bruto~Nível_Instrução+Raça+Idade, design=Base)
summary(object=modeloLin)

Call:
svyglm(formula = Rendimento_bruto ~ Nível_Instrução + Raça + 
    Idade, design = Base)

Survey design:
subset(base, vars[[nome_coluna]] %in% c(0, 1))

Coefficients:
                                                     Estimate Std. Error
(Intercept)                                          -658.631     93.557
Nível_InstruçãoFundamental incompleto ou equivalente  420.194     36.784
Nível_InstruçãoFundamental completo ou equivalente    856.624     44.061
Nível_InstruçãoMédio incompleto ou equivalente       1050.529     50.686
Nível_InstruçãoMédio completo ou equivalente         1311.015     55.812
Nível_InstruçãoSuperior incompleto ou equivalente    1790.372     70.420
Nível_InstruçãoSuperior completo                     4077.126    106.676
RaçaPreta                                            -670.738     46.133
RaçaAmarela                                          -117.818    214.938
RaçaParda                                            -608.849     41.372
RaçaIndígena                                         -746.442     97.356
RaçaIgnorado                                          363.825   1324.524
Idade                                                  39.837      1.872
                                                     t value Pr(>|t|)    
(Intercept)                                           -7.040 2.10e-12 ***
Nível_InstruçãoFundamental incompleto ou equivalente  11.423  < 2e-16 ***
Nível_InstruçãoFundamental completo ou equivalente    19.442  < 2e-16 ***
Nível_InstruçãoMédio incompleto ou equivalente        20.726  < 2e-16 ***
Nível_InstruçãoMédio completo ou equivalente          23.490  < 2e-16 ***
Nível_InstruçãoSuperior incompleto ou equivalente     25.424  < 2e-16 ***
Nível_InstruçãoSuperior completo                      38.220  < 2e-16 ***
RaçaPreta                                            -14.539  < 2e-16 ***
RaçaAmarela                                           -0.548    0.584    
RaçaParda                                            -14.716  < 2e-16 ***
RaçaIndígena                                          -7.667 1.98e-14 ***
RaçaIgnorado                                           0.275    0.784    
Idade                                                 21.285  < 2e-16 ***
---
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: 2
##Intervalos de confiança 95%

ic_linear<-confint(modeloLin)

coef_df <- data.frame(
 Variavel = names(coef(modeloLin)),
 Coeficiente = as.numeric(coef(modeloLin)),
 IC_Inf = ic_linear[, 1],
 IC_Sup = ic_linear[, 2]
)

coef_df
                                                                                                 Variavel
(Intercept)                                                                                   (Intercept)
Nível_InstruçãoFundamental incompleto ou equivalente Nível_InstruçãoFundamental incompleto ou equivalente
Nível_InstruçãoFundamental completo ou equivalente     Nível_InstruçãoFundamental completo ou equivalente
Nível_InstruçãoMédio incompleto ou equivalente             Nível_InstruçãoMédio incompleto ou equivalente
Nível_InstruçãoMédio completo ou equivalente                 Nível_InstruçãoMédio completo ou equivalente
Nível_InstruçãoSuperior incompleto ou equivalente       Nível_InstruçãoSuperior incompleto ou equivalente
Nível_InstruçãoSuperior completo                                         Nível_InstruçãoSuperior completo
RaçaPreta                                                                                       RaçaPreta
RaçaAmarela                                                                                   RaçaAmarela
RaçaParda                                                                                       RaçaParda
RaçaIndígena                                                                                 RaçaIndígena
RaçaIgnorado                                                                                 RaçaIgnorado
Idade                                                                                               Idade
                                                     Coeficiente      IC_Inf
(Intercept)                                           -658.63051  -842.02949
Nível_InstruçãoFundamental incompleto ou equivalente   420.19439   348.08821
Nível_InstruçãoFundamental completo ou equivalente     856.62422   770.25288
Nível_InstruçãoMédio incompleto ou equivalente        1050.52940   951.17111
Nível_InstruçãoMédio completo ou equivalente          1311.01518  1201.60830
Nível_InstruçãoSuperior incompleto ou equivalente     1790.37207  1652.32816
Nível_InstruçãoSuperior completo                      4077.12608  3868.01043
RaçaPreta                                             -670.73813  -761.17246
RaçaAmarela                                           -117.81785  -539.15746
RaçaParda                                             -608.84930  -689.95080
RaçaIndígena                                          -746.44221  -937.28735
RaçaIgnorado                                           363.82458 -2232.61649
Idade                                                   39.83676    36.16794
                                                         IC_Sup
(Intercept)                                          -475.23153
Nível_InstruçãoFundamental incompleto ou equivalente  492.30056
Nível_InstruçãoFundamental completo ou equivalente    942.99555
Nível_InstruçãoMédio incompleto ou equivalente       1149.88770
Nível_InstruçãoMédio completo ou equivalente         1420.42206
Nível_InstruçãoSuperior incompleto ou equivalente    1928.41597
Nível_InstruçãoSuperior completo                     4286.24173
RaçaPreta                                            -580.30380
RaçaAmarela                                           303.52177
RaçaParda                                            -527.74780
RaçaIndígena                                         -555.59706
RaçaIgnorado                                         2960.26564
Idade                                                  43.50557

Modelagens: Regressão Logística (Parte 1)

#REGRESSÃO LOG

modeloLog <- svyglm(formula=Doença~Sexo+Raça+Idade, design=Base, family="binomial")
summary(object=modeloLog)

Call:
svyglm(formula = Doença ~ Sexo + Raça + Idade, design = Base, 
    family = "binomial")

Survey design:
subset(base, vars[[nome_coluna]] %in% c(0, 1))

Coefficients:
               Estimate Std. Error  t value Pr(>|t|)    
(Intercept)  -3.1782216  0.0290104 -109.554  < 2e-16 ***
SexoMulher    0.3486645  0.0158377   22.015  < 2e-16 ***
RaçaPreta     0.0013958  0.0277866    0.050    0.960    
RaçaAmarela   0.2138631  0.2078184    1.029    0.303    
RaçaParda    -0.0919759  0.0189782   -4.846 1.28e-06 ***
RaçaIndígena -0.0042507  0.1099412   -0.039    0.969    
RaçaIgnorado -0.0345987  0.6965349   -0.050    0.960    
Idade         0.0559756  0.0004777  117.176  < 2e-16 ***
---
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: 4
#Probabilidades:
#Intercept:exp(Estimate)/(1+exp(Estimate))
#Outros: 1 - exp(Estimate)
#Idade: Incremento por ano

Modelagens: Regressão Logística (Parte 2)

# Obter os coeficientes (log-odds) e seus intervalos de confiança
log_odds <- coef(modeloLog)
ic_log_odds <- confint(modeloLog)

# Exponenciar para obter os Odds Ratios e seus ICs, evitando valores negativos
resultado_or <- data.frame(
 Variavel = names(log_odds),
 OR = round(exp(log_odds), 3),
 IC_Inferior = round(exp(ic_log_odds[, 1]), 3),
 IC_Superior = round(exp(ic_log_odds[, 2]), 3)
)

print(resultado_or)
                 Variavel    OR IC_Inferior IC_Superior
(Intercept)   (Intercept) 0.042       0.039       0.044
SexoMulher     SexoMulher 1.417       1.374       1.462
RaçaPreta       RaçaPreta 1.001       0.948       1.057
RaçaAmarela   RaçaAmarela 1.238       0.824       1.861
RaçaParda       RaçaParda 0.912       0.879       0.947
RaçaIndígena RaçaIndígena 0.996       0.803       1.235
RaçaIgnorado RaçaIgnorado 0.966       0.247       3.784
Idade               Idade 1.058       1.057       1.059

Modelagens: Regressão Logística (Gráfico)

Pacote Convey

Pacote Convey: Gini

##Índice de Gini
##FAZER ESTE POR ULTIMO
dadosPNS <- convey_prep(design=Base)
giniBR <- svygini(formula=~Rendimento_bruto, design=dadosPNS, na.rm=TRUE)
giniBR
                    gini     SE
Rendimento_bruto 0.50802 0.0056
cv(object=giniBR)
                 Rendimento_bruto
Rendimento_bruto       0.01094457
giniUF <- svyby(formula=~Rendimento_bruto, by=~V0001, design=dadosPNS, FUN=svygini, na.rm=TRUE) #esse demora um pouco..
giniUF
                                  V0001 Rendimento_bruto          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.011241314

Pacote Convey: Curva de Lorenz

# Curva de Lorenz
curvaLorenz <- svylorenz(~E01602, design = dadosPNS, quantiles = seq(0, 1, 0.05), na.rm = TRUE)