Category

Hackeando o R

Hackeando o R: introdução à estatística bayesiana

By | Hackeando o R

No Hackeando o R de hoje, vamos fazer uma introdução à estatística bayesiana, passando pelos seus conceitos básicos e realizando um exemplo no R. Esse post será o primeiro de uma série sobre o tema, abordando temas mais avançados a cada semana.

A ideia por trás da estatística bayesiana é de que, a partir de uma crença prévia sobre os parâmetros que queremos estimar, é possível melhorar nosso conhecimento sobre seus valores verdadeiros através de uma amostra. Para realizar essa melhora, utilizamos o famoso Teorema de Bayes:


Pense no caso de uma variável com distribuição binomial, como uma amostra de pessoas que responderam a uma pergunta de 'sim' ou 'não'. Tal distribuição pode ser descrita como:

Note que a distribuição da amostra é dependente do parâmetro da binomial, que é o que queremos encontrar. Para resolvermos isso, basta traduzir nosso problema para o Teorema de Bayes: a probabilidade do parâmetro ser de certo valor, a partir dos dados, é igual à probabilidade dos dados serem gerados por esse valor, multiplicada pela probabilidade do parâmetro, dividia pela probabilidade total da amostra.

A probabilidade do parâmetro é escolhida previamente, pois é a nossa crença a priori quanto ao parâmetro verdadeiro, sobrando apenas a integral para calcularmos. De modo geral, a estimação dessa integral é o grande problema para a solução analítica de aplicações desse tipo, porém no caso que estamos utilizando de exemplo, essa questão se simplifica ao utilizarmos uma crença a priori uniforme entre 0 e 1. Considerando nosso exemplo de amostras de 'sim' ou 'não', e, dado que não temos nenhum estudo prévio sobre a proporção das respostas, tal crença parece válida, pois dá peso igual a todos os valores válidos para a proporção, que é o parâmetro de uma binomial.

Utilizando o Teorema de Bayes (e alguns resultados fora do escopo do nosso post), podemos chegar em uma fórmula fechada para a probabilidade a posteriori - isto é, a nossa nova crença, após olharmos os dados - do parâmetro estimado:

Com isso, a partir de uma amostra, podemos encontrar uma distribuição estimada para o valor verdadeiro da proporção. A distribuição uniforme para a crença a priori é um tipo particular da distribuição Beta, que é a posteriori acima. Quando as crenças são da mesma família, dizemos que elas são conjugadas.

Digamos que possuímos duas amostras diferentes de tamanho 5 e 85, e ambas apresentam proporção de 'sim' de 20%. Utilizando a regra acima, podemos encontrar como fica a distribuição após computarmos a nova distribuição. Abaixo, temos um gráfico com a crença original (uniforme), e as outras duas. Como podemos ver, elas têm pico em torno da proporção da amostra, porém a distribuição com maior número de observações é muito mais estreita, indicando o maior nível de certeza que temos sobre o valor real do parâmetro.

library(ggplot2)

df = data.frame(x = seq(0.01, 1, by = 0.01), y = 1)
ggplot(df, aes(x=x, y = y)) + geom_line() +
stat_function(fun = function(x) dbeta(x, 2, 5), color = "red",
size = 1) +
stat_function(fun = function(x) dbeta(x, 20, 67), color = "darkred",
size = 1) +
xlab("Proporção") +
ylab('Densidade') +
theme_minimal(
)

Hackeando o R: diminuindo seu código com purrr

By | Hackeando o R

No Hackeando o R de hoje, vamos mostrar como podemos facilitar a repetição de funções em larga escala com o auxílio do pacote purrr. Além de possuir melhor performance do que loops básicos, as funções disponíveis nesse pacote permitem a geração de objetos em larga escala de modo conciso. Inicialmente, vamos calcular a média de cada coluna de um dataframe, comparando um loop com o purrr:

library(purrr)

dados = as.data.frame(matrix(rexp(200, rate=.1), ncol=20))

# for loop
media <- vector('double', ncol(dados))

for (i in 1:ncol(dados)) {
media[i] <- mean(dados[,i])
}

# purrr
media <- map_dbl(dados, mean)
<pre>
media
V1 V2 V3 V4 V5 V6 V7 V8 V9 V10 V11 V12 V13 
15.506445 9.066860 9.606222 7.875641 8.021313 14.424925 7.983791 8.276200 5.910139 8.545459 12.199477 9.404228 5.866048 
V14 V15 V16 V17 V18 V19 V20 
8.440641 10.624942 9.263813 5.243025 6.235834 14.552733 12.405802 

Os dois códigos geram o mesmo resultado, porém no primeiro precisamos gerar um objeto e assinalar manualmente onde deve ir cada resultado, enquanto a função map_dbl do purrr retorna diretamente o vetor de resultados. Sua interpretação é simples: para o objeto dados, calculamos a função mean() sobre cada coluna, e guardamos seu valor em um vetor do tipo double. Se quiséssemos aplicar sobre as linhas, basta transpor a matriz que gerou os dados, ou utilizar a função pmap.

Considere agora o caso de uma iteração múltipla. Digamos que você tem o dataset do Gapminder, que contém informações como o PIB per capita para diversos países separados por continente e ano, e quer gerar um gráfico para cada continente, em cada ano. No total, teríamos 60 gráficos distintos, que exigiriam a realização de um for para uma lista de continentes, e outro for dentro do primeiro com a lista de anos. Ao invés disso, podemos utilizar a map2, que itera uma função sobre as mesmas linhas de vetores distintos. O código abaixo mostra como isso fica:

library(tidyverse)
library(ggplot2)

gapminder <- read.csv("https://raw.githubusercontent.com/swcarpentry/r-novice-gapminder/gh-pages/_episodes_rmd/data/gapminder-FiveYearData.csv")

continent_year <- gapminder %>% distinct(continent, year)
continents <- continent_year %>% pull(continent) %>% as.character
years <- continent_year %>% pull(year)

plot_list <- map2(.x = continents,
.y = years,
.f = ~{
gapminder %>%
filter(continent == .x,
year == .y) %>%
ggplot() +
geom_point(aes(x = gdpPercap, y = lifeExp)) +
ggtitle(glue::glue(.x, " ", .y))
})

A partir daqui, podemos acessar um gráfico qualquer verificando sua numeração no objeto continent_year.

plot_list[15]

 

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

 

 

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

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