Momento R do Dia – Replicando um exercício…com outra base de dados

Uma das maneiras mais simples de se aprender qualquer coisa – inclusive R – é ganhar confiança replicando um exemplo, ainda que você não tenha condições, no início, de entender todas as implicações do que está fazendo. Funciona bem com tudo na vida (ou você sabia andar de bicicleta antes de subir nela pela primeira vez?) e, claro, também no R.

Depois de replicar, claro, é legal mudar a base de dados e treinar novamente. Assim, após replicar o exemplo de cross-validation do BlogR, resolvi, como economista, fazer um exercício de, digamos, estática comparativa. Tudo o mais constante, mudei a base de dados. Mas o que vamos fazer? O que é cross-validation, ou melhor, o que é leave-one-out cross-validation? Citando o autor:

Leave-one-out is a type of cross validation whereby the following is done for each observation in the data:

  • Run model on all other observations
  • Use model to predict value for observation

This means that a model is fitted, and a predicted is made n times where n is the number of observations in your data.

Seu exemplo usa a (já) clássica base de dados do R, mtcars. Eu resolvi usar uma do pacote AER , a base OECDGrowth.

Um ponto importante é que a amostra precisa ser formatada para uso com os comandos do exemplo original. Veremos isso a seguir.

A regressão que vou replicar é uma equação de crescimento que Kleiber e Zeileis estimam na seção 4.4 de seu ótimo livro (obviamente, o Applied Econometrics with R), quando discutem resistant regressions.  Trata-se desta equação:

solow_lm <- lm(log(gdp85/gdp60) ~ log(gdp60) + log(invest) + log(popgrowth + .05), data = OECDGrowth)

Vejamos o que fiz.

Primeiro criei as variáveis já em logaritmos para evitar possíveis problemas com o uso de fórmulas dentro das funções utilizadas no exemplo (por via das dúvidas, prefiro me precaver). Eis o código.

library(AER)

data("OECDGrowth")
# vamos criar a base das variaveis usadas na solow_lm
solow_lm &lt<- lm(log(gdp85/gdp60) ~ log(gdp60) +
                    + log(invest) + log(popgrowth + .05), data = OECDGrowth)

gdp85<-OECDGrowth$gdp85
gdp60<-OECDGrowth$gdp60
invest<-OECDGrowth$invest
popgrowth<-OECDGrowth$popgrowth
lgdp85_60<-log(gdp85/gdp60)
lgdp60<-log(gdp60)
linves<-log(invest)
lpopgrowthplus005<-log(popgrowth+0.05)
data1<-data.frame(lgdp85_60,lgdp60,linvest,lpopgrowthplus005)

Pronto para a replicação? É bom ler o código original com atenção. Ao fazer isso, descobre-se facilmente que é preciso alterar o nome da base de dados original – e o da variável dependente – para não ter surpresas desagradáveis. O trecho do código abaixo pode ser comparado com o original citado lá no início deste pequeno texto.

library(pipelearner)
pl <- pipelearner(data1, lm, lgdp85_60 ~ .)
pl <- learn_cvpairs(pl, k = nrow(data1))
pl <- learn(pl)

# ou
pl <- pipelearner(data1, lm, lgdp85_60 ~ .) %>% 
  learn_cvpairs(k = nrow(data1)) %>% 
  learn()

O que está acontecendo? Não é esta a minha preocupação, neste estágio do exercício (como disse no início do texto). Mas você pode descobrir um pouco mais aqui. Basicamente, estamos usando um algoritmo de machine learning na regressão de crescimento original.

Os passos seguintes não necessitam alterações. Para não ficar tão igual, vou traduzir os elementos textuais do gráfico (ok, eu sei que é preguiça).

library(tidyverse)

# Extrai os valores observados e previstos da variavel dependente para cada observação
pl <- pl %>% 
  mutate(true = map2_dbl(test, target, ~as.data.frame(.x)[[.y]]),
         predicted = map2_dbl(fit, test, predict)) 

# Resumo dos resultados
results <- pl %>% 
  summarise(
    sse = sum((predicted - true)^2),
    sst = sum(true^2)
  ) %>% 
  mutate(r_squared = 1 - sse / sst)

