Nota prévia de leitura

Antes que você comece a ler essa minha pequena aventura, acho que é muito importante ressaltar que todos os posts aqui no blog são escritos diretamente no R, usando o pacote RMarkdown - mesmo quando usamos algo de python, Julia ou, nesse caso, Stata. O Daniel tem um post bom explicando o nosso workflow de maneira detalhada disponível preeliminarmente aqui.

O paper

Dia desses saiu a edição de Novembro da American Economic Review e nela um paper me chamou muito à atenção: The Costs of Patronage: Evidence from the British Empire, de Guo Xu. O resumo me deixou dando pulinhos de alegria:

“I combine newly digitized personnel and public finance data from the British colonial administration for the period 1854-1966 to study how patronage affects the promotion and incentives of governors. Governors are more likely to be promoted to higher salaried colonies when connected to their superior during the period of patronage. Once allocated, they provide more tax exemptions, raise less revenue, and invest less. The promotion and performance gaps disappear after the abolition of patronage appointments. Patronage therefore distorts the allocation of public sector positions and reduces the incentives of favored bureaucrats to perform.”

Tem tudo aí. Uma base de dados inédita, uma questão interessante, uma abordagem claramente econômica (incentivos, regras) e a quantificação do impacto de um comportamento. Se alguém arranjar um anão e uma garrafa de champagne vira uma festa daquelas. Pois, como eu tenho muita coisa pra fazer, prova para estudar e trabalho para entregar, decidi replicar nem que seja uma regressão desse paper tão bonito - e mais importante ainda - deixar uma base pronta para quem quiser replicar o resto.

Extraindo os dados

Se você está interessado em ter a base de dados maneira que o autor desse paper conseguiu, é fácil. O código a seguir puxa eles diretamente do repositório da AER.

library(readstata13) # importante para ler o arquivo com a base
library(tibble) # não tão importante, mas tibbles > data.frames

temp = tempfile() ## geramos um arquivo temporário
download.file("https://www.aeaweb.org/doi/10.1257/aer.20171339.data",
              destfile = temp, mode = "wb") # puxamos diretamente os dados do repositório da AEA

dados = unzip(temp, file = "analysis.dta") # descomprimimos
dados = as.tibble(read.dta13(dados)) #lemos como um tibble

dim(dados)

## [1] 4687   55

dados #exploramos para garantir que está tudo certo

## # A tibble: 4,687 x 55
##       sid  year                      aid log_salary_governor_gbp
##  * <fctr> <dbl>                   <fctr>                   <dbl>
##  1   Aden  1938    Bernard Rawdon Reilly                8.086411
##  2   Aden  1937    Bernard Rawdon Reilly                8.086411
##  3   Aden  1939    Bernard Rawdon Reilly                8.086411
##  4   Aden  1944        John Hathorn Hall                8.086411
##  5   Aden  1941        John Hathorn Hall                8.086411
##  6   Aden  1942        John Hathorn Hall                8.086411
##  7   Aden  1940        John Hathorn Hall                8.086411
##  8   Aden  1943        John Hathorn Hall                8.086411
##  9   Aden  1946 Reginald Stuart Champion                8.086411
## 10   Aden  1947 Reginald Stuart Champion                8.086411
## # ... with 4,677 more rows, and 51 more variables:
## #   log_rev_total_gbp <dbl>, log_rev_customs_gbp <dbl>,
## #   log_rev_internal_gbp <dbl>, log_exp_total_gbp <dbl>,
## #   log_exp_tax_gbp <dbl>, log_exp_pubworks_gbp <dbl>,
## #   log_population_imputed <dbl>, shared_ancestry <dbl>, both_arist <dbl>,
## #   both_oxbridge <dbl>, both_eton <dbl>, connected <dbl>,
## #   age_entry <dbl>, cambridge <int>, oxford <int>, eton <int>,
## #   peerage <dbl>, civilservant <dbl>, politician <dbl>, military <dbl>,
## #   ordinance_total <dbl>, ordinance_health <dbl>,
## #   ordinance_education <dbl>, ordinance_pubworks <dbl>,
## #   ordinance_welfare <dbl>, ordinance_tax <dbl>, ordinance_trade <dbl>,
## #   customs_exemptions <int>, social_unrest <dbl>, hansard_mention <dbl>,
## #   polarity <dbl>, award_highest <dbl>, log_dist_london <dbl>,
## #   area_tropics <dbl>, logem4 <dbl>, sum_d_pre_post <dbl>,
## #   post1930 <dbl>, tenure_aid <dbl>, landlocked <int>,
## #   quinquennial <dbl>, decade <dbl>, bilateral <dbl>, no_colonies <dbl>,
## #   duration <dbl>, state_aid <dbl>, min_year <dbl>, min_revenue <dbl>,
## #   full <dbl>, post1930_connected <dbl>, connected_year1930 <dbl>,
## #   tenure <dbl>

