Cachorros artificiais: agentes de agent-based modeling usando R (ABM – Parte 1)

(Primeiro Post da série. Veja os outros: Parte 2, Parte 3, Parte 4 e Parte 5)


robot-dog-meets-a-real-doberman-dog

Modelos baseados em agentes (ou Agent-Based Models – ABM) já foram um tópico discutido no blog outras vezes (por exemplo, aqui e aqui). Mas julgo que esse é um assunto “quente”, que sempre vai acabar voltando.

Para aqueles que não sabem ao que estou me referindo, vou dar uma ideia simples e rápida. Trata-se de simulações de algo como “sociedades artificiais” através de modelos computacionais. O propósito NÃO É construir a Matrix ou qualquer coisa do tipo… ao, contrário, é algo bem mais “humilde” e focado. Funciona assim: o pesquisador tem hipóteses sobre como indivíduos ou grupos se comportam, cria um ambiente simples que contém apenas algumas características fundamentais e indispensáveis para a caracterização de uma situação típica e então toca o play. Com isso, pode investigar se os mecanismos interacionais supostos de fato produzem, no nível macro ou agregado, os padrões efetivamente observáveis. Trata-se de um exercício lógico, antes de qualquer outra coisa.

Economistas neoclássicos, trabalhos que envolvem modelagem formal e escolha racional fazem exercícios lógicos de natureza semelhante o tempo todo: “suponha indivíduos homogêneos (idênticos em todas as suas características), com determinado tipo de preferência e orçamento; de fronte a um mercado que se lhes apresente certas quantidades e preços, eles se comportarão desta e daquela maneira, gerando um estado de equilíbrio tal e qual…”. No entanto, nem todo tipo de modelagem deve envolver escolha racional… e nem matemática. Certos problemas são complexos o suficiente para não permitir “forma fechada” (i.e. uma equação bem definida, cuja solução informa um resultado predito). Os ABMs entram em cena exatamente aí.


Um agente é uma entidade conceitual, que pode ser representada de diversas formas e nas mais diversas linguagens de programação. Mas duas propriedades fundamentais devem estar presentes: atributos e ações. Um indivíduo, por exemplo, tem diversas características: idade, peso, altura, raça, sexo, ocupação etc… Esses são os atributos. E também é capaz de agir no mundo de diversas formas: correr, gritar, socializar, procurar emprego etc. Um “agente”, nesse sentido a que me refiro, seria a representação computacional dos aspectos analiticamente relevantes de um indivíduo para determinados fins de pesquisa. Ou seja, ele não contém todos os atributos e ações possíveis de um indivíduo, mas apenas aqueles imprescindíveis para a análise. Isso vai ficar mais claro adiante.

O exemplo que vou dar será utilizando o R, minha “linguagem nativa” (que uso no dia a dia para fazer análises estatísticas). Mas reconheço que ele não é o melhor ambiente pra isso… Existem linguagens de programação específicas para ABM (a mais conhecida é a LOGO, que se tornou bastante popular devido ao software NetLogo, que implementa uma de suas versões). Também bastante utilizada é a linguagem Python — que é uma linguagem de programação completa (não específica para ABM) e tão flexível quando o R (e muito parecida, inclusive). ABMs escritos em Python têm scripts mais simples, diretos e rodam mais rápido; de acordo com minha experiência (mas esse último ponto é controverso).

As estruturas de dados mais simples e comuns no R são vetores, matrizes, listas e data.frames. Basicamente, todos os outros tipos de objetos existentes são algum tipo de combinação ou transformação desses… Além disso, os comandos que executamos também são objetos e são chamados funções. Criar um agente é guardar seus atributos em algum desses tipos de estrutura de dados e representar suas ações por funções. Uma forma conveniente de mesclar atributos e funções num único objeto-agente é fazer uso de um tipo de estrutura de dados pouco utilizada no R, chamada Reference Class (veja mais coisas sobre isso aqui).

Dou um exemplo. Desejamos representar um cachorro como agente. Os atributos relevantes serão a idade e a raça. E as ações possíveis serão latir e nos dizer qual é a sua raça e a sua idade, se lhe perguntarmos (sim! é um cachorro falante! e daí!?). Criaremos uma Reference Class do tipo “Cachorro”

Cachorro = setRefClass("Cachorro",
             fields = c("raca",
                        "idade"),

             methods = list(
                 late = function(){
                     print("Au! Au!")
                 },

                 diz_raca = function(){
                     print(paste("Eu sou um", .self$raca))
                 },

                 diz_idade = function(){
                    print(paste("Tenho", .self$idade,
                                "anos, o que significaria",
                                 7*.self$idade,
                                 "para um humano" ))
                 }
) )

Acredito que para quem conhece um pouco de R, mesmo sem jamais ter visto uma Reference Class antes, o procedimento acima é mais ou menos inteligível. Com o comando setRefClass estamos criando esse novo tipo de objeto — tipo que será chamado “Cachorro”. Ele terá dois atributos, cujos nomes estão declarados no argumento  fields (como vetores character). Depois especificamos as ações que os Cachorros são capazes de realizar: latir, dizer a raça e dizer a idade. Observem uma coisa importante, que é específica das Reference Classes: a expressão .self. Quando .self$raca estamos dizendo que o atributo raça está contido dentro do próprio objeto Cachorro, não é preciso procurá-lo no ambiente principal do R. Essa é a principal característica das Reference Classes: a capacidade de se auto-referenciar. Qualquer pessoa familiarizada com C++,  Java ou Python logo identifica que essa é uma das principais características de qualquer “linguagem orientada a objetos“. Essa expressão simplesmente significa que os objetos criados são auto-contidos: trazem consigo informações, dados e também os métodos para acessá-los, modificá-los e utilizá-los.

Tá… e agora? Agora que definimos as propriedades de um cachorro em abstrato, o próximo passo é efetivamente criar um (é como se tivéssemos definido a “idéia” de cachorro e agora fossemos de fato criar um “cachorro empírico”).

Jimmy = Cachorro$new(raca = "Vira-latas", idade = 4)

