Could you please check if I specify the model by brms correctly?

Interfaces brms
1

Hello,

I want to specify two model, one is logistic model,
(1) y \sim Binominal(p)
logit(p)=log(p/1-p)=\beta_0+X\beta+\alpha
where \beta \sim N(0,10000), \alpha \sim N(0,0.001)
and I specify by:

prior_logit<-prior(normal(0,10000),class=b,coef=f1)+prior(normal(0,10000),class=b,coef=f2)+
  prior(normal(0,10000),class=b,coef=f3)+prior(normal(0,10000),class=b,coef=f4)+
  prior(normal(0,10000),class=b,coef=f5)+prior(normal(0,10000),class=b,coef=f6)+
  prior(normal(0,0.001),class=Intercept)

model_logit<-brms::brm(y~f1+f2+f3+f4+f5+f6,
            data=s,family=bernoulli,warmup=1000,iter=5000,chains=1,prior=prior_test)

(2) model 2 is logistic random effect model,
y \sim Binominal(p)
logit(p)=log(p/1-p)=\beta_0+X\beta+u_j(i)
where \beta \sim N(0,10000), u_j(i) is grouping factor, u_j \sim N(0,\tau), \tau \sim gamma(0.001,0.001)

and I specify in brms:

prior_re<prior(normal(0,10000),class=b,coef=f1)+prior(normal(0,10000),class=b,coef=f2)+
prior(normal(0,10000),class=b,coef=f3)+prior(normal(0,10000),class=b,coef=f4)+
prior(normal(0,10000),class=b,coef=f5)+prior(normal(0,10000),class=b,coef=f6)+
prior(gamma(0.001,0.001),class=sd)

model_re<-brms::brm(y~f1+f2+f3+f4+f5+f6+(1|u),
            data=s,family=bernoulli,warmup=1000,iter=5000,chains=1,prior=prior_re)

Am I specify the model correctly? Thanks!

1 Like
2

Generally speaking having such wide priors is not good for many reasons (yours are basically flat). Please see here for some background info:

In the end, you should always simulate from your priors, so you see that they are sane.

1 Like
3

As @torkar said, you shouldn’t be using such flat priors, especially in a logistic regression setting. What those priors imply is that a standardized regression coefficient as large as 1000 (in absolute value) has 30% probability of occurring:

2 * pnorm(1000, 0, 1000, lower.tail=FALSE)
[1] 0.3173105

Even a coefficient of at least 10 occurs almost certainly (99%), and note that the corresponding odds ratio would be exp(10) = 22026. In Statistical Rethinking, @richard_mcelreath has nice plots that show that even a normal prior with standard deviation of 10 would create spikes at 0 and 1 in logistic regression, which is totally the opposite of flat.

It’s also quite strange that you set very tight priors on the intercept. Would you have a justification for that?

Coming back to @torkar’s suggestion, prior predictive checks would give you a good idea on how to set priors. For that, Statistical Rethinking 2nd edition is a great resource.

3 Likes
4

Hello,

Regarding the models, there are some small things:

  • You want a binomial but use a bernoulli. While the bernoulli is a special case binomial, you don’t make clear that which one you want (or I am missing something).
  • Your prior for the first model has two two different names (prior_logit vs prior_test)
  • Think about running multiple chains for better rhat estimation. Using options(mc.cores = parallel::detectCores()) will let R run the chains in parallel if you have multiple cores.

And I agree with the priors being very wide (and rather tight for the intercept) as stated before.

2 Likes
5

Hello everyone,

I have read your suggestions…Is it more appropriate to set parameter prior N(0,10e6) and N(0,1) for the intercept? Or use the default student prior student_t(3, 0, 10) for Intercept?
And I will use invgamma(0.001,0.001) for the random term in model2.

Thank for all of your help!

6

No, I don’t think those are appropriate. To see how N(0, 10e6) behaves, just try the following in R, which extracts 10 random numbers from that distribution:

rnorm(10, 0, 10e6)
 [1]   4016654   8621579  -7582142  -3323129  20135648 -23859572  -1489552
 [8] -18364477 -18275152  -7689779

Do coefficients like these make sense for your model? Presumably not. If you try something like rnorm(10, 0, 10) you would get

 [1]  -6.4048409 -20.3726069   8.3001658   3.0496285  -0.1736941   9.4504703
 [7]   6.9147799   6.5376653  -5.9216552   7.3235563

Are these values closer to what you would expect? That’s a way to start choosing priors.

Also, we’ve concentrated mainly on the standard deviation term as it was the wildest one. However, think also about the mean of the normal: if your coefficients are centred (so they have mean zero), then your mean is fine. But if that’s not the case, and you expect them to have a mean of say 50, also the corresponding prior should have a mean of 50 and not 0.

1 Like
7

I think I have understood what you mean…The prior set by me is very arbitraty which can have effects on estimation results? Thanks!

1 Like