All Posts By

Vitor Ostrensky

Extraindo dados eleitorais no R

By | Notas

O pacote electionsbr oferece um conjunto de funções que possibilitam a extração de dados eleitorais diretamente da base de dados do Tribunal Superior Eleitoral (TSE). Com ele é possível obter dados de todas as eleições desde 1994 desagregados até por zona eleitoral. Além disso, é possível obter informações sobre filiação partidária, características dos eleitores de cada zona eleitoral e declaração patrimonial de todos os candidatos.

Para mostrar a funcionalidade do pacote, iremos utilizar os dados da última eleição federal, em 2018. Iremos separar apenas os votos para presidente e selecionar o partido - no caso a chapa - que recebeu mais votos em cada município.


elec_2018 = party_mun_zone_fed(2018, br_archive = TRUE) %>%
filter(DESCRICAO_CARGO == "Presidente" & NUM_TURNO == 1) %>%
group_by(CODIGO_MUNICIPIO) %>%
top_n(1, QTDE_VOTOS_NOMINAIS) %>%
mutate(CODIGO_MUNICIPIO = as.double(CODIGO_MUNICIPIO)) %>%
select(CODIGO_MUNICIPIO, SIGLA_PARTIDO)

O código utilizado como identificador de cada município é específico do TSE e é diferente do código IBGE padrão. Por isso, como iremos colocar em um mapa utilizando o pacote geobr, precisamos de uma tabela para conversão.

</pre>
tse_ibge_id <- read_csv("https://raw.githubusercontent.com/betafcc/Municipios-Brasileiros-TSE/master/municipios_brasileiros_tse.csv")

elec_2018 = left_join(elec_2018, tse_ibge_id, by = c("CODIGO_MUNICIPIO" = "codigo_tse"))
<pre>

 

Assim, fazemos o download dos mapas por município e estado e juntamos com a nossa base.

</pre>
mapa_muni = read_municipality(showProgress = FALSE)
mapa_state = read_state(showProgress = FALSE)

df <- left_join(mapa_muni, elec_2018, by = c("code_muni" = "codigo_ibge"))

Agora, podemos fazer o mapa com os dados, que mostram grande concentração espacial do voto.

 


ggplot() +
geom_sf(df, mapping = aes(fill = SIGLA_PARTIDO), colour = NA) +
geom_sf(mapa_state, mapping = aes(fill = NA)) +
scale_fill_manual(values=c( "#fff200", "#263272", "#c4122d")) +
theme_minimal() +
labs(title = "Eleição em 2018",
subtitle = "Partido mais votado por município na eleição presidencial em 2018")

Análise estatística espacial no R (parte 2)

By | Data Science, Economia

Semana passada mostramos como estimar a dependência espacial de variáveis no R. Hoje, iremos continuar a falar sobre econometria espacial, apresentando como rodar modelos de regressão que levem em consideração o espaço. Para isso, utilizaremos um exemplo simples em que estimamos a relação  da taxa de homicídios com o PIB per capita.

O argumento mais importante para se estimar modelos de econometria espacial é se a suposição de independência entre as observações não é mais válida. Assim, características de um local podem influenciar nos resultados de outro.


library(geobr)
library(sidrar)
library(spdep)

library(stargazer)

set.ZeroPolicyOption(TRUE)

Utilizaremos basicamente o mesmo processo de extração de dados do post anterior, com a utilização de microrregiões como unidade geográfica. A diferença será a inclusão da variável de com os dados de taxa de homicídio.

#Dados de PIB
pib = get_sidra(5938,
variable = 37, 
geo = "MicroRegion")
#Dados de população
pop= get_sidra(202,
variable = 93,
geo = "MicroRegion") %>%
filter(`Sexo (Código)` == 0 & `Situação do domicílio (Código)` == 0)

#Juntando e criando a variável de pib per capita
df <- left_join(pib, pop, by = "Microrregião Geográfica (Código)") %>%
mutate(pibpc = (Valor.x/Valor.y) - mean(Valor.x/Valor.y)) #centralizar no 0

# colocando a taxa de homicídios
tx_homic <- ipeadatar::ipeadata("THOMIC") %>%
filter(uname == "Microregions" &
date == "2018-01-01") %>%
mutate("Microrregião Geográfica (Código)" = as.character(tcode)) %>%
rename(tx_homic = value) %>%
select(tx_homic, "Microrregião Geográfica (Código)" )

df <- left_join(df, tx_homic, by = "Microrregião Geográfica (Código)")

