Parece que Sérgio Almeida e Mauro Rodrigues, do Economistas X, estão com um paper bacana no forno: coletar os próprios dados não é tarefa fácil, confiram no post algumas das agruras pelas quais os dois passaram!
PS: sou partidário da idéia de que a coleta de dados interessantes vale um paper por si só. E, claro, que os dados sejam abertos ao público!
Depois que o Nate Silver saiu do NYT, o jornal montou um novo time de previsão para eleições. O novo modelo, denominado LEO, está no ar e com uma descrição bem amigável de seu funcionamento. A parte gráfica dos resultados também está bacana, inclusive com uma seção comparando o modelo do NYT com seus concorrentes (vale lembrar que o site conta com a ajuda do Michael Bostock, um dos caras que desenvolveu o fantástico D3 para JavaScript)
Mas, tem algo ainda melhor. O NYT liberou os dados e os códigos do modelo no github. E o modelo é em R. Ou seja, agora, para replicar e adaptar o modelo à realidade brasileira só faltam duas coisas: tempo e vontade.
Suponha que você queira medir o efeito de X em Y (isto é, o quanto uma variação de X afeta Y – uma relação causal) e que você tenha mais duas variáveis que podem ser incluídas como controle, Z1 e Z2. Suponha ainda que você saiba que o modelo é linear, isto é, não há nenhuma incerteza com relação à especificação. Quais variáveis você incluiria no seu modelo?
Hoje, provavelmente você diria o seguinte: depende da significância! São Z1 e Z2 significantes? Se sim, eles devem ser incluídos.
Vejamos um exemplo de uma simulação. O código em R está ao final do post.
Vamos rodar as três regressões: uma só com X, outra incluindo Z1 e por fim uma com todos os controles. Os resultados foram os seguintes:
Equação 1: Y = -10 + 43X ***
Equação 2: Y = -7 + 13X * + 107Z1 ***
Equação 3: Y = -5 – 9X * + 46Z1 *** + 37Z2 ***
Pelos resultados, tanto Z1 quanto Z2 são significantes, então preferimos a equação 3. Concluímos que, na média, uma variação de 1 unidade de X reduz Y em 9 unidades. Certo?
***
Errado. O modelo ideal neste caso seria a equação 2. O efeito real de X sobre Y é de 10 (veja que valor estimado foi 13, bem próximo). O problema aqui é que a significância estatística não vai te responder sobre a pertinência de incluir ou não uma variável para estimar o efeito de X sobre Y. Infelizmente, não há almoço grátis. Como diz Judea Pearl, sem saber a estrutura do problema, não é possível determinar quais variáveis devem ser incluídas. Agora pense. Como é a lógica de trabalho dos artigos aplicados hoje? *** A simulação A nossa simulação tem a seguinte estrutura (U1 e U2 dizem respeito a duas variáveis não observadas, só observamos Y, X, Z1 e Z2): 
O código em R para gerar os resultados é:
gen_data <- function(N=200,s=2,beta1=10, beta2=100){
Z1 <- rnorm(N,0,s)
U2 <- rnorm(N,0,s) + Z1
U1 <- rnorm(N,0,s) + Z1
Z2 <- rnorm(N,0,s) + U2 + U1
X <- rnorm(N,0,s) + U1
Y <- rnorm(N,0,s) + beta1*X + beta2*U2
data.frame(Y,X,Z1,Z2)
}
set.seed(100)
data <- gen_data()
summary(lm(Y~X, data))
summary(lm(Y~X + Z1, data))
summary(lm(Y~X + Z1 + Z2, data))
Você pode brincar mais com o paradoxo de Simpson aqui; e o gráfico você pode fazer aqui.
O famoso site de programação StackOverflow (SO) ganhou uma versão tupiniquim.
O SO é um excelente site de perguntas e respostas. Seu diferencial é ser direto: as perguntas têm que ser bem definidas e as respostas têm de resolver diretamente o problema. Quer saber, por exemplo, como agregar uma base de dados no R? Pergunte lá e surgirão várias respostas diferentes de como se fazer isso.
Ainda há poucos usuários ativos no R do SO em português. Mas estamos fazendo um esforço para popular o site com perguntas e respostas. Você pode fazer perguntas sobre problemas que está enfrentando atualmente ou, inclusive, registrar perguntas e respostas que você já sabe, como, por exemplo, aqui (gráfico em 3d), aqui (barplot) ou aqui (contar ocorrências em um vetor) – alguém certamente passará pela mesma dificuldade e a solução que você encontrou para o problema pode ser útil. Ou, ainda, outro usuário pode ter uma solução mais interessante do que a que você propôs. De uma olhada nas perguntas que já foram feitas sobre R aqui.
Se você usa R (Python), cadaste-se no StackOverflow em Português e ajude o site a crescer! Podemos torná-lo um ótimo ambiente para a comunidade brasileira de R, tal como é hoje o SO em inglês.
Peter Singer, Paul Bloom e Dan Ariely irão discutir agora, ao vivo, suas pesquisas sobre desonestidade, moralidade e ética.
Dilbert, resumindo Os Números (não) mentem
Ok, tenho resgatado uns cartoons velhos. Mas valem a pena, não valem!?
Recentemente conheci um site com uma iniciativa bem bacana chamado Onde Fui Roubado. Lá qualquer pessoa pode reportar um crime especificando local, hora, objetos roubados e inclusive fornecer um relato. Há mais de 16 mil registros para várias cidades do país, e resolvi fazer um webscraping para ver como são estes dados.
Especificamente para Brasília, infelizmente, existem apenas cerca de 200 registros. A maioria na Asa Sul, Asa Norte e Sudoeste, com mais de 100. A ideia aqui será montar um mapa de calor, ou de densidade, dos roubos no Plano Piloto.
Temos, entretanto, dois problemas que valem ser ressaltados: (i) a amostra é pequena; e, (ii) possivelmente viesada. Isto é, como o site ainda não parece ser muito conhecido, não necessariamente o público que está informando é representativo da população do local. Ainda assim, tendo em mente essas ressalvas, vamos brincar um pouco com a visualização dos dados!
Primeiro, vejamos um mapa com todos os casos – note que, quanto mais vermelho, maior a concentração de roubos reportados na região. A maior parte dos registros foram na Asa Sul e Asa Norte. Na Asa Norte, em especial, a região próxima à UnB tem destaque. Lembre que talvez isto seja decorrência, por exemplo, de pessoas mais jovens conhecerem o site e reportarem mais casos.
Vamos dividir agora o mapa por horário do roubo, entre manhã, tarde, noite e madrugada. A maior parte dos roubos registrados ocorreu durante a noite, com focos na Asa Norte e início da Asa Sul.
Vejamos, ainda, uma divisão por dias da semana. De maneira consistente com os mapas anteriores, aparece um foco nas sextas, na região próxima à UnB.
Poderíamos fazer um mapa cruzando dias da semana e hora, mas temos poucos dados para isso. A ideia aqui é mostrar como podem ser poderosas essas visualizações! Se a Secretaria de Segurança Pública liberar os microdados dos BO’s (se alguém tiver estes dados, por favor, entre em contato), seria possível montar mapas bem acurados. E imagine cruzá-los com as informações de imóveis – poderíamos medir o impacto da criminalidade nos preços imobiliários.
Por fim, reforço a divulgação do Onde Fui Roubado, é uma iniciativa louvável!
***
A quem interessar, seguem os códigos para a construção dos mapas. Os dados podem ser baixados aqui.
library(ggmap)
library(dplyr)
### carrega dados
dados <- readRDS("roubo2.rds")
### Pega mapa de Brasília
q<-qmap("estadio mane garrincha, Brasilia", zoom=13, color="bw")
### transformando data em POSIXlt e extraindo hora
dados$hora <- as.POSIXlt(dados$data)$hour
### selecionando a base de dados do plano piloto, criando semanas e horários
bsb <- filter(na.omit(dados), cidade=="Brasília/DF",
lon > -47.95218, lon < -47.84232,
lat > -15.83679, lat < -15.73107)%.%
mutate(semana = weekdays(data),
hora = cut(hora,
breaks=c(-1,6,12,18,25),
labels=c("Madrugada", "Manhã", "Tarde", "Noite")))
### reordenando os dias da semana
bsb$semana <- factor(bsb$semana, levels = c("segunda-feira", "terça-feira",
"quarta-feira", "quinta-feira",
"sexta-feira", "sábado", "domingo"))
### estrutura básica do gráfico
map <- q + stat_density2d(
aes(x = lon, y = lat, fill = ..level.., alpha = ..level..),
size = 2, bins = 4, data = bsb,
geom = "polygon")
### mapa geral
map + scale_fill_gradient(low = "black", high = "red", guide=FALSE)+
scale_alpha(guide=FALSE)
### mapa por dia da semana
map+scale_fill_gradient(low = "black", high = "red", guide=FALSE)+
facet_wrap(~ semana)+scale_alpha(guide=FALSE)
### mapa por horário
map+scale_fill_gradient(low = "black", high = "red", guide=FALSE)+
facet_wrap(~ hora) + scale_alpha(guide=FALSE)
Hoje, com 30 dias de coleta e mais de 60.000 observações de preços de aluguéis de Brasília, resolvi explorar um pouco os dados.
Será que, como nos diz a intuição, quanto mais tempo o imóvel passou ofertado, maiores as reduções observadas do preço do aluguel?
Vejamos com o gráfico abaixo.
No eixo x temos quantos dias o imóvel ficou ofertado durante os 30 dias de coleta e, no eixo y, a soma da variação percentual do valor do aluguel no período:
Parece que os preços de oferta são relativamente rígidos, mas depois de 20 dias sem alugar começam a ceder. Vejamos se o padrão se mantém e como isso se comporta mais para frente!
PS: como muitos imóveis podem estar no mesmo ponto – por exemplo, a maioria tem variação zero no preço – o gráfico pode dar a impressão de que há poucas observações com poucos dias de anúncio. Na verdade há muitos pontos ali, o problema é que eles estão um em cima do outro, o que chamamos de overplotting. Uma outra forma de visualizar a distribuição tentando suavizar o overplotting é com um pouco de jitter (desvios aleatórios na posição dos pontos), você pode ver o mesmo gráfico com jitter aqui. Outra coisa que vale a pena ser novamente ressaltada é que o gráfico não é uma série temporal! Ele relaciona a quantidade de dias que um anúncio ficou no ar com a variação percentual do preço deste anúncio.
***
Resolvi também testar outra forma de visualização espacial dos dados. No mapa abaixo, quanto mais vermelho, mais caro o aluguel e, quanto maior a bola, maior o apartamento (em metros quadrados). Os dados são de hoje.
Análise Real 