Econometria: heteroscedasticidade: salários de advogados - Wooldridge(exemplo 7.8)

Abstract

This is an undergrad student level exercise for class use.

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: heteroscedasticidade: salários de advogados - Wooldridge(exemplo 7.8). Campo Grande-MS, Brasil: RStudio/Rpubs, 2020. Disponível em https://rpubs.com/amrofi/wooldridge_ex7_8_heterosc.

1 Introdução - Primeiros passos

Wooldridge Source: Collected by Kelly Barnett, an MSU economics student, for use in a term project. The data come from two sources: The Official Guide to U.S. Law Schools, 1986, Law School Admission Services, and The Gourman Report: A Ranking of Graduate and Professional Programs in American and International Universities, 1995, Washington, D.C. Data loads lazily.

O dataset do Wooldridge contém dados sobre a mediana dos salarios iniciais dos formados em Direito. Uma variável explicativa é a classificação da faculdade cursada. O grupo base é das faculdades com rank maior que 100.

O script foi adaptado de https://www.r-econometrics.com/reproduction/wooldridge/wooldridge07/.

No presente post faremos a análise específica da heteroscedasticidade complementando o post anterior em https://rpubs.com/amrofi/wooldridge_ex7_8 e https://rpubs.com/amrofi/wooldridge_ex7_8_spec, como também ao video da discussão de especificação em https://youtu.be/J7NK23J54bI.

Temos um data.frame com 156 observações para 21 variáveis:

rank: law school ranking    
salary: median starting salary    
cost: law school cost    
LSAT: median LSAT score    
GPA: median college GPA    
libvol: no. volumes in lib., 1000s    
faculty: no. of faculty    
age: age of law sch., years    
clsize: size of entering class    
north: =1 if law sch in north    
south: =1 if law sch in south    
east: =1 if law sch in east    
west: =1 if law sch in west    
lsalary: log(salary)    
studfac: student-faculty ratio    
top10: =1 if ranked in top 10    
r11_25: =1 if ranked 11-25    
r26_40: =1 if ranked 26-40    
r41_60: =1 if ranked 41-60    
llibvol: log(libvol)    
lcost: log(cost)    

2 Modelagem

A parte inicial reproduz o post de https://rpubs.com/amrofi/wooldridge_ex7_8_spec.

2.1 Carregar os dados

> library(wooldridge)
> data("lawsch85")

Vamos gerar dummies (não era necessário, pois já estavam no dataset) para os rankings das faculdades. Criamos abaixo para ver como ficam, pois o dataset contém uma variável rank para os salarios (salary). O leitor pode investigar a forma como criamos dummies usando a variável rank e critérios de classes.

> lawsch85$top10 <- as.numeric(lawsch85$rank <= 10)
> lawsch85$r11.25 <- as.numeric(lawsch85$rank >= 11 & lawsch85$rank <= 25)
> lawsch85$r26.40 <- as.numeric(lawsch85$rank >= 26 & lawsch85$rank <= 40)
> lawsch85$r41.60 <- as.numeric(lawsch85$rank >= 41 & lawsch85$rank <= 60)
> lawsch85$r61.100 <- as.numeric(lawsch85$rank >= 61 & lawsch85$rank <= 100)

As dummies foram criadas e já incorporadas ao dataset lawsch85. Podemos conferir, olhando pelo pacote knitr e função kable.

> library(knitr)
> kable(lawsch85[1:15,c(16:19,22:25)],caption="Dummies do dataset: recorte de 15 linhas")

Dummies do dataset: recorte de 15 linhas

top10	r11_25	r26_40	r41_60	r11.25	r26.40	r41.60	r61.100
0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0
0	0	1	0	0	1	0	0
0	0	0	1	0	0	1	0
0	0	0	0	0	0	0	1
0	0	0	0	0	0	0	1
0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	1
0	0	0	1	0	0	1	0
0	1	0	0	1	0	0	0
0	0	0	0	0	0	0	1
1	0	0	0	0	0	0	0
0	0	1	0	0	1	0	0

