Blog

Introdução e objetivo

Este tutorial tem por objetivo revisitar alguns conceitos de estatística, assim como a implementação de funções no R que executam os cálculos relacionados a cada um destes. Será tratado aqui: medidas de tendência central, medidas de dispersão, medidas de associação, distribuições de probabilidade e intervalos de confiança.

O pacote dplyr, tidyr e ggplot2 serão utilizados para manipulação e visualização dos dados. No entanto, não é escopo deste tutorial detalhar a documentação destes.

key words: Estatística descritiva, análise de dados, dplyr, tidyr

Medidas de tendência central

Medidas de tendência central referem-se a métricas que sumarizam regiões aonde os dados provenientes de uma variável aleatória costumam se agrupar. Não deve-se atribuir, entretanto, a esses valores centrais uma maior probabilidade de ocorrência (isso é válido para algumas distribuições de probabilidade específicas, e não uma regra geral).

As medidas de tendência central permitem sumarizar um conjuntos de dados. Veremos como calcular algumas destas. Seja uma amostra de observações \(X_i\) com \(i=1,...,N\).

1. Média Aritmética: \(\bar{X}=\frac{1}{N}\sum_{i=1}^{N}{X_i}\)

2. Média Geométrica: \(X^G=\sqrt[N]{\prod_{i=1}^{N}{X_i}}\)

3. Mediana: É o valor ao qual metade dos dados são iguais ou inferiores e a outra metade é igual ou superior.

4. Moda: É o valor cuja frequência é máxima na amostra.

Vamos calcular essas estatísticas para os dados de preço médio semanal da energia elétrica contratada para diferentes regiões do Brasil.

library(readxl)

dados <- read_excel("Dados//1_preco_eletricidade.xlsx", col_names = TRUE)

dados %>% summarise(MEDIA_PRECO = mean(PRECO),               ## calcula a média aritmética do preço
                    MEDIA_G_PRECO = exp(mean(log(PRECO))),   ## calcula a média grométrica do preço
                    MEDIANA_PRECO = median(PRECO),           ## calcula a mediana grométrica do preço
                    MODA_PRECO = DescTools::Mode(PRECO))     ## calcula a moda grométrica do preço

## # A tibble: 1 x 4
##   MEDIA_PRECO MEDIA_G_PRECO MEDIANA_PRECO MODA_PRECO
##         <dbl>         <dbl>         <dbl>      <dbl>
## 1        167.          78.5          96.6       18.6

Podemos, facilmente, calcular essas estatísticas por região e para um período específico utilizando as funções do dplyr. Suponha que queiramos calcular as médias mensais, por região, para o ano de 2021.

##Calculando a média mensal por região
dados %>% filter(ANO==2021) %>% group_by(MES, REGIAO) %>% 
  summarise(MEDIA_PRECO = mean(PRECO)) %>%
  pivot_wider(names_from = REGIAO, values_from = MEDIA_PRECO)

## # A tibble: 11 x 5
## # Groups:   MES [11]
##      MES NORDESTE NORTE SUDESTE   SUL
##    <dbl>    <dbl> <dbl>   <dbl> <dbl>
##  1     1    227.  229.    231.  228. 
##  2     2    148.  147.    162.  163. 
##  3     3     79.4  53.3   104.  104. 
##  4     4     89.7  79.2   136.  140. 
##  5     5    196.  196.    223.  229. 
##  6     6    358.  365.    366.  366. 
##  7     7    584.  584.    584.  584. 
##  8     8    584.  584.    584.  584. 
##  9     9    568.  575.    575.  575. 
## 10    10    220.  220.    220.  220. 
## 11    11     91.7  91.8    91.8  91.8

Você é desafiado para tentar reproduzir o gráfico acima utilizando o ggplot.

Medidas de dispersão e associação

Medidas de dispersão dizem respeito a forma com que, dados gerados a partir de um mesmo processo estocástico, variam. De forma simples, essas métricas iram captar o quão disperso os dados se encontram.

Veremos como calcular algumas dessas medidas. Seja uma amostra de observações \(X_i\) com \(i=1,...,N\), com média aritmética dada por \(\bar{X}\).

1. Variância: \(VAR(X)=\frac{1}{N-1}\sum_{i=1}^{N}{(X_i-\bar{X})^2 }\)

2. Desvio-Padrão: \(SD(X)=\sqrt{\frac{1}{N-1}\sum_{i=1}^{N}{(X_i-\bar{X})^2 }}\)

3. Amplitude: \(A(X)=\max{(X_i)}-\min{(X_i)}\)

