Como executar scripts demorados em segundo plano no R

By | Data Science

Algumas tarefas possuem um custo computacional maior do que outras, como, por exemplo, renderizar um gráfico de alta resolução, executar um loop ou estimar modelos mais robustos. Para estes casos é útil executar o código em "segundo plano" mantendo o Console do RStudio livre para trabalhar. Neste exercício demonstramos como fazer isso utilizando o pacote job.

Exemplo 1: renderização de gráficos

A depender da quantidade de dados, a renderização de um gráfico pode se tornar muito lenta. Vamos simular a renderização de um gráfico do ggplot2 com um dataset de 1 milhão de linhas:

# Dataset simulado
n <- 1e6
dataset <- data.frame(x = seq_len(n), y = rnorm(n))
# Criar objeto com gráfico ggplot2
grafico <- ggplot2::ggplot(dataset, ggplot2::aes(x, y)) +
ggplot2::geom_point(size = 0.1, alpha = 0.05)

Com os dados e a "gramática" do gráfico criados, podemos executar sua renderização como uma tarefa (job) na aba Jobs do RStudio. Neste caso iremos realizar a tarefa de renderizar o gráfico salvando o mesmo como um arquivo de imagem no computador, bastando passar o código dentro da função job::job({script aqui}):

# Instalar pacote
if(!require("job")) install.packages("job")
# Renderizar gráfico em um job
job::job({
ggplot2::ggsave("grafico_pontos.png", plot = grafico, width = 7, height = 3)
})

Pode-se acompanhar a execução da tarefa em segundo plano na aba Jobs do RStudio.

Exemplo 2: modelo de regressão

Outro exemplo bastante comum é a estimação de modelos de regressão. Aqui iremos demonstrar a criação de jobs para dois modelos bayesiano simples, usando o pacote brms:

# Carregar pacote
if(!require("brms")) install.packages("brms")
library(brms)
# Preparar dados e criar fórmulas
dados <- mtcars[mtcars$hp > 100, ]
modelo_1 <- mpg ~ hp * wt
modelo_2 <- mpg ~ hp + wt
# Criar jobs para estimar os modelos
job::job({
fit_1 = brms::brm(modelo_1, dados) # job do modelo 1
})
job::job({
fit_2 = brms::brm(modelo_2, dados) # job do modelo 2
})
# Continuar trabalhando no Console
cat("Estou livre! \n\nAtenciosamente, \nConsole do RStudio.")
# Estou livre!
# Atenciosamente, 
# Console do RStudio.

Como visto, o pacote é um forte candidato a entrar no toolkit de trabalho de qualquer pessoa que utilize intensivamente o R e RStudio. Confira estes e outros exemplos na documentação do pacote.

Hackeando o R: visualizando o efeito de variáveis em um modelo linear

By | Hackeando o R

No Hackeando o R de hoje, vamos mostrar como fazer a visualização do impacto das variáveis de um modelo linear com o pacote Effects. Esse tipo de visualização é interessante para facilitar a comunicação de resultados estatísticos, garantindo a interpretação correta de seus modelos. Vamos iniciar nosso exemplo gerando um modelo linear usual:

library(car)

Prestige$type = factor(Prestige$type, levels=c("bc", "wc", "prof"))
lm1 = lm(prestige ~ education + poly(women, 2) +
log(income)*type, data=Prestige)

summary(lm1)

lm(formula = prestige ~ education + poly(women, 2) + log(income) * 
type, data = Prestige)

Residuals:
Min 1Q Median 3Q Max 
-12.1070 -3.8277 0.2736 3.8382 16.4393

Coefficients:
Estimate Std. Error t value Pr(&amp;gt; |t|) 
(Intercept) -137.5002 23.5219 -5.846 8.18e-08 ***
education 2.9588 0.5817 5.087 2.01e-06 ***
poly(women, 2)1 28.3395 10.1900 2.781 0.00661 ** 
poly(women, 2)2 12.5663 7.0954 1.771 0.07998 . 
log(income) 17.5135 2.9159 6.006 4.06e-08 ***
typewc 0.9695 39.4947 0.025 0.98047 
typeprof 74.2759 30.7357 2.417 0.01771 * 
log(income):typewc -0.4661 4.6200 -0.101 0.91986 
log(income):typeprof -7.6977 3.4512 -2.230 0.02823 * 
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.199 on 89 degrees of freedom
(4 observations deleted due to missingness)
Multiple R-squared: 0.8793, Adjusted R-squared: 0.8685 
F-statistic: 81.08 on 8 and 89 DF, p-value: &amp;lt; 2.2e-16

Dentre os regressores do modelo, apenas education possui uma interpretação direta, de que uma unidade adicional aumenta o valor de prestige em 2.95. Para as outras variáveis, temos efeitos que variam de magnitude, como no caso de women, e transformações de escala misturadas com interações, fazendo com que a compreensão do modelo não seja muito intuitiva. Para resolver isso, vamos utilizar a função plot() do pacote effects, que permite visualizar o efeito de uma das variáveis. Abaixo, o efeito de education:

