Regressão Logística Bayesiana

Modelos Lineares Generalizados – Binomial

Jose Storopoli https://scholar.google.com/citations?user=xGU7H1QAAAAJ&hl=en (UNINOVE)https://www.uninove.br
August 1, 2021

Saindo do universo dos modelos lineares, começamos a nos aventurar nos modelos linares generalizados (generalized linear models – GLM). O primeiro deles é a regressão logística (também chamada de regressão binomial).

Uma regressão logística 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 a função logística desse valor:

\[ \text{Logística}(x) = \frac{1}{1 + e^{(-x)}} \]

Usamos regressão logística quando a nossa variável dependente é binária. Ela possui apenas dois valores distintos, geralmente codificados como \(0\) ou \(1\). Veja a figura 1 para uma intuição gráfica da função logística:

ggplot(data = tibble(
  x = seq(-10, 10, length.out = 100),
  y =  1 / (1 + exp(-x))
  ),
  aes(x, y)) +
  geom_line(size = 2, color = "steelblue") +
  ylab("Logística(x)")
Função Logística

Figure 1: Função Logística

Como podemos ver, a função logística é basicamente um mapeamento de qualquer número real para um número real no intervalo entre 0 e 1:

\[ \text{Logística}(x) = \{ \mathbb{R} \in [- \infty , + \infty] \} \to \{ \mathbb{R} \in (0, 1) \} \]

Ou seja, a função logística é a candidata ideal para quando precisamos converter algo contínuo sem restrições para algo contínuo restrito entre 0 e 1. Por isso ela é usada quando precisamos que um modelo tenha como variável dependente uma probabilidade (lembrando que qualquer numero real entre 0 e 1 é uma probabilidade válida). No caso de uma variável dependente binária, podemos usar essa probabilidade como a chance da variável dependente tomar valor de 0 ou de 1.

Comparativo com a Regressão Linear

A regressão linear segue a seguinte formulação matemática:

\[ \text{Linear} = \alpha + \beta_1 x_1 + \beta_2 x_2 + \dots + \beta_k x_k \]

Se você implementar uma pequena gambiarra matemática, você terá a regressão logística:

Exemplo:

\[\text{Previsão de Morte} = \text{Logística} \big(-10 + 10 \cdot \text{cancer} + 12 \cdot \text{diabetes} + 8 \cdot \text{obesidade} \big)\]

Regressão Logística Bayesiana

Podemos modelar regressão logística de duas maneiras. A primeira opção com uma função verossimilhança Bernoulli e a segunda opção com uma função verossimilhança binomial.

Com a verossimilhança Bernoulli modelamos uma variável dependente binária \(\boldsymbol{y}\) que é o resultado de um experimento de Bernoulli com uma certa probabilidade \(p\).

Já na verossimilhança binomial modelamos uma variável dependente contínua \(\boldsymbol{y}\) que é o número de sucessos de \(n\) experimentos Bernoulli independentes.

Usando Verossimilhança Bernoulli

\[ \begin{aligned} \boldsymbol{y} &\sim \text{Bernoulli}\left( p\right) \\ p &\sim \text{Logística}(\alpha + \mathbf{X} \boldsymbol{\beta}) \\ \alpha &\sim \text{Normal}(\mu_\alpha, \sigma_\alpha) \\ \boldsymbol{\beta} &\sim \text{Normal}(\mu_{\boldsymbol{\beta}}, \sigma_{\boldsymbol{\beta}}) \end{aligned} \]

Sendo que:

Usando Verossimilhança Binomial

\[ \begin{aligned} \boldsymbol{y} &\sim \text{Binomial}\left(n, p\right) \\ p &\sim \text{Logística}(\alpha + \mathbf{X} \boldsymbol{\beta}) \\ \alpha &\sim \text{Normal}(\mu_\alpha, \sigma_\alpha) \\ \boldsymbol{\beta} &\sim \text{Normal}(\mu_{\boldsymbol{\beta}}, \sigma_{\boldsymbol{\beta}}) \end{aligned} \]

Sendo que:

Em ambas famílias de verossimilhança, o que falta é especificar quais são as prioris dos parâmetros do modelo:

Stan instanciará um modelo bayesiano com os dados fornecidos (\(\boldsymbol{y}\) e \(\mathbf{X}\)) e encontrará a distribuição posterior dos parâmetros de interesse do modelo (\(\alpha\) e \(\boldsymbol{\beta}\)) calculando a distribuição posterior completa de:

\[ P(\boldsymbol{\theta} \mid \boldsymbol{y}) = P(\alpha, \boldsymbol{\beta} \mid \boldsymbol{y}) \]

Regressão Logística com o rstanarm e brms

O rstanarm e brms podem tolerar qualquer modelo linear generalizado e regressão logística não é uma exceção. Não há função de verossimilhança Bernoulli em Stan (e, consequentemente, em nenhuma de suas interfaces) pois ela é apenas um caso especial da função de verossimilhança binomial com apenas um único experimento Bernoulli. Ou seja:

\[ \text{Bernouli}(p) = \text{Binomial}(1, p) \]

Para designar um modelo binomial 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 - Propensão a mudar de poço de água contaminado

Para exemplo, usaremos um dataset chamado wells (Gelman & Hill, 2007) que está incluído no rstanarm. É uma survey com 3.200 residentes de uma pequena área de Bangladesh na qual os lençóis freáticos estão contaminados por arsênico. Respondentes com altos níveis de arsênico nos seus poços foram encorajados para trocar a sua fonte de água para uma níveis seguros de arsênico.

Possui as seguintes variáveis:

options(mc.cores = parallel::detectCores())
options(Ncpus = parallel::detectCores())

library(rstanarm)
data(wells)

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 (Waring et al., 2021) com a função skim():

library(skimr)

skim(wells)
Table 1: Data summary
Name wells
Number of rows 3020
Number of columns 5
_______________________
Column type frequency:
numeric 5
________________________
Group variables None

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
switch 0 1 0.58 0.49 0.00 0.00 1.0 1.0 1.0 ▆▁▁▁▇
arsenic 0 1 1.66 1.11 0.51 0.82 1.3 2.2 9.7 ▇▂▁▁▁
dist 0 1 48.33 38.48 0.39 21.12 36.8 64.0 339.5 ▇▂▁▁▁
assoc 0 1 0.42 0.49 0.00 0.00 0.0 1.0 1.0 ▇▁▁▁▆
educ 0 1 4.83 4.02 0.00 0.00 5.0 8.0 17.0 ▇▇▅▁▁

Modelo Binomial

O modelo possuirá as seguintes variáveis como independentes: dist,arsenic, assoc e educ. Para todas essas variáveis, como os coeficientes estarão em uma escala logit (que é o inverso da transformação logística), usarei uma priori de uma distribuição normal com média 0 e desvio padrão 2.5. O que resulta em prior = normal(0, 2.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.

model_binomial <- stan_glm(
  switch ~ dist + arsenic + assoc + educ,
  data = wells,
  family = binomial,
  prior = normal(0, 2.5),
  prior_intercept = normal(0, 2.5)
)

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

summary(model_binomial)

Model Info:
 function:     stan_glm
 family:       binomial [logit]
 formula:      switch ~ dist + arsenic + assoc + educ
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 3020
 predictors:   5

Estimates:
              mean   sd   10%   50%   90%
(Intercept) -0.2    0.1 -0.3  -0.2   0.0 
dist         0.0    0.0  0.0   0.0   0.0 
arsenic      0.5    0.0  0.4   0.5   0.5 
assoc       -0.1    0.1 -0.2  -0.1   0.0 
educ         0.0    0.0  0.0   0.0   0.1 

Fit Diagnostics:
           mean   sd   10%   50%   90%
mean_PPD 0.6    0.0  0.6   0.6   0.6  

The mean_ppd is the sample average posterior predictive distribution of the outcome variable (for details see help('summary.stanreg')).

MCMC diagnostics
              mcse Rhat n_eff
(Intercept)   0.0  1.0  3655 
dist          0.0  1.0  3518 
arsenic       0.0  1.0  2853 
assoc         0.0  1.0  3460 
educ          0.0  1.0  3419 
mean_PPD      0.0  1.0  3312 
log-posterior 0.0  1.0  1781 

For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).

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 binomial na figura 2:

pp_check(model_binomial)
*Posterior Preditive Check* do modelo binomial

Figure 2: Posterior Preditive Check do modelo binomial

Este é um ótimo exemplo de um Posterior Predictive Check no qual não conseguimos identificar qualquer anomalia no modelo.

Interpretação dos coeficientes