Acredito que o código acima seja mais o menos intuitivo. Mas é bom discuti-lo. Observem que estamos utilizando o $ para acessar uma função que está “dentro” do objeto Cachorro, chamada new(). Essa não era uma das ações definidas no argumento methods (late, diz_idade, diz_raca). Trata-se de uma função presente em todo e qualquer objeto do tipo Reference Class, que significa, em termos simples, “crie pra mim uma espécie ou instancia empírica dessa classe”. Informamos as características daquele exemplar particular de cachorro (Vira-latas, com 4 anos). E pronto, fiat Jimmy. Esse é o Jimmy “por dentro”:

jimmy1

Ele é capaz de latir, dizer a idade e a sua raça:

jimmy2

Vamos agora criar novamente a classe Cachorro, mas com duas novas capacidades de ação: dar a pata e abaixá-la. Além dessas duas novas funções, que deverão ser informada dentro da lista passada para o argumento methods, temos que também definir um novo atributo, o “status” da pata (i.e., se está levantada ou abaixada). Se já estiver levantada e pedirmos mesmo assim para que ele a levante, o cachorro vai nos informar. O mesmo se já estiver abaixada e ainda assim pedirmos a ele para abaixá-la.

Cachorro = setRefClass("Cachorro",
             fields = c("raca",
                        "idade",
                        "pata"),

             methods = list(
                 initialize = function(pata = 0, ...){
                     .self$pata = pata
                     callSuper(...)
                 },
                 late = function(){
                     print("Au! Au!")
                 },

                 diz_raca = function(){
                     print(paste("Eu sou um", .self$raca))
                 },

                 diz_idade = function(){
                    print(paste("Tenho", .self$idade,
                                "anos, o que significaria",
                                 7*.self$idade,
                                 "para um humano" ))
                 },

                 da_pata = function(){
                     if( .self$pata == 0 ){
                         .self$pata = 1
                         print("Levantei a pata")
                     }else{
                         print("Já te dei a pata, oras...!")
                     }
                 },

                 abaixa_pata = function(){
                     if( .self$pata == 0 ){
                         print("Minha pata já está abaixada!")
                     }else{
                         .self$pata = 0
                         print("Abaixei a pata")
                 }
}
) )

Observem que uma outra função teve que ser definida dentro de methods, denominada initialize. Essa função simplesmente especifica valores-padrão que serão passados para todos os agentes daquela classe no momento de sua criação. Neste caso, todos os cachorros começam com pata = 0 (com a pata abaixada). Vejamos agora o que o Tobby, nosso novo cachorro, faz:

jimmy3

Tobby não é bobo. Não peça pra que ele dê a pata ou a abaixe duas vezes seguidas. E Tobby tem a capacidade de “modificar-se a si mesmo”. Quando levanta ou abaixa a pata, ele altera o atributo pata, que existe dentro de si.

O ponto aqui é que Reference Classes são adequadas para fazer uma aproximação operacional dos conceitos de atributos e ações. Poderíamos fazer isso de inúmeros outros modos, sem apelar para programação orientada a objeto. Mas há simplicidade nos códigos acima. Agentes agem (ou, nos termos de Taylor Swift, “players are gonna play”). Essa representação sugere uma analogia com a realidade. Como diria meu amigo Davoud, especialista em ABM, “it’s all about ontology”.

 

 

 

Anúncios

Decompondo as desigualdades: material para a replicação completa de “Os impactos da geração de empregos…”

a06graf01

Publiquei na RBCS: “Os impactos da geração de empregos sobre as desigualdades de renda: uma análise da década de 2000“. O link na página do Scielo é esse aqui. Pra quem prefere o PDF (que ficou bem bonitinho), é esse aqui.

O texto é em co-autoria com Flavio Carvalhaes, Pedro Herculado F. de Souza e Carlos Costa Ribeiro. Estou muito feliz. Mas produzir todas as análises foi um processo longo e de muito aprendizado…

Nosso trabalho trata do seguinte: houve uma grande geração de empregos na década de 2000 e, simultaneamente, grande queda da desigualdade de rendimentos… Perguntamos então: como esses dois fenômenos se relacionam? Noutros termos: a mudança composicional da força de trabalho (da distribuição dos indivíduos entre as ocupações) exerceu influência sobre o verificado movimento de queda?

Bem… descobrimos que sim, mas esse não foi o fator principal. A geração de empregos foi um fenômeno positivo, que trouxe melhoria dos postos de trabalho existentes (o que pode trazer efeitos de mais longo prazo); mas sua contribuição imediata para o saldo de queda das desigualdades na última década foi de 18%. Outros fatores — principalmente relacionados à educação — foram mais importantes (o que corrobora outras pesquisas já realizadas sobre o assunto). No entanto, justamente porque o “componente ocupacional” não caiu tão depressa quanto os demais, hoje em dia, sua participação na parcela restante de desigualdade de renda (ainda muito alta!) se tornou mais importante. Noutros termos, o movimento de queda trouxe mudança qualitativa da composição das desigualdades.

Para aqueles que se interessam, esse aqui é o link para o material completo de replicação do nosso texto. Tudo foi feito no R. Na pasta principal, há dois scripts que executam toda a análise (não é preciso acessar os demais, localizados dentro das outras pastas). É preciso apenas que o usuário mude, dentro dos scripts, o nome das pastas onde os arquivos estão localizados.

As bases de dados utilizadas podem ser baixadas no site do Centro de Estudos da Metrópole.

Investimento Estrangeiro Direto no Brasil (mapa por País de Origem Imediata)

Um uso muito legal de coisas sobre as quais já conversamos aqui no Sociais & Métodos.
PS.: Conheçam o Análise Real, Blog de Carlos Cinelli

Análise Real

Que tal visualizar os dados do Censo de Capitais Estrangeiros de uma maneira diferente?

Abaixo, mapa com a distribuição do Investimento Estrangeiro Direto (IED) no Brasil, critério participação no capital, em 2010, segundo o país de origem imediata. O mapa foi feito no R. Quanto mais escuro, maior o investimento daquele país em empresas brasileiras.
IED_Pais

PS: agradeço ao Rogério pelo didático post ensinando o caminho das pedras.

Ver o post original

Milhões de casos em segundos: os Censos no R