library(effects)

e1.lm1 = predictorEffect("education", lm1)

plot(e1.lm1)

O gráfico gerado apresenta uma reta cuja angulação é o coeficiente do regressor no modelo, e o valor da função de efeito é prestige em função de education, com os outros regressores fixos em valores padrões, como a média deles, sendo assim o efeito parcial de education. A banda desenhada é o intervalo de confiança para a estimação desse valor, se baseando na matriz de covariâncias dos regressores da amostra. Para um parâmetro simples, não há grandes ganhos sobre a interpretação, porém no caso da variável income, que entra no modelo em logaritmo e tem interação com dummies, o efeito é mais complicado, e o gráfico se torna mais interessante:

plot(predictorEffect("income", lm1),
lines=list(multiline=TRUE))

 

No caso da própria variável type, que é categórica, o efeito depende da categoria, e do valor de income. Para entendermos como funciona o modelo em níveis distintos de income, são gerados pontos para os 5 principais quantis:

plot(predictorEffect("type", lm1, xlevels = 5), lines=list(multiline=TRUE))

 

 

Inflação no Atacado e o avanço das commodities

By | Comentário de Conjuntura

A inflação segue sendo uma preocupação importante para a conjuntura brasileira, a despeito de um hiato do produto ainda bastante aberto. O IPCA-15 de julho, por exemplo, teve avanço de 0,72%, reagindo, sobretudo, à alta da energia elétrica. Julho é um mês que, sazonalmente, a inflação mensal deveria começar a apresentar algum arrefecimento. A verdade, como temos discutido nesse espaço já há alguns meses, é que a inflação atual está longe de ser fruto de apenas um ou outro choque. Os núcleos e o índice de difusão mostram, de fato, que a inflação atual já é um processo bastante consolidado e que, se houve choques primários no passado, isso se reverberou sobre os demais preços da economia.

Um importante fruto de tensão para os preços ao consumidor, nesse aspecto, como notado pelo Diogo Wolff, no Relatório AM publicado ontem nesse espaço, vem do atacado. O Índice de Preços ao Produtor Amplo (IPA) tem mostrado uma inflação acumulada em 12 meses próxima a 50% nos últimos meses.

Por trás desse aumento consistente do IPA está um avanço forte nos preços das commodities em reais, seja porque houve de fato um aumento de diversas matérias-primas ao longo da pandemia, seja porque o câmbio R$/US$ se depreciou no período.

Tanto o IPA de produtos agrícolas quanto o IPA de produtos industrais apresentaram no período avanços consistentes, fruto em grande parte do aumento das commodities.

A reversão da inflação, por suposto, parece ainda bastante distante do cenário base, ao considerar toda essa pressão vinda do atacado. O hiato do produto aberto, diga-se, dificulta o repasse integral para o varejo, mas é ingênuo imaginar que isso não está ocorrendo.

Sobre isso, em breve, divulgaremos estudo sobre efeito de choques no atacado na inflação ao consumidor, no âmbito do Clube AM.

________________

Os códigos do exercício estarão disponíveis na plataforma do Clube AM.

Relatório AM #11 - IGP-M

By | Indicadores

O IGP-M, índice divulgado mensalmente pela FGV, é um importante indicador da inflação no país, pois também captura a inflação por parte da oferta, causada por pressões nos insumos da indústria. Como podemos ver abaixo, o IGP e suas desagregações estão tendo uma forte pressão desde o ano passado:

Os índices IGP são compostos por 3 sub-índices: o Índice de Preços por Atacado (IPA), o Índice de Preços ao Consumidor (IPC) e o Índice de Preços de Custo da Construção civil (INCC). O IGP-M dá a cada sub-índice pesos 0.6, 0.3 e 0.1, respectivamente. Nos gráficos acima, vemos que o IPC está em valores alinhados com o IPCA atual, e, como o INCC tem o menor peso, a alta do IGP deve ser explicada pela variação no IPA. De fato, podemos ver a abaixo que tal índice se descolou fortemente do IPCA nos últimos meses:

O que explica tais resultados? Há diversos fatores em pauta, como a escalada do câmbio no final do ano passado (que parece ter se estabilizado mais recentemente), aumentando o custo de empresas que dependem de importação de insumos, e o próprio aumento do valor de commodities, fatores essenciais para as diversas cadeias de produção.

Datação de recessões e ciclos econômicos no R com o algoritmo de Harding-Pagan

By | Artigos de Economia, Dados Macroeconômicos, Data Science, Economia

Ao longo do tempo a economia apresenta o que se chama de ciclos econômicos, ou seja, períodos de expansão e recessão. Mas de que forma podemos saber em qual ponto do ciclo econômico a economia se encontra? Como sabemos se a economia está em recessão? Estas são perguntas de grande interesse para acadêmicos e profissionais da área, e neste breve exercício demonstramos como replicar a datação de ciclos econômicos que instituições como NBER (EUA) e CODACE (Brasil) tradicionalmente publicam.