4. Coeficiente de Variação: \(CV(X)=\hat{\sigma}/\bar{X}\), em que \(\hat{\sigma}\) é a estimativa do desvio-padrão da variável.

Vamos a alguns exemplos.

#Para todos os dados disponíveis
dados %>% summarise(VAR_PRECO = var(PRECO),             ## calcula a variância do preço
                    DP_PRECO = sd(PRECO),               ## calcula o desvio-padrão do preço
                    AMP_PRECO = max(PRECO)-min(PRECO),  ## calcula a amplitude do preço
                    CV_PRECO = sd(PRECO)/mean(PRECO))   ## calcula o coeficiente de variação do preço

## # A tibble: 1 x 4
##   VAR_PRECO DP_PRECO AMP_PRECO CV_PRECO
##       <dbl>    <dbl>     <dbl>    <dbl>
## 1    36759.     192.      819.     1.15

#Por região para o ano de 2021
dados %>% filter(ANO==2021) %>% group_by(REGIAO) %>% 
  summarise(VAR_PRECO = var(PRECO),            
            DP_PRECO = sd(PRECO),               
            AMP_PRECO = max(PRECO)-min(PRECO),
            CV_PRECO = sd(PRECO)/mean(PRECO))

## # A tibble: 4 x 5
##   REGIAO   VAR_PRECO DP_PRECO AMP_PRECO CV_PRECO
##   <chr>        <dbl>    <dbl>     <dbl>    <dbl>
## 1 NORDESTE    41633.     204.      534.    0.699
## 2 NORTE       43375.     208.      534.    0.718
## 3 SUDESTE     38372.     196.      513.    0.645
## 4 SUL         38177.     195.      513.    0.642

Lembre-se que a unidade de medida da variância é sempre o quadrado da unidade de medida original. Neste caso seria \((R\$/MWh)^2\).

Medidas de associação dizem respeito a forma com que, duas ou mais variáveis, variam de forma conjunta. Isto é, possuem trajetórias semelhante. Dentre as medidas mais comuns estão a covariância e correlação. Seja uma amostra de observações \(X_i\) com \(i=1,...,N\), com média aritmética dada por \(\bar{X}\) e uma amostra \(Y_i\) com \(i=1,...,N\), com média aritmética dada por \(\bar{Y}\).

1. Covariância: \(COV(X,Y)=\frac{1}{N-1}\sum_{i=1}^{N}{(X_i-\bar{X})(Y_i-\bar{Y})}\)

2. Correlação: \(\rho(X,Y)=\frac{\sum_{i=1}^{N}{(X_i-\bar{X})(Y_i-\bar{Y})}}{\sqrt{\sum_{i=1}^{N}{(X_i-\bar{X})^2 }}\sqrt{\sum_{i=1}^{N}{(Y_i-\bar{Y})^2 }}}\)

A correlação, medida muito utilizada para medição de associação, corresponde a um índice cujo intervalo de valores possíveis está restrito entre -1 e 1. Quanto mais próximo de -1, as duas variáveis possuem trajetórias conjuntas que divergem de direção, enquanto o contrário é válido para valores próximos a 1. Todavia, é importante distinguir o conceito de correlação com o de causalidade. Correlação não deve ser interpretado como causalidade entre variáveis. Causalidade é um conceito muito mais forte e complexo de ser averiguado empiricamente do que correlação. Correlação elevada entre duas variáveis não é condição suficiente e nem mesmo necessária para que haja causalidade.

Agora vamos ver como podemos calcular essas métricas utilizando o R. Utilizamos os dados da região Sudeste e Sul para tanto.

#Para todos os dados disponíveis
dados %>% filter(REGIAO %in% c("SUL", "SUDESTE")) %>% 
  pivot_wider(names_from = REGIAO, values_from = PRECO) %>% 
  summarise(COV_PRECO = cov(SUDESTE,SUL),      ## calcula a covariância do preço das regiões
            CORREL_PRECO = cor(SUDESTE,SUL))   ## calcula a correlação do preço das regiões      

## # A tibble: 1 x 2
##   COV_PRECO CORREL_PRECO
##       <dbl>        <dbl>
## 1    33260.        0.913

Como era de se esperar pela ilustração gráfica dos dados, a correlação entre as variáveis é bastante alta.

Distribuição de probabilidades

Uma distribuição de probabilidade é a função matemática que representa o processo gerador de dados. Ou seja, ela descreve o comportamento que um evento aleatório deve possuir. Para ilustrar essas definição um pouco mais técnica, imagine que, ao questionar um meteorologista a respeito das chuvas para o mês de janeiro, ele responda que elas devem permanecer entre 50 mm e 350 mm e que, em média, costuma ser de 200 mm. Todas essas afirmações são derivadas das características do evento aleatório “chuva em janeiro” e são consequência do processo gerador das chuvas que pode ser ilustrado por uma função de distribuição de probabilidades.

