Que variáveis incluir na regressão? Ou, por que grande parte dos trabalhos aplicados está errada.


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 Pearlsem 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): dagitty-model 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.

Sobre a acurácia das variáveis econômicas – Dilbert Style.


Dilbert, resumindo Os Números (não) mentem em uma tirinha.
dilbert

Ok, tenho resgatado uns cartoons velhos. Mas valem a pena, não valem!?

 

Mapas de roubos em Brasília?


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.

crimes_geral

 

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.

hora

 

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.

semana

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)

 

Quanto mais tempo sem alugar, maior a variação do preço do aluguel? E mais um mapa.


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:

variacao

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.

mapa_aluguel

Cartoon: não olhe somente a média!


 

Uma antiga, mas excelente, do SMBC:
fotoPara ver outras relacionadas a economia ou estatística, clique aqui.

 

Distribuindo pacotes no R. Qual o alcance?


Faz mais ou menos 1 mês que o pacote benford.analysis 0.1 foi disponibilizado no CRAN. 

Achei que valeria o esforço adicional de criar o pacote por alguns motivos e, entre eles, dois se destacam: (i) pacotes deixam os arquivos fontes mais estruturados, facilitam o uso das funções, forçam a criar uma documentação e passam por uma bateria de sanity tests que ajudam a criar boas práticas de programação; (ii) pacotes tornam o compartilhamento do código extremamente simples, ainda mais se o pacote estiver no CRAN, pois, basta rodar

 install.packages("benford.analysis") 

e qualquer pessoa de qualquer lugar do mundo terá o pacote instalado em sua máquina.

Sobre este último ponto, infelizmente, não é possível ter dados de download do CRAN de uma maneira consolidada, pois há diversos espelhos do site ao redor do mundo e nem todos guardam informações de acesso. Entretanto, o CRAN do RStudio faz esse registro. Assim resolvi baixar os dados de lá e ver se alguma outra alma além de mim baixou o benford.analysis.

Sinceramente, eu achei que encontraria uns 10 ou 11 downloads registrados – no máximo -, pois estamos com dados de apenas um espelho do CRAN e estamos falando de um pacote simples e relativamente desconhecido. Ocorre que neste 1 mês de existência o benford.analysis foi baixado 190 vezes em mais de 40 países diferentes considerando apenas o espelho do RStudio.  Um número pequeno quando comparado com pacotes como ggplot2 (que deve estar em virtualmente quase toda máquina de usuário do R), mas, ainda assim, grande o suficiente para me surpreender!

E também para me preocupar. Nesse meio tempo encontrei dois pequenos bugs no pacote. E se antes achava que não deveria ter pressa para corrigi-los, agora esperem uma atualização em breve (mas, claro, depois do carnaval)!

Mapa de aluguel em Brasília (Plano Piloto)


Em post anterior fizemos uma breve análise dos dados de aluguel no plano piloto.

Agora, que tal navegar por todos imóveis em um mapa da cidade, vendo a localização, tamanho, número de quartos e valor do aluguel? Clique aqui ou na mapa abaixo para navegar.

Atenção,  ainda é um protótipo!

Se o mapa não aparecer na sua tela, provavelmente o seu navegador bloqueou a execução do javaScript. Procure por um cadeado no navegador (canto superior direito ou esquerdo, geralmente) e autorize o carregamento do site.

Captura de Tela 2014-02-23 às 21.13.59

PS: agora já estamos coletando diariamente e automaticamente preços online de imóveis dos principais sites e das principais capitais do país. Ainda estamos testando métodos de análise e visualização.

Existe um pacote para isso.


Por que eu uso R, em 30 segundos:

O que o Facebook diz sobre o seu relacionamento?


O time de análise de dados do Facebook fez uma série de 6 posts sobre o valentine’s day (dia dos namorados) nos Estados Unidos.

Recomendo fortemente a leitura de todos. O posts tratam dos seguintes temas:

  • O primeiro post trata de amor e religião e constata que há poucos casais de religiões diferentes, mesmo em países com alta diversidade religiosa.
  • O segundo post trata da diferença de idade entre casais. Na média, homens são mais de dois anos mais velhos do que suas  parceiras.
  • O terceiro post trata da duração dos relacionamentos. Um dos resultados: quanto mais tempo de relacionamento, menor a chance de o casal se separar.
  • O quarto post trata das “melhores” cidades para os solteiros (como são cidades dos EUA, provavelmente não interessará muito os leitores deste blog).
  • O quinto post trata da mudança de comportamento dos casais antes e depois do relacionamento. Esse é um dos mais bacanas. Para quem quiser ler algo em português, a Folha fez uma matéria. Vale reproduzir um gráfico, relacionando a quantidade de posts com palavras positivas e os dias antes/após o início do namoro:

1898250_10152219519288415_127545461_n

Os dados confirmam aquilo que você já percebia: casais recém formados postam sobre unicórnios vomitando arco-iris e o efeito pode durar muito, muito tempo (destaque para o gráfico feito com ggplot2).

  • Por fim, o último post trata do que acontece após o término do relacionamento. As interações, principalmente de apoio dos amigos, aumentam bastante.

O Facebook é, muito provavelmente, a organização com a maior base de dados sobre informações pessoais do mundo. O potencial disso é inimaginável. No final do ano passado, eles contrataram o professor da NYU Yann LeCun para liderar o departamento de inteligência articial da empresa – parece que ainda há muita coisa interessante por esperar.

Mais sobre análise de dados do Facebook neste blog, aqui (analise seus próprios dados) e aqui (descubra características  da pessoa – como a orientação sexual – com base no que ela curte).

Analise a pesquisa mensal de emprego com R


Mais um post no excelente analyze survey data for free, agora com a PME.

Parabéns ao Anthony Damico e ao Djalma Pessoa pela iniciativa!