USP com a mão na massa!


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! 

 

 

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)

 

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.

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!

 

Valores de aluguel em Brasília (plano piloto)


Está pesquisando apartamento para alugar em Brasília?  Um pouco de web scrapingmanipulação e visualização de dados com os valores (de oferta) dos aluguéis de 1.030 imóveis (Asa Sul, Asa Norte e Sudoeste) do site wimoveis pode ajudar a responder algumas perguntas interessantes.

A primeira delas: qual o bairro mais caro para se alugar, hoje, no plano piloto? Esta é uma pergunta que, veremos, depende do ponto de vista. Veja a tabela abaixo (versão ampliada aqui). M2 quer dizer metro quadrado e pm2 preço por metro quadrado.

Captura de Tela 2014-01-30 às 00.08.44

Na média e mediana – em conformidade com a impressão pessoal de muitos – a Asa Sul é o bairro mais caro para se alugar dos três. Entretanto, note que isso ocorre porque há mais apartamentos maiores para aluguel na Asa Sul, e não porque o valor por metro quadrado é mais caro. Na verdade, o valor por metro quadrado, na média, é maior na Asa Norte e, na mediana, maior no Sudoeste.

Podemos agora decompor a tabela acima não somente por bairro, mas por bairro e número de quartos (versão ampliada aqui) . Na média, o bairro mais barato/caro para morar não é o mesmo a depender de quantos cômodos você quer no apartamento. E, uma curiosidade: na amostra, a média do tamanho dos apartamentos da Asa Norte, em todos os grupos de números de quartos, é menor do que a média do tamanho da Asa Sul.

tabela_wi_2014_01_24

Uma última forma de visualizar as diferenças de preços pode ser com um gráfico de densidade (versão ampliada aqui):

teste

Veja que o pico do Sudoeste (em verde) é em valores mais altos do que na Asa Norte e na Asa Sul. Entretanto, a Asa Sul tem a “cauda” mais pesada em valores próximos a R$ 5.000.

Uma outra pergunta que podemos tentar responder é a seguinte: dos anúncios que temos hoje, na média, os preços daqueles atualizados em 2014 são maiores do que aqueles cuja última atualização foi feita em dezembro de 2013? Pelo quadro abaixo (versão ampliada aqui), infelizmente, sim, e por mais ou menos RS$100,00.

Captura de Tela 2014-01-30 às 00.24.58

E como é a concentração da oferta dos anúncios por corretora? A distribuição de anúncios por imobiliária é homogênea? 

Aparentemente, não. Veja o gráfico abaixo (versão ampliada aqui).

corretoras

Enquanto algumas imobiliárias têm 30 a 40 apartamentos listados, muitas outras têm apenas 1 ou 5.

Isso quer dizer que os anúncios são concentrados? Não necessariamente. Note que apesar de a distribuição de anúncios não ser homogênea, a concorrência é bem grande, e usando como exemplo o índice de Herfindahl–Hirschman chegamos a um valor de 0.013, comumente considerado indicador de alta competitividade.

Há mais que poderíamos ver sobre aluguel. Mas deixemos para depois.  No próximo (em algum próximo) post veremos os dados de valor de venda!

PS: iremos acompanhar regularmente esses preços. E não somente para Brasília. Uma área específica do blog será criada para isso. 

benford.analysis 0.1


O pacote benford.analysis (versão 0.1) está disponível no CRAN e você já pode instalar no R com o comando:

 install.packages("benford.analysis") 

O objetivo do pacote é prover algumas funções que facilitem a validação de dados utilizando a Lei de Benford (para saber mais sobre a lei, veja aqui e aqui).

Validar como e para quê?

Um dos objetivos pode ser o auxilío na detecção de manipulações contábeis. Dados financeiros (como pagamentos) tendem a seguir a Lei de Benford e tentativas de manipulação podem acabar sendo identificadas.

Por exemplo, a lei 8.666/93 estabelece que o limite para se fazer uma licitação na modalidade convite é de R$80.000,00. Será que os valores de licitações seguiriam a Lei de Benford? Pode ser que sim. E, caso haja a tendência, uma tentativa de manipular artificialmente valores licitados para algo pouco abaixo de R$80 mil geraria um “excesso” de dígitos iniciais 7. Restaria verificar uma amostra desses registros para confirmar a existência ou não de manipulação indevida.

