Error applying horseshoe prior in Poisson glm

Interfaces brms
1

Please also provide the following information in addition to your question:

  • Operating System: R 3.4.4 on Ubuntu 18.04.1
  • brms Version: 2.4.4

Hello,

I’m trying to fit a Poisson glm to my data, and I would like to apply a normal prior to some parameters in my model, and a horseshoe prior to others.
My model is parameterized as follows:

 hsForm = bf(count ~ alpha + eta1 + eta2 + log(nucGen), 
             alpha ~ 1,  
             eta1 ~ 0 + <myNormalPriorParams>,
             eta2 ~ 0 + <myHorseshoePriorParams>,
             nl=TRUE,
             family=poisson())

pr = c(prior(normal(0,3), nlpar="eta1"), 
       prior(horseshoe(), nlpar="eta2"), 
       prior(normal(-23,5), nlpar="alpha"))

When I try to fit this model to my data, I get the following error:

SYNTAX ERROR, MESSAGE(S) FROM PARSER:

variable "sigma" does not exist.
  error in 'modele86151a69e5_filee865f8cf93' at line 55, column 105
  -------------------------------------------------
    53: } 
    54: transformed parameters { 
    55:   vector[K_eta2] b_eta2 = horseshoe(zb_eta2, hs_local_eta2, hs_global_eta2, hs_scale_global_eta2 * sigma, hs_scale_slab_eta2^2 * hs_c2_eta2); 
                                                                                                                ^
    56: } 
  -------------------------------------------------

Error in stanc(model_code = paste(program, collapse = "\n"), model_name = model_cppname,  : 
  failed to parse Stan model 'filee865f8cf93' due to the above error.

Of course if I use the gaussian family, the model compiles, because sigma exists. I really do need to use Poisson, however. How can I parameterize my model to apply the horseshoe shrinkage prior on some parameters, but a normal on others, while using a Poisson glm?

Thanks!

2

That looks like a bug. Can you provide a reproducible example?

3

I would be happy to do that. I’m pretty new to this stuff, though. When you ask for a reproducible example, do you just need more elaboration on my model’s parameters and a couple example rows in a typical dataset?

If so, here is some example data:

background binVar nucGen count Signal1 Signal2 Signal3
wt 0 365690700 1 1.96040799009779 1.2120092926712 1.27138013020913
wt 0 310059900 0 1.60797118707689 1.11098600125775 0.363046518815955
wt 1 58903200 2 1.19423577475165 0.853740870198064 2.46228178036029
wt 1 303515100 0 1.72044334660999 1.48068239325765 1.89201291098705
wt 0 67084200 0 1.48112731648809 1.25606462797278 0.83583370217392
mutant 0 522765900 1 2.47680830484492 1.59257122660975 2.04069252530003
mutant 1 31087800 1 0.66422015479385 0.12771625892178 1.4971909598719
mutant 0 672478200 0 2.07548017856752 1.76626031532154 1.51622851733695
mutant 1 258519600 0 0.065752965294899 -0.698184227708412 0.581839196256855
mutant 0 400050900 0 3.11353662066563 2.61093865000118 1.97723216922194

And a more detailed model description:

hsForm = bf(count ~ alpha + eta1 + eta2 + log(nucGen), # nucGen is exposure
             alpha ~ 1,  
             eta1 ~ 0 + background + binVar + Signal1 + Signal2 + Signal3,
             eta2 ~ 0 + background:binVar + background:Signal1 + background:Signal2 + background:Signal3,
             nl=TRUE,
             family=poisson())

pr = c(prior(normal(0,3), nlpar="eta1"), 
       prior(horseshoe(), nlpar="eta2"), 
       prior(normal(-23,5), nlpar="alpha"))

fit = brm(hsForm, 
          data=myDat, 
          prior=pr, cores=1, chains=1, 
          iter=2000, control = list(adapt_delta=0.99))

Please let me know if you need anything else from me.

4

Thanks! I need some code, I can just run after pasting it into my editor. Actually, I would like a minimal reproducible example that makes it as simple as possible, but still triggers the error.

For you code to be reproducible I actually needed some dummy data defined directly in the R code, but for this example, I will just try to go from what I see.

I report back as soon as I know what’s going on.

5

The problem is now fixed in the dev version of brms to be installed via

if (!requireNamespace("devtools")) {
  install.packages("devtools")
}
devtools::install_github("paul-buerkner/brms")
1 Like
6

You’re awesome! I updated to the latest dev version, and my model now compiles and is gently sampling away!

Thanks!