funilCensos Demográficos são bancos de dados muito pesados! Milhões de casos… Operações simples, como frequências, médias e proporções podem demorar muitos minutos. Modelos estatísticos complexos podem demorar horas… ou dias.  O  R convencional (assim como o Stata) carrega todos as informações com as quais está trabalhando na memória RAM. Ou seja, as análises são realizadas muito rapidamente, mas não é possível abrir um banco que seja maior do que a memória disponível. O SPSS e o SAS executam as análises a partir dos arquivos no HD – assim, suportam bancos grandes, mas são muito lentos.

Neste em 2012 e 2013 trabalhei muito com os Censos brasileiros (cf. Projeto Censo). Era necessário descobrir uma maneira de agilizar as análises. Foi quando descobri uma versão do feita exatamente para lidar com grandes bancos de dados, produzida pela Revolution AnalyticsEssa  empresa construiu um software (não aberto) em cima do R convencional — o RevolutionR. Ele funciona com a mesma linguagem, suporta os mesmos pacotes e funcionalidades — mas com alguns adicionais (como dizem, “100% R and more“). Uma versão gratuita do RevolutionR para uso individual ou acadêmico pode ser baixada aqui. Uma das principais diferenças é a presença de um pacote para “big data”, que não pode ser instalado no R convencional, chamado RevoScaleR. Quais são suas vantagens?

  • É possível trabalhar com dados de qualquer tamanho, pois ele acessa informações a partir do HD (assim como o SPSS e o SAS), superando os limites de memória RAM.
  • O RevolutionR permite salvar dados num formato próprio, com extensão XDF. Nesses arquivos, grandes bancos de dados são fragmentados e salvos em blocos separados (como uma planilha com múltiplas abas), chamados chunks ou blocks. Processar cada bloco de uma vez é menos pesado do que trabalhar com o banco todo (esse é o princípio já adotado, no R, por pacotes como ff, ffbase e biglm).
  • As funções de análise já acessam e integram os resultados dos múltiplos blocks. E tudo isso já fazendo uso de processamento paralelo — usando todos os cores de um computador ou vários computadores ao mesmo tempo.

Resumo: uma regressão múltipla com 25 milhões de casos pode ser feita em menos de 10 segundos.

Vamos por a mão na massa então. Vou usar os dados do Censo de 2010, que podem ser baixados no site do Centro de Estudos da Metrópole, como já dissemos aqui (lembrem-se dessas instruções, antes de usar os microdados).

Em primeiro lugar, é necessário dizer que vamos usar computação paralela, todos os cores do computador:

# Define um contexto de computação paralela
parallelContext = RxLocalParallel()
rxSetComputeContext(parallelContext)

Depois temos que salvar o arquivo .SAV (formato SPSS) como XDF (formato RevolutionR). Essa é a única parte demorada: 40 minutos ou mais… Mas só precisa ser feita uma vez.

#Indicando o local do arquivo
censo2010spss = RxSpssData(file.path("C:/dados/censo2010.sav"))

#"Salvando como"
rxImportToXdf(
	inSource = censo2010spss,
	outSource = "C:/dados/censo2010.xdf",
	reportProgress = 1,overwrite = TRUE)

Então apontamos o local do arquivo XDF no HD:

censo2010 = file.path("C:/dados/censo2010.xdf")

Agora é moleza. Que tal um gráfico da média de renda por idade?

# Tabela: Média de renda por idade -- 10 a 80 anos  (com peso)
censoCube = rxCube(renda.trab ~ F(idade.anos), data=censo2010, 
   rowSelection = (idade.anos >= 10) & (idade.anos <= 80), 
   fweights = "peso_pessoa") 

# Tranformando a idade de factor para numeric
censoCube$idade.anos = as.integer(levels(censoCube$F_idade.anos))

# Transformando a tabela num data.frame
censoCube = as.data.frame(censoCube) 

# Pronto: gráfico de renda por idade
rxLinePlot(renda.trab ~ idade.anos, data=censoCube,
   title="Relação entre Renda e Idade",
   xlim=c(0,85), ylim=c(0,2000))

O que acham?

grafico

Bonitinho, né!? E o mais interessante: Computation time: 3.697 seconds. Mais? Vamos repetir o mesmo gráfico, mas com recorte por sexo.

# Uma regressão para estimar as médias de renda
# por idade e sexo (efeito interativo). A opção
# "cube=TRUE" indica que salvaremos o cruzamento
linMod = rxLinMod(renda.trab~sexo:F(idade.anos), data=censo2010,
   pweights="peso_pessoa",  cube=TRUE, 
   rowSelection = (idade.anos >= 16) & (idade.anos <= 80)) 

# Extraio a tabela de cruzamento
plotData = linMod$countDF

# Converto a idade de factor para numeric
plotData$idade.anos = as.integer(levels(plotData$F.idade.anos.))[plotData$F.idade.anos.]

# Gráfico de idade e renda por sexo
rxLinePlot(renda.trab~idade.anos, groups=sexo, data=plotData,
   xlim=c(0,85), ylim=c(0,2500))

Computation time: 4.053 seconds Nada mal, né!? Uma regressão, transformações em variáveis e um gráfico…

grafico2

Pois é… seu laptop era um super computador e você nem sabia…

Nota:  As flutuações uma flutuações na renda do trabalho nas idades mais avançadas provavelmente são devidas ao baixo número de casos de trabalhadores nessas faixas. Logo, outliers passam a influenciar mais as observações. Além disso, podem ter ocorrido problemas na captação da informação sobre renda e há limites intrínsecos à própria amostra (os microdados do Censo abarcam +/- 10% da população, cerca de 21 milhões de casos). Mas esta não é uma particularidade do Censo de 2010. Em outros censos pesquisas (como PNADs, PMEs, PEDs etc) encontramos o mesmo.

Como medir desigualdades: Slide e materiais

Compartilho aqui os slides e os materiais utilizados em duas seções das quais participei na disciplina “Desigualdades: conceito, mensuração e novas abordagens”, ofertada por Márcia Lima no curso de Ciências Sociais da USP.

Achei que o diálogo foi bem legal. Os materiais contém bancos de dados utilizados como exemplos, scripts em R que ilustram e aplicam os conceitos apresentados, bem como dois textos pequenos e introdutórios. O link para baixar tudo é esse aqui. Mas não sei se essas coisas bastam por si mesmas, sem a parte expositiva.