Outro objetivo pode ser acadêmico: a validação de dados de pesquisas e censos. Por exemplo, dados de população de municípios, ou dados de renda dos indivíduos tendem a ter distribuição conforme a lei de Benford. Assim, desvios dos valores observados em relação aos valores esperados podem ajudar a identificar e corrigir dados anômalos, melhorando a qualidade da estatística.

Vejamos rapidamente alguns exemplos das funções básicas do pacote.

O benford.analysis tem 6 bases de dados reais, retiradas do livro do Mark Nigrini, para ilustrar as análises. Aqui vamos utilizar 189.470 registros de pagamentos de uma empresa no ano de 2010. Os valores vão desde lançamentos negativos (estornos) até valores na ordem de milhões de dólares.

Primeiramente, precisamos carregar o pacote (se você já o tiver instalado) e em seguida carregar os dados de exemplo:

library(benford.analysis) #carrega pacote

data(corporate.payment) #carrega dados

Para analisar os dados contra a lei de benford, basta aplicar a função benford nos valores que, no nosso caso, estão na coluna ‘Amount’.

 bfd.cp <- benford(corporate.payment$Amount)

Com o comando acima criamos um objeto chamado “bfd.cp” contendo os resultados da análise para os dois primeiros dígitos dos lançamentos positivos, que é o padrão. Caso queira, você também pode mudar quantos digitos deseja analisar, ou se quer analisar os dados negativos e positivos juntos, entre outras opções Para mais detalhes, veja a ajuda da função:

 ?benford

Com a análise feita, vejamos os principais gráficos com o comando:

 plot(bfd.cp) 

Os gráficos resultantes se encontram abaixo. Os dados da empresa estão em azul e os valores esperados pela lei de benford em vermelho.

plot_cp

O primeiro gráfico diz respeito à contagem de observações com relação aos seus dois primeiros dígitos, comparando-a com o valor esperado pela Lei de Benford.  Percebe-se que os dados da empresa se ajustam à Lei, mas, também, que há um salto claro no dígito 50!

O segundo gráfico é análogo ao primeiro, mas faz esta contagem para a diferença dos dados ordenados. Como nossos dados são discretos, este saltos em 10, 20, 30,  são naturais e não devem ser encarados como algo suspeito. Por fim, o terceiro gráfico tem um objetivo diferente e, em geral, você também não deve esperar encontrar um bom ajuste dos dados à reta vermelha, principalmente com dados de cauda pesada. Ali se encontra a soma dos valores das observações agrupadas por primeiros dígitos e a  intenção é identificar grupos de valores influentes (que, se estiverem errados, podem afetar bastante uma estatística).

Vejamos agora os principais resultados da análise com o comando

 print(bfd.cp) 

ou somente

bfd.cp

results_cp

Primeiramente são mostrados dados gerais da análise, como o nome da base de dados, o número de observações e a quantidade de primeiros dígitos analisados.

Logo em seguida têm-se as principais estatísticas da mantissa do log das observações. Se um conjunto de dados segue a lei de benford esses valores deveriam ser próximos de:

  • média: 0.5;
  • variância: 1/12 (0.08333…);
  • curtose: 1.2;
  • assimetria: 0.

Que são, de fato, similares aos dados de pagamento da empresa, confirmando a tendência.

Após isso, temos um ranking com os 5 maiores desvios, que é o que mais nos interessa aqui. Veja que o primeiro grupo é o dos números que começam com o dígito 50, como estava claro no gráfico, e o segundo grupo é o dos números que começam com 11. Esses registros são bons candidatos para uma análise mais minuciosa.

Por fim, temos um conjunto de estatísticas de grau de ajuste – que não irei detalhar agora para não prolongar muito este post. Tomemos como exemplo o teste de Pearson, que é bem conhecido. Veja que o p-valor do qui-quadrado é praticamente zero, sinalizando um desvio em relação ao esperado. Mas, como já dissemos várias vezes neste blog, o mais importante não é saber se os dados seguem ou não a lei de benford exatamente. Para isso você não precisaria sequer testá-los. O mais importante é verificar qual o tamanho do desvio e a sua importância prática. Assim, há um pequeno aviso ao final: Real data will never conform perfectly to Benford’s Law. You should not focus on p-values!

