0.1 Dados

library(tidyverse)

msp <- read_rds("C:\\Users\\rauld\\Google Drive\\meu_projeto\\dados e scripts\\tabelas_output\\tab_FINAL.rds")

msp <- msp   %>% 
  select(ano,
         distrito,
         pop_jovem_onu,
         homicidio,
         populacao,
         area_constr_resid_baixop,
         renda_domic_media,
         renda_domic_dp,
         emprego_formal,
         armas_apreendidas,
         favela)                                                %>%
  filter(ano == "2003")                                         %>%
  mutate(tx_homicidio = (homicidio/populacao)*100000,
         eformais     = (emprego_formal/populacao)*100000,
         tx_belica    = (armas_apreendidas/populacao)*100000)   %>% 
  rename(jovens1524   = pop_jovem_onu)


# o modelo estimado por MQO
ols <- lm(tx_homicidio ~ jovens1524 + area_constr_resid_baixop + renda_domic_media + renda_domic_dp + favela + eformais + tx_belica, data = msp)
 
msp <- msp %>% mutate(residuos = ols$residuals) %>% glimpse()
## Observations: 80
## Variables: 15
## $ ano                      <int> 2003, 2003, 2003, 2003, 2003, 2003, 2...
## $ distrito                 <chr> "SÉ", "BOM RETIRO", "CAMPOS ELÍSIOS",...
## $ jovens1524               <dbl> 3927, 5061, 8225, 20268, 10795, 4926,...
## $ homicidio                <dbl> 53, 6, 45, 4, 20, 13, 18, 23, 15, 17,...
## $ populacao                <dbl> 21079, 28656, 50595, 120821, 64310, 3...
## $ area_constr_resid_baixop <dbl> 115205, 223771, 365532, 344136, 16226...
## $ renda_domic_media        <dbl> 9.155, 11.321, 11.291, 15.672, 14.405...
## $ renda_domic_dp           <dbl> 9.114, 8.559, 8.564, 9.052, 8.689, 8....
## $ emprego_formal           <dbl> 79883, 40115, 76795, 114601, 24619, 2...
## $ armas_apreendidas        <dbl> 110, 68, 103, 61, 73, 94, 122, 57, 13...
## $ favela                   <dbl> 0, 400, 0, 0, 0, 0, 65, 0, 36, 5170, ...
## $ tx_homicidio             <dbl> 251.435078, 20.938023, 88.941595, 3.3...
## $ eformais                 <dbl> 378969.591, 139988.135, 151783.773, 9...
## $ tx_belica                <dbl> 521.84639, 237.29760, 203.57743, 50.4...
## $ residuos                 <dbl> 38.589179, -59.712277, 16.604226, 4.8...

0.2 Introdução

No contexto da heterocedasticidade, as variâncias para todas as observações não são as mesmas. Quebra-se a hipótese de homocedasticidade. Quais as suas consequências? Como detectá-la? E como corrigir?

0.2.1 Quais as suas consequências?

No contexto da homocedasticidade, temos violada a condição \(var(e_i)=\sigma^2\):

\[var(e_i)=\sigma^{2}_{i}\]

  • O estimador por MQO ainda é enviesado, mas não é mais o melhor, pois o melhor é o que tem a menor variância.

  • Os erros padrão obtidos pelo estimador MQO é incorreto, afetando os intervaos de confiânça e, consequentemente, os testes de hipóteses.

A heterocedásticidade ocorre com frequência. Mas o MQO sem correção para heterocedásticidade ainda pode ser um estimador razoável no contexto das grandes amostras, onde as variâncias do MQO são pequenas o suficiente para estimação de parãmetros mais precisos.

0.3 Detectando:

0.3.1 Análise de resíduos

É uma maneira informal de detecção: faça um plot nos resíduos dos etimadores por MQO.

