Say I have a set of data for three quantities $x$, $z$, and $y$ at various locations, $i$. I then fit a linear model,
$$ y_i = a x_i + b z_i + \epsilon_i $$
to this data, assuming $\epsilon_i \sim N(0, \sigma)$, with the value of $\sigma$ unknown.
I'm then interested in predicting $y$ at some new locations $j \in \{1,2,...n\}$, and most importantly I'd like to quantify uncertainty in that prediction. We have $x$ and $z$ values in these new locations. The question is how do I generate something like a histogram of $y$ at each of these locations that account for all uncertainty. Presumably, simply saying
$$y_j \sim N( \hat a x_j + \hat b z_j, \hat \sigma ),$$
where $\hat a$ and $\hat b$ are the estimates of the coefficients from the linear regression model and $\hat \sigma$ is the residual standard error, is not account for all of the uncertainty, because it assumes the estimated coefficients are true.
So for example how might I plot a histogram for each $y_j$, given the model that accounts for all uncertainty. I think this is related to prediction, vs. confidence intervals for which there are a lot of tutorials online, but I'm interested in the full distribution. The R code below does roughly what I want to do, but I don't think it accounts for enough uncertainty.
## data for fitting
n = 100;
x = rnorm(n, 1, .5);
z = rnorm(n, .5, .1);
y = rnorm(n, mean = 2*x - 6*z, sd = 1)
mod = lm(y~x+z)
## new data
xn = rnorm(n, 1, .5);
zn = rnorm(n, .5, .1);
ytrue = rnorm(n, mean = 2*xn - 6*zn, sd = 1)
## predict with uncertainty
j=2; #choose a site to plot uncertainty
y_pred_samp = rnorm(500,
mean = mod$coef[1]*xn[j] + mod$coef[2]*zn[j],
sd = sigma(mod))
hist(y_pred_samp, main = paste('true value is: ', ytrue[j]))