All Posts By

diogowolff

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(> |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: < 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))

 

 

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.

Hackeando o R: calculando o carry-over estatístico de uma variável

By | Hackeando o R

No Hackeando o R de hoje, vamos mostrar como fazer a visualização do carry-over estatístico de uma série temporal. Essa estatística pode ser uma importante ferramenta para a análise de agregações de dados ao longo do tempo pois facilita identificar a variação que ocorreu apenas no período corrente, extraindo a variação que é apenas resíduo do período anterior, chamada de carry-over. Abaixo, visualizamos os dois efeitos teóricos com uma simulação de dados

library(RcppRoll)
library(ggplot2)

set.seed(1)
x = data.frame(valor = (1:100) + rnorm(100), t = 1:100)

ggplot(x[49:72,], aes(x = t, y = valor))+geom_bar(stat = 'identity') +
coord_cartesian(ylim = c(45, 73))+
geom_segment(aes(x=49, xend=60, y = 54.68361, yend=54.68361), size = 1.2)+
geom_label(aes(x=49, y=54.68361, label = 'A'))+
geom_segment(aes(x=61, xend=72, y = 59.86495, yend=59.86495), size = 1.2)+
geom_label(aes(x=61, y=59.86495, label = 'B'))+
geom_segment(aes(x=61, xend=72, y = 66.85647, yend=66.85647), size = 1.2)+
geom_label(aes(x=61, y=66.85647, label = 'C'))+
labs(x='', y = '')+
theme_bw()

No exemplo acima, A é a média do ano anterior, C a média do ano corrente, e B é o valor da última observação do ano anterior, repetido para o ano corrente, ou seja, a média do ano corrente caso não houvesse crescimento. Ao compararmos a variação interanual dos dois períodos, podemos decompor esse valor em duas partes: a variação percentual de A a B, chamada de carry-over, e a variação percentual de B a C (mensurada no nível de A), que é o crescimento que ocorreu apenas a partir da última observação do ano anterior. A função abaixo calcula tais valores para uma variável mensal qualquer:


calcula_carry_over_anual <- function(data) {
A <- dplyr::lag(RcppRoll::roll_meanr(data, n=12), n=12)
B <- dplyr::lag(data, n = 12)
C <- RcppRoll::roll_meanr(data, n=12)

carry_over <- (B-A)/A
cresc_real_do_periodo <- (C-B)/A

lista = data.frame(carry_over, cresc_real_do_periodo, carry_over+cresc_real_do_periodo)
return(lista)
}

 

Então, vamos fazer a decomposição da série de nível do IBC como exemplo:

library(BETS)
library(tidyverse)
library(ggplot2)
library(scales)

ibc = BETSget(24363, data.frame=TRUE)

tibble(ibc$date, calcula_carry_over_anual((ibc$value))*100) %>%
magrittr::set_colnames(c('date', 'carry_over', 'cresc_real', 'soma')) %>%
pivot_longer(-date, names_to = 'var', values_to = 'val') %>%
filter(date>as.Date('2018-01-01') & var != 'soma') %>%
mutate(idk = RcppRoll::roll_sumr(val, n=2),
idk = ifelse(rep(c(FALSE, TRUE), times = 39), idk, NA)) %>%
ggplot(aes(x=date, y = val, fill = var))+geom_bar(stat = 'identity')+
scale_x_date(breaks = date_breaks('3 months'),
labels = date_format("%b/%Y"))+
scale_fill_manual(labels = c('Carry over', 'Crescimento real'), values = c('#244747', '#9ae5de'))+
geom_line(aes(x=date,y=idk, color = 'Agregado'), size= 1.2, linetype='solid')+
scale_color_manual(values = c('Agregado' = '#e89835'))+
geom_hline(yintercept=0, colour='black', linetype='dashed')+
labs(title='Decomposição da variação do nível do IBC', y = '%',
caption='Fonte: IBGE')+
theme(panel.background = element_rect(colour = 'white', fill='white'),
legend.position = 'right',
strip.text = element_text(size=8, face='bold'),
axis.text.x = element_text(angle = 45, hjust=1),
plot.title = element_text(size=10, face='bold'),
legend.title = element_blank(),
plot.caption.position = 'plot',
axis.title.x = element_blank())

Relatório AM #10 - IPCA-15

By | Indicadores

No Relatório AM dessa semana, vamos comentar brevemente o comportamento do IPCA-15, divulgado na metade de todo mês. Esse índice é importante pois é um previsor do IPCA cheio, logo suas características tendem a ser semelhantes. Uma delas é a sazonalidade, como podemos ver no gráfico abaixo:

