Predictions from Bayesian posteriors

Question

I am reading the answer to this post (Prediction based on bayesian model) about new predictions, and I have a question. Could I find beta0 and beta1 distributions in regular way:

 beta0 ~ dnorm(0, 10)
 beta1 ~ dnorm(0, 10)
 sigma ~ dunif(0, 50)
 for (i in 1:N) {
    y[i] ~ dnorm(beta0 + beta1 * x[i], sigma)
 }

And then in R do something like this (beta0 and beta1 would be taken from the posterior distribution (regression) and x_new (1:40) would be an array where I need extrapolation):

beta0 beta1 x_new y_new
1.09  0.01  1     1.1
0.98  0.015 1     0.995
1.08  0.012 1     1.092
...
...
...
1.09  0.01  2     1.1
0.98  0.015 2     1.01
1.08  0.012 2     1.11
...
...
...
1.09  0.01  40     1.49
0.98  0.015 40     1.59
1.08  0.012 40     1.56
...
...
...

Can I use this approach to create Bayes inference? Is it wrong? Why do I need "sigma" in the answer from the link for y_new?

Thank you very much for explanation

Answer 1

Writing your code move verbosely:

for (i in 1:N) {
    mu[i] <- beta0 + beta1 * x[i]
    y[i] ~ dnorm(mu[i], sigma)
}

Your model is a Bayesian linear regression, so as standard linear regression, it predicts the mean $E[y_i|x_i] = \mu_i = \beta_0 + \beta_1 x_i$. If you are interested in the mean $\mu_i$, sure, take the predictions for the mean. However if you want to look at the samples from the conditional distribution of $y_i|x_i$ itself, then you sample from the posterior predictive distribution. In many cases, saving both the mean and samples from the posterior predictive distribution would be useful.