8-Regressao_Poisson.Rmd

---
title: "Regressão de Poisson Bayesiana"
description: "Modelos Lineares Generalizados -- Poisson"
author:
  - name: Jose Storopoli
    url: https://scholar.google.com/citations?user=xGU7H1QAAAAJ&hl=en
    affiliation: UNINOVE
    affiliation_url: https://www.uninove.br
    orcid_id: 0000-0002-0559-5176
date: August 1, 2021
citation_url: https://storopoli.github.io/Estatistica-Bayesiana/8-Regressao_Poisson.html
slug: storopoli2021regressaopoissonbayesR
bibliography: bib/bibliografia.bib
csl: bib/apa.csl
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
library(tibble)
theme_set(theme_minimal())
options(brms.backend = "rstan",
        brms.normalize = TRUE)
set.seed(123)
```

<link rel="stylesheet" href="https://cdn.rawgit.com/jpswalsh/academicons/master/css/academicons.min.css"/>

Saindo do universo dos modelos lineares, começamos a nos aventurar nos modelos linares generalizados (*generalized linear models* - GLM). O segundo deles é a regressão de Poisson.

Uma regressão de Poisson se comporta exatamente como um modelo linear: faz uma predição simplesmente computando uma soma ponderada das variáveis independentes $\mathbf{X}$ pelos coeficientes estimados $\boldsymbol{\beta}$, mais uma constante $\alpha$. Porém ao invés de retornar um valor contínuo $\boldsymbol{y}$, como a regressão linear, retorna o logarítmo natural desse valor:

$$
\log(\boldsymbol{y})= \alpha \cdot \beta_1 x_1 \cdot \beta_2 x_2 \cdot \ldots \cdot \beta_k x_k
$$

que é o mesmo que:

$$
\boldsymbol{y} = e^{(\alpha + \beta_1 x_1 + \beta_2 x_2 + \ldots + \beta_k x_k)}
$$

* $\boldsymbol{\theta} \in \mathbb{R}^p$ - parâmetros do modelo
  * $\theta_0$ - constante
  * $\theta_1, \theta_2, \dots, \theta_p$ - coeficientes das variáveis independentes $x_1, x_2, \dots, p$

* $n$ - número de variáveis independentes

A função $e^x$ é chamada de função exponencial. Veja a figura \@ref(fig:exponential-function) para uma intuição gráfica da função exponencial:

```{r exponential-function, fig.cap='Função Exponencial'}
ggplot(data = tibble(
  x = seq(0, 10, length.out = 100),
  y =  exp(x)
  ),
  aes(x, y)) +
  geom_line(size = 2, color = "steelblue") +
  ylab("Exponencial(x)")
```


Regressão de Poisson é usada quando a nossa variável dependente só pode tomar valores positivos, geralmente em contextos de dados de contagem.

## Regressão de Poisson Bayesiana

Podemos fazer uma regressão de Poisson se a variável dependente $\boldsymbol{y}$ for uma variável com dados de contagem, ou seja, $\boldsymbol{y}$ somente toma valores positivos. A função de verossimilhança de Poisson usa uma constante $\alpha$ e os coeficientes $\boldsymbol{\beta}$ porém estes são exponenciados ($e^x$):

$$
\begin{aligned}
\boldsymbol{y} &\sim \text{Poisson}\left( e^{(\alpha +  \mathbf{X} \boldsymbol{\beta})} \right) \\
\alpha &\sim \text{Normal}(\mu_\alpha, \sigma_\alpha) \\
\boldsymbol{\beta} &\sim \text{Normal}(\mu_{\boldsymbol{\beta}}, \sigma_{\boldsymbol{\beta}})
\end{aligned}
$$

Como podemos ver, o preditor linear $\alpha + \mathbf{X} \boldsymbol{\beta}$ é o logaritmo do valor de $\boldsymbol{y}$. Portanto precisamos exponenciar os valores do preditor linear:

$$
\begin{aligned}
\log(y) &= \alpha + \mathbf{X} \boldsymbol{\beta} \\
\boldsymbol{y} &= e^{ \left( \alpha + \mathbf{X} \boldsymbol{\beta} \right) } \\
\boldsymbol{y} &= e^{\alpha} \cdot e^{\left( \mathbf{X} \boldsymbol{\beta} \right)}
\end{aligned}
$$

A constante $\alpha$ e os coeficientes $\boldsymbol{\beta}$ são apenas interpretados após a exponenciação.

## Regressão de Poisson com o `rstanarm` e `brms`

O `rstanarm` e `brms` podem tolerar qualquer modelo linear generalizado e regressão de Poisson não é uma exceção. Para designar um modelo de Poisson no `rstanarm` e `brms` é preciso simplesmente alterar o argumento `family` da função `stan_glm()` ou `brm()`. Isso faz com que a família da função de verossimilhança do modelo seja modificada. Você pode controlar também a função de ligação com argumento `link` mas acredito que para a vasta maioria dos casos não é necessário alterar o padrão que Stan usa para as diferentes famílias de funções de verossimilhança.

### Exemplo Poisson - Exterminação de baratas

Para exemplo, usaremos um *dataset* chamado `roaches` [@gelmanDataAnalysisUsing2007] que está incluído no `rstanarm`. É uma base de dados com 262 observações sobre a eficácia de um sistema de controle de pragas em reduzir o número de baratas (*roaches*) em apartamentos urbanos.

Possui as seguintes variáveis:

-   `y`: variável dependente - número de baratas mortas
-   `roach1`: número de baratas antes da dedetização
-   `treatment`: *dummy* para indicar se o apartamento foi dedetizado ou não
-   `senior`: *dummy* para indicar se há apenas idosos no apartamento
-   `exposure2`: número de dias que as armadilhas de baratas foram usadas

```{r rstanarm}
# Detectar quantos cores/processadores
options(mc.cores = parallel::detectCores())
options(Ncpus = parallel::detectCores())

