Estatística Bayesiana com R e Stan: rstanarm e brms

Jose Storopoli

A principal ferramenta para computação Bayesiana é a linguagem probabilística Stan (Carpenter et al., 2017). O problema do Stan é que ele é uma linguagem de programação e, portanto, possui um acesso dificultado a não-programadores. Abaixo um código que mostra como é um programa escrito em Stan¹:

data {
  int<lower=0> N;
  vector<lower=0, upper=200>[N] kid_score;
  vector<lower=0, upper=200>[N] mom_iq;
}
parameters {
  real alpha;
  real beta;
  real<lower=0> sigma;
}
model {
  sigma ~ cauchy(0, 2.5);
  kid_score ~ normal(alpha + beta * mom_iq, sigma);
}

Para remediar essa barreira de acesso ao Stan, temos interfaces abstratas que interpretam a intenção do usuário e lidam com a parte mais obral de codificação. As principais são rstanarm (Goodrich, Gabry, Ali, & Brilleman, 2020) e brms (Bürkner, 2017). Ambos usa a mesma síntaxe e as mesmas famílias de funções de verossimilhança, mas com algumas diferenças:

rstanarm todos os modelos são pré-compilados e brms não possui os modelos pré-compilados então os modelos devem ser todos compilados antes de serem rodados. A diferença prática é que você irá esperar alguns instantes antes do R começar a amostrar do modelo.
Como brms compila os modelos conforme a especificação do usuário (não possui modelos pré-compilados) ele pode criar modelos um pouco mais eficientes que o rstanarm. Caso o tempo de compilação dos modelos brms seja menor que ganho de velocidade em amostragem do modelo, é vantajoso especificar seu modelo com brms.
brms dá maior poder e flexibilidade ao usuário na especificação de funções de verossimilhança e também permite modelos mais complexos que o rstanarm (um exemplo notório é que rstanarm não permite usarmos distribuição \(t\) de Student como função de verossimilhança enquanto que isso é possível no brms. Portanto a Aula 9 - Regressão Robusta usará exclusivamente o brms).

Fórmulas

Todos os modelos especificados pelo rstanarm e brms usam uma fórmula com a seguinte síntaxe:

dependente ~ independente_1 + independente_2 + ...

Para quem conhece R essa síntaxe é a mesma utilizada em regressões lineares frequentista com o lm() ou em modelos lineares generalizados frequentistas com glm().

`family`

Todo modelo especificado pelo rstanarm e brms devem especificar qual família da função de verossimilhança (family) respectivamente com a função de ligação (‘link’) que fará o mapeamento dos parâmetros condicionados nos dados para a variável dependente. Caso o usuário não designe nenhum valor para esses dois parâmetros, rstanarm e brms usarão a verossimilhança Gaussiana (family = gaussian) e a função de identidade como função de ligação (link = "identity").

Estes argumentos family e link são conhecidos para quem usa R para estatística frequentista pois são os mesmos da função glm() de R para modelos lineares generalizados frequentistas. Algumas das principais famílias com suas funções de ligação padrões são:

Gaussiana – family = gaussian – link = "identity"
Log-Normal – family = lognormal – link = "log"
Binomial – family = binomial – link = "logit"
Poisson – family = poisson – link = "log"
Binomial Negativa – family = negbinomial – link = "log"
\(t\) de Student – family = student – link = "identity"
Exponencial – family = exponential – link = "log"

`rstanarm`

O rstanarm é a porta de entrada para estatística Bayesiana com Stan. O nome rstanarm é:

r: pacote para R
stan: usa a linguagem probabilística Stan
arm: acrônimo para Applied Regression Modeling

Ele possui as seguintes funções para especificação de modelos Bayesianos:

stan_glm() – modelos lineares generalizados (generalized linear model)
stan_lm() – modelos lineares regularizados (linear model)
stan_aov() – modelo ANOVA (analysis of variance)
stan_glmer() – modelos linares generalizados multiníveis
stan_lmer() – modelos linares regularizados multiníveis
stan_jm() – modelos longitudinais e de sobrevivência
stan_nlmer() – modelos não-lineares multiníveis (non-linear model)
stan_polr() – modelos ordinais
stan_gamm4() – modelos aditivos linares multiníveis

Neste curso usaremos apenas stan_glm e stan_glmer, mas saiba que você possui uma vasta categoria de modelos bayesianos à disposição.