Com essas variáveis, criamos um modelo inicial (mod1) do tipo log-log:

\[ \operatorname{log(salary)} = \alpha + \beta_{1}(\operatorname{log(rank)}) + \beta_{2}(\operatorname{log(LSAT)}) + \beta_{3}(\operatorname{log(GPA)}) + \beta_{4}(\operatorname{log(libvol)}) + \beta_{5}(\operatorname{log(cost)}) + \epsilon \] Depois um (mod2) linear:

\[ \operatorname{(salary)} = \alpha + \beta_{1}(\operatorname{rank}) + \beta_{2}(\operatorname{LSAT}) + \beta_{3}(\operatorname{GPA}) + \beta_{4}(\operatorname{libvol}) + \beta_{5}(\operatorname{cost}) + \epsilon \] e um modelo log-log (mod3) com as dummies:

\[ \begin{array}{l} \operatorname{log(salary)} = \alpha + \beta_{1}(\operatorname{top10}) + \beta_{2}(\operatorname{r11\_25}) + \beta_{3}(\operatorname{r26\_40}) +\\ \beta_{4}(\operatorname{r41\_60}) + \beta_{5}(\operatorname{r61.100}) + \beta_{6}(\operatorname{log(LSAT)}) + \beta_{7}(\operatorname{log(GPA)}) +\\ \beta_{8}(\operatorname{log(libvol)}) + \beta_{9}(\operatorname{log(cost)}) + \epsilon \end{array} \]

2.2 Primeiros resultados

Portanto, os modelos são estimados e os resultados estão abaixo.

> mod1 <- lm(log(salary) ~ log(rank) + log(LSAT) + log(GPA) 
+            + log(libvol) + log(cost), data = lawsch85)
> mod2 <- lm((salary) ~ (rank) + (LSAT) + (GPA) 
+            + (libvol) + (cost), data = lawsch85)
> mod3 <- lm(log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + r61.100 
+            + log(LSAT) + log(GPA) + log(libvol) + log(cost), data = lawsch85)
> mod1$AIC <- AIC(mod1)
> mod2$AIC <- AIC(mod2)
> mod3$AIC <- AIC(mod3)
> mod1$BIC <- BIC(mod1)
> mod2$BIC <- BIC(mod2)
> mod3$BIC <- BIC(mod3)
> library(stargazer)
> star.1 <- stargazer(mod1, mod2, mod3,
+                     title="Título: Resultados das Regressões",
+                     align=TRUE,
+                     type = "text",
+                     keep.stat=c("aic","bic","rsq", "adj.rsq","n")
+ )

Título: Resultados das Regressões
=========================================================
                             Dependent variable:         
                    -------------------------------------
                    log(salary)   (salary)    log(salary)
                        (1)          (2)          (3)    
---------------------------------------------------------
log(rank)            -0.227***                           
                      (0.019)                            
                                                         
top10                                          0.700***  
                                                (0.053)  
                                                         
r11_25                                         0.593***  
                                                (0.039)  
                                                         
r26_40                                         0.374***  
                                                (0.034)  
                                                         
r41_60                                         0.262***  
                                                (0.028)  
                                                         
r61.100                                        0.131***  
                                                (0.021)  
                                                         
log(LSAT)              0.726                    0.885*   
                      (0.559)                   (0.472)  
                                                         
log(GPA)               0.191                     0.057   
                      (0.282)                   (0.239)  
                                                         
log(libvol)            0.017                     0.036   
                      (0.033)                   (0.026)  
                                                         
log(cost)              0.008                     0.001   
                      (0.029)                   (0.025)  
                                                         
rank                             -115.815***             
                                  (16.153)               
                                                         
LSAT                               155.286               
                                  (192.598)              
                                                         
GPA                             14,500.240***            
                                 (4,357.070)             
                                                         
libvol                            12.956***              
                                   (3.261)               
                                                         
cost                               0.366**               
                                   (0.146)               
                                                         
Constant             7.383***    -33,124.550    5.562**  
                      (2.577)   (25,210.760)    (2.141)  
                                                         
