urgently need your help! In doing bayesian modelling, how do we set prior for the standard deviation of regression slope, so-called hyper-parameter. For instance, y = b0 + b1x1+b2x2+error, where b2 ~ N(0, tau), and tau ~ U(1,1). so, how to set the prior for the tau? Many thanks!
In the following, I modified the example codes in ?rstan so that the modified codes use a hyper parameter as follows
Y \sim \text{Normal}(\mu,1),
\mu \sim \text{Normal}(hyper, 10),
hyper \sim \text{Normal}(-1,1)
So, we can describe the prior of hyper parameter in the ususal way.
stanmodelcode <- "
data {
int<lower=0> N;
real y[N];
}
parameters {
real mu;
real hyper;
}
model {
target += normal_lpdf(hyper | -1, 1);
target += normal_lpdf(mu | hyper, 10);
target += normal_lpdf(y | mu, 1);
}
"
y <- rnorm(20)
dat <- list(N = 20, y = y);
fit <- stan(model_code = stanmodelcode, model_name = "example",
data = dat, iter = 2012, chains = 3, verbose = FALSE,
sample_file = file.path(tempdir(), 'norm.csv'))
print(fit)
1 Like
Many thanks @Jean_Billie! Is it possible to set those in the function ‘brm’? I would like to set prior for hyper-parameter the brm. Actually I am not familiar with rstan syntax. Thanks.
Sorry, I cannot. I am not familiar with brm package, but only Stan.
You can do this in brms by passing Stan code to brm() via two arguments: prior and stanvars (which adds tau as a parameter). Effectively this is exactly what @Jean_Billie recommends, but leans on brms to create as much of the Stan code as possible.
bpriors <- prior(normal(0, 2), class = "Intercept") +
prior(normal(0, 4), class = "b", coef = "x1") +
prior(normal(0, tau), class = "b", coef = "x2")+ # here is where we add tau as a hyperparameter for the standard deviation of x2
prior("target += exponential_lpdf(tau | 1)", check = FALSE) + # here is where we define the prior for tau; I picked exponential(1)
prior(exponential(1), class = "sigma")
stanvars <- stanvar(scode = " real<lower=0> tau;", # here is where we add the parameter for tau
block = "parameters")
bmod <- brm(y ~ 1 + x1 + x2, family = gaussian(),
data = dat, chains = 2,
prior = bpriors, stanvars = stanvars) # this line passes the prior & modified parameter blocks into brm()
summary(bmod)
(ht: @paul.buerkner for posting a similar example https://github.com/paul-buerkner/brms/issues/459#issuecomment-399041581)
If your actual model is more complicated than the example you provided, you may find it easier to use the rethinking package. Its syntax is very similar to how you wrote your question. rethinking::ulam is the workhorse function, and within that alist wraps the model (<-) & priors (~)
rmod <- ulam(
alist(
y ~ dnorm(mu, sig),
mu <- a + b*x1 + c*x2,
a ~ dnorm(0, 2),
b ~ dnorm(0, 4),
c ~ dnorm(0, tau),
tau ~ dexp(1),
sig ~ dexp(1)
),
data = dat, chains = 2
)
precis(rmod)
Finally, this code wouldn’t be complete without the data needed to make it reproducible. Here’s one way to generate a dummy dataset:
b0 <- 2
b1 <- 0.5
b2 <- 0.75
sigma <- 1
n <- 50
x1 <- runif(n)
x2 <- runif(n)
lp <- b0 + b1 * x1 + b2 * x2
y <- rnorm(n, lp, sigma)
dat <- data.frame(x1, x2, y)
2 Likes