Econometria com Complex Survey: Resolução exercício 7.6 Heeringa et al (2010)

Abstract

This is an undergrad student level instruction for class use. EXAMPLE 7.6 WITH COMPLEX SAMPLE ADJUSTMENT AND WEIGHTS USING SVYGLM.

Licença

This work is licensed under the Creative Commons Attribution-ShareAlike 4.0 International License. To view a copy of this license, visit http://creativecommons.org/licenses/by-sa/4.0/ or send a letter to Creative Commons, PO Box 1866, Mountain View, CA 94042, USA.

License: CC BY-SA 4.0

Citação

Sugestão de citação: FIGUEIREDO, Adriano Marcos Rodrigues. Econometria com Complex Survey: Resolução exercício 7.6 Heeringa et al (2010). Campo Grande-MS,Brasil: RStudio/Rpubs, 2020 (updated links in 2025). Disponível em http://rpubs.com/amrofi/Econometrics_survey_heeringa76.

1 Introdução

Fonte: HEERINGA et al (2010, p.224-226).

Seja o dataset 2005–2006 NHANES data (<chapter_exercises_nhanes0506.rdata> em https://smponline.isr.umich.edu/asda/R%20data%20chapter_exercises.zip), pede-se estimar o modelo linear de regressão assumindo: variável dependente é a medida de pressão sanguínea (BPXSY1), contra os regressores BMI (índice de massa corporal - BMXBMI), idade do respondente (RIDAGEYR), Gênero (RIAGENDR: 1 = Male, 2 = Female), Raça (RIDRETH1: 1 = Mexican, 2 = Other Hispanic, 3 = White, 4 = Black, 5 = Other), e índice de pobreza (INDFMPIR). Certifique de usar o plano amostral complexo. Lembrar que gênero e raça são variáveis qualitativas (usar factor). Qual é a evidência de que os preditores sócio-demográficos (idade, gênero, raça e pobreza) melhoram o ajuste do modelo comparativamente a um modelo apenas com o preditor BMI? SUGESTÃO: olhar o livro do Heeringa et al (2010) e o arquivo de análise do capítulo 7 em https://smponline.isr.umich.edu/asda/ASDA3_Datasets/Analysis%20Examples%20and%20Chapter%20Exercises%20Data%20Sets%20ASDA3%20%20Stata%20SAS%20and%20R%20Formats.zip

Portanto, os passos para desenvolvimento são:

1.    acessar livro    
2.    baixar dados    
3.    processar variáveis    
4.    rodar modelo    
5.    interpretar saída.    

2 Puxar dados

Os dados, conforme enunciado, estão em (<chapter_exercises_nhanes0506.rdata>) e foram previamente baixados do arquivo zip em (http://www.isr.umich.edu/src/smp/asda/R%20data%20chapter_exercises.zip) e colocados no diretório de trabalho do projeto.

# chamar os dados
load("chapter_exercises_nhanes0506.rdata")
nhanesdata <- chapter_exercises_nhanes0506
rm(chapter_exercises_nhanes0506)
print(nhanesdata)

   seqn riagendr ridageyr ridreth1 indfmpir wtint2yr wtmec2yr sdmvpsu sdmvstra
1 31127        1        0        3     0.75 6434.950 6571.396       2       44
2 31128        2       11        4     0.77 9081.701 8987.042       1       52
3 31129        1       15        4     2.71 5316.895 5586.719       1       51
  bpxpls bpxpuls bpxsy1 bpxdi1 bpxsy2 bpxdi2 bpxsy3 bpxdi3 bpxsy4 bpxdi4 bmxbmi
1     NA       1     NA     NA     NA     NA     NA     NA     NA     NA     NA
2     84       1    100     62     NA     NA     NA     NA     NA     NA  17.45
3     96       1    104     76    106     72    108     70     NA     NA  26.53
  lbxtc lbdhdd irregular edcat marcat bpxdi1_1 age18p age51
1    NA     NA         0    NA     NA       NA      0     0
2   129     55         0     1     NA       62      0     0
3   170     46         0     1      3       76      0     0
 [ reached 'max' / getOption("max.print") -- omitted 10345 rows ]

names(nhanesdata)

 [1] "seqn"      "riagendr"  "ridageyr"  "ridreth1"  "indfmpir"  "wtint2yr" 
 [7] "wtmec2yr"  "sdmvpsu"   "sdmvstra"  "bpxpls"    "bpxpuls"   "bpxsy1"   
[13] "bpxdi1"    "bpxsy2"    "bpxdi2"    "bpxsy3"    "bpxdi3"    "bpxsy4"   
[19] "bpxdi4"    "bmxbmi"    "lbxtc"     "lbdhdd"    "irregular" "edcat"    
[25] "marcat"    "bpxdi1_1"  "age18p"    "age51"    

# 'sdmvpsu' (which contains masked versions, or approximations, of the true
# primary sampling unit codes for each respondent for the purposes of variance
# estimation; see Section 4.3.1) and 'sdmvstra' (which contains the estimation
# purposes)
summary(nhanesdata)

      seqn          riagendr        ridageyr     ridreth1        indfmpir    
 Min.   :31127   Min.   :1.000   Min.   : 0   Min.   :1.000   Min.   :0.000  
 1st Qu.:33714   1st Qu.:1.000   1st Qu.: 9   1st Qu.:1.000   1st Qu.:1.000  
 Median :36301   Median :2.000   Median :19   Median :3.000   Median :1.950  
    wtint2yr         wtmec2yr         sdmvpsu         sdmvstra    
 Min.   :  1225   Min.   :     0   Min.   :1.000   Min.   :44.00  
 1st Qu.:  7247   1st Qu.:  7019   1st Qu.:1.000   1st Qu.:47.00  
 Median : 18971   Median : 18037   Median :2.000   Median :51.00  
     bpxpls          bpxpuls          bpxsy1         bpxdi1      
 Min.   : 40.00   Min.   :1.000   Min.   : 74    Min.   :  0.00  
 1st Qu.: 66.00   1st Qu.:1.000   1st Qu.:106    1st Qu.: 56.00  
 Median : 74.00   Median :1.000   Median :116    Median : 66.00  
     bpxsy2          bpxdi2           bpxsy3          bpxdi3      
 Min.   : 72.0   Min.   :  0.00   Min.   : 78.0   Min.   :  0.00  
 1st Qu.:106.0   1st Qu.: 58.00   1st Qu.:106.0   1st Qu.: 58.00  
 Median :114.0   Median : 66.00   Median :114.0   Median : 66.00  
     bpxsy4          bpxdi4           bmxbmi           lbxtc      
 Min.   : 76.0   Min.   :  0.00   Min.   : 11.74   Min.   : 78.0  
 1st Qu.:108.0   1st Qu.: 56.00   1st Qu.: 19.29   1st Qu.:155.0  
 Median :116.0   Median : 66.00   Median : 24.31   Median :179.0  
     lbdhdd         irregular           edcat           marcat     
 Min.   : 15.00   Min.   :0.00000   Min.   :1.000   Min.   :1.000  
 1st Qu.: 43.00   1st Qu.:0.00000   1st Qu.:1.000   1st Qu.:1.000  
 Median : 52.00   Median :0.00000   Median :1.000   Median :2.000  
    bpxdi1_1          age18p           age51       
 Min.   :  4.00   Min.   :0.0000   Min.   :0.0000  
 1st Qu.: 56.00   1st Qu.:0.0000   1st Qu.:0.0000  
 Median : 66.00   Median :1.0000   Median :0.0000  
 [ reached 'max' / getOption("max.print") -- omitted 4 rows ]

3 Organizando os dados

As variáveis, conforme sugestão do enunciado do exercício e descrito no arquivo de análise do capítulo 7 de Heeringa et al (2010) em http://www.isr.umich.edu/src/smp/asda/mysite/R%20analysis%20examples/Chapter%207%20ANALYSIS%20EXAMPLES%20REPLICATION%20FALL%202011%20R.pdf é fazer uma organização para retirar os dados ausentes, NAs, e para especificar os rótulos de gênero (Masculino e Feminino) e criar as variáveis de fatores para raça (Mexicano, outros hispânicos, brancos, pretos e outros), estado civil (casado, casado anteriormente, nunca casado), educação etc

nhanesdata$bpxdi1_1 <- nhanesdata$bpxdi1
# especificar bpxdi1 zeros para missing (NA)
nhanesdata$bpxdi1_1[nhanesdata$bpxdi1_1 == 0] <- NA
nhanesdata$genderc <- factor(nhanesdata$riagendr, levels = 1:2, labels = c("M", "F"))
# create factor variables
nhanesdata$racec <- factor(nhanesdata$ridreth1, levels = 1:5, labels = c("Mexican",
    "Other Hispanic", "White", "Black", "Other"))
nhanesdata$marcatc <- factor(nhanesdata$marcat, levels = 1:3, labels = c("Married",
    "Previously Married", "Never Married"))
nhanesdata$edcatc <- factor(nhanesdata$edcat, levels = 1:4, labels = c("0-11", "12",
    "13-15", "16+"))
names(nhanesdata)

 [1] "seqn"      "riagendr"  "ridageyr"  "ridreth1"  "indfmpir"  "wtint2yr" 
 [7] "wtmec2yr"  "sdmvpsu"   "sdmvstra"  "bpxpls"    "bpxpuls"   "bpxsy1"   
[13] "bpxdi1"    "bpxsy2"    "bpxdi2"    "bpxsy3"    "bpxdi3"    "bpxsy4"   
[19] "bpxdi4"    "bmxbmi"    "lbxtc"     "lbdhdd"    "irregular" "edcat"    
[25] "marcat"    "bpxdi1_1"  "age18p"    "age51"     "genderc"   "racec"    
[31] "marcatc"   "edcatc"   

4 Plano amostral

Aqui deve-se cuidar para especificar o strata, o id (indicador do psu), pesos e dataset.

library(survey)
nhanessvy2 <- svydesign(strata = ~sdmvstra, id = ~sdmvpsu, weights = ~wtmec2yr, data = nhanesdata,
    nest = T)
# subnhanes <- subset(nhanessvy2 , age18p == 1)

5 Modelagem: “EXAMPLE 7.6 WITH COMPLEX SAMPLE ADJUSTMENT AND WEIGHTS USING SVYGLM”

Portanto, agora vou rodar o modelo. Sejam as variáveis de análise:

# medida de pressão sanguínea (BPXSY1), variável dependente

contra os regressores 

# BMI (índice de massa corporal - BMXBMI),
# idade do respondente  (RIDAGEYR), 
# Gênero (RIAGENDR: 1 = Male, 2 = Female), 
# Raça (RIDRETH1: 1 = Mexican, 2 = Other Hispanic, 3 = White, 4 = Black, 5 = Other), 
# Índice de pobreza (INDFMPIR).

ex76_svyglm <- svyglm(bpxsy1 ~ bmxbmi + ridageyr + factor(riagendr) + factor(ridreth1) +
    indfmpir, design = nhanessvy2)
summary(ex76_svyglm)

Call:
svyglm(formula = bpxsy1 ~ bmxbmi + ridageyr + factor(riagendr) + 
    factor(ridreth1) + indfmpir, design = nhanessvy2)

Survey design:
svydesign(strata = ~sdmvstra, id = ~sdmvpsu, weights = ~wtmec2yr, 
    data = nhanesdata, nest = T)

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       95.73268    1.20124  79.695 1.29e-11 ***
bmxbmi             0.39411    0.05326   7.399 0.000149 ***
ridageyr           0.44837    0.01368  32.783 6.36e-09 ***
factor(riagendr)2 -3.45526    0.48273  -7.158 0.000184 ***
factor(ridreth1)2 -0.69765    1.63213  -0.427 0.681895    
factor(ridreth1)3  0.53653    0.56510   0.949 0.374013    
factor(ridreth1)4  3.27957    0.55792   5.878 0.000613 ***
factor(ridreth1)5 -0.35906    1.12646  -0.319 0.759217    
indfmpir          -0.92678    0.15280  -6.065 0.000508 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 279.5849)

Number of Fisher Scoring iterations: 2

plot(ex76_svyglm)

6 Opção

Este procedimento utiliza as variáveis criadas anteriormente e faz o mesmo modelo. É possível ver que resulta no mesmo modelo anterior.

ex76_svyglm2 <- svyglm(bpxsy1 ~ bmxbmi + ridageyr + genderc + racec + indfmpir, design = nhanessvy2)
summary(ex76_svyglm2)

Call:
svyglm(formula = bpxsy1 ~ bmxbmi + ridageyr + genderc + racec + 
    indfmpir, design = nhanessvy2)

Survey design:
svydesign(strata = ~sdmvstra, id = ~sdmvpsu, weights = ~wtmec2yr, 
    data = nhanesdata, nest = T)

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)         95.73268    1.20124  79.695 1.29e-11 ***
bmxbmi               0.39411    0.05326   7.399 0.000149 ***
ridageyr             0.44837    0.01368  32.783 6.36e-09 ***
gendercF            -3.45526    0.48273  -7.158 0.000184 ***
racecOther Hispanic -0.69765    1.63213  -0.427 0.681895    
racecWhite           0.53653    0.56510   0.949 0.374013    
racecBlack           3.27957    0.55792   5.878 0.000613 ***
racecOther          -0.35906    1.12646  -0.319 0.759217    
indfmpir            -0.92678    0.15280  -6.065 0.000508 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for gaussian family taken to be 279.5849)