Como podemos ver, o componente sazonal nos meses da metade do ano é baixo, indicando que as variações recentes que estamos vendo podem ser atribuídas em quase sua totalidade a choques exógenos que estão afetando a economia. Existem diversas fontes para tais choques, mas a mais importante do momento é o setor de energia, com a iminente crise hídrica que estamos sofrendo. Com isso, o IPCA-15 do último mês fechou em forte alta:

 

Hackeando o R: visualizando dados categóricos com mosaicos

By | Hackeando o R

No Hackeando o R de hoje, vamos mostrar métodos diferentes de visualizar dados distribuídos em categorias. Os dados utilizados no exemplo serão do dataset Titanic, disponível no R base.

ftable(Titanic)

A tabela acima é perfeitamente válida para acessarmos os dados conforme a necessidade. Apesar disso, a comparação entre linhas e suas subdivisões pode levar algum tempo, de modo que a criação de gráficos se justifica como método de facilitar a compreensão dos dados. Uma visualização inicial que podemos fazer é conferir o número de pessoas no navio por sexo. Para isso, um gráfico de barras simples é válido:


library(ggplot2)

df_titanic &lt;- as.data.frame(Titanic)

ggplot(df_titanic, aes(x=Sex, y= Freq)) + geom_bar(stat = 'identity') +
labs(x='Sexo', y = 'Número de pessoas') + theme_minimal()

A partir desse gráfico, uma expansão simples é dividir as pessoas entre quem sobreviveu ou não. Para isso, basta adicionar um fill:



ggplot(df_titanic, aes(x=Sex, y= Freq, fill = Survived)) + 
geom_bar(stat = 'identity', position = position_dodge()) +
labs(x='Sexo', y = 'Número de pessoas') + theme_minimal()

A inclusão das barras separadas já traz resultados interessantes, mostrando que a maior parte das mulheres sobreviveram, enquanto a taxa de sobrevivência para homens ficou abaixo de 25%. No código do gráfico, utilizamos o argumento position_dodge, que deixa as colunas de cada grupo organizadas horizontalmente, tornando a comparação entre número de sobreviventes para cada sexo rápida, pois basta comparar o nível no eixo y para cada cor. Agora, vamos separar os grupos entre classes, para verificar disparidades entre grupos diferentes de pessoas no navio:


ggplot(df_titanic, aes(x=Sex, y= Freq, fill = Survived)) + geom_bar(stat = 'identity', position = position_dodge()) +
labs(x='Sexo', y = 'Número de pessoas') + facet_wrap(~Class) + theme_minimal()

A separação indica que quase nenhuma mulher na primeira classe morreu, e quase nenhum homem da segunda classe sobreviveu. A escala de todos os gráficos é idêntica por padrão, o que pode ou não ser justificável, dependendo do tipo dos dados utilizados. No nosso caso, é importante manter tal configuração, pois permite a comparação entre classes diferentes. A última informação que podemos adicionar é a separação entre idades. Para fazer isso, vamos adicionar linhas que indicam a idade (criança ou adulto), gerando subdivisões das divisões originais. A função utilizada está disponível no pacote ggpattern.


#remotes::install_github("coolbutuseless/ggpattern")
library(ggpattern)

ggplot(df_titanic, aes(x=Sex, y= Freq, fill = Survived)) + geom_bar(stat = 'identity', position = position_dodge()) +
geom_col_pattern(
aes(Sex, Freq, pattern_fill = Age, fill = Survived),
color = 'black'
) +
labs(x='Sexo', y = 'Número de pessoas') + facet_wrap(~Class) + theme_minimal()

O resultado indica que quase todas as crianças foram salvas. Podemos variar as opções do geom_col_pattern, porém a visualização já se torna complicada pois é difícil incluir tantos detalhes em um gráfico de barras. Outra opção seria quebrar o gráfico em múltiplas categorias com o facet_wrap, porém rapidamente temos um número grande de gráficos pequenos, difíceis de comparar entre si. Uma solução que iremos apresentar aqui é a introdução de gráficos de mosaico, com o pacote vcd. A ideia de um gráfico desse tipo é utilizar os 4 lados dele como eixos, permitindo a análise de múltiplas categorias de modo conciso. Para utilizarmos a função mosaic(), os dados devem ser um array de categorias.


library(vcd)

mosaic(Titanic, shade = TRUE)


Cada retângulo do gráfico acima é facilmente identificado analisando cada um dos 4 eixos, e a comparação de tamanho entre os grupos é facilitada pois estão próximos em um mesmo gráfico. Ademais, as cores geradas são o resultado de um teste estatístico que verifica se a distribuição da amostra é independente dos atributos, sendo setores azuis estatisticamente acima do esperado, e vermelhos abaixo. O resultado indica que há muito mais tripulantes que não sobreviveram do que ocorreria se fossem salvas pessoas aleatórias, assim como muito mais mulheres foram salvas. Por outro lado, menos homens da primeira classe foram salvos do que esperado.

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