Regressão robusta, erro de medida e preços de imóveis


Um amigo estava tendo problemas ao analisar sua base de dados e pediu ajuda — ao olhar alguns gráficos o problema parecia claro: erro de medida. Resolvi revisitar um post antigo e falar um pouco mais sobre como poucas observações influentes podem afetar sua análise e como métodos robustos podem te dar uma dica se isso está acontecendo.

Voltemos, então, ao nosso exemplo de uma base de dados de venda de imóveis online:

arquivo <- url("https://dl.dropboxusercontent.com/u/44201187/dados/vendas.rds")
con <- gzcon(arquivo)
vendas <- readRDS(con)
close(con)

Suponha que você esteja interessado na relação entre preço e tamanho do imóvel. Basta um gráfico para perceber que a base contém alguns dados muito corrompidos:

with(vendas, plot(preco ~ m2))

unnamed-chunk-15-1

Mas, não são muitos pontos. Nossa base tem mais de 25 mil observações, será que apenas essas poucas observações corrompidas podem alterar tanto assim nossa análise? Sim. Se você rodar uma regressão simples, ficará desapontado:

summary(m1 <- lm(preco ~ m2, data = vendas))
##
## Call:
## lm(formula = preco ~ m2, data = vendas)
##
## Residuals:
##       Min        1Q    Median        3Q       Max
##  -6746423   -937172   -527498     99957 993612610
##
## Coefficients:
##                Estimate  Std. Error t value             Pr(>|t|)
## (Intercept) 1386226.833   18826.675  73.631 < 0.0000000000000002 ***
## m2               18.172       3.189   5.699         0.0000000121 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9489000 on 254761 degrees of freedom
## Multiple R-squared:  0.0001275,  Adjusted R-squared:  0.0001235
## F-statistic: 32.48 on 1 and 254761 DF,  p-value: 0.00000001208

A regressão está sugerindo que cada metro quadrado extra no imóvel corresponde, em média, a um aumento de apenas 18 reais em seu preço! Como vimos no caso do post anterior, limpar um percentual bem pequeno da base é suficiente para estimar algo que faça sentido.

Mas, suponha que você não tenha noção de quais sejam os outliers da base e também que, por alguma razão, você não saiba que 18 reais o metro quadrado é um número completamente absurdo a priori. O que fazer? (Vale fazer um parêntese aqui – se você está analisando um problema em que você não tem o mínimo de conhecimento substantivo, não sabe julgar sequer se 18 é um número grande ou pequeno, plausível ou não, isso por si só é um sinal de alerta, mas prossigamos de qualquer forma!)

Um hábito que vale a pena você incluir no seu dia-a-dia é rodar regressões resistentes/robustas, que buscam levar em conta a possibilidade de uma grande parcela dos dados estar corrompida.

Vejamos o que ocorre no nosso exemplo de dados online:

library(robust)
summary(m2 <- lmRob(preco ~ m2, data = vendas))
##
## Call:
## lmRob(formula = preco ~ m2, data = vendas)
##
## Residuals:
##         Min          1Q      Median          3Q         Max
## -3683781389     -202332      -23119       64600   994411077
##
## Coefficients:
##               Estimate Std. Error t value            Pr(>|t|)
## (Intercept) -15926.247    589.410  -27.02 <0.0000000000000002 ***
## m2            9450.762      5.611 1684.32 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 171800 on 254761 degrees of freedom
## Multiple R-Squared: 0.4806
##
## Test for Bias:
##             statistic p-value
## M-estimate     502.61       0
## LS-estimate     86.91       0

Agora cada metro quadrado correponde a um aumento de R$9.450,00 no preço do imóvel! A mensagem aqui extrapola dados online, que são notórios por terem observações com erros de várias ordens de magnitude. Praticamente toda base de dados que você usa está sujeita a isso, mesmo de fontes oficiais. No post anterior vimos um exemplo em que pesquisadores não desconfiaram de uma queda de 36% (!!!) do PIB na Tanzânia.

Por fim, vale fazer a ressalva de sempre: entender o que está acontencedo nos seus dados — por que os valores são diferentes e a razão de existir de alguns outliers  — é fundamental. Dependendo do tipo de problema, os outliers podem não ser erros de medida, e você não quer simplesmente ignorar sua influência. Na verdade, há casos em que outliers podem ser a parte mais interessante da história.

Anúncios

Google Trends no R


O pacote gtrendsR está passando por uma reformulação e parece que vai ficar ainda mais fácil analisar dados do Google Trends no R. A nova versão ainda não está no CRAN, mas já pode ser testada pelo github. Para instalar:

install.packages("devtools")
devtools::install_github('PMassicotte/gtrendsR', ref = 'new-api')

A grande novidade dessa versão é que não será mais preciso fazer login no google trends para ter acesso. Para brasileiros, outra novidade é que os bugs com problema de encoding parecem estar diminuindo.

Vejamos um exemplo simples, pegando dados das buscas pelos nomes dos candidatos nas eleições de 2014 no Brasil:

library(gtrendsR)
eleicoes2014 <- gtrends(c("Dilma Rousseff", "Aécio Neves", "Marina Silva"), geo = c("BR"), time = "2014-01-01 2014-12-31")
plot(eleicoes2014)

rplot01

Para ilustrar novamente, vejamos um exemplo mais recente — as buscas pelos nomes dos candidatos das eleições norte-americanas:

USelections2016 <- gtrends(c("Donald Trump", "Hillary Clinton"), geo = c("US"), time = "2016-01-01 2016-12-31")
plot(USelections2016)

rplot

Personalizando seu gráfico do ggplot2 – Exports and Imports, William Playfair


O ggplot2 é muito bom para explorar visualmente, de forma dinâmica, sua base de dados.  Mas às vezes queremos editar cada detalhe do gráfico para uma publicação, é possível fazer isso?

Como, por exemplo, reproduzir o famoso gráfico de exportações e importações do William Playfair?

Playfair-bivariate-area-chart

Hoje resolvi testar o quão difícil seria gerar uma imagem parecida e, brincando um pouco com os parâmetros, cheguei na figura abaixo. É um pouco trabalhoso – pois temos que colocar cada texto separadamente – mas não é difícil, nem tão demorado assim.

playfair

Se você tiver um pouco mais de paciência para ajustar detalhes talvez consiga tornar a reprodução ainda mais fiel. E, caso o faça, favor compartilhar o código com todos por aqui!

***

Segue abaixo o código para gerar o gráfico acima. Os dados bem como o próprio código também estão no github.

 

# load packages -----------------------------------------------------------
library(reshape2)
library(ggplot2)

# prepare data for ggplot2 ------------------------------------------------
## data extracted from https://plot.ly/~MattSundquist/2404/exports-and-imports-to-and-from-denmark-norway-from-1700-to-1780/#plot
playfair <- readRDS("william_playfair.rds")

## create min for geom_ribbon
playfair$min <- with(playfair, pmin(exp, imp))
year <- playfair$year

## melt data
molten_data <- melt(playfair, id.vars = c("year", "min"))

# ggplot2 -----------------------------------------------------------------
ggplot(molten_data, aes(x = year, y = value)) +
geom_line(aes(col = variable), size = 1) +
geom_ribbon(aes(ymin = min, ymax = value, fill = variable), alpha = 0.3) +
scale_color_manual(values = c("darkred", "gold3"), guide = F) +
scale_fill_manual(values = c("#90752d", "#BB5766"), guide = F) +
theme_bw() +
annotate("text", x = year[5], y = 100000, label = "Line", angle = 25, size = 3, family = "Garamond") +
annotate("text", x = year[6] - 100, y = 104000, label = "of", angle = 0, size = 3, family = "Garamond") +
annotate("text", x = year[7], y = 101000, label = "Imports", angle = 340, size = 3, family = "Garamond") +
annotate("text", x = year[5] + 400, y = 73000, label = "Line", angle = 345, size = 3, family = "Garamond") +
annotate("text", x = year[6], y = 70000, label = "of", angle = 330, size = 3, family = "Garamond") +
annotate("text", x = year[7] - 200, y = 64000, label = "Exports", angle = 335, size = 3, family = "Garamond") +
annotate("text", x = year[8], y = 83000, label = "italic('BALANCE AGAINST')", angle = 0, family = "Garamond", parse = TRUE) +
annotate("text", x = year[16] + 400, y = 110000, label = "italic('BALANCE in\nFAVOUR of\nENGLAND')", angle = 0, family = "Garamond", parse = TRUE) +
annotate("text", x = year[16], y = 82000, label = "Imports", angle = 30, size = 3, family = "Garamond") +
annotate("text", x = year[14] + 200, y = 131000, label = "Exports", angle = 65, size = 3, family = "Garamond") +
ggtitle("Exports and Imports to and from DENMARK & NORWAY from 1700 to 1780") +
scale_x_date(breaks = seq(year[1], year[18], by = "10 years"),
labels = format(seq(year[1], year[18], by = "10 years"), "%Y")) +
scale_y_continuous(breaks = seq(0, 190e3, by = 10e3),
labels = seq(0, 190, by = 10)) +
theme(title = element_text(size = 8, face = 'bold', family = "Garamond"),
axis.title = element_blank(),
axis.text = element_text(family = "Garamond"),
panel.grid.minor = element_blank())