results

pl %>% 
  ggplot(aes(true, predicted)) +
  geom_point(size = 2) +
  geom_abline(intercept = 0, slope = 1, linetype = 2)  +
  theme_minimal() +
  labs(x = "Valor observado", y = "Valor previsto") +
  ggtitle("Valores Observados contra previstos baseados\n em leave-one-one cross validation")

Agora, vejamos o gráfico.

crossvalidation_2017

Comparando com o gráfico do autor para seu exercício em mtcars confesso que não consigo dizer se há alguma superioridade no ajuste. Isto é bastante razoável, dado que o R² obtido por ele (0.95) não é tão distante do que o obtido aqui (0.94).

Ok, concluímos o nosso objetivo de refazer os passos do exercício original com outra base de dados. Vale ressaltar que o pipelearner tem funções mais interessantes e talvez você queira estudar outros exemplos criados pelo autor do pacote (aliás, caso você queira realmente entender o que foi feito aqui, é bom dar uma olhada lá). Aliás, após o sucesso de um exercício como este, a gente ganha confiança para, justamente, seguir em busca de maior compreensão, não é?

Até a próxima.

Dica R do Dia – Penn World Tables 9.x

O prof. Zeileis tem feito ótimos pacotes para o R. O último é o que contém os dados da última versão da Penn World Tables (ele tem também pacotes para as versões anteriores). Claro, deve-se ter muito cuidado com esta base de dados.

Ela já foi mais simples. Hoje, está longe de ter seu uso feito sem um bocado de restrições. Leia, por exemplo, a documentação que se encontra aqui.

Um exemplo bem simples de como aprender a analisar dados em R

É um blog cheio de propaganda, mas vale a pena. Didático e simples. A dica é: copie e cole os comandos e replique lendo o texto. Depois, pense no que você pode fazer com isso (ou seja, com os dados que te interessam). Pense também nas limitações do exemplo (o que poderia ser feito para minimizar seu trabalho).

R é tudo, heim? ^_^

Momento R do Dia – Gastos per capita e Gastos em Percentual do PIB

O livro de finanças públicas de Joseph Stiglitz, em sua segunda edição, apresenta um gráfico interessante em seu capítulo 2. Ele diz respeito à correlação entre o percentual do gasto do governo no PIB e a renda per capita para um conjunto de doze países para 1982. Usando os dados do Banco Mundial para o ano de 2014, resolvi ampliar a base de dados.

Primeiramente, tive que excluir alguns subconjuntos de países que o Banco Mundial usa para manter uma amostra apenas de países e não de grupos como “alta renda”, “baixa renda”, etc. Em segundo lugar, para a renda per capita, utilizei a renda por população empregada calculada na PPP de 2011.

Claro, usei o R.

# reproduz a figura 2.8 do livro de Stiglitz de setor público.
data <- read.delim("clipboard")
data<-as.data.frame(data)
data
library(ggplot2)
ggplot(data, aes(x=GDP_per_person_employed_2011_PPP, y=General_gov_final_cons_expenditure_perc_of_GDP))+
  stat_smooth(method="auto",level=0.99)+
  geom_text(aes(label=Code), size=3)+
  ggtitle("Breve correlação - 2014") +
  xlab("GDP per person employed - 2011 PPP") +
  ylab("General gov. final cons. exp. in percent of GDP") 

ggplot(data, aes(x=log(GDP_per_person_employed_2011_PPP), y=log(General_gov_final_cons_expenditure_perc_of_GDP)))+
  stat_smooth(method="auto",level=0.99)+
  geom_text(aes(label=Code), size=3)+
  ggtitle("Breve correlação - 2014 (ambas em log)") +
  xlab("GDP per person employed - 2011 PPP") +
  ylab("General gov. final cons. exp. in percent of GDP")

O gráfico obtido é pouco elucidativo (no código acima, repare que copiei os dados do Excel e colei no R com o bom e velho “clipboard”).

stiglitz1
Assim, passei ambas para a escala logarítmica para, ao menos, visualizar melhor os países que apareciam bem misturados no início do gráfico.

