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())

Seminário – Ciência de Dados e Sociedade, dia 15 de Junho às 19h, no Auditório do Instituto de Ciência Política da UnB


Inscrições e mais informações aqui. Estaremos no seminário eu, Daniel Marcelino e Rommel Carvalho.