Number of Fisher Scoring iterations: 2

plot(ex76_svyglm2)

7 Resultados pelo jtools

O pacote jtools é uma opção para visualizar alguns resultados e pode ser interessante para alguns usuários, de acordo com o que deseja observar nos resultados, por exemplo, \(R^2\), \(R^2_ajustado\), intervalo de confiança, correlação parcial, VIFs.

ex76_svyglm2 <- svyglm(bpxsy1 ~ bmxbmi + ridageyr + genderc + racec + indfmpir, design = nhanessvy2)
jtools::summ(ex76_svyglm2, vifs = TRUE, part.corr = TRUE, confint = TRUE)

Observations	6297
Dependent variable	bpxsy1
Type	Survey-weighted linear regression

R²	0.31
Adj. R²	-0.81

	Est.	2.5%	97.5%	t val.	p	VIF
(Intercept)	95.73	92.89	98.57	79.70	0.00	NA
bmxbmi	0.39	0.27	0.52	7.40	0.00	1.71
ridageyr	0.45	0.42	0.48	32.78	0.00	1.68
gendercF	-3.46	-4.60	-2.31	-7.16	0.00	1.40
racecOther Hispanic	-0.70	-4.56	3.16	-0.43	0.68	2.55
racecWhite	0.54	-0.80	1.87	0.95	0.37	2.55
racecBlack	3.28	1.96	4.60	5.88	0.00	2.55
racecOther	-0.36	-3.02	2.30	-0.32	0.76	2.55
indfmpir	-0.93	-1.29	-0.57	-6.07	0.00	1.37
Standard errors: Robust