Voltando, portanto, à identificação dos desvios, você pode pegar os conjuntos dos dados “suspeitos” para análise com a função getSuspects.

 suspeitos <- getSuspects(bfd.cp, corporate.payment)

Isto irá gerar uma tabela com os dados dos 2 grupos de dígitos com maior discrepância (pela diferença absoluta), conforme ilustrado abaixo. Veja que são exatamente os dados que começam com 50 ou 11. Você pode personalizar qual a métrica de discrepância utilizar e também quantos grupos analisar. Para mais detalhes veja a ajuda da função:

 ?getSuspects

suspects

Note que nossa base de dados é de mais de 189 mil observações. Verificar todos os dados seria infactível. Poderíamos analisar uma amostra aleatória desses dados. Mas, não necessariamente isso seria eficiente.  Veja, assim, que a análise de benford, com apenas os dois primeiros dígitos, nos deu um grupo de dados suspeitos com menos de 10% das observações, permitindo um foco mais restrito e talvez mais efetivo para análise.

Há outras funcionalidades no pacote e na ajuda há exemplos com dados reais. O pacote é bem simples, o intuito é fornecer um mínimo de funções que automatize os procedimentos, facilite a vida e minimize a quantidade de caracteres digitados de quem queira fazer a análise. Algumas funcionalidades que serão adicionadas no futuro são: melhoria na parte gráfica com o Lattice, inclusão de comparação dos dados com a lognormal e inclusão de mais dados de exemplo.

Se encotrar algum bug, tiver alguma dúvida ou se quiser deixar alguma sugestão, comente aqui!

PS: com relação às sugestões, tanto o R quanto este pacote são livres e com código aberto. Então, sinta-se à vontade não somente para sugerir, mas principalmente para escrever novas funções e funcionalidades!

Investimento Estrangeiro Direto no Brasil por Estado (Indústria)


Os dados do Censo de Capitais Estrangeiros no País, em 2010, trouxeram a distribuição do Investimento Estrangeiro Direto (IED) na indústria por Unidade da Federação (UF).

Somente da Indústria? E como foi feita a distribuição? Aqui voltamos ao que já dissemos sobre erro de medida (ver aqui, aqui, aqui e aqui, por exemplo). Distribuir o estoque investimento estrangeiro por UF é algo complicado, sujeito a erros diversos, tanto ao se definir a metodologia, quanto ao se mensurar o valor. No censo de 1995, por exemplo, os dados foram distribuídos por estado “[…] tomando por base o endereço da sede da empresa”. Será que essa é uma boa medida? Depende.

Percebe-se que uma indústria que concentre o grosso da sua estrutura produtiva no Pará, mas que tenha sede em São Paulo, será considerada um investimento nesta última UF.  Se a intenção é medir onde se encontra o centro administrativo, esta medida poderá ser boa. Todavia, se intenção é medir onde se encontram as unidades produtivas, esta medida terá, talvez, distorções significativas. Qual a melhor forma, então, de se distribuírem os investimentos por estado? Pela localização da sede? Pela localização do ativo imobilizado? Pela distribuição dos funcionários?  Particularmente, acho que não existe uma métrica única que se sobressaia às demais – a melhor opção depende do uso que você irá fazer da estatística.

Voltando ao Censo, a pesquisa passou a considerar a distribuição do ativo imobilizado como critério para alocação do IED – e apenas para a indústria . Os declarantes distribuíram percentualmente o seu imobilizado pelos diferentes estados e isso foi utilizado para ponderar o investimento direto pelas UF’s.

Segue abaixo mapa do Brasil com a distribuição do IED da indústria por Unidade Federativa:

IED_UF

 Para aprender a fazer o mapa, veja aqui.

Analisando microdados do IBGE com o R


Os materiais do Seminário de Metodologia do IBGE de 2013 estão disponíveis para download. Dentre eles, destaco o do mini-curso Introdução à análise de dados amostrais complexos. Lá você vai aprender a replicar os resultados da POF, da PNAD e amostra do Censo levando em conta o desenho amostral das pesquisas (que é necessário para se calcular corretamente medidas de precisão, como a variância). O material é bastante focado no blog de Anthony Damico, Analyze Survey Data for Free. O blog é fantástico, com diversos exemplos de como baixar e analisar dados de pesquisas públicas levando em conta o plano amostral, tudo com ferramentas gratuitas como o R.