rstanarm e brms

As interfaces amigáveis do Stan

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

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 Stan1:

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:

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:

rstanarm

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

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

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 estimado2 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 deviationMAD_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:

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

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 estimado3 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

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

  1. 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↩︎

  2. geralmente chamamos objetos stanreg que são oriundos das funções de modelos Bayesianos do rstanarm como fit.↩︎

  3. geralmente chamamos objetos brmsfit que são oriundos das funções de modelos Bayesianos do brms como fit.↩︎

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: rstanarm e brms. Retrieved from https://storopoli.io/Estatistica-Bayesiana/3-rstanarm.html

BibTeX citation

@misc{storopoli2021priorbayesR,
  author = {Storopoli, Jose},
  title = {Estatística Bayesiana com R e Stan: rstanarm e brms},
  url = {https://storopoli.io/Estatistica-Bayesiana/3-rstanarm.html},
  year = {2021}
}