paul-buerkner / brms

brms R package for Bayesian generalized multivariate non-linear multilevel models using Stan
https://paul-buerkner.github.io/brms/
GNU General Public License v2.0
1.28k stars 184 forks source link

Priors on sigma not getting passed to Stan code #213

Closed mvuorre closed 7 years ago

mvuorre commented 7 years ago

There appears to be a problem when trying to specify priors on the standard deviation of a normal response distribution. It looks like the priors aren't getting passed to the right parameters, but instead the effects on sigma get the priors that were assigned to effects on the mean of the distribution. For details see below:

Example data and model from https://cran.r-project.org/web/packages/brms/vignettes/brms_distreg.html

library(brms)
## Loading required package: Rcpp
## Loading required package: ggplot2
## Loading 'brms' package (version 1.6.0.9000). Useful instructions 
## can be found by typing help('brms'). A more detailed introduction 
## to the package is available through vignette('brms_overview').
set.seed(1)
group <- rep(c("treat", "placebo"), each = 30)
symptom_post <- c(rnorm(30, mean = 1, sd = 2), rnorm(30, mean = 0, sd = 1))
dat1 <- data.frame(group, symptom_post)
head(dat1)
##   group symptom_post
## 1 treat   -0.2529076
## 2 treat    1.3672866
## 3 treat   -0.6712572
## 4 treat    4.1905616
## 5 treat    1.6590155
## 6 treat   -0.6409368

I want to implement priors on all four parameters, and treat the intercept as an ordinary pop. level effect.

Formula <- bf(symptom_post ~ 0 + intercept + group, 
              sigma ~ 0 + intercept + group)
get_prior(Formula, dat1, family=gaussian())
##   prior class       coef group nlpar bound
## 1           b                             
## 2           b grouptreat                  
## 3           b  intercept                  
## 4           b                  sigma      
## 5           b grouptreat       sigma      
## 6           b  intercept       sigma

I then specify some arbitrary priors to illustrate the problem. The call to get_prior() suggests that I should specify the priors on class = "b", nlpar = "sigma".

Prior <- c(
    set_prior("student_t(10, 0, 1)", class = "b", coef = "intercept"),
    set_prior("student_t(20, 0, 2)", class = "b", coef = "grouptreat"),
    set_prior("student_t(30, 0, 3)", class = "b", coef = "intercept", nlpar = "sigma"),
    set_prior("student_t(40, 0, 4)", class = "b", coef = "grouptreat", nlpar = "sigma")
)
make_stancode(Formula, family=gaussian(), data = dat1, 
              prior = Prior, sample_prior = T)
## // generated with brms 1.6.0.9000
## functions { 
## } 
## data { 
##   int<lower=1> N;  // total number of observations 
##   vector[N] Y;  // response variable 
##   int<lower=1> K;  // number of population-level effects 
##   matrix[N, K] X;  // population-level design matrix 
##   int<lower=1> K_sigma;  // number of population-level effects 
##   matrix[N, K_sigma] X_sigma;  // population-level design matrix 
##   int prior_only;  // should the likelihood be ignored? 
## } 
## transformed data { 
## } 
## parameters { 
##   vector[K] b;  // population-level effects 
##   vector[K_sigma] b_sigma;  // population-level effects 
## } 
## transformed parameters { 
## } 
## model { 
##   vector[N] mu; 
##   vector[N] sigma; 
##   mu = X * b; 
##   sigma = X_sigma * b_sigma; 
##   for (n in 1:N) { 
##     sigma[n] = exp(sigma[n]); 
##   } 
##   // prior specifications 
##   b[1] ~ student_t(10, 0, 1); 
##   b[2] ~ student_t(20, 0, 2); 
##   b_sigma[1] ~ student_t(10, 0, 1); 
##   b_sigma[2] ~ student_t(20, 0, 2); 
##   // likelihood contribution 
##   if (!prior_only) { 
##     Y ~ normal(mu, sigma); 
##   } 
## } 
## generated quantities { 
##   real prior_b_1; 
##   real prior_b_2; 
##   real prior_b_sigma_1; 
##   real prior_b_sigma_2; 
##   // additionally draw samples from priors 
##   prior_b_1 = student_t_rng(10,0,1); 
##   prior_b_2 = student_t_rng(20,0,2); 
##   prior_b_sigma_1 = student_t_rng(10,0,1); 
##   prior_b_sigma_2 = student_t_rng(20,0,2); 
## }

It seems that the sigmas are not getting the priors assigned on them, but instead they get the priors for effects on the means. Estimation gives an error:

fit <- brm(formula = Formula, 
          family = gaussian(),
          data = dat1, 
          prior = Prior,
          sample_prior = T,
          cores = 4)
## Compiling the C++ model
## Start sampling
## Error in warning(..., call. = FALSE): formal argument "call." matched by multiple actual arguments
stancode(fit)
## Error in stancode(fit): object 'fit' not found

I've tried this with the current CRAN version of brms, and also the github version:

sessionInfo()
## R version 3.3.3 (2017-03-06)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
## Running under: macOS Sierra 10.12.4
## 
## 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   base     
## 
## other attached packages:
## [1] brms_1.6.0.9000 ggplot2_2.2.1   Rcpp_0.12.10   
## 
## loaded via a namespace (and not attached):
##  [1] knitr_1.15.19        magrittr_1.5         rstantools_1.2.0    
##  [4] munsell_0.4.3        lattice_0.20-35      colorspace_1.3-2    
##  [7] R6_2.2.0             dplyr_0.5.0          stringr_1.2.0       
## [10] plyr_1.8.4           tools_3.3.3          bayesplot_1.1.0     
## [13] parallel_3.3.3       grid_3.3.3           nlme_3.1-131        
## [16] gtable_0.2.0         loo_1.1.0            coda_0.19-1         
## [19] DBI_0.6-1            matrixStats_0.52.1   StanHeaders_2.14.0-1
## [22] htmltools_0.3.5      assertthat_0.1       yaml_2.1.14         
## [25] abind_1.4-5          lazyeval_0.2.0       rprojroot_1.2       
## [28] digest_0.6.12        tibble_1.3.0         rstan_2.14.2        
## [31] gridExtra_2.2.1      inline_0.3.14        evaluate_0.10       
## [34] rmarkdown_1.4.0.9001 stringi_1.1.5        scales_0.4.1        
## [37] backports_1.0.5      stats4_3.3.3
paul-buerkner commented 7 years ago

You are right, I will fix that! Thanks!

paul-buerkner commented 7 years ago

On the dev version of brms on github, it should now be working correctly! :-)