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

1 Base de dados

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

2 Regressão

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

3 Análise de Resíduos

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?

4 Teste para heterocedasticidade

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

5 Correção de White: Erros padrões robustos

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.

5.1 Fazendo manualmente

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

5.2 Usando o pacote

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 !