As distribuições de probabilidade possuem parâmetros que definem como a realização de observações daqueles dados deve ocorrer. Vamos utilizar três distribuições distintas para ilustrar isso: i) distribuição uniforme e ii) distribuição normal.

1. Distribuição Uniforme A distribuição uniforme, como o próprio nome indica, possui densidade de probabilidade uniforme para todos os valores de seu suporte. A distribuição de densidade de probabilidade neste caso é:

• \(\mathbb{P}(X)=\frac{1}{b-a}\textbf{I}_{[a,b]}(X)\). Os parâmetros da distribuição são \(a\) e \(b\).

Todos os valores entre \(a\) e \(b\) são igualmente prováveis de ocorrer nesses caso. Graficamente, uma amostra (ou população) que possui essa distribuição deve possuir frequências relativas semelhantes a figura abaixo. Neste caso, estamos assumindo \(a=10\) e \(b=100\)

2. Distribuição Normal

A distribuição normal é uma das mais populares na estatística. Suas propriedades (simetria, curtose, simplicidade de interpretação dos parâmetros, etc) a tornam extremamente conveniente, além de que graças ao teorema do limite central é possível assumir, em condições bastante gerais, que diversas variáveis irão convergir em distribuição para a normal. A distribuição de densidade de probabilidade neste caso é:

• \(\mathbb{P}(X)=\frac{1}{\sqrt{2\pi\sigma^2}}\textrm{e}^{-\frac{1}{2}\big{(}\frac{X-\mu}{\sigma}\big{)}^2}\textbf{I}_{[-\infty,\infty]}(X)\). Os parâmetros da distribuição são \(\mu\) e \(\sigma\).

Nesse caso, os valores mais próximos a \(\mu\) são os mais prováveis de se ocorrer. O valor de \(\sigma\) irá definir o quanto essas observações podem estar dispersas em relação a média (quanto maior \(\sigma\), maior o número de observações extremas). Quando trabalhamos com uma distribuição normal, os valores de \(\mu\) e \(\sigma\) podem ser estimados pela média aritmética de amostras \(\bar{X}\) e pelo estimador da variância \(\hat{\sigma}\) apresentado nas seções anteriores. Conforme o tamanho dessa amostra aumenta a estimativa dada por \(\bar{X}\) e \(\hat{\sigma}\) convergem para o verdadeiro valor dos parâmetros dessa distribuição. Graficamente, uma amostra (ou população) que possui essa distribuição deve possuir frequências relativas semelhantes a figura abaixo. Neste caso, estamos assumindo \(\mu=0\) e \(\sigma=1\).

Modelos de previsão e distribuição de probabilidades

O conceito de distribuição de probabilidade é importante pois diversos modelos possuem algum pressuposto quanto a distribuição dos dados. Portanto, a validade dessas hipóteses quando utilizamos essas técnicas é fundamental para que sejam considerados os procedimentos mais adequados de análise. Por exemplo, imagine que uma variável de interesse \(Y_t\), em que \(t\) indica a indexação no tempo desses valores, tenha distribuição normal, de modo que \(Y_t \sim \mathcal{N}(\mu,\sigma^{2})\) para todo \(t\). Se isso for verdade, a melhor previsão que poderíamos fazer para \(Y_t\) seria a média \(\mu\).

Uma forma conveniente de representar este modelo, cuja premissa é que \(Y_t \sim \mathcal{N}(\mu,\sigma^{2})\), é utilizando a equação:

\[ Y_t = \mu + \epsilon_t | \epsilon_t \sim \mathcal{N}(0,\sigma^{2})\]

Note que podemos interpretar que a variável \(Y_t\) é a soma de um componente determinístico (isto é, não aleatório) \(\mu\) e de um termo aleatório \(\epsilon_t\).

No entanto, na prática, só possuímos um conjunto de observações. Nesse caso, como descobrir qual a distribuição dos dados observados? Embora hajam testes e procedimentos mais formais, um bom primeiro passo é representar o histograma com a frequência relativa dos dados. Essa frequência relativa é equiparável com as densidades de probabilidade mencionadas anteriormente. Vamos utilizar o R para obter um histograma para os dados de preços que vinhamos trabalhando até então.