Ao vermos a fórmula de regressão logística percebemos a interpretação dos coeficientes requer uma transformação. A transformação que precisamos fazer á a que inverte a função logística. Mas antes preciso falar sobre qual a diferença matemática entre probabilidade (em inglês probability) e chances (em inglês odds). Probabilidade já discutimos na Aula 0 - O que é Estatística Bayesiana?: um número real entre 0 e 1 que representa a certeza de que um evento irá acontecer por meio de frequências de longo-prazo (probabilidade frequentista) ou níveis de credibilidade (probabilidade Bayesiana). Chances é um número positivo real (\(\mathbb{R}^+\)) que mensura também a certeza de um evento. Mas essa certeza não é expressa como uma probabilidade (algo entre 0 e 1), mas como uma razão entre a quantidade de resultados que produzem o evento desejado e a quantidade de resultados que não produzem o evento desejado:

\[ \text{Chances} = \frac{p}{1-p} \]

onde \(p\) é a probabilidade. Chance com o valor de \(1\) é uma chance neutra, algo como uma moeda justa \(p = \frac{1}{2}\). Chances abaixo de \(1\) decrescem a probabilidade de vermos um certo evento e chances acima de \(1\) aumentam a probabilidade do evento. É um sistema de quantificação de incerteza muito usado em casas de apostas. Chances 2:1 quer dizer que a casa de aposta pagará R$ 2 para cada R$ 1 apostado caso o evento apostado ocorra.

Se você revisitar a função logística, verá que ela os coeficientes de \(\boldsymbol{\beta}\) são literalmente o logarítmo da chance (em inglês logodds):

\[ \begin{aligned} p &\sim \text{Logística}(\alpha + \mathbf{X} \boldsymbol{\beta} ) \\ p &\sim \text{Logística}(\alpha) + \text{Logística}( \mathbf{X} \boldsymbol{\beta}) \\ \boldsymbol{\beta} &= \frac{1}{1 + e^{(-\boldsymbol{\beta})}}\\ \boldsymbol{\beta} &= \log(\text{Chance}) \end{aligned} \]

Portanto, os coeficientes de uma regressão logística são expressados em logodds no qual \(0\) é o elemento neutro e qualquer número acima ou abaixo aumenta ou diminui as chances de obtermos um “sucesso” de \(\boldsymbol{y}\). Para termos uma interpretação mais intuitiva (igual a das casas de apostas) precisamos converter as logodds em chances revertendo a função \(\log\). Para isso basta “exponenciar” os valores de \(\boldsymbol{\beta}\):

\[ \text{Chances} = e^{\boldsymbol{\beta}} \]

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

coeff <- exp(model_binomial$coefficients)
coeff
(Intercept)        dist     arsenic       assoc        educ 
       0.85        0.99        1.60        0.88        1.04 

Atividade Prática

Dois datasets estão disponíveis no diretório datasets/:

  1. Titanic Survival: datasets/Titanic_Survival.csv
  2. IBM HR Analytics Employee Attrition & Performance: datasets/IBM_HR_Attrition.csv
###

Ambiente

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] gapminder_0.3.0      ggExtra_0.9          dplyr_1.0.6         
 [4] rstan_2.21.2         StanHeaders_2.21.0-7 MASS_7.3-54         
 [7] ggforce_0.3.3        gganimate_1.0.7      plotly_4.9.3        
[10] carData_3.0-4        DiagrammeR_1.0.6.1   brms_2.15.0         
[13] rstanarm_2.21.1      Rcpp_1.0.6           skimr_2.1.3         
[16] readr_1.4.0          readxl_1.3.1         tibble_3.1.2        
[19] ggplot2_3.3.3        patchwork_1.1.1      cowplot_1.1.1       

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

References

Corrections

If you see mistakes or want to suggest changes, please create an issue on the source repository.

Reuse

Text and figures are licensed under Creative Commons Attribution CC BY-SA 4.0. Source code is available at https://github.com/storopoli/Estatistica-Bayesiana, unless otherwise noted. The figures that have been reused from other sources don't fall under this license and can be recognized by a note in their caption: "Figure from ...".

Citation

For attribution, please cite this work as

Storopoli (2021, Aug. 1). Estatística Bayesiana com R e Stan: Regressão Logística Bayesiana. Retrieved from https://storopoli.github.io/Estatistica-Bayesiana/7-Regressao_Logistica.html

BibTeX citation

@misc{storopoli2021regressaologisticabayesR,
  author = {Storopoli, Jose},
  title = {Estatística Bayesiana com R e Stan: Regressão Logística Bayesiana},
  url = {https://storopoli.github.io/Estatistica-Bayesiana/7-Regressao_Logistica.html},
  year = {2021}
}