Um adendo sobre os scripts em R:

1 – Há a replicação da regressão de prestigio sócio-ocupacional de Duncan (1961), que serviu, posteriormente de base para a construção de outras escalas contínuas. Inclui um plot em 3d, hehe.

2 – Aplicação de diversas medidas de desigualdade de renda a duas edições das PNADs (acompanha banco de microdados — com apenas algumas variáveis).

3 – Simulações para ilustrar a sensibilidade de cada indicador de desigualdade de renda.

Aluguéis em SP: o viés das imobiliárias on-line

Procurar um lugar pra morar é bastante estressante — principalmente tendo em vista que tempo e dinheiro são bens escassos. Contudo, aparentemente, essa tarefa se tornou mais fácil, com o advento dos sites que conjugam num só lugar diversas imobiliárias. Pretensamente, esses sites também acirrariam a competição, por diminuírem o custo de obtenção da informação, “empoderando” o potencial cliente. Ou seja, com mais competição, se espera menor preço. Mas recentemente descobri que isso pode não ser bem verdade…

Li o post São Paulo, uma cidade gentrificada, do Arquitetura da Gentrificação e fiquei pensando se eu conseguiria fazer um levantamento de preços de aluguéis a partir desses sites de imóveis e então calcular a média do que se costuma cobrar por metro quadrado. Construí dois scripts de webscrapping em R: um para fazer download das informações do ImóvelWeb e outro para o ZapImóveis. Mas neste post, vou explorar apenas as informações do ImóvelWeb.

Baixei informações sobre apartamentos para locação (excluindo casas, lofts, studios, kitnets e imóveis comerciais). Montei um banco de dados com informações sobre preço, área (metros quadrados), bairro, nº de banheiros, quartos, suítes e vagas em garagem. Sei que esses sites tem entradas repetidas de imóveis (às vezes devido ao fato de estarem simultaneamente no cadastro de duas ou mais imobiliárias). Exclui os duplicados — certamente jogando alguns fora indevidamente e ainda deixando passar outros (que apesar de aparecerem mais de uma vez, não tem registros exatamente idênticos).

Depois obtive a latitude e longitude dos pontos médios de cada bairro onde estavam localizados com uso da função geocode (sobre a qual já tratei em post anterior). Sim… seria mais legal se fosse possível geocodificar os endereços exatos dos apartamentos, mas não há essa informação para todos os casos. No final das contas fiquei com dados para 1426 apartamentos.

Baixei um mapa do GoogleMaps no R e plotei a mediana do valor por metro quadrado para cada bairro. O resultado é a figura abaixo:

precos

Algumas observações

  1. De forma nada surpreendente (para alguém que vive em São Paulo), os lugares mais caros são: Itaim, Moema, Brooklin, Jardins, Morumbi… Mas esses dados são um pouco diferentes daqueles do Arquitetura da Gentrificação: Pinheiros, Perdizes e Pompéia são sim caros… mas não estão no mesmo páreo que aqueles.
  2. Reparem nos preços: grande parte das bolinhas azuis são indicadoras de valores entre 25 e 45 R$/m². Isso parece bem acima dos valores médios identificados por eles.

Obviamente há diferenças importantes na composição da amostra e nas fontes de dados. Não são dados completamente comparáveis. Creio que as duas pesquisas, a minha e a deles, dizem respeito a coisas diferentes. Assumo que os dados deles informam “de fato” sobre a tendência dos aluguéis em SP… Os meus informam sobre o características do ImóvelWeb — que podem ser avaliadas tomando os dados deles como parâmetros.

Plotei também uma mancha de densidade dos imóveis para locação, para ter uma noção da cobertura e dos lugares mais anunciados. No mapa abaixo, quanto mais azul claro, maior a frequência de imóveis anunciados na região:

cobertura

As conclusões (tentativas) são simples:

  1. Os valores de aluguel informados pelas imobiliárias online estão acima da média da cidade para quase todos os bairros (as bolinhas, afinal, indicam preços comparativamente mais altos por toda extensão do mapa). Alugar pela internet parece ser mais caro!
  2. O ImóvelWeb (e possivelmente outros sites desse tipo — ainda tenho que verificar) possui um viés de cobertura em favor dos bairros ricos. O número de imóveis notificados em determinadas localidades é mínimo. Não se reduz o custo informacional igualmente para todos os lugares.

Arrisco a razão para esses fatos: provavelmente apenas as maiores imobiliárias anunciam nesses sites. Assim, despontam na frente das pequenas, em termos de visibilidade, e podem cobrar um preço maior. A competição talvez ocorra sim, mas apenas entre essas grandes. Suspeito que talvez encontraríamos preços mais baratos batendo na porta de alguns escritórios, à moda antiga mesmo. Se a razão que apontei está correta, fica também fácil compreender o viés de cobertura: as maiores imobiliárias são aquelas que atendem os ricos.

O script para a replicação completa está disponível neste link. Não está muito bonitinho, nem muito formatado. É a pressa.

Mapas no R, parte 2: utilizando shapes

brasilTrabalhar com mapas ou dados georreferenciados frequentemente significa fazer uso de shapes. Em linhas gerais, shapes são de arquivos de mapas formados por polígonos geocodificados (ie., suas coordenadas são latitudes e longitudes) que delineiam os contornos de unidades espaciais de interesse, como limites municipais, estaduais etc. Basicamente o que se pretende nesses casos é “colorir” os polígonos de acordo com alguma característica ou atributo relevante. Este é o tipo de mapa que Leonardo utilizou em seus posts sobre médicos no Brasil e sobre transportes na cidade de SP.

Trabalhar com mapas como polígonos significa manusear (ao menos) três arquivos conjuntamente: um que possui a extensão .shp (que traz o desenho do mapa propriamente), outro com extensão .shx (que é um indexador de informações, para facilitar buscas) e um terceiro com extensão .dbf (que é um arquivo de banco de dados, trazendo informações e atributos das unidades espaciais – p.ex.: população do município, renda per capita etc.) O que liga esses três arquivos é um ID, único para cada unidade espacial, geralmente chamado de Geocódigo. Há uma explicação mínima sobre shapefiles na Wikipédia em português e uma explicação mais completa na versão em inglês.