library(rstanarm)
data(roaches)
```

#### Descritivo das variáveis

Antes de tudo, analise **SEMPRE** os dados em mãos. Graficamente e com tabelas. Isso pode ajudar a elucidar *prioris* além de embasar melhor conhecimento de domínio do fenômeno de interesse.

Pessoalmente uso o pacote `skimr` [@skimr] com a função `skim()`:

```{r skimr}
library(skimr)

skim(roaches)
```

#### Modelo de Poisson

O modelo possuirá as seguintes variáveis como independentes: `roach1`, `treatment` e `senior`. Para todas essas variáveis, como os coeficientes estarão em um escala logarítma usarei uma *priori* de uma distribuição normal com média 0 e os respectivos desvios padrões de 0.1, 5 e 5^[depois de diversos *prior predictive checks* cheguei a essas *prioris* que não são tão diferentes das padrões do `rstanarm` para modelos de Poisson.]. O que resulta em `prior = normal(c(0, 0, 0), c(0.1, 5, 5))`. Vou colocar a *priori* da constante $\alpha$ como uma normal centrada em 0 e com um desvio padrão 2.5, basicamente traduzindo o que usamos como *priori* em modelos gaussianos (regressão linear). Isto resulta em `prior_intercept = normal(0, 2.5)`. Notem que aqui não temos erro do modelo, pois a função de verossimilhança não faz pressupostos sobre desvios como era o caso na regressão linear com uma verossimilhança Gaussiana/Normal. Portanto não há o que se especificar para `prior_aux`.

```{r model_poisson}
model_poisson <- stan_glm(
  y ~ roach1 + treatment + senior,
  data = roaches,
  family = poisson,
  prior = normal(c(0, 0, 0), c(0.1, 5, 5)),
  prior_intercept = normal(0, 2.5)
)
```

Vamos ver como ficaram as estimação dos parâmetros de interesse do modelo com `summary()`:

```{r summary-model_poisson}
summary(model_poisson)
```

Não tivemos nenhum problema nas correntes Markov pois todos os `Rhat` estão bem abaixo de `1.01`.

Vamos verificar o *Posterior Predictive Check* do modelo de Poisson na figura \@ref(fig:pp-check-poisson):

```{r pp-check-poisson, fig.cap='*Posterior Preditive Check* do modelo de Poisson'}
pp_check(model_poisson) + xlim(0, 200)
```

Este é um bom exemplo de um *Posterior Predictive Check* no qual há algo errado com o modelo. Isto indica que devemos pensar melhor nas variáveis nas *prioris*, ou até incorporar mais variáveis no modelo. A função de verossimilhança faz um bom trabalho em se moldar à densidade da variável dependente $\boldsymbol{y}$ mas ainda é necessário alguns ajustes finos no modelo.

## Interpretação dos coeficientes

Ao vermos a fórmula de regressão de Poisson percebemos a interpretação dos coeficientes requer uma transformação. A transformação que precisamos fazer é a que inverte o logarítmo de $\boldsymbol{y}$, no caso a função de exponenciação:

$$
\begin{aligned}
\log(\boldsymbol{y}) &= \alpha + \mathbf{X} \boldsymbol{\beta} \\
\boldsymbol{y} &= e^{ \left( \alpha + \mathbf{X} \boldsymbol{\beta} \right) } \\
\boldsymbol{y} &= e^{\alpha} \cdot e^{\left( \mathbf{X} \boldsymbol{\beta} \right) }
\end{aligned}
$$

Fazemos isso com a função `exp()` do `R` nos coeficientes do `model_poisson`:

```{r coefficients}
coeff <- exp(model_poisson$coefficients)
coeff
```

-   `(Intercept)`: a taxa basal de exterminação das baratas $\boldsymbol{y}$, $23$
-   `roach1`: a cada uma barata antes da exterminação há um **incremento** de 1% de baratas exterminadas
-   `treatment`: se o apartamento foi dedetizado há uma **redução** de 40% de baratas exterminada
-   `senior`: se o apartamento possui somente idosos há uma **redução** de 30% de barata exterminadas

## Inflação de Zeros^[também chamados de superdispersão.] - Modelo Binomial Negativo

Às vezes temos dados de contagem na qual há uma super-inflação de zeros na variável dependente $\boldsymbol{y}$. Nestes casos os modelos de Poisson podem não ser a melhor opção, pois ele força que a média e variância da variável dependente sejam iguais. Para isso é melhor usar função de verossimilhança binomial negativa. Lembrando que a distribuição negativa binomial tem dois parâmetros:

* Número de Sucessos ($k$)
* Probabilidade de Sucessos ($p$)

Qualquer fenômeno que pode ser modelo com uma distribuição de Poisson, pode ser modelo com uma distribuição binomial negativa [@gelman2013bayesian; @gelman2020regression].

Usando a função de verossimilhança binomial negativa nosso modelo se torna:

$$
\begin{aligned}
\boldsymbol{y} &\sim \text{Binomial Negativa} \left( e^{(\alpha + \mathbf{X} \boldsymbol{\beta})}, \phi \right) \\
\phi &\sim \text{Gamma}(0.01, 0.01) \\
\alpha &\sim \text{Normal}(\mu_\alpha, \sigma_\alpha) \\
\boldsymbol{\beta} &\sim \text{Normal}(\mu_{\boldsymbol{\beta}}, \sigma_{\boldsymbol{\beta}})
\end{aligned}
$$

Veja que, quando comparamos com o modelo de Poisson, temos uma novo parâmetro $\phi$ que parametriza a distribuição binomial negativa. Esse parâmetro é a probabilidade de sucessos $p$ da distribuição binomial negativa e geralmente usamos uma distribuição gamma como *priori* para que $\phi$ cumpra a função de um parâmetro de "dispersão recíproca".

Para usar funções de verossimilhança binomial negativa basta adicionar:

* `rstanarm`: `family = neg_binomial_2`
* `brms`: `family = negbinomial`

## Inflação de Zeros - Mistura Binomial Negativa

Mesmo usando uma função de verossimilhança binomial negativa, caso a superdispersão seja muito acentuada, o seu modelo ainda pode resultar em patologias. Uma outra sugestão é usar uma mistura com uma função de verossimilhança binomial negativa. Precisa de apenas uma única mudança na especificação do modelo:

$$
\boldsymbol{y} \begin{cases}
= 0, &\text{ se } S_i = 0 \\
\sim \text{Binomial Negativa} \left( e^{(\alpha + \mathbf{X} \boldsymbol{\beta})}, \phi \right), &\text{ se } S_i = 1
\end{cases}
$$

Aqui, $S_i$ é uma variável binária (*dummy*) indicando se a observação $i$ tem valor diferente de zero ou não. $S_i$ pode ser modelado usando uma regressão logística:

$$
\begin{aligned}
P(S_i = 1) &= \operatorname{logit}^{-1}(\mathbf{X} \boldsymbol{\gamma}) \\
\boldsymbol{\gamma} &\sim \text{Beta}(1, 1)
\end{aligned}
$$

onde $\boldsymbol{\gamma}$ é um novo conjunto de coeficientes para essa parte do modelo com *prioris* uniformes de $\text{Beta} (1, 1)$. Esse modelo de duas etapas (*two-stage*) pode ser feito apenas especificando `family = zero_inflated_negbinomial` no `brms` ou codificando diretamente no Stan. `rstanarm` não possui suporte para esse tipo de modelo baseado em uma mistura.

### Exemplo Binomial Negativo - Revisitando a exterminação de baratas

Vamos dar uma olhada melhor no histograma da variável dependente `y` (figura \@ref(fig:hist-y)). Como vocês podem ver, temos muitos zeros como observações de `y`:

```{r hist-y, fig.cap='Histograma da variável `y`'}
ggplot(roaches, aes(y)) +
  geom_histogram(bins = 15, fill = "steelblue") +
  labs(
    title = "Histograma da variável y",
    y = "Frequência"
  )
