Negative prior distribution

I’m working my way through the exercises in Richard McElreath’s Statistical Rethinking. In one exercise, the ask is to replicate a model from the chapter and then re-estimate with new priors. The original model can be fit in {brms} as follows:

# remotes::install_github("rmcelreath/rethinking")
library(rethinking)
library(brms)

data(tulips)
d <- tulips

d$blooms_std <- d$blooms / max(d$blooms)
d$water_cent <- d$water - mean(d$water)
d$shade_cent <- d$shade - mean(d$shade)

mod <- brm(blooms_std ~ 1 + water_cent + shade_cent + water_cent:shade_cent,
           data = d, family = gaussian,
           prior = c(prior(normal(0.5, 0.25), class = Intercept),
                     prior(normal(0, 0.25), class = b, coef = water_cent),
                     prior(normal(0, 0.25), class = b, coef = shade_cent),
                     prior(normal(0, 0.25), class = b, coef = "water_cent:shade_cent"),
                     prior(exponential(1), class = sigma)),
           iter = 4000, warmup = 2000, chains = 4, cores = 4, seed = 1234)

summary(mod)
#>  Family: gaussian 
#>   Links: mu = identity; sigma = identity 
#> Formula: blooms_std ~ 1 + water_cent + shade_cent + water_cent:shade_cent 
#>    Data: d (Number of observations: 27) 
#> Samples: 4 chains, each with iter = 4000; warmup = 2000; thin = 1;
#>          total post-warmup samples = 8000
#> 
#> Population-Level Effects: 
#>                       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
#> Intercept                 0.36      0.03     0.30     0.41 1.00     8764
#> water_cent                0.21      0.03     0.14     0.27 1.00     9244
#> shade_cent               -0.11      0.03    -0.18    -0.05 1.00    10097
#> water_cent:shade_cent    -0.14      0.04    -0.22    -0.06 1.00     9670
#>                       Tail_ESS
#> Intercept                 5307
#> water_cent                6043
#> shade_cent                5108
#> water_cent:shade_cent     5938
#> 
#> Family Specific Parameters: 
#>       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
#> sigma     0.14      0.02     0.11     0.19 1.00     6132     5473
#> 
#> 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).

Created on 2020-09-20 by the reprex package (v0.3.0)

In the exercise, we are asked to use new priors that constrain water_cent to be positive and shade_cent to be negative. I can constrain water_cent by using a lognormal distribution. However, I have not been able to figure out how to define a prior that will constrain shade_cent to be negative. I have tried setting boundaries:

mod <- brm(blooms_std ~ 1 + water_cent + shade_cent + water_cent:shade_cent,
           data = d, family = gaussian,
           prior = c(prior(normal(0.5, 0.25), class = Intercept),
                     prior(normal(0, 0.25), class = b, coef = water_cent, lb = 0),
                     prior(normal(0, 0.25), class = b, coef = shade_cent, ub = 0),
                     prior(normal(0, 0.25), class = b, coef = "water_cent:shade_cent"),
                     prior(exponential(1), class = sigma)),
           iter = 4000, warmup = 2000, chains = 4, cores = 4, seed = 1234)
#> Error: Argument 'coef' may not be specified when using boundaries.

But as you can see, I can’t specify priors for specific coefficients while setting boundaries.

Following this discussion, I tried to use a positive lognormal prior and multiply by -1, but that fails also:

mod <- brm(blooms_std ~ 1 + water_cent + (-1 * shade_cent) + water_cent:shade_cent,
           data = d, family = gaussian,
           prior = c(prior(normal(0.5, 0.25), class = Intercept),
                     prior(lognormal(0, 1), class = b, coef = water_cent),
                     prior(lognormal(0, 1), class = b, coef = shade_cent),
                     prior(normal(0, 0.25), class = b, coef = "water_cent:shade_cent"),
                     prior(exponential(1), class = sigma)),
           iter = 4000, warmup = 2000, chains = 4, cores = 4, seed = 1234)
#> Error: The following priors do not correspond to any model parameter: 
#> Intercept ~ normal(0.5, 0.25)
#> b_shade_cent ~ lognormal(0, 1)
#> Function 'get_prior' might be helpful to you.