O objetivo desse post é replicar este mapa simples sobre PIBs municipais, disponível no site da Wikipédia. O exemplo é de Diogo Ferrari.

Precisamos então seguir os seguintes passos:

  1. Fazer download de um shape de estados brasileiros
  2. Fazer um webscrapping para obter as informações dispostas na tabela do Wiki (e “limpar” tais informações)
  3. Definir aspectos da formatação do mapa (cores, rótulos da legenda)
  4. Parear os dados da tabela e as informações do arquivo de atributos (.dbf)
  5. Plotar

Primeiramente, carregamos os pacotes necessários:

require(XML)
require(RCurl)
require(maptools)
require(RColorBrewer)

Em seguida, para guardar os arquivos do mapa que serão baixados, criamos uma pasta e a definimos como diretório de trabalho (esse passo é opcional):

dir.create("c:/mapas/")
setwd("c:/mapas/")

Agora utilizaremos a função download.file para baixar um arquivo compactado (.zip) que contém o mapa (.shp, .shx e .dbf).  Isto está disponível no site Gismaps. Em seguida aplicamos a função unzip para extrair o conteúdo compactado. Por fim, para ler o mapa e guardá-lo num objeto do R, aplicamos a função readShapePoly do pacote maptools. Lendo o arquivo .shp já estamos implicitamente fazendo referência ao .shx e ao .dbf.

# Carregando dados do Mapa estadual
download.file("https://dl.dropboxusercontent.com/u/35964572/estados_2010.zip",destfile="estados_2010.zip")

unzip("estados_2010.zip")

mapaUF = readShapePoly("estados_2010.shp")

Nesse ponto, o mapa já é “plotável”:

plot(mapaUF)

plot_zoom

Mas só conseguimos ver seus formatos… não há cores que indiquem quaisquer informações sobre cada polígono (unidade da federação). O arquivo de atributos (.dbf) não contém ainda dados interessantes. Por isso, vamos ao passo 2: fazer o download das informações sobre PIBs estaduais.

Vamos usar a função readHTMLTable, já apresentada num post anterior.

PIBEstadual = readHTMLTable(getURL("https://pt.wikipedia.org/wiki/Anexo:Lista_de_estados_do_Brasil_por_PIB", ssl.verifypeer = FALSE),which=2)

Como praticamente toda informação baixada na internet, os dados da tabela chegam “sujos”  ou desformatados no R. É preciso fazer alguns procedimentos antes de utilizá-la:

# Mantemos apenas as colunas 3 e 4, que contém o nome da UF e seu PIB
PIBEstadual = PIBEstadual[2:nrow(PIBEstadual),3:4]

#renomeando as colunas
names(PIBEstadual)=c("UF","PIB")

#retira o símbolo '.' , que estava dividindo os milhares
PIBEstadual$PIB = gsub("\\.", "",PIBEstadual$PIB)

#substitui o nome pela sigla das UFs
PIBEstadual$UF = c('SP', 'RJ', 'MG', 'RS', 'PR', 'BA', 'SC', 'DF', 'GO', 'PE', 'ES', 'CE', 'PA', 'AM', 'MT', 'MA', 'MS', 'RN', 'PB', 'AL', 'SE', 'RO', 'PI', 'TO', 'AC', 'AP', 'RR')

# informa que os dados contidos na coluna 'PIB' são números
PIBEstadual$PIB = as.numeric(PIBEstadual$PIB)

# Divide o valor do PIB por 10.000, para que seja expresso em bilhões de reais
# (já está em R$ 1.000)
PIBEstadual$PIB = PIBEstadual$PIB/10^6

# Transforma os dados do PIB em uma variável categórica.
PIBEstadual$PIB_cat = cut(PIBEstadual$PIB, breaks=c(0,30,60,90,120,500,5000),
 labels=c('até 30','+ 30', '+ 60', '+ 90', '+ 120', '+ 500'))

Passamos ao passo 3: definir a formatação.

# Selecionamos algumas cores de uma paleta de cores do pacote RColorBrewer
paletaDeCores = brewer.pal(9, 'OrRd')
paletaDeCores = paletaDeCores[-c(3,6,8)]

# Agora fazemos um pareamento entre as faixas da variável sobre PIB (categórica) e as cores:
coresDasCategorias = data.frame(PIB_cat=levels(PIBEstadual$PIB_cat), Cores=paletaDeCores)
PIBEstadual = merge(PIBEstadual, coresDasCategorias)

Agora o passo 4: fazer o pareamento entre os dados da tabela e o mapa:

# Primeiramente, guardamos os dados do .dbf num objeto do R.
# Ele é um atributo do objeto mapaUF
mapaData = attr(mapaUF, 'data')

# Guardamos o número das linhas numa nova variável
# Esse passo é necessário pois o pareamento entre esses dados e a tabela do PIB
# muda a ordem os casos, o que prejudica, depois, a construção do mapa
mapaData$Index = row.names(mapaData)

# Mudando o nome da variável que indica a sigla dos estados
names(mapaData)[3] = "UF"

# Fundimos então as duas tabelas:
mapaData = merge(mapaData, PIBEstadual, by="UF")

# Reordenamos os dados do mapa
mapaData = mapaData[order(as.numeric(mapaData$Index)),]

# E guardamos essas informações novamente como sendo um atributo do arquivo de mapa.
attr(mapaUF, 'data') = mapaData

Agora, o passo final, plotar:


# Configurando tela (reduzindo as margens da figura)
parDefault = par(no.readonly = T)
layout(matrix(c(1,2),nrow=2),widths= c(1,1), heights=c(4,1))
par (mar=c(0,0,0,0))

# Plotando mapa
plot(mapaUF, col=as.character(mapaData$Cores))
plot(1,1,pch=NA, axes=F)
legend(x='center', legend=rev(levels(mapaData$PIB_cat)),
 box.lty=0, fill=rev(paletaDeCores),cex=.8, ncol=2,
 title='Mapa dos estados brasileiros segundo o PIB em\n2010. Em bilhões de reais:')

Pronto:

mapa.pib.estados

Aprendendo R – Tutorial interativo + Vídeos