```

`roaches` é um bom candidato para um modelo binomial negativo. Vamos então estimá-lo com o `rstanarm`. Aqui vou usar as mesmas *prioris* do modelo de Poisson. `rstanarm` por padrão usa a distribuição $\text{Exponencial} (1)$ como *priori* de $\phi$^[`rstanarm` não possui suporte para especificar *prioris* com distribuição gamma.]. Caso queira mudar é só especificar a sua *priori* desejada como argumento do `prior_aux`:

```{r model_negbinomial}
model_negbinomial <- stan_glm(
  y ~ roach1 + treatment + senior,
  data = roaches,
  family = neg_binomial_2,
  prior = normal(c(0, 0, 0), c(0.1, 5, 5)),
  prior_intercept = normal(0, 2.5),
  prior_aux = rstanarm::exponential(1)
)
```

Vamos ver como ficaram as estimação dos parâmetros de interesse do modelo com `summary()`:

```{r summary-model_negbinomial}
summary(model_negbinomial)
```

Não tivemos nenhum problema nas correntes Markov pois todos os `Rhat` estão bem abaixo de `1.01`. Vejam que as inferências sobre os parâmetros estão diferentes do modelo de Poisson. Notem também que o número de amostras efetivas (*Effective Sample Size* - ESS) para todos os parâmetros do modelo negativo binomial estão bem mais elevadas que o modelo de Poisson, indicando que com a parametrização binomial negativa o amostrador MCMC de Stan conseguiu explorar muito melhor a posterior do que na parametrização de Poisson.

Vamos verificar o *Posterior Predictive Check* do modelo binomial negativo na figura \@ref(fig:pp-check-negbinomial):

```{r pp-check-negbinomial, fig.cap='*Posterior Preditive Check* do modelo Binomial Negativo'}
pp_check(model_negbinomial) + xlim(0, 200)
```

Vemos também que *Posterior Preditive Check* do modelo binomial negativo possui um aspecto muito melhor que o modelo de Poisson.

### Exemplo Mistura Binomial Negativo - $\textbf{Revisitando}^2$ a exterminação de baratas

Vamos tentar aprimorar o modelo binomial negativo incorporando a mistura sugerida acima. Para isso teremos que usar o `brms` já que o `rstanarm` não possui suporte à misturas. Vou manter todas as *prioris* dos modelos anteriores. Mas para $\phi$ usarei a distribuição $\text{Gamma}(0.01, 0.01)$ (`prior(gamma(0.01, 0.01), class = shape)`) que possui suporte no `brms` e para os coeficientes $\boldsymbol{\gamma}$ da segunda etapa do modelo usarei a *priori* uniforme $\text{Beta}(1, 1)$ (`prior(beta(1,1), class = zi)`):

```{r model_zero_inflated}
library(brms)
model_zero_inflated <- brm(
  y ~ roach1 + treatment + senior,
  data = roaches,
  family = zero_inflated_negbinomial,
  prior = c(
    prior(normal(0, 0.1), class = b, coef = roach1),
    prior(normal(0, 5), class = b, coef = treatment),
    prior(normal(0, 5), class = b, coef = senior),
    prior(normal(0, 2.5), class = Intercept),
    prior(gamma(0.01, 0.01), class = shape),
    prior(beta(1, 1), class = zi)
    )
)
```


Vamos ver como ficaram as estimação dos parâmetros de interesse do modelo com `summary()`:

```{r summary-model_zero_inflated}
summary(model_zero_inflated)
```

Não tivemos nenhum problema nas correntes Markov pois todos os `Rhat` estão bem abaixo de `1.01`. Vejam que as inferências sobre os parâmetros estão similares ao modelo negativo binomial, mas com menores erros (`Est.Error` no `brms` vs `sd` no `rstanarm`).

Vamos verificar o *Posterior Predictive Check* da mistura binomial negativa na figura \@ref(fig:pp-check-zero-inflated):

```{r pp-check-zero-inflated, fig.cap='*Posterior Preditive Check* da mistura Binomial Negativa'}
pp_check(model_zero_inflated, nsamples = 40) + xlim(0, 200)
```

As mistura negativa binomial tem quase o mesmo *Posterior Preditive Check* do modelo binomial, porém com menores erros de mensuração dos parâmetros de interesse. Isto nos satisfaz que das três alternativas aqui apresentadas a mistura negativa binomial é a melhor candidata para modelar `roaches`. O conteúdo auxiliar sobre [Comparação de Modelos](aux-Model_Comparison.html) mostra como comparar modelos Bayesianos de maneira objetiva que é um aprimoramento sobre comparação subjetivas de gráficos de *Posterior Predictive Check*.

## Atividade Prática

Veja o *dataset* `NYC_bicycle` que está disponível no diretório `datasets/`:

-  [New York City - East River Bicycle Crossings](https://www.kaggle.com/new-york-city/nyc-east-river-bicycle-crossings): `datasets/NYC_bicycle.csv`

```{r atividade}
###
```

## Ambiente

```{r SessionInfo}
sessionInfo()
```