Temos um painel que cobre boa parte do século XIX e início do XX. As informações são bem detalhadas, desde nome do governador da colônia, onde ele estudou, se era militar e até distância da colônia até Londres.

library(dplyr)

dados %>%
  group_by(decade) %>%
  summarise(idade = round(mean(age_entry, na.rm = TRUE), digits = 2), 
            politico = round(mean(politician, na.rm = TRUE), digits = 2), 
            militar = round(mean(military, na.rm = TRUE), digits = 2))

## # A tibble: 13 x 4
##    decade idade politico militar
##     <dbl> <dbl>    <dbl>   <dbl>
##  1   1854 43.52     0.19    0.45
##  2   1860 41.96     0.20    0.38
##  3   1870 41.92     0.16    0.37
##  4   1880 43.97     0.12    0.39
##  5   1890 45.61     0.13    0.44
##  6   1900 45.26     0.07    0.44
##  7   1910 46.38     0.07    0.44
##  8   1920 49.82     0.06    0.40
##  9   1930 50.36     0.01    0.35
## 10   1940 48.55     0.04    0.38
## 11   1950 47.01     0.00    0.29
## 12   1960 48.65     0.02    0.25
## 13     NA 45.48     0.16    0.54

Queria, em particular, reproduzir a tabela 2 em que salários de governadores de colônias são explicadas por modelos com alguns termos de efeitos fixos e variáveis de conectividade com o encarregado dos apontamentos. Essas variáveis incluem ancentrais comuns, ser aristocrata e ter estudado em Eton, Oxford ou Cambridge junto com quem fez os apontamentos.

A loucura do Stata, ou como eu virei um herege

Acabei fuçando os arquivos para reprodução que a própria AER fornece e achei o código original em Stata - que por sinal poderia ser executado no próprio RMarkdown com a ajuda do pacote Statamarkdown e vou mostrar como fazer isso.

Primeiro carregamos a biblioteca, ainda em R:

library(Statamarkdown)

Depois armazenamos um objeto com nome stataexe o endereço do Stata no seu computador. Se tiver instalado, pode usar a função Statamarkdown::find_stata() para isso.

stataexe <- "C:/Users/Pedro/Desktop/Meus Dados/Stata14/StataMP-64.exe"

Agora você avisa o knitr que está realmente fazendo essa loucura de rodar Stata no RStudio:

knitr::opts_chunk$set(engine.path = stataexe)

Voi lá:

use "analysis.dta", replace

sum log_salary_governo, detail

                  ln governor salary (gbp)
-------------------------------------------------------------
      Percentiles      Smallest
 1%     5.703783       5.521461
 5%      6.39693       5.521461
10%     6.684612       5.521461       Obs               4,325
25%     7.495542       5.521461       Sum of Wgt.       4,325

50%     8.086411                      Mean           7.957734
                        Largest       Std. Dev.      .7848033
75%     8.517193       9.545025
90%     8.853665       9.546813       Variance       .6159162
95%     8.922658       9.546813       Skewness      -.7752303
99%      9.21034       9.546813       Kurtosis       3.013498

Para de fato rodar o código original, precisamos do pacote reghdfe, do Sergio Correia - disponível no Github. Então vamos instala-lo com o código sugerido pelo próprio criador:

* Install ftools (remove program if it existed previously)
cap ado uninstall ftools
net install ftools, from("https://raw.githubusercontent.com/sergiocorreia/ftools/master/src/")

* Install reghdfe 4.x
cap ado uninstall reghdfe
net install reghdfe, from("https://raw.githubusercontent.com/sergiocorreia/reghdfe/master/src/")

* Install boottest for Stata 11 and 12
if (c(version)<13) cap ado uninstall boottest
if (c(version)<13) ssc install boottest

* Install moremata (sometimes used by ftools but not needed for reghdfe)
cap ssc install moremata

ftools, compile
reghdfe, compile