codeschoolO R aparece muito por aqui… Um dos objetivos deste blog é popularizá-lo.  Fato: o custo inicial de aprendizagem é um pouco elevado, se comparado aos de outros softwares. Mas vencido o atrito inercial, as possibilidades analíticas são muito maiores.

Este breve post tem o objetivo de apresentar dois tutoriais. O primeiro deles é o tutorial interativo da Code School. Esse tutorial é composto de uma série de explicações e exercícios de nível básico, que podem ser concluídos muito rapidamente. O próprio site simula o prompt de comando do R — nem é preciso instalar nada. Mas tem duas limitações: 1) não chega em conteúdos de nível intermediário e 2) está em inglês (mas não há nada muito avançado). Essa segunda limitação fere um pouco do propósito deste blog, que é popularizar e apresentar materiais em português…

A segunda indicação é uma série de vídeos que está sendo produzida por um usuário do Youtube, Henrique Gomide. Parece legal também. A partir do vídeo abaixo, é possível acessar o restante da lista de reprodução. Até o momento de publicação deste post, já eram 18 pequenos vídeos.

Mapas e dados georreferenciados no R – parte 1

Existem diversos softwares especializados para trabalhar com dados geográficos, que possuem funções das mais diversas, dedicadas desde desenhar e/ou corrigir polígonos, passando por gerenciamento de bancos dados até executar complexos modelos de estatística espacial. No entanto, pouco a pouco, a comunidade de usuários do R vem implementando pacotes, funções e interfaces para lidar com sistemas de informação geográficas (GIS). Talvez ainda hoje não possibilitam todos os usos que um software como o ArcGIS. Mas conhecendo um pouco da grande velocidade e gana com que se desenvolvem pacotes para o R, é bem provável daqui a pouco isso seja revertido (se é que é verdade!).

Neste post ilustro uma forma de lidar com dados georreferenciados (latitude/longitude) e plotar mapas. O pacote a ser usado é o ggmap, que possui uma função chamada geocode. Seu uso é bem simples: basta informar um endereço postal e recebemos de volta sua latitude e longitude. Detalhe: não pode haver acentos nas vogais, nem “ç”. Os caracteres aceitos são aqueles do teclado americano. :-/

Vamos obter a latitude e a longitude do MASP, Museu de Arte de São Paulo:

require(ggmap)
geocode("Av. Paulista, 1578 - Bela Vista, Sao Paulo, Brasil")

Information from URL : http://maps.googleapis.com/maps/api/geocode/json?address=Av.+Paulista,+1578+-+Bela+Vista,+Sao+Paulo,+Brasil&sensor=false
Google Maps API Terms of Service : http://developers.google.com/maps/terms
lon lat
1 -46.65726 -23.56049

A função geocode acessa o API do Google Maps e faz uma pesquisa pelo endereço informado. Obviamente, nem todos os endereços são unívocos. Qualquer um que já utilizou o Google Maps sabe que esse aplicativo retorna um rol de alternativas em caso de múltiplas possibilidades. O resultado acima é a alternativa mais provável (e, na realidade, é a correta). Mas se for o caso, é possível acessar, dentro do objeto que é resultado dessa pesquisa, as demais alternativas. É possível também procurar vários endereços de uma vez só. Basta que cada endereço esteja gravado dentro de um elemento de um vetor de caracteres (string ou character). O uso do API, no entanto, tem uma limitação: só é possível fazer 2500 consultas por dia. Limites são muito comuns em APIs de acesso gratuito. Há planos pagos que permitem elevar aquele número.

O que fazer com a latitude e a longitude? É possível plotá-las num mapa!

No exemplo abaixo, vou até o site da prefeitura de SP, baixo o código fonte da página que contém uma lista de endereços de 13 unidades do SESC na região metropolitana (existem mais unidades… mas para o exemplo, essas bastam).

# Baixando, na página da prefeitura de SP, os endereços dos SESCs
pag_pref = readLines("http://www9.prefeitura.sp.gov.br/secretarias/smpp/sites/saopaulomaisjovem/cultura/index.php?p=25")[71]

Reparem que já extraio diretamente a linha de número 71. É exatamente nela que estão todos os endereços. Então baixar o resto era desnecessário. Mas esse conteúdo não está formatado:

pag

É preciso fazer uma limpeza do código HTML e extrair apenas o conteúdo desejado:

# Separando as linhas e removendo conteúdos desnecessários
pag_pref = unlist(strsplit(pag_pref,"<br />|/ Tel"))
pag_pref =gsub("<strong>|</strong>","",pag_pref)

# Mantém apenas as linhas que contêm a expressão "Endereço"
pag_pref = pag_pref[grep("Endereço",pag_pref)]

# Remove a expressão "Endereço"
pag_pref = gsub("Endereço: ","",pag_pref)

# Retira todos os caracteres especiais
pag_pref = gsub("[[:punct:]]", ",", pag_pref)

# Remove conteúdo desnecessário da linha 1
pag_pref = gsub("esquina com a Rua 24 de Maio","Sao Paulo, SP",pag_pref)

# Adiciona a cidade à linha 8
pag_pref[8] = paste(pag_pref[8],", Sao Paulo, SP")

# Adiciona o país a todas as linhas
pag_pref = paste(pag_pref,", Brasil")

# Remove todos os acentos
pag_pref = iconv(pag_pref, to="ASCII//TRANSLIT")

Pronto, esse é o resultado, uma lista de endereços guardada num vetor de caracteres:

pag

Agora é só aplicar a função geocode para obter todas as latitudes e longitudes:

latlon = geocode(pag_pref)

Falta o mapa.

A função get_map do pacote ggmap acessa um repositório público de mapas, o stamen, copia a imagem (geocodificada) da região desejada e a retorna como um objeto do R. Há várias opções de formatação, cores etc. E há também outros repositórios de mapas (inclusive o próprio Google Maps).

# Baixa o mapa de SP (centrado na Sé - isto pode ser alterado)
sp.map =get_map(location="Sao Paulo", zoom = 11,
 source = "stamen", maptype = "toner", color = "color")

# Transforma o arquivo de mapa em um grafico (ggplot2)
sp.map.2012 <- ggmap(sp.map, base_layer =
 ggplot(aes(x = lon, y = lat), data = latlon), extent = "device")

#Plota o resultado (os pontos dos endereços)
sp.map.2012 + geom_point(size = I(4),colour="red", alpha = 2/3)