# mapa
mapa_micro = read_micro_region(showProgress = F)
mapa_micro$code_micro = as.character(mapa_micro$code_micro)

merged = dplyr::left_join(mapa_micro, df,
by = c("code_micro" = "Microrregião Geográfica (Código)"))

 

Também iremos produzir a matriz de vizinhança da mesma forma que antes. Iremos retirar da análise a microrregião de Fernando de Noronha, pois ela não apresenta vizinhos, o que impossibilita a estimação.

nb <- poly2nb(merged, queen=TRUE)
lw <- nb2listw(nb, style="W", zero.policy=TRUE)
df = df  %>%  filter(`Microrregião Geográfica.x` != "Fernando de Noronha - PE")

 

Existem vários modelos que se encaixam nesse arcabouço, mas os dois principais são o Autorregressivo Espacial (SAR) e o de Erros Espaciais (SEM). Ambos são modificações do modelo clássico OLS:

(1)   \begin{equation*} y=X \beta+\epsilon \end{equation*}

O modelo SAR introduz um processo autoregressivo, semelhante ao que ocorre nos modelos de séries temporais:

(2)    \begin{equation*} y=\rho W y+X \beta+\epsilon \end{equation*}

Já o modelo SEM introduz o componente do processo autoregressivo no termo de erro:

(3)    \begin{equation*} \begin{gathered} y=X \beta+\epsilon \\ \epsilon=\lambda W \epsilon+u \end{gathered} \end{equation*}

Para determinar qual dos dois modelos rodar, é preciso estimar o teste do tipo multiplicador de Lagrange para avaliar se o Lambda e o W são diferentes de zero.

ols <- lm(tx_homic ~ pibpc, data = df)

lmLMtests <- lm.LMtests(ols,lw, test=c("RLMerr", "RLMlag"))
lmLMtests

Como em ambos os casos a hipótese nula foi rejeitada, o procedimento padrão é estimar o modelo que tenha retornado a estatística mais significativa. No nosso caso, estimaremos o modelo SAR.

 

sar = lagsarlm(tx_homic ~ pibpc, data = df, lw, zero.policy = TRUE)
stargazer(ols,sar)

Veja que na especificação OLS temos que a há uma relação significativa entre as duas variáveis. Quando há a inclusão do componente espacial, o parâmetro do PIB per capita deixa de ser significativo. Isso mostra que, nesse caso, estimar o modelo OLS pode gerar o chamado viés da variável omitida.

 \begin{table}[!htbp] \centering \caption{} \label{} \begin{tabular}{@{\extracolsep{5pt}}lcc} \\[-1.8ex]\hline \hline \\[-1.8ex] & \multicolumn{2}{c}{\textit{Dependent variable:}} \\ \cline{2-3} \\[-1.8ex] & \multicolumn{2}{c}{tx\_homic} \\ \\[-1.8ex] & \textit{OLS} & \textit{spatial} \\ & \textit{} & \textit{autoregressive} \\ \\[-1.8ex] & (1) & (2)\\ \hline \\[-1.8ex] pibpc & $-$0.128$^{***}$ & $-$0.017 \\ & (0.038) & (0.033) \\ & & \\ Constant & 25.402$^{***}$ & 11.914$^{***}$ \\ & (0.740) & (1.380) \\ & & \\ \hline \\[-1.8ex] Observations & 557 & 557 \\ R$^{2}$ & 0.020 & \\ Adjusted R$^{2}$ & 0.018 & \\ Log Likelihood & & $-$2,322.888 \\ $\sigma^{2}$ & & 230.258 \\ Akaike Inf. Crit. & & 4,653.776 \\ Residual Std. Error & 17.460 (df = 555) & \\ F Statistic & 11.043$^{***}$ (df = 1; 555) & \\ Wald Test & & 126.609$^{***}$ (df = 1) \\ LR Test & & 118.879$^{***}$ (df = 1) \\ \hline \hline \\[-1.8ex] \textit{Note:} & \multicolumn{2}{r}{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} \\ \end{tabular} \end{table}

 

 

Análise estatística espacial no R

By | Data Science

Trabalhos aplicado em economia regional e urbana dependem muito de dados com referência a locais medidos como pontos em espaço, o que gera consequências importantes na análise estatística. Isso acontece pois há dependência espacial entre as observações e pois há heterogeneidade nas relações que estamos modelando. Modelos que levam em conta esses fatores entram no campo de estatística/econometria espacial.