Exemplo usando o `mtcars`

Vamos estimar modelos Bayesianos usando o dataset já conhecido mtcars da Aula 1 - Comandos Básicos de R. A idéia é usarmos como variável dependente a autonomia do carro (milhas por galão – mpg) e como independentes o peso do carro (wt) e se o carro é automático ou manual (am). A fórmula então fica:

mpg ~ wt + am

Note que também devemos especificar a localização dos nossos dados com o argumento data. Para garantir que toda a funcionalidade do rstanarm esteja disponível para seu modelo, recomendo que especifique sempre o valor de data como um objeto data.frame ou tibble.

library(rstanarm)
rstanarm_fit <- stan_glm(mpg ~ wt + am, data = mtcars)

Podemos ver o sumário do modelo estimado² com a função print() ou summary():

print(rstanarm_fit)

stan_glm
 family:       gaussian [identity]
 formula:      mpg ~ wt + am
 observations: 32
 predictors:   3
------
            Median MAD_SD
(Intercept) 37.3    3.2  
wt          -5.3    0.8  
am           0.0    1.5  

Auxiliary parameter(s):
      Median MAD_SD
sigma 3.2    0.4   

------
* For help interpreting the printed output see ?print.stanreg
* For info on the priors used see ?prior_summary.stanreg

summary(rstanarm_fit)


Model Info:
 function:     stan_glm
 family:       gaussian [identity]
 formula:      mpg ~ wt + am
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 32
 predictors:   3

Estimates:
              mean   sd   10%   50%   90%
(Intercept) 37.2    3.2 33.2  37.3  41.1 
wt          -5.3    0.8 -6.3  -5.3  -4.3 
am           0.0    1.6 -2.0   0.0   2.1 
sigma        3.2    0.4  2.7   3.2   3.8 

Fit Diagnostics:
           mean   sd   10%   50%   90%
mean_PPD 20.1    0.8 19.1  20.1  21.1 

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.1  1.0  1945 
wt            0.0  1.0  2085 
am            0.0  1.0  1949 
sigma         0.0  1.0  2793 
mean_PPD      0.0  1.0  3434 
log-posterior 0.0  1.0  1394 

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).

A interpretação e significado da saída dos modelos rstanarm serão explicadas nas aulas seguintes. A função print() é mais concisa que a função summary(). Para quem já rodou modelos de regressão, o output da função print() mostra a mediana (Median) dos coeficientes e erro (sigma) do modelo junto com o desvio absoluto médio (mean absolute deviation – MAD_SD). No output da função summary(), a tabela Estimates é a tabela que mostra a média (mean) dos coeficientes e erro (sigma) do modelo junto com o desvio padrão (sd) e os percentis 10%, 50% (mediana) e 90% baseados na mediana e desvio absoluto médio.

Podemos também especificar os percentis desejados no summary():

summary(rstanarm_fit, probs = c(0.025, 0.975))


Model Info:
 function:     stan_glm
 family:       gaussian [identity]
 formula:      mpg ~ wt + am
 algorithm:    sampling
 sample:       4000 (posterior sample size)
 priors:       see help('prior_summary')
 observations: 32
 predictors:   3

Estimates:
              mean   sd   2.5%   97.5%
(Intercept) 37.2    3.2 30.9   43.4   
wt          -5.3    0.8 -6.9   -3.7   
am           0.0    1.6 -3.1    3.3   
sigma        3.2    0.4  2.5    4.2   

Fit Diagnostics:
           mean   sd   2.5%   97.5%
mean_PPD 20.1    0.8 18.5   21.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.1  1.0  1945 
wt            0.0  1.0  2085 
am            0.0  1.0  1949 
sigma         0.0  1.0  2793 
mean_PPD      0.0  1.0  3434 
log-posterior 0.0  1.0  1394 

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).

`brms`

O brms alia toda a comodidade do rstanarm com o poder e flexibilidade do Stan. O nome brms quer dizer:

b: Bayesian
r: regression
m: models
s: usando Stan

Ao invés de possuir diversas funções para diferentes tipos de modelo, brms tem apenas uma única função para especificar modelos: brm() – acrônimo para Bayesian Regression Model. O usuário consegue especificar qualquer modelo que quiser a partir da função brm() apenas mudando seus parâmetros internos. Os parâmetros da função brm() mais importantes são (como mencionado acima):