---------------------------------------------------------
Observations            136          136          136    
R2                     0.870        0.806        0.911   
Adjusted R2            0.865        0.798        0.905   
Akaike Inf. Crit.    -227.844     2,737.692    -271.086  
Bayesian Inf. Crit.  -207.455     2,758.080    -239.046  
=========================================================
Note:                         *p<0.1; **p<0.05; ***p<0.01

Conforme o post https://rpubs.com/amrofi/wooldridge_ex7_8_spec, o melhor modelo foi identificado para o modelo 3 (como indicado no video da discussão de especificação em https://youtu.be/J7NK23J54bI).

Daremos prosseguimento então para a avaliação da heteroscedasticidade neste modelo 3.

2.3 Heteroscedasticidade no modelo 3

2.3.1 Teste de White sem termos cruzados

Optamos por realizar o teste de White com a função bptest e a especificação sem termos cruzados em virtude de julgarmos pequeno o número de observações do modelo. A hipótese nula é:
H0: resíduos são homoscedásticos. A alternativa será então H1: resíduos heteroscedásticos.

> lmtest::bptest(mod3,~ log(LSAT) + log(GPA) + log(libvol) + log(cost)+
+                    I(log(LSAT)^2)+I(log(GPA)^2)+
+                    I(log(libvol)^2)+I(log(cost)^2), data = lawsch85)

    studentized Breusch-Pagan test

data:  mod3
BP = 15.673, df = 8, p-value = 0.04731

Concluo pela presença de heteroscedasticidade dos resíduos em mod3, uma vez que o p-value = 0.047 é menor que 0.10 e rejeita-se a hipótese nula de resíduos homoscedásticos. Uma alternativa direta seria aplicar a correção pela matriz de variância-covariância pela função hccm , do pacote car .

2.3.2 Correção de Var-cov conforme White

A correção pela função car::hccm (heteroskedasticity-corrected covariance matrix) admite diferentes formatos de especificação da var-cov (“hc0”,“hc1”,“hc2”,“hc3”,“hc4”), e trabalharemos aqui com a opção hc1, em seu formato mais tradicional. Para a matriz de var-cov dos parâmetros no formato

