Aula 2 - Problemas de Especificação
Isis Michelan
Econometria II
Temos visto ao longo das aulas de econometria que rodar um modelo de regressão exige bem mais que jogar Y ~ X no R. Os problemas que lidaremos hoje se referem ao refinamento do econometrista quanto ao conhecimento do seu banco de dados e do modelo que busca realizar inferência (suas implicações econômicas).
São alguns dos problemas de especificação que iremos abordar aqui brevemente:
- Omissão de uma ou mais variáveis relevantes.
- Adoção da forma funcional errada.
- Erros de medida.
Vamos começar abrindo o banco de dados (disponibilizado no classroom):
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(dplyr)
pnad <- read_csv("pnad_pv.csv") %>%
select(-1)## New names:
## Rows: 376821 Columns: 9
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (3): uf, sexo, cor dbl (6): ...1, ano, trim, idade, escolaridade, salario
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`glimpse(pnad)## Rows: 376,821
## Columns: 8
## $ ano <dbl> 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 2022, 202…
## $ trim <dbl> 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 4, 1, 1, …
## $ uf <chr> "Rondônia", "Rondônia", "Rondônia", "Rondônia", "Rondônia…
## $ sexo <chr> "Homem", "Mulher", "Mulher", "Mulher", "Mulher", "Homem",…
## $ idade <dbl> 30, 29, 49, 24, 2, 34, 72, 54, 33, 37, 10, 32, 45, 58, 54…
## $ cor <chr> "Parda", "Parda", "Parda", "Parda", "Parda", "Branca", "B…
## $ escolaridade <dbl> 12, 12, 16, 12, NA, 16, 16, 16, 8, 12, 4, 12, 12, 5, 12, …
## $ salario <dbl> 2500, 9999, 3500, 9999, 9999, 9999, 9999, 3500, 9999, 250…Como vimos nas últimas aulas, precisamos antes garantir que nossas variáveis categóricas estejam classificadas como “factor”:
pnad$sexo <- as.factor(pnad$sexo)
pnad$cor <- as.factor(pnad$cor)
summary(pnad)## ano trim uf sexo
## Min. :2022 Min. :1.000 Length:376821 Homem :183018
## 1st Qu.:2022 1st Qu.:1.000 Class :character Mulher:193803
## Median :2022 Median :2.000 Mode :character
## Mean :2022 Mean :2.481
## 3rd Qu.:2022 3rd Qu.:3.000
## Max. :2022 Max. :4.000
##
## idade cor escolaridade salario
## Min. : 0.00 Amarela : 2636 Min. : 0.000 Min. : 4
## 1st Qu.: 18.00 Branca :147984 1st Qu.: 5.000 1st Qu.: 1800
## Median : 36.00 Ignorado: 87 Median : 9.000 Median : 9999
## Mean : 37.01 Indígena: 2084 Mean : 8.425 Mean : 6796
## 3rd Qu.: 54.00 Parda :186809 3rd Qu.:12.000 3rd Qu.: 9999
## Max. :114.00 Preta : 37221 Max. :16.000 Max. :240000
## NA's :22569Veja que agora temos uma visualização das estatísticas de cada variável. Para a variável de cor, devido à baixa ocorrência de algumas categorias, podemos agrupá-las em duas para facilitar a comparação:
pnad <- pnad %>%
filter(cor != "Ignorado") %>%
mutate(cor = ifelse(cor %in% c("Indígena","Parda","Preta"),"Pretos ou Pardos",cor)) %>%
mutate(cor = ifelse(cor %in% c(1,2), "Brancos",cor))
#observando a nova classificação de cor
pnad$cor <- as.factor(pnad$cor)
summary(pnad)## ano trim uf sexo
## Min. :2022 Min. :1.000 Length:376734 Homem :182982
## 1st Qu.:2022 1st Qu.:1.000 Class :character Mulher:193752
## Median :2022 Median :2.000 Mode :character
## Mean :2022 Mean :2.481
## 3rd Qu.:2022 3rd Qu.:3.000
## Max. :2022 Max. :4.000
##
## idade cor escolaridade salario
## Min. : 0.00 Brancos :150620 Min. : 0.000 Min. : 4
## 1st Qu.: 18.00 Pretos ou Pardos:226114 1st Qu.: 5.000 1st Qu.: 1800
## Median : 36.00 Median : 9.000 Median : 9999
## Mean : 37.01 Mean : 8.424 Mean : 6796
## 3rd Qu.: 54.00 3rd Qu.:12.000 3rd Qu.: 9999
## Max. :114.00 Max. :16.000 Max. :240000
## NA's :22567#categoria de referência
levels(pnad$sexo) ## [1] "Homem" "Mulher"#para alterar a categoria de referência
#pnad$sexo <- relevel(pnad$sexo, ref = "Homem")
levels(pnad$cor)## [1] "Brancos" "Pretos ou Pardos"#retirando os NAs
anyNA(pnad)## [1] TRUEpnad <- na.omit(pnad)Vamos focar agora nosso olhar para a variável de salário, veja que sua média é bem alta e sua mediana um número relativamente estranho. Olhando para a distribuição da variável então temos:
#DISTRIBUIÇÃO AMOSTRAL
hist(pnad$salario)
grafico <- pnad %>%
ggplot(aes(x= salario)) +
geom_histogram() +
xlim(NA,30000)+
theme_light()+
ggtitle("Histograma")
grafico## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
O que acontece neste caso é que os valores nulos (NAs) de não respondentes da variável de salário foram reportados como “9999”, gerando uma distorção dos valores e devemos portanto retirá-los do banco de dados ou transformá-los em NAs. Os valores extremos de R$240.000,00, entretanto, não nos permite saber se este valor foi corretamente reportado ou não (erro de medida), dessa forma, não podemos removê-lo.
#ajustando os dados
pnad <- pnad %>%
filter(salario != 9999)
grafico <- pnad %>%
ggplot(aes(x= salario)) +
geom_histogram() +
xlim(NA,30000)+
theme_light()+
ggtitle("Histograma")
grafico## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
Agora, vamos dar uma olhada em um gráfico de “nuvem de pontos” do salário em relação a idade, que chamamos de “scatterplot”:
#forma funcional
grafico2 <- pnad %>%
filter(salario <= 10000) %>%
ggplot(aes(x= idade, y = salario))+
scale_x_continuous()+
geom_point()+
geom_smooth(method=lm , color="red", fill="#69b3a2", se=TRUE) +
theme_classic()
grafico2## `geom_smooth()` using formula = 'y ~ x'
Veja que por mais suja que a imagem esteja por termos muitas observações, há claramente um padrão não linear na distribuição de salários contra idade. Podemos ver pelo ajuste linear que se regredirmos linearmente Y~X teremos um problema de especificação; na equação abaixo veremos um dos jeitos de lidar com isso, no caso, adicionando a variável \(idade^2\).
Vamos rodar novamente a equação de Mincer, adicionando da nossa primeira estimação agora o termo de interação entre sexo e escolaridade. A ausência do termo pode viesar os parâmetros, uma vez que o efeito da educação sobre o salário é diferente para cada um dos sexos (por exemplo, se o coeficiente for positivo em relação à categoria mulher, significa que,em média as mulheres se beneficiam mais que os homens do aumento da educação em relação ao ganho salarial).
library(lmtest)## Carregando pacotes exigidos: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(sandwich)
#Criando a variável de idade^2
pnad <- pnad %>%
mutate(idade2 = idade^2)
#rodando a regressão
mqo <- lm(log(salario) ~ escolaridade + idade + idade2 + sexo + cor + sexo*escolaridade , data = pnad)
summary(mqo)##
## Call:
## lm(formula = log(salario) ~ escolaridade + idade + idade2 + sexo +
## cor + sexo * escolaridade, data = pnad)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.3931 -0.3939 0.0354 0.4497 5.1698
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.156e+00 1.766e-02 292.01 <2e-16 ***
## escolaridade 9.805e-02 6.057e-04 161.89 <2e-16 ***
## idade 6.123e-02 7.982e-04 76.71 <2e-16 ***
## idade2 -5.552e-04 9.299e-06 -59.71 <2e-16 ***
## sexoMulher -5.789e-01 1.132e-02 -51.15 <2e-16 ***
## corPretos ou Pardos -2.813e-01 4.048e-03 -69.48 <2e-16 ***
## escolaridade:sexoMulher 1.867e-02 9.632e-04 19.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7819 on 159318 degrees of freedom
## Multiple R-squared: 0.298, Adjusted R-squared: 0.298
## F-statistic: 1.127e+04 on 6 and 159318 DF, p-value: < 2.2e-16#erro padrão robusto pois já vimos que há heterogeneidade
coeftest(mqo, vcov = vcovHC(mqo, type="HC0"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 5.1555e+00 1.8805e-02 274.153 < 2.2e-16 ***
## escolaridade 9.8052e-02 6.8738e-04 142.645 < 2.2e-16 ***
## idade 6.1227e-02 8.8873e-04 68.892 < 2.2e-16 ***
## idade2 -5.5522e-04 1.0826e-05 -51.285 < 2.2e-16 ***
## sexoMulher -5.7895e-01 1.3463e-02 -43.002 < 2.2e-16 ***
## corPretos ou Pardos -2.8127e-01 4.0141e-03 -70.072 < 2.2e-16 ***
## escolaridade:sexoMulher 1.8674e-02 1.1159e-03 16.734 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Sabemos ainda que outras variáveis como UF (estado) e Zona (rural ou urbana), principalmente com dados brasileiros, tem grande importância nessa estimação e podem ser incluídos facilmente para melhorar o ajuste da regressão. Entretanto, existe porém uma categoria de variável não observada que é relevante para explicar o salário – a habilidade – e neste caso, com os dados que posuímos não é possível realizar uma estimação sem que haja viés. Para isso, muitos economistas utilizam a variável QI como proxy para a habilidade do trabalhador, apesar de sabermos hoje que está não é a melhor aproximação e existem outras formas mais modernas – que incluem outras áreas do conhecimento – para observar a habilidade dos trabalhadores no mercado.
Como saber se temos uma boa especificação?
Podemos realizar o teste RESET para compreender se há alguma má especificação do modelo, testando se há significância do nosso y estimado com relação à formas quadráticas e cúbicas. Seja \(\hat{y}\) nossa variável dependente estimada, temos que:
\[y = \beta_0+\beta_1*x_1+...+\beta_k*x_k+\delta_1*\hat{y}^2+\delta_2*\hat{y}^3+ erro\] Nossa hipótese nula será: \(H_0: \delta_1=\delta_2=0\). Isto é, caso o p-valor seja menor que os valores usuais de significância \((5\%,10\%)\) rejeitamos a hipótese nula e portanto teremos um problema de especificação.
#-----------------------------------------------------------------------------------------
#TESTE RESET
#-----------------------------------------------------------------------------------------
resettest(mqo)##
## RESET test
##
## data: mqo
## RESET = 1426.6, df1 = 2, df2 = 159316, p-value < 2.2e-16Para os modelos não aninhados (não hierárquicos, para simplificação), temos um teste que nos permite testar qual modelo de especificação se ajusta melhor aos dados. Vejamos o caso da variável de idade:
Modelo 1: \(log(salario) = \beta_0 + \beta_1*escolaridade + \beta_2*idade + \beta_3*idade2 + \beta_4*sexo + \beta_5*cor + \beta_6*sexo*escolaridade + u\)
Modelo 2: \(\beta_0 + \beta_1*escolaridade + \beta_2*log(idade) + \beta_3*sexo + \beta_4*cor + \beta_5*sexo*escolaridade + u\)
Modelo completo (não restrito): \(log(salario) = \beta_0 + \beta_1*escolaridade + \beta_2*idade + \beta_3*idade2 + \beta_4*log(idade) + \beta_5*sexo + \beta_6*cor + \beta_7*sexo:escolaridade + \varepsilon\)
O teste de Encompassing segue o modelo de Davidson & MacKinnon apresentado em aula; o teste tem por hipótese nula que os valores dos parâmetros do modelo irrestrito comparado ao modelo restrito devem ser iguais a zero, i.e, não existe ganho em adicionar outro tipo de especificação.
#-----------------------------------------------------------------------------------------
#TESTE PARA MODELOS NAO ANINHADOS
#-----------------------------------------------------------------------------------------
#MODELOS 1 E 2 CONTRA O MODELO ABRANGENTE
modelo1 <- lm(log(salario) ~ escolaridade + idade + idade2 + sexo + cor + sexo*escolaridade , data = pnad)
modelo2 <- lm(log(salario) ~ escolaridade + log(idade) + sexo + cor + sexo*escolaridade , data = pnad)
encomptest(modelo1,modelo2)## Encompassing test
##
## Model 1: log(salario) ~ escolaridade + idade + idade2 + sexo + cor + sexo *
## escolaridade
## Model 2: log(salario) ~ escolaridade + log(idade) + sexo + cor + sexo *
## escolaridade
## Model E: log(salario) ~ escolaridade + idade + idade2 + sexo + cor + log(idade) +
## escolaridade:sexo
## Res.Df Df F Pr(>F)
## M1 vs. ME 159317 -1 205.05 < 2.2e-16 ***
## M2 vs. ME 159317 -2 852.37 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Como podemos ver pelo resultado do teste (p-valor baixo tanto para o modelo 1 quanto para o modelo 2) o modelo irrestrito tem uma melhora no ajuste sobre os dois modelos, o que significa que nenhum dos dois é preferido e tem boa especificação.
Por hoje é só, espero que aproveitem este material. Recomendo leitura do capítulo 13 - Gujarati, Damodar e Porter, Dawn. Econometria básica - 5ª edição.
(Fonte: JEFFREY, M. Wooldridge, introductory econometrics—A modern approach. 2018.)