aux-Model_Comparison.Rmd

---
title: "Comparação de Modelos"
description: "Como comparar modelos Bayesianos usando métricas objetivas"
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/aux-Model_Comparison.html
slug: storopoli2021bayescompararmodelosR
bibliography: bib/bibliografia.bib
csl: bib/apa.csl
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
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"/>

Depois de estimarmos um modelo Bayesiano, muitas vezes queremos medir sua precisão preditiva, por si só ou para fins de comparação, seleção ou cálculo de média do modelo [@geisser1979predictive].

Nas aulas desta disciplina nos debruçamos sobre diferentes gráficos de *Posterior Predictive Check* de diferentes modelos. Esta é uma maneira subjetiva e arbitrária de analisarmos e compararmos modelos entre si usando sua precisão preditiva. Há uma maneira objetiva de compararmos modelos Bayesianos com uma métrica robusta que nos ajude a selecionar qual o melhor modelo dentre o rol de modelos candidatos. Ter uma maneira objetiva de comparar modelos e escolher o melhor dentre eles é muito importante pois no *workflow* Bayesiano geralmente temos diversas iterações entre *prioris* e funções de verossimilhança o que ocasiona na criação de diversos modelos diferentes [@gelmanBayesianWorkflow2020].

## Técnicas de Comparação de Modelos

Temos diversas técnicas de comparação de modelos que usam a precisão preditiva, sendo as principais:

* *Leave-one-out cross-validation* (LOO) [@vehtariPracticalBayesianModel2015]
* *Deviance Information Criterion* (DIC) [@spiegelhalter2002bayesian],  mas sabe-se que tem alguns problemas, que surgem em parte por não ser totalmente Bayesiano, pois se baseia em uma estimativa pontual [@van2005dic]
* *Widely Applicable Information Criteria* (WAIC) [@watanabe2010asymptotic], totalmente Bayesiano no sentido de que usa toda a distribuição posterior, e é assintoticamente igual ao LOO [@vehtariPracticalBayesianModel2015]

Destes, já descartamos o DIC por termos problemas e ser baseado em uma estimativa pontual (afinal somos Bayesianos, se estivéssemos interessados em estimativas pontuais estaríamos ainda maximizando funções de verossimilhança e achando a moda como os frequentistas fazem).

LOO é computacionalmente intensivo, afinal na validação cruzada (*cross-validation*) re-estimamos o modelo para cada partição dos dados. *Leave-one-out* quer dizer que para um *dataset* de tamanho $N$ estimaremos $N-1$ modelos para $N-1$ partições do *dataset*. Ou seja, deixamos uma observação para fora e estimamos o modelo usando a partição de dados sem essa observação e repetimos para todas as observações do *dataset*. Para superar essa dificuldade, LOO pode ser aproximado por amostragem de importância [@gelfand1996model]. Mas tal estimativa pode ser imprecisa. Usando uma amostragem de importância com suavização de Pareto (*Pareto smoothed importance sampling* -- PSIS) podemos aproximar uma estimativa confiável ajustando uma distribuição de Pareto à cauda superior da distribuição dos pesos de importância. PSIS nos permite calcular LOO usando pesos de importância que de outra forma seriam instáveis. Tal aproximação é cunhada de PSIS-LOO [@vehtariPracticalBayesianModel2015].

WAIC pode ser visto como uma melhoria do DIC para modelos bayesianos. Apesar de WAIC ser assintoticamente igual ao LOO, PSIS-LOO é mais robusto no quando usamos *prioris* não-informativas ou na presença observações influentes (*outliers*).

## Como mensuramos precisão preditiva?

Bayesianos mensuram precisão preditiva usando simulações da distribuição posterior $\tilde{y}$ do modelo. Para isso temos a distribuição preditiva posterior (*predictive posterior distribution*):

$$
p(\tilde{y} \mid y) = \int p(\tilde{y}_i \mid \theta) p(\theta \mid y) d \theta
$$

Onde $p(\theta \mid y)$ é a distribuição posterior do modelo (aquela que o `rstanarm` e `brms` estima para nós). A fórmula acima significa que calculamos a integral de toda a probabilidade conjunta da distribuição posterior preditiva com a distribuição posterior do nosso modelo: $p(\tilde{y}_i \mid \theta) p(\theta \mid y)$. Quanto maior a distribuição preditiva posterior $p(\tilde{y} \mid y)$ melhor será a precisão preditiva do modelo. Para mantermos comparabilidade com um dado *dataset*, calculamos a esperança dessa medida (do inglês *expectation* que pode ser também interpretada como a média ponderada) para cada uma das $N$ observações do *dataset*:

$$
\operatorname{elpd} = \sum_{i=1}^N \int p_t(\tilde{y}_i) \log p(\tilde{y}_i \mid y) d \tilde{y}
$$

onde $\operatorname{elpd}$ é esperança do log da densidade preditiva pontual (*expected log pointwise predictive density*) e $p_t(\tilde{y}_i)$ é a distribuição representando o verdadeiro processo generativo dos dados para $\tilde{y}_i$. Os $p_t(\tilde{y}_i)$ são desconhecidos e geralmente usamos validação cruzada ou WAIC para aproximar a estimação da $\operatorname{elpd}$.

### *Leave-one-out Cross-Validation* (LOO)

Podemos calcular a $\operatorname{elpd}$ usando LOO:

$$
\operatorname{elpd}_{\text{loo}} = \sum_{i=1}^N \log p(y_i \mid y_{-i})
$$

onde

$$
p(y_i \mid y_{-i}) = \int p(y_i \mid \theta) p(\theta \mid y_{-i}) d \theta
$$

que é a densidade preditiva com uma observação a menos condicionada nos dados sem a observação $i$ ($y_{-i}$). Quase sempre usamos a aproximação PSIS-LOO pela sua robustez e baixo custo computacional.

### *Widely Applicable Information Criteria*^[também chamado às vezes de *Watanabe-Akaike Information Criteria*.] (WAIC)

WAIC [@watanabe2010asymptotic], assim como o LOO também é uma abordagem alternativa para calcularmos a $\operatorname{elpd}$ e é definida como:

$$
\widehat{\operatorname{elpd}}_{\text{waic}} = \widehat{\operatorname{lpd}} - \widehat{p}_{\text{waic}}
$$

onde $\widehat{\operatorname{lpd}}$ é a estimativa do log da densidade preditiva pontual (*log pointwise predictive density*):

$$
\widehat{\operatorname{lpd}} = \sum_{i=1}^N \log p(y_i \mid y) = \sum_{i=1}^N \log \int p(y_i \mid \theta) p(\theta \mid y) d \theta
$$

e $\widehat{p}_{\text{waic}}$ é o número estimado efetivo de paramêtros e calculado com base em:

$$
\widehat{p}_{\text{waic}} = \sum_{i=1}^N \operatorname{var}_{\text{post}} (\log p(y_i \mid \theta))
$$

que conseguimos calcular usando a variância posterior do log da densidade preditiva para cada observação $y_i$:

$$
\widehat{p}_{\text{waic}} = \sum_{i=1}^N V^S_{s=1} (\log p(y_i \mid \theta^s))
$$

onde $V^S_{s=1}$ representa a variância da amostra:

$$
V^S_{s=1} a_s = \frac{1}{S-1} \sum^S_{s=1} (a_s - \bar{a})^2
$$

### *$K$-fold Cross-Validation*

Da mesma maneira que conseguimos cacular a $\operatorname{elpd}$ usando LOO com $N-1$ partições do *dataset* podemos também calcular com qualquer número de partições que quisermos. Tal abordagem é chamada de validação cruzada usando $K$ partições (*$K$-fold Cross-Validation*, encurtado para *$K$-fold CV*). Ao contrário de LOO, não conseguimos aproximar a $\operatorname{elpd}$ usando *$K$-fold CV* e precisamos fazer a computação atual da $\operatorname{elpd}$ sobre $K$ partições que quase sempre envolve um alto custo computacional.

## Implementações em Stan, `rstanarm` e `brms`