Nessa área é muito comum a utilização de softwares especializados, como o GeoDa. Entretanto, vamos mostrar que é possível realizar esse tipo de operação no R. No post de hoje iremos apresentar como fazer a estimação do Índice de Moran, que permite fazer a clusterização de observações no espaço. Na semana que vem mostraremos como estimar modelos de econometria espacial. Para isso, utilizaremos o pacote spdep.

library(tidyverse)
library(geobr)
library(sidrar)
library(spdep) 

Para esse exemplo, iremos usar dados de PIB por microrregião, obtidos pelo pacote sidrar. Veja que estamos centralizando a variável em 0. Então, utilizamos o pacote geobr para fazer o download do mapa e o unimos com os dados.

 

#Dados de PIB
pib = get_sidra(5938,
variable = 37, 
geo = "MicroRegion")</pre>
#Dados de população
pop= get_sidra(202,
variable = 93,
geo = "MicroRegion") %>%
filter(`Sexo (Código)` == 0 & `Situação do domicílio (Código)` == 0)

#Juntando e criando a variável de pib per capita
df <- left_join(pib, pop, by = "Microrregião Geográfica (Código)") %>%
mutate(pibpc = (Valor.x/Valor.y) - mean(Valor.x/Valor.y)) #centralizar no 0

#
mapa_micro = read_micro_region()
mapa_micro$code_micro = as.character(mapa_micro$code_micro)

merged = dplyr::left_join(mapa_micro, df, by = c("code_micro" = "Microrregião Geográfica (Código)"))

 

Para estimar os modelos da área de econometria espacial é preciso construir a chamada matriz de vizinhança. Essa matriz tem dimensão n x X, sendo n o número de unidades geográficas. No caso deste exemplo, ela mostra quais microrregiões tem vizinhança com quais. Assim, para aquelas em que há vizinhança, há um valor 1. Para todas as outras, 0.

Em sequência é preciso transforma-la em uma matriz de pesos. Uma microrregião com mais vizinhos tem menor peso para cada um deles. Já uma unidade geográfica que tenha apenas 1 vizinho terá peso 1 para este. Assim, utilizamos essa matriz de peso para estimar o chamado I de Moran, que indica o grau de associação espacial da variável que estamos analisando.


nb <- poly2nb(merged, queen=TRUE)
lw <- nb2listw(nb, style="W", zero.policy=TRUE)

moran <- localmoran(merged$pibpc, lw)

ggplot(merged) +
geom_point(aes(x = pibpc, y =lag_pibpc)) +
theme_minimal() +
xlab("Pib per capita") +
ylab("Pib per capita medio dos vizinhos")

Com isso, separamos as observações com relação estatística significante em quatro grupos a depender se elas estão acima ou abaixo da média de renda e se seus vizinhos estão a baixo e acima da média.

merged$moransig <- moran[,5] #p valor
merged$lag_pibpc<- lag.listw(lw, merged$pibpc)</pre>
merged <- merged %>%
mutate(quadrante = case_when(pibpc >= 0 & lag_pibpc >= 0 & moransig < 0.05~"Alto-Alto",
pibpc <= 0 & lag_pibpc <= 0 & moransig < 0.05~"Baixo-Baixo",
pibpc >= 0 & lag_pibpc <= 0 & moransig < 0.05~"Alto-Baixo",
pibpc <= 0 & lag_pibpc >= 0 & moransig < 0.05~"Baixo-Alto"))
<pre>

Assim, podemos plotar o mapa dos clusters que estimamos. Veja que apenas dois tipos de clusters estão representados. Isso significa que não há microrregiões ricas com vizinhos pobres ou pobres com vizinhos ricos.


ggplot(merged) +
geom_sf(aes(fill = quadrante), color=NA) +
scale_fill_manual(values = c("blue", "red"), na.value = "lightgrey") +
theme_minimal()


Importando planilhas mal formatadas no R

By | Data Science, Hackeando o R

Quando trabalhamos em uma equipe, nem sempre podemos contar que todos usem as mesmas ferramentas. No ramo de dados isso pode ser um complicador importante, dado a necessidade de padrões que possam ser utilizado por diversos meios. Assim, quem está acostumado a utilizar R para análises deve estar habituado com o padrão tidy, que segue três regras principais: 1) Cada variável é uma coluna; 2) Cada observação é uma linha e 3) Cada célula tem apenas uma medida. No entanto, quando esse padrão não é seguido, pode ser difícil seguir com o trabalho. Frequentemente, as pessoas utilizam o Excel para criar análises, gerando planilhas muito elaboradas, mas que se tornam impossíveis de serem importadas por qualquer outra ferramenta.