family – especifica a família da função de verossimilhança do modelo (padrão gaussian)
link – especifica a função de ligação que fará o mapeamento dos parâmetros condicionados nos dados para a variável dependente do modelo (padrão varia conforme o family, mas para family = gaussian, a função identidade é a função de ligação padrão link = "identity")

Vamos usar o mesmo modelo que usamos para o rstanarm acima. Note que a fórmula não muda e usaremos a mesma:

mpg ~ wt + am

Note que também devemos especificar a localização dos nossos dados com o argumento data. Para garantir que toda a funcionalidade do brms esteja disponível para seu modelo, recomendo que especifique sempre o valor de data como um objeto data.frame ou tibble.

library(brms)
brms_fit <- brm(mpg ~ wt + am, data = mtcars)

Podemos ver o sumário do modelo estimado³ com a função print() ou summary(). No caso do brms não há diferença entre elas e elas literalmente produzem o mesmo output:

print(brms_fit)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: mpg ~ wt + am 
   Data: mtcars (Number of observations: 32) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    37.27      3.21    30.88    43.79 1.00     2340     2513
wt           -5.34      0.83    -6.98    -3.68 1.00     2355     2334
am           -0.01      1.63    -3.18     3.15 1.00     2474     2639

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.21      0.44     2.50     4.18 1.00     2760     2835

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

summary(brms_fit)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: mpg ~ wt + am 
   Data: mtcars (Number of observations: 32) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    37.27      3.21    30.88    43.79 1.00     2340     2513
wt           -5.34      0.83    -6.98    -3.68 1.00     2355     2334
am           -0.01      1.63    -3.18     3.15 1.00     2474     2639

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.21      0.44     2.50     4.18 1.00     2760     2835

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Note que temos quase que o mesmo output, sendo que brms por padrão mostra a média (Estimate) dos coeficientes e erro (sigma) do modelo junto com o desvio padrão (Est.Error) e os percentis 5% (l-95%) e 95% (u-95%) baseados na média e desvio padrão. Caso queira mediana e desvio absoluto médio, forneça o argumento robust = TRUE para a função print():

summary(brms_fit, robust = TRUE)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: mpg ~ wt + am 
   Data: mtcars (Number of observations: 32) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept    37.34      3.18    30.88    43.79 1.00     2340     2513
wt           -5.36      0.82    -6.98    -3.68 1.00     2355     2334
am           -0.04      1.62    -3.18     3.15 1.00     2474     2639

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma     3.17      0.42     2.50     4.18 1.00     2760     2835

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

Note que agora temos valores diferentes para o output de summary() com robust = TRUE, mas as colunas são as mesmas. Não se engane, agora temos a mediana (Estimate) dos coeficientes e erro (sigma) do modelo junto com o desvio absoluto médio (Est.Error) e os percentis 5% (l-95%) e 95% (u-95%) baseados na mediana e desvio absoluto médio.

Podemos também especificar os percentis desejados no summary(). Aqui a síntaxe é um pouco diferente da síntaxe do rstanarm:

summary(brms_fit, prob = 0.9)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: mpg ~ wt + am 
   Data: mtcars (Number of observations: 32) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Population-Level Effects: 
          Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
Intercept    37.27      3.21    31.94    42.42 1.00     2340     2513
wt           -5.34      0.83    -6.69    -3.98 1.00     2355     2334
am           -0.01      1.63    -2.70     2.69 1.00     2474     2639

Family Specific Parameters: 
      Estimate Est.Error l-90% CI u-90% CI Rhat Bulk_ESS Tail_ESS
