# Bayesian Linear Regression - Creating a distribution for a new prediction

I'm using MCMC to fit a linear regression model with the end goal of making predictions for new observations. See reproducible example below:

library(ggplot2)
library(MCMCpack)
library(dplyr)
data(mpg)
attach(mpg)

# Variance of hwy for each trans category
# Some categories e.g. auto(s5) have a bigger variance
# How can I make the distribution for a new prediction reflect this?
vars <- mpg %>%
group_by(trans) %>%
summarise(hwy_var=var(hwy)) %>%
ungroup()

# Fit model
bm <- MCMCregress(hwy ~ year + cyl + trans + cty)

# Want to find distribution of hwy for new data point
new_obs <- data.frame(year=1999, cyl=4, trans="auto(l6)", cty=16)


What I want to capture in predicting hwy for the new_obs is a predictive distribution, taking into account its predictor values. For example, auto(l6) is a category with a larger hwy variance (as seen in vars) so I'd like the predictive distribution for this point to be wider. Is this possible with Bayesian statistics? Or any other method?

You are probably looking for the posterior predictive distribution. So, what you could do is the following. Suppose that $$(\beta^{(s)},\sigma^{(s)})$$ is the $$s$$th posterior sample of your parameter vector, consisting of the regression coefficients including the intercept and the residual standard deviation. Further, let $$\mathbf x$$ be the set of predictor values for which you want to have the prediction (e.g., (1999, 4,...)), where the first element of $$\mathbf x$$ is a one. Using your $$s=1,2,...,S$$ posterior samples of $$(\beta^{(s)},\sigma^{(s)})$$, calculate $$\hat y^{(s)} = \mathbf x^\top \beta^{(s)}$$ and then draw $$\tilde y^{(s)}$$ from $$\text{Normal}(\hat y^{(s)}, \sigma^{(s)})$$. This will give you $$S$$ draws from the posterior predictive distribution at $$\mathbf x$$.