\(Var\left( {\hat \beta |X} \right) = {\left( {X'X} \right)^{ - 1}}X'\Omega X{\left( {X'X} \right)^{ - 1}}\)

têm-se então diferentes possibilidades para \[\Omega\].

• “hc0”:  matriz clássica de correção de Eicker (1963) e popularizada por White (1980);

• “hc1”, “hc2” e “hc3”: matrizes de correção sugeridas por MacKinnon e White (1985) e aperfeiçoadas para pequenas amostras conforme Long e Ervin (2000);

• “hc4”: matriz de correção conforme Cribari-Neto (2004) para aperfeiçoar a performance em pequenas amostras com presença de observações influentes.

As matriz modificam conforme a composição da matriz de ponderação:

Diferentes possibilidades para as matrizes de variância-covariância dos resíduos. Fonte: KLEIBER e ZEILEIS (2008).

> library(car) 
> #possibilidades: hccm(regressao1,type=c("hc0","hc1","hc2","hc3","hc4"))
> vcov.white0<-hccm(mod3,type=c("hc1"))
> # pelo coeftest
> lmtest::coeftest(mod3,vcov.white0)

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) 5.5619076  2.1432862  2.5950   0.01058 *  
top10       0.7002290  0.0572244 12.2365 < 2.2e-16 ***
r11_25      0.5931739  0.0467471 12.6890 < 2.2e-16 ***
r26_40      0.3744336  0.0524192  7.1431 6.463e-11 ***
r41_60      0.2623891  0.0263407  9.9613 < 2.2e-16 ***
r61.100     0.1311293  0.0181522  7.2239 4.248e-11 ***
log(LSAT)   0.8847047  0.4854581  1.8224   0.07076 .  
log(GPA)    0.0570718  0.2579492  0.2213   0.82525    
log(libvol) 0.0362878  0.0259046  1.4008   0.16372    
log(cost)   0.0010463  0.0293381  0.0357   0.97161    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> # ou no sumário
> summary(mod3)

Call:
lm(formula = log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + 
    r61.100 + log(LSAT) + log(GPA) + log(libvol) + log(cost), 
    data = lawsch85)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.294920 -0.040379 -0.001164  0.044715  0.277004 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 5.561908   2.140706   2.598   0.0105 *  
top10       0.700229   0.052517  13.333  < 2e-16 ***
r11_25      0.593174   0.039208  15.129  < 2e-16 ***
r26_40      0.374434   0.034036  11.001  < 2e-16 ***
r41_60      0.262389   0.027979   9.378 3.57e-16 ***
r61.100     0.131129   0.021049   6.230 6.42e-09 ***
log(LSAT)   0.884705   0.471736   1.875   0.0630 .  
log(GPA)    0.057072   0.239032   0.239   0.8117    
log(libvol) 0.036288   0.025998   1.396   0.1652    
log(cost)   0.001046   0.025055   0.042   0.9668    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.08558 on 126 degrees of freedom
  (20 observations deleted due to missingness)
Multiple R-squared:  0.9111,    Adjusted R-squared:  0.9047 
F-statistic: 143.4 on 9 and 126 DF,  p-value: < 2.2e-16

Será preciso pequenos passos para expressar adequadamente na tabela do stargazer. O modelo 3 inicial está colocado ao lado do modelo 3 robusto, ou seja, com correção para heteroscedasticidade.

> # fazer stargazer com dois modelos
> cov <- vcov.white0
> robust.se <- sqrt(diag(cov))
> 
> stargazer(mod3, mod3,
+           se=list(NULL, robust.se),
+           column.labels=c("MQO-mod3","mod3 robusto"), 
+           title="Título: Resultado da Regressão",
+           align=TRUE,
+           type = "text", style = "all",
+           keep.stat=c("aic","bic","rsq", "adj.rsq","n"))

Título: Resultado da Regressão
================================================
                        Dependent variable:     
                    ----------------------------
                            log(salary)         
                       MQO-mod3    mod3 robusto 
                         (1)            (2)     
------------------------------------------------
top10                  0.700***      0.700***   
                       (0.053)        (0.057)   
                      t = 13.333    t = 12.237  
                      p = 0.000      p = 0.000  
r11_25                 0.593***      0.593***   
                       (0.039)        (0.047)   
                      t = 15.129    t = 12.689  
                      p = 0.000      p = 0.000  
r26_40                 0.374***      0.374***   
                       (0.034)        (0.052)   
                      t = 11.001     t = 7.143  
                      p = 0.000      p = 0.000  
r41_60                 0.262***      0.262***   
                       (0.028)        (0.026)   
                      t = 9.378      t = 9.961  
                      p = 0.000      p = 0.000  
r61.100                0.131***      0.131***   
                       (0.021)        (0.018)   
                      t = 6.230      t = 7.224  
                      p = 0.000      p = 0.000  
log(LSAT)               0.885*        0.885*    
                       (0.472)        (0.485)   
                      t = 1.875      t = 1.822  
                      p = 0.064      p = 0.069  
log(GPA)                0.057          0.057    
                       (0.239)        (0.258)   
                      t = 0.239      t = 0.221  
                      p = 0.812      p = 0.825  
log(libvol)             0.036          0.036    
                       (0.026)        (0.026)   
                      t = 1.396      t = 1.401  
                      p = 0.166      p = 0.162  
log(cost)               0.001          0.001    
                       (0.025)        (0.029)   
                      t = 0.042      t = 0.036  
                      p = 0.967      p = 0.972  
Constant               5.562**       5.562***   
                       (2.141)        (2.143)   
                      t = 2.598      t = 2.595  
                      p = 0.011      p = 0.010  
------------------------------------------------
Observations             136            136     
R2                      0.911          0.911    
Adjusted R2             0.905          0.905    
Akaike Inf. Crit.      -271.086      -271.086   
Bayesian Inf. Crit.    -239.046      -239.046   
================================================
Note:                *p<0.1; **p<0.05; ***p<0.01

Embora do ponto de vista de resultados tenhamos as mesmas variáveis significativas que antes (no modelo 3), sabemos que foi feita a correção para heteroscedasticidade no modelo 3 robusto. Não cabe agora refazer o teste de heteroscedasticidade sobre o modelo 3 robusto.

Entretanto, levantamos algumas indagações na próxima seção, para o leitor compreender um pouco mais sobre esse problema da heteroscedasticidade. No exercício de Figueiredo (2019), concluiu-se que os resíduos do modelo 3 não são normalmente distribuídos.

> ## Teste de Jarque-Bera para normalidade
> 
> u.hat<-resid(mod3)
> library(tseries)
> JB.mod3<-jarque.bera.test(u.hat)
> JB.mod3

    Jarque Bera Test

data:  u.hat
X-squared = 20.538, df = 2, p-value = 3.47e-05

Fez-se na sequência a análise de presença de outliers.

2.4 Investigação de outliers - teste de Bonferroni para outlier

> outlierTest(mod3)

   rstudent unadjusted p-value Bonferroni p
3 -3.850327          0.0001873     0.025473

> qqPlot(mod3)

[1]  3 47

> vif(mod3)

      top10      r11_25      r26_40      r41_60     r61.100   log(LSAT) 
   3.488650    2.635882    1.730471    1.589780    1.542470    3.478267 
   log(GPA) log(libvol)   log(cost) 
   3.693818    2.225695    1.653360 

As observações 3 e 47 sugerem outliers a serem controlados. Colocarei uma dummy para ambos.

> attach(lawsch85)
> dados<-as.data.frame(cbind(salary,rank,LSAT,GPA,libvol,cost,top10,r11_25,r26_40,r41_60,r61.100))
> detach(lawsch85)
> attach(dados)
> dados$dummy<-0
> dados$dummy[3]=1
> dados$dummy[47]=1
> dados$dummy[79]=1 # detectei no outro post que o 79 também é outlier

Farei novas regressões partindo de mod3 para inclusão de dummies alterando, para os outliers, intercepto e inclinações em diferentes variáveis.

> attach(dados)
> mod4 <- lm(log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + r61.100 
+            + log(LSAT) + log(GPA) + log(libvol) + log(cost)+dummy+
+              I(dummy*LSAT), data = dados)
> mod4$AIC <- AIC(mod4)
> mod4$BIC <- BIC(mod4)
> mod5 <- lm(log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + r61.100 
+            + log(LSAT) + log(GPA) + log(libvol) + log(cost)+dummy+
+              I(dummy*GPA), data = dados)
> mod5$AIC <- AIC(mod5)
> mod5$BIC <- BIC(mod5)
> mod6 <- lm(log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + r61.100 
+            + log(LSAT) + log(GPA) + log(libvol) + log(cost)+dummy+
+              I(dummy*libvol), data = dados)
> mod6$AIC <- AIC(mod6)
> mod6$BIC <- BIC(mod6)
> mod7 <- lm(log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + r61.100 
+            + log(LSAT) + log(GPA) + log(libvol) + log(cost)+dummy+
+              I(dummy*cost), data = dados)
> mod7$AIC <- AIC(mod7)
> mod7$BIC <- BIC(mod7)
> library(stargazer)
> star.2 <- stargazer(mod4, mod5,mod6,mod7,
+                     title="Título: Resultados das Regressões",
+                     column.labels=c("MQO-mod4","MQO-mod5","MQO-mod6","MQO-mod7"),
+                     align=TRUE,
+                     type = "text",
+                     keep.stat=c("aic","bic","rsq", "adj.rsq","n")
+ )

Título: Resultados das Regressões
==========================================================
                             Dependent variable:          
                    --------------------------------------
                                 log(salary)              
                    MQO-mod4  MQO-mod5 MQO-mod6  MQO-mod7 
                       (1)      (2)       (3)       (4)   
----------------------------------------------------------
top10               0.717***  0.713*** 0.752***  0.745*** 
                     (0.051)  (0.052)   (0.047)   (0.046) 
                                                          
r11_25              0.604***  0.602*** 0.630***  0.624*** 
                     (0.038)  (0.039)   (0.035)   (0.035) 
                                                          
r26_40              0.356***  0.356*** 0.412***  0.389*** 
                     (0.033)  (0.034)   (0.031)   (0.030) 
                                                          
r41_60              0.268***  0.267*** 0.285***  0.280*** 
                     (0.027)  (0.027)   (0.025)   (0.025) 
                                                          
r61.100             0.134***  0.134*** 0.148***  0.143*** 
                     (0.020)  (0.021)   (0.019)   (0.018) 
                                                          
log(LSAT)            0.853*    0.898*   0.923**   0.774*  
                     (0.461)  (0.468)   (0.417)   (0.417) 
                                                          
log(GPA)             -0.003    0.005    -0.090    -0.074  
                     (0.231)  (0.235)   (0.211)   (0.210) 
                                                          
log(libvol)           0.035    0.034     0.023     0.031  
                     (0.025)  (0.026)   (0.023)   (0.023) 
                                                          
log(cost)            -0.005    -0.005   -0.017    -0.011  
                     (0.025)  (0.025)   (0.022)   (0.022) 
                                                          
dummy               -5.833*** -2.467** 8.269***  -8.352***
                     (2.031)  (1.211)   (1.368)   (1.385) 
                                                          
I(dummy * LSAT)     0.038***                              
                     (0.013)                              
                                                          
I(dummy * GPA)                0.775**                     
                              (0.365)                     
                                                          
I(dummy * libvol)                      -0.020***          
                                        (0.003)           
                                                          
I(dummy * cost)                                  0.0005***
                                                 (0.0001) 
                                                          
Constant            5.852***  5.623*** 5.772***  6.410*** 
                     (2.087)  (2.118)   (1.890)   (1.889) 
                                                          
----------------------------------------------------------
Observations           136      136       136       136   
R2                    0.919    0.917     0.933     0.934  
Adjusted R2           0.912    0.909     0.927     0.928  
Akaike Inf. Crit.   -280.148  -275.942 -305.474  -306.833 
Bayesian Inf. Crit. -242.284  -238.077 -267.610  -268.968 
==========================================================
Note:                          *p<0.1; **p<0.05; ***p<0.01

O modelo 7 aparenta melhor ajustamento e continuarei análises a partir dele. Testarei novamente para normalidade:

> dados$obs<-1:156
> u.hat<-resid(mod7)
> library(tseries)
> JB.mod7<-jarque.bera.test(u.hat)
> JB.mod7

    Jarque Bera Test

data:  u.hat
X-squared = 0.44485, df = 2, p-value = 0.8006

> outlierTest(mod7)

No Studentized residuals with Bonferroni p < 0.05
Largest |rstudent|:
   rstudent unadjusted p-value Bonferroni p
96 2.444258           0.015933           NA

> qqPlot(mod7)

[1]  96 134

Agora sim, colocando 3, 47 e 79 como dummies consegui resíduos normais no modelo 7. Vou olhar o teste RESET para especificação.

> TesteRESET.power<-lmtest::resettest(mod7, power = 2:3)
> TesteRESET.power

    RESET test

data:  mod7
RESET = 2.2131, df1 = 2, df2 = 122, p-value = 0.1137

Embora ainda com probabilidade baixa, vou considerar o modelo 7 bem especificado e investigar outros itens como a heteroscedasticidade.

2.5 Teste de heteroscedasticidade no modelo 7

Vou checar a heterocedasticidade. A dummy somente entra em nível, não adiantando colocar ao quadrado.

> lmtest::bptest(mod7,~ log(LSAT) + log(GPA) + log(libvol) + log(cost)
+ 
+                  +dummy+I(dummy*cost)+
+ 
+                    I(log(LSAT)^2)+I(log(GPA)^2)+
+ 
+                    I(log(libvol)^2)+I(log(cost)^2), data = dados)

    studentized Breusch-Pagan test

data:  mod7
BP = 14.633, df = 10, p-value = 0.146

Não detecto heterocedasticidade. OU seja, ao controlar pelos outliers, resolvi o problema de heterocedasticidade. Como a probabilidade foi próxima a 0.10, sugiro investigar outras possibilidades. Farei o teste padrão de BP.

> BP_t <- lmtest::bptest(mod7)
> print(BP_t)

    studentized Breusch-Pagan test

data:  mod7
BP = 29.246, df = 11, p-value = 0.002079

Ou seja, o modelo apresenta heteroscedasticidade sem controlar outliers e como visto, pode obter outro resultado se optar por controlar para os outliers. O pesquisador deve estar atento pois cada caso pode levar a alternativas diferentes. Como o teste BP acusou heterocedasticidade, faz-se a correção e estimação robusta. Agora utilizaremos o pacote sandwich e a função vcovHC. Esta função opera de modo semelhante ao do pacote car e também se pode especificar tipos de ponderação. Usarei a opção HC1 para ficar como na forma do hccm.

> V_HC_1 <- sandwich::vcovHC(mod7, type = "HC1")

Vamos apresentar a tabela de coeficientes com erros robustos:

> print(lmtest::coeftest(mod7, vcov. = V_HC_1))

t test of coefficients:

                   Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)      6.4102e+00  1.8605e+00  3.4454 0.0007788 ***
top10            7.4497e-01  4.7493e-02 15.6859 < 2.2e-16 ***
r11_25           6.2351e-01  4.1367e-02 15.0726 < 2.2e-16 ***
r26_40           3.8881e-01  4.0996e-02  9.4841 2.243e-16 ***
r41_60           2.7962e-01  2.1732e-02 12.8667 < 2.2e-16 ***
r61.100          1.4318e-01  1.4309e-02 10.0062 < 2.2e-16 ***
log(LSAT)        7.7380e-01  4.1461e-01  1.8663 0.0643582 .  
log(GPA)        -7.3788e-02  2.4278e-01 -0.3039 0.7616933    
log(libvol)      3.1176e-02  2.1683e-02  1.4378 0.1530143    
log(cost)       -1.1273e-02  2.4602e-02 -0.4582 0.6476093    
dummy           -8.3518e+00  9.5763e-01 -8.7214 1.509e-14 ***
I(dummy * cost)  4.9380e-04  5.7356e-05  8.6095 2.782e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

E agora com os erros viesados pela heteroscedasticidade:

> print(lmtest::coeftest(mod7))

t test of coefficients:

                   Estimate  Std. Error t value  Pr(>|t|)    
(Intercept)      6.4102e+00  1.8893e+00  3.3928 0.0009291 ***
top10            7.4497e-01  4.6412e-02 16.0514 < 2.2e-16 ***
r11_25           6.2351e-01  3.4590e-02 18.0257 < 2.2e-16 ***
r26_40           3.8881e-01  3.0177e-02 12.8843 < 2.2e-16 ***
r41_60           2.7962e-01  2.4595e-02 11.3690 < 2.2e-16 ***
r61.100          1.4318e-01  1.8482e-02  7.7472 2.894e-12 ***
log(LSAT)        7.7380e-01  4.1695e-01  1.8559 0.0658484 .  
log(GPA)        -7.3788e-02  2.0985e-01 -0.3516 0.7257230    
log(libvol)      3.1176e-02  2.2902e-02  1.3612 0.1759067    
log(cost)       -1.1273e-02  2.2222e-02 -0.5073 0.6128562    
dummy           -8.3518e+00  1.3855e+00 -6.0281 1.760e-08 ***
I(dummy * cost)  4.9380e-04  8.0869e-05  6.1062 1.210e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Embora tenhamos alteracoes nos desvios padrões e consequentes estatísticas t e probabilidades, o resultado por MQG - HCE (heteroskedastic consistent estimator) neste caso não alterou a interpretação relativamente ao modelo por MQO.

2.6 Outra estimação robusta pelo estimatr

> library(estimatr)
> library(car)
> # Erros padrões robustos HC1 – padrão do Stata
> model <- lm_robust(log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + r61.100 
+            + log(LSAT) + log(GPA) + log(libvol) + log(cost)+dummy+
+              I(dummy*cost), data = dados,
+                    se_type = "stata")
> summary(model)

Call:
lm_robust(formula = log(salary) ~ top10 + r11_25 + r26_40 + r41_60 + 
    r61.100 + log(LSAT) + log(GPA) + log(libvol) + log(cost) + 
    dummy + I(dummy * cost), data = dados, se_type = "stata")

Standard error type:  HC1 

Coefficients:
                  Estimate Std. Error t value  Pr(>|t|)   CI Lower   CI Upper
(Intercept)      6.4101833  1.861e+00  3.4454 7.788e-04  2.728e+00 10.0926428
top10            0.7449703  4.749e-02 15.6859 3.173e-31  6.510e-01  0.8389722
r11_25           0.6235082  4.137e-02 15.0726 8.242e-30  5.416e-01  0.7053849
r26_40           0.3888092  4.100e-02  9.4841 2.243e-16  3.077e-01  0.4699515
r41_60           0.2796213  2.173e-02 12.8667 1.375e-24  2.366e-01  0.3226353
r61.100          0.1431794  1.431e-02 10.0062 1.222e-17  1.149e-01  0.1715011
log(LSAT)        0.7738031  4.146e-01  1.8663 6.436e-02 -4.683e-02  1.5944390
log(GPA)        -0.0737878  2.428e-01 -0.3039 7.617e-01 -5.543e-01  0.4067455
log(libvol)      0.0311757  2.168e-02  1.4378 1.530e-01 -1.174e-02  0.0740928
log(cost)       -0.0112726  2.460e-02 -0.4582 6.476e-01 -5.997e-02  0.0374214
dummy           -8.3518485  9.576e-01 -8.7214 1.509e-14 -1.025e+01 -6.4564349
I(dummy * cost)  0.0004938  5.736e-05  8.6095 2.782e-14  3.803e-04  0.0006073
                 DF
(Intercept)     124
top10           124
r11_25          124
r26_40          124
r41_60          124
r61.100         124
log(LSAT)       124
log(GPA)        124
log(libvol)     124
log(cost)       124
dummy           124
I(dummy * cost) 124

Multiple R-squared:  0.9336 ,   Adjusted R-squared:  0.9277 
F-statistic:   169 on 11 and 124 DF,  p-value: < 2.2e-16

3 Referências

FIGUEIREDO, Adriano Marcos Rodrigues. Econometria: Salários de advogados - Wooldridge(exemplo 7.8) adaptado de MOHR (2018). Campo Grande-MS,Brasil: RStudio/Rpubs, 2019. Disponível em <http://rpubs.com/amrofi/wooldridge_ex7_8%3E.

HEISS, Florian. Using R for Introductory Econometrics. 2.ed. Florian Heiss, 2020. Recurso online. Disponível em: http://www.urfie.net/.

KLEIBER, Christian;·ZEILEIS, Achim. Applied Econometrics with R. New York-NY: Springer Science+Business Media, LLC, 2008.

SHEA, Justin M. wooldridge: vignette. 2020. Disponível em: https://justinmshea.github.io/wooldridge/articles/Introductory-Econometrics-Examples.html.

SHEA, Justin M. wooldridge: 111 Data Sets from “Introductory Econometrics: A Modern Approach, 6e” by Jeffrey M. Wooldridge. R package version 1.3.1. 2018. Disponível em: https://CRAN.R-project.org/package=wooldridge.

WOOLDRIDGE, J.M. Introdução à Econometria: uma abordagem moderna. São Paulo: Pioneira Thomson Learning, 2006.(tradução da segunda edição americana).

WOOLDRIDGE, Jeffrey M. Introductory econometrics: A modern approach. Nelson Education, 2016.