Antes de iniciar, vamos instalar os pacotes necessários (comando install.packages) e depois requisitar eles (comando library). Depois de instalado em seu computador, não precisa instalar mais, só requisitar.
Por isso, depois de instalar, coloque um # em frente ao install.packages para que não instale novamente. (Quando se coloca um # em frente ao comando, esse comando se torna um comentário)
#install.packages("dplyr") #pacote para limpar as bases
#install.packages("ggplot2") #pacote para plotar gráficos
#install.packages("lmtest") #pacote para fazer testes de regressão multipla
#install.packages("sandwich") #pacote para calcular erros padrões robustos
library(dplyr)
library(ggplot2)
library(lmtest)
library(sandwich)
options(scipen = 999) #Desligamos a notação cientifica
A base de dados utilizada será a PNAD 2015-Pesquisa Nacional por Amostra de Domicílios. Essa base de dados é uma base real do Brasil, e as variáveis contidas nelas são as seguintes:
X | Variável de contagem |
---|---|
id_dom | Variável identificadora do domicílio |
uf | Codigo do estado |
renda | Rendimento mensal em dinheiro no trabalho principal da semana de referência |
raca | Dummy indicando 1=Branco,0=Negro |
anos_estudo | Número de anos de estudo |
idade | Idade |
sexo | Dummy indicando 1=Homens, 0=Mulheres |
Foram feitos alguns ajustes para facilitar o exercício: foram removidos observações que não continham alguma dessas variáveis, e a variável de anos de estudo varia de 0 a 17 anos.
Vamos abrir a base de dados e nomea-la como “df” (abreviação de data frame). Lembre-se que para abrir a base de dados, ela deve estar na pasta do seu projeto. Caso não esteja na pasta, você pode colocar o endereço completo de onde esteja.
df<-read.csv("PNAD-2015-ajustada.csv")
Vamos ver como a base de dados está e algumas descritivas das variáveis.
names(df)
## [1] "X" "id_dom" "uf" "renda" "raca"
## [6] "anos_estudo" "idade" "sexo"
summary(df$anos_estudo)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.000 6.000 11.000 9.263 11.000 17.000
table(df$anos_estudo)
##
## 0 1 2 3 4 5 6 7 8 9 10 11
## 7813 1847 2865 4149 11313 6687 4523 5422 15638 4679 5307 49069
## 12 13 14 15 16 17
## 2954 3372 3564 14341 6166 2586
summary(df$renda)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 5 788 1100 1786 2000 200000
Vamos criar uma variável de idade ao quadrado
df$idade2<-(df$idade)^2
Vamos fazer uma regressão múltipla baseada na regressão de Mincer (1974).
MQO<-lm(renda~anos_estudo+raca+idade+idade2+sexo, data=df)
summary(MQO)
##
## Call:
## lm(formula = renda ~ anos_estudo + raca + idade + idade2 + sexo,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5045 -960 -315 391 196639
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3352.08771 55.74775 -60.13 <0.0000000000000002 ***
## anos_estudo 231.90138 1.59958 144.98 <0.0000000000000002 ***
## raca 481.14486 13.38146 35.96 <0.0000000000000002 ***
## idade 76.82222 2.68913 28.57 <0.0000000000000002 ***
## idade2 -0.40543 0.03224 -12.57 <0.0000000000000002 ***
## sexo 817.82548 13.33011 61.35 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2525 on 152289 degrees of freedom
## Multiple R-squared: 0.1671, Adjusted R-squared: 0.1671
## F-statistic: 6110 on 5 and 152289 DF, p-value: < 0.00000000000000022
Vamos plotar os resíduos para entender a intuição do problema da heterocedasticidade. Para isso, primeiro, vamos primeiro salvar o resíduo de cada observação em nosso data frame, e também salvar nosso valor previsto para a renda (\(\hat{y}\)).
df$residuos<-MQO$residuals
df$ychapeu<-MQO$fitted
Agora sim vamos plotar:
graf<-ggplot(data=df,
mapping=aes(y=residuos,x=ychapeu))+
geom_point(col="purple")
graf
Você acha que tem heterocedasticidade? Porque?
Vamos usar o teste de Breusch Pagan (versão LM do teste de heterocedasticidade).
\[H0 : Var(\mu \mid x_1,x_2,...,x_k)=\sigma^2 \text{ (Não há heterocedasticidade)}\] \[H1: Var(\mu \mid x_1,x_2,...,x_k) \neq \sigma^2 \text{ (Há heterocedasticidade)} \]
bptest(MQO)
##
## studentized Breusch-Pagan test
##
## data: MQO
## BP = 413.61, df = 5, p-value < 0.00000000000000022
P valor <0,05, rejeita-se a hipótese nula
Sabemos que um estimador válido da \(Var(\hat{\beta})\) na regressão mútipla é dado por:
\[\widehat{Var(\hat{\beta})}=\frac{\sum_{i=1}^{n}\hat{r_{ij}}^2\hat{\mu_i}^2}{SQR_j^2}\] (Wooldridge, 2011, p.250)
Em que \(r_{ij}\) representa o i-ésimo resíduo da regressão \(x_j\) sobre todas as outras variáveis independentes (regressão auxiliar).
E em que \(SQR_j= \sum_{i=1}^{n}r_i^2\) , ou seja, é a soma dos resíduos da regressão auxiliar.
Vamos calcular então o estimador da variância do \(\beta\) para a variável anos_estudo. Vamos fazer de duas formas: i) manualmente e ii) com o pacote. Depois vamos comparar os resultados.
Primeiro vamos calcular o \(r_{ij}\) em que \(j=\) anos_estudo. Vamos fazer a regressão auxiliar, e depois, salvar os resíduos.
Reg_aux=lm(anos_estudo~raca+idade+idade2+sexo,data=df)
df$residuos_aux<-Reg_aux$residuals
Vamos calcular o \(SQR_j= \sum_{i=1}^{n}r_i^2\)
SQR_j<-sum((df$residuos_aux)^2)
print(SQR_j)
## [1] 2491088
Vamos estimar o erro padrão robusto do estimador do coeficiente anos_estudo:
ep_rob_anos_estudo=sqrt(sum((df$residuos_aux)^2*(df$residuos)^2)/SQR_j^2)
print(ep_rob_anos_estudo)
## [1] 2.646141
Usando o pacote, calculamos todos os erros padrões robustos de uma vez só. Na forma manual, fizemos somente para anos_estudo (mas você pode fazer para as outras variáveis)
coeftest(MQO, vcov = vcovHC(MQO, type="HC0"))
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3352.087705 72.412373 -46.2916 < 0.00000000000000022 ***
## anos_estudo 231.901376 2.646141 87.6376 < 0.00000000000000022 ***
## raca 481.144862 12.779603 37.6494 < 0.00000000000000022 ***
## idade 76.822222 4.637187 16.5666 < 0.00000000000000022 ***
## idade2 -0.405426 0.063019 -6.4334 0.0000000001252 ***
## sexo 817.825479 13.982571 58.4889 < 0.00000000000000022 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Compare os resultados da forma manual e com o pacote !