Os pacotes tidyxl e unpivotr cobrem essa dificuldade, facilitando an importação de planilhas mal estruturadas no R. O primeiro serve para ler a planilha de excel de maneira diferente do usual, célula por célula. O segundo serve para estruturar essa informação em uma tabela em formato tidy.

Utilizaremos um exemplo simples para mostrar as funcionalidades destes pacotes. Veja que a mesma planilha traz diversas informações financeiras para três produtos, dois clientes e dois períodos, que são arranjados em duas tabelas diferentes.

 


library(tidyverse)
library(readxl)
library(tidyxl)
library(unpivotr)
library(knitr)

Veja que se tentarmos ler o arquivo pelo modo tradicional, pelo pacote readxl, não há possibilidade nenhuma de prosseguirmos.


kable(readxl::read_xlsx("exemplo.xlsx")[1:5,1:5])


 

APÊNDICE II ...2 ...3 ...4 ...5
VENDAS TOTAIS DA EMPRESA NA NA NA NA
NA NA NA NA NA
NA Período 1 VENDAS NA NA
NA NA Vendas Vendas Faturamento Bruto
NA NA Peso Unidades NA

 

 

Agora importando pelo pacote tidyxl


cells = xlsx_cells("exemplo.xlsx")

A partir daqui a funções utilizadas dependem muito da estrutura da planilha. No caso estamos filtrando todas as células em branco. Além disso, estamos retirando as linhas iniciais e as linhas apenas com os nomes das colunas na segunda tabela. Isso vai ser útil para trabalharmos apenas com uma tabela.


extracted <-
cells %>%
filter(!is_blank,
row > 4 & !(row %in% c(16,17,18,19))) %>%
select(row, col, data_type, numeric, character)

Agora utilizaremos o comando "behead" do unpivotr para extraírmos as colunas. O comando é dividido em duas partes. Na primeira dizemos a direção em que está o título da coluna/linha e na segundo o nome que daremos a esta coluna. Na primeira a direção é para cima ("up"), pois é onde está a linha com os nomes das colunas (Vendas, Faturamento, etc.). Na segunda vamos de novo para cima, para selecionar a linha que diferencia as duas colunas de vendas, entre peso e unidade. Na terceira, a direção é esquerda-baixo ("left-down"). Esse é um caso importante, pois estamos tratando de células mescladas. Neste caso, os títulos estão à esquerda e abaixo os valores que ela descreve. Por último, criamos a coluna de produtos, também à esquerda.

Além disso, criamos uma coluna com o período. Se a linha for menor ou igual a 15, as informações são sobre o período 1, caso contrário, do período 2. Também juntamos a coluna de variável com a de unidade, para os casos das variáveis de vendas, onde há a separação. Perceba que agora temos uma tabela "tidy".

tidied <- extracted %>%
behead("up", "variable") %>%
behead("up", "tipo") %>%
behead("left-down", "cliente") %>%
behead("left", "produto") %>%
mutate(periodo = ifelse(row <= 15, 1,2),
variable = ifelse(is.na(tipo), variable, paste0(variable,"_",tipo))) %>%
select(variable, cliente, periodo, produto, numeric) %>%
drop_na(cliente)

kable(head(tidied,10))

variable cliente periodo produto numeric
Vendas_Peso Cliente A 1 Produto 1 2
Vendas_Unidades Cliente A 1 Produto 1 4
Faturamento Bruto (em R$) Cliente A 1 Produto 1 50
IPI Cliente A 1 Produto 1 8
ICMS Cliente A 1 Produto 1 10
PIS Cliente A 1 Produto 1 14
COFINS Cliente A 1 Produto 1 16
Total de Impostos Cliente A 1 Produto 1 48
Receita Operacional Líquida (R$) Cliente A 1 Produto 1 2
Vendas_Peso Cliente B 1 Produto 2 1

 

 

Para torna-la no formato semelhante ao da tabela inicial, basta utilizarmos a função pivot_wider do tidyr.



kable(tidied %>%
pivot_wider(names_from = variable, values_from = numeric) %>%
select(1:6))