Stan e suas interfaces (`rstanarm` e `brms`) conseguem calcular ou aproximar a $\operatorname{elpd}$ com o pacote [`loo`](https://mc-stan.org/loo) [@loo]. Conseguimos aproximar com LOO e calcular com WAIC e *$K$-fold CV*. As respectivas funções são:

* `loo()`
* `waic()`
* `kfold()` -- o padrão são 10 partições (`K = 10`)

### Revisitando os modelos de Poisson

Na [Aula 8 - Regressão de Poisson](8-Regressao_Poisson.html) tínhamos 3 modelos que estimamos usando o *dataset* `roaches` [@gelmanDataAnalysisUsing2007] incluído no `rstanarm`:

* Regressão de Poisson usando o `rstanarm`
* Regressão binomial negativa usando o `rstanarm`
* Mistura binomial usando o `brms`

Nesses modelos fizemos uma verificação subjetiva e arbitrária com o gráfico *Posterior Predictive Check* que compara o histograma da variável dependente $y$ contra o histograma de variáveis dependentes simuladas pelo modelo $y_{\text{rep}}$ após a estimação dos parâmetros. A ideia é que os histogramas reais e simulados se misturem e não haja divergências. Fizemos isso com a função `pp_check()`.

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

library(rstanarm)
data(roaches)
```

O modelo de Poisson no `rstanarm`:

```{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),
  seed = 123
)
```

O modelo negativo binomial no `rstanarm`:

```{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),
  seed = 123
)
```

A mistura negativo binomial no `brms`:

```{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)
    ),
  seed = 123
)
```

Então aproximamos a ${\operatorname{elpd}}$ de todos os modelos usando PSIS-LOO com o pacote [`loo`](https://mc-stan.org/loo) [@loo] e a função `loo()`, que possui o mesmo nome do pacote:

```{r loo}
library(loo)

loo_poisson <- loo::loo(model_poisson)
loo_negbinomial <- loo::loo(model_negbinomial)
loo_zero_inflated <- loo::loo(model_zero_inflated)
```

Conseguimos comparar os modelos agora de maneira objetiva com o `loo_compare()` inserindo como argumento os outputs da função `loo()` de cada modelo desejado:

```{r loo_compare}
loo::loo_compare(loo_poisson, loo_negbinomial, loo_zero_inflated)
```

`elpd_diff` é a diferença da $\operatorname{elpd}_{\text{loo}}$ para os modelos relativo sempre ao modelo que tem o maior valor de $\operatorname{elpd}_{\text{loo}}$. Ou seja, o modelo com maior $\operatorname{elpd}_{\text{loo}}$ recebe sempre o valor de 0 na comparação.

`se_diff` é o erro padrão das diferenças dos modelos comparados ao melhor modelo. A fórmula é:

$$
\operatorname{se}(\widehat{\operatorname{elpd}}_{\text{loo}}^A - \widehat{\operatorname{elpd}}_{\text{loo}}^B) = \sqrt{n V^n_{i=1}(\widehat{\operatorname{elpd}}_{\text{loo}, i}^A - \widehat{\operatorname{elpd}}_{\text{loo}, i}^B)}
$$

onde $\widehat{\operatorname{elpd}}_{\text{loo}}^A$ é a $\operatorname{elpd}$ estimada por PSIS-LOO do modelo $A$.

No nosso caso, podemos sem dúvida afirmar que o modelo de Poisson é bem inferior na precisão preditiva. Afinal sua  $\operatorname{elpd}^{\text{loo}}$ é `-5356.3` o que é no mínimo 7 vezes maior que o erro padrão `720.7` quando comparado com o modelo com maior $\operatorname{elpd}^{\text{loo}}$: o modelo negativo binomial `model_negbinomial`.

Já a diferença entre o modelo negativo binomial e a mistura negativa binomial a diferença é negligente pois além de pequena (`0.7`) está à aproximadamente 1 erro padrão `0.5` de diferença. Caso queira, podemos escolher o modelo negativo binomial devido à sua menor complexidade mensurado pelo número de parâmetros do modelo (menos parâmetros, menos complexo, menos pressupostos etc.).

## Comentários Finais

Comparação de modelos Bayesianos é um campo de pesquisa ainda em movimento. LOO, em especial PSIS-LOO, reflete o estado da arte sobre como comparar modelos Bayesianos de maneira robusta e com baixo custo computacional. Caso o leitor se interesse, recomendo se aprofundar no pacote pacote [`loo`](https://mc-stan.org/loo) [@loo] e também no artigo que descreve a sua fórmulação matemática e estatística [@vehtariPracticalBayesianModel2015]. Por fim, sugiro que vejam as publicações e os materiais disponíveis do [Aki Vehtari](https://avehtari.github.io/modelselection/), um dos desenvolvedores Stan e a principal referência em comparação de modelos Bayesianos.

## Ambiente

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