Agora com erros robustos, observe a informação do pacote survey - Robust standard errors are reported by default in the survey package:

# ex76_svyglm2 <- svyglm(bpxsy1 ~ bmxbmi+ridageyr+genderc + racec + indfmpir,
# design=nhanessvy2)
library(jtools)
jtools::summ(ex76_svyglm2, vifs = TRUE, robust = "HC3", part.corr = TRUE, confint = TRUE)

Observations	6297
Dependent variable	bpxsy1
Type	Survey-weighted linear regression

R²	0.31
Adj. R²	-0.81

	Est.	2.5%	97.5%	t val.	p	VIF
(Intercept)	95.73	92.89	98.57	79.70	0.00	NA
bmxbmi	0.39	0.27	0.52	7.40	0.00	1.71
ridageyr	0.45	0.42	0.48	32.78	0.00	1.68
gendercF	-3.46	-4.60	-2.31	-7.16	0.00	1.40
racecOther Hispanic	-0.70	-4.56	3.16	-0.43	0.68	2.55
racecWhite	0.54	-0.80	1.87	0.95	0.37	2.55
racecBlack	3.28	1.96	4.60	5.88	0.00	2.55
racecOther	-0.36	-3.02	2.30	-0.32	0.76	2.55
indfmpir	-0.93	-1.29	-0.57	-6.07	0.00	1.37
Standard errors: Robust

