I’ve run into a weird issue while designing teaching materials. Prior predictive simulation in brms is not demonstrating the behaviour that I expected. I am looking to better understand what brms is doing or to clear up a conceptual misunderstanding on my end.
The example I’m creating concerns coding categorical variables with dummy coding and how that can create problems when setting priors. To keep things concrete, I will walk through a reproducible example. Below I create a variable ‘x’ that has two levels: 0 and 1. We then create
y as a continuous variable from a normal distribution where the mean is conditional on the level of x.
library(brms) # Simulating data set.seed(123) x <- rep(c(0, 1), times = 50) y <- rnorm(100, 10+(x*2), 1) # Creating dataframe and changing x to factor df <- data.frame(x, y) df$x <- as.factor(df$x)
As the level ‘0’ will be represented by the Intercept, and the estimated effect of ‘1’ will be the difference from the Intercept, I would always expect any set of priors to have greater uncertainty regarding the average of y when x=1, because by definition it has the uncertainty of the prior on the Intercept and the uncertainty of the difference between them. I checked that was the case with a simple simulation below:
# Running a simple prior simulation x0 <- rnorm(1e6, mean = 8, sd = 3) x1 <- rnorm(1e6, mean = 8, sd = 3) + rnorm(1e6, mean = 0, sd = 2) # Means are equal as expected mean(x0) # 7.998458 mean(x1) # 7.995305 # SD of x1 is more uncertain, as expected sd(x0) # 2.999817 sd(x1) # 3.605881
However, when I run a prior check using the
brm() function, I don’t see that the uncertainty regarding the averages is different. I run the prior check using the code below.
prior_check <- brm(y ~ x, data = df, prior = c( prior(normal(8,3), class = "Intercept"), prior(normal(0,2), class = "b"), prior(normal(0,3), class = "sigma")), sample_prior = "only")
First I visualized the means using the following code, which produced the plot below it.
The uncertainty regarding the means looks about the same, which I didn’t understand, so I looked at the draws myself, which led to the same result:
# Getting a dataframe of draws prior_df <- as_draws_df(prior_check) # Getting the means - similar as expected mean(prior_df$b_Intercept) # 7.97645 mean(prior_df$b_Intercept+prior_df$b_x1) # 8.072909 # Getting the sds, similar contrary to expectations sd(prior_df$b_Intercept) # 3.175866 sd(prior_df$b_Intercept+prior_df$b_x1) # 3.11842
All of the materials I’ve read, as well as many posts on this website and others, refer to this problem of setting priors on dummy variables that I don’t seem to be encountering in brms. Any insight into why this is happening in brms and/or what I do not understand about setting priors would be awesome!
sessionInfo() R version 4.1.1 (2021-08-10) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045) Matrix products: default locale:  LC_COLLATE=English_Canada.1252 LC_CTYPE=English_Canada.1252  LC_MONETARY=English_Canada.1252 LC_NUMERIC=C  LC_TIME=English_Canada.1252 attached base packages:  stats graphics grDevices utils datasets methods base other attached packages:  brms_2.18.0 Rcpp_1.0.10 loaded via a namespace (and not attached):  nlme_3.1-152 matrixStats_0.60.1 xts_0.12.1  threejs_0.3.3 rstan_2.21.2 tensorA_0.36.2  tools_4.1.1 backports_1.2.1 utf8_1.2.2  R6_2.5.1 DT_0.19 mgcv_1.8-36  projpred_2.0.2 DBI_1.1.1 colorspace_2.0-2  withr_2.5.0 tidyselect_1.2.0 gridExtra_2.3  prettyunits_1.1.1 processx_3.5.2 Brobdingnag_1.2-6  emmeans_1.6.3 curl_4.3.2 compiler_4.1.1  cli_3.2.0 shinyjs_2.0.0 sandwich_3.0-1  labeling_0.4.2 colourpicker_1.1.0 posterior_1.2.1  scales_1.2.1 dygraphs_220.127.116.11 checkmate_2.0.0  mvtnorm_1.1-2 ggridges_0.5.3 callr_3.7.0  stringr_1.5.0 digest_0.6.27 StanHeaders_2.21.0-7  minqa_1.2.4 base64enc_0.1-3 pkgconfig_2.0.3  htmltools_0.5.2 lme4_1.1-29 fastmap_1.1.0  htmlwidgets_1.5.3 rlang_1.0.6 rstudioapi_0.13  shiny_1.6.0 farver_2.1.0 generics_0.1.0  zoo_1.8-9 jsonlite_1.7.2 crosstalk_1.1.1  gtools_3.9.2 dplyr_1.0.8 distributional_0.2.2  inline_0.3.19 magrittr_2.0.3 loo_2.4.1  bayesplot_1.8.1 Matrix_1.5-3 munsell_0.5.0  fansi_0.5.0 abind_1.4-5 lifecycle_1.0.3  stringi_1.6.2 multcomp_1.4-17 MASS_7.3-54  pkgbuild_1.2.0 plyr_1.8.6 grid_4.1.1  parallel_4.1.1 promises_18.104.22.168 crayon_1.4.1  miniUI_0.1.1.1 lattice_0.20-44 splines_4.1.1  ps_1.6.0 pillar_1.8.1 igraph_1.2.6  boot_1.3-28 markdown_1.1 estimability_1.3  shinystan_2.5.0 reshape2_1.4.4 codetools_0.2-18  stats4_4.1.1 rstantools_2.1.1 glue_1.6.2  V8_3.4.2 RcppParallel_5.1.4 nloptr_22.214.171.124  vctrs_0.5.2 httpuv_1.6.2 gtable_0.3.0  assertthat_0.2.1 ggplot2_3.4.1 mime_0.11  xtable_1.8-4 coda_0.19-4 later_1.3.0  survival_3.2-11 rsconnect_0.8.24 tibble_3.1.4  shinythemes_1.2.0 gamm4_0.2-6 TH.data_1.0-10  ellipsis_0.3.2 bridgesampling_1.1-2 ```