Analisando seu histórico de pesquisas do Google


Hoje descobri que é possível fazer o download de todo seu histórico de buscas no Google. TODO seu histórico de TUDO o que você busca no Google. Já que a opção está disponível, por que não dar uma olhada nos dados?

Por alguma razão meu histórico só vai até 2014 — acredito que tenha deletado o histórico anterior — então no meu caso temos apenas dois anos de dados para analisar (não vou considerar 2016 aqui pois o ano ainda não terminou). Além disso, esses dados certamente não contemplam tudo o que pesquisei na internet neste período, porque: (i) além do Google eu uso o DuckDuckGo; e, (ii) muitas vezes não estou logado quando faço pesquisas no próprio Google.

Feitas as ressalvas anteriores, a primeira coisa que tentei montar foi uma nuvem com as palavras mais utilizadas nas buscas. Em 2014 e 2015, segundo o registro do google, fiz aproximadamente 19 mil buscas, utilizando aproximadamente 69 mil palavras-chave. Após remover algumas “stopwords” em inglês e português — isto é, preposições, artigos etc — fiz uma nuvem com aquelas palavras que representam cerca de 20% da frequência total, e o resultado foi o seguinte:

Não tem muita surpresa aí. Previsivelmente, “R” foi a palavra chave mais utilizada, seguida de “package”, “statistics”, “Mac”, “Data”, “Los Angeles”, “UCLA” entre outras.

Após verificar as palavras mais utilizadas, procurei ver se encontrava alguns padrões nos meus hábitos de busca. Primeiramente, calculei a média de buscas por dia da semana. Nesses dois anos, as buscas parecem ter alcançado seu pico de segunda a quarta:

por_semana

Em seguida calculei a média por hora. Tirando a madrugada e o início da manhã, não parece existir diferença significativa entre os horários.  Há, contudo, um problema com essa informação: elas estão no horário brasileiro. Como estive fora do país em certas datas, isso distorce o horário original de algumas pesquisas — e ainda não descobri como consertar esse problema de maneira automática.

por_hora

Essa questão das viagens para fora do país suscitou outra pergunta: o total de buscas no Google Maps altera quando estou viajando? A princípio, diria que sim, e é isso o que o gráfico a seguir mostra, com algumas viagens destacadas:

Isto é, pelo menos neste caso, é muito fácil identificar viagens utilizando apenas a série histórica do total de buscas do Google Maps.

Para finalizar, montei um gráfico com a média de pesquisas por hora, separados por dia da semana e ano, mas não parece ter havido mudança relevante entre os padrões de 2014 e 2015.

hora_semana_ano
Quer analisar seus dados também?

Para fazer o download dos dados, basta seguir essas instruções. Os dados virão em um arquivo zip com vários arquivos no formato JSON. Para tratá-los, você pode se basear no script de R que coloquei aqui.

PS: É um pouco assustador perceber que, com análises bastante simples de dados de busca, já é possível inferir bastante coisa sobre os hábitos de uma pessoa.

Replicação em economia


John Cochrane soltou um post bacana sobre replicação em economia. Vale a pena conferir.

On replication in economics. Just in time for bar-room discussions at the annual meetings.

“I have a truly marvelous demonstration of this proposition which this margin is too narrow to contain.” -Fermat

“I have a truly marvelous regression result, but I can’t show you the data and won’t even show you the computer program that produced the result” – Typical paper in economics and finance.

The problem 

Science demands transparency. Yet much research in economics and finance uses secret data. The journals publish results and conclusions, but the data and sometimes even the programs are not available for review or inspection.  Replication, even just checking what the author(s) did given their data, is getting harder.

Quite often, when one digs in, empirical results are nowhere near as strong as the papers make them out to be.

I have seen many examples of these problems, in papers published in top journals. Many facts that you think are facts are not facts. Yet as more and more papers use secret data, it’s getting harder and harder to know.

The solution is pretty obvious: to be considered peer-reviewed “scientific” research, authors should post their programs and data. If the world cannot see your lab methods, you have an anecdote, an undocumented claim, you don’t have research. An empirical paper without data and programs is like a theoretical paper without proofs.

(continue lendo no blog do Cochrane)

Dados de pesquisas eleitorais de 1989 a 2015


Neale Ahmed El-Dash, do Polling Data (que já mencionamos aqui algumas vezes, como no modelo de impeachment), acabou de divulgar dados de pesquisas eleitorais brasileiras publicadas entre 1989 a 2015. Você pode acessar os dados clicando em  “Acervo/Past Elections”.