Pronto, o mapa final é esse aqui:

plot_zoom

Mas o assunto sobre dados georreferenciados e mapas está longe de ser esgotado aqui…

Webscrapping II – Baixando tabelas com o readHTMLTable

Neste segundo post sobre webscrapping, vamos apresentar um modo mais simples de extrair informações de tabelas em sites.

No exemplo anterior, acessamos o código fonte da página e copiamos seu texto através da função readLines. Desta vez, utilizaremos um comando do pacote XML, que é bastante utilizado para coleta de dados na internet com o R.

A linguagem HTML, básica para a escrita do conteúdo web, é organizada de forma “hierárquica”. No código-fonte, todo conteúdo de uma página se localiza entre as expressões: <HTML> e </HTML>. Essas duas expressões marcam o início e o fim de uma área dentro da qual todas as demais informações estarão contidas. Podemos entendê-las como demarcando as fronteiras de um grande conjunto. O conteúdo exibido no corpo principal da página é delimitado por <BODY> e </BODY>. Dentro desse “sub-conjunto”, é possível criar um parágrafo, cujo início e fim são assinalados por <p> e </p>. Vejam esse exemplo aqui.

Ou seja, todo conteúdo HTML se localiza dentro de “nós” ou nodes. Há um nó principal que contém todos os demais — e que marca o início e o fim de qualquer site. O que nos interessa é que o conteúdo de uma tabela é delimitado por <TABLE> e </TABLE>. Neste link, vocês encontram exemplos.

Algumas funções do pacote XML conseguem compreender essa estrutura da linguagem HTML e acessar o conteúdo dos nodes determinados. A função readHTMLTable é especifica para tabelas.

Vejamos como podemos simplificar código do exemplo anterior, através do uso desse comando:

require(XML)

# Endereço básico
baseurl <- "http://www.portaldatransparencia.gov.br/servidores/OrgaoExercicio-ListaServidores.asp?CodOS=25201&DescOS=BANCO%20CENTRAL%20DO%20BRASIL&CodOrg=25201&DescOrg=BANCO%20CENTRAL%20DO%20BRASIL&Pagina=XX"

#Loop para fazer download de uma sequencia de páginas
data <- data.frame()
for (i in 1:278) { #contador: vai da página 1 à pág. 278
  print(i) #imprime na tela a página que está sendo acessada
  url <- gsub("XX", i, baseurl) #substitui a expressão XX no endereço pelo índice da página
  x <- readHTMLTable(url)[[2]]

  data <- rbind(data, x) #salva os resultados num vetor
}

data #resultado final

O comando readHTMLTable extrai todas as tabelas da página. Aplicando um índice depois do comando, é possível acessar apenas a tabela de número desejado. Nesse caso, desejávamos a tabela de número “2”, por isso, readHTMLTable(url)[[2]].

Reparem que não foi preciso aplicar funções para “limpar” o conteúdo e retirar caracteres especiais.

Bem mais simples agora.

Webscrapping I: Fazendo download de dados de tabelas em sites

[POST ATUALIZADO EM 03/09/2013 – Correções apontadas por Thiago Silva] Webscrapping é a tarefa de coletar dados pela internet e prepará-los para a análise de forma automatizada. Dados de diversos tipos: tabelas/bancos de dados, textos, mapas, informações geocodificadas, softwares… Neste primeiro post sobre o assunto, trago um exemplo bastante simples: baixar dados de uma tabela e transformá-los em um banco de dados em R. O script abaixo foi produzido por uma amiga há bastante tempo — e levemente modificado por mim, depois. O site-alvo é o seguinte: http://www.portaldatransparencia.gov.br/servidores/OrgaoExercicio-ListaServidores.asp?CodOS=25201&DescOS=BANCO%20CENTRAL%20DO%20BRASIL&CodOrg=25201&DescOrg=BANCO%20CENTRAL%20DO%20BRASIL&Pagina=1 A consulta gerou 278 páginas — 278 tabelas! Seria possível, claro, copiar e colar o resultado de todas as páginas em planilhas e formatá-las manualmente. Mas que trabalhinho chato! tabela A primeira dica é: percebam o padrão contido no endereço da página, reparem que ele termina com: …&Pagina=1 . Se substituirmos isso por &Pagina=2, somos direcionados para a próxima lista de resultados. Então o primeiro passo é tentar automatizar essa mudança de páginas. Em segundo lugar, as informações contidas na tabela são, na realidade, um texto HTML formatado. Ou seja, de alguma forma, estão contidas dentro do código-fonte do site. Para visualizar isso, basta clicar com o botão direito em qualquer lugar da página e selecionar a opção “Ver código fonte” (ou algo parecido, a depender do navegador). E lá está! Na linha 285: tabela Então como fazemos? Guardamos num objeto o endereço básico da HTML, substituindo o número-índice da página (&Pagina=1, &Pagina=2 etc.) por XX: 

# Endereço básico</span>
baseurl <- "http://www.portaldatransparencia.gov.br/servidores/OrgaoExercicio-ListaServidores.asp?CodOS=25201&DescOS=BANCO%20CENTRAL%20DO%20BRASIL&CodOrg=25201&DescOrg=BANCO%20CENTRAL%20DO%20BRASIL&Pagina=XX"

Depois, criamos uma função simples para retirar caracteres irrelevantes presentes no código-fonte que vai ser baixado: 

# Cria um função simples para retirar tabs e carateres especiais da página
cleanupHTML <- function(z) {
  z <- gsub("\t","",z) #retira tabs
  z <- gsub("  ", "",z) #retira espaços duplos
  z <- gsub("^\\s+", "", z) # retira espaços no início das linhas
  z <- gsub("\\s+$", "", z) # retira espaços no fim das linhas
  return(gsub("<(.|\n)*?>","",z)) #retira caracteres especiais e retorna o resultado
}

Agora, para acessar todas as 278 páginas, baixar o código-fonte, aplicar a função que “limpa” e retira os caracteres indesejados em todo material, podemos criar um loop. Esse loop substitui a expressão “XX” pelo número da página, acessa o endereço, copia o conteúdo num objeto de tipo texto (string / character), seleciona apenas as linhas que contém informações relevantes (onde estão os dados da tabela), aplica a função de limpeza e guarda os resultados:

#Loop para fazer download de uma sequencia de páginas
data <- NULL
    for (i in 1:278) { #contador: vai da página 1 à pág. 278
    print(i) #imprime na tela a página que está sendo acessada
    url <- gsub("XX", i, baseurl) #substitui a expressão XX no endereço pelo índice da página
    x <- readLines(url) #lê o código fonte da página número 'i'
    x <- iconv(x) #converte o sistema de codificação de caracteres para UTF-8
    # Este último comando funciona apenas para Windows.
    # No Linux, troque pelo seguinte comando:
    # x <- iconv(x, "latin1", "utf-8")

    #delimina as linhas inicial e final (onde estão contidos os dados desejados)
    start <- grep("<table summary=\"Lista de servidores por Nível da Função\">", x)
    end <- grep("<script type=\"text/javascript\" language=\"javascript\" src=\"../../Recursos/helper.js\"></script>", x)[1]

    x <- cleanupHTML(x[start:end]) #aplica a função que retira tabs e caracteres especiais
    data <- c(data, x) #salva os resultados num vetor
 }

O que produzimos até agora é um vetor de várias linhas, que contém, em cada uma delas, não somente os dados, como também alguns espaços vazios, células em branco e “restos” que não foram limpados pela função que montamos. Por isso, um trabalhinho adicional ainda é necessário. Além disso, em cada uma das páginas de resultado, o cabeçalho da tabela se repete. É melhor retirá-lo completamente e depois re-adicioná-lo, mas já como nome das variáveis do banco de dados que vamos compor:


temp <- data

#Retira espaços vazious e marcas de comentários
 temp <- temp[-which(temp %in% c("", " ", " ", " ", " ", " ", " "," ", " <!--", " -->"))]

#Retira cabeçalhos (pois foram coletados várias vezes -- uma vez em cada página)
 temp <- temp[-which(temp %in% c("CPF", "Nome do servidor", "Órgão de lotação"))]

#organiza os resultados como um Dataframe
 temp_m <- as.data.frame(matrix(temp, byrow=T, ncol=3))

#Dá nome às variáveis
 names(temp_m) = c("CPF", "Nome do servidor", "Órgão de lotação")

View(temp_m)

O resultado é uma banco de dados (dataframe) como este: tabela Existem outros meios de fazer esse mesmo procedimento no R (alguns mais simples, inclusive). Pacotes como o XML, RJSONIO, rjson, RCurl e RHTMLForms possuem funções que facilitam muito o webscrapping. No script acima usamos apenas as funções contidas nos pacotes básicos, que já vem com o R.

R+Google Ngram Viewer = Análise de milhões de livros

O Google Ngram Viewer é uma ferramenta pouco conhecida. O que ela faz é simplesmente contar palavras, mas em milhões de livros ao mesmo tempo — que são parte do acerco do Google Books. Ou seja, faz, de modo quase instantâneo, um procedimento simples de Análise de Conteúdo.

Os seus criadores são pessoas de diversas áreas de formação e que pretendem dar início a um estudo massivo sobre cultura e comportamento através de dados quantitativos presentes na internet, iniciando pelo conteúdo dos livros. Chamam essa iniciativa de CulturonomicsNão sei se é possível dizer que já foi de todo bem sucedida, mas fato é que possuem artigo na Science e na Nature, ilustrando formas de uso da ferramenta e algumas consequências analíticas.

Mas até recentemente, o uso desses dados por parte da maioria dos usuários era feito via consultas no próprio site ou então baixando uma quantidade muito massiva de dados, partidos em inúmeros arquivos… Mas recentemente foi lançado um pacote para R que acessa, busca e retorna resultados do Ngram Viewer, permitindo, depois, analisá-los normalmente com as demais ferramentas e pacotes do R.

O código abaixo instala o pacote ngramr e faz uma busca pelos termos “science” e “religion” no corpus de livros em inglês publicados entre 1800 e 2008. Como o buscador é “case-sensitive”, vários formatos de escrita foram utilizados.

#Limpando o ambiente
rm(list=ls())

# Instalando o ngramr
install.packages("ngramr")

# Carregando os pacotes necessários (disponíveis no CRAN)
require(tseries)
require(ngramr)
require(ggplot2)
#########################################################

#Delimitando um período inicial e final para a pesquisa
start = 1800
end=2008

# Faz as pesquisas no Google Ngram Viewer e salva no objeto ng
ng = ngram(c("science+Science+SCIENCE+sciences+Sciences+SCIENCES",
      "religion+Religion+RELIGION+religions+Religions+RELIGIONS"),
      year_start = start,
      year_end=end,
      smoothing=1)

# Gráfico dos resultados
ggplot(ng, aes(x=Year, y=Frequency, colour=Phrase)) +
      geom_line(size=1.1) + opts(legend.position="bottom")

O resultado é um gráfico como esse:

Rplot

Ora… os resultados são como uma série temporal. Podemos aplicar técnicas específicas pra isso então.

# Transformando o objeto em Dataframe...
ng = as.data.frame(ng)

# Transformando as variáveis em um objeto do tipo Time Series
science = ts(ng[as.numeric(ng$Phrase)==1,]$Frequency,
      start=c(start,1),end=c(end,1),freq=1)

religion = ts(ng[as.numeric(ng$Phrase)==2,]$Frequency,
      start=c(start,1),end=c(end,1),freq=1)

######################################################
# Fazendo algumas análises:

# Aparentemente, as duas séries estão bastante correlacionadas
# O que visualmente percebemos, quando observamos as tendências
# divergentes, crescente e decrescente.

# O coeficiente de pearson corrobora isso:
cor(science,religion)
[1] -0.8703351

# Mas o Dickey-Fuller mostra que as duas séries não são estacionárias
adf.test(science)
adf.test(religion)

# Mas quando first differenced, são...
adf.test(diff(science))
adf.test(diff(religion))

# Fazendo a correlação das variáveis diferenciadas:
cor(diff(science),diff(religion))
[1] -0.05208317

É… o problema “ciência-religião” não é tão trivial… 

Ferramenta legal, não? Dá pra pensar em muitos usos — e as sugestões dos seus criadores já são bastante amplas.