#Para todos os dados disponíveis
dados %>% filter(REGIAO =="SUDESTE" & ANO>=2015) %>% 
  ggplot(aes(x=PRECO, y=stat(count)/sum(count))) + 
  geom_histogram(colour="black", fill="gray") +
  labs(title = "Histograma do Preço", x="Resíduos", y="Frequência Relativa")

A princípio o histograma não possui uma aparência que nos leve a acreditar que possamos representar esse processo com uma distribuição normal, por exemplo. Porém, e se alterarmos o modelo utilizado para:

\[ Y_t = \mu + \phi Y_{t-1} + \epsilon_t | \epsilon_t \sim \mathcal{N}(0,\sigma^{2})\] Agora \(Y_t\) é a soma de um termo determinístico e dois termos aleatórios (\(Y_{t-1}\) e \(\epsilon_t\)). Em algumas lições mais adiante mostraremos algumas técnicas para identificar a melhor especificação de modelos, assim como estimá-los. Por hora, estamos interessados apenas em ilustrar como as distribuições de probabilidade estão relacionadas a construção de modelos e como estas auxiliam na verificação do quão adequado é o procedimento utilizado.

Neste exemplo, caso o modelo proposto esteja correto, é esperado que nossas estimativas sobre \(\mu\) e \(\phi\) façam com que \(\epsilon_t\) possua uma distribuição que se assemelhe a uma normal. Vamos checar essa afirmação.

Utilizando o pacote fabletools podemos estimar esse modelo. Abaixo você visualizará o código para estimação e o histograma dos resíduos, isto é, os valores estimados de \(\epsilon_t\).

library(fabletools)
library(fable)
library(tsibble)

#Estimando o modelo para a região SUDESTE com dados a partir de 2015
dados %>% filter(REGIAO =="SUDESTE" & ANO>=2015) %>% select(SEMANA,REGIAO,PRECO) %>%
  as_tsibble(index = SEMANA, key=REGIAO) %>%
  model(ar=ARIMA(log(PRECO))) %>%
  augment() %>% ggplot(aes(x = .resid, y=stat(count)/sum(count))) +
  geom_histogram(colour="black", fill="gray") +
  labs(title = "Histograma dos resíduos", x="Resíduos", y="Frequência Relativa")

O histograma se assemelha ao de uma distribuição normal com média zero. Vale ressaltar que a visualização gráfica é um procedimento com baixo rigor para validação do modelo. No entanto, o objetivo aqui é ilustrar, de uma forma intuitiva, como as hipóteses acerca de distribuição de probabilidades se refletem nos modelos estimados.

Intervalos de confiança

O intervalo de confiança, ou estimativa por intervalo, tem por objetivo incorporar as incertezas associadas a uma estatística calculada. Por exemplo, em pesquisas de intenção de voto é comum que anunciem a chamada “margem de erro de 95%” que gira em torno de 2%. Portanto, se um candidato tem uma intenção de, digamos, 15%, então com 95% de confiança a verdadeira intensão de voto neste candidato está no intervalo de 13% e 17%. Esse é chamado de intervalo de confiança (IC).

Do ponto de vista prático o IC permite criar cenários baseados na distribuição de probabilidade dos dados e estimativas. Vamos ilustrar a utilização deste conceito realizando a previsão do preço da energia elétrica para a região sudeste, utilizando o modelo estimado no tópico anterior.

dados %>% filter(REGIAO =="SUDESTE" & ANO>=2015) %>% select(SEMANA,REGIAO,PRECO) %>%
  as_tsibble(index = SEMANA, key=REGIAO) %>%
  model(ar=ARIMA(log(PRECO))) %>%
  forecast(h=3) %>%
  autoplot(dados %>% filter(REGIAO =="SUDESTE" & ANO>=2021) %>% as_tsibble(index = SEMANA, key=REGIAO)) +
  labs(title="Histórico e previsão de preços - Sudeste",
        y="R$/MWh", x="Semana")

No gráfico ilustramos a previsão para às próximas 3 semanas. A linha azul sólida representa a previsão pontual, enquanto as áreas em azul correspondem ao IC com 80% e 95%, respectivamente. Por exemplo, a previsão fornecida para o modelo para a próxima semana dos dados foi de R$102/MWh. Porém, com 95% de confiança, o valor verdadeiro do preço deve estar no intervalo entre R$59/MWh e R$165/MWh.

Além disso, perceba que conforme o horizonte de previsão aumenta o IC fica mais amplo. Intuitivamente (pois existe uma justificativa formal para isso), esse fato é consequência da maior incerteza associada a previsões mais distantes.

Referências sugeridas

HOFFMANN, R. Estatística para Economistas. 4ª edição. São Paulo, Cengage Learning, 2013.