Posso normalizar os coeficientes para 1 desvio-padrão a fim de melhor visualizar os efeitos de cada variável sobre a medida de pressão sanguínea (BPXSY1), variável dependente. Depois plotar os resultados padronizados. Isso no caso de mudar para scale = T. Abaixo, não farei essa padronização.

jtools::summ(ex76_svyglm2, vifs = TRUE, robust = "HC3", scale = F, part.corr = TRUE,
    confint = TRUE, transform.response = F)

Observations	6297
Dependent variable	bpxsy1
Type	Survey-weighted linear regression

R²	0.31
Adj. R²	-0.81

	Est.	2.5%	97.5%	t val.	p	VIF
(Intercept)	95.73	92.89	98.57	79.70	0.00	NA
bmxbmi	0.39	0.27	0.52	7.40	0.00	1.71
ridageyr	0.45	0.42	0.48	32.78	0.00	1.68
gendercF	-3.46	-4.60	-2.31	-7.16	0.00	1.40
racecOther Hispanic	-0.70	-4.56	3.16	-0.43	0.68	2.55
racecWhite	0.54	-0.80	1.87	0.95	0.37	2.55
racecBlack	3.28	1.96	4.60	5.88	0.00	2.55
racecOther	-0.36	-3.02	2.30	-0.32	0.76	2.55
indfmpir	-0.93	-1.29	-0.57	-6.07	0.00	1.37
Standard errors: Robust

coef_names <- c(`BMI (índice de massa corporal)` = "bmxbmi", `idade do respondente` = "ridageyr",
    `Gênero feminino` = "gendercF", `Raça (outro hispânico)` = "racecOther Hispanic",
    `Raça (Branca)` = "racecWhite", `Raça (Preta)` = "racecBlack", `Raça (Outra)` = "racecOther",
    `Índice de pobreza` = "indfmpir", Constant = "(Intercept)")
coef_names <- coef_names[1:8]  # retirei intercept do plot
plot_summs(ex76_svyglm2, coefs = coef_names)

Podemos observar o efeito preditor de uma variável com uso da função jtools::effect_plot, por exemplo, para o BMI (bmxbmi). Neste caso ele plota os dados brutos das variáveis:

effect_plot(ex76_svyglm2, pred = bmxbmi, interval = TRUE, plot.points = TRUE)