In frequentist linear regression and generalized linear models, we typically check model assumptions using residual diagnostics.
I’m aware that metrics like ELPD (expected log predictive density) can be used to compare models, and that posterior predictive checks are a standard way to assess model fit and conformance to data. However, do these approaches specifically confirm whether the assumption of normally distributed residuals holds in Bayesian linear regression, or whether the selected distribution is appropriate for residuals in generalized linear models? How should these assumption checks be performed and interpreted in PyMC?
Having said that, the expectation part of a GLM is usually much more important for downstream inference and you can often think of the error term as just greasing the wheels (because you can’t do a linear fit of an overdetermined system). Violations of normality are usually not a big deal unless there’s something like residual correlation among errors (e.g., in unmodeled time series dependence) or there is really bad misspecification in the tails or if you care about error terms in individual predictions going forward (the latter is where you really want to get error terms right, especially when there are lives or even just money or time on the line).
In the frequentist case, we can go further and fully decouple the computations. First compute the MLE for \beta assuming \sigma = 1, then estimate \sigma^2 from the squared residuals given the MLE \widehat{\beta}.