sigma     3.21      0.44     2.59     4.02 1.00     2760     2835

Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

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] brms_2.15.0     rstanarm_2.21.1 Rcpp_1.0.6      skimr_2.1.3    
 [5] readr_1.4.0     readxl_1.3.1    tibble_3.1.2    ggplot2_3.3.3  
 [9] 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         repr_1.1.3           splines_4.1.0       
  [7] crosstalk_1.1.1      rstantools_2.1.1     inline_0.3.19       
 [10] digest_0.6.27        htmltools_0.5.1.1    magick_2.7.2        
 [13] rsconnect_0.8.18     fansi_0.5.0          magrittr_2.0.1      
 [16] RcppParallel_5.1.4   matrixStats_0.59.0   xts_0.12.1          
 [19] prettyunits_1.1.1    colorspace_2.0-1     textshaping_0.3.4   
 [22] xfun_0.23            dplyr_1.0.6          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] gtable_0.3.0         V8_3.4.2             pkgbuild_1.2.0      
 [34] rstan_2.21.2         abind_1.4-5          scales_1.1.1        
 [37] mvtnorm_1.1-1        DBI_1.1.1            miniUI_0.1.1.1      
 [40] xtable_1.8-4         stats4_4.1.0         StanHeaders_2.21.0-7
 [43] DT_0.18              htmlwidgets_1.5.3    threejs_0.3.3       
 [46] RColorBrewer_1.1-2   ellipsis_0.3.2       pkgconfig_2.0.3     
 [49] loo_2.4.1            farver_2.1.0         sass_0.4.0          
 [52] utf8_1.2.1           tidyselect_1.1.1     labeling_0.4.2      
 [55] rlang_0.4.11         reshape2_1.4.4       later_1.2.0         
 [58] munsell_0.5.0        cellranger_1.1.0     tools_4.1.0         
 [61] cli_2.5.0            generics_0.1.0       ggridges_0.5.3      
 [64] evaluate_0.14        stringr_1.4.0        fastmap_1.1.0       
 [67] yaml_2.2.1           ragg_1.1.2           processx_3.5.2      
 [70] knitr_1.33           purrr_0.3.4          nlme_3.1-152        
 [73] mime_0.10            projpred_2.0.2       xml2_1.3.2          
 [76] compiler_4.1.0       bayesplot_1.8.0      shinythemes_1.2.0   
 [79] rstudioapi_0.13      curl_4.3.1           gamm4_0.2-6         
 [82] png_0.1-7            bslib_0.2.5.1        stringi_1.6.2       
 [85] highr_0.9            ps_1.6.0             Brobdingnag_1.2-6   
 [88] lattice_0.20-44      Matrix_1.3-3         nloptr_1.2.2.2      
 [91] markdown_1.1         shinyjs_2.0.0        vctrs_0.3.8         
 [94] pillar_1.6.1         lifecycle_1.0.0      jquerylib_0.1.4     
 [97] bridgesampling_1.1-2 httpuv_1.6.1         R6_2.5.0            
[100] bookdown_0.22        promises_1.2.0.1     gridExtra_2.3       
[103] codetools_0.2-18     distill_1.2          boot_1.3-28         
[106] colourpicker_1.1.0   MASS_7.3-54          gtools_3.8.2        
[109] assertthat_0.2.1     rprojroot_2.0.2      withr_2.4.2         
[112] shinystan_2.5.0      mgcv_1.8-35          parallel_4.1.0      
[115] hms_1.1.0            grid_4.1.0           tidyr_1.1.3         
[118] coda_0.19-4          minqa_1.2.4          rmarkdown_2.8       
[121] downlit_0.2.1        shiny_1.6.0          lubridate_1.7.10    
[124] base64enc_0.1-3      dygraphs_1.1.1.6

Bürkner, P.-C. (2017). brms: An R package for Bayesian multilevel models using Stan. Journal of Statistical Software, 80(1), 1–28. https://doi.org/10.18637/jss.v080.i01

Carpenter, B., Gelman, A., Hoffman, M. D., Lee, D., Goodrich, B., Betancourt, M., … Riddell, A. (2017). Stan : A Probabilistic Programming Language. Journal of Statistical Software, 76(1). https://doi.org/10.18637/jss.v076.i01

Goodrich, B., Gabry, J., Ali, I., & Brilleman, S. (2020). Rstanarm: Bayesian applied regression modeling via Stan. Retrieved from https://mc-stan.org/rstanarm

notem que a síntaxe é bem similar à C++ com uma diferença notória que pontos flutuantes double são real e o ~ designa uma distribuição probabilística a uma variável↩︎
geralmente chamamos objetos stanreg que são oriundos das funções de modelos Bayesianos do rstanarm como fit.↩︎
geralmente chamamos objetos brmsfit que são oriundos das funções de modelos Bayesianos do brms como fit.↩︎

rstanarm e brms

Fórmulas

`family`

`rstanarm`

Exemplo usando o `mtcars`

`brms`

Ambiente

References

Corrections

Reuse

Citation

rstanarm e brms

Fórmulas

family

rstanarm

Exemplo usando o mtcars

brms

Ambiente

References

Corrections

Reuse

Citation

`family`

`rstanarm`

Exemplo usando o `mtcars`

`brms`