Is there any documentation I have missed for specifying a prior distribution that will constrain a specific coefficient to be negative?


  • Operating System: macOS 10.14.6
  • brms Version: 2.13.9
sessioninfo::session_info()
#> ─ Session info ─────────────────────────────────────────────────────────────────────────────
#>  setting  value                       
#>  version  R version 4.0.2 (2020-06-22)
#>  os       macOS Mojave 10.14.6        
#>  system   x86_64, darwin17.0          
#>  ui       RStudio                     
#>  language (EN)                        
#>  collate  en_US.UTF-8                 
#>  ctype    en_US.UTF-8                 
#>  tz       America/Chicago             
#>  date     2020-09-20                  
#> 
#> ─ Packages ─────────────────────────────────────────────────────────────────────────────────
#>  package        * version    date       lib source                                
#>  abind            1.4-5      2016-07-21 [1] CRAN (R 4.0.2)                        
#>  assertthat       0.2.1      2019-03-21 [1] CRAN (R 4.0.2)                        
#>  backports        1.1.10     2020-09-15 [1] CRAN (R 4.0.2)                        
#>  base64enc        0.1-3      2015-07-28 [1] CRAN (R 4.0.2)                        
#>  bayesplot        1.7.2      2020-05-28 [1] CRAN (R 4.0.2)                        
#>  bridgesampling   1.0-0      2020-02-26 [1] CRAN (R 4.0.2)                        
#>  brms           * 2.13.9     2020-09-20 [1] Github (paul-buerkner/brms@b4702f3)   
#>  Brobdingnag      1.2-6      2018-08-13 [1] CRAN (R 4.0.2)                        
#>  callr            3.4.4      2020-09-07 [1] CRAN (R 4.0.2)                        
#>  cli              2.0.2      2020-02-28 [1] CRAN (R 4.0.2)                        
#>  clipr            0.7.0      2019-07-23 [1] CRAN (R 4.0.2)                        
#>  coda             0.19-3     2019-07-05 [1] CRAN (R 4.0.2)                        
#>  codetools        0.2-16     2018-12-24 [1] CRAN (R 4.0.2)                        
#>  colorspace       1.4-1      2019-03-18 [1] CRAN (R 4.0.2)                        
#>  colourpicker     1.1.0      2020-09-14 [1] CRAN (R 4.0.2)                        
#>  crayon           1.3.4      2017-09-16 [1] CRAN (R 4.0.2)                        
#>  crosstalk        1.1.0.1    2020-03-13 [1] CRAN (R 4.0.2)                        
#>  curl             4.3        2019-12-02 [1] CRAN (R 4.0.1)                        
#>  desc             1.2.0      2018-05-01 [1] CRAN (R 4.0.2)                        
#>  devtools       * 2.3.2      2020-09-18 [1] CRAN (R 4.0.2)                        
#>  digest           0.6.25     2020-02-23 [1] CRAN (R 4.0.2)                        
#>  dplyr            1.0.2      2020-08-18 [1] CRAN (R 4.0.2)                        
#>  DT               0.15       2020-08-05 [1] CRAN (R 4.0.2)                        
#>  dygraphs         1.1.1.6    2018-07-11 [1] CRAN (R 4.0.2)                        
#>  ellipsis         0.3.1      2020-05-15 [1] CRAN (R 4.0.2)                        
#>  evaluate         0.14       2019-05-28 [1] CRAN (R 4.0.1)                        
#>  fansi            0.4.1      2020-01-08 [1] CRAN (R 4.0.2)                        
#>  fastmap          1.0.1      2019-10-08 [1] CRAN (R 4.0.2)                        
#>  fs               1.5.0      2020-07-31 [1] CRAN (R 4.0.2)                        
#>  generics         0.0.2      2018-11-29 [1] CRAN (R 4.0.2)                        
#>  ggplot2        * 3.3.2      2020-06-19 [1] CRAN (R 4.0.2)                        
#>  ggridges         0.5.2      2020-01-12 [1] CRAN (R 4.0.2)                        
#>  glue             1.4.2      2020-08-27 [1] CRAN (R 4.0.2)                        
#>  gridExtra        2.3        2017-09-09 [1] CRAN (R 4.0.2)                        
#>  gtable           0.3.0      2019-03-25 [1] CRAN (R 4.0.2)                        
#>  gtools           3.8.2      2020-03-31 [1] CRAN (R 4.0.2)                        
#>  highr            0.8        2019-03-20 [1] CRAN (R 4.0.2)                        
#>  htmltools        0.5.0      2020-06-16 [1] CRAN (R 4.0.2)                        
#>  htmlwidgets      1.5.1      2019-10-08 [1] CRAN (R 4.0.2)                        
#>  httpuv           1.5.4      2020-06-06 [1] CRAN (R 4.0.2)                        
#>  igraph           1.2.5      2020-03-19 [1] CRAN (R 4.0.2)                        
#>  inline           0.3.16     2020-09-06 [1] CRAN (R 4.0.2)                        
#>  jsonlite         1.7.1      2020-09-07 [1] CRAN (R 4.0.2)                        
#>  knitr            1.29       2020-06-23 [1] CRAN (R 4.0.2)                        
#>  later            1.1.0.1    2020-06-05 [1] CRAN (R 4.0.2)                        
#>  lattice          0.20-41    2020-04-02 [1] CRAN (R 4.0.2)                        
#>  lifecycle        0.2.0      2020-03-06 [1] CRAN (R 4.0.2)                        
#>  loo              2.3.1      2020-07-14 [1] CRAN (R 4.0.2)                        
#>  magrittr         1.5        2014-11-22 [1] CRAN (R 4.0.2)                        
#>  markdown         1.1        2019-08-07 [1] CRAN (R 4.0.2)                        
#>  MASS             7.3-51.6   2020-04-26 [1] CRAN (R 4.0.2)                        
#>  Matrix           1.2-18     2019-11-27 [1] CRAN (R 4.0.2)                        
#>  matrixStats      0.56.0     2020-03-13 [1] CRAN (R 4.0.2)                        
#>  memoise          1.1.0      2017-04-21 [1] CRAN (R 4.0.2)                        
#>  mime             0.9        2020-02-04 [1] CRAN (R 4.0.2)                        
#>  miniUI           0.1.1.1    2018-05-18 [1] CRAN (R 4.0.2)                        
#>  munsell          0.5.0      2018-06-12 [1] CRAN (R 4.0.2)                        
#>  mvtnorm          1.1-1      2020-06-09 [1] CRAN (R 4.0.2)                        
#>  nlme             3.1-148    2020-05-24 [1] CRAN (R 4.0.2)                        
#>  pillar           1.4.6      2020-07-10 [1] CRAN (R 4.0.2)                        
#>  pkgbuild         1.1.0      2020-07-13 [1] CRAN (R 4.0.2)                        
#>  pkgconfig        2.0.3      2019-09-22 [1] CRAN (R 4.0.2)                        
#>  pkgload          1.1.0      2020-05-29 [1] CRAN (R 4.0.2)                        
#>  plyr             1.8.6      2020-03-03 [1] CRAN (R 4.0.2)                        
#>  prettyunits      1.1.1      2020-01-24 [1] CRAN (R 4.0.2)                        
#>  processx         3.4.4      2020-09-03 [1] CRAN (R 4.0.2)                        
#>  promises         1.1.1      2020-06-09 [1] CRAN (R 4.0.2)                        
#>  ps               1.3.4      2020-08-11 [1] CRAN (R 4.0.2)                        
#>  purrr            0.3.4      2020-04-17 [1] CRAN (R 4.0.2)                        
#>  R6               2.4.1      2019-11-12 [1] CRAN (R 4.0.2)                        
#>  Rcpp           * 1.0.5      2020-07-06 [1] CRAN (R 4.0.2)                        
#>  RcppParallel     5.0.2      2020-06-24 [1] CRAN (R 4.0.2)                        
#>  remotes        * 2.2.0      2020-07-21 [1] CRAN (R 4.0.2)                        
#>  reprex           0.3.0.9001 2020-09-20 [1] Github (tidyverse/reprex@d3fc4b8)     
#>  reshape2         1.4.4      2020-04-09 [1] CRAN (R 4.0.2)                        
#>  rethinking     * 2.13       2020-08-27 [1] Github (rmcelreath/rethinking@3b48ec8)
#>  rlang            0.4.7      2020-07-09 [1] CRAN (R 4.0.2)                        
#>  rmarkdown        2.3.8      2020-09-20 [1] Github (rstudio/rmarkdown@fbd0b1f)    
#>  rprojroot        1.3-2      2018-01-03 [1] CRAN (R 4.0.2)                        
#>  rsconnect        0.8.16     2019-12-13 [1] CRAN (R 4.0.2)                        
#>  rstan          * 2.21.2     2020-07-27 [1] CRAN (R 4.0.2)                        
#>  rstantools       2.1.1      2020-07-06 [1] CRAN (R 4.0.2)                        
#>  rstudioapi       0.11       2020-02-07 [1] CRAN (R 4.0.2)                        
#>  scales           1.1.1      2020-05-11 [1] CRAN (R 4.0.2)                        
#>  sessioninfo      1.1.1      2018-11-05 [1] CRAN (R 4.0.2)                        
#>  shape            1.4.5      2020-09-13 [1] CRAN (R 4.0.2)                        
#>  shiny            1.5.0      2020-06-23 [1] CRAN (R 4.0.2)                        
#>  shinyjs          2.0.0      2020-09-09 [1] CRAN (R 4.0.2)                        
#>  shinystan        2.5.0      2018-05-01 [1] CRAN (R 4.0.2)                        
#>  shinythemes      1.1.2      2018-11-06 [1] CRAN (R 4.0.2)                        
#>  StanHeaders    * 2.21.0-6   2020-08-16 [1] CRAN (R 4.0.2)                        
#>  stringi          1.5.3      2020-09-09 [1] CRAN (R 4.0.2)                        
#>  stringr          1.4.0      2019-02-10 [1] CRAN (R 4.0.2)                        
#>  styler           1.3.2      2020-02-23 [1] CRAN (R 4.0.2)                        
#>  testthat       * 2.3.2      2020-03-02 [1] CRAN (R 4.0.2)                        
#>  threejs          0.3.3      2020-01-21 [1] CRAN (R 4.0.2)                        
#>  tibble           3.0.3      2020-07-10 [1] CRAN (R 4.0.2)                        
#>  tidyselect       1.1.0      2020-05-11 [1] CRAN (R 4.0.2)                        
#>  usethis        * 1.6.3      2020-09-17 [1] CRAN (R 4.0.2)                        
#>  V8               3.2.0      2020-06-19 [1] CRAN (R 4.0.2)                        
#>  vctrs            0.3.4      2020-08-29 [1] CRAN (R 4.0.2)                        
#>  withr            2.2.0      2020-04-20 [1] CRAN (R 4.0.2)                        
#>  xfun             0.17       2020-09-09 [1] CRAN (R 4.0.2)                        
#>  xtable           1.8-4      2019-04-21 [1] CRAN (R 4.0.2)                        
#>  xts              0.12.1     2020-09-09 [1] CRAN (R 4.0.2)                        
#>  yaml             2.2.1      2020-02-01 [1] CRAN (R 4.0.2)                        
#>  zoo              1.8-8      2020-05-02 [1] CRAN (R 4.0.2)                        
#> 
#> [1] /Library/Frameworks/R.framework/Versions/4.0/Resources/library

Not the ideal fix, but if you literally just wanted to get something so that you can move on with exercises, you could specify the dreaded uniform prior, between some negative value and 0. Then it would be bounded at zero. This is almost certainly not what Richard would like you to do, but perhaps it could be a quick fix to move forward, and you can come back when someone more experienced has a better solution!?

I don’t think you can apply an upper-bound to a coefficient parameter in brms (@paul.buerkner might have some good insight here), but could you instead use the lognormal prior and multiply the data by -1 before passing it to the model:

mod <- 
    d %>% mutate(shade_cent = -1*shade_cent) %>%
    brm(blooms_std ~ 1 + water_cent + shade_cent + water_cent:shade_cent,
          data = ., family = gaussian,
           prior = c(prior(normal(0.5, 0.25), class = Intercept),
                     prior(lognormal(0, 1), class = b, coef = water_cent),
                     prior(lognormal(0, 1), class = b, coef = shade_cent),
                     prior(normal(0, 0.25), class = b, coef = "water_cent:shade_cent"),
                     prior(exponential(1), class = sigma)),
           iter = 4000, warmup = 2000, chains = 4, cores = 4, seed = 1234)