cliente periodo produto Vendas_Peso Vendas_Unidades Faturamento Bruto (em R$)
Cliente A 1 Produto 1 2 4 50
Cliente B 1 Produto 2 1 2 30
Cliente B 1 Produto 3 1 2 20
Cliente B 1 Total (A) 4 8 100
Cliente B 1 Produto 1 2 4 40
Total (A) + (B) 1 Produto 2 1 2 35
Total (A) + (B) 1 Produto 3 1 2 25
Total (A) + (B) 1 Total (B) 4 8 100
Total (A) + (B) 1 NA 8 16 200
Cliente A 2 Produto 1 2 4 50
Cliente B 2 Produto 2 1 2 30
Cliente B 2 Produto 3 1 2 20
Cliente B 2 Total (A) 4 8 100
Cliente B 2 Produto 1 2 4 40
Total (A) + (B) 2 Produto 2 1 2 35
Total (A) + (B) 2 Produto 3 1 2 25
Total (A) + (B) 2 Total (B) 4 8 100
Total (A) + (B) 2 NA 8 16 200

 

Utilizando o pacote basedosdados no R

By | Data Science, R News

Neste post vamos mostrar como utilizar o recém lançado pacote basedosdados, que fornece um jeito simples para acessar o datalake da organização Base dos Dados. São centenas de dados disponíveis, já tratados e de fácil compatibilização entre si. Entre as bases disponíveis nesse datalake estão RAIS, CAGED, comércio exterior, dados eleitorais e dados de CNPJ.

Para acessar os dados é necessário ter uma conta e um projeto no Google Cloud. Assim, tendo um projeto, é preciso colocar a sua chave identificadora utilizando a função "set_billing_id".

 

library(tidyverse)
library(basedosdados)

basedosdados::set_billing_id(XXXXXX) # trocar para o seu identificador

Iremos utilizar um exemplo simples para mostrar a facilidade do cruzamento de dados. Nota-se que o objetivo não é fazer nenhum tipo de inferência, mas apenas mostrar a funcionalidade do pacote.  Iremos cruzar três variáveis a nível municipal: População, valor adicionado pela indústria e óbitos por doenças respiratórias (CID-J). As duas primeiras são provenientes do IBGE, já a última vem do Sistema de Informações sobre Mortalidade (SIM), do Datasus. Cada "query", ou seja, a seleção dos dados, é feita por meio de SQL.

#óbitos
query1 <- "SELECT ano, id_municipio, SUM(numero_obitos) AS obitos
FROM `basedosdados.br_ms_sim.municipio_causa`
WHERE LEFT(causa_basica,1) = 'J' # Apenas doenças respiratórias
GROUP BY ano, id_municipio"

obitos <- read_sql(query1)

#população
query2 <- "SELECT *
FROM `basedosdados.br_ibge_populacao.municipios`"

pop <- read_sql(query2)

#PIB
query3 <- "SELECT id_municipio, ano, VA_industria
FROM `basedosdados.br_ibge_pib.municipios`"

pib <- read_sql(query3)

Uma grande facilidade trazida pelo Base dos Dados é fornecer centralização e padronização. Por exemplo, nesse caso, podemos juntar as três tabelas pelo Código IBGE de cada município e pelo ano, que já estão com o mesmo nome em todas elas. Quem já trabalhou com dados municipais sabe que os identificadores dos municípios podem estar em formatos diferentes ou até não estarem presentes, dificultando bastante o tratamento dos dados.

Assim, juntando os três data frames e filtrando apenas para valores de 2018, podemos mostrar a relação entre óbitos por doenças respiratórias e o valor per capita adicionado pela indústria.

 

df <- left_join(obitos, pop, by = c("id_municipio", "ano"))
df <- left_join(df, pib, by = c("id_municipio", "ano"))

df_18 <- df %>%
filter(ano == 2018) %>%
mutate(obitos_pc = obitos*100000/populacao,
industria_pc = VA_industria/populacao,
log_industria_pc = log(industria_pc),
pc = predict(prcomp(~log_industria_pc+obitos_pc, .))[,1])


ggplot(data = df_18, aes(x = log_industria_pc, y = obitos_pc, color = pc)) +
geom_point(show.legend = FALSE, shape = 16, size = 2, alpha = .5) +
theme_minimal() +
ylab("Óbitos/100 mil habitantes por doenças respiratórias") +
xlab("Valor adicionado pela indústria (em log)") +
scale_color_gradient(low = "#0091ff", high = "#f0650e")


 

 

Conheça o Curso de Avaliação de Políticas Públicas usando o R

Receba diretamente em seu e-mail gratuitamente nossas promoções especiais
e conteúdos exclusivos sobre Análise de Dados!

Assinar Gratuitamente
{"cart_token":"","hash":"","cart_data":""}