De maneira prática, neste exercício replicamos o algoritmo de Harding & Pagan (2002) para datar o ciclos de negócios do Produto Interno Bruto (PIB) brasileiro. Em resumo, o método considera algumas regras impostas ao comportamento de uma série temporal para classificar picos e vales. Recessão é o período entro o pico da atividade econômica e seu subsequente vale, ou ponto mínimo. Entre o vale e o pico, diz-se que a economia está em expansão.

O método é bastante simples e poderoso, conseguindo praticamente replicar a cronologia de recessões desenvolvidas pelas instituições mencionadas acima.

Pacotes

Para aplicar o algoritmo utilizaremos o pacote BCDating na linguagem R, criado por Majid Einian (Central Bank of Islamic Republic of Iran) e Franck Arnaud (National Institute of Statistics and Economic Studies, France). Outros pacotes são utilizados para coleta, tratamento e visualização de dados:


# Instalar/carregar pacotes
if(!require("pacman")) install.packages("pacman")
pacman::p_load(
"magrittr",
"BCDating",
"sidrar",
"dplyr",
"lubridate",
"timetk",
"zoo",
"ggplot2",
"ggthemes"
)

Dados

Neste exercício utilizaremos a série do PIB a preços de mercado (série encadeada do índice de volume trimestral com ajuste sazonal, média de 1995 = 100), disponível no SIDRA/IBGE. Para coletar os dados via API pode-se usar o pacote sidrar, especificando o código de coleta. Além disso realizamos a preparação dos dados para utilização posterior:


# Coleta e tratamento de dados
pib <- sidrar::get_sidra(api = "/t/1621/n1/all/v/all/p/all/c11255/90707/d/v584%202") %>%
dplyr::select("date" = `Trimestre (Código)`, "value" = `Valor`) %>%
dplyr::mutate(value = value, date = lubridate::yq(date)) %>%
dplyr::as_tibble()

Algoritmo de Harding & Pagan (2002)

Para aplicar o algoritmo e obter as datações de ciclo de negócios, primeiro transformamos o objeto pro formato de série temporal e, em seguida, utilizamos a função BBQ() do pacote BCDating. Optamos por deixar com os valores predefinidos os demais argumentos da função, que servem para definir os valores mínimos de duração do ciclo (pico ao pico ou vale ao vale) e da fase do ciclo (pico ao vale ou vale ao pico).


# Obter datação de ciclo de negócios
bc_dates <- pib %>%
timetk::tk_ts(select = value, start = c(1996, 1), frequency = 4) %>%
BCDating::BBQ(name = "Ciclo de Negócios do PIB do Brasil")

Resultados

Como pode ser visto abaixo, o objeto retornado traz como resultado as datas (trimestres) de picos e vales, assim como a duração do ciclo.


# Exibir resultados
show(bc_dates)

## Peaks Troughs Duration
## 1 2001Q1 2001Q4 3
## 2 2002Q4 2003Q2 2
## 3 2008Q3 2009Q1 2
## 4 2014Q1 2016Q4 11
## 5 2019Q4 2020Q2 2

Outras informações podem ser obtidas com a função summary(), porém o mais interessante é avaliar o resultado visualmente através de um gráfico. Para tal, fazemos um tratamento dos dados retornados pela função BBQ() e utilizamos o ggplot2 para gerar o gráfico com as áreas sombreadas referente às datas de recessão que foram identificadas pelo algoritmo, acompanhadas do comportamento do PIB no período:

Comparação com cronologia do CODACE/FGV

Por fim, vamos comparar os resultados aqui encontrados com a Cronologia de Ciclos de Negócios Brasileiros elaborada pelo Comitê de Datação de Ciclos Econômicos (CODACE). A última reunião do comitê foi em 29 de junho de 2020, na qual reportou a seguinte situação do ciclo de negócios:

Fonte: CODACE/FGV

Percebe-se que a série utilizada pelo comitê inicia-se em 1980, mas se analisarmos a partir de 1996 (período de início da série utilizada em nosso exercício), verificamos que 5 de 6 recessões datadas pelo CODACE são identificadas pelo algoritmo de Harding & Pagan (2002). Apenas a recessão do 1º trimestre de 1998 ao 1º trimestre de 1999 não foi detectada. Apesar disso, o resultado é empolgante!

Por fim, vamos comparar os resultados de ambas as datações mais a fundo na tabela a seguir, na qual contabilizamos o período de recessão partindo do trimestre imediatamente posterior ao pico até o subsequente vale:

Perceba que ambas as datações são idênticas! A única diferença está na última datação, a qual o CODACE ainda não definiu o próximo vale. Dessa forma, fica demonstrado o poder e facilidade de uso do algoritmo de Harding & Pagan para datação de ciclos econômicos.

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":""}