msp %>%
  mutate(residuos = ols$residuals) %>% 
  ggplot(data = ., aes(y = residuos, x = tx_homicidio)) + 
   geom_point() + 
   geom_abline(slope = 0) +
   theme_classic()

Taxas de homicídio mais altas não apresentam resíduos negativos. Taxas de homicidio abaixo de (aprox.) 80 sem padrão evidente.

0.3.2 Como detectá-la? O teste Breusch-Pagan

Aplicamos um teste \(\chi^2\) sobre a função de variãncia:


-\(H_{0}\): não há heterocedasticidade, logo \(\alpha_2=\alpha_3=\ldots\ =\alpha_s=0\)
-\(H_{1}\): c.c. (pelo menos um \(\alpha \neq 0\))


  • Função de variância: \[ \begin{aligned} var(y_i)=\sigma_i^2=E(e_i^2)=h(\alpha_1+\alpha_2z_{i2}+\alpha_3z_{i3}+ \ldots\ +\alpha_sz_{is}) \end{aligned}\]


Neste caso, \(z_{i1},z_{i2},\ldots\ ,z_{is}\) é um conjunto de variáveis explicativas da função de variância que são diferentes de \(x_{i1},x_{i2},\ldots\ ,x_{is}\).


  • se \(\chi^2≥\chi_{(1-\alpha,S-1)}^2\) rejeitamos \(H_{0}\).


Para obter a estatística:


\[\begin{aligned} e_i^2=E(e_i^2)=\alpha_1+\alpha_2z_{i2}+ \ldots\ +\alpha_sz_{is} \end{aligned}\]


Assumimos \(v_i=e_i^2-E(e_i^2)\), portanto:


\[\begin{aligned} e_i^2=E(e_i^2)+v_i=\alpha_1 + \alpha_2z_{i2} + \ldots\ +\alpha_sz_{is} + v_i \end{aligned}\]


Estimando \(\hat{e}_i^2\), temos:


\[\begin{aligned} \hat{e}_i^2=\alpha_1 + \alpha_2z_{i2} + \ldots + \alpha_sz_{is} + v_i \end{aligned}\]


Se \(H_{0}\) é vrdadeira, então o número de observações \(N\) multiplicado pelo \(R^{2}\) (que mede a proporção da variância em \(\hat{e}_i^2\) que é explicada por \(z\)). Este produto seque uma distribuição \(\chi^2\) com \(S-1\) graus de liberdade. Ou seja:


\[\begin{aligned} \chi^2=N \times R^2 \sim \chi_{S-1}^2 \end{aligned}\]


0.3.3 Aplicando

var_func <- lm(residuos^2 ~ tx_homicidio, data = msp)
summary(var_func)
## 
## Call:
## lm(formula = residuos^2 ~ tx_homicidio, data = msp)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -813.9 -310.6 -234.6   62.0 3261.3 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   217.114    106.013   2.048   0.0439 *
## tx_homicidio    4.163      1.711   2.433   0.0173 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 665.6 on 78 degrees of freedom
## Multiple R-squared:  0.07055,    Adjusted R-squared:  0.05863 
## F-statistic:  5.92 on 1 and 78 DF,  p-value: 0.01726

\(R^{2} = 0.05863\), com \(N=80\), o que nos da um valor \(\chi_{(1-\alpha,S-1)}^2=\) = 4.6904. Vamos extrair agora o valor crítico (tabelado), com nível de significância a 5%:

qchisq(.95, df = 1)
## [1] 3.841459


Como \(\chi^2 < \chi_{(1-\alpha,S-1)}^2\), a hipotese nula é rejeitada e é constatada a presença de heterocedasticidade.

De um modo mais fácil, o R faz o teste com ajuda da função lmtest::bptest():

library(lmtest)
bptest(ols)
## 
##  studentized Breusch-Pagan test
## 
## data:  ols
## BP = 19.314, df = 7, p-value = 0.007259

