Comparação de Modelos

Como comparar modelos Bayesianos usando métricas objetivas

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 (Geisser & Eddy, 1979).

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 (Gelman et al., 2020).

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) (Vehtari, Gelman, & Gabry, 2015)
• Deviance Information Criterion (DIC) (Spiegelhalter, Best, Carlin, & Van Der Linde, 2002), mas sabe-se que tem alguns problemas, que surgem em parte por não ser totalmente Bayesiano, pois se baseia em uma estimativa pontual (Van Der Linde, 2005)
• Widely Applicable Information Criteria (WAIC) (Watanabe & Opper, 2010), totalmente Bayesiano no sentido de que usa toda a distribuição posterior, e é assintoticamente igual ao LOO (Vehtari et al., 2015)

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 (Gelfand, 1996). 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 (Vehtari et al., 2015).

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

WAIC (Watanabe & Opper, 2010), 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 (Vehtari et al., 2020). 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 tínhamos 3 modelos que estimamos usando o dataset roaches (Gelman & Hill, 2007) 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().

# Detectar quantos cores/processadores
options(mc.cores = parallel::detectCores())
options(Ncpus = parallel::detectCores())

library(rstanarm)
data(roaches)

O modelo de Poisson no rstanarm:

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:

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:

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 (Vehtari et al., 2020) e a função loo(), que possui o mesmo nome do pacote:

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:

loo::loo_compare(loo_poisson, loo_negbinomial, loo_zero_inflated)

                    elpd_diff se_diff
model_negbinomial       0.0       0.0
model_zero_inflated    -0.7       0.5
model_poisson       -5356.3     720.0

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 (Vehtari et al., 2020) e também no artigo que descreve a sua fórmulação matemática e estatística (Vehtari et al., 2015). Por fim, sugiro que vejam as publicações e os materiais disponíveis do Aki Vehtari, um dos desenvolvedores Stan e a principal referência em comparação de modelos Bayesianos.

Ambiente

sessionInfo()

R version 4.1.0 (2021-05-18)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.1/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods  
[7] base     

other attached packages:
 [1] loo_2.4.1            DescTools_0.99.41    gapminder_0.3.0     
 [4] ggExtra_0.9          dplyr_1.0.6          rstan_2.21.2        
 [7] StanHeaders_2.21.0-7 MASS_7.3-54          ggforce_0.3.3       
[10] gganimate_1.0.7      plotly_4.9.3         carData_3.0-4       
[13] DiagrammeR_1.0.6.1   brms_2.15.0          rstanarm_2.21.1     
[16] Rcpp_1.0.6           skimr_2.1.3          readr_1.4.0         
[19] readxl_1.3.1         tibble_3.1.2         ggplot2_3.3.3       
[22] patchwork_1.1.1      cowplot_1.1.1       

loaded via a namespace (and not attached):
  [1] utf8_1.2.1           tidyselect_1.1.1     lme4_1.1-27         
  [4] htmlwidgets_1.5.3    grid_4.1.0           munsell_0.5.0       
  [7] codetools_0.2-18     ragg_1.1.2           distill_1.2         
 [10] DT_0.18              gifski_1.4.3-1       miniUI_0.1.1.1      
 [13] withr_2.4.2          Brobdingnag_1.2-6    colorspace_2.0-1    
 [16] highr_0.9            knitr_1.33           rstudioapi_0.13     
 [19] stats4_4.1.0         bayesplot_1.8.0      labeling_0.4.2      
 [22] repr_1.1.3           mnormt_2.0.2         polyclip_1.10-0     
 [25] farver_2.1.0         bridgesampling_1.1-2 rprojroot_2.0.2     
 [28] coda_0.19-4          vctrs_0.3.8          generics_0.1.0      
 [31] xfun_0.23            R6_2.5.0             markdown_1.1        
 [34] isoband_0.2.4        gamm4_0.2-6          projpred_2.0.2      
 [37] assertthat_0.2.1     promises_1.2.0.1     scales_1.1.1        
 [40] rootSolve_1.8.2.1    gtable_0.3.0         downlit_0.2.1       
 [43] processx_3.5.2       lmom_2.8             rlang_0.4.11        
 [46] systemfonts_1.0.2    splines_4.1.0        lazyeval_0.2.2      
 [49] checkmate_2.0.0      inline_0.3.19        yaml_2.2.1          
 [52] reshape2_1.4.4       abind_1.4-5          threejs_0.3.3       
 [55] crosstalk_1.1.1      backports_1.2.1      httpuv_1.6.1        
 [58] rsconnect_0.8.18     tools_4.1.0          bookdown_0.22       
 [61] ellipsis_0.3.2       jquerylib_0.1.4      RColorBrewer_1.1-2  
 [64] proxy_0.4-25         ggridges_0.5.3       plyr_1.8.6          
 [67] base64enc_0.1-3      visNetwork_2.0.9     progress_1.2.2      
 [70] purrr_0.3.4          ps_1.6.0             prettyunits_1.1.1   
 [73] zoo_1.8-9            here_1.0.1           magrittr_2.0.1      
 [76] data.table_1.14.0    magick_2.7.2         colourpicker_1.1.0  
 [79] tmvnsim_1.0-2        mvtnorm_1.1-1        matrixStats_0.59.0  
 [82] hms_1.1.0            shinyjs_2.0.0        mime_0.10           
 [85] evaluate_0.14        xtable_1.8-4         shinystan_2.5.0     
 [88] gridExtra_2.3        rstantools_2.1.1     compiler_4.1.0      
 [91] V8_3.4.2             crayon_1.4.1         minqa_1.2.4         
 [94] htmltools_0.5.1.1    mgcv_1.8-35          later_1.2.0         
 [97] tidyr_1.1.3          expm_0.999-6         Exact_2.1           
[100] RcppParallel_5.1.4   lubridate_1.7.10     DBI_1.1.1           
[103] tweenr_1.0.2         boot_1.3-28          Matrix_1.3-3        
[106] cli_2.5.0            parallel_4.1.0       igraph_1.2.6        
[109] pkgconfig_2.0.3      xml2_1.3.2           dygraphs_1.1.1.6    
[112] bslib_0.2.5.1        stringr_1.4.0        callr_3.7.0         
[115] digest_0.6.27        rmarkdown_2.8        cellranger_1.1.0    
[118] gld_2.6.2            curl_4.3.1           shiny_1.6.0         
[121] gtools_3.8.2         nloptr_1.2.2.2       lifecycle_1.0.0     
[124] nlme_3.1-152         jsonlite_1.7.2       viridisLite_0.4.0   
[127] fansi_0.5.0          pillar_1.6.1         lattice_0.20-44     
[130] fastmap_1.1.0        httr_1.4.2           pkgbuild_1.2.0      
[133] survival_3.2-11      glue_1.4.2           xts_0.12.1          
[136] png_0.1-7            shinythemes_1.2.0    class_7.3-19        
[139] stringi_1.6.2        sass_0.4.0           textshaping_0.3.4   
[142] e1071_1.7-7         

---

1. também chamado às vezes de Watanabe-Akaike Information Criteria.↩︎