Should my residuals have the same distribution as my likelihood function in a Bayesian linear regression? If I have a standard linear regression model:
$$
y \sim \mathcal{N}(\beta_0 + \beta_1X + ..., \sigma)
$$
With whatever priors on $\beta$'s and $\sigma$. Does that imply that my residuals $r = y - (\beta_0 + \beta_1X + ...)$ should be $r \sim \mathcal{N}(0, \sigma)$? I know that in a traditional least-squares setting, this is a common check. Is it relevant in a Bayesian setting? Further, in traditional least-squares, not meeting the check means prediction intervals would be off as well as confidence intervals for parameters. I can understand how the prediction interval issue would be relevant in the Bayesian context as well, but does the parameter issue apply as well? Are the posterior estimates for the parameters potentially off because my data doesn't match my assumed likelihood distribution?
Ultimately, I'm curious about the relationship between residuals and the likelihood function in a Bayesian context.

In principle what Tim says; however, there are two subtleties
a) Given that your parameter estimates are distributions, what's your definition of "residual"? To get a "point" residual, you could take the model predictions at the MAP or the posterior mean, or else you evaluate the entire posterior (then, your residual is a distribution though). Alternatively, you can calculate quantile residuals from posterior predictive simulations, which creates a point-residual accounting for the posterior uncertainty.
b) The second point is the influence of the priors - if priors dominate, you could force your model to have pretty much any residual distribution.
See comments in https://cran.r-project.org/web/packages/DHARMa/vignettes/DHARMaForBayesians.html