Como 0.0072591 < 0.05, nós rejeitamos a hipótese nula de homocedasticidade.

0.3.4 E como corrigir? Regressão com erros robustos

lembre-se que as consequencias são duas: 1) os estimadores MQO não são o de menor variãncia; 2) Erros padrões não são confiáveis para teste de hipótese.

0.3.4.1 heteroskedasticity-consistent standard errors or simply robust standard errors

Obtemos, primeiro, a variancia do estimador MQO \(b_{2}\):


\[\begin{aligned} var(b_2)=\frac{\sum_{i=1}^{N}\big[(x_i-\bar{x})^2\sigma_i^2\big]}{\big[\sum_{i=1}^{N}(x_i-\bar{x}^2)\big]^2} \end{aligned}\]


Na equação acima, devemos substituir o \(\sigma_i^2\) pelo quadrado dos resíduos MQO \(\hat{e}_i^2=y_i-b_1-b_2x_i\), incluindo graus de liberdade dado por \(N/N-K\), onde \(K\) representa o número de parâmetros no modelo. No nosso modelo temos \(K\) = 8.


\[\begin{aligned} \hat{var(b_2)}=\frac{N}{N-2} \frac{\sum_{i=1}^N\big[(x_i-\bar{x})^2\hat{e}_i^2\big]}{\big[\sum_{i=1}^N(x_i-\bar{x})^2\big]^2} \end{aligned}\]


No modelo por MQO:

summary(ols)
## 
## Call:
## lm(formula = tx_homicidio ~ jovens1524 + area_constr_resid_baixop + 
##     renda_domic_media + renda_domic_dp + favela + eformais + 
##     tx_belica, data = msp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -59.712  -9.297   0.200  10.516  46.901 
## 
## Coefficients:
##                            Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               3.014e+01  3.359e+01   0.897  0.37254    
## jovens1524                6.535e-05  3.976e-04   0.164  0.86991    
## area_constr_resid_baixop -1.275e-05  1.040e-05  -1.226  0.22415    
## renda_domic_media        -6.457e+00  8.211e-01  -7.864 2.75e-11 ***
## renda_domic_dp            4.734e+00  3.124e+00   1.516  0.13401    
## favela                    1.288e-03  8.222e-04   1.566  0.12162    
## eformais                  1.792e-04  6.032e-05   2.971  0.00404 ** 
## tx_belica                 2.529e-01  2.734e-02   9.251 7.13e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.1 on 72 degrees of freedom
## Multiple R-squared:  0.7882, Adjusted R-squared:  0.7676 
## F-statistic: 38.27 on 7 and 72 DF,  p-value: < 2.2e-16

Com a função sandwich::coeftest() obtemos o erro padrão robusto:

## 
## t test of coefficients:
## 
##                             Estimate  Std. Error t value  Pr(>|t|)    
## (Intercept)               3.0142e+01  3.2752e+01  0.9203   0.36049    
## jovens1524                6.5347e-05  4.3329e-04  0.1508   0.88054    
## area_constr_resid_baixop -1.2750e-05  1.0614e-05 -1.2012   0.23362    
## renda_domic_media        -6.4571e+00  1.0243e+00 -6.3037 2.066e-08 ***
## renda_domic_dp            4.7339e+00  2.9979e+00  1.5791   0.11870    
## favela                    1.2880e-03  8.0128e-04  1.6075   0.11233    
## eformais                  1.7919e-04  9.7142e-05  1.8446   0.06921 .  
## tx_belica                 2.5291e-01  2.9616e-02  8.5396 1.507e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Os erros padrões:

coeficiente mqo white
(Intercept) 33.59 32.75
jovens1524.00 0.00 0.00
area_constr_resid_baixop 0.00 0.00
renda_domic_media 0.82 1.02
renda_domic_dp 3.12 3.00
favela 0.00 0.00
eformais 0.00 0.00
tx_belica 0.03 0.03