stiglitz2
Os códigos dos países e as variáveis você encontra aqui. Há, certamente, muita coisa para se incomodar nestes gráficos.

Primeiro, há a questão de se entender a proxy de renda per capita utilizada pois países como Brunei (BRN) aparecem em uma – intrigante – alta posição no que diz respeito à esta variável, por exemplo.

Em segundo lugar, não há uma relação clara entre as duas variáveis e talvez isso possa querer dizer alguma coisa. Volto a isto adiante.

Em terceiro lugar, claro, existe a questão de se olhar para um fenômeno apenas em um período de tempo. Um painel seria interessatne, mas a visualização ficaria mais complicada.

Em quarto lugar, há quem goste de brigar por conta do uso de um “gasto de X em percentual do PIB” em uma suposta contraposição ao “gasto de X per capita”. Os gráficos acima não fazem um contraponto a este suposto problema, é verdade. Mas o interessante é que ele faz uma correlação que deve incomodar alguns: compara o gasto governamental em percentual do PIB com o PIB per capita para diversos países (alguma discussão interessante sobre isto aqui).

A idéia era apenas a de ver o que Stiglitz parecia apontar: países com renda per capita mais alta tenderiam a apresentar um gasto do governo mais alto (em percentual do PIB). Bem, o próprio Stiglitz chama a atenção para a cautela de se comparar países distintos e acaba limitando seu gráfico para uma amostra de países desenvolvidos. Em princípio, eu esperava obter grupos de países separados nos gráficos acima. Contudo, não encontrei nada nem perto disto. Talvez os problemas acima sejam importantes para se prosseguir com o tema, mas a idéia era apenas a de fazer um exercício em R.

Até a próxima!

Econometria aplicada…em português

Já se vão anos desde minha apostila sobre o Eviews que, segundo me dizem, já ajudou muita gente por este país. Recentemente publiquei a apostila para o R (está aqui) e o Vitor Wilher publicou a dele (que não é gratuita e está aqui) e não é que tenho o prazer de saber, por meio do prof. Gabrielito, da FURG, que o pessoal da FEE (órgão estadual do RS) acaba de publicar uma apostila gratuita para o Eviews?

Outro dia o Leo Monasterio me instigou a publicar um post retrospectivo sobre a blogosfera no debate. Não sei se tenho ânimo para tanto, mas eis aqui uma evidência de que a internet e a blogosfera mudaram, inclusive, a forma como aprendemos e também como ensinamos. Apostilas de pacotes econométricos em português a um preço baixo (ou mesmo gratuitas)? Foram anos vendo materiais em inglês e, aos poucos, vemos a mudança aparecer por aqui.

Que ótima notícia!

Nova Apostila do R para Econometria Aplicada !

apostilar_capaTenho a satisfação de anunciar a disponibilização da nossa apostila de econometria aplicada em R. Sim, antes que você pergunte, ela é gratuita.

Ela não está ainda do jeito que eu gostaria. Ao longo dos anos, vários co-autores foram convidados mas, no final, ficamos apenas eu e o prof. Rodrigo. Portanto, embora versões anteriores tenham circulado com nomes de outros colegas, esta versão  – totalmente repaginada para uso em RStudio leva apenas nossos nomes, ok?

Confesso que tentei algumas editoras, mas todas elas foram taxativas: não insista, não queremos saber de um livro de apoio para econometria aplicada em programas gratuitos. Diante desta resposta, pensei em arriscar um ebook, mas é tanto trabalho com formatações chatas que, sozinho, eu não conseguiria.

Além disso, embora eu não acredite que exista uma horda de alunos famintos por econometria aplicada no Brasil (muitas vezes eles sequer lêem uma PEC de duas páginas antes de palpitarem como especialistas em finanças públicas…), acredito que há sempre um aluno aqui ou ali, sem muitos recursos para aprender econometria e que, infelizmente, ainda não venceu a barreira do inglês. É para ele que dedico a apostila.

Ao longo dos anos, muitos amigos e colegas me incentivaram a usar o R. São tantos que não vou nomeá-los para não cometer injustiças. De qualquer forma, que fiquem registrados meus agradecimentos a eles.

O endereço da página oficial da apostila é este.