Agora sim podemos estimar os modelos com o código original em Stata, que estima todas as especificações da tabela 2. Vou colocar * antes das linhas que seguem a primeira para que o documento não renderize todas as tabelas.

use "analysis.dta"

reghdfe log_salary_governor_gbp no_colonies shared_ancestry if full==1, absorb(aid year duration) vce(cluster bilateral)

* reghdfe log_salary_governor_gbp no_colonies both_arist if full==1, absorb(aid year duration) vce(cluster bilateral)

* reghdfe log_salary_governor_gbp no_colonies both_eton if full==1, absorb(aid year duration) vce(cluster bilateral)

* reghdfe log_salary_governor_gbp no_colonies both_oxbridge if full==1, absorb(aid year duration) vce(cluster bilateral)

* reghdfe log_salary_governor_gbp no_colonies shared_ancestry both_arist both_eton both_oxbridge if full==1, absorb(aid year duration) vce(cluster bilateral)

* reghdfe log_salary_governor_gbp no_colonies connected if full==1, absorb(aid year duration) vce(cluster bilateral)

> b(aid year duration) vce(cluster bilateral)
(MWFE estimator converged in 26 iterations)

HDFE Linear regression                            Number of obs   =      3,510
Absorbing 3 HDFE groups                           F(   2,   1517) =      22.90
Statistics robust to heteroskedasticity           Prob > F        =     0.0000
                                                  R-squared       =     0.9253
                                                  Adj R-squared   =     0.9108
                                                  Within R-sq.    =     0.0962
Number of clusters (bilateral) =      1,518       Root MSE        =     0.2376

                          (Std. Err. adjusted for 1,518 clusters in bilateral)
------------------------------------------------------------------------------
             |               Robust
log_salary~p |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
 no_colonies |   .2207025   .0350703     6.29   0.000     .1519111    .2894939
shared_anc~y |   .1034715   .0470474     2.20   0.028     .0111866    .1957563
       _cons |   7.496656   .0664925   112.74   0.000     7.366229    7.627083
------------------------------------------------------------------------------

Absorbed degrees of freedom:
-----------------------------------------------------+
 Absorbed FE | Categories  - Redundant  = Num. Coefs |
-------------+---------------------------------------|
         aid |       456           0         456     |
        year |       110           1         109     |
    duration |         7           1           6    ?|
-----------------------------------------------------+
? = number of redundant parameters may be higher

> id year duration) vce(cluster bilateral)
> d year duration) vce(cluster bilateral)
> b(aid year duration) vce(cluster bilateral)
> _eton both_oxbridge if full==1, absorb(aid year duration) vce(cluster bilater
> al)
> d year duration) vce(cluster bilateral)

Encontrei duas dificuldades:

• O console sai no documento então para ver o resultado de uma linha de código, é preciso compilar o documento todo e conferir a saída lá

• Cada chunk de código em Stata gera uma sessão nova, então é preciso importar os dados sempre que abrir um chunk novo

E como fazer isso no R?

Vamos usar o pacote lfe, que implementa um método muito similar com a função felm. Como não conheço nem localizei implementação funcional para R de erros-padrão agrupados então a estimação não é estritamente a mesma, mas o espírito está lá.

library(lfe)

modelo1 = felm(log_salary_governor_gbp ~ no_colonies + shared_ancestry | aid + year + duration, 
               data = dados,
               subset = (full == 1))

summary(modelo1)

Call:
   felm(formula = log_salary_governor_gbp ~ no_colonies + shared_ancestry |      aid + year + duration, data = dados, subset = (full == 1)) 

Residuals:
     Min       1Q   Median       3Q      Max 
-1.26531 -0.07417  0.00110  0.08623  0.96266 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
no_colonies      0.22070    0.01289  17.123  < 2e-16 ***
shared_ancestry  0.10347    0.02425   4.267 2.04e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.2376 on 2937 degrees of freedom
Multiple R-squared(full model): 0.9253   Adjusted R-squared: 0.9108 
Multiple R-squared(proj model): 0.09618   Adjusted R-squared: -0.07985 
F-statistic(full model):63.63 on 572 and 2937 DF, p-value: < 2.2e-16 
F-statistic(proj model): 156.3 on 2 and 2937 DF, p-value: < 2.2e-16 
*** Standard errors may be too high due to more than 2 groups and exactDOF=FALSE

Observe que o código que estima o modelo1 é praticamente o mesmo que estima as outras especificações, basta somente adicionar ou retirar variáveis explicativas na fórmula.