1 Load packages
library(tidyverse) # data wrangling
library(rstanarm) # Bayes modelling
library(palmerpenguins) # data
library(easystats) # unified API
library(tictoc) # measure computation time2 Motivation
2.1 Bayes and Frequentis
Bayes inference comforts its user with an interpretation that many practitioners seem to prefer (and to understand more quickly), ie., Bayes posterior distribution provides the probability of the hypothesis given the data at hand, \(Pr(H|D)\). Classical Frequentist inference does not provide this probability, but rather the long-run probability of finding an at least as extreme empirical result, given that the Null hypothesis is true. This is all well-known.
To be clear, both approaches have their appeal and may be optimal in different sitaution.
2.2 Technical setup for Bayes analysis provides a barrier
However, a disadvantage of Bayes analysis, at least at its current state, is that it has higher technical and computational demands. For beginners in particular, this may present a substantial (entry) burden. Teaching statistics, I have found that students (and many colleagues) have had difficulties installing Stan (particularly the C++ compiler needed in order to run Stan); Stan is the probabilistic programming language which many front-end Bayes engines use such as brms in R.
Thus, the installation process being not so user-friendly, a burden is placed for beginners which may prevent using Bayes methods.
In that light, this post explores the numerical simarilities of Bayes regression models and Frequentis models. The idea is to use a Frequentist regression model as a proxi for a full Bayesian analysis. The value added is the quick computation and the simple technical setup.
3 Numerical convergence of Bayes and Frequentist approaches
To be clear, the philosophical differences and, consequently, the interpretation of Bayes and Frequentist models are very different.
However, under some instances, the numerical results (for point estimates) may be very similar comparing Bayesian and Frequentist models, and the Frequentist results may thus be interpreted in a Bayesian way (see for a proof https://arxiv.org/abs/1411.5018; for an example, see this paper).
More precisely, a linear regression model will provide the same numerical results as a Bayesian model (ceteris paribus) given that
- the likelihood is normal
- the priors are flat
- data are informative
Note: Flat priors are not non-informative in general.
It is not generally advisable to use flat priors as explained in this post.
In rstanarm, an R package for applied Bayes regression modelling, flat priors (for the regression coefficient) can be specified in this following way:
flat_prior_test <- stan_glm(mpg ~ hp, data = mtcars, prior = NULL, refresh = 0)Check:
prior_summary(flat_prior_test)
#> Priors for model 'flat_prior_test'
#> ------
#> Intercept (after predictors centered)
#> Specified prior:
#> ~ normal(location = 20, scale = 2.5)
#> Adjusted prior:
#> ~ normal(location = 20, scale = 15)
#>
#> Coefficients
#> ~ flat
#>
#> Auxiliary (sigma)
#> Specified prior:
#> ~ exponential(rate = 1)
#> Adjusted prior:
#> ~ exponential(rate = 0.17)
#> ------
#> See help('prior_summary.stanreg') for more details4 Example 1: mctcars
Let’s use mtcars to compute the following linear model: mpg ~ hp.
4.1 Frequentist model
lm_mtcars_freq <- lm(mpg ~ hp, data = mtcars)
parameters(lm_mtcars_freq)
#> Parameter | Coefficient | SE | 95% CI | t(30) | p
#> ------------------------------------------------------------------
#> (Intercept) | 30.10 | 1.63 | [26.76, 33.44] | 18.42 | < .001
#> hp | -0.07 | 0.01 | [-0.09, -0.05] | -6.74 | < .0014.2 Bayesian numerically equivalent model
lm_mtcars_bayes <- stan_glm(mpg ~ hp, data = mtcars, prior = NULL, refresh = 0)
parameters(lm_mtcars_bayes)
#> Parameter | Median | 95% CI | pd | Rhat | ESS | Prior
#> ----------------------------------------------------------------------------------------
#> (Intercept) | 30.08 | [26.79, 33.40] | 100% | 1.000 | 3460.00 | Normal (20.09 +- 15.07)
#> hp | -0.07 | [-0.09, -0.05] | 100% | 1.000 | 3628.00 | Uniform4.3 Conclusion
There’s a small difference due to the fact that the Bayes model involves random numbers. Beside that, the numerical results (of the point estimates) are (virtually) identical.
We have not put flat priors on the intercept and on sigma, hence those parameters may differ from the frequentist ones.
5 Example 2: penguins
data(penguins)5.1 Frequentist model
lm_penguins_freq <- lm(body_mass_g ~ bill_length_mm, data = penguins)
parameters(lm_penguins_freq)
#> Parameter | Coefficient | SE | 95% CI | t(340) | p
#> ---------------------------------------------------------------------------
#> (Intercept) | 362.31 | 283.35 | [-195.02, 919.64] | 1.28 | 0.202
#> bill length mm | 87.42 | 6.40 | [ 74.82, 100.01] | 13.65 | < .0015.2 Bayesian numerically equivalent model
lm_penguins_bayes <- stan_glm(body_mass_g ~ bill_length_mm, data = penguins, prior = NULL, refresh = 0)
parameters(lm_penguins_bayes)
#> Parameter | Median | 95% CI | pd | Rhat | ESS | Prior
#> ----------------------------------------------------------------------------------------------------
#> (Intercept) | 357.45 | [-180.90, 905.48] | 90.08% | 0.999 | 4579.00 | Normal (4201.75 +- 2004.89)
#> bill_length_mm | 87.55 | [ 75.17, 99.76] | 100% | 0.999 | 4586.00 | Uniform5.3 Conclusion
Again, the numerical results are very similar.
6 Example 3: diamonds
data(diamonds)6.1 Frequentist model
lm_diamonds_freq <- lm(price ~ carat, data = diamonds)
parameters(lm_diamonds_freq)
#> Parameter | Coefficient | SE | 95% CI | t(53938) | p
#> ----------------------------------------------------------------------------
#> (Intercept) | -2256.36 | 13.06 | [-2281.95, -2230.77] | -172.83 | < .001
#> carat | 7756.43 | 14.07 | [ 7728.86, 7784.00] | 551.41 | < .0016.2 Bayesian numerically equivalent model
tic()
lm_diamonds_bayes <- stan_glm(price ~ carat, data = diamonds, prior = NULL, refresh = 0)
toc()
#> 24.06 sec elapsed
parameters(lm_diamonds_bayes)
#> Parameter | Median | 95% CI | pd | Rhat | ESS | Prior
#> ----------------------------------------------------------------------------------------------------
#> (Intercept) | -2255.92 | [-2283.22, -2229.19] | 100% | 1.000 | 2911.00 | Normal (3932.80 +- 9973.60)
#> carat | 7756.16 | [ 7728.01, 7784.72] | 100% | 1.000 | 2730.00 | Uniform6.3 Conclusion
Note the prolonged computation time for the Bayes model due to the larger sample size (and possibly due to the unconstrained, ie., flat, prior).
Again, the regression coefficients are very similar.
7 Caveats
The less data you have, the more the Posterior distribution will depend on your priors. In other words, make sure that your sample size is “large enough”. This post provides some deeper insights.
We have not specified the prior_intercept in our Stan models, as we were not interested in that. However, using the same argument as above, we could set the intercept prior to a flat prior, too. Read here more about setting priors in rstanarm.
8 Concluding remarks
As can be seen (and as it is stated in the literature), the point estimates of Bayesian and Frequentist models converge given that a flat prior is used and that the likelihood is normal. Furthermore, the data must be informative enough to “convince” the posterior distribution.
Hence, using lm can be seen as second-best option for computing simple Bayesian models.
Remember that the idea of this post was to provide a simple soluttion practitioners who do not want to install additional softare (such as Stan). Of course, using a second-best option is not optimal by definition.
9 Further reading
This post comments on the numerical equivalence of flat prior Bayesian models with Frquentis models.
10 Reproducibility
#> ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
#> setting value
#> version R version 4.2.1 (2022-06-23)
#> os macOS Big Sur ... 10.16
#> system x86_64, darwin17.0
#> ui X11
#> language (EN)
#> collate en_US.UTF-8
#> ctype en_US.UTF-8
#> tz Europe/Berlin
#> date 2024-01-29
#> pandoc 3.1.1 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/ (via rmarkdown)
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
#> package * version date (UTC) lib source
#> abind 1.4-5 2016-07-21 [1] CRAN (R 4.2.0)
#> backports 1.4.1 2021-12-13 [1] CRAN (R 4.2.0)
#> base64enc 0.1-3 2015-07-28 [1] CRAN (R 4.2.0)
#> bayesplot 1.10.0 2022-11-16 [1] CRAN (R 4.2.0)
#> bayestestR * 0.13.1 2023-04-07 [1] CRAN (R 4.2.0)
#> blogdown 1.18 2023-06-19 [1] CRAN (R 4.2.0)
#> bookdown 0.36 2023-10-16 [1] CRAN (R 4.2.0)
#> boot 1.3-28.1 2022-11-22 [1] CRAN (R 4.2.0)
#> bslib 0.5.1 2023-08-11 [1] CRAN (R 4.2.0)
#> cachem 1.0.8 2023-05-01 [1] CRAN (R 4.2.0)
#> callr 3.7.3 2022-11-02 [1] CRAN (R 4.2.0)
#> checkmate 2.2.0 2023-04-27 [1] CRAN (R 4.2.0)
#> cli 3.6.1 2023-03-23 [1] CRAN (R 4.2.0)
#> coda 0.19-4 2020-09-30 [1] CRAN (R 4.2.0)
#> codetools 0.2-19 2023-02-01 [1] CRAN (R 4.2.0)
#> colorout * 1.3-0 2023-11-08 [1] Github (jalvesaq/colorout@8384882)
#> colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.2.0)
#> colourpicker 1.3.0 2023-08-21 [1] CRAN (R 4.2.0)
#> correlation * 0.8.4 2023-04-06 [1] CRAN (R 4.2.1)
#> crayon 1.5.2 2022-09-29 [1] CRAN (R 4.2.1)
#> crosstalk 1.2.1 2023-11-23 [1] CRAN (R 4.2.1)
#> curl 5.1.0 2023-10-02 [1] CRAN (R 4.2.0)
#> datawizard * 0.9.0 2023-09-15 [1] CRAN (R 4.2.0)
#> devtools 2.4.5 2022-10-11 [1] CRAN (R 4.2.1)
#> digest 0.6.33 2023-07-07 [1] CRAN (R 4.2.0)
#> distributional 0.3.2 2023-03-22 [1] CRAN (R 4.2.0)
#> dplyr * 1.1.3 2023-09-03 [1] CRAN (R 4.2.0)
#> DT 0.30 2023-10-05 [1] CRAN (R 4.2.0)
#> dygraphs 1.1.1.6 2018-07-11 [1] CRAN (R 4.2.0)
#> easystats * 0.7.0 2023-11-05 [1] CRAN (R 4.2.1)
#> effectsize * 0.8.6 2023-09-14 [1] CRAN (R 4.2.0)
#> ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.2.0)
#> emmeans 1.8.9 2023-10-17 [1] CRAN (R 4.2.0)
#> estimability 1.4.1 2022-08-05 [1] CRAN (R 4.2.0)
#> evaluate 0.21 2023-05-05 [1] CRAN (R 4.2.0)
#> fansi 1.0.5 2023-10-08 [1] CRAN (R 4.2.0)
#> farver 2.1.1 2022-07-06 [1] CRAN (R 4.2.0)
#> fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.2.0)
#> forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.2.0)
#> fs 1.6.3 2023-07-20 [1] CRAN (R 4.2.0)
#> generics 0.1.3 2022-07-05 [1] CRAN (R 4.2.0)
#> ggplot2 * 3.4.4 2023-10-12 [1] CRAN (R 4.2.0)
#> glue 1.6.2 2022-02-24 [1] CRAN (R 4.2.0)
#> gridExtra 2.3 2017-09-09 [1] CRAN (R 4.2.0)
#> gtable 0.3.4 2023-08-21 [1] CRAN (R 4.2.0)
#> gtools 3.9.4 2022-11-27 [1] CRAN (R 4.2.0)
#> hms 1.1.3 2023-03-21 [1] CRAN (R 4.2.0)
#> htmltools 0.5.6.1 2023-10-06 [1] CRAN (R 4.2.0)
#> htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.2.0)
#> httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.2.0)
#> igraph 1.5.1 2023-08-10 [1] CRAN (R 4.2.0)
#> inline 0.3.19 2021-05-31 [1] CRAN (R 4.2.0)
#> insight * 0.19.7 2023-11-26 [1] CRAN (R 4.2.1)
#> jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.2.0)
#> jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.2.0)
#> knitr 1.45 2023-10-30 [1] CRAN (R 4.2.1)
#> later 1.3.1 2023-05-02 [1] CRAN (R 4.2.0)
#> lattice 0.21-8 2023-04-05 [1] CRAN (R 4.2.0)
#> lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.2.1)
#> lme4 1.1-34 2023-07-04 [1] CRAN (R 4.2.0)
#> loo 2.6.0 2023-03-31 [1] CRAN (R 4.2.0)
#> lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.2.0)
#> magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.2.0)
#> markdown 1.11 2023-10-19 [1] CRAN (R 4.2.0)
#> MASS 7.3-60 2023-05-04 [1] CRAN (R 4.2.0)
#> Matrix 1.5-4.1 2023-05-18 [1] CRAN (R 4.2.0)
#> matrixStats 0.63.0 2022-11-18 [1] CRAN (R 4.2.0)
#> memoise 2.0.1 2021-11-26 [1] CRAN (R 4.2.0)
#> mime 0.12 2021-09-28 [1] CRAN (R 4.2.0)
#> miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.2.0)
#> minqa 1.2.5 2022-10-19 [1] CRAN (R 4.2.1)
#> modelbased * 0.8.6 2023-01-13 [1] CRAN (R 4.2.1)
#> multcomp 1.4-25 2023-06-20 [1] CRAN (R 4.2.0)
#> munsell 0.5.0 2018-06-12 [1] CRAN (R 4.2.0)
#> mvtnorm 1.2-2 2023-06-08 [1] CRAN (R 4.2.0)
#> nlme 3.1-162 2023-01-31 [1] CRAN (R 4.2.0)
#> nloptr 2.0.3 2022-05-26 [1] CRAN (R 4.2.0)
#> palmerpenguins * 0.1.1 2022-08-15 [1] CRAN (R 4.2.0)
#> parameters * 0.21.3 2023-11-02 [1] CRAN (R 4.2.1)
#> performance * 0.10.8 2023-10-30 [1] CRAN (R 4.2.1)
#> pillar 1.9.0 2023-03-22 [1] CRAN (R 4.2.0)
#> pkgbuild 1.4.0 2022-11-27 [1] CRAN (R 4.2.0)
#> pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.2.0)
#> pkgload 1.3.2.1 2023-07-08 [1] CRAN (R 4.2.0)
#> plyr 1.8.8 2022-11-11 [1] CRAN (R 4.2.0)
#> posterior 1.5.0 2023-10-31 [1] CRAN (R 4.2.1)
#> prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.2.0)
#> processx 3.8.2 2023-06-30 [1] CRAN (R 4.2.0)
#> profvis 0.3.8 2023-05-02 [1] CRAN (R 4.2.0)
#> promises 1.2.1 2023-08-10 [1] CRAN (R 4.2.0)
#> ps 1.7.5 2023-04-18 [1] CRAN (R 4.2.0)
#> purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.2.0)
#> QuickJSR 1.0.5 2023-09-03 [1] CRAN (R 4.2.0)
#> R6 2.5.1 2021-08-19 [1] CRAN (R 4.2.0)
#> Rcpp * 1.0.11 2023-07-06 [1] CRAN (R 4.2.0)
#> RcppParallel 5.1.7 2023-02-27 [1] CRAN (R 4.2.0)
#> readr * 2.1.4 2023-02-10 [1] CRAN (R 4.2.0)
#> remotes 2.4.2.1 2023-07-18 [1] CRAN (R 4.2.0)
#> report * 0.5.8 2023-12-07 [1] CRAN (R 4.2.1)
#> reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.2.0)
#> rlang 1.1.1 2023-04-28 [1] CRAN (R 4.2.0)
#> rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.2.0)
#> rstan 2.32.3 2023-10-15 [1] CRAN (R 4.2.0)
#> rstanarm * 2.26.1 2023-09-13 [1] CRAN (R 4.2.0)
#> rstantools 2.3.1.1 2023-07-18 [1] CRAN (R 4.2.0)
#> rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.2.0)
#> sandwich 3.0-2 2022-06-15 [1] CRAN (R 4.2.0)
#> sass 0.4.7 2023-07-15 [1] CRAN (R 4.2.0)
#> scales 1.2.1 2022-08-20 [1] CRAN (R 4.2.0)
#> see * 0.8.1 2023-11-03 [1] CRAN (R 4.2.1)
#> sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.2.0)
#> shiny 1.8.0 2023-11-17 [1] CRAN (R 4.2.1)
#> shinyjs 2.1.0 2021-12-23 [1] CRAN (R 4.2.0)
#> shinystan 2.6.0 2022-03-03 [1] CRAN (R 4.2.0)
#> shinythemes 1.2.0 2021-01-25 [1] CRAN (R 4.2.0)
#> StanHeaders 2.26.28 2023-09-07 [1] CRAN (R 4.2.0)
#> stringi 1.7.12 2023-01-11 [1] CRAN (R 4.2.0)
#> stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.2.1)
#> survival 3.5-5 2023-03-12 [1] CRAN (R 4.2.0)
#> tensorA 0.36.2 2020-11-19 [1] CRAN (R 4.2.0)
#> TH.data 1.1-2 2023-04-17 [1] CRAN (R 4.2.0)
#> threejs 0.3.3 2020-01-21 [1] CRAN (R 4.2.0)
#> tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.2.0)
#> tictoc * 1.2 2023-04-23 [1] CRAN (R 4.2.0)
#> tidyr * 1.3.0 2023-01-24 [1] CRAN (R 4.2.0)
#> tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.2.0)
#> tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.2.0)
#> timechange 0.2.0 2023-01-11 [1] CRAN (R 4.2.0)
#> tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.2.0)
#> urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.2.0)
#> usethis 2.2.2 2023-07-06 [1] CRAN (R 4.2.0)
#> utf8 1.2.3 2023-01-31 [1] CRAN (R 4.2.0)
#> V8 4.4.0 2023-10-09 [1] CRAN (R 4.2.0)
#> vctrs 0.6.4 2023-10-12 [1] CRAN (R 4.2.0)
#> withr 2.5.2 2023-10-30 [1] CRAN (R 4.2.1)
#> xfun 0.40 2023-08-09 [1] CRAN (R 4.2.0)
#> xtable 1.8-4 2019-04-21 [1] CRAN (R 4.2.0)
#> xts 0.13.1 2023-04-16 [1] CRAN (R 4.2.0)
#> yaml 2.3.7 2023-01-23 [1] CRAN (R 4.2.0)
#> zoo 1.8-12 2023-04-13 [1] CRAN (R 4.2.0)
#>
#> [1] /Users/sebastiansaueruser/Rlibs
#> [2] /Library/Frameworks/R.framework/Versions/4.2/Resources/library
#>
#> ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────